## Abstract

We construct a two-dimensional continuum model to describe the energetics of shape transitions in fully faceted epitaxial quantum dots (strained islands) via minimization of elastic energy and surface energy at fixed volume. The elastic energy of the island is based on a third-order approximation, enabling us to consider shape transitions between pyramids, domes, multifaceted domes and asymmetric intermediate states. The energetics of the shape transitions are determined by numerically calculating the facet lengths that minimize the energy of a given island type of prescribed island volume. By comparing the energy of different island types with the same volume and analysing the energy surface as a function of the island shape parameters, we determine the bifurcation diagram of equilibrium solutions and their stability, as well as the lowest barrier transition pathway for the island shape as a function of increasing volume. The main result is that the shape transition from pyramid to dome to multifaceted dome occurs through sequential nucleation of facets and involves asymmetric metastable transition shapes. We also explicitly determine the effect of corner energy (facet edge energy) on shape transitions and interpret the results in terms of the relative stability of asymmetric island shapes as observed in experiment.

## 1. Introduction

The spontaneous formation of nanoscale quantum dots (islands) during the epitaxial growth of strained semiconductor films has attracted intense interest as a route to the self-assembly of nano-structured devices for electronic and optoelectronic devices (see [1,2] for reviews). Since the electronic and optical properties are determined by a combination of the elastic strain and the quantum dot shape, a key element in developing self-assembled nanostructures is to understand how the island size and shape are controlled by the material properties and growth parameters.

The self-assembly process during epitaxial growth is driven by reduction in the elastic energy associated with the misfit between the crystal lattice of the film and that of the substrate. Initially, a film grows in a layer-by-layer manner until it reaches a critical layer thickness; then the planar film becomes unstable due to the misfit strain, and self-assembled strained islands begin to form on the wetting layer (the Stranksi–Krastanov growth mode). In the Ge/Si(001) heteroepitaxy system, islands undergo a series of shape transformations as they grow [3–10]. Small unfaceted prepyramid islands grow into a pyramid shape which is bounded by four shallow facets (11° slope, {105} crystallographic facet planes). At larger volume, the pyramid islands transform into a faceted dome shape which includes additional steeper facets with slopes of 22° ({113} facet planes) and 33° ({15 3 23} facet planes). At even larger volume, they will change into other multifaceted shapes which include more steeper facets. These shapes are typically symmetric; however, asymmetric shapes are also observed in experiments [6–10] and simulations [11,12] under some circumstances. Thus, the faceted island shape progression in Ge/Si has two interesting mathematical features: (i) the island shape changes are discontinuous, involving the creation of new crystalline facets, and (ii) the possible existence of asymmetric island shapes.

The important physical components controlling the size and shape of epitaxial islands are elastic strain energy, surface energy and wetting interactions of the thin film with the substrate. Island shapes can be viewed as a manifestation of the stress-driven elastic instability [13–17] for a thin film with a strong wetting effect [18]. Essentially, the formation of an island reduces the elastic energy of the film but increases the surface area of the exposed film surface, resulting in an increase in the surface energy. Competition between the elastic energy and surface energy sets a length scale for the island geometry, and the wetting of the film to the substrate results in compact, spatially localized island shapes. Mathematically, equilibrium shapes can be viewed as a free boundary problem arising from the minimization of a shape-dependent energy functional, in which the elastic energy contribution depends non-locally on the shape through the equations of elastic equilibrium. An example of such shapes for a two-dimensional (2D) model (i.e. an island ridge geometry) with isotropic surface energy is given in [18], showing that small island shapes are governed by a linear integro-differential equation, with numerically determined island shapes at larger volume approaching the shape of a ball sitting atop the substrate. A related body of work in the literature [19–29] has focused on determining how the island shape is affected by other effects such as anisotropy of surface energy, non-ridge (three-dimensional, 3D) geometries, non-uniform film composition and film substrate wetting as discussed below. Rigorous existence and regularity results for the free boundary problem appear in [30].

A key difficulty in determining equilibrium island shapes theoretically is the nonlinear, non-local character of the mathematical problem due to the coupling between elasticity and the island shape. As such, nearly all analytical results involve some type of approximation to solve the elasticity problem for an island of a given shape to determine the elastic energy. For example, the small-volume results in [18] use a small-slope approximation for the island shape to find a two-term asymptotic approximation for the elastic energy as described in [19]. The two-term approximation is also used in [20,21,23,24], while a third-order elastic energy approximation was developed for smooth surfaces in 2D in [31], for axisymmetric nanoring structures in [29], and for smooth surfaces in 3D in [27,28]. Other approaches involve numerical solution of the elasticity problem using a boundary integral method for the 2D geometry [18,22] or an analytic fitting model of the elastic relaxation function for the 3D axisymmetric geometry [25,26]. The results we detail here develop the third-order elastic approximation for faceted 2D island geometries. While third-order results have appeared previously for smooth surfaces [27,28,31], we have carefully extended the results to faceted surfaces which have corners.

An interesting feature of strained island morphologies is that the shape is typically size-dependent. This is in contrast to the classic behaviour of a liquid drop on a planar substrate for negligible gravity, for which the drop shape is independent of its size. Broadly speaking, small islands are short and wide, while large islands have a taller aspect ratio. The precise shape of the island is largely determined by the anisotropy of the surface energy. In the limiting case of isotropic surface energy, all surface orientations have the same surface energy and the resulting island shapes are smooth and progress from a sinusoidal-like bump at small volume to a ball-like shape at large volume [18]. Similar analyses for the shape at small volume for the isotropic case have been published for the axisymmetric island [23], and including compositional stress effects due to non-uniform composition within an alloy axisymmetric island [24]. At a strong level of anisotropy, some ranges of orientations are energetically forbidden, leading to the presence of corners. The resulting shape progression starts with shallow smooth islands as in the isotropic case, but at larger volume there is a discontinuous shape transition to an island with corners (see [22] for theory using a regularization term for the corner).

The case of extreme anisotropy corresponds to a faceted surface with orientations permitted at only a number of discrete orientations. There have been a number of studies which have investigated the energy minimizing shape of such piecewise-linear island shapes. Tersoff & Tromp [19] considered the shape transition of a faceted rectangular pyramid island, and they found that the symmetric island transforms to an elongated ‘quantum wire’ during island growth. Daruka and colleagues [20,21] used a 2D model analogous to that in [19] to study the energy of 2D faceted islands as a function of the shape, facet angles and facet energies to obtain a phase diagram for the island type as a function of the volume and surface energy parameters. A similar but more general analysis was completed by Gill & Cocks [25,26], who constructed a 3D model for the energetics of axisymmetric faceted islands, which is valid for high-aspect ratio islands and allows different elastic properties of the film and substrate. They determined similar equilibrium phase diagrams as in [20] but for 3D faceted axisymmetric cone islands, and interpret the results in terms of 3D faceted islands in [26].

The general picture that emerges from these studies of faceted equilibrium shapes is that the small islands are composed of small-angle facets and thus have relatively small aspect ratio, while as the island becomes larger its equilibrium shape incorporates more facets with larger angle and has larger aspect ratio. But what is not clear from these equilibrium theories is *how* new facets are created on the island as it grows. Put another way, there is an apparent first-order phase transition from one island shape to another at a critical island size. There are dynamical simulations of thin-film growth which display island formation [27,28] using a 3D continuum model with a third-order elastic approximation and weakly anisotropic surface energy giving rise to pyramid-shaped islands, but the surface energy model is not a strictly faceted model, and it does not allow for the possibility of higher angle orientations on the island, so is not currently viable for analysing the transition from pyramids to dome islands.

Recently, in [32], a 2D continuum model with a second-order elastic energy approximation has been used to determine the energy of a more general class of faceted island shapes which include both symmetric and asymmetric island shapes. These results show that the transition from a pyramid to a dome (which involves the creation of a new steep facet on each side of the island) generally occurs through a *metastable asymmetric half-dome state* which has a steep facet on one side only, and there is a critical facet nucleus and associated energy barrier for passing into the metastable state that can be found as a saddle point in a parametrization of the energy surface in terms of the island shape.

In this article, we use a 2D model similar to [32] but use a third-order approximation to the elastic energy to study the faceted island shape transition. In particular, the higher order elastic energy approximation permits us to determine not only the characteristics of the transition from pyramids to domes, but also the transformation from domes to multifaceted domes. As mentioned above, this third-order elastic energy approximation is a generalization of the smooth-surface results in [27,28,31] to explicitly allow for corners on the island shape. In addition to permitting us to describe the transition to multifacet domes, the higher order elastic energy approximation also gives an improved description of the pyramid to dome transition which gives an important qualitative difference in the details of this transition at the critical transition volume. A general summary of our results is that the pyramid-to-dome and dome-to-multifacet dome transitions have the same behaviour and a strikingly similar structure for the shape bifurcation diagram in the vicinity of the transition. During both of these two transitions, there exist metastable asymmetric shapes denoted as ‘half dome’ and ‘half multifaceted dome’ which include a ‘steep facet’ on one side only. Then by considering the effect of corner energy on the energies associated with these shape transitions, we conclude that *the lowest barrier transition pathway must necessarily pass through such asymmetric metastable shapes.* Finally, by constructing the transition pathways and calculating the corresponding energy barriers, we show that half-multifaceted domes are more likely to be observed in experiments compared to half domes, which agrees with the experimental results in [6].

The rest of the article is organized as follows. In §2, we construct our isotropic linear elastic model and find the third-order approximation to the elastic energy by analytically solving this elasticity problem. In §3, we explain our numerical methods to calculate the energies of islands for all shapes as a function of island volume. Then we construct the transition pathways between island types and determine the corresponding energy barriers, and, in addition, compare our results with those of [32] obtained using a second-order approximation to the elastic energy. In §4, we discuss the corner energy effect and compare our numerical results with experimental results and discuss other aspects of our theoretical model. In §5, we summarize our findings.

## 2. Mathematical model

### (a) Isotropic linear elastic model

In this article, we consider a 2D, isolated, single-component, completely faceted, epitaxially strained island on a substrate. Denoting each of the *n* facets on an island by index *i*, we take the facet orientation angle *θ*_{i} as a multiple of *θ*_{0}=11.2°; i.e. *θ*_{i}=*mθ*_{0} for some *m*=0,±1,±2,±3, where *θ* is measured counterclockwise from the vertical axis to the exterior facet normal. This choice of facet orientations is a 2D analogue for the {105}(*m*=±1), {113}(*m*=±2) and {15 3 23}(*m*=±3) facets on SiGe islands grown on Si(001). Thus, we can represent the facet orientations of an island as a vector of integers with *i*th element corresponding to the orientation *m* of the *i*th facet. Moreover, we can classify different island types that have differently oriented facets by the orientation vectors. For example, a pyramid corresponds to (1,0,−1), a dome corresponds to (1,2,1,0,−1,−2,−1) and a multifaceted dome corresponds to (1,2,3,2,1,0,−1,−2,−3,−2,−1). In figure 1, we illustrate different island types such as pyramid, dome and multifaceted dome. Note that pointed-top pyramids with facets (−1,1) are a special case of our pyramids corresponding to a top facet of length zero; we find this special shape is always energetically unfavourable in comparison to the flat-top pyramid.

We describe the island shape by *y*=*h*(*x*), where *h*(*x*) is a continuous piecewise-linear function defined in a finite domain. Note that in the Stranski–Krastanov growth mode, there is a wetting layer that covers the substrate before islands form. Here we use the glued wetting layer model [33] in which the wetting layer is assumed to be infinitesimally thin and glued to the substrate. Thus, we can extend *h*(*x*) to the entire domain to describe the film surface. In particular, *h*(*x*)>0 for the island and *h*(*x*)=0 for the wetting layer outside the island. Now *h*(*x*) is a continuous piecewise-linear function in the domain

The energy of the film-substrate system consists of two contributions, *E*_{tot}=*E*_{el}+*E*_{surf}. The surface energy is
*γ*(*θ*) is the strongly anisotropic surface energy density, which is dependent on the surface orientation *θ*. Since our 2D islands only have facets in several specific orientations, for simplicity (as in [32]), we assume that these facets have the same surface energy density, i.e. *γ*(*θ*_{i})=*γ* for *i*=1,2,…,*n*. We discuss the generalizations of the model to different facet energies in §4c.

The elastic energy is
*S*(*x*,*y*) is the elastic energy density.

We use isotropic linear elasticity to determine the elastic energy for the 2D system. We use a similar formulation to that in [18] with an asymptotic expansion approach similar to that in [31]. Defining *σ*^{F} and *σ*^{S} as the stress tensors in the film and substrate, respectively, the stresses satisfy the equilibrium equations
*y* direction. Moreover, the stress decays to zero at infinity in the substrate,
*e*_{0}=(*a*_{s}−*a*_{f})/*a*_{s} due to the different lattice parameters *a*_{f} in the film and *a*_{s} in the substrate. Therefore, a jump condition in the strain occurs at the interface of the film and the substrate
*e*^{F} and *e*^{S} are the strain tensors in the film and substrate, respectively, and
*ν* is the Poisson ratio. For isotropic linear elasticity, we have the constitutive equations to relate stress to strain in each of the film and substrate
*i*,*j*,*k* range over *x*,*y*, repeated indices imply summation, *E* is the Young modulus, *δ*_{ij} is the Kronecker delta, and where *E* and *ν*, to be the same in both film and substrate. This assumption lets one reduce the elasticity problem for the film/substrate system to an equivalent problem in a single semi-infinite domain; it is also a reasonable simplifying assumption because the elastic constants for Si and Ge are similar. We discuss general elastic constants in §4c. Finally, the elastic energy density in each of the film and substrate is

Considering a planar film *y*=*h*_{0} on the substrate, where *h*_{0}≥0 is a constant, the elasticity problem has the exact solution *σ*_{0}=(*E*/(1−*ν*^{2}))*e*_{0} is the corresponding stress due to the misfit strain *e*_{0}. The elastic energy density for the planar film is

Now for a non-planar film *y*=*h*(*x*), we define new variables for the stress tensor ** σ** and strain tensor

**for the entire domain such that**

*e**y*=0 in (2.6) and (2.8) reduce to the continuity of both stress and strain. Then

**satisfies**

*σ**y*=

*h*(

*x*) obtained from (2.5) is

**and**

*σ***reduce to a single-material elasticity problem for the semi-infinite domain**

*e**y*<

*h*(

*x*).

### (b) Perturbation solution and energies

To solve this reduced elasticity problem, we use a perturbation method based on a thin-film and small-slope approximation where the surface height *h*(*x*) is assumed to be relatively small, i.e. *h*(*x*) is of order *ε*, where *ε*≪1 . Expand the stress in powers of *ε*
*k*=1,2,…. Substituting the perturbation expansion of *σ*_{ij} into equations (2.16)–(2.19) and expanding each order term *h* in (2.17) and (2.18), we have at *O*(*ε*)
*O*(*ε*^{2})

Now we have reduced the elasticity problem in the solid *y*<*h*(*x*) to a sequence of problems for an elastic half-space with known surface tractions. In particular, both problems (2.21) and (2.22) are of the general form

Using these solutions for equation (2.21), we directly obtain the solutions at *O*(*ϵ*). Then we can determine the boundary conditions in (2.22). For example, the right side of the third line of equation (2.22) contains the derivative of *G*(*ξ*;*x*,*y*) satisfies the following properties:
*δ*(*x*) is the Dirac delta function. So by taking the limit of equation (2.27), we get the required boundary term appearing in (2.22)
*y*=0 up to *O*(*ϵ*^{2}) are
*H*(*f*) is the Hilbert transform of *f*(*x*), which is defined using the Cauchy principal value, denoted by p.v.,

With this perturbation solution, we determine the total elastic energy in the film and substrate, *S*(*x*,*y*) valid for small *h*(*x*)=*O*(*ϵ*)
*S*(*x*,*y*) obtained from equation (2.11) in the last step. In the substrate, we use the divergence theorem

where ∂*S* consists of the *x*-axis, denoted *C*_{1}, and an arc at infinity, denoted *C*_{2}. We use the relation between the strain and displacement
_{j}*σ*_{ij}=0 for *i*=*x*,*y* in the last step.

Considering the integral on *C*_{1}, we use the perturbation expansion of displacements and stresses and obtain
*y*=0. To obtain *C*_{1}

For the integral on *C*_{2}, we claim that *π*≤*θ*≤0, and get
*k*,*l*=1,2, and *G*_{ij}(*x*,*y*,*ξ*),*F*_{ij}(*x*,*y*,*ξ*) are the kernel functions and *f*^{(l)},*g*^{(l)} are the boundary conditions for *O*(*ε*^{l}). The kernel functions have the general form
*a*,*b*=0,1,2,3 such that *a*+*b*=3. So we get *G*_{ij}(*x*,*y*,*ξ*), *F*_{ij}(*x*,*y*,*ξ*)=*O*(1/*r*) independent of *ξ*. Combining with the fact that our boundary conditions *f*^{(l)},*g*^{(l)} (*l*=1,2) are integrable functions with compact support, we obtain *i*) for *x*_{j}=*x*,*y*. Using the same strategy with equation (2.42), we can show that
*C*(*r*) is some function of *r*. Moreover, the boundary conditions *f*^{(l)},*g*^{(l)} (*l*=1,2) satisfy

With the results of *k*,*l*=1,2. Thus, we complete our proof for

Combining the results for *E*^{F}_{el} and

Our results for the remainder of the paper are shown in terms of characteristic scales for length *l*_{0}=*γ*/*S*_{0}, volume *E*_{0}=*γ*^{2}/*S*_{0}. Using the parameters for SiGe films (for Si: *a*=5.4310×10^{−8} cm, *E*=13.0×10^{11} erg cm^{−3}, *ν*=0.278; for Ge: *a*=5.6575×10^{−8} cm, *γ*=1927 erg cm^{−2}) these characteristic scales are *l*_{0}∼20 nm, *E*_{0}∼3×10^{−3} erg cm^{−1} ∼20 eV Å^{−1}.

## 3. Numerical calculations and results

### (a) Numerical method

We find the minimum energy shape *h*(*x*) of an island with a fixed volume

Given an island type, as a sequence of facets with prescribed orientations, the shape function *h*(*x*) of the island depends only on the lengths of the facets since the number of facets and the slope of each facet are already determined. Suppose an island has *n* facets and *L*_{i} is the length for the *i*th facet, then we need to minimize *x*_{0} and *x*_{n} are the left and right edges of the island, respectively, and the height constraint guarantees that the island is continuous with the wetting layer. Note equation (3.3) gives the domain of the variables: if *L*_{i}<0 for some *i*, *h*(*x*) is a non-physical shape that is out of our consideration; if *L*_{i}=0 for some *i*, *h*(*x*) corresponds to another island type which has fewer facets. Since the objective function represents the energy of an island, the constrained optimization problem is expected to be well-behaved in the domain: it is smooth and has at most one local minimum away from the boundary (all *L*_{i}>0); local minima at the boundary (*L*_{i}=0 for some *i*) can also exist but correspond to the minimum energy shapes for other island types that have fewer facets. In other words, we have only one minimum energy shape, if exists, for each island type with a fixed volume. Note that, for a given island type with *n* facets, the minimum away from the boundary is a stationary point in the *n*-dimensional space while a minimum at the boundary may not be a stationary point in this *n*-dimensional space but corresponds to the stationary minimum in a lower dimensional space.

Using the first two constraints, we can express two of our variables in terms of the other *n*−2 variables so that we reduce the dimension of the variables to *n*−2. For symmetric islands, we can further reduce the dimension of the variables to (*n*−1)/2. To explicitly avoid non-physical shapes with *L*_{i}<0, we impose a large energy penalty to such non-physical shapes. In this way, the *n*-dimensional constrained minimization problem is reduced to a (*n*−2)-dimensional (or (*n*−1)/2 for symmetric islands) unconstrained minimization problem. To solve this minimization problem, we use a quasi-Newton functional minimization method [35] to find the possible stationary local minimum away from the boundary. The quasi-Newton functional minimization method will not converge to the non-stationary minima at the boundary. For minima close to the boundary, to avoid numerical difficulties when *L*_{i}<0 due to the discontinuity, we use the Nelder–Mead method instead of the quasi-Newton method. While we still have the possibility of nearby local minima at the boundary (*L*_{i}=0 for some *i*), in general, we can choose the initial guess away from the boundary to find any interior minima close to the boundary. For both the quasi-Newton method and the Nelder–Mead method, we choose a convergence criteria such that the iteration stops if the relative error in objective function is smaller than tolerance 10^{−8} and for the quasi-Newton method, if the gradient is smaller than tolerance 10^{−6}.

In addition to finding the local minima of the energy, we also analyse the energy surface to determine saddle points and construct transition pathways. Since we care about how the new steeper facets are formed during each stage of the transition process, we analyse the energy surface in terms of the lengths of two new facets, denoted as *W*_{1} and *W*_{2}, which enables consideration of asymmetric islands. Given an asymmetric island type and a fixed volume, we find the minimum energy shape for each pair (*W*_{1},*W*_{2}) by the same method discussed above except that here we fix the lengths of two new facets as (*W*_{1},*W*_{2}) and have only *n*−4 variables.

### (b) Numerical results

By examining all possible island types as a function of volume, we find that the pyramid is stable for small volume, dome is stable for intermediate volume and multifaceted dome is stable for large volume. It is important to note that the stable pyramid, dome and multifaceted dome, denoted as P, D_{1} and D_{m} respectively, are symmetric and that the orientations of adjacent facets of these three island types differ by *θ*_{0}.

Considering the pyramid–dome transition process, we suppose that an island at a fixed volume is changing from pyramid (P) to dome (D_{1}) by nucleating a pair of steep facets, one on each side of the island. To understand how the island forms two new facets, we plot in figure 2 the energy surface as a function of the lengths of two steep facets for some specific volumes. In figure 2, when both of *W*_{1} and *W*_{2} are zero, it corresponds to a pyramid; when one of *W*_{1} and *W*_{2} is zero, it corresponds to a half dome with only one steep facet; when *W*_{1} and *W*_{2} are equal and non-zero, it corresponds to a symmetric dome.

Figure 2 shows how the energy surface changes as island volume increases: at first, in figure 2*a* pyramid P is the only energy minimum and it is stable; then at larger volume the local minimum D_{1} and local maximum *b*D_{1} appear simultaneously in figure 2*b*; following, two asymmetric saddle points *b*A_{1} split from *b*D_{1} in figure 2*c*; next, boundary minimum H_{1} and boundary maximum bH_{1} for half dome appear in a pair in figure 2*d*; afterwards, saddles bA_{1} move to the boundary and merge with H_{1} in figure 2*e*; at last, for the largest volume, bH_{1} and uD_{1} move towards P and eventually merge with P when P becomes the maximum of the energy surface in figure 2*f*.

With understanding of the energy surfaces, we construct the lowest barrier transition pathways for each volume as marked by the arrows in figure 2*a*–*f*. For example in figure 2*d*, the island is at the volume where pyramid (P) and dome (D_{1}) have the same energy, which is the critical volume for the pyramid to dome transition. The island nucleates one steep facet first (bH_{1}) after which this facet grows and the island reaches a metastable state corresponding to a half dome (H_{1}). Then the island nucleates the second steep facet on the other side (*b*A_{1}) and finally becomes a stable dome (D_{1}). Through this transition pathway, the island needs to climb over two energy barriers at bH_{1} and bA_{1}, each for nucleating one steep facet. Note this transition pathway involving sequential nucleation of two steep facets has smaller first barrier P–bH_{1} compared to the pathway where two steep facets are nucleated simultaneously. Moreover, the transition pathway for nucleating two steep facets simultaneously, the symmetric pathway P–uD_{1}–D_{1}, has the largest energy barrier P–uD_{1} of any pathway from P to D_{1}.

More generally, we plot the energy versus volume curves for the stable, metastable, unstable and barrier states of different equilibrium island types for the pyramid-dome transition in figure 3. The energy surfaces and transition pathways in figure 2*a*–*f* are the six distinct cases associated with the volumes labelled at the top of figure 3. In particular, figure 2*d* corresponds to the volume where P and D_{1} have equal energy; figure 2*e* is at the volume where asymmetric barrier points bA_{1} merge with H_{1} and figure 2*f* indicates the volume where P becomes unstable. We can see in figure 3 that for islands at the volume range between *d* and *e* labelled at the top of figure 3, the energy barrier for the symmetric transition pathway *E*(uD_{1})–*E*(P) is much larger than both the energy barriers for the asymmetric transition, *E*(bH_{1})–*E*(P) and *E*(bA_{1})–*E*(H_{1}), so the island prefers the asymmetric transition pathway P–bH_{1}–H_{1}–bA_{1}–D_{1} as depicted by ① in figure 1. For islands at the middle volume range between *e* and *f*, an asymmetric pathway is also preferred but H_{1} is not a metastable state anymore. Since the saddle bA_{1} merges with H_{1} (figure 2*e*), there is no second energy barrier H_{1}–bA_{1} to climb. Therefore after bH_{1}, the transition path will leave the boundary in figure 2*e* before H_{1} when the steepest descent gradient points into interior. In other words, the island will not get to H_{1} but will nucleate the second steep facet before H_{1} and then end at D_{1} as path ② in figure 1. For islands at volume larger than *f*, pyramid (P) becomes an unstable state so that the island transforms to D_{1} through a monotonic decrease in energy through symmetric shapes corresponding to path ③ in figure 1.

The symmetry-breaking in the pyramid–dome transition is also demonstrated by Spencer & Tersoff [32]. They used the small-slope approximation up to second-order for the elastic energy and obtained a similar topography for the energy surface as shown in figure 2. A qualitative difference, however, is that they have the energy surface structure in figure 2*c* instead of that in figure 2*d* at the P–D_{1} critical transition volume, and hence have one extra transition pathway P–bA_{1}–D_{1} as marked by the arrows in figure 2*c*. For this path, the island has both steep facets nucleated simultaneously but with different sizes, which is an unusual phenomenon. The reason why we do not have this unusual transition path is that here we are using the small-slope elastic energy approximation up to third-order, which is more accurate than the second-order approximation in [32]. As a result, each curve in figure 3 is actually shifted relative to their bifurcation diagram so that the bifurcation structure is slightly changed. In figure 3, the curves for H_{1} and *bH*_{1} start at the volume smaller than the critical transition volume, while in [32] the starting point for these two curves is on the right side of the intersection of the curves for P and D_{1}. Therefore, we have the path P–bH_{1}–H_{1}–bA_{1}–D_{1} in figure 2*d* at the critical transition volume and the path P–bA_{1}–D_{1} in figure 2*c* exists only for subcritical volumes for which the P–D_{1} transition is not energetically favourable.

Because we have a third-order elastic energy approximation, we can use our approach to investigate the next island shape, multifaceted dome, at an even larger volume and describe the island shape transitions from dome to multifaceted dome that are not considered in [32]. For the dome to multifaceted dome transition process, the island will nucleate two more steep facets with orientation angles ±3*θ*_{0}. The topography of the energy surface in terms of the lengths of these two steeper facets has a similar structure to that for the pyramid–dome transition. In addition, the bifurcation diagram of the equilibrium states of dome and multifaceted dome, such as D_{1}, H_{m}, D_{m}, uD_{m}, bA_{m} and bH_{m}, has the same structure as in figure 3. As a consequence, the dome to multifaceted dome transition also has three different transition pathways for islands at different volume ranges (similar to paths ①–③ for pyramid to dome), which are represented by paths ④, ⑤ and ⑥ in figure 1.

In figure 4, we plot the energy curves of all stable and metastable island types to summarize all shape transitions during the growth of islands from pyramid to multifaceted dome. We can see that pyramid (P) is stable for small islands, dome (D_{1}) is stable for medium-sized islands and multifaceted dome (D_{m}) is stable for large islands. So the pyramid–dome transition takes place when small islands grow to medium size and dome to multifaceted dome transition happens when medium-sized islands become larger. We label the volume ranges with numbers ①–⑥ at the top of figure 4 that correspond to the transition pathways ①–⑥ in figure 1. Islands in volume ranges ① or ④ will have asymmetric transition pathways ① or ④ (figure 1) that pass through the metastable states half dome (H_{1}) or half multifaceted dome (H_{m}), respectively. For islands in volume ranges ② or ⑤, half island types are unstable states so that the pathways ② and ⑤ (figure 1) do not pass through half island types but they are asymmetric because of the lower energy barriers for nucleating the first facet. For islands in volume ranges ③ or ⑥, pyramid (P) or dome (D_{1}) are unstable states, so islands directly transform to stable dome/multifaceted dome states through the symmetric paths ③ or ⑥ (figure 1), respectively.

## 4. Discussion

### (a) Corner energy effects

The numerical results and conclusions above are based on the approximate energy

To estimate the magnitude of the edge/corner energy effect, we consider the corner energy regularization model discussed in [36–40] for the energy of a rounded corner between two facets. For a round corner region with size comparable to the lattice parameter *a*, we have
*κ*∝1/*a* is the curvature in the corner region, *β* is a corner energy parameter and Δ*θ* is the angle between facets. A typical magnitude for *β* can be determined from considering the extra surface energy at a right-angle corner due to broken bonds
*β*=(2/*π*)*E*_{bond}. Using the parameters for SiGe films from §2 and *E*_{bond}=42.6 kcal mole^{−1}, we can obtain our scaled corner energy between facets of orientation angle difference *θ*_{0}=11.2° as *E*_{c}=1.1191×10^{−3}*E*_{0} (≈20 meV Å^{−1}). Our estimate obtained by counting the excess bond energy is consistent with recent atomistic calculations of corner energy for pyramids in Ge/Si by Retford *et al.* [41]. They gave an average value for the (001)/(105) corner as 10 meV Å^{−1}, which is within a factor of two compared to our estimate.

We consider first how the pyramid to dome transition in figure 3 is modified by corner energy. We show results using a value for the corner energy as estimated above; different values of the corner energy result in the same qualitative behaviour. As shown in figure 5, adding corner energy will result in the corresponding upward shift for each energy curve relative to pyramid P (which we regard as the relative standard). Since the dome has four more corners relative to the pyramid (two corners on each side), the energy should be increased by an amount of 4*E*_{c} due to these four corners, i.e. *E*(cD_{1})=*E*(D_{1})+4*E*_{c}. Similarly for the half dome, *E*(*cH*_{1})=*E*(H_{1})+2*E*_{c} due to two more corners on one side. We can summarize the changes of figure 3 due to corner energy: the energy curves for D_{1}, uD_{1} and bA_{1} will be shifted upward by 4*E*_{c} and those for H_{1} and bH_{1} will be shifted upward by 2*E*_{c}. Although each curve is shifted upward by an amount corresponding to corner energy, the bifurcation structure of the solutions in the small volume range (corresponding to the original volume range between *d* and *e* in figure 3) does not change except that the critical transition volume is shifted to the right. Thus, islands in the small volume range still have the lowest barrier transition pathway P–bH_{1}–H_{1}–bA_{1}–D_{1}. But now both of the two energy barriers P–bH_{1} and H_{1}–bA_{1} increase by 2*E*_{c} due to the corner energy. In the middle volume range (corresponding to the original volume range between *e* and *f* in figure 3), since bA_{1} (not shown in figure 5) is shifted by 4*E*_{c} while H_{1} is shifted by 2*E*_{c}, the curves for bA_{1} and H_{1} will not merge anymore, but instead they will end with separation 2*E*_{c}. Similar results apply to the energy versus volume curves of P, uD_{1} and bH_{1} in the large volume range, where instead of merging, they will terminate with uniform energy difference such that *E*(uD_{1})=*E*(bH_{1})+2*E*_{c} and *E*(bH_{1})=*E*(*P*)+2*E*_{c}. In other words, the island has to conquer the corner energy barrier to form each new facet, in addition to the original energy barrier. While the presence of corner energy does not affect the bifurcation structure at small volume, its presence does result in a qualitative change in the transition sequence at larger volume. Recall from figures 2*d*–*f* and 3 that the half-dome state H_{1} could be unstable at larger volume, thus permitting transition from pyramid to dome that did not pass through a metastable half-dome state. Now, however, due to the presence of corner energy, we observe that for all volumes the half-dome state H_{1} is metastable, and correspondingly *the lowest-barrier transition from pyramid to dome must always pass through the metastable half-dome state*. To see this, considering the case in figure 2*e*, even though there was no original energy barrier from a half dome H_{1} to an asymmetric dome A_{1} without corner energy, when the corner energy is included, there is a corner energy barrier for H_{1}–A_{1} so that half dome H_{1} is now a metastable state. So the medium-sized island must also pass through H_{1} to reach D_{1}. For the case in figure 2*f*, with corner energy, there is a corner energy barrier 4*E*_{c} for the symmetric pathway P–D_{1}, while there are two corner energy barriers both equal to 2*E*_{c}, one for P–H_{1} and the other for H_{1}–A_{1}, for the asymmetric pathway. Based on the transition-state theory, even though the two transition pathways require the same total energy barrier, the asymmetric pathway is still preferred since each of its energy barriers is smaller than that for the symmetric pathway. As a consequence, in the entire volume range for pyramid to dome transition, half dome (H_{1}) is always a metastable state so that we have unique transition pathway P–bH_{1}–H_{1}–bA_{1}–D_{1} for all-sized islands. Similarly, for the dome to multifaceted dome transition with corner energy, the transition is shifted to larger volume and we have only one transition pathway D_{1}–bH_{m}–H_{m}–bA_{m}–D_{m}.

In addition to determining how corner energy affects the existing bifurcation structure, we also examine the influence of corner energy on the stability of pyramid, dome and multifaceted dome relative to other island types that have the same steepest facet orientation but fewer corners. We find that pyramid, dome and multifaceted dome are all still stable island types and no other island types become energetically favourable due to the corner energy. For example, a steep-faceted dome with orientation vector (1,3,1,0,−1,−3,−1) has four corners on each side, two with angle *θ*_{0} and two with angle 2*θ*_{0}; a multifaceted dome with the orientation vector (1,2,3,2,1,0,−1,−2,−3,−2,−1) has six corners with *θ*_{0} on each side. Then the total corner energy of a steep-faceted dome is 12*E*_{c}, where each large corner with angle 2*θ* has corner energy 2*E*_{c}; the total corner energy of a multifaceted is 12*E*_{c}, *E*_{c} for each corner. So steep-faceted domes are still unfavourable even though they have fewer corners than multifaceted domes do.

### (b) Comparing results with experiments

#### (i) Island edge details

The island shapes in figure 1 are based on our numerical calculations so we can explicitly determine the size and position of each facet during the transition. Initially, the new steeper facets are always formed near the middle of the original steep facet and then grow to extend nearly to the edge of the island, but not all the way to the edge of the island. It is energetically preferred that shallow contact facets are present at the base as the ‘feet’ of the island although they are very small relative to other facets [21,32] (figure 1). Note Sutter *et al.* [4] observed in experiment that a shallow {105} contact facet was present at the base of Si_{1−x}Ge_{x} dome, which connected the steep facet with the substrate. In addition, they even found a deep depression at the base of the barn-shaped island to provide island–substrate contact via {105} facets. This deep depression observed in experiment means that the wetting layer had disintegrated due to the effect of alloyed material, high stress at the base of the barn shape, and/or the presence of neighbouring islands [4]. These experimental results demonstrate that an island base configuration is more complicated than the intersection of a planar substrate with a steep facet. Although our 2D theoretical model assumes a single-material island resting on a planar substrate without any depression, our results also show that it is energetically favourable that a shallow contact facet appears at the island base to connect the island with the planar substrate instead of the steep facet directly contacting the substrate, similar to the results in [21,32].

#### (ii) Asymmetric islands

Based on our 2D model calculations, we find that metastable asymmetric shapes exist during both the pyramid to dome transition and the dome to multifaceted dome transition. In volume ranges ① and ④ in figure 4, the half dome H_{1} and the half multifaceted dome H_{m} are the metastable states. Considering the volume range ①, there are two energy barriers during the pyramid to dome transition, and we note that the first energy barrier P–bH_{1} is much larger than the second energy barrier H_{1}–bA_{1} as shown in figure 6. Thus, for finite fluctuations, an island transforms from H_{1} to bA_{1} much faster than it does from P to bH_{1}, which makes the lifespan of the metastable state H_{1} relatively short so that it is difficult to observe the half dome in experiment. Similarly, for the dome to multifaceted dome transition in the volume range ④, the first energy barrier D_{1}–bH_{m} is also somewhat larger than the second one H_{m}–bA_{m} in figure 6, again suggesting that the metastable state is somewhat short-lived.

We now consider these findings of metastable asymmetric states in our 2D model in view of experiments on 3D domes. In particular, Ross *et al.* [6] do observe asymmetric islands in experiment. In their experiments, the islands are found to evolve from pyramid to a ‘transitional pyramid’, which has the {113} facets present at the four lowest vertices of the pyramid shape. Then the islands transform through sequential nucleation of four pairs of the {15 3 23} facets to become a complete dome. During the transition, asymmetric dome shapes with only one or three pairs of {15 3 23} facets are observed, but they do not observe asymmetric transitional pyramids in their experiments. We suggest that this difference might be explained on the basis of the nucleation barriers found in our 2D model. In 2D, the pyramid to dome transition corresponds to the development of {113} facets and the dome to multifaceted dome transition corresponds to the development of {15 3 23} facets. Although both the asymmetric metastable states H_{1} and H_{m} are short-lived relative to the stable states, the half multifaceted dome H_{m} is relatively stable compared with the half dome H_{1} because the difference between the energy barriers D_{1}–bH_{m} and H_{m}–bA_{m} is significantly smaller than the difference between the energy barriers P–bH_{1} and H_{1}–bA_{1} as shown in figure 6. Therefore in our 2D model, half-multifaceted dome shapes H_{m} are more likely to be observed than half-dome shapes H_{1}. This result is consistent with the 3D experimental observation in [6] that asymmetric domes are observed during the transition rather than asymmetric transitional pyramids.

### (c) Model discussion

Our calculations are based on a third-order approximation to the elastic energy. While we do not have a rigorous error bound on this approximation, Gill [29] uses an analogous approximation for a 3D conical island and determines the error in the elastic energy as a function of the surface slope. At slope values corresponding to the maximum slope on our pyramid, dome and multidome these errors are less than 1%, 2% and 13%, respectively, suggesting the third-order approximation is quantitatively accurate for pyramids and domes, and reasonably accurate for multifacet domes.

The calculations here are based on our reduced 2D model. We assume an isolated, fully faceted 2D island resting on a planar substrate. Our results are based on the third-order approximation to the elastic energy in 2D and the elastic energy is different from that in 3D, where in 3D we allow the strain in the third direction. While our 2D model obviously does not give a full description of the 3D shape transitions involving the {105}, {113} and {15 3 23} facets, the model does show that each facet nucleation event seems to have a generic bifurcation structure, and moreover that the differences in energy barriers associated with the formation of pairs of steeper facets are correlated to the observations of steep faceted asymmetric states in 3D islands.

In our model, we use the simplifying assumption that different facets have same surface energy because we are considering mainly the non-trivial effect of elastic energy in promoting asymmetric shape transitions. While it is true that different facets can have different surface energies, these differences are not very large. Small differences in surface energies should result in small changes to the bifurcation diagram, so the general sequence of asymmetric transitions is expected to be robust with respect to different surface energies for the different facets. In future work, one could allow for different facet surface energies and investigate the parameter space of the different facet energy combinations (see for example the equilibrium shape diagrams for quantum dots in [20] (2D symmetric) and [26] (3D axisymmetric)). Additionally, one could also consider the effect of surface stress/strain on modifying the effective surface energy of each facet [42]. Our model could also be generalized to incorporate different elastic constants in the film and substrate as in [43] (2D) and [27] (3D), and the effect of neighbouring islands on elastic energy, either by considering interactions of island pairs or the growth of an array of strained islands. It is left for future work to understand how the above components affect the asymmetric transition of quantum dots.

## 5. Summary

We construct a 2D continuum model to describe the energetics of shape transitions in fully faceted epitaxial islands via minimization of elastic energy and surface energy at fixed volume. The elastic energy of the island is based on a third-order approximation, enabling us to consider shape transitions between pyramids, domes, multifaceted domes and asymmetric intermediate states. The energetics of the shape transitions are determined by numerically calculating the facet lengths that minimize the energy of a given island type for prescribed island volume. By comparing the energy of different island types with the same volume and analysing the energy surface as a function of the island shape parameters (allowing for explicitly asymmetric shapes), we determine the bifurcation diagram of equilibrium solutions and their stability, as well as the lowest barrier transition pathway for the island shape as a function of increasing volume.

Similar to the results in [32] (which used a second-order elastic energy approximation), we find that the pyramid to dome transition can occur through a metastable half-dome state. Here, however, because of the increased accuracy of the third-order elastic energy approximation, we find that the there is a qualitative difference in the bifurcation diagram: at the critical transition volume for the pyramid to dome transition the barrier state is the metastable half dome, whereas in [32] the second-order elastic energy gives the barrier state at the critical transition volume as an asymmetric dome. Thus, we find that sequential nucleation of facets occurs over a wider range of volumes and in particular occurs at the critical transition volume when the dome first becomes lower energy than the pyramid.

The third-order elastic energy approximation also lets us consider multifaceted dome shapes with steep facets which could not be considered in [32]. For the resulting bifurcation diagram, we find that the symmetric shapes of pyramids, domes and multifaceted domes are stable, with the asymmetric shapes of half dome and half multifaceted dome being metastable transition shapes. Moreover, we find that there is a self-similar bifurcation structure for the pyramid-to-dome and dome-to-multifaceted dome transitions.

By considering the effect of corner energy (facet edge energy) on the energetics of the shape transition, we determine how the bifurcation diagram and stability of solutions are altered. We find that while the critical volume for each shape transition is increased, the relative stability of the solutions remains the same with pyramid, dome and multifaceted dome as stable states. Shape transitions also still occur through asymmetric intermediate states, but an important finding is that the additional barrier of the corner energy means that *each shape transition between symmetric states must pass through a metastable one-sided (half dome, half multifaceted dome) state*. This is in contrast to the results without corner energy for which the one-sided state is unstable at large volume and in which case the shape transition avoids passing through a metastable state.

Finally, we relate the predictions of our 2D model to observations of 3D islands in experiment. While a precise comparison is not possible because of the different dimensionality, we suggest that the relative stability of asymmetric transition states in experiment is consistent with the relative heights of the energy barriers for the facet nucleation events in our 2D model. In particular, asymmetric {15 3 23} facet domes are observed in experiment but asymmetric {113} transitional pyramids are not [6].

## Data accessibility

The research in this paper did not involve datasets from experiments. The numerical codes are based directly on the algorithms and inputs as described in the text.

## Authors' contributions

The theoretical development was done by C.W. and B.J.S. The numerical codes were written by C.W. in consultation with B.J.S. The paper was written primarily by C.W. with assistance by B.J.S.

## Competing interests

We have no competing interests.

## Funding

Portions of this work were supported by a Simons Foundation Collaboration Grant for Mathematicians (B.J.S.).

- Received April 11, 2016.
- Accepted May 26, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.