## Abstract

We derive, in the form of coupled partial differential equations, the evolution equations for the epitaxial growth, via step flow, of a multispecies crystal on a stepped surface. Both adsorption–desorption on the terraces and attachment–detachment along the step edges are accompanied by chemical reactions and adatom diffusion. Moreover, we account for deposition from either a vacuum, e.g., in molecular beam epitaxy, or a gas, e.g., during vapour phase epitaxy (chemical or physical). Our theory (i) endows the steps with a *thermodynamic* structure whose main ingredients are a free-energy density and species edge chemical potentials, (ii) incorporates *anisotropy* into the terrace species diffusion as well as into the edge free energy, species mobilities, attachment–detachment and reaction-rate coefficients, (iii) allows for large departures from local equilibrium along the steps, and (iv) ensures the consistency of the constitutive relations for the terrace and edge chemical rates with the second law. In particular, a configurational force balance at each step yields a generalization of the classical Gibbs–Thomson relation. Finally, the special case of steady-state growth of a *binary compound* is discussed.

## 1. Introduction

A central challenge facing condensed-matter physicists, materials scientists and applied mathematicians concerns the methods for controlling the growth of single crystals (for a general introduction, cf., e.g., Saito 1996; Pimpinelli & Villain 1998). In this regard, it is well accepted that step flow, i.e. the propagation of steps on a vicinal surface at sufficiently high temperatures but below the roughening transition, constitutes the ideal growth regime since it avoids the complications of island nucleation and coalescence. In particular, recent experimental and theoretical efforts have focused on the formation of one-dimensional structures, such as quantum wires, as a result of morphological instabilities triggered during the epitaxy of self-assembling materials (for a recent review of the available literature on crystal growth and related instabilities, cf. Krug 2004).

The purpose of the present work is to contribute to a better mathematical understanding of the mechanisms, both physical and chemical, underlying the evolution of steps on a vicinal surface during *multispecies* epitaxial growth.1 The step-flow regime has been thoroughly investigated since the celebrated work of Burton *et al*. (1951). The limitations of this so-called BCF model, namely the assumption that steps act as perfect adatom sinks and the confinement to special geometries, specifically flat steps or spirals, have been relaxed in subsequent extensions to account for finite attachment–detachment kinetics and curved steps endowed with a thermodynamic structure (Bales & Zangwill 1990; Chen *et al.* 2001; Caflisch & Li 2003; Otto *et al*. 2004; see Van der Eerden 1993 for a detailed account of crystal growth mechanisms).

Owing to the inherent complexity of epitaxial growth as embodied in the presence of a multitude of competing mechanisms (e.g. terrace and edge species diffusion, terrace adsorption–desorption, step attachment–detachment, line-tension effects and elastic and entropic step–step interactions) and a spectrum of time and length-scales (from atomistic to continuum), most theoretical studies have focused on *single-component* systems.2 In contrast, we shall be concerned with *multispecies* alloys and compounds, e.g., SiGe, GaAs and YBa_{2}Cu_{3}O_{7−δ}. For such materials, the coupling between species deposition, diffusion and attachment kinetics, compounded by the presence of *chemical reactions*, renders the evolution equations significantly more intricate in structure and therefore less amenable to analysis and simulation.

Of particular relevance to our work is that of Pimpinelli & Videcoq (2000) who, in their ‘minimal’ model for a binary compound, show that chemical interactions can lead to morphological instabilities otherwise absent from single-species systems. Specifically, whereas the Ehrlich–Schwoebel barrier is stabilizing during the deposition of a single-species crystal, it can, in the presence of two chemically coupled species, yield instabilities of the kind that lead to step-bunching, i.e. the formation of macrosteps consisting of multiple monatomic steps and separated by wide terraces.3 Here we derive a thermodynamically consistent general theory that relaxes the simplifying assumptions underlying the models of Pimpinelli & Videcoq (2000) and Pimpinelli *et al*. (2003). Specifically, our framework (i) endows the steps with a thermodynamic structure whose main ingredients are a free energy and species chemical potentials, (ii) incorporates anisotropy into the terrace species diffusion and into the edge free energy, species mobilities, attachment–detachment coefficients and reaction rates, thus accounting for the material symmetry of the underlying crystal, (iii) allows for steps to evolve far from equilibrium, with the standard Gibbs–Thomson relation resulting from a restriction to local equilibrium in a single-component system and, importantly, (iv) accounts for the presence of multiple chemically reacting (atomic or molecular) species in a manner consistent with the second law.

Note that continuum models are numerically more attractive than their computationally expensive atomistic counteparts. Moreover, a theory that incorporates anisotropy is useful; Politi & Krug (2000) and Mysliveček *et al*. (2002) have shown that strongly anisotropic terrace diffusion and edge diffusion along curved steps can each cause step bunching. In addition, a framework that accounts for chemical reactions and couples nonlinearly the various species, would provide guidance towards an analytical and numerical understanding of the relevance of these effects on multispecies epitaxy.

The formulation of the evolution equations (in the form of *coupled* partial differential equations for each step and its adjacent upper and lower terraces) is based on mass balances for all chemical species present on the terraces and or along the step edges, an energy balance and a configurational force balance, in conjunction with a dissipation inequality which ensures the compatibility of the constitutive relations with the second law.4

For a system of *N* chemically reacting species, let * n* be the unit normal to the step pointing into the lower terrace (figure 1),

*s*the arc length along the step and

*ϑ*the angle between

*and a high-symmetry axis,*

**n***V*the normal velocity of the step and

*κ*its curvature,

*g*

^{±}the limiting values of

*g*, an arbitrary field defined on both the lower and upper terraces, as the step is approached from the former and the latter respectively, and 〚

*g*〛=

*g*

^{+}−

*g*

^{−}. Also, for a field

*φ*defined at the step, denotes its normal time derivative following the step. The equations that govern the evolution of a step consist of 2

*N*‘flux-matching’ conditionsaugmented by

*N*species equationsand a normal configurational force balancewhere summation over repeated indices is implied (1≤

*k*≤

*N*). Here, for the

*i*th species (1≤i≤

*N*),

*ϱ*

_{i},

*δ*

_{i}and

*ϱ*

_{i}

^{b}denote the terrace, step and bulk densities and

*μ*

_{i}and

*μ*

_{i}

^{s}the terrace and step chemical potentials. Additionally,

*ω*is the terrace grand canonical potential,

*Ψ*

^{b}and

*ψ*are the bulk and edge free energy (per unit area and length, respectively) and

*R*

_{i}

^{c}is the rate of chemical production of

*i*adatoms along the step. Finally,

**D**_{ik}denote the 2×2 binary mobility matrices associated with the (

*i*,

*k*) pair of

*terrace*species, whereas

*D*

_{ik}

^{s}and

*α*

_{ik}

^{±}are scalar mobility and attachment–detachment coefficients for the (

*i*,

*k*) pair of

*edge*species. Finally,

*β*is the step kinetic modulus. These step evolution equations are also coupled to

*N*reaction–diffusion PDEs on each of the terraceswith

*r*

_{i}

^{v}the rate of deposition of particles of the

*i*th species on to the terraces from a vacuum or vapour phase and

*r*

_{i}

^{c}is the rate of production of

*i*-adatoms by chemical reactions on the terraces.

Importantly, when one assumes that the mixtures of terrace and edge adatoms behave like two- and one-dimensional ideal lattice gases, respectively, the constitutive relations for the chemical production rates *r*_{i}^{c}, as given by (5.25) and (5.19), and *R*_{i}^{c}, as prescribed by (5.31) and (5.30), are *consistent with the second law as* expressed by the dissipation inequalities (5.10) and (5.15).

Finally, the system of coupled PDEs above being complicated, a simplification is afforded by considering the special case of a *binary* compound under the further assumptions of steady state and local equilibrium (cf. §6).5

## 2. Notation and basic kinematic identities

In the remainder, we focus on the evolution of a single step separating two adjacent terraces (figure 1). Specifically, the step is viewed as a time-dependent, smooth, simple *curve* =(*t*) in the plane, separating the upper terrace *Ω*^{−}=*Ω*^{−}(*t*) from the lower terrace *Ω*^{+}=*Ω*^{+}(*t*) (both time-dependent open domains in ) and across which jump discontinuities in the terrace, e.g. adatom densities and diffusive fluxes, can occur.6

Given a (scalar, vector or tensor) field *g* on possibly discontinuous at the step, we use the notationWe choose a parametrization of where *s* denotes the arc parameter, such that the upper terrace *Ω*^{−} lies locally to the right of the unit tangent * t* (clockwise parametrization, cf. figure 1). We define the unit normal

*to as the vector obtained by rotating*

**n***by*

**t***π*/2 counterclockwise, so that

*points into the lower terrace*

**n***Ω*

^{+}. Then the curvature

*κ*of is defined by the relation . Additionally, we denote by

*ϑ*the angle between

*and an in-plane reference crystalline axis. Moreover, write for the velocity of the step and denote by*

**n***V*=

*.*

**v***its normal component. Given a smooth time-dependent field*

**n***φ*(

*s*,

*t*) on , the

*normal*time derivative of

*φ*is given by

Finally, let be a portion of the step whose endpoints * x*(

*s*

_{0}(

*t*),

*t*) and

*(*

**x***s*

_{1}(

*t*),

*t*) satisfy

*s*

_{0}(

*t*)<

*s*

_{1}(

*t*) for all

*t*and assume that the motion of is normal, i.e. only the normal component of the velocity

*V*is non-zero. The

*transport theorem*then takes the form (Gurtin 1993)(2.1)where we have introduced the notation

## 3. Step flow equations in the absence of edge diffusion: the single-species case

Although our goal is a theory for the epitaxy of *multicomponent* alloys and compounds, we begin, for the sake of clarity, with the more classical single-species case.

### (a) The classical model

The standard models for the evolution of a step on a vicinal surface are based on three balance equations, namely, two mass balances for adatoms on each of the upper and lower terraces and a ‘global’ balance of mass that involves both terrace and crystallized adatoms. The resulting system of equations is closed by adding suitable constitutive relations for the attachment–detachment fluxes of adatoms at the step.

Specifically, let *ϱ*, * h* and

*r*be the adatom density (per unit area), diffusive flux (per unit length) and supply (per unit area) on the terraces, and denote by

*ϱ*

^{b}the density of bulk atoms (per unit area) in the crystallized upper terrace. Moreover, let

*j*

^{+}and

*j*

^{−}be the scalar fluxes (per unit length) of adatoms incorporated into the step from the lower and upper terrace, respectively. The adatom mass balances on the upper and lower terraces take the form(3.1)while the ‘global’ mass balance is expressed as(3.2)for any two-dimensional domain intersecting the step , with , , and

*the outward unit normal to ∂ (see figure 1). Localization of (3.1) away from the step yields the reaction–diffusion equation for adatoms on the terraces,(3.3)and the jump conditions at the step(3.4)with*

**ν***the unit normal to the step (pointing into the lower terrace), whereas localization of (3.2) at the step gives(3.5)Finally, combining (3.4) and (3.5), we obtain(3.6)The motion of the step can now be determined by assigning explicit forms to*

**n***j*

^{+}and

*j*

^{−}in terms of

*ϱ*

^{+}and

*ϱ*

^{−}. A simple (linear) choice frequently found in the step-flow literature (cf., e.g. Saito 1996) is(3.7)with

*ϱ*

_{eq}>0 the density of adatoms when the step is at equilibrium and

*K*

^{±}≥0 constant attachment–detachment coefficients.7

### (b) Consistency with thermodynamics

The above set of equations for the step does not involve energetic considerations and, in particular, is not necessarily consistent with the second law. Our purpose here is to show that thermodynamics does impose restrictions on constitutive relations such as (3.7) and, in doing so, to discuss the conditions under which these explicit expressions for the adatom attachment–detachment fluxes are compatible with the second law.

Let *Ψ*(*ϱ*) be the terrace adatom free energy and *Ψ*^{b} the free energy of crystallized atoms on the upper terrace (both per unit area). Denote by the terrace adatom chemical potential, by *ω*=*Ψ*−*μϱ* the adatom grand canonical potential and by *μ*^{v} the vapour chemical potential, defined as the amount of energy per unit deposited mass. For isothermal situations such as ours, the first and second law of thermodynamics combine into a single dissipation inequality of the form(3.8)for an arbitrary domain as defined above. This inequality states that the rate at which energy increases in is bounded by the rate at which energy is transported into by *diffusion* across ∂ and that supplied from the vapour through *adsorption–desorption* onto . Localizing (3.8) the step yields the conditionwhich, appealing to (3.4), reduces to(3.9)The local form (3.9) of the second law at the step is not very informative as it contains the quantities *j*^{±} and *V* which are related via the global mass balance. We therefore eliminate *V* in (3.9) by means of (3.6) obtaining(3.10)an inequality that involves only the attachment fluxes. Constitutive relations expressing *j*^{+} and *j*^{−} in terms of *μ*^{+} and *μ*^{−} must be consistent with (3.10). A simple choice consists of(3.11)with *C*^{±}≥0 prescribed constant (attachment–detachment) coefficients.

Before discussing the consistency of (3.7) with (3.10), we note that an alternative procedure may be used to express (3.9) in terms of independent variables. Indeed, (3.6) and (3.9) together are equivalent to the family of inequalities(3.12)with an *arbitrary* scalar field on . To each choice of there corresponds a specific form of the second law at the step, which is automatically consistent with (3.6) and (3.9). Imposition of *local equilibrium* (cf. the paragraph containing (3.22)) then allows the selection of on a physical basis:(3.13)so that (3.12) reduces to(3.14)As before, the simplest choice of constitutive relations for *j*^{+} and *j*^{−} that are consistent with (3.14) is(3.15)with *C*^{±} constant non-negative coefficients. Importantly, using (3.13) to eliminate in (3.15), we obtain (3.11).

In §3*c*, we show that , as determined by (3.13), can be identified with the step chemical potential.

### (c) The step chemical potential

Consider now an external mass supply *R*^{e} of adatoms to the step.8 The partial mass balances at the step remain as given by (3.4), whereas the ‘global’ mass balance now includes the external supply to the step:(3.16)or, appealing to (3.4),(3.17)

As for the second law, in order to account for the energy income due to the external supply, we must add to the right-hand side of (3.8) the term(3.18)with *μ*^{s} the *step chemical potential*, a measure of the energy inflow due to adatom incorporation into the step. Localization of (3.8), as modified by (3.18), along yields, after some simplifications, the inequality(3.19)which has the same form as (3.12) with . However, *μ*^{s} now has a well-defined physical meaning and the status of an additional *unknown* field to be determined together with the solution of the original problem. Once more, the simplest choice of constitutive relations for the attachment–detachment fluxes consistent with (3.19) is(3.20)with *C*^{±}≥0 constant coefficients, together with the *kinetic relation*9(3.21)where *β* is a non-negative constant coefficient (whose inverse can be referred to as the step mobility). In this context, (3.21) may be viewed as a closure condition by which to determine the additional unknown *μ*^{s}, but in §3*d*, we will show that (3.21) is in fact a consequence of the configurational force balance at the step.

We define local equilibrium as corresponding to growth situations in which the dissipation at the step is exclusively due to the attachment–detachment kinetics. This is consistent with the assumption that *β*≡0. Thus the kinetic relation (3.21) reduces to (cf. (3.13))(3.22)which allows us to eliminate *μ*^{s} and yields the constitutive relations (3.11). Hence the introduction of the step chemical potential is consistent with the development in the previous sections.

### (d) An alternative formulation: the kinetic relation and local-equilibrium condition as consequences of the configurational force balance

Since the local equilibrium condition (3.22), or its generalization to the kinetic relation (3.21), plays a central role, it is natural to inquire about its status within the theory. To this purpose, we reformulate the model using the notion of configurational forces, a notion which unifies a broad class of evolution equations for phase fronts and material defects (for a general introduction to the role of configurational forces, cf. Gurtin 2000).

The basic idea is to augment the total mass balance (3.2) with a configurational force balance of the form(3.23)for the arbitrary region defined above, with * C* the terrace configurational stress,

*the terrace configurational internal force and*

**l***the configurational internal force at the step. Localizing (3.23) away from the step, we find(3.24)whereas at the step we obtain(3.25)The configurational internal force*

**g***should be viewed here as indeterminate, since no structural changes occur on the terraces (away from the step). Hence, (3.24) provides no additional field equation for the terraces. Rather it serves to specify*

**l***in a way consistent with (3.23). On the other hand, as we now show, (3.25) yields directly the kinetic relation (3.21).*

**l**The role of the configurational stresses and forces are best clarified when the second law is formulated in terms of *migrating* regions, i.e. we now assume that =(*t*). We denote by * w* the velocity of ∂(

*t*) and rewrite the dissipation inequality in the form(3.26)for every migrating region (

*t*), with

*μ*

^{b}the (indeterminate) chemical potential of the crystalline phase and

_{ext}given by (3.18). Note the presence of (i) the

*flow-like*terms and associated with the motion of ∂(

*t*) and (ii) the

*power-like*term which embodies the notion that the configurational traction performs work over the velocity

*at which ∂(*

**w***t*) evolves.10

Choosing a region that does not contain the step, invariance under the choice of velocity field * w* and the arbitrariness of (

*t*) (cf. Gurtin 1995) yield the

*Eshelby representations*(3.27)with

*ω*=

*Ψ*−

*μϱ*and

*ω*

^{b}=

*Ψ*

^{b}−

*μ*

^{b}

*ϱ*

^{b}the terrace and crystallized

*grand canonical potentials*. Next, localizing (3.26) to the edge, we find thatwhich, upon using (3.4), (3.16) and (3.25), reduces to(3.28)with

*g*=

*.*

**g***, the normal component of the internal configurational force at the step.*

**n**To close the system, it is necessary to specify constitutive relations for *j*^{±} and *g* in terms of (*V*, *μ*^{s}, *μ*^{±}). We choose (see (3.20))with the attachment–detachment coefficients *C*^{±} and the kinetic modulus *β* non-negative constants. With the above constitutive relations, the normal component of the configurational balance (3.25), i.e. , yields the kinetic relation (see (3.21))(3.29)where we have used the Eshelby representations (3.27). Interestingly, (3.29) does not contain the indeterminate quantity *μ*^{b}.

### (e) Compatibility of the classical model with the second law

To compare (3.11) with (3.7), we assume that the adatom free energy *Ψ*(*ϱ*) is *convex* and vanishes at *ϱ*=0, so that the grand canonical potential is negative: *ω*(*ϱ*)<0 for all *ϱ*≠0.11 Next, we define the adatom *equilibrium density* *ϱ*_{eq} at the edge as the (unique) density for which , and we expand *μ* and *ω* about *ϱ*_{eq}:where we have used the definition of the grand canonical potential, *ω*=*Ψ*−*μϱ*, that of the adatom chemical potential, , andInserting these relations into (3.11) we obtain(3.30)with .

Finally, assume that the adatoms are sufficiently *rarified* that the approximationholds, so that the term can be neglected. Thus, we obtain the classical constitutive relations (3.7) as a first-order approximation of (3.11).

## 4. The single-species case with edge diffusion and anisotropic step free energy

We now extend the approach outlined above to encompass adatom diffusion along a step endowed with a free energy density (per unit length). Throughout this and the following sections we *neglect external supplies* since the step chemical potential can now be associated with (in the sense that its derivative with respect to arc length is *conjugate* to) the diffusive flux along the step.

### (a) Mass and configurational force balances

Let *δ* be the adatom density (per unit length) and *k*=* k*.

*the (tangential) diffusive flux along the step. The ‘partial’ mass balances for adatoms on the upper and lower terraces still have the local form (3.3) but the ‘global’ mass balance (3.2) is now modified to account for edge diffusionfor the arbitrary*

**t***time-dependent*two-dimensional control volume (

*t*) as defined earlier. Localizing to the step, and appealing to (2.1) and (3.4), we obtain(4.1)Moreover, since the step is now endowed with an internal (energetic) structure, we introduce a

*step configurational traction*

*, and modify the configurational force balance (3.23) by addition on the left-hand side of the term . Subsequently, localization at the step yieldsthe normal component of which reduces to(4.2)where*

**c***g*=

*.*

**n***and*

**g***=*

**c***γ*

*+*

**t***σ*

*, with*

**n***γ*usually referred to as the line tension and

*σ*the (scalar) configurational shear.

### (b) The second law

We now write the dissipation inequality in the formShrinking the region (*t*) to the edge and assuming that the evolution of is *normal*, we find(4.3)Next, since (*t*), and by consequence ∂*Γ*(*t*) and , may be arbitrarily chosen, we obtain the *identity between the line tension and the grand canonical potential of the step*, , which, in turn, implies the *Eshelby representation*,for the edge configurational traction. Finally, appealing to the adatom flux-matching conditions (3.4), the configurational force balance (4.2), the mass balance (4.1) and the kinematic identity (cf., e.g. Gurtin 1993), the dissipation inequality (4.3) reduces to(4.4)with *G* the thermodynamic driving force at the step, given by

### (c) Constitutive theory

For specificity, we assume that each terrace behaves like a two-dimensional ideal lattice gas, i.e. the terrace-adatom free energy density has the formwith *μ*^{*} and *Γ* positive constants. Hence, the terrace-adatom chemical potential is given bywith *μ*_{0}=*μ*^{*}+*Γ*, and the terrace grand canonical potential reduces toIn the same vein, we assume that the step behaves like a one-dimensional anisotropic ideal lattice gas so that the step free energy is given bywith and *Λ*(ϑ) positive functions of the step orientation. Thus, the step chemical potential reduces towith , whereas the step configurational shear is given byand the step grand canonical potential takes the formFinally, the attachment–detachment fluxes, the (thermodynamic) driving force and the diffusive flux along the step are assumed to be given bywith the scalar coefficients *α*^{±}, *β* and *D*^{s} positive functions of *ϑ*, and *s* the arc parameter along the step.12

### (d) Step evolution equations and the Gibbs–Thomson relation

The evolution equations for the step consist of the ‘partial’ mass balances(4.5)and the ‘global’ mass balance(4.6)augmented by the normal configurational force balance:(4.7)Importantly, the assumption of local equilibrium, i.e. when *β*=0, (4.7) yields a *generalized* Gibbs–Thomson relation(4.8)This condition reduces to the *standard* Gibbs–Thomson when the step density *δ* is neglected and the adatom densities are small with respect to *ϱ*^{b}, as in the paragraph containing (3.30). Indeed, formally take *δ*=0 in (4.8), assume that and view *μ*^{s} as an indeterminate field to be specified by the local equilibrium condition. Hence, within the limit of rarefied terrace adatoms, i.e. when *ϱ*_{eq}≪*ϱ*^{b}, the same reasoning as in §3 gives(4.9)which yields first-order approximations for *j*^{±} analogous (modulo anisotropy) to the constitutive relations of Bales & Zangwill (1990):(4.10)with .

## 5. The multispecies case with edge diffusion and anisotropic step free energy. The role of chemical reactions

Here, we extend the framework established in §4 to account for the presence of a multitude of chemically reacting atomic (molecular) species.13

### (a) Balance laws

Consider a vicinal surface consisting of *N* species and let *ϱ*_{i} and **h**_{i} be the density and diffusive flux of adatoms of the *i*th species on the terraces. The species equations have the form (1≤*i*≤*N*)(5.1)with *r*_{i}^{v} the rate of adsorption–desorption of particles of the *i*th species from the vapour onto the terraces and *r*_{i}^{c} the rate of production–consumption of *i*-adatoms due to the chemical reactions that occur on the terraces (both per unit area), whereas(5.2)where *j*_{i}^{±} are the attachment–detachment (scalar) fluxes of *i*-adatoms from the lower and upper terraces, respectively, and(5.3)with *δ*_{i} and *k*_{i}=**k**_{i}.* t* the edge

*i*-adatom density and (scalar) diffusive flux,

*ϱ*

_{i}

^{b}the bulk density of crystallized particles of the

*i*th species on the upper terrace and

*R*

_{i}

^{c}the rate of production–consumption of edge

*i*-adatoms by chemical reactions occurring

*along the step*(per unit length).

Finally, the normal step configurational force balance remains unaffected (cf. (4.2)).

### (b) Dissipation inequality

Consistent with the general setting in the presence of configurational forces, we write the dissipation inequality for arbitrary migrating control surfaces (*t*):(5.4)where (and hereafter) summation over repeated indices is implied, and, for the *i*th species, *μ*_{i}, *μ*_{i}^{b}, *μ*_{i}^{v}, and *μ*_{i}^{s} are the terrace, bulk, vapour and edge chemical potentials.

The same arguments as before yield the generalized Eshelby representations * C*=

*ω*

**1**on the lower terrace,

*=(*

**C***ω*+

*ω*

^{b})

**1**on the upper terrace, and

*=*

**c***σ*

*+*

**n***ω*

^{s}

*along the step with*

**t***ω*=

*Ψ*−

*μ*

_{i}

*ϱ*

_{i}, , and the terrace, step, and crystallized grand canonical potentials, respectively.

Away from the step, localization of (5.4) yields(5.5)whereas, along the step, we obtain(5.6)with .

### (c) Constitutive theory

Assume that the terrace free energy depends on the terrace species densities:with . From rather standard variational arguments (cf. Fried & Gurtin 2003, 2004), it follows that(5.7)and the dissipation inequality (5.5) reduces to(5.8)The constitutive relations(5.9)where the *N*×*N* matrix of 2×2 mobility matrices **D**_{ik} and the *N*×*N* matrix of scalar adsorption–desorption coefficients *γ*_{ik} are both *positive semi-definite*, and satisfy the inequality (sum over repeated indices)(5.10)are *sufficient* for the second law as expressed in (5.8) to be satisfied.14

Finally, let the free energy density of the step depend on its chemical composition and orientation:with . As before, an extension of the Coleman–Noll procedure yields(5.11)Hence, the step dissipation inequality (5.6) takes the reduced form(5.12)The constitutive relations(5.13)augmented by(5.14)where the orientation-dependent species edge mobilities *D*_{ik}^{s} and attachment–detachment coefficients *α*_{ik}^{±} each form an *N*×*N* positive semi-definite matrix, *β* is a non-negative kinetic modulus, and satisfy the inequality(5.15)and are *sufficient* to ensure consistency with the second law (5.12).

### (d) The terrace adatoms as a two-dimensional ideal lattice gas: thermodynamic compatibility

Assuming that the terrace behaves like an ideal lattice gas (in the dilute limit), its free energy density is given by(5.16)with and *Γ* constant. By (5.7),(5.17)with .

Moreover, we consider situations in which *M* chemical reactions (involving only terrace adatoms) occur on any individual terrace. The *α*-reaction (1≤*α*≤*M*) is such that reactant species _{i} combine to form product species _{j} according to the reversible mechanism:(5.18)where the positive integer constant denotes the stochiometric coefficient of the *i*th reactant (*j*th product) species. Associated with the forward (backward) reaction is a *positive* scalar coefficient such that the rate of the *α*-reaction is given by(5.19)Appealing to (5.17), (5.19) can be rewritten as(5.20)

Assume now that chemical equilibrium can be attained by the reaction (5.18), i.e. there exist equilibrium chemical potentials such that(5.21)It is straightforward to check that a necessary and sufficient condition for chemical equilibrium to be attained or, equivalently, that there exist solutions to (5.21), is that(5.22)When substituted into (5.20), (5.22) yields the following expression for the *α*-reaction rate:(5.23)where(5.24)

Finally, the chemical production rate of *i*-adatoms on the (upper or lower) terrace can be expressed as:(5.25)with *r*^{α} given by (5.23). Thus, by (5.24), the terrace chemical dissipation takes the formso that, appealing to (5.23), the inequality (5.10) reduces to(5.26)It is trivial to see that this inequality is automatically satisfied, therefore ensuring the consistency, when the terrace free-energy density is given by (5.16), of the constitutive relations (5.25) and (5.19) with the second law.

### (e) The step as a chemically reacting one-dimensional anisotropic ideal lattice gas: consistency with the second law

Assuming that the mixture of edge adatoms behaves like an *anisotropic* one-dimensional ideal lattice gas (in the dilute limit), the dependence of the step free-energy density on its chemical composition and orientation is given by(5.27)Appealing to (5.11)_{1}, it follows that the edge-adatom chemical potentials and the step configurational shear take the form(5.28)with .

Let *K* be the number of chemical reactions that occur along the step (and involve only edge adatoms) and let the *β*-reaction be such a reaction(5.29)Denote by *R*^{β} the corresponding reaction rate:(5.30)Proceeding in a manner analogous to that leading to (5.25), it is easy to show that the chemical production rate of edge *i*-adatoms can be reduced to(5.31)where *R*^{β} is given by(5.32)with(5.33)A reasoning identical to the one of the previous section ensures that the constitutive relations (5.31), (5.32) and (5.33) are *sufficient* for the edge dissipation inequality (5.15) to be satisfied.

## 6. The case of a binary system

Assume that only two species of adatoms, _{1}= and _{2}=, are deposited on the terraces. After diffusion and subsequent attachment, they react along the *isotropic* step to form the species _{3}= (of which the bulk crystalline phase is constituted), according to the reversible chemical reaction(6.1)and let *λ*_{b} be the rate coefficient associated with the backward mechanism in (6.1). Moreover, (i) assume that adatom diffusion on the terraces proceeds within the steady-state regime, (ii) neglect edge adatom densities and diffusion and (iii) approximate (5.9)_{1} with the dilute-limit approximation of Fick's law, that is, and . Then, the terrace diffusion (5.1) reduces to(6.2)with Δ being the Laplacian in , and and the deposition fluxes.15 Additionally, consistent with the steady-state assumption, we neglect the convective terms and in the ‘partial’ species equations along the step (5.2)(6.3)with and being the normal derivatives of and , respectively, as the step is approached from the lower and upper terraces, respectively, and the attachment–detachment fluxes and are given by the relations (see (5.13))(6.4)where we have implicitly assumed that the - and - adatoms along the edge are sufficiently dilute that the species attachment–detachment kinetics along the step are *uncoupled*. The ‘global’ species (5.3) reduce to(6.5)with *ϱ*^{b} the bulk density, and , and the rates of production of , and edge adatoms, respectively, due to the step chemical reaction (6.1). Appealing to (5.31) and (5.32), we have(6.6)Finally, closure is provided by the configurational balance (5.14)_{1}, which, under the further assumption of *local equilibrium*, reduces to(6.7)

In summary, the basic unknowns of the problem are the adatoms chemical potentials and defined on the terraces together with the step chemical potentials , and , and the step normal velocity *V*. The corresponding free-boundary problem consists of the diffusion (6.2) on the upper and lower terraces, the ‘flux-matching’ conditions (6.3), the ‘global’ species mass balances (6.5), and the generalized Gibbs–Thomson condition (6.7), augmented by the constitutive relations (6.4) and(6.6).

## Acknowledgments

P.C. was supported by the Italian M.I.U.R. project ‘Modelli matematici per la scienza dei materiali’ and M.J. was supported by the NSF under grant no. DMS-0204939 and the KSEF under grant no. 148-502-02-15.

## Footnotes

↵Epitaxy designates growth during which the film inherits, at least in the early stages, the crystalline structure of the substrate. If the film and substrate are made of distinct materials, the mismatch between their lattice parameters generates a misfit stress in the film. Here the resulting strain is ignored as the emphasis is on the interplay between diffusion and chemistry. During step flow, adsorption–desorption occurs mainly on the terraces, followed by the diffusion of adatoms, i.e. adsorbed atoms or molecules, until the steps are reached where attachment–detachment and edge diffusion take place. Additionally, for multicomponent systems, chemical reactions are present, both on the terraces and along the steps, with the latter leading to the incorporation of particles into the bulk, by which the step edges evolve laterally.

↵Notable exceptions include Jabbour & Bhattacharya (1999), Pimpinelli & Videcoq (2000), Vladimirova

*et al*. (2000) and Pimpinelli*et al*. (2003).↵As pointed out by these authors, the possibility of step bunching in binary materials offers an alternative mechanism for the development of instabilities that does not recourse to the controversial assumption of an

*inverse*Ehrlich–Schwoebel barrier (cf. Vladimirova*et al*. 2000; Krug in press). See also Jabbour & Bhattacharya (1999), who show that, in the presence of*competition*between the vapour-phase precursors for open adsorption sites, even slight deviations from vapour stochiometry can lead to a drastic drop in the growth rate during chemical vapour deposition.↵Within continuum physics, configurational forces are associated with the evolution of grain boundaries, phase interfaces, film surfaces and disclinations in liquid crystals, etc. (cf. Gurtin 2000; Cermelli & Fried 2002; Fried & Gurtin 2003, 2004; Jabbour & Bhattacharya 2003). In the context of epitaxy on vicinal surfaces, these forces can be viewed as

*driving the lateral motion of steps*.↵The steady-state approximation is valid for growth situations in which the time-scale associated with the motion of steps is negligible in comparison with the time-scale for diffusion of adatoms on terraces.

↵Note that, in accordance with the above continuum description, the length-scale considered herein, although submicroscopic, is large enough that a distinction between kinks, corners and straight segments along a step is unwarranted. Instead, it is implicit that attachment of adatoms to the steps occurs preferentially at kinks.

↵The classical BCF model is based on the assumption of

*infinite*attachment–detachment kinetics, i.e.*K*^{±}→∞, so that the step acts as a perfect sink:The step velocity is then determined by (3.5) which now reduces towith*ϱ*^{b}assumed constant and the adatom flux given by a*linearized*Fick's law . Finally, note that*j*^{±}are now constitutively indeterminate, and are specified a posteriori by (3.4) once the adatom densities and step velocity have been resolved.↵External supplies can be assigned

*arbitrarily*and, although this may seem artificial, the availability of these supplies provides us with a useful conceptual tool to investigate the constitutive dependence of*j*^{±}on the variables of the theory.↵See, for example, Heidug & Lehner (1985), Truskinowsky (1987) and Abeyaratne & Knowles (1990).

↵The terrace and edge configurational forces

and**l**do not expend power on (**g***t*) as they act*internally*to the migrating control surface.↵This convexity condition is automatically satisfied for single- and multi-component ideal lattice gases (see §§4 and 5).

↵To illustrate the anisotropy of the step, assume that lies on the [001]-plane of a

*cubic*crystal. Hence, the step free-energy density has square symmetry, i.e. , and , where , and are smooth periodic functions with period 2*π*. Additionally, we require that , and have minima at the high-symmetry directions*ϑ*=0,*π*/2,*π*and 3*π*/4 (cf. Van der Eerden 1993). Analogous constitutive relations can be postulated for*α*^{±},*β*and*D*^{s}.↵We still allow for species diffusion along the step which, as before, is endowed with an anisotropic free-energy density.

↵The positive-definiteness of the matrix of matrices

**D**_{ik}signifies that for all . (An analogous definition holds for the matrix of scalar coefficients*α*_{ik}.) Note that (5.9)_{1}has the classical Fickean form for multicomponent systems, while (5.9)_{2}generalizes the standard expression for the deposition rate to the case of multispecies crystals. If the diffusion of adatoms is*isotropic*for all species, then the mobility matrices have the form**D**_{ik}=*D*_{ik}**1**. The explicit form of (5.9)_{3}will be discussed more explicitly below in the context of ideal lattice gas.↵In the context of MBE, and are assumed constant. This is consistent with the second law (5.8). Indeed, when adsorption of

*i*-particles occurs,*F*_{i}>0 and so that .- Received October 20, 2004.
- Accepted April 5, 2005.

- © 2005 The Royal Society