Thin Film Evolution Equation for a Strained Anisotropic Solid Film on a Deformable Isotropic Substrate

Wondimu Tekalign, Agegnehu Atena

Department of Mathematics, Savannah State University, Savannah, GA, USA.

DOI: 10.4236/jamp.2018.64074                

Abstract

We consider a continuum model for the evolution of an epitaxially-strained dislocation-free anisotropic thin solid film on isotropic deformable substrate in the absence of vapor deposition. By using a thin film approximation we derived a nonlinear evolution equation. We examined the nonlinear evolution equation and found that there is a critical film thickness below which every film thickness is stable and a critical wave number above which every film thickness is stable.

Keywords

Thin Films, Evolution Equation, Anisotropy, Hooke’s Law

1. Introduction

Spontaneously formed periodic domain structures of nanoscale islands (quantum dots) in epitaxially strained thin solid films have become a subject of intense theoretical and experimental study. These islands have unique, optical, electronic and magnetic properties which signify their importance in quantum dot applications [1] [2] . These islands are small (a few nanometers) in size and hence difficult to prepare by standard lithographic techniques. One promising way is the formation of islands by a Stranski-Krastanow growth process whereby the planar film undergoes a morphological instability [3] [4] [5] [6] . During heteroepitaxial growth, the instability of surfaces under strain and subsequent island formation is caused by the competition between the surface free energy and the strain energy of the system [7] [8] [9] [10] [11] .

In experimental study of the nonlinear evolution of the stress-driven instability of thick films the formation of deep, cusp-like grooves was observed [12] [13] . In [14] this instability was studied numerically and they showed that the surface instability creates a groove that sharpens as it grows deeper. In [15] a fully nonlinear bifurcation analysis was performed and tracking the branch of steady state solutions numerically they found that the steady state solution branch terminates as the solutions form a cusp singularity.

In a thin film, however, cusp formation is suppressed as the surface approaches the film substrate interface. The different stress fields and the different surface energies of the film and the substrate affect the surface morphology and the film-substrate interface is prevented from being exposed when the wetting criterion is satisfied. Stranski-Krastanow islands will be formed in this case [16] . The steady states of island shapes were studied by many researchers [17] - [22] .

While the understanding of some of the theoretical and modeling issues is well developed, the implementation of the models as large-scale numerical simulations is not yet feasible. The central issue is that the dynamics of the surface morphology is coupled to the elastic strain in the system, so dynamic models require solving the elasticity problems throughout the film and substrate at each time step and are limited by storage limitations for 3-dimensional problem. Simulations involving the full elasticity problem are limited to one or few islands [23] [24] . Only one recent work [25] explored a large number of islands using a large-scale 3-dimensional calculation. This work uses the same elastic constants for both the film and the substrate. In [26] the evolution of a large number of islands was obtained using small slope approximation on a rigid substrate.

In [27] we developed an approximate solution to the elasticity problem which is valid when the film is thin. This elasticity solution is valid for arbitrary elastic constants in the film and substrate. The resulting elasticity solution then removes the necessity for solving the full 3-dimensional elasticity problem numerically, and may provide a means for implementing large-scale simulations. Our work here is to include anisotropic properties of the film to the evolution equation to study the formation of islands. Within this framework, in [28] a non-linear evolution equation with a second-order approximation for the stress field and a nonlinear wetting potential for the interface was considered and it was claimed that the combined effect of nonlinear stress and wetting can terminate the coarsening process and lead to the formation of arrays of equal-sized islands. And, [29] , found that wetting interaction can damp the long-wave perturbations and lead to Turing-type instability, further a weakly nonlinear analysis showed a possibility for spatially periodic arrays of quantum dots which are unstable.

The rest of the paper is organized as follows. In Section 2 we present the full nonlinear model for morphological evolution in thin solid films. In Sections 3 and 4 we describe the thin-film scalings and derive a systematic approximation to find the dominant terms in thin film evolution. In Section 5 we analyze the stability properties of this reduced equation and finally in Section 6 we summarize our results.

2. Model for Film Growth

We begin with the model from [6] . This model describes the evolution of a strained film due to surface diffusion in response to the driving forces of elastic strain energy and surface energy. We consider here the simple case of annealing of a film (no vapor deposition). The film lies in $0<z<h\left(x,y,t\right)$ , the vapor in $z>h\left(x,y,t\right)$ , and the substrate occupies $z<0$ . The vapor is taken to be at zero pressure The states of stress and strain in a deformed crystal being idealized as a continuum are characterized by symmetric second-rank tensors ${\sigma }_{ij}$ and $E_{ij}$ , respectively, each comprising six independent components. Hooke’s law of linear elasticity for the most general anisotropic solid expresses each component of the stress tensor linearly in terms of all components of the strain tensor in the form

${\sigma }_{ij}=c_{ijkl}{E}_{ij}$ (1)

where $c_{ijkl}$ is the array of elastic stiffness constants. Each of the nine equations for a stress component involves nine material parameters. The fourth-order tensor $c_{ijkl}$ comprises 81 components. The symmetry of the stress and strain tensors further imply that the components of the stiffness tensor must satisfy $c_{ijkl}=c_{ijlk}=c_{jikl}$ . As a consequence, the number of independent elastic constants is reduced from 81 to 36. Using this fact and rewriting the stress and strain using their symmetric properties, 1 can be written in a simpler format

${\sigma }_i=c_{ij}{E}_i.$ (2)

Cubic symmetry is a property of crystals that possess three fourfold axes of rotational symmetry, the cube axes, and four threefold axes of rotational symmetry, the cube diagonals. Alternatively, cubic symmetry may be described as invariance of material structure under a translation of a certain distance in any of three mutually orthogonal directions; these directions are usually identified as the cube axes. Consider a cubic material for which the [100], [010] and [001] cube axes are parallel to the axes of an underlying rectangular x; y; z-coordinate system. For this case, it is evident that

$c_{11}=c_{22}=c_{33};c_{12}=c_{23}=c_{31};c_{44}=c_{55}=c_{66}$ (3)

All the other elastic constants vanish because of the fourfold rotational symmetry of the reference axes. Hence elastic response of any cubic crystal is characterized by three independent elastic constants and the stress-strain relationship is given by

${\sigma }_{11}=c_{11}{E}_{11}+c_{12}{E}_{22}+c_{12}{E}_{33}$ (4)

${\sigma }_{22}=c_{12}{E}_{11}+c_{11}{E}_{22}+c_{12}{E}_{33}$ (5)

${\sigma }_{33}=c_{12}{E}_{11}+c_{12}{E}_{22}+c_{11}{E}_{33}$ (6)

${\sigma }_{23}=2c_{44}{E}_{23}$ (7)

${\sigma }_{31}=2c_{44}{E}_{31}$ (8)

${\sigma }_{12}=2c_{44}{E}_{12}$ (9)

where $c_{11}\mathrm{,\; <mi>\; </mi>}{c}_{12}$ and $c_{44}$ are the elastic stiffnesses of the material and

$E_{ij}=\frac{1}{2}\left({\partial }_i{u}_j+{\partial }_j{u}_i\right).$ (10)

Here $u_i$ is the i^th Cartesian component of the displacement vector where the index $i=1,2,3$ corresponds to the $x\mathrm{,}y\mathrm{,}z$ coordinates respectively, and ${\partial }_j$ indicates partial derivative with respect to the j^th coordinate. The quantities $u_i\mathrm{,}{\sigma }_{ij}\mathrm{,}{E}_{ij}$ are defined separately in the film (F) and substrate (S). Since mechanical equilibrium exists within the film and the substrate,

${\partial }_j{\sigma }_{ij}=0\text{in}\text{F},\text{S}.$ (11)

Upon substituting the formulas for stress and strain in these equations we obtain Navier’s equations for the equilibrium displacements, which is valid both in the substrate and film:

$c_{11}\frac{{\partial }^2{u}_1}{\partial x^2}+\left(c_{12}+c_{44}\right)\frac{{\partial }^2{u}_2}{\partial x\partial y}+c_{44}\frac{{\partial }^2{u}_1}{\partial y^2}+\left(c_{12}+c_{44}\right)\frac{{\partial }^2{u}_3}{\partial x\partial z}+c_{44}\frac{{\partial }^2{u}_1}{\partial z^2}=0$ (12)

$c_{44}\frac{{\partial }^2{u}_2}{\partial x^2}+\left(c_{12}+c_{44}\right)\frac{{\partial }^2{u}_1}{\partial x\partial y}+c_{11}\frac{{\partial }^2{u}_2}{\partial y^2}+\left(c_{12}+c_{44}\right)\frac{{\partial }^2{u}_3}{\partial y\partial z}+c_{44}\frac{{\partial }^2{u}_2}{\partial z^2}=0$ (13)

$c_{44}\frac{{\partial }^2{u}_3}{\partial x^2}+\left(c_{12}+c_{44}\right)\frac{{\partial }^2{u}_1}{\partial x\partial z}+c_{44}\frac{{\partial }^2{u}_3}{\partial y^2}+\left(c_{12}+c_{44}\right)\frac{{\partial }^2{u}_2}{\partial y\partial z}+c_{11}\frac{{\partial }^2{u}_3}{\partial z^2}=0$ (14)

The stress balance boundary conditions at the film free surface,

${\sigma }_{ij}^F{\hat{n}}_j+\frac{\partial \gamma }{\partial h}{\delta }_{i3}=0\text{on}z=h\left(x,y,t\right)$ (15)

where

$\hat{n}=\frac{\left(-h_x,-h_y,1\right)}{\sqrt{1+h_x^2+h_y^2}}$ (16)

is the unit normal to the film surface and at the film-substrate interface read

${\sigma }_{ij}^f{\hat{n}}_j-{\sigma }_{ij}^s{\hat{n}}_j+\frac{\partial \gamma }{\partial h}{\delta }_{i3}=\mathrm{0\; <mi>\; </mi>\; <mi>\; </mi>}\text{on}\mathrm{<mi>\; </mi>}z=0$ (17)

where $\partial \gamma /\partial h$ is the wetting stress.

The substrate is taken to be semi-infinite, and so the strains vanish far beneath the film,

$E_{ij}^S\to 0\text{as}z\to -\infty .$ (18)

Finally on $z=0$ (the film/substrate interface), continuity of displacement taking into account the lattice mismatch ò is

$u_i^F=u_i^S+ϵ\left[\begin{array}{c}x\\ y\\ 0\end{array}\right].$ (19)

The evolution equation is given by surface diffusion in response to a chemical potential μ,

$\frac{\partial h}{\partial t}=D\sqrt{1+{\left|\nabla h\right|}^2}{\nabla }_S^2\left(\mu \right)$ (20)

where ${\nabla }_S^2$ the surface Laplacian,

$\begin{array}{c}{\nabla }_S^2=\frac{1}{1+h_x^2+h_y^2}[\left(1+h_y^2\right){\partial }_x^2-2h_x{h}_y{\partial }_x{\partial }_y+\left(1+h_x^2\right){\partial }_y^2\\ -\frac{\left(1+h_y^2\right)h_{xx}-2h_x{h}_y{h}_{xy}+\left(1+h_x^2\right)h_{yy}}{1+h_x^2+h_y^2}\left(h_x{\partial }_x+h_y{\partial }_y\right)]\end{array}$ (21)

D is a constant related to surface diffusion, and the surface chemical potential is

$\mu =\mathcal{E}+\gamma \kappa +\omega \left(h\right)$ (22)

where E is the elastic energy density, $\gamma \kappa$ represents the surface energy, and $\omega \left(h\right)$ is the wetting energy. In the above,

$\mathcal{E}=\frac{1}{2}{\sigma }_{ij}^F{E}_{ij}^F\text{on}z=h\left(x,y,t\right)$ (23)

and the curvature of the film κ is given by

$\kappa =-\frac{\left(1+h_y^2\right)h_{xx}-2h_x{h}_y{h}_{xy}+\left(1+h_x^2\right)h_{yy}}{{\left(1+h_x^2+h_y^2\right)}^{3/2}}.$ (24)

For the wetting energy $\omega \left(h\right)$ we use the two-layer wetting model where the surface energy depends on the film thickness according to

$\gamma \left(h\right)={\gamma }_F+\left({\gamma }_S-{\gamma }_F\right){\text{e}}^{-h/{\delta }_w}$ (25)

The model for the wetting energy $\omega \left(h\right)$ is from [30] , based on a surface energy which depends on the film thickness and undergoes a rapid transition from ${\gamma }_F$ to ${\gamma }_S$ over a length scale δ:

$\gamma \left(h\right)=\frac{1}{2}\left({\gamma }_F+{\gamma }_S\right)+\frac{1}{2}\left({\gamma }_F-{\gamma }_S\right)f\left(h/\delta \right),$ (26)

which gives the wetting term,

$\omega \left(h\right)=n_y{\gamma }^{\prime }\left(h\right)$ (27)

where

$f\left(h/\delta \right)=\frac{2}{\text{π}}\arctan \left(h/\delta \right)$ (28)

and

$\Delta \gamma ={\gamma }_S-{\gamma }_F.$ (29)

Hence we get

$\omega \left(h\right)=\frac{-1}{\sqrt{1+{\left|\nabla h\right|}^2}}\frac{\Delta \gamma }{\text{π}}\frac{\delta }{{\delta }^2+h^2}.$ (30)

Equation (20) is a non-linear moving boundary problem coupled to partial differential equations for the elasticity problem (12)-(19).

3. Steady State Solutions

The governing equations in Section 2 describe the stress state and surface evolution of an epitaxially strained film. They have a basic-state solution corresponding to a completely relaxed, stress-free substrate,

${\bar{u}}_i^s=0,\mathrm{<mi>\; </mi>\; <mi>\; </mi>}{\bar{\sigma }}_{ij}^s=\mathrm{0\; <mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>}\text{for}i,j=1,2,3$ (31)

and a planar film with spatially uniform stress and strain,

${\bar{u}}_1^f=ϵx,\mathrm{<mi>\; </mi>\; <mi>\; </mi>}{\bar{u}}_2^f=ϵy,\mathrm{<mi>\; </mi>\; <mi>\; </mi>}{\bar{u}}_3^f=\frac{1}{c_{11}}\left(-2c_{12}ϵ-\frac{\partial \gamma }{\partial h}\right)z$ (32)

${\bar{\sigma }}_{11}^f={\bar{\sigma }}_{22}^f=M_{001}ϵ-\frac{c_{12}}{c_{11}}\frac{\partial \gamma }{\partial h},\mathrm{<mi>\; </mi>\; <mi>\; </mi>}{\bar{\sigma }}_{33}^f=-\frac{\partial \gamma }{\partial h}$ (33)

where

$M_{001}=\left(c_{11}+c_{12}-2\frac{c_{12}^2}{c_{11}}\right)$ (34)

which is the biaxial modulus in the plane with normal in the (001) material direction (in our case (001) material direction is the direction parallel to the z-axis).

The total elastic energy store in the film due to epitaxial and wetting stresses is

${\mathcal{E}}_0=\frac{1}{2}{\bar{\sigma }}_{ij}^f{\bar{E}}_{ij}^f=2{ϵ}^2{M}_{001}+\frac{1}{c_{11}}{\left(\frac{\partial \gamma }{\partial h}\right)}^2$ (35)

In Section 5 we perform a linear stability analysis of this basic state of an epitaxial film.

4. Nondimensional Anisotropic Governing Equations

Here we derive the evolution equation based on approximation that wavelength of surface undulations is large compared to the characteristic film thickness $H_0$ . Define

$\alpha =\frac{H_0}{l}\ll 1$ (36)

where l is characteristic length scale in $\left(x\mathrm{,}y\right)$ . Let us use the following scalings:

$\begin{array}{l}h=\alpha Hl\\ x=lX\\ y=lY\\ z=\alpha lZ\\ t=\tau T\\ u_i\left(x,y,z\right)=l{U}_i\left(X,Y,Z\right)\\ \delta =l\tilde{\delta }\end{array}\}$ (37)

Next let us write the elastic stresses, strains and elastic energy density as expansions in α:

${\sigma }_{ij}={\left({\sigma }_{ij}\right)}_0+\alpha {\left({\sigma }_{ij}\right)}_1+{\alpha }^2{\left({\sigma }_{ij}\right)}_2+\cdots$ (38)

$E_{ij}={\left(E_{ij}\right)}_0+\alpha {\left(E_{ij}\right)}_1+{\alpha }^2{\left(E_{ij}\right)}_2+\cdots$ (39)

$\mathcal{E}={\mathcal{E}}_0+\alpha {\mathcal{E}}_1+{\alpha }^2{\mathcal{E}}_2+\cdots$ (40)

Expanding the other quantities in (9) we obtain

$\begin{array}{c}\kappa =-\frac{1}{l}(\alpha \left(H_{XX}+H_{YY}\right)+{\alpha }^3(H_Y^2{H}_{XX}-2H_X{H}_Y{H}_{XY}+H_X^2{H}_{YY}\\ -3/2{\left|\nabla H\right|}^2\left(H_{XX}+H_{YY}\right)))+\mathcal{O}\left({\alpha }^5\right)\end{array}$ (41)

and the surface Laplacian has the expansion,

$\begin{array}{c}{\nabla }_S^2=\frac{1}{l^2}[\left({\partial }_X^2+{\partial }_Y^2\right)+{\alpha }^2(H_Y^2{\partial }_X^2-2H_X{H}_Y{\partial }_X{\partial }_Y+H_X^2{\partial }_Y^2\\ +\left(H_{XX}+H_{YY}\right)\left(H_X{\partial }_X+H_Y{\partial }_Y\right)-{\left|\nabla H\right|}^2\left({\partial }_X^2+{\partial }_Y^2\right))]+\mathcal{O}\left({\alpha }^4\right).\end{array}$ (42)

We choose the time scale τ as

$\tau =\frac{l^4}{D\gamma }$ (43)

and length scale l as

$l=\frac{\gamma }{{\mathcal{E}}_0}.$ (44)

We also assume that the characteristic film thickness is much larger than the wetting layer thickness, so

$\alpha \gg \tilde{\delta }\mathrm{.}$ (45)

Hence (30) can be written as:

$\omega \left(\alpha lH\right)~\frac{-\Delta \gamma }{l\text{π}}\frac{\tilde{\delta }}{{\alpha }^2{H}^2}\left(1-\frac{1}{2}{\alpha }^2{\left|\nabla H\right|}^2-\frac{{\tilde{\delta }}^2}{{\alpha }^2{H}^2}+\cdots \right).$ (46)

To balance the wetting energy term $\mathcal{O}\left(\tilde{\delta }/{\alpha }^2\right)$ with the surface energy term $\mathcal{O}\left(\alpha \right)$ in (20) we choose

$\tilde{\delta }=\mathcal{O}\left({\alpha }^3\right)$ (47)

and define

$\tilde{\delta }={\delta }^{*}{\alpha }^3$ (48)

where ${\delta }^{\mathrm{<mo>\; *\; </mo>}}=\mathcal{O}\left(1\right)$ . Then defining

$\tilde{\omega }=\frac{1}{\alpha }\frac{l}{\gamma }\omega \left(\alpha lH\right),$ (49)

we obtain

$\tilde{\omega }={\tilde{\omega }}_0+{\alpha }^2{\tilde{\omega }}_2+\mathcal{O}\left({\alpha }^4\right)$ (50)

where

${\tilde{\omega }}_0=\frac{-\Delta \gamma }{\text{π}\gamma }\frac{{\delta }^{*}}{H^2}$ (51)

${\tilde{\omega }}_2=\frac{1}{2}\frac{-\Delta \gamma }{\text{π}\gamma }\frac{{\delta }^{*}}{H^2}{\left|\nabla H\right|}^2.$ (52)

We substitute these expansions into the evolution Equation (20) to obtain

$\frac{\partial H}{\partial T}=\left({\nabla }^2\right)\left(E_1-{\nabla }^{2}H+{\tilde{\omega }}_0\right)+\alpha \left({\nabla }^2\right)E_2+\mathcal{O}\left({\alpha }^2\right)$ (53)

where

${\tilde{\mathcal{E}}}_1=\frac{l\mathcal{E}}{\gamma }=\frac{{\mathcal{E}}_1}{{\mathcal{E}}_0}$ (54)

and

${\tilde{\mathcal{E}}}_2=\frac{l{\mathcal{E}}_2}{\gamma }=\frac{{\mathcal{E}}_2}{{\mathcal{E}}_0}.$ (55)

Equation (53) is the thin film evolution equation. It depends on elastic response through ${\tilde{\mathcal{E}}}_1$ and ${\tilde{\mathcal{E}}}_2$ , which we now determine.

4.1. Elastic Response of Film

Now we find ${\tilde{\mathcal{E}}}_1$ by solving the elasticity problem. Using scalings given in (37), (12)-(14) can be written as:

$\begin{array}{l}{c}_{11}\frac{{\partial }^2{U}_1}{\partial X^2}+c_{44}\frac{{\partial }^2{U}_1}{\partial Y^2}+c_{44}\frac{1}{{\alpha }^2}\frac{{\partial }^2{U}_1}{\partial Z^2}+\left(c_{12}+c_{44}\right)\frac{{\partial }^2{U}_2}{\partial X\partial Y}\\ +\left(c_{12}+c_{44}\right)\frac{1}{\alpha }\frac{{\partial }^2{U}_3}{\partial X\partial Z}=0,\end{array}$ (56)

$\begin{array}{l}{c}_{44}\frac{{\partial }^2{U}_2}{\partial X^2}+c_{11}\frac{{\partial }^2{U}_2}{\partial Y^2}+c_{44}\frac{1}{{\alpha }^2}\frac{{\partial }^2{U}_2}{\partial Z^2}+\left(c_{12}+c_{44}\right)\frac{{\partial }^2{U}_1}{\partial X\partial Y}\\ +\left(c_{12}+c_{44}\right)\frac{1}{\alpha }\frac{{\partial }^2{U}_3}{\partial Y\partial Z}=0,\end{array}$ (57)

$\begin{array}{l}{c}_{44}\frac{{\partial }^2{U}_3}{\partial X^2}+c_{44}\frac{{\partial }^2{U}_3}{\partial Y^2}+c_{11}\frac{1}{{\alpha }^2}\frac{{\partial }^2{U}_3}{\partial Z^2}+\left(c_{12}+c_{44}\right)\frac{1}{\alpha }\frac{{\partial }^2{U}_1}{\partial X\partial Z}\\ +\left(c_{12}+c_{44}\right)\frac{1}{\alpha }\frac{{\partial }^2{U}_2}{\partial Y\partial Z}=0.\end{array}$ (58)

The same way we can write stress and strain using (37) and from the boundary condition (15) we obtain:

$\begin{array}{l}\alpha \left(c_{11}\frac{\partial U_1}{\partial X}+c_{12}\frac{\partial U_2}{\partial Y}+\frac{1}{\alpha }{c}_{12}\frac{\partial U_3}{\partial Z}\right)H_X+\alpha c_{44}\left(\frac{\partial U_1}{\partial Y}+\frac{\partial U_2}{\partial X}\right)H_Y\\ -c_{44}\left(\frac{1}{\alpha }\frac{\partial U_1}{\partial Z}+\frac{\partial U_3}{\partial X}\right)=0,\end{array}$ (59)

$\begin{array}{l}\alpha c_{44}\left(\frac{\partial U_1}{\partial Y}+\frac{\partial U_2}{\partial X}\right)H_X+\alpha \left(c_{12}\frac{\partial U_1}{\partial X}+c_{11}\frac{\partial U_2}{\partial Y}+c_{12}\frac{1}{\alpha }\frac{\partial U_3}{\partial Z}\right)H_Y\\ -c_{44}\left(\frac{1}{\alpha }\frac{\partial U_2}{\partial Z}+\frac{\partial U_3}{\partial Y}\right)=0,\end{array}$ (60)

$\begin{array}{l}\alpha c_{44}\left(\frac{1}{\alpha }\frac{\partial U_1}{\partial Z}+\frac{\partial U_3}{\partial X}\right)H_X+\alpha c_{44}\left(\frac{1}{\alpha }\frac{\partial U_2}{\partial Z}+\frac{\partial U_3}{\partial Y}\right)H_Y\\ -\left(c_{12}\frac{\partial U_1}{\partial X}+c_{12}\frac{\partial U_2}{\partial Y}+\frac{1}{\alpha }{c}_{11}\frac{\partial U_3}{\partial Z}\right)=0\end{array}$ (61)

on $Z=H\left(X,Y,T\right)$ .

We proceed by finding elasticity solutions in the film satisfying the boundary condition on $Z=H\left(X,Y,T\right)$ . These solutions have unknown constants to be determined from the boundary conditions on $z=0$ . Then we construct elasticity solutions in the substrate to determine these constants. Finally we substitute the final elasticity solutions into evolution equation.

Now let us expand the displacements U in a,

$U_i=U_{i0}+\alpha U_{i1}+{\alpha }^2{U}_{i2}+\cdots$ (62)

We substitute (62) in (56)-(58) and compare by order in a.

At $\mathcal{O}\left(1\right)$ we obtain:

$\begin{array}{l}{c}_{44}\frac{{\partial }^2{U}_{10}}{\partial Z^2}=0\hfill \\ c_{44}\frac{{\partial }^2{U}_{20}}{\partial Z^2}=0\hfill \\ c_{11}\frac{{\partial }^2{U}_{30}}{\partial Z^2}=0\hfill \end{array}\}$ (63)

with the boundary conditions,

$\begin{array}{l}-c_{44}\frac{\partial U_{10}}{\partial Z}=0\text{on}Z=H\hfill \\ -c_{44}\frac{\partial U_{20}}{\partial Z}=0\text{on}Z=H\hfill \\ -c_{11}\frac{\partial U_{30}}{\partial Z}=0\text{on}Z=H\hfill \end{array}\}$ (64)

We expect the $\mathcal{O}\left(1\right)$ strain tensor is

$E^F=\left[\begin{array}{ccc}ϵ& 0& 0\\ 0& ϵ& 0\\ 0& 0& -2ϵc_{12}/c_{11}\end{array}\right]$ (65)

and solving (63) and using (64) we obtain $\mathcal{O}\left(1\right)$ solution:

$\begin{array}{l}{U}_{10}=ϵX\hfill \\ U_{20}=ϵY\hfill \\ U_{30}=0\hfill \end{array}\}$ (66)

Solving the $\mathcal{O}\left(\alpha \right)$ differential equations we obtain

$U_{i1}=B_{i1}\left(X,Y\right)Z+A_{i1}\left(X,Y\right),i=1,2,3.$ (67)

Using $\mathcal{O}\left(\alpha \right)$ boundary conditions,

$B_{11}=0,$ (68)

$B_{21}=0$ (69)

and

$B_{31}=\frac{-2ϵc_{12}}{c_{11}}.$ (70)

The same way using the $\mathcal{O}\left({\alpha }^2\right)$ differential equations we obtain

$U_{i2}=B_{i2}\left(X,Y\right)Z+A_{i2}\left(X,Y\right)\text{for}i=1,2,3.$ (71)

And the $\mathcal{O}\left({\alpha }^2\right)$ boundary conditions gives us,

$B_{12}=\frac{ϵ}{c_{44}}\left(c_{11}+c_{12}-\frac{2c_{12}^2}{c_{11}}\right)H_X-A_{31,X},$ (72)

$B_{22}=\frac{ϵ}{c_{44}}\left(c_{11}+c_{12}-\frac{2c_{12}^2}{c_{11}}\right)H_Y-A_{31,Y},$ (73)

and

$B_{32}=-\frac{c_{12}}{c_{11}}\left(A_{11,X}+A_{21,Y}\right).$ (74)

Returning to our expansion for $\mathcal{E}$ we obtain

${\mathcal{E}}_0={ϵ}^2\left(c_{11}+c_{12}-2\frac{c_{12}^2}{c_{11}}\right)\equiv {ϵ}^2{M}_{001}$ (75)

and

${\mathcal{E}}_1=\frac{ϵ}{c_{11}}\left(c_{11}-c_{12}\right)\left(c_{11}+2c_{12}\right)\left(\frac{\partial A_{11}}{\partial X}+\frac{\partial A_{21}}{\partial Y}\right).$ (76)

4.2. Elastic Response of Substrate

The functions $A_{11}$ and $A_{21}$ appearing in ${\mathcal{E}}_1$ are determined by the response of the substrate. Since the substrate is a semi-infinite domain, we can seek solutions in terms of Fourier transforms in $\left(X\mathrm{,}Y\right)$ . In the substrate let

$\begin{array}{l}\left(x\mathrm{,}y\mathrm{,}z\right)=l\left(X\mathrm{,}Y\mathrm{,}\tilde{Z}\right)\hfill \\ u_i=l{U}_i,i=\mathrm{1,2,3}\hfill \\ a=\sqrt{a_X^2+a_Y^2}\hfill \end{array}\}$ (77)

The elasticity Equations (12)-(14) and far field condition (18) become

$\begin{array}{l}{c}_{11}^S\frac{{\partial }^2{U}_1}{\partial X^2}+c_{44}^S\frac{{\partial }^2{U}_1}{\partial Y^2}+c_{44}^S\frac{{\partial }^2{U}_1}{\partial Z^2}+\left(c_{12}^S+c_{44}^S\right)\frac{{\partial }^2{U}_2}{\partial X\partial Y}\\ +\left(c_{12}^S+c_{44}^S\right)\frac{{\partial }^2{U}_3}{\partial X\partial Z}=0,\end{array}$ (78)

$\begin{array}{l}{c}_{44}^S\frac{{\partial }^2{U}_2}{\partial X^2}+c_{11}^S\frac{{\partial }^2{U}_2}{\partial Y^2}+c_{44}^S\frac{{\partial }^2{U}_2}{\partial Z^2}+\left(c_{12}^S+c_{44}^S\right)\frac{{\partial }^2{U}_1}{\partial X\partial Y}\\ +\left(c_{12}^S+c_{44}^S\right)\frac{{\partial }^2{U}_3}{\partial Y\partial Z}=0,\end{array}$ (79)

$\begin{array}{l}{c}_{44}^S\frac{{\partial }^2{U}_3}{\partial X^2}+c_{44}^S\frac{{\partial }^2{U}_3}{\partial Y^2}+c_{11}^S\frac{{\partial }^2{U}_3}{\partial Z^2}+\left(c_{12}^S+c_{44}^S\right)\frac{{\partial }^2{U}_1}{\partial X\partial Z}\\ +\left(c_{12}^S+c_{44}^S\right)\frac{{\partial }^2{U}_2}{\partial Y\partial Z}=0.\end{array}$ (80)

and

$U_i\to 0\text{as}\tilde{Z}\to -\infty \mathrm{.}$ (81)

Using the Fourier transform to define $\hat{H}\left(a_X\mathrm{,}{a}_Y\right)$ and ${\hat{U}}_i\left(a_X\mathrm{,}{a}_Y\mathrm{,}Z\right)$ we have,

$\begin{array}{l}H={\displaystyle {\int }_{-\infty }^{\infty }}{\displaystyle {\int }_{-\infty }^{\infty }}\hat{H}\exp \left(-i{a}_{X}X-i{a}_{Y}Y\right)\text{d}{a}_X\text{d}{a}_Y\\ U_i^S={\displaystyle {\int }_{-\infty }^{\infty }}{\displaystyle {\int }_{-\infty }^{\infty }}{\hat{U}}_i^S\exp \left(-i{a}_{X}X-i{a}_{Y}Y\right)\text{d}{a}_X\text{d}{a}_Y\\ {\hat{U}}_i^S=\frac{1}{{\left(\text{2π}\right)}^2}{\displaystyle {\int }_{-\infty }^{\infty }}{\displaystyle {\int }_{-\infty }^{\infty }}{U}_i^S\exp \left(i{a}_{X}X+i{a}_{Y}Y\right)\text{d}X\text{d}Y\end{array}\}$ (82)

Then using (82) we write (78)-(80) as

$\begin{array}{l}{c}_{44}^S{\partial }_{\tilde{Z}}^2{\hat{U}}_1-\left(a_X^2{c}_{11}^S+a_Y^2{c}_{44}^S\right){\hat{U}}_1-\left(c_{12}^S+c_{44}^S\right)a_X{a}_Y{\hat{U}}_2\\ -i{a}_X\left(c_{12}^S+c_{44}^S\right){\partial }_{\tilde{Z}}{\hat{U}}_3=\mathrm{0,}\end{array}$ (83)

$\begin{array}{l}{c}_{44}^S{\partial }_{\tilde{Z}}^2{\hat{U}}_2-\left(a_X^2{c}_{44}^S+a_Y^2{c}_{11}^S\right){\hat{U}}_2-\left(c_{12}^S+c_{44}^S\right)a_X{a}_Y{\hat{U}}_1\\ -i{a}_Y\left(c_{12}^S+c_{44}^S\right){\partial }_{\tilde{Z}}{\hat{U}}_3=0,\end{array}$ (84)

and

$c_{11}^S{\partial }_{\tilde{Z}}^2{\hat{U}}_3-a^2{c}_{44}^S{\hat{U}}_3-i\left(c_{12}^S+c_{44}^S\right)\left(a_X{\partial }_{\tilde{Z}}{\hat{U}}_1+a_Y{\partial }_{\tilde{Z}}{\hat{U}}_2\right)=0$ (85)

in $\tilde{Z}<0$ . With the condition that

${\hat{U}}_i^S\to 0\text{as}\tilde{Z}\to -\infty \mathrm{.}$ (86)

We expand displacement in the substrate using powers in a,

$U_i^S=U_{i0}^S+\alpha U_{i1}^S+\cdots$ (87)

with

$U_{ij}^S={\displaystyle {\int }_{-\infty }^{\infty }}{\displaystyle {\int }_{-\infty }^{\infty }}{\hat{U}}_{ij}^S\exp \left(-i{a}_{X}X-i{a}_{Y}Y\right)\text{d}{a}_X\text{d}{a}_Y.$ (88)

Using (87) and the solution for $U_{i0}^F$ we obtain,

$U_{i0}^S=0,i=1,2,3.$ (89)

Using (86), (89) and (19), an $\mathcal{O}\left(\alpha \right)$ comparison gives us

${U_{i1}^F|}_{Z=0}={U_{i1}^S|}_{\tilde{Z}=0}\text{for}i=1,2,3.$ (90)

The $\mathcal{O}\left(\alpha \right)$ solution of the form (82) satisfying (83)-(85) and (86) has the form

$\left[\begin{array}{c}{\hat{U}}_{11}^S\\ {\hat{U}}_{21}^S\\ {\hat{U}}_{31}^S\end{array}\right]=\left[\begin{array}{c}{\alpha }_1\\ {\alpha }_2\\ {\alpha }_3\end{array}\right]{\text{e}}^{a\tilde{Z}}+\left[\begin{array}{c}{\beta }_1\\ {\beta }_2\\ {\beta }_3\end{array}\right]\tilde{Z}{\text{e}}^{a\tilde{Z}}$ (91)

where ${\beta }_i$ is obtained from (83)-(85) in terms of ${\alpha }_i$ and ${\alpha }_i$ is obtained from the boundary conditions at the film/substrate interface. Then for

$A_{i1}={\displaystyle {\int }_{-\infty }^{\infty }}{\displaystyle {\int }_{-\infty }^{\infty }}{\hat{A}}_{i1}\exp \left(-i{a}_{X}X-i{a}_{Y}Y\right)\text{d}{a}_X\text{d}{a}_Y,$ (92)

using (88) and (90)

${\hat{A}}_{i1}={\alpha }_i,i=1,2,3.$ (93)

The second boundary condition is (17). At $\mathcal{O}\left(\alpha \right)$ we obtain

${\left({\sigma }_{31}^F\right)}_1={\left({\sigma }_{31}^S\right)}_1,$ (94)

${\left({\sigma }_{32}^F\right)}_1={\left({\sigma }_{32}^S\right)}_1$ (95)

and

${\left({\sigma }_{33}^F\right)}_1={\left({\sigma }_{33}^S\right)}_1$ (96)

which gives us a linear system of equations for ${\alpha }_1\mathrm{,\; <mi>\; </mi>}{\alpha }_2$ and ${\alpha }_3$ . Solving this system we obtain

$-ϵM_{001}i{a}_X\hat{H}={\mu }^s\left(-i{a}_X{u}_3^1+a{u}_1^1+a_X{\delta }_3^1\right)$ (97)

$-ϵM_{001}i{a}_Y\hat{H}={\mu }^s\left(-i{a}_Y{u}_3^1+a{u}_2^1+a_Y{\delta }_3^1\right)$ (98)

and

${\nu }^S\left(-i{a}_X{U}_1^1-i{a}_Y{U}_2^1\right)+\left(1-{\nu }^S\right)\left(a{U}_3^1+ia{\delta }_3^1\right)=0.$ (99)

Solving (97), (98) and (99) we obtain

$U_1^1=M_{001}\frac{a_X\hat{H}ϵi}{a}\frac{{\nu }^s-1}{{\mu }^s}$ (100)

$U_2^1=M_{001}\frac{a_Y\hat{H}ϵi}{a}\frac{{\nu }^s-1}{{\mu }^s}$ (101)

$U_3^1=M_{001}\frac{\hat{H}ϵ}{2}\frac{2{\nu }^s-1}{{\mu }^s}$ (102)

Hence the Fourier transform ${\hat{\mathcal{E}}}_1=\mathcal{F}\left[{\mathcal{E}}_1\right]$ can be written as

${\hat{\mathcal{E}}}_1=-{\mathcal{E}}_{0}Ea\hat{H}$ (103)

where

$E=M_{001}\frac{1-{\nu }^S}{{\mu }^S}.$ (104)

Hence in nondimensional form,

${\widehat{\stackrel{˜}{\mathcal{E}}}}_1=-Ea\hat{H}$ (105)

Note that the constant E contains the interaction of the elastic response of the anisotropic film and isotropic substrate.

The evolution equation at $\mathcal{O}\left(1\right)$ is thus

$\frac{\partial H}{\partial T}={\nabla }^2\left({\tilde{\mathcal{E}}}_1-{\nabla }^{2}H-\frac{r}{H^2}\right)$ (106)

where

${\mathcal{E}}_1={\mathcal{F}}^{-1}\left[{\hat{\mathcal{E}}}_1\right]={\mathcal{F}}^{-1}\left[-aE\hat{H}\right]$ (107)

with the strength of the wetting energy quantified by

$r=\frac{\Delta \gamma {\delta }^{*}}{\text{π}\gamma }.$ (108)

Equation (106) contains the dominant effects for thin films, derived in a self-consistent thin film approximation. The three terms on the right hand side, correspond to elastic energy ( ${\tilde{\mathcal{E}}}_1$ ), surface energy ( ${\nabla }^{2}H$ ) and wetting energy ( $r/H^2$ ). Of particular note is that in this equation elasticity is linear in H, but remains a non-local term in the evolution equation. In addition, the only nonlinear contribution is from the wetting energy.

5. Linear Stability Analysis

To determine the stability of planar films of thickness $\bar{H}$ we consider normal-mode perturbation with wavenumbers $a_X$ and $a_Y$ :

$H=\bar{H}+\hat{H}\exp \left(\sigma T+i{a}_{X}X+i{a}_{Y}Y\right).$ (109)

We consider $\hat{H}\ll \bar{H}$ , substitute (109) into (106), and linearize to obtain the characteristic equation

$\sigma =-a^4+E{a}^3-\frac{2r}{{\bar{H}}^3}{a}^2.$ (110)

When the film wets the substrate, $r>0$ and a typical graph of (110) is as shown in Figure 1.

In Figure 2 we plot the neutral stability condition $\sigma =0$ for the film thickness $\bar{H}$ versus a. We see that there is a critical film thickness, ${\bar{H}}_c$ below

Figure 1. Linear stability growth rate σ versus the wave number a.

Figure 2. Schematic of the stability diagram for the film thickness $\bar{H}$ versus perturbation wavenumber a.

which every film thickness is stable,

${\bar{H}}_c={\left(\frac{8r}{E^2}\right)}^{1/3}$ (111)

and there is a critical wave number, $a_c$ above which every film thickness is stable,

$a_c=E.$ (112)

The evolution equation thus has the property that sufficiently thin films are stabilized by the wetting effect, but thicker films are unstable to the stress driven morphological instability.

6. Summary

We derived a self-contained evolution equation where the film thickness is smaller than wavelength of surface variations. Our evolution equation includes effects of anisotropic elastic constants for cubic symmetry in the film and isotropic elastic constants in the substrate, isotropic surface energy, and wetting energy. This evolution equation possesses steady state solutions corresponding to island formation, and is a possible candidate for use in large scale simulations of island systems.

Acknowledgements

Sincere thanks to the members of JAMP for their professional performance. And W. Tekalign would like to thank his Ph.D. adviser Dr. Brian J. Spencer for his valuable inputs and support throughout the years.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Kastner, M.A. (1996) Mesoscopic Physics with Artificial Atoms. Proceedings of the 23rd International Conference on Physics of Semiconductors, 1, 27.
[2]	Chen, M. and Porod, W. (1995) Design of Gate-Confined Quantum-Dot Structures in the Few-Electron Regime. Journal of Applied Physics, 78, 1050. https://doi.org/10.1063/1.360339
[3]	Asaro, R.J. and Tiller, W.A. (1972) Interface Morphology Development during Stress Corrosion Cracking: Part I. via Surface Diffusion. Metallurgical and Materials Transactions B, 3, 1789-1796. https://doi.org/10.1007/BF02642562
[4]	Grinfeld, M.A. (1986) Instability of the Interface between a Non-Hydrostatically Stressed Elastic Body and a Melt. Doklady Akademii Nauk SSSR Sov.Phys. Dokl. 31, 831.
[5]	Srolovitz, D.J. (1989) On the Stability of Surfaces of Stressed Solids. Acta Metallurgica, 37, 621-625. https://doi.org/10.1016/0001-6160(89)90246-0
[6]	Spencer, B.J., Voorhees, P.W. and Davis, S.H. (1991) Morphological Instability in Epitaxially Strained Dislocation-Free Solid Films. Physical Review Letters, 67, 3696. https://doi.org/10.1103/PhysRevLett.67.3696
[7]	Daruka, I., Tersoff, J. and Barbasi, A.L. (1999) Shape Transition in Growth of Strained Islands. Physical Review Letters, 82, 2753. https://doi.org/10.1103/PhysRevLett.82.2753
[8]	Nix, W.D. (1989) Mechanical Properties of Thin Films. Metallurgical and Materials Transactions A, 20A, 2217. https://doi.org/10.1007/BF02666659
[9]	Gao, H. (1991) Stress Concentration at Slightly Undulating Surfaces. Journal of the Mechanics and Physics of Solids, 39, 443-458. https://doi.org/10.1016/0022-5096(91)90035-M
[10]	Gao, H. (1991) A Boundary Perturbation Analysis for Elastic Inclusions and Interfaces. International Journal of Solids and Structures, 28, 703-725. https://doi.org/10.1016/0020-7683(91)90151-5
[11]	Grinfeld, M.A. (1993) The Stress Driven Instability in Elastic Crystals: Mathematical Models and Physical Manifestations. Journal of Nonlinear Science, 3, 35-83. https://doi.org/10.1007/BF02429859
[12]	Berrehar, J., Caroli, C., Lapersonne-Meyer, C. and Schott, M. (1992) Surface Patterns on Single-Crystal Films under Uniaxial Stress: Experimental Evidence for the Grinfeld Instability. Physical Review B, 46, 13487. https://doi.org/10.1103/PhysRevB.46.13487
[13]	Jesson, D.E., Pennycook, S.J., Baribeau, J.M. and Houghton, D.C. (1993) Direct Imaging of Surface Cusp Evolution during Strained-Layer Epitaxy and Implications for Strain Relaxation. Physical Review Letters, 71, 1744. https://doi.org/10.1103/PhysRevLett.71.1744
[14]	Yang, W.H. and Srolovitz, D.J. (1994) Surface Morphology Evolution in Stressed Solids: Surface Diffusion Controlled Crack Initiation. Journal of the Mechanics and Physics of Solids, 42, 1551-1574. https://doi.org/10.1016/0022-5096(94)90087-6
[15]	Spencer, B.J. and Meiron, D.I. (1994) Nonlinear Evolution of the Stress Driven Morphological Instability in a Two Dimensional Semi-Infinite Solid. Acta Metallurgica et Materialia, 37, 621.
[16]	Xiang, Y. and Weinan, E. (2002) Nonlinear Evolution Equation for the Stress-Driven Morphological Instability. Journal of Applied Physics, 91, 9414. https://doi.org/10.1063/1.1477259
[17]	Spencer, B.J. and Tersoff, J. (1996) Equilibrium Shapes of Small Strained Islands. Materials Research Society Symposia Proceedings, 399, 283-288. https://doi.org/10.1557/PROC-399-283
[18]	Spencer, B.J. and Tersoff, J. (1997) Equilibrium Shapes and Properties of Epitaxially Strained Islands. Physical Review Letters, 79, 4858-4861. https://doi.org/10.1103/PhysRevLett.79.4858
[19]	Kukta, R.V. and Freund, L.B. (1997) Minimum Energy Configuration of Epitaxial Material Clusters on a Lattice-Mismatched Substrate. Journal of the Mechanics and Physics of Solids, 45, 1835-1860. https://doi.org/10.1016/S0022-5096(97)00031-8
[20]	Rudin, C.D. and Spencer, B.J. (1999) Equilibrium Island Ridge Arrays in Strained Solid Films. Journal of Applied Physics, 86, 5530-5536. https://doi.org/10.1063/1.371556
[21]	Shanahan, L.L. and Spencer, B.J. (2002) Interfaces and Free Boundaries. Interfaces and Free Boundaries, 4, 1-25. https://doi.org/10.4171/IFB/50
[22]	Tekalign, W.T. and Spencer, B.J. (2007) Thin-Film Evolution Equation for a Strained Solid Film on a Deformable Substrate: Numerical Steady States. Journal of Applied Physics, 102, Article ID: 073503. https://doi.org/10.1063/1.2785024
[23]	Zhang, Y.W. and Bower, A.F. (1999) Numerical Simulations of Island Formation in a Coherent Strained Epitaxial Thin Film System. Journal of the Mechanics and Physics of Solids, 47, 22730-2297. https://doi.org/10.1016/S0022-5096(99)00026-5
[24]	Bergamaschini, R., Salvalaglio, M., Backofen, R., Voigt, A. and Montalenti, F. (2016) Continuum Modelling of Semiconductor Heteroepitaxy: An Applied Perspective. Advances in Physics, 1, 331-367.
[25]	Liu, P., Zhang, Y.W. and Lu, C. (2003) Coarsening Kinetics of Heteroepitaxial Islands in Nucleationless Stranski-Krastanov Growth. Physical Review B, 68, Article ID: 035402. https://doi.org/10.1103/physrevb.68.035402
[26]	Golovin, A.A., Davis, S.H. and Voorhees, P.W. (2003) Self-Organization of Quantum Dots in Epitaxially Strained Solid Films. Physical Review E, 68, Article ID: 056203. https://doi.org/10.1103/PhysRevE.68.056203
[27]	Tekalign, W.T. and Spencer, B.J. (2004) Evolution Equation for a Thin Epitaxial Film on a Deformable Substrate. Journal of Applied Physics, 96, 5505-5512. https://doi.org/10.1063/1.1766084
[28]	Pang, Y. and Huang, R. (2006) Nonlinear Effect of Stress and Wetting on Surface Evolution of Epitaxial Thin Films. Physical Review B, 74, Article ID: 075413. https://doi.org/10.1103/PhysRevB.74.075413
[29]	Levine, M.S., Golovin, A.A., Davis, S.H. and Voorhees, P.W. (2007) Self-Assembly of Quantum Dots in a Thin Epitaxial Film Wetting an Elastic Substrate. Physical Review B, 75, Article ID: 205312. https://doi.org/10.1103/PhysRevB.75.205312
[30]	Spencer, B.J. (1999) Asymptotic Derivation of the Glued-Wetting-Layer Model and Contact-Angle Condition for Stranski-Krastanow Islands. Physical Review B, 59, 2011. https://doi.org/10.1103/PhysRevB.59.2011