format: text

text:
# A continuum description of the energetics and evolution of stepped surfaces in strained nanostructures

## 1 Introduction

Self-organized semiconductor nanostructures hold the promise for manufacture of micro-electronic devices with unprecedented performance characteristics. Because of the close connection between size and electronic characteristics, these structures can be tuned to very specific requirements. Potential applications include field effect transistors, quantum memory devices, and solid-state lasers.

Strain-driven nucleation, growth and coarsening of epitaxial islands offers a versatile approach to manufacturing nanoscale devices. To exploit this phenomenon in manufacturing, a fundamental understanding of the role of strain, surface energies, and the kinetics of transport on the formation and evolution of material structures is required.

Surface energy is an important concept in considering the evolution of microstructure in small-scale material systems. The free energy of the bounding surface of a crystal is a macroscopic quantity representing, in some sense, the net work that had to be done to create that surface. This energy is distributed over a specific mathematical surface that approximates the physical boundary between the crystal and its surroundings. From this definition of surface energy it is clear that the reference level for energy of a free surface is the state of the material on that same crystallographic surface when it is embedded deep within a perfect crystal. A number of analyses of surface stability and evolution have been reported. Asaro and Tiller (1972), Grinfeld (1986) and Srolovitz (1989) independently showed that the flat surface of a stressed solid under two-dimensional plane strain conditions is unstable if the surface energy is independent of surface orientation. The same physical model was extended to three dimensions, but restricted to small amplitude surface fluctuations, by Freund (1995). A numerical method for handling large amplitude surface fluctuations was developed by Zhang and Bower (1998). Gao and Nix (1999) adopted the same physical framework for describing the breakup of an unstable film into islands. In all of these studies, the surface energy was assumed to be isotropic, independent of surface strain and to have no connection to the discrete nature of the crystalline materials being modeled. This study is focused on a departure from this point of view.

The value of surface energy per unit area of a given crystallographic surface orientation is determined by the fine scale structure of that surface. For a high symmetry orientation in a crystal, such as a {100} surface of a cubic crystal, the surface is atomically flat. For other orientations close to this surface, the structure usually consists of flat terraces with well defined local surface energies, separated by atomic scale ledges or steps as illustrated schematically in Fig. 1. The steps alter the macroscopic surface energy by an amount corresponding to their energy of formation in the configuration relevant to the structure. Furthermore, the energies of individual features and the interactions of these features are influenced by the presence of strain in the crystal. Macroscopically, the surface is assumed to be smooth but to have a local surface energy at a point on the surface determined by the orientation of the tangent plane at that point and the level of elastic strain in the crystal. A quantitative interpretation of the free energy of a strained crystal in terms of the physics of crystallographic steps is developed in the sections that follow. The driving force for alterations in surface morphology is the surface chemical potential, defined in terms of variations of system free energy with surface shape. The surface chemical potential of a stepped surface is also derived, and implications for morphology evolution by surface diffusion are examined. The actual values of the parameters involved in the characterization can be determined only through atomistic simulation or experiment; such work focused on the SiGe material systems is being carried out in parallel with this study.

The present work is motivated by recent experimental observations reported by Sutter and Lagally (2000) and Tromp et al. (2000) on the growth of Si1-xGex films on Si(001) substrates with Ge concentration $`x`$ in the range $`0.1x0.4`$. The experiments have clearly revealed that, during early stages of growth, shallow stepped mounds whose side walls are made up of widely spaced steps emerge as an inherent morphological instability of the film. As more material is deposited, the spacing between the steps gradually decreases until the sidewalls reach a certain crystallographic orientation. This type of nucleationless growth of epitaxial islands is not predicted by a widely used nucleation barrier model for the energetics of faceted islands (Tersoff and LeGoues 1994), which is based on the competition between gain in elastic energy and the energy spent in creating the surface of the island. Approaches to surface instability that account for neither orientation dependence nor strain dependance of surface energy have been unable to explain the observed growth mode. A similar growth mode has also been observed in very recent experiments on Ge films grown on Si(001) substrates (Vailionis et al. 2000, Rastelli et al. 2001), where small islands with widely spaced steps called ”pre-pyramids” are the three dimensional features that first appear during growth; these pre-pyramids eventually evolve to faceted islands. In this paper, we will show that the emergence of stepped islands, without a nucleation barrier, can be explained by incorporating the physics of steps, in particular their interactions and the dependance of their formation energies on the mismatch strain.

Before we proceed with development of the continuum description, we recall a few results regarding stepped crystal surfaces that are well-established in surface science. Below the characteristic roughening temperature, which is in excess of 1000 C for most semiconductors, a nominally flat surface of a crystal that is misoriented by a small angle from a high-symmetry direction consists of a train of straight parallel steps, as illustrated schematically in Fig. LABEL:vicinalsurf. The surface energy of such a “vicinal” surface is given by

$$\gamma (\theta )=\gamma _0\mathrm{cos}\theta +\beta _1|\mathrm{sin}\theta |+\beta _3\frac{|\mathrm{sin}\theta |^3}{\mathrm{cos}^2\theta },$$
(1)
where $`\theta `$ is the misorientation angle, $`\gamma _0`$ the energy density of the atomically flat surface and the parameters $`\beta _1`$ and $`\beta _3`$ are related to step creation and interaction energies respectively. Specifically, $`\beta _1=\beta /h_s`$, where $`\beta `$ is the energy to create a unit length of an isolated step and $`h_s`$ is the height of an atomic step. Like any other defect, say a dislocation or a vacancy, steps give rise to long-range stress fields in the solid; the steps interact with each other, as well as with other defects in the material, through these fields. The amplitude of the stress field produced by a step decays with the inverse square of the distance from the step (Marchenko and Parshin 1980). The last term on the right in (1) arises from the interactions between steps through their stress fields.

The discussion proceeds in the following way. In Sec. 2 we will derive the energy of a stepped surface for small strains but arbitrary orientations, starting from the surface energy of a vicinal surface given in (1). To facilitate calculations, this energy will then be simplified so that it applies to surfaces whose deviation from a flat orientation is small. The surface chemical potential which represents the driving force for change in shape due to mass transport on the surface will be derived in Sec. 3. In Sec. 4, we develop a variational formulation of surface evolution and apply it to study the growth of a thin film on a lattice mismatched substrate under the influence of a constant deposition flux. Our focus will be on the evolution and interactions of stepped mounds that appear on the substrate surface as soon the deposition flux is turned on. It will be shown that the spacing between the steps on the sidewalls of these mounds decreases continuously until it reaches an optimum value, determined by the competition between repulsive step interactions and a strain induced lowering of step formation energies. We provide a summary of the key results along with future directions of research in Sec. 5. The connections between the continuum description developed in this paper and discrete elastic models of steps will be discussed in Appendix A. In Appendix B, we calculate the energy of a strained island using both the continuum and discrete descriptions of stepped surfaces. Here, we show that the mismatch strain can lead to lowering of the surface energy of the vicinal surfaces and can lead to elimination of the barrier to nucleation of these islands.

## 2 Energy of a strained solid with a stepped surface

Consider a strained crystal that occupies the region $`R`$ and that undergoes two-dimensional generalized plane strain deformation. The region is bounded by a free surface $`S`$ as shown in Fig. 2; the remainder of the boundary of $`R`$ is “workless” due to symmetry constraints or its remoteness. The free surface is viewed as a being made up of an array of infinitely long, straight monatomic steps of height $`h_s`$. Initially, the surface $`S`$ is atomically flat with its normal vector in the $`x_2`$direction. The goal of this section is to represent the change in energy associated with departure from flatness in terms of the evolving surface shape.

When the surface $`S`$ is flat, the crystal is subjected to a spatially uniform strain $`ϵ_{ij}^0`$. To ensure that $`S`$ is stress-free in this configuration, this initial strain must satisfy

$$C_{i2kl}ϵ_{kl}^0=0,$$
(2)
where $`C_{ijkl}`$ is the array of elastic constants. The total strain in the crystal is then written as $`ϵ_{ij}=\frac{1}{2}(u_{i,j}+u_{j,i})+ϵ_{ij}^0`$, so that $`u_i`$ represents the additional displacement field that arises as the surface deviates from its flat orientation. The material stretch along the evolving surface can be expressed in terms of the limiting value of the bulk strain on the surface as $`ϵ=m_iϵ_{ij}m_j`$, where $`m_i`$ is a unit vector that is locally tangent to $`S`$. For the surface of the stressed solid under consideration, as illustrated schematically in Fig. 2, we can generalize (1) to include the effects of surface stretch $`ϵ`$ to express the surface energy density as

$$\gamma [h(s),ϵ(s)]=\left(\gamma _0+\tau _0ϵ(s)\right)\sqrt{1h^{}(s)^2}+\left(\beta _1+\stackrel{~}{\beta _1}ϵ(s)\right)\left|h^{}(s)\right|+\beta _3\frac{\left|h^{}(s)\right|^3}{1h^{}(s)^2},$$
(3)
where $`h(s)`$ is the deviation of the height of the surface from its flat orientation as a function of arclength $`s`$ along $`S`$, $`\tau _0=\gamma _0+d\gamma _0/dϵ`$ is the surface stress of the flat surface at the current level of stretch, and $`\stackrel{~}{\beta _1}=d\beta _1/dϵ`$ is a measure of the sensitivity of the formation energy of a step on a strained surface. The last term has the same interpretation as in (1) and its inclusion is essential. In principle, we should also expand the coefficient of $`\left|h^{}(s)\right|^3`$ in terms of surface stretch. However, it will be shown in the appendix, on the basis of direct comparison with discrete step models, that ignoring the higher order terms in the expansion provides a reasonable approximation. Consequently, only the leading order term represented by $`\beta _3`$ is retained.

The total energy of the solid is the sum of the elastic energy and the surface energy,

$$E[h(s),ϵ_{ij}(𝐱)]=\frac{1}{2}_RC_{ijkl}ϵ_{ij}(𝐱)ϵ_{kl}(𝐱)𝑑𝐱+_S\mathrm{\Sigma }(s)𝑑s+_S\tau (s)ϵ(s)𝑑s,$$
(4)
where $`𝐱`$ is any material point in $`R`$, and $`\mathrm{\Sigma }(s)`$ and $`\tau (s)`$ are the local surface energy and surface stress at zero surface stretch. The latter quantities are given by

$$\mathrm{\Sigma }(s)=\gamma _0\sqrt{1h^{}(s)^2}+\beta _1\left|h^{}(s)\right|+\beta _3\frac{\left|h^{}(s)\right|^3}{1h^{}(s)^2}$$
(5)
and

$$\tau (s)=\tau _0\sqrt{1h^{}(s)^2}+\stackrel{~}{\beta _1}\left|h^{}(s)\right|,$$
(6)
respectively. Because mechanical equilibrium is achieved very quickly on the time scale associated with diffusive mass transport on crystal surfaces, we will first derive conditions for equilibrium at fixed surface shape $`h(s)`$. Then, by using the equilibrium fields obtained, we can proceed to derive an expression for the change in the energy of the system in terms of surface shape change, which leads naturally to a corresponding definition of surface chemical potential as an energetic driving force for surface shape change.

Mechanical equilibrium requires that the variation in the free energy $`E`$ at fixed surface shape $`h(s)`$ due to a small perturbation in displacement field $`\delta u_i`$ from equilibrium values of $`u_i`$ must vanish to linear order in $`\delta u_i`$. The first order variation in energy in the present case is

$$\delta E=_RC_{ijkl}u_{l,kj}\delta u_i𝑑𝐱+_SC_{ijkl}\left(u_{l,k}+ϵ_{kl}^0\right)n_j\delta u_i𝑑s+_S\tau (s)m_i\delta u_{i,j}m_j𝑑s,$$
(7)
where $`n_i`$ is the local outward unit normal vector to $`S`$. The last term can be conveniently rewritten as

$$_S\tau (s)m_i\delta u_{i,j}m_j𝑑s=_S\tau (s)m_i\frac{d\delta u_i}{ds}𝑑s.$$
(8)
Because the bounding surface of $`R`$ is workless, we can rewrite this expression in the form

$$_S\tau (s)m_i\frac{d\delta u_i}{ds}𝑑s=_S\frac{d(\tau (s)m_i)}{ds}\delta u_i𝑑s.$$
(9)
Noting that

$$\frac{d(\tau (s)m_i)}{ds}=\frac{d\tau (s)}{ds}m_i+\tau (s)\kappa (s)n_i,$$
(10)
where $`\kappa (s)`$ is the curvature of the surface, which is taken as positive if the center of curvature is outside the material, we can rewrite the variation in the total energy as

$$\delta E=_RC_{ijkl}u_{l,kj}\delta u_i𝑑𝐱+_S\left[C_{ijkl}\left(u_{l,k}+ϵ_{kl}^0\right)n_j\frac{d\tau (s)}{ds}m_i\tau (s)\kappa (s)n_i\right]\delta u_i𝑑s.$$
(11)
Because the variation must vanish for arbitrary $`\delta u_i`$, we obtain the usual equilibrium condition that the divergence of stress must vanish at each material point, $`C_{ijkl}u_{l,kj}=0`$ in the region $`R`$. In addition, traction boundary conditions reflecting the structure of the surface emerge in the form

$`C_{ijkl}u_{l,k}n_jm_i`$ $`=`$ $`C_{ijkl}ϵ_{kl}^0n_jm_i+\stackrel{~}{\beta _1}h^{\prime \prime }(s)\mathrm{Sgn}[h^{}(s)]{\displaystyle \frac{\tau _0h^{}(s)h^{\prime \prime }(s)}{\sqrt{1h^{}(s)^2}}},`$
$`C_{ijkl}u_{l,k}n_jn_i`$ $`=`$ $`C_{ijkl}ϵ_{kl}^0n_jn_i+\left(\tau _0\sqrt{1h^{}(s)^2}+\stackrel{~}{\beta _1}\left|h^{}(s)\right|\right)\kappa (s),`$ (12)
where $`\mathrm{Sgn}[p]=d|p|/dp`$ for any quantity $`p0`$. These must be satisfied on $`S`$ by the elastic field. The first terms on the right sides of (12) represent, respectively, the normal and shear tractions applied to the surface to compensate for the relaxation of the bulk stress $`C_{ijkl}ϵ_{ij}^0`$ due to surface reorientation; these tractions ensure that the surface $`S`$ remains traction free. The other terms on the right side of (12)1 represent the shear traction due to variation of surface stress with position along the surface, whereas the second term in (12)2 represents the normal traction on $`S`$ due to the so-called Laplace pressure which arises from surface stress. Note that these conditions arise only as a result of the assumption that energies of surface features may depend on strain. They involve no other ad hoc assumptions.

Up to this point, there has been no restriction on the magnitude of the surface shape change. From this point onward, however, attention will be restricted to surface shapes that are locally misoriented by only a small angle from the flat orientation, that is, the condition $`|h^{}(s)|1`$ prevails everywhere on $`S`$. Using the fact that $`m_1=1h^{}(x_1)^2/2=n_2`$ and $`m_2=h^{}(x_1)=n_1`$ to second-order in the slope, the boundary conditions become

$`C_{12kl}u_{l,k}`$ $`=`$ $`C_{11kl}ϵ_{kl}^0h^{}(x_1)+\stackrel{~}{\beta _1}h^{\prime \prime }(x_1)\mathrm{Sgn}\left[h^{}(x_1)\right]f_1(x_1),`$
$`C_{22kl}u_{l,k}`$ $`=`$ $`\tau _0h^{\prime \prime }(x_1)f_2(x_1),`$ (13)
where $`f_1`$ and $`f_2`$ are the components of surface traction in the coordinate directions that give rise to stress fields in the crystal associated with shape change. The displacement field at any point in the crystal can now be obtained using the elastic half-space Green’s function $`G_{ij}(𝐱)`$ as

$$u_i(𝐱)=_{\mathrm{}}^{\mathrm{}}G_{ij}(𝐱x^{}𝐞_1)f_j(x^{})𝑑x^{}.$$
(14)
We can now recast the minimum free energy at fixed surface shape in terms of the surface tractions introduced in (13). Applying the divergence theorem to the first term in (4), recalling the workless nature of the boundary constraints, and using the fact that

$`\left(\tau _0+\stackrel{~}{\beta }_1\left|h^{}\right|\right)m_iϵ_{ij}^0m_j`$ $`=`$ $`\left(\tau _0+\stackrel{~}{\beta }_1\left|h^{}\right|\right)ϵ_{11}^0+2\tau _0ϵ_{12}^0h^{}`$ (15)
$`+\left(\tau _0(ϵ_{22}^0ϵ_{11}^0)+2\stackrel{~}{\beta }_1ϵ_{12}^0\mathrm{Sgn}\left[h^{}\right]\right)\left(h^{}\right)^2+O(h^3)`$
we find that

$`E[h(x)]`$ $`=`$ $`E_0+{\displaystyle _{S_1}}\left(\beta _1+\stackrel{~}{\beta _1}ϵ_{11}^0+2\tau _0ϵ_{12}^0\mathrm{Sgn}\left[h^{}(x)\right]\right)\left|h^{}(x)\right|𝑑x+{\displaystyle _{S_1}}\beta _3\left|h^{}(x)\right|^3𝑑x`$ (16)
$`+`$ $`{\displaystyle _{S_1}}\left(\tau _0(ϵ_{22}^0ϵ_{11}^0)+2\stackrel{~}{\beta _1}ϵ_{12}^0\mathrm{Sgn}\left[h^{}(x)\right]\right)\left(h^{}(x)\right)^2𝑑x{\displaystyle \frac{1}{2}}{\displaystyle _{S_1}}f_i(x)u_i(x)𝑑x,`$
where $`E_0`$ is the total energy when the surface is flat, $`S_1`$ refers to the projection of $`S`$ in the $`x_1`$-direction and the coordinate $`x_1`$ has been replaced by $`x`$ for simplicity. This is the total free energy of the system with a perturbed surface shape.

Before proceeding to study the evolution of stepped surfaces, we make a few observations about the continuum representation of system free energy change given by (4) and (16):

1. All the terms in (16) except the last are local, in the sense that the energy density at any point depends only on the local surface slope. The last term, however, represents the interaction of the surface force distributions given in (13); in the jargon of surface science, these represent force monopoles. Since these monopole interactions are long-ranged, the energy density at a point on the surface depends on the surface shape at all points on the surface.
2. The free energy in the small slope approximation is correct to $`O(h^2)`$; only one term of $`O(h^3)`$ has been retained. It will be shown in the Appendix A that this approximation where all the other $`O(h^3)`$ terms are ignored is a reasonable one, by means of comparison of the continuum description with discrete step models.
3. The small slope approximation is not essential to the development. If this approximation is not made, the total energy can be evaluated by solving the boundary value problem, which consists of the bulk equilibrium equations subject to boundary conditions (12), using a numerical technique such as the finite element or boundary element method.
4. The quantity $`\stackrel{~}{\beta _1}`$ must be either positive or negative, if the step formation energy is strain sensitive. It follows that $`\beta _1+\stackrel{~}{\beta _1}ϵ_{11}^0<0`$ for compressive or tensile strain $`ϵ_{11}^0`$, respectively. This implies that the energy of the stepped surface becomes lower than the energy of the surface without steps for strain of one sign or the other. A plot of the surface energy of a vicinal surface for which $`ϵ_{12}^0=0`$ is shown in Fig. 3. In this case, when strains are sufficiently compressive with $`\stackrel{~}{\beta _1}>0`$, the variation of surface energy with surface orientation $`\theta `$, as given by

$$\gamma (\theta )=\gamma _0+\tau _0ϵ_{11}^0+\left(\beta _1+\stackrel{~}{\beta _1}ϵ_{11}^0\right)|\theta |+\tau _0(ϵ_{22}^0ϵ_{11}^0)\theta ^2+\beta _3|\theta |^3,$$
(17)
develops minima away from $`\theta =0`$. At small misorientation angles, the step formation energy dominates the repulsive interaction energy between the steps, while at larger angles, the step interactions are larger in magnitude. The competition between these two opposing effects results in optimum misorientation angles denoted by $`\theta =\pm \theta ^{}`$ in Fig. 3. It follows that $`\stackrel{~}{\beta _1}`$ is a key quantity in determining the morphology of strained crystals111The value of $`\stackrel{~}{\beta _1}`$ cannot be deduced from a continuum analysis. However, by means of atomistic simulations of stepped surfaces on strained crystals of Si being conducted in parallel with this study, we have concluded that the anticipated minima do indeed arise within the strain range of interest.. Interestingly, since the force monopole in (13) has a contribution that depends on $`\stackrel{~}{\beta _1}`$, this key quantity also enters the non-local part of the total energy.
5. If the energy of a stepped surface under strain becomes lower than the energy of a flat surface, there is no barrier to nucleation of epitaxial islands whose sidewalls are made up of surfaces with widely spaced steps. In Appendix B, we consider the energetics of such stepped mounds using both the discrete models of steps and the present continuum description.

## 3 Surface chemical potential

The evolution of surface morphology of strained crystals as a result of mass transport by surface diffusion or other mechanism has been studied extensively on the basis of continuum models, beginning with the introduction of a surface chemical potential by Herring (1953). The essential ingredients for describing evolution are the chemical potential field defined as a function of position over the evolving surface of the solid and a kinetic relationship between surface mass flux and the gradient in chemical potential. The rate of change of surface shape then follows from conservation of mass. The kinetic relationship must be defined in such a way as to ensure that the free energy decreases as surface evolution progresses. In this section, the surface chemical potential for a stepped surface will be obtained from the free energy expression given in (4). The result is written first for arbitrary magnitude of surface slope, and it is then reduced to a simpler form applicable for the case of surface perturbations with slopes of small magnitude.

In deriving the mechanical equilibrium equations in Section 2, it was assumed that the surface was fixed with respect to the material. This was justified by the observation that mechanical fields equilibrate very fast on the timescale of surface evolution. We now assume that the surface shape depends on time $`t`$, and we consider the associated rate of change the free energy of the system as the surface shape changes, with the system being continuously in mechanical equilibrium. The rate of change of this free energy measure then has the form

$$\dot{E}(t)=_S\mu [h(s,t)]v_n(s,t)𝑑s,$$
(18)
where $`v_n(s,t)`$ is the outward normal velocity of the surface $`S`$ with respect to the material instantaneously on it and $`\mu [h(s,t)]`$ is the local chemical potential. At any point along the surface, the quantity $`v_nds`$ is the rate of addition of material volume to the surface. It follows that its coefficient is the change in free energy per unit volume of material added, that is, the coefficient $`\mu `$ represents the surface chemical potential field (Freund 1998).

We can write the rate of change of free energy (4), following application of the divergence theorem to the first term, as

$$\dot{E}=_S\left(C_{ijkl}ϵ_{kl}\dot{u}_in_j+Uv_n\right)𝑑s+_S\left(\frac{\mathrm{\Sigma }}{\theta }\frac{d\theta }{dt}+\frac{\tau }{\theta }\frac{d\theta }{dt}ϵ+\tau \frac{dϵ}{dt}\right)𝑑s_S\left(\mathrm{\Sigma }+\tau ϵ\right)\kappa v_n𝑑s,$$
(19)
where $`U=\frac{1}{2}C_{ijkl}ϵ_{ij}ϵ_{kl}`$ is the elastic strain energy density, a superposed dot denotes a material time derivative, and $`d/dt`$ denotes the total time derivative whenever it must be distinguished from the material derivative. For example, the rate of change of surface stretch

$$\frac{dϵ}{dt}=\dot{ϵ}+m_im_jϵ_{ij,k}n_kv_n$$
(20)
includes a convective term, and the rate of change of surface orientation following the non-material surface requires that

$$\frac{d\theta }{dt}=\frac{v_n}{s}\text{and}\frac{dm_i}{dt}=n_i\frac{v_n}{s}.$$
(21)
The first term in (19) accounts for the change in strain energy in the crystal due to surface evolution, whereas the remaining terms account for the change in surface energy. Recalling that the system is enclosed within workless boundaries, we can rewrite (19) in the form (18) with

$$\mu [h(s,t)]=U\kappa \left[\mathrm{\Sigma }+ϵ\tau +\frac{^2\mathrm{\Sigma }}{\theta ^2}+ϵ\frac{^2\tau }{\theta ^2}\right]\frac{\tau }{\theta }\frac{ϵ}{s}2\frac{}{s}(m_in_jϵ_{ij}\tau )+m_im_jϵ_{ij,k}n_k\tau .$$
(22)
In the limit of small slope $`|h,_x(x,t)|1`$, where the subscript $`x`$ following the comma indicates partial differentiation with respect to $`x`$, this expression can be further simplified to the form

$`\mu [h(x,t)]`$ $`=`$ $`U_0h,_{xx}[(\beta _1+\stackrel{~}{\beta }_1ϵ_{11}^0)\delta \left(h,_x\right)+2\tau _0(ϵ_{22}^0ϵ_{11}^0)+4\stackrel{~}{\beta }_1ϵ_{12}^0\mathrm{Sgn}\left[h,_x\right]+6\beta _3|h,_x\left|\right]`$ (23)
$`+C_{11kl}ϵ_{kl}^0u_{1,x}\stackrel{~}{\beta }_1\mathrm{Sgn}\left[h,_x\right]u_{1,xx}\tau _0u_{2,xx},`$
where $`U_0`$ is the value of $`U`$ when the surface is flat and $`\delta ()`$ is the Dirac delta function.

A vexing feature of the chemical potential function obtained in (23) is the appearance of the Dirac delta function in the term representing step creation energy. We expect this term to play a role in evolution of the surface of the strained crystal because the step creation energy can become negative for strain of one sign or the other, if the magnitude is sufficiently large. While it is possible to introduce any number of ad hoc “regularizations” to deal with this singularity, it is not evident that the results obtained will be independent of the particular choice made. Therefore, we adopt a variational framework for describing surface evolution which is based on the chemical potential derived but which circumvents the problem posed by the singular step creation energy term.

## 4 A variational approach for strain-driven surface evolution

In this section, we apply the continuum description of the energy derived in (16) to study the growth of a strained film bonded to a lattice-mismatched substrate. The constraint of the substrate maintains the mismatch stress in the film. As noted before, recent experiments (Sutter and Lagally 2000 and Tromp et al. 2000) have observed that growth of the deposited film proceeds via formation of shallow mounds whose stepped side-walls ultimately evolve to a faceted orientation. The evolution takes place through diffusion of atoms across the terraces that separate the steps on the side-walls of these mounds.

In the case of terrace diffusion-limited kinetics, the mass flux on the surface can be related to the gradient of the chemical potential derived in (23) through

$$j(x,t)=c\mu _{,x},$$
(24)
where $`c`$ is a coefficient representing surface mobility of diffusing species. If there is a constant deposition flux $`f`$, mass conservation can be invoked to connect the mass flux to the surface shape through the relation

$$h_{,t}+j_{,x}=f.$$
(25)
Our goal is to develop a variational framework for modeling the evolution of these stepped surfaces. If we focus attention on surface profiles that are periodic with wavelength $`\lambda `$, the thickness of the film $`h(x,t)`$, which is illustrated in Fig. 4, can be expressed as

$$h(x,t)=\underset{n=0}{\overset{N_f}{}}a_n(t)\mathrm{cos}(nkx),\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0}x\lambda ,$$
(26)
where $`k=2\pi /\lambda `$. The wavelength should be large enough to accommodate several stepped mounds and other features of interest. The number of terms $`N_f`$ included in the series is finite but otherwise unrestricted. The variational formulation provides first order coupled ordinary differential equations for the Fourier coefficients $`a_n(t)`$; these equations will be numerically integrated to obtain the evolving surface shape for $`t>0`$ using information on the surface shape at $`t=0`$.

Using (26) in (25), we can express the mass flux in the form

$$j(x,t)=\underset{n=1}{\overset{N_f}{}}\dot{a}_n(t)\frac{\mathrm{sin}(nkx)}{nk},$$
(27)
and the mean height of the surface as $`a_0(t)=ft`$. While (24) and (25) constitute a closed set of equations, they cannot be directly integrated to obtain the surface shape because of the singular nature of the chemical potential in (23), as noted in the preceding section. We now turn to the variational approach which provides well-behaved evolution equations for the Fourier coefficients in (26) by using the weak form of (24).

The variational approach can be formulated (Suo 1997) by noting that the functional

$$\mathrm{\Phi }[j(x,t)]=_0^\lambda \mu ,_xj(x,t)dx+_0^\lambda \frac{j^2(x,t)}{2c}dx$$
(28)
attains a minimum value when $`j(x,t)`$ satisfies the kinetic relation (24). The idea is to express the functional in terms of the time derivatives of the Fourier coefficients $`\dot{a}_n(t)`$ and to then minimize the functional with respect to variations in these rates. The first of these tasks can be accomplished by integrating the first term by parts and by using (25) and (27) to obtain

$$\mathrm{\Phi }[\dot{a}_1,\dot{a}_2\mathrm{}]=\underset{n=1}{\overset{N_f}{}}\frac{E[a_0,a_1\mathrm{}]}{a_n}\dot{a}_n+\underset{n=1}{\overset{N_f}{}}\frac{\pi }{2cnk^3}\dot{a}_n^2,$$
(29)
where the total free energy $`E[a_0,a_1\mathrm{}]`$ of the film-substrate system can be directly evaluated from (16) with one minor modification. Since the mismatch strain is present only in the film but not in the substrate (refer to Fig. 4), $`ϵ_{ij}^0`$ in (16) has to replaced by a spatially dependent function $`ϵ_{ij}^0(x_2)=ϵ_{ij}^0H(x_2)`$, where $`H(x_2)`$ is the Heaviside function which has the value of zero when $`x_2<0`$ and the value one when $`x_2>0`$. Minimizing the functional in (29) with respect to the $`\dot{a}_n`$’s, we get the evolution equations

$$\dot{a}_n(t)=\frac{cnk^3}{\pi }\frac{E}{a_n},$$
(30)
which are unambiguous and well-behaved. We now proceed to integrate these rate equations to determine the time evolution of the deposited film. At each time step of integration, the quantity $`E/a_n`$ is calculated numerically. The important role of the nonlinear contribution to (30) for the convergence of this approach is illustrated in Appendix C by considering only a single Fourier mode.

The numerical integration of (30) was carried out using the fourth-order Runge-Kutta procedure with adaptive step-size control as described in detail by Press et al. (1992). The parameters in the continuum description of the energy were chosen as follows: $`\beta _1=0.03J/m^2`$, $`\stackrel{~}{\beta _1}=15J/m^2`$, $`\tau _0=1J/m^2`$, $`ϵ_{11}^0=0.01`$ and $`\beta _3=2.86J/m^2`$. It is assumed that the substrate and the film are isotropic with similar elastic properties, and that the Young’s modulus and the Poisson’s ratio are $`10^{11}N/m^2`$ and $`0.3`$, respectively. Using these parameters, we find that $`ϵ_{22}^0=0.0023`$ and $`ϵ_{12}^0=0`$. Also, for the compressively strained film, the surface energy of the film (sketched schematically in Fig. 3) attains a minimum when $`\theta ^{}=0.12`$, which implies that the sidewalls of the stepped mounds would eventually evolve to this angle.

The initial profile of the film was chosen to be a sinusoid of wavelength $`\lambda =400`$nm, with an amplitude of $`0.4`$ nm, so that the only non-vanishing Fourier components at $`t=0`$ are $`a_0`$ and $`a_1`$. In our calculations, sixteen Fourier coefficients were used to keep track of the surface shape. The evolution of the deposited film is shown in Fig. 5. As the deposition flux is turned on, the material on the surface very quickly gathers into five stepped mounds with slopes much smaller than the optimum value of $`0.12`$ as shown in Fig. 5(a). The insert in Fig. 5(a-c) tracks the evolution of one of these mounds. As more material is deposited, the side-walls of the mound become steeper until they reach the optimum angle $`\theta =0.12`$ in Fig. 5(c). It can also be observed that the center of mass of the mound shifts gradually to the right in going from Fig. 5(a) to Fig. 5(c). This can be understood by looking at the interactions between the islands. Since elastic relaxation is achieved for widely spaced islands, the islands tends to repel each other (Floro et al. 2000). It can be seen that the island in Fig. 5(a) is located closer to the island on its left and would therefore tend to shift towards the right via diffusion of atoms from the side-wall on the left to the one on the right. Once the side-walls of the islands reach the optimum orientation, they grow in a self-similar fashion, until they come in contact with their neighbors as shown in Fig. 5(c). At this point, self-similar coarsening is initiated, which leads to a decrease in areal coverage (fraction of the surface of the substrate covered by the film) as is evident in Fig. 5(e). Here there are two islands, with side walls at $`\theta =0.12`$, separated by 50nm. As more material is deposited, these further grow in size and, eventually, one of the islands grows at the expense of the other island as shown in Fig. 5(f)-(h).

The evolution of the film in the early stages of island growth are in close agreement with the recent experiments of Sutter and Lagally (2000) and Tromp et al. (2000). The key result of these experiments is that the islands evolve as a natural instability without any nucleation barrier. As we show in the Appendix, there is no barrier for nucleating stepped islands if the strains are compressive. Islands with orientations below the optimum facet orientation lower the energy of the film and provide a kinetic pathway for obtaining faceted islands that is free of any nucleation barrier. This is indeed what is seen in experiments and during early stages of growth in Fig. 5(a-c).

It can be seen from Appendix B that the energy of the islands is dominated by the surface energy at small island volumes and by the elastic energy at large volumes. Using the parameters adopted in our calculations, we can use (43) to show that the crossover between these two regimes takes place when the base width of the island is about 200-300nm. When the island sizes are smaller than this value, the sidewalls are oriented at the optimum angle that minimizes the surface energy. With increasing base width, it is known that the islands undergo shape transitions (Medeiros-Ribeiro et al. 1998), where the sidewall angles change to a steeper orientation, which is usually a low-energy crystallographic orientation that makes a larger angle with the substrate. Since such orientations are not yet included in our model, we do not observe such transitions. We do, however, find that the slope of the sidewalls of the large islands in Fig. 5(e) is about $`0.15`$, which is about 25% bigger than the optimum slope. This indicates that elastic energy of these larger islands is becoming comparable to the surface energy.

There is no fundamental impediment to including additional low energy surface orientations in order to study steeper sidewall facets. For example, these might appear as additional relative minima in the variation of $`\gamma (\theta )`$ with orientation $`\theta `$ at angles greater in magnitude than $`\theta ^{}`$. This will be subject of future investigations.

## 5 Summary

In summary, we have developed a continuum framework to model the energetics and evolution of stepped surfaces of nanostructures. The surface energy of strained surfaces is obtained by a generalization of the energy of a vicinal surface to account for the effects of stretching of the surface caused by the mismatch strain. We find that the dependence of the formation energy of steps on the sign of mismatch strain plays a crucial role in establishing the morphology of evolving nanostructures. While this dependence has been seen in atomic-scale simulations, we have developed a consistent continuum framework that naturally accounts for this crucial effect in modeling surface evolution. Furthermore, we have shown in Appendix A that our continuum description provides a direct way to analyze stepped surfaces, without evaluating complicated sums involving discrete steps.

The continuum framework was applied to study the evolution of an epitaxial film bonded to a lattice-mismatched substrate. Our simulations show that the deposited material initially gathers in shallow stepped mounds whose side-wall angles eventually evolve to an orientation that is determined by the competition between strain-induced lowering of the step formation energy and the repulsive interactions between the steps. This kinetic pathway has no nucleation barrier and occurs as a natural instability, in agreement with the recent experimental observations. The simulations were also used to look at self-similar growth of islands leading to their impingement with neighboring islands and subsequent coarsening.

In this article, discussions were limited to two-dimensional deformation fields and evolution of surfaces in one dimension through the diffusion atoms, all of the same type. Work is currently in progress to extend the current formulation to study the evolution of two-dimensional surfaces. Also, since the films of interest in device applications are typically alloys of two different semiconductors, further insight into surface evolution can be obtained by including the possibility of alloy segregation in binary materials.

## Acknowledgments

The research support of the National Science Foundation through grant CMS-0093714 and the Brown University MRSEC Program, under award DMR-0079964, is gratefully acknowledged.

Appendices

## Appendix A Discrete surface step models

In developing a continuum description of stepped surfaces, we started with the surface energy of a strained vicinal surface given by (3), which is a generalization of the surface energy of the vicinal surface to include the effects of surface stretch. The key idea of this section is to demonstrate that if one chooses to start with elastic models of discrete steps and proceeds to construct a homogenized continuum description, then the elastic displacement fields (given in (14)) and the energy of stepped surfaces (given in (16)) are recovered. This procedure is very instructive, as it will allow us to give physical interpretation for parameters such as $`\stackrel{~}{\beta }_1`$ and $`\beta _3`$ in the continuum description in terms of quantities that determine the elastic fields of discrete atomic steps. As we will see, these quantities depend on bonding of the atoms near the step edges and therefore have to be determined using an atomic level calculation or, alternatively, from experiments. In what follows, we will first look at the elastic displacement fields and the corresponding strain energy stored in the crystal with a traction-free stepped surface. As noted earlier, the contribution of the strain energy to the total free energy enters through the non-local contribution in (16). We will then proceed to analyze the local part of the free energy that depends on the local slope of the surface.

### A.1 Non-local parts of the free energy and elastic displacement fields

In their pioneering work on step interactions, Marchenko and Parshin (1980) analyzed the strain field of a surface step by viewing the step as a point dipole on the surface of strength $`𝐃=D_1𝐞_1+D_2𝐞_2`$, with components in both the coordinate directions perpendicular to the step which lies along the $`𝐞_3`$-direction in the present case. Since the $`𝐞_3`$ component of the moment should vanish in equilibrium, they concluded that $`D_2=\pm \tau _0h_s`$, where $`\tau _o`$ is the surface stress and $`+()`$ denotes a positive(negative) step222For a positive step located at $`x_0`$, the surface height satisfies $`h(x_0^+)h(x_0^{})=h_s`$, where $`h(x_0^+)`$ ($`h(x_0^{})`$) is the surface height slightly to the right (left) of the step, while for a negative step the relation $`h(x_0^+)h(x_0^{})=h_s`$ is satisfied. respectively. This value of $`D_2`$ ensures that the moment created by the surface stress is compensated. The other component, $`D_1`$, depends on details of atomic bonding around the step edge and has to be obtained from an atomic level calculation or from experiments. The 1-component is a dipole without moment, whereas 2-component is a dipole with moment, as noted. If we focus attention on crystal structures that are invariant under rotations of 180 about the $`𝐞_2`$ axis, it is clear that a positive step becomes a negative step under such rotations. While the 2-component flips its sign under such a rotation, the 1-component does not. Therefore, on such high symmetry surfaces (for example, the (100) surfaces of cubic materials commonly employed in strained epitaxy), both positive and negative steps have identical values for $`D_1`$. In addition to these components of the force dipole, for a strained crystal equilibrium requirements can be invoked to show that a monopole $`M=h_sC_{11kl}ϵ_{kl}^0`$ in the $`𝐞_1`$-direction must be included, in addition to the point dipole, to model the stress fields of the step associated with the mismatch strain (Tersoff et al. 1995). We can now use the surface Green’s function to write the elastic displacement field of an isolated step located at the origin as

$$u_i(𝐱)=\pm MG_{i1}(𝐱)+D_1\frac{G_{i1}(𝐱)}{x_1}\pm D_2\frac{G_{i2}(𝐱)}{x_1}$$
(31)
for a positive/negative step.

If we consider the surface shown in Fig. 2, with a continuous distribution of steps, the surface forces can be modeled using continuum monopole and dipole densities

$$m(x)=\frac{M}{h_s}\frac{dh}{dx}\mathrm{and}𝐝(x)=\frac{D_1}{h_s}\left|\frac{dh}{dx}\right|𝐞_1+\frac{D_2}{h_s}\frac{dh}{dx}𝐞_2$$
(32)
respectively. The displacement fields in the crystal due to of distribution of steps can be written as a superposition of the displacement fields of individual steps, or

$`u_i(𝐱)`$ $`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}[{\displaystyle \frac{M}{h_S}}G_{i1}(𝐱x^{}𝐞_1){\displaystyle \frac{dh(x^{})}{dx^{}}}{\displaystyle \frac{D_1}{h_s}}{\displaystyle \frac{G_{i1}(𝐱x^{}𝐞_1)}{x^{}}}|{\displaystyle \frac{dh(x^{})}{dx^{}}}|`$ (33)
$`\tau _0{\displaystyle \frac{G_{i2}(𝐱x^{}𝐞_1)}{x^{}}}{\displaystyle \frac{dh(x^{})}{dx^{}}}]dx^{}.`$
Once the second and third terms are integrated by parts, it is easily seen that (33) can be written as

$`u_i(𝐱)`$ $`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}[{\displaystyle \frac{M}{h_S}}G_{i1}(𝐱x^{}𝐞_1){\displaystyle \frac{dh(x^{})}{dx^{}}}+{\displaystyle \frac{D_1}{h_s}}G_{i1}(𝐱x^{}𝐞_1){\displaystyle \frac{d^2h(x^{})}{dx^2}}\mathrm{Sgn}\left[{\displaystyle \frac{dh(x^{})}{dx^{}}}\right]`$ (34)
$`+\tau _0G_{i2}(𝐱x^{}𝐞_1){\displaystyle \frac{d^2h(x^{})}{dx^2}}]dx^{},`$
which is in the form given in (14) with the surface traction $`𝐟`$ precisely as given by (13) if we identify $`\stackrel{~}{\beta _1}`$ with $`D_1/h_s`$. It can also be verified that the contribution to the free energy from this displacement field is

$$E_{nonlocal}=\frac{1}{2}_{S_1}𝐟𝐮𝑑x_1,$$
(35)
which is nothing but the last term in (16) which represents the non-local contribution to the free energy.

To summarize, an interesting outcome of the simple homogenization process that we have carried out is that even though individual steps are modeled using both monopoles and dipoles, the displacement fields are determined solely in terms of a distribution of monopole-like surface tractions given in (13). The physical reason for this can be understood by appealing to the analogous homogenization problem of dipoles in electrostatics or magnetostatics. In a polarized dielectric with dipole moment density $`𝐏(𝐱)`$, the homogenized electric fields can be obtained by replacing the spatially varying dipole density with an equivalent charge density $`\rho (𝐱)=𝐏(𝐱)`$; this means that a spatial variation of the dipole density physically corresponds to charge accumulation. In the case of steps, the force monopole distribution that originated from the dipoles in (13) is seen to be nothing but $`d𝐝(x)/dx`$, where $`𝐝(x)`$ is the dipole density introduced earlier. It is also clear that a homogenous distribution of dipoles does not give rise to any long-range elastic displacement fields; if we consider a vicinal surface shown in Fig. 1, for regions with size comparable to or larger than the step spacing, symmetry arguments can be invoked to show that the average or homogenized displacements in such regions vanish. Of course, on length scales much smaller than the step spacing there can be local field fluctuations due to individual dipoles. However, at scales larger than the step spacing, these fluctuations average out to zero and the homogenized field is determined solely by the spatial variations in the step density.

### A.2 Local parts of the free energy

The local part of the free energy has contributions from two sources, one from the self-energy of the force dipoles at the steps and the other due to the interaction of these dipoles with the epitaxial mismatch stress. In each of these cases, the interactions arise from the local displacement fields near the steps as opposed to the homogenized displacement fields calculated in (34). Below, we consider each of these contributions separately.

The self-energy of a force dipole distribution can be obtained from the interaction energy between force dipoles on a vicinal surface. If we consider the vicinal surface in Fig. LABEL:vicinalsurf, the interaction energy of one of the dipoles with all the others can be expressed as

$$E_{dipole}=\frac{\pi ^2}{6}\left[D_1^2\frac{^2}{^2x}G_{11}(x𝐞_1)|_{x=1}+D_2^2\frac{^2}{^2x}G_{22}(x𝐞_1)|_{x=1}\right]\frac{1}{d^2}\frac{G}{d^2},$$
(36)
where $`d=h_s/\mathrm{tan}\theta `$ is the spacing between steps and the factor $`\pi ^2/6`$ comes from evaluating the infinite sum $`_{i=1}^{\mathrm{}}i^2`$. The self-energy of a continuous distribution of dipoles can be calculated using (36) as

$$E_{self}=_{S_1}\frac{G}{d^2(x_1)}\frac{|h^{}(x_1)|}{h_s}𝑑x_1=_{S_1}\frac{G}{h_s^3}|h^{}(x_1)|^3𝑑x_1,$$
(37)
where the $`local`$ step spacing $`d(x_1)`$ is expressed in terms of the surface slope as $`d(x_1)=h_s/|h^{}(x_1)|`$. The reason that the local approximation works for the self-energy is that the inverse square decay of interaction between dipoles guarantees that the contribution of steps that are far from a given step are much smaller in magnitude than the contributions from those that are close by. In distinct contract, the logarithmic nature of the interaction between force monopoles (given in (13)) that arise from the spatial variation of the dipoles, requires us to assign equal importance to steps that are both near and far. Hence, this contribution must be treated in the non-local approximation.

The coupling of the strain field of a step with the mismatch strain gives rise to another local contribution to the free energy, which can also be evaluated by considering the vicinal surface in Fig. LABEL:vicinalsurf. The interaction energy between a positive step and the mismatch strain is $`D_1u_{1^{},1^{}}+|D_2|u_{2^{},1^{}}`$, where $`u_{1^{},1^{}}`$ and $`u_{1^{},2^{}}`$ are displacement gradients in the transformed coordinate system shown in Fig. 6. The coordinate axes in this transformed system are chosen so that the $`1^{}`$-axis lies along the vicinal surface and the $`2^{}`$-axis coincides with the normal to the vicinal surface. Using the standard rules for transformation of strains, we have

$`u_{1^{},1^{}}`$ $`=`$ $`ϵ_{11}^0+2ϵ_{12}^0\theta +O(\theta ^2)`$
$`u_{2^{},1^{}}`$ $`=`$ $`u_{2,1}^0+(ϵ_{22}^0ϵ_{11}^0)\theta +O(\theta ^2).`$ (38)
For a continuous surface, the total contribution to the local part of the free energy can be obtained by summing the interaction energies of each of the steps on the surface with the mismatch strain, or

$$E_{dipolestrain}=_S\left[D_1\frac{|h^{}(x_1)|}{h_s}(ϵ_{11}^0+2ϵ_{12}^0h^{}(x_1))+|D_2|\frac{h^{}(x_1)}{h_s}(u_{2,1}^0+(ϵ_{22}^0ϵ_{11}^0)h^{}(x_1))\right]𝑑x_1.$$
(39)
Making use of the fact that $`D_1=\stackrel{~}{\beta _1}h_s`$ and $`|D_2|=\tau _0h_s`$, it can be verified that the terms of $`O(h^{})`$ and $`O(h^2)`$ in (16) add up to the local contribution evaluated in (39) save the factor proportional to $`u_{2,1}^0`$. We do not have an explanation for this difference.

Within the continuum description the parameter $`\stackrel{~}{\beta }_1`$ was introduced to include the effect of surface stretch on the step formation energy. Comparing the continuum description with the discrete step models, we conclude that this parameter is closely related to the intrinsic dipole at the step. If there are significant modifications in atomic bonding at the step edges, the intrinsic dipole will be large in magnitude. For such steps, strain can significantly reduce their formation energy. Atomistic calculations on reconstructed Si steps (Xie et al. 1994, Roland 1995) have shown large reduction in step formation energies due to strain. The results, however, were neither interpreted in terms of the intrinsic dipole nor used in a continuum framework to describing evolution of surfaces. Our continuum formalism accounts this crucial effect in the description of evolution of strained material surfaces.

## Appendix B Energetics of strained islands

### B.1 Continuum description

We will now apply the continuum model to calculate the energy of strained islands. The island energies will also be obtained directly by summing up the interactions of individual steps. It will be shown that the continuum model not only provides an easy way of calculating the energetics but also provides a direct of identifying the contributions to the total energy arising from edges and corners. As we will see, extracting these quantities from the discrete step picture is much more involved and tedious.

Let us consider the island made up of $`N`$ positive and negative steps with base width $`2L`$ as shown in Fig. 7; the angle made by the facet is given by $`\theta Nh_s/L`$ and the two-dimensional volume of the island is $`A=LNh_s`$. For the sake of simplicity, the surface Green’s function of an isotropic material will be employed in evaluating the non-local contribution to the energy of the island.

Within the continuum description, the local contribution to energy of the island is

$$E_L=2A^{1/2}\theta ^{1/2}\left[(\beta _1+\stackrel{~}{\beta }_1ϵ_{11}^0)+\tau _0(ϵ_{22}^0ϵ_{11}^0)\theta +\beta _3\theta ^2\right],$$
(40)
while the non-local contributions to the energy can be calculated from the surface tractions given in (13). Since the slope on each of the side-walls of the island is a constant, the components of the surface tractions are given by

$`f_1`$ $`=`$ $`{\displaystyle \frac{M}{h_s}}\theta \left[\mathrm{\Theta }(xL)+\mathrm{\Theta }(x+L)2\mathrm{\Theta }(x)\right]+\stackrel{~}{\beta }_1\theta \left[\delta (xL)\delta (x+L)\right],`$
$`f_2`$ $`=`$ $`\tau _0\theta [2\delta (x)+\delta (xL)+\delta (x+L)]`$ (41)
where $`\mathrm{\Theta }()`$ is the unit step function and $`\delta ()`$ is the Dirac delta function. The first term in (41)1 comes from the force monopoles at the steps, while the other terms that act at the apex and the edges of the island, as shown in Fig. 7, come from the force dipoles at the step. The different components of the non-local part of the island energy can be written as

$`E_{MM}`$ $`=`$ $`{\displaystyle \frac{\alpha M^2\theta ^2}{h_s^2}}\left[2{\displaystyle _0^L}(Lx)\mathrm{log}xdx{\displaystyle _L^L}(L|x|)\mathrm{log}(L+x)𝑑x\right]=\alpha A\theta (C_{11ij}ϵ_{ij}^0)^2\mathrm{log}4,`$
$`E_{MD}`$ $`=`$ $`{\displaystyle \frac{2\alpha M\stackrel{~}{\beta }_1\theta ^2}{h_s}}{\displaystyle _0^L}\mathrm{log}(2L/x1)𝑑x=2\alpha \stackrel{~}{\beta }_1A^{1/2}\theta ^{3/2}(C_{11ij}ϵ_{ij}^0)\mathrm{log}4,`$
$`E_{DD}`$ $`=`$ $`\alpha \stackrel{~}{\beta }_{1}^{}{}_{}{}^{2}\theta ^2\mathrm{log}(2L/h_s)\alpha \tau _{0}^{}{}_{}{}^{2}\theta ^2[4\mathrm{log}(L/h_s)\mathrm{log}(2L/h_s)]`$ (42)
$``$ $`\left[\alpha \stackrel{~}{\beta }_{1}^{}{}_{}{}^{2}+3\alpha \tau _{0}^{}{}_{}{}^{2}\right]\theta ^2\mathrm{log}(2L/h_s),`$
where the subscripts $`M`$ and $`D`$, respectively, indicate terms arising from the monopoles and dipoles at the steps, and $`\alpha =\frac{2(1\nu ^2)}{\pi E}`$.

The relative importance of the edge contribution to the island energy, given by (42)2 plus (42)3, can be determined by comparing its magnitude with the local contribution given by (40). It can easily be verified that the ratios $`E_{MD}/E_L`$ and $`E_{DD}/E_L`$ are equal to $`\theta `$ and $`\mathrm{log}(\stackrel{~}{A})/\stackrel{~}{A}`$ within factors of order unity, respectively, where $`\stackrel{~}{A}=A/h_s^2`$. Since our focus is on shallow islands ($`\theta 1`$) whose dimensions are large compared to atomic dimensions, it is clear that the edge contribution can be safely ignored and the energy of the island becomes

$$E_{island}=2A^{1/2}\theta ^{1/2}\left[(\beta _1+\stackrel{~}{\beta }_1ϵ_{11}^0)+\tau _0(ϵ_{22}^0ϵ_{11}^0)\theta +\beta _3\theta ^2\right]\alpha A\theta (C_{11ij}ϵ_{ij}^0)^2\mathrm{log}4.$$
(43)
If the strain dependence of the first term is ignored, there is an energy barrier to nucleation of epitaxial islands (Tersoff and Le Goues 1994). However if the condition $`\beta _1+\stackrel{~}{\beta }_1ϵ_{11}^0<0`$ is met, for side-wall angles in the range $`0\theta \theta ^{}`$ (see Fig. 3 ) the nucleation barrier is absent. Indeed, experimental observations (Sutter and Lagally 2000, Tromp et al. 2000) confirm that such shallow stepped mounds islands form as a natural instability during strained heteroepitaxy.

### B.2 Discrete description

The energy of the strained island in Fig. 7 can be also obtained considering the interaction energy of discrete steps that make up the sidewalls of the islands. Using italicized symbols to denote energies, the local part of the energy that includes that step formation energy and the interaction of the force dipole with the mismatch strain can be expressed as

$$_L=2N\left[(\beta _1+\stackrel{~}{\beta }_1ϵ_{11}^0)+\tau _0(ϵ_{22}^0ϵ_{11}^0)\theta \right].$$
(44)
If we compare this expression with the local contribution in the continuum description, we find that the self-energy of the force dipoles does not appear in the discrete framework. The reason is that within the discrete model, both the self-energy and the non-local part due the spatial variation of the step density have to be evaluated from the dipole-dipole interactions between the steps that constitute the island. The total interaction energy, which is a sum of interactions between pairs of steps can be written as

$$_{int}=\frac{1}{2}\underset{i,j}{}\left[\left(M_iM_j(M_iM_j)D_1\frac{}{x_{ij}}+D_1^2\frac{^2}{^2x_{ij}}\right)G_{11}(x_{ij}𝐞_1)+D_2^2\frac{^2}{^2x_{ij}}G_{22}(x_{ij}𝐞_1)\right],$$
(45)
where the steps are labeled with indices $`i`$ and $`j`$ and $`x_{i,j}=x_ix_j`$. Below, we will evaluate the total energy for the island in Fig. 7 by using the elastic Green’s function of the isotropic solid.

The double sum in (45) can be simplified by rearranging the terms in (45) (Shenoy et al. 1998, 2000), so that sum over one of the indices can be conveniently evaluated; the step interaction energies are then given by

$`_{MM}`$ $`=`$ $`\alpha M^2\left[2{\displaystyle \underset{i=1}{\overset{N1}{}}}(Ni)\mathrm{log}i{\displaystyle \underset{N+1}{\overset{N1}{}}}(N|i|)\mathrm{log}(N+i)\right],`$
$`_{MD}`$ $`=`$ $`2\alpha M\stackrel{~}{\beta }_1\theta {\displaystyle \underset{i=N+1}{\overset{N1}{}}}{\displaystyle \frac{N|i|}{N+i}},`$
$`_{DD}`$ $`=`$ $`\alpha \stackrel{~}{\beta }_{1}^{}{}_{}{}^{2}\theta ^2{\displaystyle \underset{i=1}{\overset{2N1}{}}}{\displaystyle \frac{2Ni}{i^2}}+\alpha ^2\tau _0^2\theta ^2\left[2{\displaystyle \underset{i=1}{\overset{N1}{}}}{\displaystyle \frac{Ni}{i^2}}+{\displaystyle \underset{i=N+1}{\overset{N1}{}}}{\displaystyle \frac{N|i|}{(N+i)^2}}\right].`$ (46)
When $`N`$ becomes large, the discrete sums in (46)1 and (46)2 can be evaluated by converting to integrals. The results so obtained are identical to the continuum results in (42). The dipole-dipole term given by (46)3, in the limit of large $`N`$, yields

$$_{DD}=\frac{\alpha \theta ^2\pi ^2}{3}[\stackrel{~}{\beta }_{1}^{}{}_{}{}^{2}+\tau _0^2]\left[\alpha \stackrel{~}{\beta }_{1}^{}{}_{}{}^{2}+3\alpha \tau _{0}^{}{}_{}{}^{2}\right]\theta ^2\mathrm{log}(2L/h_s).$$
(47)
The first term is identical to the self-energy term of the dipoles in the continuum description, while the next two terms agree with the non-local contribution given in (42)3. Thus, the expressions for energy of the island computed using both the discrete and continuum descriptions are in agreement.

The above exercise shows that the continuum model provides a direct route to computing the energy of stepped surfaces without recourse to tedious evaluation of discrete sums. It also provides a direct way of looking at the energies associated with corners and edges. It is difficult to obtain these contributions form the discrete approach because the self-energy and the long-range part of the dipole interactions are not clearly separated.

## Appendix C Convergence of the variational formulation

The variational formulation is based upon the expression of the surface shape in terms of a Fourier series with a finite number of expansion coefficients; refer to (26). In this appendix we investigate the convergence of the expansion, focusing on the growth a single Fourier mode. The goal of this exercise is to show that the amplitudes of the modes with small wavelengths do not grow without bound.

In order to perform the stability analysis of short wavelength modes, the highly nonlinear nature of the evolution equations behooves us to consider not just the leading order contribution of the energy represented by the negative step formation energy, but also the repulsive step interaction energy. It can be shown that the nonlocal contribution, given by the last term in (16), becomes significant only when the wavelengths are large. If only the second and third terms in (16) are retained, the evolution of the surface shape consisting of a single mode,

$$h(x)=a_n(t)\mathrm{cos}(nkx),$$
(48)
can be written as

$$\dot{a}_n(t)=\frac{4cn^2k^3}{\pi }\left[\beta _1+\stackrel{~}{\beta }_1ϵ_{11}^0+2\beta _3a_n^2k^2n^2\right]\mathrm{Sgn}[a_n].$$
(49)
If $`\beta _1+\stackrel{~}{\beta }_1ϵ_{11}^0<0`$, the growth of this mode is determined by a competition between the negative step formation energy, which promotes the growth of small wavelength modes, and the repulsive step interactions. It is evident that $`a_n`$ achieves the stable amplitude

$$|a_n|=\frac{1}{nk}\sqrt{\frac{(\beta _1+\stackrel{~}{\beta }_1ϵ_{11}^0)}{2\beta _3}}.$$
(50)
Because $`|a_n|`$ vanishes as $`n\mathrm{}`$, the use of a finite number of terms in the expansion of surface shape in (26) is reasonable. This expectation is borne out by numerical experiments in the full problem.

In the above analysis, we did not include the effects of the substrate. Since the mismatch strain is present only in the deposited film, the step formation energy in the substrate is always positive. This implies that the amplitude obtained in (50) represents an upper bound of sorts; if the effect of the substrate is included, the stable amplitude would be smaller.

## References

* R. J. Asaro and W. A. Tiller, Interface morphology development during stress corrosion cracking: part I. via diffusion, Metall. Trans. 3, 1789-1796 (1972).
* J. A. Floro, M. B. Sinclair, E. Chason, L. B. Freund, R. D. Twesten, R. Q. Hwang and G. A. Lucadamo, Novel SiGe island coarsening kinetics: Ostwald ripening and elastic interactions, Phys. Rev. Lett. 84, 701-704 (2000).
* L. B. Freund, Evolution of waviness of the surface of a strained elastic solid due to stress-driven diffusion, Int. J. Solids Structures 32, 911-923 (1995).
* L. B. Freund, A surface chemical potential for elastic solids, J. Mech. Phys. Sols. 46, 1835-1844 (1998).
* H. Gao and W. D. Nix, Surface Roughening of heteroepitaxial thin films, Annu. Rev. Mater. Sci. 29, 173-209 (1999).
* M. A. Grinfeld, Instability of the separation boundary between non-hydrostatically stressed elastic solid and melt, Sov. Phys. Dokl. 31, 831-834 (1996).
* C. Herring, The use of classical macroscopic concepts in surface energy problems, Structure and Properties of Solid Surfaces, edited by R. Gomer and C. S. Smith, University of Chicago Press, 5-72 (1953).
* V. I. Marchenko and A. Ya. Parshin, Elastic properties of crystal surfaces, Sov. Phys. JETP 52, 129-131 (1980).
* G. Medeiros-Ribeiro, A. M. Bratkovski, T. I. Kamins, D. A. A. Ohlberg and R. S. Williams, Shape transition of Ge nanocrystals on Si (001) surface from pyramids to domes, Science 279, 353-355 (1998).
* W. H. Press, S. A. Teukolsky and W. T. Vetterling and B. P. Flannery, ”Numerical recipes: The art of scientific computing”, (Cambridge University Press, Cambridge, 1992) pages 708-715.
* A. Rastelli, M. Kummer and H. von Kanel, Reversible shape evolution of Ge islands on Si(001), Phys. Rev. Lett. 87, 6101-6104 (2001).
* C. Roland, Effect of stress on step energies and surface roughness, MRS Bulletin 21, 27-30 (1996).
* V. B. Shenoy, S. Zhang and W. F. Saam, Bunching transitions on vicinal surfaces and quantum-$`n`$ mers, Phys. Rev. Lett. 81, 3475-3478 (1998).
* V. B. Shenoy, S. Zhang and W. F. Saam, Step bunching transitions on vicinal surfaces with attractive step interactions, Surf. Sci. 467 58-84 (2000).
* D. J. Srolovitz, On the stability of surfaces of stressed solids, Acta. Metall. 37, 621-625 (1989).
* Z. Suo, Motion of microscopic surfaces in materials, Adv. Appl. Mech. 33, 193-294 (1997).
* P. Sutter and M. G. Lagally, Nucleationless three-dimensional island formation in low-misfit heteroepitaxy, Phys. Rev. Lett. 84, 4637-4640 (2000).
* J. Tersoff and F. K. Le Goues, Competing relaxation mechanisms in strained layers, Phys. Rev. Lett. 72, 3570-3573 (1994).
* J. Tersoff, Y. H. Phang, Z. Y. Zhang and M. G. Lagally, Step-bunching instability of vicinal surfaces under stress, Phys. Rev. Lett. 75, 2730-2733 (1995).
* R. M. Tromp, F. M. Ross and M. C. Reuter, Instability-driven SiGe island growth, Phys. Rev. Lett. 84, 4641-4644 (2000).
* A. Vailionis, B. Cho, G. Glass, P. Desjardins, D. G. Cahill and J. E. Greene, Pathway to strain-driven two-dimensional to three-dimensional transition during growth of Ge on Si(001), Phys. Rev. Lett. 85, 3672-3675 (2000).
* Y. H. Xie, G. H. Gilmer, C. Roland, P. J. Silverman, S. K. Buratto, J. Y. Cheng, E. A. Fitzgerald, A. R. Kortan, S. Schuppler, M. A. Marcus and P. H. Citrin, Semiconductor surface roughness - dependence on sign and magnitude of bulk strain, Phys. Rev. Lett. 73, 3006-3009 (1994).
* Y. W. Zhang and A. F. Bower, Numerical simulations of island formation in a coherent strained epitaxial thin film system, J. Mech. Phys. Sols. 47, 2273-2297 (1998).


text:
# Search for 𝑹_𝒑-violating SUSY and excited fermions at HERA

## 1 Introduction

At the HERA collider, electrons (positrons) and protons collide at a center of mass energy of about $`\sqrt{s}=318`$ GeV (300 GeV before 1998). Since 1994, integrated luminosities of about $`=110\mathrm{pb}^1`$ in $`e^+p`$ and $`=15\mathrm{pb}^1`$ in $`e^{}p`$ scattering have been collected by each of the two experiments H1 and ZEUS. These data are used to search for new physics beyond the SM. Recent results on searches for excited fermions and for R-parity-violating SUSY are presented in this contribution.

## 2 Excited fermions

Fermionic substructure as predicted in compositeness models naturally leads to excitations of fermions. At HERA, such excited fermions could be singly produced via the t-channel exchange of a gauge boson, and would subsequently decay into a SM fermion and a boson. In the 1994–97/1999–2000 $`e^+p`$ and the 1998–99 $`e^{}p`$ data, both collaborations have found no evidence for a deviation from the SM due to the production of an excited electron111If not particularly emphasized, electron can mean either an electron or a positron., neutrino or quark, which decay into $`\gamma `$, $`Z`$ or $`W`$ and a SM fermion, followed by the subsequent decay of the boson into $`e`$, $`\mu `$, $`\nu `$ or hadrons.

To determine the experimental sensitivity, the data have been interpreted in the framework of a compositeness model where the excited fermion $`f^{}`$ is a composite of a scalar and a spin $`\frac{1}{2}`$ constituent . If $`f^{}`$ carries spin and isospin $`\frac{1}{2}`$, the effective Lagrangian reads

$`_{eff}={\displaystyle \frac{1}{2\mathrm{\Lambda }}}\overline{F}_{R}^{}{}_{}{}^{}\sigma ^{\mu \nu }(gf{\displaystyle \frac{\tau ^a}{2}}W_{\mu \nu }^a+g^{}f^{}{\displaystyle \frac{Y}{2}}B_{\mu \nu }+g_sf_s{\displaystyle \frac{\lambda _a}{2}}G_{\mu \nu }^a)F_L+h.c.,`$ (1)
with the standard electroweak and strong couplings denoted by $`g`$, $`g^{}`$ and $`g_s`$ and the SM gauge field strength tensors $`W_{\mu \nu }^a`$ (SU(2)), $`B_{\mu \nu }`$ (U(1)) and $`G_{\mu \nu }^a`$ ($`\mathrm{SU}(3)_\mathrm{C}`$). $`\tau ^a`$, $`Y`$ and $`\lambda _a`$ are the Pauli matrices, the weak hypercharge operator and the Gell-Mann matrices, respectively. The new couplings between the weak isodoublets $`F_R^{}`$, $`F_L^{}`$ of the excited fermion fields and the SM fermion fields $`F_R`$, $`F_L`$ are modified by $`f`$, $`f^{}`$ and $`f_s`$. The compositeness scale is $`\mathrm{\Lambda }`$. The coupling between e.g. a photon, an excited fermion and a SM fermion is given by $`c_{\gamma f^{}f}=\frac{1}{2}(fI_3+f^{}\frac{Y}{2})`$, where $`I_3`$ is the third component of the isospin of the fermion. Assuming relations between the new couplings, upper limits on $`f`$ (or $`f/\mathrm{\Lambda }`$) at the 95% confidence level (CL) have been obtained.

For excited electrons, the ZEUS upper limits have been improved recently by a factor of 2 – 3 and are competitive with those from LEP for larger masses $`m_e^{}`$  as can be inferred from Figure 1a. Figure 1b shows that, for small values of $`f_s`$, the H1 upper limits on $`f`$ ($`f=f^{}`$; $`\mathrm{\Lambda }=M(q^{})`$) for excited quarks are more restrictive than those of the TeVatron for not too large masses $`M(q^{})`$. The HERA results obtained on $`q^{}`$ production via the electroweak couplings $`f`$ and $`f^{}`$ are complementary to the results obtained at the TeVatron, where $`q^{}`$ production occurs via the strong coupling $`f_s`$.

From Figure 1c it can be seen that the H1  and ZEUS  upper limits on $`f/\mathrm{\Lambda }`$ ($`f=f^{}`$) for excited neutrinos deliver a unique sensitivity beyond the LEP center of mass energy.222A similar picture can be found for the upper limits for excited quarks. In addition, H1 has obtained less model-dependent upper limits on $`f/\mathrm{\Lambda }`$ for arbitrary ratios $`f^{}/f`$ . This is presented in Figure 1d. Since $`\nu ^{}`$ results were obtained from a relatively small amount of $`e^{}p`$ 333In this case, $`e^{}p`$ scattering is superior to $`e^+p`$ because the dominant valence quark interaction is helicity suppressed and has smaller parton densities in $`e^+p`$ ($`d`$) compared to $`e^{}p`$ ($`u`$). data, new $`e^{}p`$ luminosity from 2001 (see Section 4) will give substantially higher sensitivity to excited neutrinos and are expected to improve on the LEP limits, even for masses below 200 GeV.

## 3 Squarks in R-Parity-violating SUSY

In the Minimal Supersymmetric Standard Model (MSSM), there is a multiplicative quantum number called R-Parity which is defined by $`R_p(1)^{3B+L+2S}`$, with the baryon number is denoted $`B`$, the lepton number $`L`$ and the spin quantum number $`S`$. For SM particles (sparticles, the supersymmetric partners of the particles), this quantum number acquires the value $`R_p=1`$ ($`R_p=1`$). The $`R_p`$-conserving MSSM has been investigated at HERA in the associated $`\stackrel{~}{e}`$$`\stackrel{~}{q}`$ (selectron-squark) production .

However, the most general supersymmetric and gauge-invariant superpotential contains three R-Parity-violating Yukawa coupling terms, two operators mediating $`L`$ violation and one $`B`$ violation . Of particular interest for HERA are the terms $`\lambda _{ijk}^{}L_iQ_j\overline{D}_k`$ with $`L_i`$ ($`Q_j`$) the lepton (quark) $`SU(2)_L`$ doublet superfields and $`\overline{D}_k`$ the down-quark $`SU(2)_L`$ singlet superfields . The R-Parity-violating dimensionless coupling constants are called $`\lambda _{ijk}^{}`$ , where $`i`$, $`j`$, $`k`$ denote the generation. For non-vanishing $`\lambda _{1jk}^{}`$ , the resonant production of single squarks is possible in $`ep`$ scattering. At HERA with an $`e^+`$ beam, predominantly the $`\lambda _{1j1}^{}`$ couplings are probed.

Both H1 and ZEUS have searched for squark production using $`e^+p`$ data from 1994–97 at a center-of-mass energy $`\sqrt{s}=300`$ GeV distinguishing between R-Parity-violating $`\mathrm{}q`$ ($`\mathrm{}=e,\nu `$) decays leading to lepton + jet topologies and decays into the superpartner of a gauge boson (gaugino) leading to lepton(s) + multi-jets topologies in the final state. The gaugino can decay to a lepton ($`e`$ or $`\nu _e`$) and two quarks with the same $`\lambda ^{}`$ coupling. It also can decay to a lighter gaugino and two fermions. H1 and ZEUS found no evidence for squark production. The results translate to constraints in SUSY parameter space in various scenarios.

Limits on R-Parity-violating couplings have been derived in the unconstrained MSSM where the sfermion and gaugino sector are assumed to be unrelated. The sfermion masses are free parameters. A scan in the SUSY parameter space ($`M_2`$,$`\mu `$), $`M_2`$ being the soft SUSY breaking mass term and $`\mu `$ the mixing mass term for 2 Higgs doublets, has been performed for two different values of $`\mathrm{tan}\beta `$, the ratio of the vacuum expectation values of the two neutral scalar Higgs fields. As a result, Figure 2a presents the ZEUS upper limits on $`\lambda _{1j1}^{}`$ as a function of the squark mass . They are widely SUSY parameter independent. For $`\stackrel{~}{c}`$ and $`\stackrel{~}{t}`$ production the constraints on $`\lambda _{1j1}^{}`$ are more stringent than the indirect bounds from atomic parity violation measurements .

H1 also sets limits  in constrained models where the sfermion and gaugino sector are related at the GUT scale and evolved down to the electroweak scale by the Renormalization Group Equations. Minimal Supergravity (mSUGRA) models assume, in addition, that the electroweak symmetry breaking is driven by radiative corrections. This leads to 5 independent parameters: $`\mathrm{tan}\beta `$, the sign of $`\mu `$, the common scalar mass $`m_0`$, the common gaugino mass $`m_{1/2}`$ and the common trilinear coupling at the GUT scale $`A_0`$. Figure 2b shows the H1 limits for fixed values of $`\lambda ^{}`$ and $`\mathrm{tan}\beta `$ in the ($`m_0`$,$`m_{1/2}`$) plane for $`\stackrel{~}{u}`$, $`\stackrel{~}{c}`$ production. Here, for reasonably large values of the Yukawa coupling (even as small as 0.1), H1 probes a domain which is not ruled out by the ($`\lambda ^{}`$ independent) TeVatron limits. For $`\lambda _{1j1}^{}`$ values of the electromagnetic strength ($`0.3`$), squark masses below 260 GeV are ruled out. However, it should be noted that LEP has the highest sensitivity for $`\stackrel{~}{u}`$, $`\stackrel{~}{c}`$ production, while for $`\stackrel{~}{t}`$ production for intermediate values of $`m_0`$ the H1 sensitivity is comparable with LEP. Note that the results are from 1994–97 data, and there are further data taken since 1998 at a higher center-of-mass energy, 318 GeV. There are also decays not yet investigated, for example $`\stackrel{~}{t}\stackrel{~}{b}+W^+`$, which can produce an event topology with an isolated lepton and missing transverse momentum, accompanied by a hadronic system with large transverse momentum. It is noteworthy that H1 observes an excess above the SM prediction in this topology .

## 4 Summary and prospects

In the HERA searches for R-Parity-violating SUSY and for excited fermions, no evidence for physics beyond the SM has been found. New constraints on squark production in R-Parity-violating SUSY and for excitations of electrons, neutrinos and quarks have been obtained. The HERA limits complement and are competitive with those of the LEP and TeVatron searches. After the HERA detector and luminosity upgrade in autumn 2001, HERA2 will accumulate approximately 1 $`\mathrm{fb}^1`$ of luminosity. Thus, a remarkably enhanced discovery potential will be available for the analyses presented here. In particular, the polarized $`e^\pm `$ beams will lead to an increased squark production cross section and sensitivity.

## 5 Acknowledgements

I wish to thank my colleagues from H1 and ZEUS collaborations for their help in preparing this talk.