Continuum modelling of semiconductor heteroepitaxy: an applied perspective

Semiconductor heteroepitaxy involves a wealth of qualitatively diﬀerent, competing phenomena. Examples include three-dimensional island formation, injection of dislocations, mixing between ﬁlm and substrate atoms. Their relative importance depends on the speciﬁc growth conditions, giving rise to a very complex scenario. The need for an optimal control over heteroepitaxial ﬁlms and/or nanostructures is widespread: semicon- ductor epitaxy by molecular beam epitaxy or chemical vapour deposition is nowadays exploited also in industrial environments. Simulation models can be precious in limiting the parameter space to be sampled while aiming at ﬁlms/nanostructures with the desired properties. In order to be appealing (and useful) to an applied audience, such models must yield predictions directly comparable with experimental data. This implies matching typical time scales and sizes, while oﬀering a satisfactory description of the main physical driving forces. It is the aim of the presentreviewtoshowthatcontinuummodelsofsemiconductorheteroepitaxyevolvedsigniﬁcantly,providingapromisingtool(evenaworkingtool,insomecases)tocomplywiththeabove requirements. Several examples, spanning from the nanometre to the micron scale, are illustrated. Current limitations and future research directions are also discussed.


Introduction
Heteroepitaxy refers to any deposition process leading to a crystalline film composed of a given material, grown over a crystalline substrate made of another material. As such, this definition embraces several different systems. While most of the aspects discussed in this review are also general, we find it natural to mainly focus on semiconductor heteroepitaxy, a research field involving a huge academic and industrial community. Even more specifically, semiconductor heteroepitaxy on silicon seems to capture ever-growing attention. The possibility to improve both electronic and optical properties of the material while still exploiting mainstream silicon technology, indeed, attracted important investments as a new generation of high-performance devices ready for highvolume industrial production is envisaged. Reviews focusing on both present and future applications exploiting heteroepitaxy of Ge (or, SiGe alloys), GaAs, GaN can be found in Refs. [1][2][3][4][5][6][7][8][9][10].
Despite considerable progress, the possibility to finely tune the deposition conditions in order to obtain the desired film quality, composition, and morphology (both low-roughness 2D films and 3D structures can be appealing, depending on the application) still calls for improvements. The task is of formidable complexity, due to the amount of competing physical phenomena.
The lattice mismatch between film and substrate induces a biaxial stress (typically compressive, if deposition on Si is considered). Two alternative paths are available to release the excess elastic energy: plastic relaxation via the introduction of misfit dislocations or elastic relaxation via formation of 3D islands [11][12][13][14][15]. 3D growth can be preceded by the formation of a thin wetting layer (Stranski-Krastanow growth) depending both on the deposition conditions and on the (film-thickness dependent) difference between film and substrate surface energies [9,16,17].
If the misfit is low, the film tends to grow flat (provided that the interfacial energy is not too high) until dislocations are injected. This is the typical behaviour observed in Si 1−x Ge x alloys grown on Si(001), with x 0.2 [18,19]. The opposite behaviour is observed at high misfits where the appearance of 3D nanostructures precedes the onset of plastic relaxation, as observed for high-x Si 1−x Ge x /Si(001), or for InGaAs/GaAs(001).
Elastic and plastic relaxation often act simultaneously. Coherent islands tend to enlarge their height-to-base aspect ratio (AR) during growth, as this allows for a better strain relaxation [9]. Once a critical volume is reached, however, the introduction of misfit dislocations become more convenient from the energetic point of view [20]. At this stage, the AR vs. volume curve suddenly changes, almost assuming the form of a plateau [21,22]. A close inspection of the island shape reveals fascinating oscillations [23]. Interestingly, shape oscillations were also obtained in the absence of dislocations, as due to intermixing [24]. Actually, despite being favoured by entropy, intermixing is yet another channel leading to strain relaxation as alloying lowers the effective lattice mismatch between film and substrate [25,26]. Intermixing becomes of particular importance at high temperatures, when deposited atoms can easily exchange position with the substrate ones [27,28]. Possible consequences involve the formation of a wetting layer not uniquely composed of the deposited material, and/or 3D islands hosting a significant concentration of substrate atoms. When deposition of pure Ge on Si is considered, the effect can be dramatic. The enthalpy of mixing of Ge/Si is indeed negligible, so that entropy of mixing forces alloying [19,29]. At deposition temperatures around 700 • C or higher, 3D islands formed upon deposition of Ge on Si(001) by Molecular Beam Epitaxy (MBE) can contain less than 50% of Ge [30]. Intermixing takes place also in InGaAs/GaAs [31], but the effect is less dramatic due to the non-negligible enthalpy of mixing. It is important to point out that, at the typical temperatures of interest, intermixing occurs only at the surface (or, at the subsurface) [32][33][34], as bulk diffusion is kinetically hindered. This observation is key, as, for any miscible system, deposition of a thin film on a thick substrate would lead to full dilution of the former in the latter, as determined by entropy maximization.
What makes the study of heteroepitaxy at the same time challenging and fascinating is that the aforementioned processes influence each other, their relative importance depending on the deposition conditions. For instance, 3D islanding following Ge deposition on Si(001) can be suppressed by growing outof-equilibrium. High deposition rates, and/or low deposition temperatures (as an alternative to surfactants), indeed, lead to fast film thicknening without allowing for the formation of islands [35][36][37]. In such a metastable state, the system relaxes by injecting dislocations suppressing further islanding as the effective lattice mismatch is reduced. This is of paramount importance if one is willing to grow flat substrates, as needed for many applications. On the opposite side, the onset of plastic relaxation can be delayed by depositing on suitably patterned substrates [38][39][40][41].
To make things even more complex, strain also influences surface energies [42]. The most striking example is given by Ge{105} whose surface energy is stabilized by compression [43][44][45] to the point that the (001) wetting layer can spontaneously break into {105} faceting [7,46]. Micron-long, horizontal Ge wires delimited by {105} facets, have been recently observed [24,47] and shown to display peculiar electronic properties.
The aforementioned examples should make it clear that a subtle competition between elastic and plastic relaxation and intermixing is often determining, together with surface energy, the final properties of heteroepitaxial films and nanostructures. Given this complex scenario, it is clear that reliable models tackling semiconductor growth dynamics could greatly help in determining optimal deposition conditions leading to a film with the desired properties.
Atomistic approaches can only tackle small systems and short time scales. Worst, diffusion processes at semiconductor surfaces are typically influenced Downloaded by [187.160.112.235] at 22:00 10 June 2016 by local changes in orbital hybridization, so that reliable estimates of diffusion rates call for refined and computationally demanding ab initio calculations [48]. Nevertheless, it is worth mentioning some significant progress based on empirical approaches exploiting Kinetic Monte-Carlo schemes [49][50][51][52].
Continuum models are a possible solution to close both spatial and temporal gaps between simulations and experiments [53][54][55][56]. Based by definition on a coarse-grained description of the actual system, any continuum model of growth will necessarily rely on some simplifying assumptions. With this respect, only a final comparison with experiments can fully establish the validity of the specific approach. If successful, a continuum model becomes extremely appealing as it allows for a deep understanding of the main physical driving forces underlying some complex evolution (much harder to be inferred, e.g. from the observation of billions of atomic trajectories).
It is the aim of the present work to offer a survey of the present status of continuum models tackling semiconductor heteroepitaxy. Importantly, authors wish to share an applied perspective. This means that attention will be focused only on the deposition conditions commonly employed when growing highquality semiconductor films, nanostructures, or heterostructures, i.e. relatively high temperatures (above the roughening transition [57] and sufficiently high to overcome typical diffusion barriers) and/or moderate fluxes. Low-temperature growth and, more in general, strongly out-of-equilibrium conditions are complex and fascinating regimes already in the homoepitaxial cases. The interested reader can find excellent treatments in several books [58][59][60]. Epitaxial breakdown in semiconductors, also not treated here, is nicely reviewed in Ref. [61].
The manuscript is organized as follows. In Section 2 we recall the basic thermodynamic and kinetic concepts behind the definition of an evolution model in the continuum. This allows us to review the seminal works of Asaro, Tiller, and Grinfeld (ATG) [62,63], and to point out the role played by non-linear effects. Then in Section 3 we discuss various numerical methods to implement the model with a dedicated focus (Section 3.2) on the diffuse-interface approach based on a phase-field description. In Section 4, we examine the major extensions of the simple ATG model towards a more realistic description of heteroepitaxy by considering the interaction with the substrate, responsible for wetting effects (Section 4.1), accounting for the mechanisms of intermixing at the surface (Section 4.3) and discussing the simulation of growth on patterned substrates (Section 4.4). The effects of plastic relaxation by dislocations are also considered in Section 4.5. Finally, in Section 5 we draw the conclusions and outline future perspectives.

Thermodynamics and evolution model
Heteroepitaxial growth by MBE or chemical vapour deposition, is generally characterized by close-to-equilibrium processes, ensured by high temperatures Downloaded by [187.160.112.235] at 22:00 10 June 2016 and/or low deposition fluxes [64]. Therefore, the dynamics is mainly driven by the tendency to decrease the system free energy 1 where U is the internal energy, S the entropy, V the total volume, T the temperature, N i the number of the i-type particles with a chemical potential defined by μ i = δG/δN i | p,T,N j =i [65]. It must be pointed out that growth conditions far-from-equilibrium might lead to very different behaviours compared with thermodynamic expectations [66] and eventually induce other instabilities (see Ref. [54], and references therein).

Evolution law
From a microscopic point of view the heteroepitaxial growth results from two main processes: deposition of adatoms from the vapour phase (i.e. the MBE beam) and redistribution of atoms by diffusion along the surface. Bulk processes can be considered to play no role due to their higher activation barrier [32,58].
Atomic diffusion by site-to-site hopping is a stochastic process, but in the thermodynamic limit it results into a net transport of material. This can be modelled in terms of surface diffusion [67][68][69]. According to the Onsager's linear law [70,71], the material flux j is driven by the gradient of the chemical potential along the surface profile, so that j = −M∇ μ with M the surface mobility and ∇ the surface-gradient (i.e. ∇ evaluated with respect to the local surface coordinates). As diffusion is a thermally activated process, M is expected to be described by an Arrhenius law: M ∼ exp ( − E d /k b T) with E d an effective energy barrier and k b the Boltzmann constant. The resulting evolution of the surface profile is conveniently inferred from the local normal velocity vn, defined at each point x of the surface asn[(dx/dt) ·n] withn =n(x) the outer surface normal. By considering the continuity equation we obtain with a source term corresponding to the external deposition flux. It is assigned as an external contribution independent of the surface energetics so that it can drive the system out-of-equilibrium. This is a simplified treatment of deposition neglecting a proper description of local fluctuations [72] expected however to provide a fair description following our high-temperature hypothesis. An illustration of the quantities which determine vn is shown in Figure 1(a).
Other transport mechanisms could also be involved, in particular the so-called evaporation/condensation process (also referred to as attachment/detachment) [67][68][69], where the local velocity is directly proportional to the chemical potential, i.e. vn = k(μ 0 − μ), with k a kinetic factor. For a conservative dynamics accounting for an "instantaneous" material redistribution along the profile, the condition μ 0 = μ = μdx/ dx must be fulfilled. Differently, μ 0 can be identified with the chemical potential of the vapour phase. In such a case a nonconservative dynamics is obtained, resulting in both deposition and desorption at the surface according to its local stability relative to the vapour. This contribution to mass transport can also be added to the surface diffusion Equation (1), providing a more general description of the profile evolution [73]. Deposition under equilibrium conditions can be modelled by replacing the external flux with the evaporation/condensation term [74,75]. It could be also useful to account for desorption effects which are, for instance, negligible for SiGe but not for GaAs [76]. It must be pointed out that the continuum dynamics here considered rigorously applies only to vicinal surfaces, where any arbitrary variation of the profile can be obtained by (barrier-less) step motion. On the contrary, on singular surfaces such as the typical {0 0 1} or {1 1 1} facets, a barrier for step nucleation must be surmounted [54]. Such barrier, however, vanishes above the roughening transition, i.e. under the typical growth conditions [57] we are interested in.

Energy contributions
The free energy for a heteroepitaxial system is essentially given by elastic G ε and surface (and/or interface) G γ energy: G = G ε + G γ . Here a single component system is assumed so that entropic contributions, mainly related with intermixing, are not considered. These will be taken into account in Section 4.3. The chemical potential, required in Equation (1) is obtained as the variation of G by adding one atom: μ = δG/δN = V a δG/δV , where V a is the atomic volume.
G ε is the energy stored in the volume of the solid because of deformations in the lattice structure, here determined by the misfit strain ε m between the lattice parameter of the film a f and of the substrate a s , i.e. ε m = (a s −a f )/a f . In the linear elasticity approximation (hereafter adopted) such deformations are described by the strain tensor ε = 1/2[∇u + (∇u) T ] where u(x) is the displacement of each point of the solid with respect to its position x in the reference state [77].
Forces produced by these deformations are defined by the stress tensor σ where the component σ ij corresponds to the force acting along the j-th direction, on the i-th face of an infinitesimal cubic volume within the solid. Stress and strain are related by the Hooke's law σ = C : ε, with C the fourth-order tensor of the elastic constants. This reduces to two independent values for an isotropic solid: C ijkl = λδ ij δ kl + μ(δ il δ jk + δ ik δ jl ) where λ and μ are the Lamé constants (equivalently, Young modulus E and Poisson ratio ν can be used), and δ ij is the Kronecker delta.
The strain/stress field of a deformed solid can be defined by solving the mechanical equilibrium equation ∇ · (σ − σ * ) = 0 with σ * the eigenstress tensor given by σ * ij = − k,l C ijkl ε * kl , where ε * ij = ε m δ ij is the eigenstrain, i.e. the reference strain corresponding to the in-plane lattice mismatch [77,78]. The additional boundary condition (σ − σ * ) ·n = 0 must be included to account for free surfaces. When considering dynamical conditions, i.e. the evolution of the profile in time, in principle, surface diffusion from Equation (1) and mechanical equilibrium should be solved simultaneously. However, the time scale for elastic relaxation is very quick compared to the time scale of diffusive motion, so that mechanical equilibrium can be assumed to be always instantaneously achieved for any profile configuration.
The elastic energy can be defined as G ε = ρ ε dV with the (volumetric) elastic energy density ρ ε = 1/2(σ − σ * ) : (ε − ε * ). The corresponding contribution to the local chemical potential is then In Figure 1(b) a map of μ ε at the surface of an island-like profile under compressive stress is reported. A cross section showing the ε xx field is shown in Figure  1(c), indicating a better strain relaxation at the top and more stressed regions at the borders. G γ is the excess of energy due to the presence of surfaces (compared to the bulk material) and can be computed as G γ = γ (n)dA, where γ (n) is the orientation-dependent surface energy density and dA is the infinitesimal surface element. In the absence of strain, γ (n) determines the equilibrium crystal shape (ECS), which can be derived by the well-known Wulff construction [79][80][81]. A dependence of surface energy on the local deformation can be included by means of the surface stress tensor s, defined by d(γ A) = A ij s ij de ij with A the surface area and e the surface strain tensor [82,83] related to γ by s ij = γ δ ij + ∂γ /∂e ij [84,85]. In solids, ∂γ /∂e ij = 0 and the ECS in presence of strain can be described by the Wulff-Kaishew theorem [86,87]. The surface energy contribution to μ is: where ξ is the Cahn-Hoffman vector [88,89], defined by ξ (r) = ∇(rγ (n)) with r = rn the position with respect to the origin. For isotropic surface energy ξ = γn, so that μ γ = V a γ κ where κ = ∇ ·n is the local curvature, corresponding to the sum of the principal curvatures κ 1 + κ 2 [89]. A map of μ γ ∼ κ is reported, for example, in Figure 1(d), on the same geometry of Figure 1(c). For anisotropic γ (n), preferential orientations are introduced into the ECS, as shown for instance in Figure 1(e). Two different conditions can be identified: weak anisotropy, when the ECS contains all the possible orientations, and strong anisotropy, when some orientations are missing in the ECS and sharp corners appear [80]. To distinguish among them a general criterion for missing orientations, based on the definition of ξ , can be found in Ref. [90]. In 2D, Equation (3) can be rewritten as μ γ = V aγ κ withγ = γ (θ) + γ (θ ) called surface stiffness, defined with respect to the profile orientation θ . In Figure 2, weak (panel a) and strong (panel b) anisotropy regimes are shown for a suitable choice of γ (θ). Whenγ > 0 for any orientation weak anisotropy is obtained, while strong anisotropy occurs whenγ < 0 for some θ . The equivalent of this latter case for a 3D structure is shown in Figure 2(c). When considering strong anisotropy, the surface diffusion equation becomes backward parabolic for any missing orientation θ in the ECS, giving rise to an ill-posed behaviour. Regularization procedures are then needed. A possibility consists in the convexification of the Wulff shape [91,92]. In 2D, this can be obtained by imposing the constraintγ ≥ 0 (also known as Frank's convexification). Another method consists of adding an energy term penalizing high curvatures (∼ κ 2 ), thus leading to rounded corners with no missing orientations [80,93,94], as sketched in Figure 2(d). The ECSs in Figure 1(e) have been obtained by surface diffusion evolution of a sphere, with strongly anisotropic γ (n) regularized by this corner term (see Ref. [95]).
Notice that, in the present approach, faceted shapes are not exactly represented by straight intersecting planes, but a continuum variation in the orientation is always present within each apparent facet. Moreover, their boundaries are always smeared over a short distance (at least because of the regularization) yielding a smooth transition in between. An alternative approach able to manage fully faceted shapes with selected orientations only, yet in a continuum description, has been proposed by Carter et al. in Ref. [96] and also applied to model heteroepitaxy [97,98].

Asaro-Tiller-Grinfeld instability
The main physics behind heteroepitaxy can be led back to the competition between elastic and surface contributions, the former favouring the growth of high AR structures and the latter opposing to any corrugation of a flat geometry, as illustrated in Figure 1(b), (d). Based on this concept, Asaro and Tiller [62], and independently Grinfeld [63,99] and Srolovitz [100], proposed a model to investigate the stability of the surface of a semi-infinite film under stress, by a linear stability analysis. Their theory, named as ATG instability, can be considered as the most essential explanation of the dynamics of heteroepitaxy [53][54][55]. Despite its simplicity, the ATG model was successfully applied to the study of the film instability observed experimentally for a wide variety of materials, such as helium crystals [101], polymers [102] and semiconductors [53,[103][104][105]. Surface waviness is typically observed for low misfit systems (provided that dislocations are not present), while for larger misfits, three-dimensional structures are generally obtained [11,12,106,107] and their formation is affected by nucleation process [108]. However, the instability defined by the ATG model can be viewed as the trigger for the formation of the initial nuclei leading to these 3D islands.
Let us consider a cosine perturbation of the surface of a (2D) film h(x) = a 0 cos (qx). Following Refs. [62,100], for isotropic elastic constants, to linear order in q the stress distribution at the surface is Summing the two contributions, the total chemical potential is obtained: μ ≈ V a U ε /2 + V a (U ε q − γ q 2 )a 0 cos (qx). The profile evolution by surface diffusion then reads Deposition is not included since in the linear regime it would just consist of a rigid shift in the vertical direction. Setting the solution for Equation (4) in the form of h(x, t) = a(t) cos (qx) it is possible to define the profile evolution in terms of the amplitude of the perturbation only: The sign of the amplification factor A(q) depends on q − q c . If q > q c the amplitude decays so that the flat surface is stable. If q < q c , the film is unstable as the amplitude of the perturbation grows exponentially. The behaviour of A(q) as function of q is shown in Figure 3(a), where each curve corresponds to a different q c . Notice that the fastest perturbation (where A is maximum) exists at The typical length scale for the evolution of the perturbation is = q −1 c . The stability analysis also shows a characteristic time scale τ = 4 /(MV a γ ), strongly dependent on the misfit (∼ ε −8 m ). As M is expected to obey an Arrhenius law, τ is extremely sensitive to temperature and decreases exponentially when T is raised.
The previous arguments are developed for a single q mode. However, a generic perturbation can always be decomposed in a Fourier series and, in the limit of linear elasticity, each component q is expected to behave independently from the others according to Equation (5). The resulting profile is then dominated by the fastest q and for a white noise (ideally comprising all q values) the rise of the fastest perturbation with q = q max (i.e. the fastest wavelength λ = λ max ) is expected. With this respect, the ATG model states that a stressed film is always unstable with respect to long-enough perturbation wavelengths λ > λ c (with λ c = 2π/q c the critical wavelength). However, when considering low misfit ε m and/or high surface energy cost γ , only very long λ are expected to grow, with extremely slow evolution (long τ ). In real systems, this condition might be physically unreasonable, and the flat film morphology can be considered practically metastable.
The description here reported for a 2D system directly applies also in 3D for the case of uniaxial stress, while for biaxial stress σ 0 = E/(1 − ν)ε m and U ε = 2(1 + ν)σ 0 ε m , meaning that the instability is shifted to longer q c by a factor (1 + ν) 2 . In this latter case, the perturbation should be considered with respect to both in-plane directions thus leading to cross-patterns dominated by the same q max perturbation.

Non-linear effects
The ATG model provides a good description of surface evolution for qa 1. In order to extend the approach to large perturbation amplitudes, it is necessary to solve Equation (1) by evaluating all derivatives along the surface (∇ ), consider the proper profile curvature κ = −h (1 + h 2 ) −3/2 in Equation (3), and more importantly, extend the calculation of the elastic contribution in Equation (2) to higher order in q.
Weakly, non-linear regimes have been explored by means of the approximation of the surface grooves with prescribed functions [104,109,110]. In the special case of cycloid-like profiles, an exact solution for the elastic field exists [111]. By means of numerical approaches, it is possible to analyse more in general the evolution as in Refs. [112,113], where the formation of cusps separated by deep trenches was also linked to crack propagation. Several other studies are reported in literature (e.g. in Refs. [53,[114][115][116][117][118][119]). More accurate calculations involving the numerical solution of the mechanical equilibrium equations by Finite Element Method (FEM) have also been exploited [120,121]. All of these approaches, however, agree in the tendency to form cusp-like structures with rounded tops and sharp trenches, yielding the best tradeoff between strain relaxation and total surface area. Indeed, the cycloid surface is found to be the most efficient stress concentrator at a fixed wavelength [111], yielding an almost strain-free surface except at the cusp-point. The formation of cusps ends into stress singularities, leading to a divergent propagation of the trenches within the solid. In the case of tensile stresses, this corresponds to the initiation of surface cracks, leading to spontaneous Griffith brittle fracture [111]. Nucleation of defects (dislocations) at the cusp tips, favoured by the high local stresses [122], might contrast material failure by enabling the alternative channel of plastic deformation for strain release (ductile fracture). In Figure 3(b) the evolution of a surface profile, including non-linear contributions, is shown. Evidently, the initial cosine profile is firstly amplified according to the linear theory (here λ > λ c ). Then, for larger amplitudes, deviations occur and lead to an asymmetric profile with cusp-like structures separated by deep trenches, eventually diverging. This feature is well consistent with experimental observations [53,102,104,105,123]. An example is reported in Figure 3(c) for a surface morphology obtained by annealing a SiGe/Si film with low Ge content [124].

Computational methods
In order to investigate the evolution of heteroepitaxial systems, several computational methods have been developed to solve Equations (1)-(3) or approximations thereof. We here provide a general overview of different approaches and describe one, the diffuse-interface approach, in detail.

Comparison of computational approaches
Mathematically, Equations (1)-(3) define a free boundary value problem, where the surface profile is unknown and has to be determined. Only recently, first analytic results for the highly non-linear sixth-order partial differential Equations (1)-(3) (including corner rounding regularization) have been obtained [125,126]. Most numerical treatments consider only specific aspects of the whole problem. In Refs. [121,127,128], Zhang et al. developed a model, based on FEM, to solve both differential problems of mechanical equilibrium and surface diffusion, but only for weak surface anisotropy. A more advanced treatment of weak surface anisotropy using parametric finite elements is considered in Ref. [129]. For strong surface anisotropy, we refer to Refs. [130,131] for parametric finite element approaches and to Ref. [132] for a finite difference method.
If the surface morphology can be described by a height function, as in Equation (4), the situation simplifies. For shallow profiles, elasticity can be generally treated by means of the so-called flat-island approximation [133,134]. It consists of describing uniformly strained layers with a thickness h(x) above the substrate, in terms of a stress field σ ij = σ 0 h(x)δ ij at the surface. Then, the linear elastic response of the medium at a point x due to a force located at another point x 0 is generally given by the displacement u i (x) = G ij (x − x 0 )∂σ jk /∂x k where G is the Green-function tensor. For a semi-infinite film the analytic solution for G exists [77]. It can also be corrected in order to account for more accurate force distributions [133,134]. From the resulting displacement field u, strain is calculated and hence μ ε , according to Equation (2). Differently, an Airy function formulation [77] can be exploited. Spectral methods have been widely employed to solve the resulting system e.g. in Refs. [117][118][119]135,136], profiting of highly optimized algorithms, such as the Fast Fourier Transform. However, height function methods are not applicable to complex structures with overhangs, e.g. mushroom-like morphologies or wetting angles above 90 • .
The methods described so far treat the surface profile explicitly. An alternative approach consists of an implicit description of the geometry, where the surface profile is located as a level-set of an introduced function. Despite the increased computational cost resulting from the reformulation of the problem in a higher dimensional space, this method offers several advantages. Indeed, the reformulated problem can be solved in a fixed domain, using standard discretization techniques, and allows one to naturally tackle topological changes. Strongly anisotropic surface diffusion has been modelled e.g. in Refs. [137][138][139], using a level-set method [140], but the elastic field was not taken into account therein. Also diffuse-interface or phase-field-type models describe the surface profile implicitly. Here, Equations (1)-(3) are only approximated by assuming a finite thickness of the surface region of order . In the limit → 0 the original, sharp-interface problem can be recovered [141][142][143][144][145]. Phase-field approaches provide a robust numerical tool to solve surface evolution problems, and open the possibility to add more and more of the relevant phenomena in a consistent way. The basis is the classical Cahn-Hilliard model [146], extended in order to include the major traits of heteroepitaxial growth, namely anisotropy in the surface free energy density, elastic misfit strain, diffusion restricted at the surface and deposition flux, as illustrated in Section 3.2. Further extensions will be also addressed in Section 4.

Diffuse-interface approach
Let us consider the two-phase system formed by the solid and the surrounding vapour. An order parameter ϕ is introduced to identify the two regions, respectively, as ϕ=1 and ϕ=0. A smooth transition in ϕ is assumed at the boundaries between the two phases, thus providing a diffused description of the solid surface (transition layer). The classical formulation of the diffuse-interface equations for isotropic systems is based on the free energy where is a small parameter that is a measure of the interface thickness and f (ϕ) = 18ϕ 2 (1 − ϕ) 2 is a double-well potential. An example of a surface profile represented within the phase-field framework is illustrated in Figure 4. Approximations for surface diffusion based on this energy have been considered in Refs. [141,147,148]. The evolution law reads where M(ϕ) is a mobility function, localized near the interface and vanishing in the solid and in the vapour phase. For the present parameter choice, M(ϕ) = (36/ )ϕ 2 (1−ϕ) 2 , if the diffusion constant is one. The chemical potential is defined as μ = δG[ϕ]/δϕ and f (ϕ) is the partial derivative of the double-well f (ϕ) with respect to ϕ. The order of convergence can be improved if μ is replaced by g(ϕ)μ, with g(ϕ) localized near the interface, as g(ϕ) = 30ϕ 2 (1 − ϕ 2 ) (see Refs. [141,145]). In all cases the profile across the interface in equilibrium is well described by ϕ (r) = 1/2 1 − tanh 3r/ , with r = r(x) the signed distance function to the 0.5 level-set of ϕ.
To model crystal growth, deposition should be considered, which here is rendered by a prescribed material flux from the vapour phase = −|∇ϕ|n · .   Figure 1(b)). On the left: the ϕ continuous field; on the right: mesh discretization with local refinement at the interface.
|∇ϕ| thereby restricts the incoming flux to the interface region and integrates in the whole domain to the size of the interface. The projection onto the outer pointing normaln = −∇ϕ/|∇ϕ| results in a net income which enters Equation (7) as an additional term for ∂ϕ/∂t. The evaluation of might be a delicate task as it may account for vapour-phase details, such as mechanisms of fluid-dynamic transport or flux shielding due to the surface profile itself. Numerical approaches for the whole system of equations have been considered for example in Refs. [141,149].

Anisotropic interface energy
Following Ref. [150], anisotropy can be considered in the surface energy as The first term is the interface energy equivalent to G γ introduced in Section 2.2. The additional term is a De Giorgi-type phase-field approximation of a Willmore regularization [151][152][153], with √ β a small-length scale over which corners are smeared out. For → 0 the energy reduces to the surface energy [γ (n) + (β/2)κ 2 ]dA, which can be viewed as the next higher order approximation of a more general surface free energyγ =γ (n, κ, . . . ) = γ (n) + (β/2)κ 2 + · · · [80,93,154]. Such a regularization prevents the model to become unstable for negative surface stiffness as discussed in Section 2.2. The corresponding evolution equations read Downloaded by [187.160.112.235] at 22:00 10 June 2016 where P = 1 −n ⊗n is the projector operator. The approximation results from f (ϕ) ≈ 2 /2|∇ϕ| 2 , which is valid in the asymptotic limit of → 0.
In contrast with other approaches [155], the above formalism ensures that the interface thickness remains independent of orientation. Furthermore, it retains the same form of the isotropic case, allowing for the combination with the Willmore regularization term. A fully customizable formulation for γ (n) has been proposed in Refs. [95,156] to account for the complexity of realistic morphologies.

Elastic contribution
As discussed in Section 2.2, surface motion should account for the additional driving force towards the reduction of elastic energy in stressed films. To avoid solving the mechanical equilibrium equation with the free-surface boundary condition (see Section 2.2) in a time varying domain, a diffuse-interface approach can be conveniently used [141,157,158]. Thus, the problem is extended to the whole domain by introducing material properties dependent on ϕ: is an interpolating function, with value 0 in the vapour (ϕ = 0) and 1 in the solid (ϕ = 1). At the interface H(ϕ) gradually changes from 1 to 0. In the diffuse-interface approach, the elastic energy density for an isotropic material reads following the definition introduced in Equation (2) with ε * (ϕ) = ε m (ϕ). Mechanical equilibrium is assumed to be instantaneously achieved. From the computed displacement vector u, the elastic energy density (ϕ, u) is obtained. Then, its derivative with respect to ϕ can be computed: and included as an additional term in the g(ϕ)μ definition in Equation (9). The resulting model is related to the Cahn-Larché system [159] and has been intensively studied, but mostly with isotropic surface energies only [160][161][162][163].

Numerical issues
The derived set of equations for the phase-field variable ϕ and the displacement field u provides a complete continuum model of heteroepitaxy. A wide variety of similar models, some accounting for further details, can be found in the literature, e.g. in Refs. [115,[164][165][166][167][168][169][170][171]. Different numerical methods can be used to solve them. FEM is a very effective way to obtain an accurate solution of the partial differential equations here defined, profiting of their representation in a weak form. However, adaptive mesh refinement [172] (as shown in Figure 4(b)) is highly desirable. Indeed, the mesh has to resolve the transition layer, which is of order , as well as the corner length scale √ β when considering strong anisotropy. These requirements of spatial resolutions also set restrictions on the time stepping, with an upper bound, such that the interface never evolves over a whole mesh cell within one time step. Furthermore, solving the full model in 3D requires high performance computing and highly scalable algorithms.

Towards realistic modelling
The simple ATG model presented in Section 2.3 provides an effective explanation for the occurrence of island growth. However, it lacks several physical aspects which play a crucial role in real systems. In this section, we review the main contributions that should be incorporated in the model for a more reliable description of heteroepitaxial growth.

Wetting effects
One of the major simplifications of the model discussed in Section 2.3 is that a semi-infinite film is assumed and the existence of a substrate underneath is accounted only for the definition of the lattice mismatch. However, real films extend for finite thicknesses h on top of the substrate. If the elastic constants differ from those of the film, a different redistribution of stresses can be achieved, thus affecting the surface relaxation and hence the film stability [173][174][175][176]. If the substrate is stiffer, the film stability is enhanced and the critical wavelength is larger. In the limiting case of a rigid substrate it has been shown that there exists a critical thickness h c = (1 − ν)/q c (with ν the Poisson ratio of the epilayer and q c given in Section 2.3) below which the film is stable for any perturbation wavelength [110,176,177]. On the contrary, when the substrate is softer than the film, shorter critical wavelengths are obtained and the profile roughening is favoured. This is the case of compliant substrates [178] or, in the extreme case, viscous substrates [179]. These effects are all dependent on the film thickness and are effective only for thin enough films, as reported in Refs. [173,174]. Even if a sharp interface between film and substrate is assumed, both elastic and surface properties are modified in the proximity of the interface. In particular, it has been shown by considering non-linear elasticity, that the relaxation of a film can be enhanced very close to the substrate [180][181][182], recovering the constant value predicted by linear elasticity only a few monolayers (MLs) above the interface. The effect of such a contribution on the profile evolution is discussed by Eisenberg and Kandel in Refs. [181,182]. Surface energy might depend as well on the film thickness. Indeed, ab initio calculations for Ge/Si(001) [45,183,184] have shown that the surface energy density of the film undergoes an exponential decay from the substrate value (γ of Si) to the one expected for the thick-film limit (γ of Ge), within 3-5 MLs from the interface. A γ (h) function fitting such a behaviour has been introduced in the evolution equations by Chiu and Gao [185] and others [119,186]. Smooth transitions between substrate and film properties through a small, finite region at the interface have also been considered [118,175,176] and represent the natural choice when considering diffuse-interface models [141,165,168]. These approaches, named as boundary layer transition models, might be physically related with the existence of a mixing region between the two materials (see Section 4.3), where their properties combine to intermediate values. The asymptotic case of vanishing thickness of the transition region, named glued-wetting-layer model, has been investigated by Spencer in Ref. [187].
Let us consider a generic dependence of the free energy G on the film thickness h: G = G(h), mimicking the behaviour of surface energy γ = γ (h) and/or elastic energy ρ ε = ρ ε (h). An additional contribution must then be added to the local chemical potential as μ wet = ∂G/∂h ≈ W 0 + W(h − h 0 ), where the linear expansion around the average film thickness h 0 = h has been considered, with W 0 = ∂G ∂h h 0 and W = ∂ 2 G ∂h 2 h 0 . Correspondingly, the Equation (4) describing the profile evolution becomes [188]: The wetting contribution behaves as a stabilizing term provided that W > 0. Furthermore, if the condition γ W > U 2 ε /4 is fulfilled, the flat profile is stable against any perturbation wavelength. By considering that W = W(h 0 ), the effect of wetting is expected to vary with the film thickness and become negligible in a thick-film limit. In particular, if W > 0 and decreasing with h 0 , there exists a critical film thickness h c above which the film becomes unstable for some q. For a better understanding, let us look in more details to the case of γ = γ (h). In such a case, W = ∂ 2 γ ∂h 2 h 0 provides a film stabilization if the surface energy of the film γ f is lower than the substrate one γ s (i.e. γ f < γ s ). By assuming an exponential law for γ (h) = γ f + (γ s − γ f ) exp ( − h/d), the critical thickness can be defined as h c = −d ln (U 2 ε d 2 )/(4γ f (γ s − γ f )) . As shown in Figure 5(a), for h 0 < h c the film is stable for any wavelength while for h 0 h c the typical curve for the ATG instability is recovered. Interestingly, for values of h 0 h c , only wavelengths in a finite range are unstable, i.e. q min < q < q max , so that the wetting contribution has the effect of stabilizing the long wavelengths.
Based on this analysis, a profile perturbation is expected to rise when annealing a flat film above the critical thickness h 0 > h c . As for the standard ATG model (Section 2.3.1), the evolution is dominated by a fastest wavelength, which however depends on the film thickness itself, increasing for thin films as evident in Figure 5(a). The presence of the wetting term also modifies the evolution of the perturbation as soon as it starts to dig trenches towards the substrate. Indeed, these regions are progressively stabilized by the increasing values of W and γ and stop propagating. A wetting layer of thickness h wet < h c is then preserved while the perturbation breaks up into islands on top of it [119,189]. The wetting contribution, coupled with non-linear stress, thus permits to avoid the formation of the cusp-singularities observed in the non-linear ATG models (see Section 2.3.1), allowing for the description of long-time dynamics, including island coarsening (see below). This evolution matches the description of the Stranski-Krastanow growth mode.
When considering weak or even non-wetting conditions W < 0 (i.e. γ f > γ s ), the formation of a wetting layer is no more favoured. The presence of the unstressed substrate ensures that the film perturbation does not penetrate inside it but breaks up into islands directly on the substrate surface, which can be partly exposed [165]. The Volmer-Weber growth mode is then obtained.
More in general, the balance between γ f , γ s and the interfacial energy γ sf between the substrate and the film should be considered when addressing the wetting conditions. This is discussed, for instance, in Ref. [167] where a phasefield approach is used to describe the three-phase system consistent in substrate, film and vapour, explicitly accounting for all the relative interfacial energies.

Growth
As far as a semi-infinite strained film is considered, as in the ATG model, the addition of a deposition term dh 0 /dt ≈ to the evolution Equation (4) is trivial since elastic and surface energy are invariant under the Galilean transformation h → h + h 0 (t). This holds true also when the film and substrate have the same elastic constants, with no wetting contributions. In the more general case, the dependence of the film instability on its thickness h 0 , discussed in the previous Section, breaks this translational invariance so that the variation of h 0 during the deposition process dynamically affects the energy balance [173,174,190]. In the linear approximation, for a small profile perturbation, the time evolution of the amplitude h of the Fourier component of wavenumber q can be written as: Downloaded by [187.160.112.235] at 22:00 10 June 2016 where A(q, t) = A(q, h 0 (t)) is the amplification factor for the q component, at the average film thickness h 0 resulting after deposition, at time t. In particular, in the case of wetting contributions, as the film grows thicker, the system continuously moves from one curve to the other of Figure 5(a). This implies the passage from stable conditions, for h 0 < h c , to the rise of the instability, with q max getting larger while approaching the thick-film limit of the ATG model, as shown in the evolution sequence of Figure 5(b). Similar effects are expected when different elastic constants are considered, especially for the case of rigid substrate, as discussed in Ref. [173].
Kinetics plays a crucial role in controlling whether the instability is going to become apparent or not with respect to the growth process. Indeed, deposition sets a reference time scale −1 to be compared with the instability time scale τ ∼ A −1 . The relative growth rate = A − /h 0 (t) can be considered [173,174] to identify the kinetic onset of the instability: for < 0, the perturbation grows too slowly and its amplitude shrinks relative to the growth of the film thickness, while it becomes dominating for > 0. Another criterion consists in the definition of an arbitrary, small threshold below which the amplitude is assumed to be not appreciable [19,190]. A dynamical critical thickness h c can then be identified as the average height h 0 reached by the film before the perturbation amplitude exceeds this threshold. h c is proportional to the /M ratio, so it is expected to be very sensitive to the growth temperature, due to the exponential dependence of the surface mobility M.
It is important to notice that, in the case of wetting effects, the existence of a dynamical critical thickness can delay the onset of (apparent) island growth to larger thicknesses with respect to the thermodynamic critical value h c , defined in the previous Section. This kinetic stabilization may permit the growth of flat films, thick enough to overtake the critical thickness for plastic relaxation prior to islanding. Planar, dislocated films, hundreds of nm thick, are indeed fabricated by exploiting low temperature and high deposition rate to skip the thermodynamic regime of island growth [191][192][193].

Coarsening dynamics
Thanks to the wetting contribution, healing the divergent behaviour of the nonlinear ATG instability (see Section 2.3.1), island growth can be explored over long time scales, as shown in Figure 5(c). This permits not only to observe the onset of their formation but also the subsequent exchange of material between them. In agreement with experimental observations [194], island coarsening occurs, following Ostwald ripening dynamics, i.e. bigger islands grow at the expenses of the smaller ones. These are then observed to disappear, leaving an almost flat wetting layer at the thickness h wet in between the larger ones. This behaviour is illustrated by the late simulation stages reported in Figure 5(c).
As the island volume increases the chemical potential monotonously decreases [119], so that a non-interrupted coarsening is expected, eventually leading to the formation of a single island [119,186,189,195]. This feature is inherent to isotropic islands, for which a continuum of island shapes is allowed [196], and does not hold when faceted structures are considered (see Section 4.2). Moreover, island chemical potentials can be altered by the island-island interactions [197], corresponding to the strength of the elastic field induced by one island at the location of the other. Such an effect is, however, significant only at short distance d and for large cross-sectional island area A. Indeed, its strength scales as A/d 2 , and it is responsible for mutual island repulsion. Coarsening dynamics occurs on a long time scale so that it can be strongly affected by deposition, especially for high growth rates [128,190].

Faceting and shape-transitions
Heteroepitaxial islands are typically characterized by well-defined faceted shapes [12,15], as shown for example in Figure 6(a). The existence of these preferential orientations is generally due to anisotropic behaviours of the system. Simulations including elastic anisotropy are not observed to produce relevant morphological changes compared to the isotropic system [165] while surface anisotropy is found to yield faceted shapes, exposing those planes corresponding to minima in the surface energy density function γ . Such facets typically include those expected from the ECS of the unstressed material (see discussion in Section 2.2) but strain effects might affect their relative extensions [198] and eventually induce the formation of additional facets (e.g. the common {105} facets observed on Ge islands [12], whose presence in the crystal shape is possible only because of a strong stabilization effect due to strain [43][44][45]).
The effect of surface anisotropy on the stability of the stressed film can be easily inferred in 2D by replacing the isotropic γ used in the linear stability analysis in Equation (4), with the surface stiffnessγ [199]. Hence, if the orientation of the flat surface (θ = 0) is a stable one, i.e. a minimum of γ ,γ is maximum so that the stability is enhanced, the critical wavelength becomes larger and the critical thickness of the wetting layer increases. More precisely, this stabilization is just a metastable condition [181] as perturbations of large enough amplitudes can introduce local orientations far from the θ = 0 minimum, for which the stiffness is lower and such to destabilize the film, even below the critical thickness defined from the linear stability analysis usingγ . On the other hand, if the flat orientation is not stable, it is strongly destabilized and tends to spontaneously evolve into a faceted profile. The critical thickness is lowered or even nullified as the strength of the anisotropy is increased [200].
Once the instability is triggered, growth of faceted islands is observed [127,156,[201][202][203][204], as shown, for example, in Figure 6(b) and (c). The existence of preferential orientations forces the possible island geometries to be bounded by such planes. Because of this, islands are expected to grow self-similar [201] and eventually convert into different shapes as additional facets are formed [156]. Coarsening between close-by islands is also affected by this shape stabilization. At the early stages, when facets are not fully developed, Ostwald ripening causes the disappearance of the smaller islands, but these dynamics is almost interrupted after fully faceted structures are formed [119], in contrast with the isotropic case discussed in Section 4.1.2.

Intermixing
Heteroepitaxy inherently involves at least two different chemical species forming the substrate and the film. Alloy materials (e.g. Si 1−x Ge x , In 1−x Ga x As, etc.) are very common, allowing one to finely control the lattice mismatch ε m by tuning the difference χ = c f − c s between the composition of the film c f and of the substrate c s , i.e. ε m = (a(c s )−a(c f ))/a(c f ) ≈ χε 0 (with ε 0 the misfit between the pure components). According to the ATG model (Section 2.3), this results into different critical wavelengths for the film instability: the smaller is χ the longer is the unstable wavelength, in agreement with the experimental observation of larger islands when considering diluted alloys [28]. A variation in the surface energy density γ is also expected when changing the film composition but the effect of the misfit is assumed to be dominant (q c ∼ ε 2 m /γ ). This trend is the same illustrated by the different curves in Figure 3.
Experimental evidences [28,29] clearly show the occurrence of a complex intermixing dynamics during the growth process which implies both non-uniform composition profiles [205] and time dependence of the composition field. In order to account for these effects, an excess contribution must be included in the definition of the system free energy: G mix = H mix − TS mix . An enthalpic contribution H mix arises from the atomic size differences of the chemical species, responsible of solute stresses in the solid, while the mixing entropy S mix is the configurational entropy −k b ν c ν ln c ν , where c ν = c ν (x) is the local composition with respect to component ν (i.e. the atomic fraction N ν / ν N ν ) and k b is the Boltzmann constant. By assuming that each component behaves independently, according to Equation (1) its local flux is given by j ν = −c ν M ν ∇ μ ν where M ν is the mobility of the ν-th component and μ ν is the corresponding local chemical potential. It consists of: Surface μ γ and elastic μ ε contributions are the same defined in Section 2, provided that γ = γ (c ν ) is considered and that the problem of mechanical equilibrium for the strain field is solved by considering the non-uniform misfit strain map ε m (c ν ) within the whole solid. The third term accounts for the solute stress ∼ ∂ρ ε /∂c ν (here given for a binary alloy) [206,207], due to a relative difference in atomic sizes η ν = (a ν − aν)/a(c ν ) with a ν , aν and a(c ν ) the lattice parameter of the pure components ν andν and of the alloy at composition c ν . trace(σ ) is the local hydrostatic strain. The fourth contribution comes from the entropic term and favours alloying.
As demonstrated by ab initio calculations [32,33], all mechanisms of atomic exchange responsible for intermixing are possible only within the first few atomic layers at the surface. In the bulk below, the activation barriers for atomic motions are too high and bulk diffusion is hindered. In order to model this condition, the system can be separated into a mobile region at the surface, extending for a small depth w (few atomic MLs), and the bulk underneath, characterized by a "frozen-in" composition profile, as sketched in Figure 7(a). The evolution of the system then results from the coupled dynamics of morphology and composition at the surface: The first equation defines the local normal velocity as resulting from the cumulative fluxes of each component ν. The second equation defines the local variation of the surface composition, resulting from the balance between the surface flux of component ν and the material exchange with the bulk region immediately below, at local composition c b , represented by the last term. In a first approximation [206,207], it can be assumed that c b ≈ c ν . However, as discussed in Ref. [208], this is true only where the surface is locally growing (vn > 0), leaving new material in the bulk. On the contrary, where the surface is retreating (vn < 0), material is extracted from the bulk, at a composition c b of the material buried immediately below the surface, possibly determined from previous growth stages. A similar approach accounting for a fully faceted description has been proposed in Ref. [97]. In Ref. [209] a Ginzburg-Landau free energy is introduced to describe the thermodynamic state of the alloy and a continuum composition field, extending from the bulk to the surface, is considered with mobility restricted at the surface to mimic surface diffusion. Attempts to couple both morphology and composition evolution in a phase-field description are proposed in Refs. [166,210]. According to Equation (16), deposition of material at composition c d on a flat substrate (at initial composition c s ) leads to a continuous variation c(h) = c d + (c s − c d ) exp ( − h/w), with h the height above the substrate interface. This provides an explanation to the transition-layer models introduced for the wetting energy in Section 4.1. Since h = · t the thickness of this diffused interface is proportional to the ratio w/ .
If a slight perturbation of the profile is considered, the film instability is expected to rise during deposition. However, the characteristic time scale for such instability strongly depends on the surface composition as τ ∼ c −8 (see Section 2.3). Therefore, at the early stages, when the composition is still diluted, it is too slow compared with the time scale of deposition ∼ w/ and the planar morphology is kinetically stabilized, as discussed in Section 4.1.1. Only by extending further the deposition, increasing the content of the strained species at the surface, it is possible to make the instability time scale comparable with deposition and hence observe its evolution into islands. As a result, there exists an apparent critical thickness for the instability to become appreciable, which strongly depends on the composition and weakly on the deposition rate [19], as illustrated in Figure 7 [19]). The solid lines, from bottom to top, correspond to increasing growth rates. The dotted line shows the behaviour when surface segregation at the top layer is considered. The inset reports a comparison with experimental data (points). (c) Simulation of island formation, growth and coarsening effects including intermixing (from Ref. [213]). A flat film (locally perturbed at the centre of the cell), with a composition profile matching the result of deposition, is set as initial condition and its evolution by annealing is simulated.
observation of increasing wetting layer thickness for highly diluted alloys [123,211]. Effects of thermodynamic surface segregation, redistributing the material within the surface layer in order to accumulate the species of lower γ (e.g. Ge in SiGe or In in InGaAs) in the top ML (see experiments in Refs. [211,212]), have been investigated as well [19] and were proved just to delay the occurrence of the instability (see dotted curve in Figure 7(b)).
Once the instability has been triggered, it grows in amplitude and islands start to form [213], as illustrated in Figure 7(c). Profile minima dig into regions of different compositions, tending to that of the substrate, and hence they gather material with lower misfit strain. As a consequence, the tendency to form deep trenches, as in Section 2.3.1, is damped and separated islands are formed. This effect is analogous to that of the wetting term discussed in Section 4.1. Low misfit material is extracted from the regions close to the bulk, redistributed along the surface and incorporated at the growing islands. The corresponding reduction of strain permits an enlargement of the island bases, reflecting an increase in the unstable wavelengths. Eventually, coarsening effects occur, favouring the bigger islands at the expenses of the smaller ones, which are pushed far away Downloaded by [187.160.112.235] at 22:00 10 June 2016 and become asymmetric both in shape and composition [214]. Based on this approach, islands are expected to have a non-uniform composition profile, with the outer region rich in the substrate component coming from trench digging. Both theoretical [26] and experimental studies [205] have shown that a typical composition map of islands is rather non-uniform with the top region enriched by the strained component and the lateral regions, close to the trenches, with a larger amount of the substrate component. This evidence is reproduced also by kinetic models as in Ref. [98,215].
The role of solute stress was studied in detail in Refs. [206,207] and it is found to enhance the growth rate of the instability, with more dramatic effects when considering intermediate deposition rates. The most unstable wavelength λ max is also predicted to shift towards shorter values. To understand this behaviour, let us consider a compressed alloy film, with η = 0. According to the ATG model, the morphological instability induces an undulation of the surface profile, with a tensile response at peaks and a further compression at valleys. Bigger atoms are then driven to accumulate at the former locations while the smaller ones are left at the latter, yielding a modulation in the composition along the surface. This redistribution enhances the strain non-uniformity, thus speeding up the growth of the instability itself. For very large values of η, the growth rate of the instability can become unbounded, even without misfit strain, producing surface decomposition of the alloy despite its stability in the bulk. This effect is known as compositional-stress instability or kinetics instability [216].
Another aspect to consider for a multi-component system is that the mobility of each species might be quite different (e.g. for SiGe, Ge diffusion is ∼10-100 times faster than Si at typical growth temperatures [217]) and depend on the local composition. This could induce kinetic segregation effects [208] and modify the film stability and evolution [218,219]. In Refs. [207,220] it is shown that a difference in surface mobilities can result into a stabilization of a strained alloy film and even suppress the stress-induced morphological instability.

Substrate patterning
The patterning of the substrate is a well-established method to direct and control the growth of islands into ordered arrays, with improved shape and size homogeneity [221][222][223][224]. From the point of view of the modelling, the effect of patterning turns into a modulation of the chemical potential. Continuum models are best suited for the study of the mechanisms involved in the growth on patterned substrates as long-range interactions, extending for distances larger than the pattern features, typically ∼100-1000 nm, are to be considered.
The most common way for substrate patterning consists in using a nonflat geometry with regular sequences of pits or mesas, produced on the surface by lithography and etching techniques. This induces both capillarity effects, towards profile flattening, and elastic effects, which alter the strain relaxation according to the relative positioning of the island on the substrate [225]. Several experiments (e.g. [223]) have demonstrated that both positional order and size uniformity can be obtained (see Figure 8(a)), provided that suitable growth conditions are matched [226]. Simulations studying the island growth on wavypatterned substrates have been reported in Refs. [225,[227][228][229][230]. Therein, the strain modulation due to the substrate undulation is accounted for by Fourier spectral approaches. When considering shallow pits with size λ s λ max [225], capillarity effects are found to first sustain the growth of the same wavelength λ s of the substrate pattern but soon after, the most unstable wavelengths around λ max emerge. As the position at the bottom is the most favourable, close to equilibrium, i.e. for low deposition flux, a single island is formed there. When higher deposition rates are considered, the island can grow slightly off-centre or multiple islands can form on the pit sidewalls, and resist to coarsening [225]. This trend is shown in Figure 8(b) by varying the deposition flux. For smaller pits, with size comparable to the most unstable wavelength λ s ≈ λ max , the perturbation is found to grow in anti-phase with the substrate profile, except for the case of very thin films, where islands localize at the substrate tops [227][228][229][230]. Similar results have been obtained for more isolated gaussian pits [231] and are consistent with simulations obtained by FEM, also including surface and elastic anisotropies, for both pits and hump geometries [75,232].
The combined effect of substrate morphology and compositional inhomogeneities has been inspected in Ref. [219], where an enhancement of the pit-filling prior to the island growth has been discussed and compared with experimental evidences (see Figure 8(c) for an example of island growth inside the pit with surface anisotropy).
The phase-field approach is well suited to account for arbitrarily complex substrate geometries using an auxiliary phase-field function to trace the film/substrate interface in a diffused way [164]. The same method has been exploited to model the effect of inclusions embedded in the substrate [233].
The modulation of the surface strain induced by buried islands [234] was proved to be an alternative way to pilot island ordering, without pre-pattetterning of the substrate. Self-aligned island superlattices have been obtained [235] by deposition of multiple film layers separated by substrate spacers. Zhang et al. in Refs. [74,195,236] reported extensive simulation results of island growth in stacks, accounting for the influence of the strain modulation due to the underlying buried layers. In particular, the authors pointed out the role of the strain energy density minima and maxima in the positioning of the islands to obtain ordering and vertical alignment. Peculiar conditions were also found to induce an alternate disposition of the islands into a body-centred tetragonal superlattice, with improved uniformity and regularity (see Figure 8(d)). Correlated growth of islands has been studied also for deposition on the top and bottom surfaces of an ultrathin substrate [237].  [223]). X and Y correspond to two 110 directions. (b) Simulation results of island growth into a cosine-patterned substrate with λ s λ max for different deposition fluxes (from Ref. [225]). (c) Simulation of Ge island formation and growth on a Si pit including faceting and intermixing, during Ge deposition (adapted from Ref. [219]). The anomalous pit-filling observed at the initial stages results from a speed up of Ge diffusion, enhanced by surface segregation as explained in Ref. [219]. (d) FEM simulation of vertical island ordering for repeated growth of five layers of strained material separated by thin substrate spacers (adapted from Ref. [74]). Notice the alternated alignment of islands between consecutive layers.

Dislocations
Plastic relaxation provides an alternative path for achieving strain release. The growth of thick films at low temperature and high deposition rate typically does not lead to misfit strain relaxation by 3D islands but via the formation of dislocation networks [191][192][193]. Even when islands are first formed, dislocations are found to appear inside those with large volumes [21][22][23]25], partially releasing the residual misfit strain not relaxed by elasticity.
Dislocations are line defects, inducing a permanent modification in the crystal lattice as sketched, for instance, in Figure 9(a), quantified by the Burgers vector b. The elastic deformation induced by the presence of a dislocation can be computed within a continuum approach by means of the linear elasticity theory. Indeed, analytic functions are known for both bulk materials [238] and flat films [239][240][241]. An example of this latter case is reported in Figure 9(b) for the ε xx component. It is worth mentioning that a singularity is present at the dislocation line, which can be regularized according to Ref. [242]. The definition of the local strain field for a generic geometry is quite more complex as it requires to numerically solve, e.g. by FEM, the problem of mechanical equilibrium in presence of the dislocation. This can be done by explicitly imposing the un-relaxed displacement field introduced by the dislocation [243,244], or, equivalently, by computing the Plots in panels c and d, here obtained by means of an improved FEM approach for the calculation of strain relaxation including dislocations, reproduce the results originally presented in Ref. [246].
changes of the analytic bulk strain/stress fields induced by the free surfaces of the actual geometry [41,245]. The amount of strain relaxed by a single dislocation is determined by b. If an array of dislocations spaced by d along x is considered in a flat film, at the interface with the substrate, the initial misfit strain ε m is reduced to an average residual strain of ε xx = ε m − b x /d [192]. A characteristic spacing of d 0 = b x /ε m is then sufficient to obtain (on average) full strain relaxation (notice that d < d 0 would imply an increase in strain, opposite to the misfit). For a typical dislocation network aligned in both x and y directions, a similar contribution has to be considered for ε yy .
In a first approximation, the effect of the dislocation array on the stability analysis of the flat film (see Section 2.3) can be inferred by recalculating the critical wavenumber q c in Equation (5) with ε xx replacing the misfit strain ε m in U ε . Evidently, the lowering of the strain has a stabilizing effect deferring the instability to longer wavelengths. For the peculiar case of d = d 0 the flat film is expected to be stable.
In a flat-island approximation, Jonsdottir and Freund in Ref. [246] showed that for dense arrays, i.e. small dislocation spacing d/d 0 , the surface profile rearranges into equilibrium shapes, as illustrated in Figure 9(c). For this lowstrain condition, the instability time scale τ is extremely slow and unstable wavelengths are very long, so that island formation is practically suppressed and the surface roughness just increases to a constant value, as shown in Figure 9(d). By increasing dislocation spacing d/d 0 the plastic relaxation becomes less efficient and the profile corrugation raises in order to release strain. Consequently, the film instability sets in and islands are formed, as shown in Figure 9(d).
The effect of a dislocation on the shape of an island has been investigated at equilibrium in Ref. [175]. As expected from experiments [21,22], dislocated islands are characterized by smaller aspect ratio and larger bases, due to the lower elastic energy. Models coupling plastic relaxation with island growth are not well Downloaded by [187.160.112.235] at 22:00 10 June 2016 established in the literature as they demand a complex description including both the nucleation process, which is still largely unknown, and dislocation gliding. In Refs. [122,247] a comprehensive phase-field model is proposed, mimicking the dislocation distribution by a dislocation-density field, eventually evolving in time. Therein, dislocations are observed to play a stabilizing role. In particular, it is shown that the propagation of grooves formed by the growth of the perturbation can be slowed down by a local accumulation of dislocations at their bottom, provided that dislocation mobility is large enough.

Conclusions
In this work, we tried to summarize the present state-of-the-art of continuum simulations of semiconductor heteroepitaxy, focusing on typical growth conditions yielding high-quality materials. The principal drawback of continuum approaches is obvious already from the name: materials are not continuous. However, most atomic-scale processes can be effectively considered, as demonstrated by the successful comparison between many of the approaches discussed here and actual experimental data. One day reliable atomistic approaches might develop to the point where growth can be reliably simulated matching both temporal and spatial laboratory time scales. Quite frankly, we do not envision this to happen in the close future. Even then, we believe that the importance of continuum approaches would not fade. As here above illustrated by the treatment of complex phenomena such as elastic and plastic relaxation, exposure of preferential surface orientations and intermixing, continuum modelling allows for a unique understanding of complex phsyics on relatively simple grounds, which would be hardly achieved by looking at billions of atomic trajectories.
Unavoidably, we focused more on the approaches and on the applications for which we can claim a direct experience. We first listed the set of competing phenomena experimentally revealed in group-IV and III/V semiconductor growth experiments and then analysed how they can be handled, often very successfully, by continuum simulations.
Among the different available methodologies, we find that diffuse-interface approaches (phase-field) are the more versatile, as they allow to combine different physical phenomena in a consistent way. As such, they are good candidates to be the heart of a general-purpose code to simulate different deposition techniques leading to heteroepitaxial growth. It is worth mentioning that phase-field continuum approaches do not necessarily work only at large scales. For instance, the Phase-Field-Crystal method [248][249][250], allows one to directly trace the crystal structure and its evolution. This approach offers several advantages over typical atomistic simulations as it filters out all vibrational dynamics while retaining atomic scale interactions yielding both elastic and even plastic deformations. However, this strength also restricts the method to relatively small scales so that we did not explicitly analyse it in the review.
Faster codes are always welcome, and computational research on efficient algorithms is going to be fundamental for further developments. Still, we believe that the major effort should be made in modelling dislocation evolution during growth. Precise defects control is nowadays the main issue in developing nextgeneration devices. This implies having information on actual dislocation positioning, on the geometry of the full dislocation line and on possible dislocation reactions. The evolution of complex dislocation networks at the level specified above can be provided by dislocation-dynamics simulations (see, e.g. Refs. [251][252][253][254][255]). A code tackling simultaneously the evolving growth morphology and the evolution of the corresponding dislocation network, considering their mutual influence, is what we envision as the future key step.
We believe that the present need to optimally control heteroepitaxy calls for a more intense synergy between experiments and theory. The parameter space is way too large to proceed by trial and errors. Actually, we wrote this review with the precise aim of stimulating further collaborations among the two communities.

Note
1. Notice that experimental conditions usually consist of constant temperature T and pressure P, so that the Gibbs free energy G is considered. However, under the assumption of very low pressure, work by pressure-volume is negligible so that G coincides with the Helmholtz free energy F.