singularity of the surface free energy at the high symmetry crystalline orientation. There are ... The BCF theory assumes that mass transfer between t...

0 downloads 3 Views 793KB Size

arXiv:cond-mat/9905031v1 [cond-mat.mtrl-sci] 4 May 1999

Navot Israeli∗ and Daniel Kandel∗∗

Department of Physics of Complex Systems, Weizmann Institute of Science, Rehovot 76100, Israel The decay of a crystalline cone below the roughening transition is studied. We consider local mass transport through surface diffusion, focusing on the two cases of diffusion limited and attachment-detachment limited step kinetics. In both cases, we describe the decay kinetics in terms of step flow models. Numerical simulations of the models indicate that in the attachmentdetachment limited case the system undergoes a step bunching instability if the repulsive interactions between steps are weak. Such an instability does not occur in the diffusion limited case. In stable cases the height profile, h(r, t), is flat at radii r < R(t) ∼ t1/4 . Outside this flat region the height profile obeys the scaling scenario ∂h/∂r = F(rt−1/4 ). A scaling ansatz for the timedependent profile of the cone yields analytical values for the scaling exponents and a differential equation for the scaling function. In the long time limit this equation provides an exact description of the discrete step dynamics. It admits a family of solutions and the mechanism responsible for the selection of a unique scaling function is discussed in detail. Finally we generalize the model and consider permeable steps by allowing direct adatom hops between neighboring terraces. We argue that step permeability does not change the scaling behavior of the system, and its only effect is a renormalization of some of the parameters. 68.35.Bs, 68.55.-a, 68.55.Jk

I. INTRODUCTION

The properties of crystalline nanostructures are of considerable interest, because of the technological importance of nanostructures in fabrication of electronic devices. Kinetic properties of nanostructures attracted particular attention, since in many cases nanostructures are thermodynamically unstable and tend to decay with time. Such decay processes have been studied both theoretically and experimentally. Considerable effort was devoted to the study of periodic structures. The decay of one and two dimensional gratings was studied extensively. The emerging experimental picture3–6 is that below the roughening temperature these structures decay in a shape preserving manner. Macroscopic facets are observed at the maxima and minima of the gratings. Although these systems are out of equilibrium, the appearance of facets is a manifestation of the cusp singularity of the surface free energy at the high symmetry crystalline orientation. There are basically two theoretical approaches to the problem of surface evolution in general and nanostructure decay in particular. On the one hand, there are phenomenological models which treat the crystal surface as a continuous medium9–14 . The evolution of the surface is then driven by the tendency of the system to lower its free energy (given in terms of continuous spatial variables). The advantage of such models is that they are relatively simple and can sometimes lead to analytical predictions of surface evolution. However, these models ignore the discrete nature of surface steps, which may become important below the roughening transition. In addition, most of them rely on assumptions of small surface slope and/or surface curvature (with the exception of Ref. 12), and are unable to properly treat the behavior of the macroscopic facets observed experimentally. On the other hand, there are models which treat surface evolution on a smaller scale. Among these are microscopic models21–24 , where the basic degrees of freedom are individual atoms, and step flow models14–18 . These models are usually solved numerically, and provide results which can be directly related to the microscopic dynamics. However, in most cases,

1

it is difficult to understand the behavior of the system on larger length scales on the basis of these results. Research efforts were also directed towards isolated surface structures. The decay of isolated step bunches, islands or hills was studied both experimentally7,8 and theoretically19,20,25 . Here too there is an apparent gap between the microscopic and macroscopic theoretical approaches. In this work we attempt to bridge over this gap in the case of a simple surface structure, i.e. an infinite crystalline cone. We give a complete account of the surface dynamics based on a step flow model, and then derive a continuum model which gives a very accurate description of the evolution of the cone and becomes exact in the long time limit. We do not make any assumptions of small slope or small surface curvature. The crystalline cone consists of an infinite number of circular concentric steps. A similar system was studied by Rettori and Villain15 who considered the decay of bidirectional surface modulations. Their results are relevant in the case of small amplitude modulations when the profile peaks and valleys affect each other. Our work addresses the opposite situation when a single peak can be considered as an isolated structure and in this sense is complimentary to theirs. Below the roughening transition atomic steps have a finite free energy. Their existence on the surface strongly affects its morphological evolution. In many cases, one can ignore the formation of islands and voids on the surface and consider only adatom diffusion and attachment and detachment processes to and from step edges. The decay of a nanostructure is then dominated by the motion of steps. In order to describe the decay process mathematically, one has to solve the diffusion equation for adatoms on the terraces between the steps, with boundary conditions at the step edges, in the spirit of the Burton-Cabrera-Frank model1 . If the geometry of the nanostructure is simple, this procedure leads to a set of coupled equations of motion for the steps. Our goal in the present work is to construct and solve these equations of motion for the simple case of an infinite crystalline cone. A partial account of this work is found in Ref. 2. The kinetic step model for the cone is derived in section II. In section III we carry out numerical simulations of the model and examine the evolution of surface morphology under various conditions. Previous experimental and theoretical research on decay of nanostructures has demonstrated that in most cases the surface reaches a scaling state where the typical length scale depends on time algebraically. Our simulations show that the cone profile also exhibits such a scaling behavior. In section IV we show analytically that the step flow model admits such solutions. We calculate the scaling exponents and derive a continuum equation for the scaling function. The properties of the scaling function are analyzed in sections V and VI. It has been suggested that in some materials, steps may be permeable to flow of atoms; i.e. atoms can hop between neighboring terraces without attaching to the step which separates them. In section VII we discuss the consequences of step permeability on the decay of an infinite cone, and show that permeability does not change the qualitative behavior of the system. II. STEP FLOW MODEL OF A CRYSTALLINE CONE

We now consider a crystalline surface which consists of flat terraces, parallel to a high symmetry plane of the crystal and separated by atomic steps. Islands and voids are ignored. The evolution of such systems of steps was treated long ago by Burton, Cabrera and Frank1 (BCF). The BCF theory assumes that mass transfer between terraces is governed by diffusion of adatoms on the terraces. These atoms are emitted and absorbed at step edges, which according to BCF act as perfect sinks and sources. This last assumption, however, is valid only when adatom diffusion is the rate limiting process. Modifications of the BCF model account for the finite rate of adatom attachment-detachment processes at step edges. The BCF model was also generalized to include elastic and entropic interactions between steps (see for example Refs. 14,26). In this section we construct such a generalized BCF model for the morphological evolution of a conic hill on a crystal surface. 2

Consider the surface of an infinite crystalline cone, which consists of circular concentric steps of radii ri (t), separated by flat terraces (Fig. 1). The index i grows in the direction away from the center of the cone. r

1

r

2

r

3

FIG. 1. Illustration of the step configuration in a crystalline cone. The steps height is exaggerated.

These steps may absorb or emit atoms, which then diffuse across the neighboring terraces with a diffusion constant Ds . Assuming no deposition of new material, no evaporation and no transport through the bulk, the adatom concentration Ci (~r) on the ith terrace satisfies the diffusion equation ∂Ci (~r) , ∂t where ~r is a two dimensional vector parallel to the high symmetry terraces. In most situations, the time scale associated with step motion is much larger than the time scale of surface diffusion. One can therefore assume that the adatom diffusion field is always in a steady state; i.e., for any step configuration, the diffusion field reaches a steady state before the steps move significantly. Within this quasi-static approximation, the r.h.s. of the diffusion equation can be neglected. Using the radial symmetry of the cone, we can write the static diffusion equation in polar coordinates as: Ds ∇2 Ci (~r) =

∂ 2 Ci (r) 1 ∂Ci (r) + =0. ∂r 2 r ∂r The general solution of this equation is Ci (r) = Ai ln r + Bi .

(1)

(2)

The coefficients Ai and Bi are determined by the boundary conditions at the step edges. To define these conditions, we assume that the flux of atoms at the two step edges bounding the ith terrace is determined by first order kinetics, characterized by an attachment-detachment rate coefficient k: ∂Ci eq Ds = k C i | ri − C i ∂r ri

∂Ci eq = k C −Ds i |ri+1 − Ci+1 , ∂r ri+1

3

(3)

where Cieq is the equilibrium concentration of atoms on the terrace adjacent to the ith step. Using these boundary conditions, we find that the constants Ai are given by Ai =

eq Cieq − Ci+1

ln ri − ln ri+1 −

Ds k

1 ri

+

1 ri+1

.

(4)

Employing mass conservation at the step we obtain the step velocity: dri = ΩDs dt

∂Ci ∂Ci−1 − ∂r ∂r

!

= ΩDs

ri

Ai − Ai−1 , ri

(5)

where Ω is the atomic area of the solid. In order to complete the solution of the diffusion problem we have to determine the equilibrium adatom concentration at step edges. This concentration depends on the local step curvature and on the radii of neighboring steps. According to the Gibbs-Thomson relation, Cieq is given by Cieq = C¯ eq exp

µi µi ≈ C¯ eq 1 + T T

,

(6)

in units where the Boltzmann constant is equal to 1. Here µi is the chemical potential associated with the addition of an atom to the solid at the ith step, T is the temperature and C¯ eq is the adatom equilibrium concentration at the edge of a straight isolated step. To evaluate the step chemical potential, µi , we take into account the step line tension, Γ, and a repulsive interaction between nearest neighbor steps. The magnitude of this interaction is inversely proportional to the square of the distance between the steps at large distances. Such a dependence is consistent with entropic as well as elastic28,29 interactions between straight steps. We follow Ref. 16 and take the interaction energy between step i + 1 and an atomic segment of the ith step to be √ G Ωri+1 U(ri , ri+1 ) = , (7) (ri + ri+1 ) (ri+1 − ri )2 where G is the interaction strength. The chemical potential of the ith step is then given by µi =

ΩΓ √ ∂ [U(ri , ri+1 ) + U(ri , ri−1 )] + Ω . ri ∂ri

(8)

This equation applies also to the chemical potential of the first step if we set r0 = 0. As we show below, in the long time limit the distance between steps is small compared with the step radius. In this limit we can approximate the step chemical potential by ΩΓ 2ri−1 1 1 2ri+1 µi = + ΩG · · 3 − ri ri+1 + ri (ri+1 − ri ) ri + ri−1 (ri − ri−1 )3

!

.

(9)

This approximation simplifies the algebra considerably, and we verified that it does not affect the analytical and numerical results described below. We are now ready to write down the step equations of motion using Eqs. (4), (5) and the Gibbs-Thomson relation (6). It is convenient to use dimensionless radii, ρi , and dimensionless time τ : T · ri , ρi = ΩΓ Ds T −1 T 2 τ = Ds C¯ eq Ω · · 1+ ·t. ΩΓ kΩΓ 4

In terms of these variables, the equations of motion take the form ρ˙ i ≡ ai =

dρi ai − ai−1 = , with dτ ρi ξi − ξi+1

ρi −q (1 − q) ln ρi+1

1 ρi

+

(10) 1 ρi+1

,

1 2ρi−1 1 2ρi+1 1 · · ξi = + g 3 − ρi ρi+1 + ρi (ρi+1 − ρi ) ρi + ρi−1 (ρi − ρi−1 )3

!

.

2

Eqs. (10) depend on two parameters: g = ΩT2 ΓG3 measures the strength of step-step interkΩΓ −1 actions G relative to the line tension Γ, while the parameter q = (1 + D ) determines sT the rate limiting process in the system. When q → 1, diffusion across terraces is fast and the rate limiting process is attachment and detachment of adatoms to and from steps. On the other hand, when q → 0, the steps act as perfect sinks and the rate limiting process is diffusion across terraces. III. RESULTS OF SIMULATIONS

We integrated Eqs. (10) numerically both in the case of diffusion limited kinetics (DL) and in the case of attachment-detachment limited kinetics (ADL). The initial configuration was a uniform step train. In principle, the initial step separation is a parameter of the model. However, for any value of this parameter one can change the units of length and time and get the same equations of motion with initial step separation of unity and different values of the parameters g and q. Thus it is sufficient to consider an initial step separation of unity. When the repulsive interactions between steps are weak (i.e. g is small), there is a striking difference between the dynamics in the DL and ADL limits. In the ADL case the system becomes unstable towards step bunching, whereas in the DL case there is no such instability. However, when g is large enough the instability disappears even in the ADL case. Let us first discuss situations where the step bunching instability does not occur. Fig. 2 shows the time evolution of the system in the ADL and the DL cases with a relatively large value of g. Each line in the figures describes the evolution of the radius of one step. We note that the innermost step shrinks while the other steps expand and absorb the atoms emitted by the first step. When the innermost step disappears, the next step starts shrinking and so on. Our observations indicate that the time at which the nth step disappears, τn , grows with n as τn ∼ n4 . This process results in a propagating front, which leaves a growing plateau or a facet at the center of the cone. At large times, the (dimensionless) position of this front behaves as ρfront (τ ) ∼ τ 1/4 . This is shown by the dashed lines in Fig. 2. This power law is an indication of a much more general and interesting phenomenon. It turns out that at large times, not only the front position but also the positions of minimal and maximal step densities scale as τ 1/4 . In fact, the step density D(ρ, τ ), defined as the inverse step separation, obeys the following scaling scenario: There exist scaling exponents α, β and γ, which define the scaled variables x ≡ ρτ −β/γ and θ ≡ τ 1/γ . In terms of these variables D(ρ, τ ) = θα F (x, θ) ,

(11)

where the scaling function F is a periodic function of θ with some period θ0 . Our ansatz is somewhat weaker than standard scaling hypotheses, which would assume the scaling function F is independent of θ. We introduce this periodic dependence because our simulations strongly indicate a periodic behavior, generated by the motion of the first step (see Fig. 2). Thus the disappearance time of step n is τn = (nθ0 )γ . An immediate consequence of the 5

scaling ansatz is that if we define θ = θ¯ + nθ0 with 0 ≤ θ¯ < θ0 , and plot θ−α D (ρ, τ ) against x, all the data with different values of n and the same value of θ¯ collapse onto a single curve, ¯ F (x, θ). b)

a) 60

60

50

50

40

40

ρ 30

ρ 30

20

20

10

10

0 0

0.5

1

1.5

τ

0 0

2

0.5

1

1.5

τ

5

x 10

2 5

x 10

FIG. 2. Time evolution of the step radii in the (a) ADL and (b) DL cases with g = 0.01. The radius of the facet edge can be fitted by a τ 1/4 power law (dashed lines).

To verify that our system obeys this scaling ansatz, we define the step density at a discrete set of points in the middle of the terraces: ρi (τ ) + ρi+1 (τ ) D ,τ 2

!

=

1 . ρi+1 (τ ) − ρi (τ )

(12)

Fig. 3 shows plots of D(ρ, τ ) as a function of x = ρτ −1/4 for a fixed value of θ¯ and 6 different values of n in the ADL and the DL cases. The excellent data collapse shows that our scaling ansatz indeed holds with α = 0, β = 1 and γ = 4 in both cases. a)

b) 1.6 1.4

1.2

1.2

1

1

F 0.8

F

0.8 0.6

0.6 0.4

0.4

0.2

0.2 0 0

1

2

3

4

5

6

0 0

7

−1/4

1

2

3

4

5

6

7

x=ρτ−1/4

x=ρτ

FIG. 3. Data collapse of the density function in the (a) ADL and (b) DL cases with g = 0.01. The values of the scaling exponents used here are α = 0, β = 1 and γ = 4. These figures show ¯ as a function of x = ρτ −1/4 . density functions with 6 different values of n and the same value of θ, Different symbols corresponds to different values of n. The unscaled data is shown in the insets.

We now turn to discuss effects of interactions between steps. As shown above, the behavior of the cone in the ADL and the DL limits is very similar when the repulsion 6

between steps is strong. We have already mentioned that when the repulsion is weak, the behavior of the system in the DL limit is very different from its behavior in the ADL case. This is not surprising, since linear stability analysis (see Appendix A) of the two cases in the absence of step-step interactions predicts that the ADL case is unstable towards step bunching while the DL case is marginal. In intermediate cases (0 < q < 1) the system is unstable when the interactions are sufficiently weak. Fig. 4 shows the time evolution and the data collapse of the density function in the DL case when step-step interactions are weak. As one can see, the behavior of the cone is qualitatively similar to the large g example (Fig. 2 (b)) and that the scaling ansatz still holds. Quantitatively, the step density near the facet edge is much higher for small values of g. Also the dependence of the scaling function on scaled time within a collapse period is much more pronounced when the interactions are weak. a)

b) 60

4 3.5

50

3 40

2.5

ρ 30

F

2 1.5

20

1 10

0 0

0.5 0.5

1

τ

1.5

0 1.5

2

2

2.5

3

−1/4

5

x 10

x=ρτ

FIG. 4. (a) Time evolution of the step radii in the DL case with g = 10−6 . (b) Data collapse of the density function in the same system. Different symbols correspond to different times.

The dependence of step kinetics on the strength of the interactions is much more complicated in the ADL case. Fig. 5 shows the evolution and when possible, the scaled density function of a series of ADL systems which differ in the value of g. The two quantitative observations we made in the DL case hold also here: The high step density region near the front becomes more dense, and the time dependence of the scaling function becomes more pronounced as the value of g is reduced (Fig. 5 (a)-(d)). In addition, we noted that the approach of the system towards the scaling state is slower. Adjacent to the dense region at the edge of the facet is another region of very low density of steps. This region becomes less dense as the value of g is reduced. Below a critical value, gc , of the interaction parameter, a few steps between the low density region and the facet edge form a bunch. At this stage the scaling ansatz breaks down. The system no longer exhibits a simple periodic nature, which results from the collapse of single steps, one at a time. Instead, it seems to adopt a complicated almost periodic pattern which involves the collapse of many steps in each period (Fig. 5 (e)). Hints for this periodicity changes are already present in Fig. 5 (a) and (c), where the steps crossing the low density region follow a threefold periodicity. If the value of g is further reduced, bunches of steps collapse together, and finally the whole system becomes unstable toward step bunching (Fig. 5 (f)). Neighboring steps merge and form a bunch, which in turn merges with another bunch and so on. Step bunches rather than isolated steps become the dynamic objects in the system.

7

a)

b) 60

2.5

50

2

40 1.5

ρ 30

F 1

20 0.5

10

0 0

0.5

1

τ

1.5

0 0

2

1

2

3

4

5

6

4

5

6

x=ρτ−1/4

5

x 10

c)

d) 60

2.5

50

2

40 1.5

ρ 30

F 1

20 0.5

10

0 0

0.5

1

τ

1.5

0 0

2

1

2

3 −1/4

5

x 10

x=ρτ

e)

f) 60

60

50

50

40

40

ρ 30

ρ 30

20

20

10

10

0 0

0.5

1

τ

1.5

0 0

2 5

x 10

0.5

1

τ

1.5

2 5

x 10

FIG. 5. Time evolution and scaled density functions of ADL systems with different interaction strengths. (a) and (b): g = 10−3 , (c) and (d): g = 5 · 10−4 , (e): g = 2.5 · 10−4 and (f): g = 10−6 . IV. SCALING ANALYSIS AND THE CONTINUUM MODEL

The above results suggest that when the step bunching instability does not occur, the time evolution of the system can be described by a step density function, which is continuous in both position and time variables. In this section we derive such a continuum model by carrying out a scaling analysis. We obtain an equation which governs the scaling function and evaluate the scaling exponents analytically. 8

Motivated by the simulation results, we assume that at long times the scaling ansatz, Eq. (11), holds. This together with conservation of the total volume of the system, already determines the two scaling exponents α and β. First, we derive a relation between these scaling exponents by considering the height profile h(ρ, τ ). Assuming steps of unit height, the height difference between two points on the surface is given by the number of steps between them. The continuous analog of this statement can be used to derive the following relation between the height profile and the step density: h(ρ, τ ) = h0 (τ ) −

Z

0

ρ

D(ρ′ , τ )dρ′ ,

(13)

where h0 (τ ) is the height at the origin. Far enough (in the limit ρ → ∞), h(ρ, τ ) does not change with time. Calculating the height change at infinity using Eq. (13) we find h0 (τ ) − h0 (0) − τ

α+β γ

Z

∞

[F (x, θ) − F (∞, 0)] dx = 0 ,

0

(14)

where we have changed the integration variable and used the definition (Eq. (11)) of the scaling function F . On the other hand, h0 (0) − h0 (τn ) = n because τn is the time of disappearance of the nth step. γ satisfies the relation τn ∼ nγ , and therefore we have h0 (0) − h0 (τn ) ∼ τn1/γ . This and the τ (α+β)/γ dependence in Eq. (14) lead to the relation α + β = 1. In addition, conservation of the total volume of the system, V, implies that V = 2π

Z

0

∞

ρh(ρ, τ )dρ

(15)

is independent of τ . Integration by parts of the derivative of this integral with respect to τ yields the following equation: Z

∞ 0

ρ2

∂D(ρ, τ ) dρ = 0 . ∂τ

(16)

Evaluating this integral in terms of the function F and the scaled variables x and θ we obtain the equation Z

0

∞

x2

"

#

αF (x, θ) ∂F (x, θ) dx = 0 . + θ ∂θ

(17)

This can be satisfied for all θ only if α = 0, since F is a periodic function of θ. Combining this result with the previously obtained relation α + β = 1, we conclude that α=0 , β=1.

(18)

Thus, γ is the only nontrivial scaling exponent in the model. To evaluate γ and the scaling function F , we continue with the equation for the full time derivative of the step density D: dD ∂D ∂D dρ = + · . dτ ∂τ ∂ρ dτ

(19)

Eq. (19) can be evaluated in the middle of the terrace between two steps (i.e. at ρ = (ρi + ρi+1 )/2). The l.h.s. of (19) is calculated by taking the time derivative of Eq. (12): 9

dD/dτ = −D 2 (ρ˙ i+1 − ρ˙ i ). This together with the fact that dρ/dτ = (ρ˙ i + ρ˙ i+1 )/2 leads to the relation ∂D ρ˙ i+1 + ρ˙ i ∂D + + D 2 (ρ˙ i+1 − ρ˙ i ) = 0 . ∂ρ 2 ∂τ

(20)

Now we change variables to θ and xi ≡ ρi θ−1 (since β = 1), and transform Eq. (20) into an equation for the scaling function F : ∂F ∂x

θ

!

θ ∂F + ρ˙ i x + − + F 2 θγ (ρ˙ i+1 − ρ˙ i ) = 0 . 2 γ γ ∂θ

γ−1 ρ˙ i+1

(21)

The step velocities ρ˙ i and ρ˙ i+1 can also be expressed in terms of the xi ’s, but we defer this algebraic manipulation to a later stage. Our next goal is to take a continuum limit of Eq. (21) in the variable x = (xi+1 + xi )/2. Such a continuum limit becomes exact in the long time limit. To see this, let us rewrite Eq. (12) in terms of xi ’s, θ and F : xi+1 − xi =

θ−1 . F ((xi+1 + xi )/2, θ)

(22)

According to this equation, the difference between successive xi ’s is of order θ−1 wherever F does not vanish. In the large θ (long time) limit these differences become vanishingly small. This justifies the continuum limit in the variable x. In practice, we take the continuum limit in the following way. We evaluate the function F at the position (xi+k + xi+k+1 )/2 by using it’s Taylor expansion F

∞ X 1 ∂ n F (x, θ) xi+k + xi+k+1 θ−1 xi+k + xi+k+1 = ,θ ≡ −x 2 xi+k+1 − xi+k n=0 n! ∂xn 2

n

.

(23)

As long as k is finite, the difference xi+k − x is small in the long time limit. It is therefore useful to expand xi+k around x xi+k = x +

∞ X

φkn θ−n ,

(24)

n=1

and insert this expansion into Eq. (23). The resulting equation is expanded as a power series in θ−1 . By requiring that the equation be obeyed to all orders in θ−1 , we can calculate the coefficients φkn for any desired k and n. These coefficients will involve the function F and its derivatives with respect to x, which are all periodic functions of θ. Next, we express the velocities ρ˙ i and ρ˙ i+1 in terms of the scaled radii xi−2 , xi−1 , . . . , xi+3 using Eq. (10) and the transformation to scaled variables. We then use Eq. (24) to expand the velocities in powers of θ−1 . The result of this expansion is ρ˙ i+1 + ρ˙ i = θ−3 (A + O(θ−1)) , ρ˙ i+1 − ρ˙ i = θ−4 (B + O(θ−1)) .

(25)

A and B are known expressions involving F , F ′ , F ′′ , F ′′′ , F ′′′′ , where the primes denote partial derivatives with respect to x. The full expressions are given in Appendix B. The existence of derivatives up to fourth order in this equation is a consequence of the fact that each step “interacts” with four other steps (two on each side) through the equations of 10

motion (10). Inserting Eqs. (25) into Eq. (21), we obtain the following differential equation for F: − F′ ·

x A θ ∂F + θγ−4 F ′ · + F 2 B + · + O(θγ−5 ) = 0 . γ 2 γ ∂θ

(26)

Consider Eq. (26) at large θ. Our expansion in the small parameter θ−1 is valid only at values of x where F does not diverge or vanish (see above). Therefore, the first term in Eq. (26) is O(1). This term has to be canceled by the second term if we require F to satisfy a single differential equation. Hence, we must have γ=4.

(27)

The fourth term vanishes as θ → ∞ since γ − 5 < 0, and the third term must vanish as well. Therefore, in the large θ limit, F is only a function of x, and we are left with an ordinary differential equation for F : F

′

A x + F 2B = 0 . − 2 4

(28)

The detailed form of this equation in the DL and ADL cases is also given in Appendix B. Let us emphasize several important properties of our scaling analysis. First, the values of the scaling exponents we calculated (α = 0, β = 1 and γ = 4) are consistent with the results of numerical simulations (see above). Secondly, our continuum model is valid for arbitrarily large surface curvature and slope (unlike other treatments14,15 ). Moreover, since our model is an expansion in the truly small parameter θ−1 (see Eq. (22)) it becomes exact in the large θ (long time) limit. Finally, note that in going to the continuum limit we lost the periodic dependence of F on θ. This periodicity is generated by the first step which follows a unique equation of motion. We did not incorporate this unique behavior into the continuum model and therefore should not be surprised that this information is lost. As we emphasized, our continuum model is an exact representation of the original discrete system. For this reason it is interesting to compare it with other continuum models, which do not emerge as limits of discrete systems of steps. Many authors use the continuity equation ∂h + ∇J~ = 0 , ∂t

(29)

to account for the surface evolution. Here h is the surface height and J~ is the current density of diffusing adatoms. This equation is of course correct and reduces the problem to ~ It is widely assumed12,14 that J~ is proportional to the gradient of the surface calculating J. chemical potential µ, as expected in diffusive systems. In the case of attachment-detachment limited systems it was suggested by Nozieres27 that the chemical potential gradient should be divided by the profile slope. Our model is consistent with the first picture in the DL case and with the second in the ADL case. To show this we take the gradient of Eq. (29). This leads to ~ ∂D ∝ ∇∇ · J~ , ∂t

(30)

~ is the gradient of the profile. We now return to Eq. (28) and transform it back where −D ′ to an equation for the step density D(ρ, τ ). Remembering that the term − xF4 in Eq. (28) arises from the time derivative of D, we find that in the general case we can write 11

∂ 1 ∂ ρ ∂µ ∂D + · ∂τ ∂ρ ρ ∂ρ 1 − q + 2qD ∂ρ ! 1 ∂D D2 µ= +g . + 3D ρ ρ ∂ρ

!!

= 0 with, (31)

Eq. (31) is nothing but the radial component of Eq. (30) written in terms of the dimensionless variables ρ and τ with J=

∂µ 1 · . 1 − q + 2qD ∂ρ

(32)

Therefore in the DL case (q = 0) the adatom current is indeed proportional to the chemical potential gradient. In the ADL case (q = 1) the chemical potential gradient is divided by the local step density as suggested by Nozieres. In addition we note that in the limit ρ → ∞ (where the steps are nearly straight), our chemical potential (Eq. (31)) becomes identical to the chemical potential of Ref. 12. V. PROPERTIES OF THE SCALING FUNCTION

In this section we use Eq. (28) to study properties of the scaling function F . We also discuss the boundary conditions necessary to solve Eq. (28). We begin by noting that according to the numerical simulations of the step model, there is a growing plateau or facet at the top of the cone. Our numerical solutions for the scaling function indicate that all physical solutions indeed have such a special point, which can be identified as the facet edge. Let us examine the behavior of the step density on the facet and at its edge. At any given time there is only one single step on the facet (the first step during its collapse towards the origin). We have seen that the size of the facet grows indefinitely, and therefore the step density on the facet vanishes in the long time scaling limit. Moreover, in Appendix C we show that in the long time limit, the step density is a continuous function even at the facet edge; i.e., it goes to zero continuously as the facet edge is approached from above. Denoting the scaled position of the facet edge by x0 , the above observations can be expressed in the form F (x) = 0 , ∀x ≤ x0 .

(33)

Note that F (x) = 0 is not a solution of Eq. (28). This does not contradict the continuum model, since the model was derived only for the case of a finite step density (see Eq. (22)). Thus the scaling function has to satisfy Eq. (28) only at x > x0 , and x = x0 is a singular point. Let us now study the nature of this singularity at the facet edge. Returning to our simulation data we note that near the facet edge, the scaling function is extremely steep (see Fig. 3). Indeed, the expansion of Eq.√(28) in small F suggests that in the vicinity of x0 , F (x) can be written as a power series in x − x0 . Thus although F is continuous at x0 , all its derivatives with respect to x diverge at the singular point. We now turn to discuss the boundary conditions necessary to solve Eq. (28). First, we examine the behavior of F at infinity. Far enough from the origin the steps do not move. Hence in the limit ρ → ∞ the step density remains at its initial value. This implies that lim F (x) = 1 .

x→∞

(34)

It turns out that this condition corresponds to two boundary conditions, since it is only possible to satisfy Eq. (34) if F ′ vanishes at infinity. In fact, by considering the asymptotic expansion of F it can be shown that for large x 12

F = 1 + a4 x−4 + O x−8 ,

(35)

where a4 = 3 (1 + g) in the DL case and a4 = 32 (1 + g) in the ADL case. In order to solve Eq. (28) we need two additional boundary conditions at x = x0 . The first of these is F (x0 ) = 0 (see Eq. (33)). Another boundary condition at x = x0 can be derived by enforcing volume conservation. From Eq. (15), the volume change during the time τ is given by 2π

Z

0

∞

ρ [h(ρ, τ ) − h(ρ, 0)] dρ .

(36)

Integrating this equation by parts and requiring volume conservation, we get ∞

Z

0

ρ2 [D(ρ, τ ) − D(ρ, 0)] dρ = 0 .

(37)

In terms of scaled variables the last equation can be written as Z

∞

0

x2 (F − 1) dx =

Z

∞

x0

x2 (F − 1) dx −

x30 =0, 3

(38)

where we have used the fact that F vanishes below x0 . The evaluation of the integral in Eq. (38) can be done by multiplying Eq. (28) by x2 and integrating it with respect to x. The resulting equation is Z

∞

2

x

x0

F ′A + F 2 B dx = 2 !

Z

∞

x0

x3 F ′ dx . 4

(39)

′

The integral M2 = x2 F 2A + F 2 B dx on the l.h.s. is carried out in appendix C, and is expressed in terms of F and its derivatives. The result of this integration combined with integration by parts of the r.h.s. of Eq. (39) leads to R

Z

∞

x0

i ∞ 1h 3 x (F − 1)dx = x (F − 1) − 4M2 . 3 x0

2

(40)

Inserting this relation in Eq. (38) we obtain the boundary condition M 2 | x0 = 0 .

(41)

We have used the facts that the surface term at infinity in Eq. (40) vanishes and F (x0 ) = 0. At this point we have four boundary conditions, two at infinity and two at x0 . We may now obtain a unique solution of Eq. (28) if we know the value of x0 . What determines x0 ? To answer this question, consider the height of the cone at the origin, h0 , at time τ = τn−1 and at time τ = τn . τn−1 and τn are the disappearance times of two successive steps, and therefore h0 (τn−1 ) − h0 (τn ) = 1. We can also calculate this difference from Eq. (14). Combining these two results we arrive at the relation

τn1/γ

−

1/γ τn−1

Z

0

∞

(F − 1) dx = −1 .

(42)

Recalling that τn = (θ0 n)γ and using the fact that F vanishes below x0 , we can rewrite the last equation as 13

∞

Z

x0

(F − 1) dx − x0 = −θ0−1 .

(43)

The integral in Eq. (43) can be evaluated by integrating Eq. (28) with respect to x. The result is: Z

∞

x0

R

(F − 1) dx = 4M0 |x0 + x0 ,

(44)

′

F A + F 2 B dx is carried out in Appendix C, and is expressed where the integral M0 = 2 in terms of F and its derivatives. Combining Eq. (43) and Eq. (44) we obtain the relation

4M0 |x0 = −θ0−1 .

(45)

The l.h.s. of this equation depends on x0 , thus relating x0 to θ0 , the scaled collapse time period of the steps. What is then the value of θ0 ? It is determined by the motion and collapse of the first step. The behavior of the first step is different from that of all the other steps, since it does not have neighboring steps with smaller radii. Thus, our continuum model, which treats all the steps on equal footing, does not contain any information on the value of θ0 . We therefore expect a family of valid scaling functions consistent with the continuum model with different values of x0 or θ0 . The unique value of x0 observed in simulations is determined by the discrete nature of the steps, and at this stage we are not able to calculate it. VI. NUMERICAL EVALUATION OF THE SCALING FUNCTION

We now find the scaling function numerically. Consider first the DL case. We choose a value of x0 and solve equation (28) √ starting at x = x0 . As explained in the previous section, F can be expanded in powers of x − x0 in the vicinity of x0 : F (x) =

∞ X

n=1

kn

√

x − x0

n

,

(46)

where we have already used the boundary condition (33). We therefore have three additional free parameters in this expansion, which should be determined by the boundary conditions (34) and (41). We choose these parameters as the first three odd coefficients in the expansion (46): k1 , k3 and k5 (the first two even coefficients vanish). Using the boundary condition (41), we can express k5 as a function of k1 and k3 through the relation k5 =

k32 k3 1 k1 − + − , 2 6x0 2k1 9x0 6gx30 k1

thus reducing the number of free parameters to two. Finally we use the boundary condition (34) to find the coefficients k1 and k3 . Since (34) addresses the value of F at infinity, it cannot be easily applied to the expansion of F in the vicinity of x0 . We therefore start from an initial guess for the values of k1 and k3 , solve Eq. (28) numerically for this choice of parameters and then tune k1 and k3 until the boundary condition (34) is satisfied. For a given choice of k1 and k3 , we found the solution of Eq. (28) in the following way. Numerical integration starting at x0 is impossible because the derivatives of F diverge there. Therefore, we first used the expansion (46), truncated at a high enough order, to evaluate the function F and its first three derivatives at x = x0 + δx, for some choice of δx. Then, using these values we integrated Eq. (28) numerically from x0 + δx. We made sure that the solution is not sensitive to the choice of δx. 14

a)

b)

1

1

0.8

0.8

F 0.6

F 0.6

0.4

0.4

0.2

0.2

0 0

2

4

6

8

0 1

10

2

3

−1/4

4

5

6

x=ρ τ−1/4

x=ρ τ

c) 1.2 1

F

0.8 0.6 0.4 0.2 0 1.5

2

2.5

3

3.5

4

4.5

x=ρ τ−1/4

FIG. 6. Numerical solutions for the DL scaling function compared with simulation data. Results for three different values of the step-step interaction parameter are shown: (a) g = 0.1, (b) g = 0.01 and (c) g = 0.001. When the step-step interaction is strong, the x∗0 solution (dashed line) agrees very well with the simulation data (circles). As the value of g is reduced, the best fit solution of Eq. (28) (solid) deviates from the x∗0 solution to higher values of x0 .

The above procedure was employed to generate the scaling function for a range of values of the scaled facet edge position, x0 , and the following picture emerged. There is a special minimal value of x0 , which we denote by x∗0 . For x0 < x∗0 there is no solution of Eq. (28) which satisfies the boundary conditions. For each value of x0 ≥ x∗0 there is a unique solution which satisfies the boundary conditions. Thus, as we anticipated, there is a one dimensional family of scaling functions parameterized by the value of x0 ≥ x∗0 . Since x0 is related to θ0 through Eq. (45), this family of solutions can alternatively be labeled by the value of θ0 . We used Eq. (45) to calculate θ0 (x0 ) and found that it is a monotonically decreasing function. In particular, θ0 (x0 ) is maximal at x∗0 . This is reasonable, since it means that a smaller facet corresponds to a longer period and thus to a slower facet edge. Despite the existence of many possible scaling functions, our simulations suggest that the system reaches a unique scaling solution independent of initial conditions. What is the selected solution? In Figs. 6 we compare the calculated scaling functions with simulation data in the DL case. We do this for three different values of g, the interaction strength parameter. When g is large there is an impressive agreement between the simulations data and the x∗0 solution (Fig 6 (a)). When g is reduced (i.e., for weaker interactions between steps) the observed scaling function deviates from the x∗0 solution (Fig. 6 (b) and (c)). 15

However, in these cases there is another solution with a larger value of x0 that best fits the simulation data. The agreement between this best fit solution and the simulation data is again excellent. The above observations suggest that x∗0 is the selected solution in the large g limit. We propose the following argument to support this scenario. Since the parameter g is a measure of the strength of the step-step interaction G relative to the step line tension Γ, the large g limit is equivalent to the small Γ limit. The collapse of the first step is driven by the step line tension. In the large g limit the collapse driving force is minimal and the collapse time period is maximal. As we mentioned above, a long collapse period is equivalent to a large value of θ0 , which corresponds to a small value of x0 . Thus the large g limit corresponds to the minimal value of x0 , i.e. x = x∗0 . a)

b) 1.5

1

0.8

1

F 0.6

F

0.4

0.5

0.2

0 0

2

4

6

8

0 1

10

−1/4

1.5

2

2.5

3

3.5

4

4.5

5

x=ρ τ−1/4

x=ρ τ

FIG. 7. Numerical solutions for the ADL scaling function compared with simulation data. Results for two different values of the step-step interaction parameter are shown: (a) g = 0.2 and (b) g = 0.01. When the step-step interaction is strong the simulation data (circles) agrees with the x∗0 = 1.09 ± 0.07 solution (dashed line). As g is reduced, we notice that the best fit solution of Eq. (28) (solid line) deviates from the x∗0 = 1.4 ± 0.07 solution to higher values of x0 .

√ Now consider the ADL case. Although in this limit we can also expand F in x − x0 , the numerical procedure described above is not an effective method of solution in this case. It turns out that the resulting scaling functions are sensitive to the choice of δx. We therefore had to use a different method to solve Eq. (28) in the ADL case. Let us denote by xpeak the minimal value of x, which corresponds to a local maximum of F . Such a point must exist for any continuous function satisfying the integral condition (38) and the boundary conditions (33) and (34). By tuning the values of F (xpeak ), F ′′ (xpeak ) and F ′′′ (xpeak ) we found a one dimensional family of scaling functions, which satisfy the boundary conditions. This family can be parameterized by the value of xpeak or alternatively by the value of the resulting x0 . As in the DL case, we find that there is a minimal value of x0 , denoted by x∗0 , below which there are no solutions satisfying the boundary conditions. Figs. 7 show the simulation data compared with calculated scaling functions for two values of the interaction strength g. We again find that for large g the simulation data agrees with the x∗0 solution (dashed line). When g is reduced, the agreement deteriorates and there is a different solution with x0 > x∗0 (solid line) which best fits the simulation data. VII. EFFECTS OF STEP PERMEABILITY

The step flow model we introduced in section II assumes that steps are impermeable; i.e., adatoms cannot hop between neighboring terraces without being incorporated in the 16

step separating them. In general, however, steps may be permeable. For example, in a recent paper16 , Tanaka et al. interpret their experimental results as evidence that steps on Si(001) are permeable. We therefore ask the following question: What is the effect of step permeability on decay of nanostructures in general, and on the decay of an infinite cone in particular? In this section we show that the scaling exponents of an infinite cone of permeable steps are identical to the exponents associated with impermeable steps. Moreover, the only effect of step permeability on the differential equation for the scaling function is a renormalization of one of its parameters. Following Ref. 16 we generalize our step flow model to include step permeability. We assume that flux of adatoms between two neighboring terraces due to direct hops is determined by first order kinetics. Introducing the permeability coefficient p, we rewrite Eq. (3) as

∂Ci eq Ds + p Ci |ri − Ci−1 |ri = k C i | ri − C i ∂r ri

∂Ci−1 eq − p Ci |ri − Ci−1 |ri . = k Ci−1 |ri − Ci −Ds ∂r ri

(47)

Step permeability does not affect the diffusion equation itself, and therefore the general form of the diffusion field, given by Eq. (2), remains unaltered. However the coefficients Ai and Bi are affected, and the equations for coefficients of different terraces become coupled. As a result, the step equations of motion cannot be written in closed form. The scaling analysis therefore becomes more cumbersome. As in the case of impermeable steps we assume that in the long time limit the scaling ansatz Eq. (11) holds, with the same definitions of the scaled variables, x and θ in terms of ρ and τ . We slightly change the definition of the dimensionless time τ to τ = t/t0 , where t0 is a time scale which we can choose to our convenience. The derivation of the exponents α and β in section IV did not involve the dynamics of the system. Thus the same derivation still holds and the values of α and β are not affected by step permeability. We proceed to evaluate the exponent γ and the scaling function F by studying Eq. (47). For convenience, we parametrize the solutions of the diffusion equation by concentration differences instead of the parameters Ai and Bi . These concentration differences, Ui and Vi , are defined as follows: Ui = Ci (ri ) − C¯ eq = Ai ln ri + Bi − C¯ eq ri Vi = Ci (ri ) − Ci (ri+1 ) = Ai ln . ri+1

(48)

In the continuum limit, Ui and Vi are continuous functions of r and t or alternatively of x and θ. We associate the values of Ui and Vi with the middle of the ith terrace, namely xi + xi+1 U , θ = Ui , 2 xi + xi+1 , θ = Vi . V 2

(49)

The scaling scenario for the functions U and V is

U (x, θ) = θµ u (x) + O θ−1

V (x, θ) = θν v (x) + O θ−1

, .

In terms of Ui , Vi and the dimensionless radii ρi , Eq. (47) takes the form: 17

(50)

T Ds V i ¯ eq ρi = k Ui − C ξi + p (Ui − (Ui−1 − Vi−1 )) , ΩΓρi ln ρi+1

−

T Ds Vi−1 eq ¯ = k U − V − C ξ i−1 i−1 i − p (Ui − (Ui−1 − Vi−1 )) . ΩΓρi ln ρi−1 ρi

(51)

Evaluating Eq. (51) at x = (xi + xi+1 )/2 we can now employ Eq. (24) to expand ξi , ρi and ρi ln ρi−1 in powers of θ−1 : ρi ln ρi+1 ρi 1 F2 ξi = θ−1 + O θ−2 , +g + 3F F′ x x ρi 1 ρi ln = − + O θ−1 , ρi+1 F ρi−1 1 ρi ln = − + O θ−1 . ρi F !!

(52)

The xi ’s dependence of U and V is unknown so we expand: ∞ n X xi−1 + xi 1 ∂ n u(x) xi−1 + xi = − x n 2 2 n=0 n! ∂x n ∞ n X 1 ∂ v(x) xi−1 + xi xi−1 + xi = −x . v n 2 2 n=0 n! ∂x

u

(53)

Since the difference between successive x’s is of order θ−1 , Eq. (53) is also an expansion in this small parameter. Using Eq. (53) together with Eqs. (49) and (50) we isolate v(x) in Eqs. (51), keeping only the lowest orders in θ−1 . θ1+ν v(x) = − θ1+ν v(x) =

ΩΓk θ1+µ u(x) − C¯ eq

ΩΓk θ

1+µ

1 x

+g

F2 x

+ 3F F ′

Ds T F + pΩΓ 2 u(x) − C¯ eq x1 + g Fx + 3F F ′ Ds T F + (p + k)ΩΓ

.

(54)

Since v(x) cannot vanish identically, the two expressions above are consistent only if ν < −1, µ = −1 and ¯ eq

u(x) = C

"

1 F2 +g + 3F F ′ x x

!#

.

(55)

Thus, both sides of Eqs. (54) decay in time. We now subtract the second line of Eq. (51) from the first line and obtain the following equation:

Ds T Vi−1 Vi = (k + 2p) (Ui − Ui−1 + Vi−1 ) . ρi−1 + ρi ΩΓ ρi ln ρi ρi ln ρi+1 Again we expand to lowest order in θ−1 and isolate v(x): 18

(56)

v(x) = −θ−(2+ν)

ΩΓ (k + 2p) u′ (x) . 2Ds T F 2 + ΩΓ (k + 2p) F

(57)

v(x) does not depend on θ. This implies that ν = −2 and v(x) is proportional to u′ (x). To finish the scaling analysis we return to Eq. (21) and find the leading orders in θ−1 of ρ˙ i+1 + ρ˙ i and ρ˙ i+1 − ρ˙ i . Using Eqs. (5) and (48) we find that

Ds T 2 V i Vi−1 ρ˙ i = . t0 ρi − 2 ΩΓ ρi ln ρi+1 ρi ln ρi−1 ρi

(58)

Putting together Eqs. (50) (52) and (58) we conclude that to lowest order in θ−1 ρ˙ i+1 + ρ˙ i = θ−3 (Ap + O(θ−1 )) , ρ˙ i+1 − ρ˙ i = θ−4 (Bp + O(θ−1 )) .

(59)

Ap and Bp are expressions involving F , F ′ , F ′′ , F ′′′ and F ′′′′ . Note that the orders of θ−1 in these expressions are identical to those in the equivalent expressions in section IV (Eq. (25)). These orders are responsible for the result γ = 4 in the case of impermeable steps. Therefore, in the permeable case we also have γ = 4. After expanding ρ˙ i+1 + ρ˙ i and ρ˙ i+1 − ρ˙ i we can use Eq. (21) to obtain the differential equation for the scaling function in the permeable case. However this exercise is not necessary. We avoid it by arguing that the resulting differential equation in the permeable case is equivalent to the one in the impermeable case with renormalization of some of the parameters. To see this we note that the treatment we presented here is valid also in the impermeable case. Thus when p = 0, the general differential equation must be equivalent to the equation derived in Section IV. In addition, the attachment/detachment rate k and the permeability p affect the step velocities only through the function v(x), which depends on k and p only trough the sum k + 2p. Therefore, k and p affect the differential equation itself only trough the sum k + 2p. We conclude that the scaling function in the case of permeable steps is identical to the scaling function of impermeable steps with k + 2p in the former replacing k in the latter. VIII. SUMMARY AND DISCUSSION

We have studied the relaxation of an infinite crystalline cone below the roughening temperature, in terms of a step flow model. The model was solved numerically and two types of dynamical evolutions were found. When the repulsive interactions are strong enough, the decay of the cone proceeds through the collapse of the innermost steps, one step at a time. Weak interactions lead to a step bunching instability (except in pure diffusion limited kinetics), and the decay process becomes much more complicated and involves collapse of bunches of steps. Focusing on stable cases, we found that in the long time limit the decaying step system obeys a scaling scenario. The step density (i.e. the slope of the height profile), defined as the α/γ −β/γ 1/γ . F is a inverse step separation, scales in time according to D(ρ, τ ) = τ F τ ρ, τ −β/γ function of the scaled position x = τ ρ and exhibits a periodic dependence on the scaled time θ = τ 1/γ . In particular, the position of the facet edge at the top of the cone grows as τ β/γ . The values of the scaling exponents which fit our simulations are α = 0, β = 1 and γ = 4. Following this observation we used a scaling ansatz to transform the discrete step flow model into a continuum description of surface evolution. The basic predictions of this continuum model are the values of the scaling exponents (which agree with the simulation results) 19

and a differential equation for the scaling function F . This continuum model becomes exact in the long time scaling limit, and it breaks down whenever the step density vanishes. This fact can be seen both from the derivation of the continuum model and from the resulting differential equation, which predicts a singular behavior at the zeros of the scaling function F. We showed that each physical solution of the continuum equation has a special point, x0 , at which the scaling function, F , vanishes. F is singular at x = x0 , and we identify this special point as the scaled position of a facet edge. Thus the edge of a macroscopic facet in the discrete system is a singular point of the scaling function F . A detailed analysis of the discrete system revealed a sufficient number of boundary conditions, which define a unique scaling function, F , for a given value of the scaled facet edge position, x0 . However, we were not able to find a unique value for x0 , and were left with a one dimensional family of solutions, parameterized by x0 . A numerical solution of the differential equation for F confirmed the existence of this family. We found that there is a minimal value of x0 , which we denoted by x∗0 . For every x0 ≥ x∗0 there is a unique scaling function, while for x0 < x∗0 there are no solutions, which satisfy the boundary conditions. A comparison of the numerical solutions of F with results of simulations of the discrete system leads to the following picture. When the step-step repulsion is strong, there is a remarkable agreement between the x0 = x∗0 solution and the simulation data. When the magnitude of the step-step interaction is reduced, the system reaches scaling solutions with x0 > x∗0 . We therefore advance the hypothesis that in the strong interaction limit the system approaches the minimal x∗0 solution. Although our work provides a detailed account of the decay process of an infinite cone, it leaves a few unresolved issues. First and most important is the fact that we were not able to uniquely determine the value of x0 . We did show that x0 depends crucially on the detailed collapse process of the innermost step, which is not included in the continuum model. Thus, in order to evaluate x0 , one has to deal directly with the discrete step system. Another open issue related to the behavior of the innermost step is the periodic behavior of the scaling function. This periodicity is observed in the kinetics of the discrete system of steps, and is absent from the solution of the continuum equation. Lastly, when the repulsive interaction between steps is weaker than a certain threshold, the system becomes unstable and step bunches are formed. The critical value of the interaction below which this instability occurs and its dependence on the kinetic parameter q have not been studied so far. We intend to address these open questions in future work. Finally, we remark that the scaling behavior of the cone profile, predicted in this work, is robust. In particular, the existence of a facet growing as τ 1/4 , as well the existence of a scaling state does not depend on the detailed form of the repulsive interactions between steps. A quantitative change of these interactions may alter the scaling function, F , but the scaling exponents do not change. Another manifestation of the robustness of the scaling solution is the effect of step permeability. In principle, step permeability could have changed the scaling behavior of the system entirely. However, we showed that its only effect on the scaling solution is to modify the attachment-detachment rate coefficient, k, to k + 2p. We are greatful to N. Bartelt, H. C. Jeong, D. J. Liu and J. D. Weeks for helpful discussions. This research was supported by grant No. 95-00268 from the United StatesIsrael Binational Science Foundation (BSF), Jerusalem, Israel. D. Kandel is the incumbent of the Ruth Epstein Recu Career Development Chair. APPENDIX A:

Here we study the linear stability of circular, uniformly spaced steps with unit step separation in the absence of step-step interactions. This configuration with ρn = n is not a steady state. However, in the large radius limit the steps move very slowly, and the uniform state is extremely close to a steady state. We therefore refer to it as a quasi steady state. We regard the quasi steady state as unstable if the growth of perturbations is faster than the motion of steps. 20

In the absence of step-step interactions Eqs. (10) simplify to ρ˙ n ≡ an =

dρn an − an−1 , with = dτ ρn 1 1 − ρn+1 ρn

ρn −q (1 − q) ln ρn+1

1 ρn

+

(A1)

1 ρn+1

(A2)

For large values of n, the velocity of the nth step in the quasi steady state is given by: ρ˙ n =

n−3 + O n−4 . 1+q

(A3)

To check the linear stability of the above configuration we perturb the step positions according to ρn = n + ∆ei(φn−ωt) .

(A4)

Equating the time derivative of this perturbation with the velocity of the nth step results in an equation for ω. We find that to the lowest orders in ∆ and n−1 ω=

4iq (1 − cos φ) n−2 . (1 + q)2

(A5)

Since the magnitude of ω is significantly larger than the step velocities in the quasi steady state, positive values of Im(ω) lead to an instability. We see that away from the DL case (q > 0) the system is unstable. The most unstable mode is φ = −π (step pairing). The φ = 0 mode (uniform translation) is marginal as is the q = 0 (DL) case. APPENDIX B:

In this appendix we give some technical details of the algebraic manipulations performed in sections IV and V. To calculate expressions A and B in Eq. (25), we first express the scaled positions of the steps as power series in θ−1 . Inserting Eq. (24) into Eq. (23) we find that to fifth order in θ−1 −1

′

−1

′

−2

24 F

24 F 9

−1

xi

9

−3 F ′2 + F F ′′ θ−3

= x−

θ 2F

21

−93 F ′3 + 61 F F ′ F ′′ − 6 F 2 F (3) θ−4

33 F ′3 − 21 F F ′ F ′′ + 2 F 2 F (3) θ−4

F θ 3θ − + − F3 2F5 12 F7 2F 4 2 2 −135 F ′ + 126 F F ′ F ′′ − 12 F 2 F ′′ − 16 F 2 F ′ F (3) + F 3 F (4) θ−5

xi−1 = x − +

5θ 3F θ − + + F3 2F5 4 F 7 2F 2 4 2 1935 F ′ − 1890 F F ′ F ′′ + 260 F 2 F ′ F (3) + F 2 180 F ′′ − 17 F F (4) θ−5

xi−2 = x − −

5 −3 F ′ 2 + F F ′′ θ−3

−2

xi+1 = x +

θ−1 2F

−3 F ′2 + F F ′′ θ−3 33 F ′3 − 21 F F ′ F ′′ + 2 F 2 F (3) θ−4 3 θ−1 F ′ θ−2 xi+2 = x + − − − F3 2F5 12 F7 2F ′4 ′ 2 ′′ 2 ′′ 2 2 ′ (3) 3 (4) θ−5 −135 F + 126 F F F − 12 F F − 16 F F F + F F − 24 F 9 ′3 ′ ′′ 2 (3) −3 ′2 ′′ −1 ′ −2 θ−4 −93 F + 61 F F F − 6 F F θ 5 −3 F + F F 5θ 3F θ xi+3 = x + − − + F3 2F5 4 F 7 2F 2 4 2 1935 F ′ − 1890 F F ′ F ′′ + 260 F 2 F ′ F (3) + F 2 180 F ′′ − 17 F F (4) θ−5 , (B1) + 24 F 9 where F and its derivatives are evaluated at x = (xi + xi+1 ) /2. Using these expressions we expand the step velocities in θ−1 and obtain Eq. (25). In the general case the expressions for A and B are too cumbersome to be written here. Instead we give these expressions in the two limiting cases. In the DL case ADL = g

2F 2 4 F ′ 10 F ′2 10 F ′′ 18 F ′ F ′′ − 2 + + + + 6 F (3) + 3 x x Fx x F F x3

BDL = g

−

!

3 5F′ 5 F ′3 7 F ′′ 10 F ′ F ′′ 9 F ′′2 F ′ 2 (5 F + 9 x2 F ′′ ) + − − + + − x4 F x3 F3 x F x2 F2 x F2 F 3 x2

5 F (3) 9 F ′ F (3) 3 F (4) + + + Fx F2 F

!

−

3 F′ − . F 2 x4 F 3 x3

(B2)

Inserting Eq. (B2) into Eq. (28) we find the differential equation which govern the scaling function in the DL case: 2

g 12 F ′ F (3) + 3 F F (4) + 9 F ′′ + −

7 F ′ 2 + F F ′′ x2

In the ADL case

15 F ′ F ′′ + 5 F F (3) x

6 F F′ 3 F2 xF′ 3 + − − =0. − x3 x4 x4 4

AADL = g

F′ 3 F ′2 3 F ′3 5 F ′′ 6 F ′ F ′′ 3 F (3) 1 − + 2 − + + + x3 F x2 F x F3 Fx F2 F

BADL = g

−

!

(B3)

1 F′ + 2 3+ 3 2 , F x F x

F′ F ′2 3 F ′3 9 F ′ 4 3 F ′′ F ′ F ′′ 21 F ′2 F ′′ 3 + − − + − + − 2 F x4 F 2 x3 F 3 x2 F4 x 2F5 F 2 x2 2 F 3 x 2F4

3 F ′′ 2 5 F (3) 3 F ′ F (3) 3 F (4) + 3 + + + F 2F2 x 2F3 2F2

!

−

3 2F′ 3 F ′2 F ′′ − − + , (B4) 2 F 3 x4 F 4 x3 2 F 5 x2 2 F 4 x2

and the differential equation for the ADL scaling function is: 3F 3F′ 3 F ′2 3 F ′3 3 F ′ 4 3 F ′′ 3 F ′ F ′′ 15 F ′2 F ′′ 3 F ′′2 5 F (3) g − 4+ − − + − 2 + − + + 2x 2 x3 2 F x2 2 F 2 x F3 x Fx 2F2 F 2x 3 F ′ F (3) 3 F (4) + + F 2

!

−

3 3F′ xF′ F ′2 F ′′ − − − + =0. 2 F x4 2 F 2 x3 4 F 3 x2 2 F 2 x2 22

(B5)

In section V we used the moments F ′A M0 = + F 2 B dx , 2 ! Z ′ 2 F A 2 M2 = x + F B dx . 2 !

Z

to set boundary conditions for the scaling function at x0 . In the DL case M0DL M2DL

F 2 2F F ′ 5F ′2 5F F ′′ 1 =g − + + + 9F ′ F ′′ + 3F F ′′′ + 3 , 3 2 x x x x x ! 2 3 3F ′ ′2 ′′ 2 ′ ′′ 2 ′′′ + . − 6F F − xF − xF F + 9x F F + 3x F F =g x x !

In the ADL case M0ADL

M2ADL

F F′ 3F ′2 3F ′3 5F ′′ 3F ′ F ′′ 3F ′′′ =g − + − + + + 2x3 2x2 2xF 2F 2 2x F 2 ′ F 1 + 3 + 2 2 , 2x F 2x F ! 5F ′ 3xF ′2 3x2 F ′3 xF ′′ 3x2 F ′ F ′′ 3x2 F ′′′ 3F − − − − + + =g 2x 2 2F 2F 2 2 F 2 ′ 3 F + + . 2xF 2F 2 !

APPENDIX C:

In this appendix we show that in the scaling state, the step density near the facet edge must vanish when the facet size diverges. This implies that F (x0 ) = 0 where x0 is the scaled position of the facet edge. We restrict ourselves to situations where steps collapse towards the origin one at a time (consistently with simulation results). It is tempting to argue that since our system is expanding and slowing down, in the long time limit it approaches the equilibrium state of straight steps in contact with a facet. In this equilibrium state the step density near the facet edge vanishes as a square root. This equilibrium state is not consistent with our density function, which also predicts a square root approach to zero, but with a time dependent coefficient. Moreover even in the long time limit steps in our system continue to collapse. In each collapse period, just before the step disappears, its chemical potential diverges. At these times the system is not close to equilibrium. Consider the velocity of the first step. It depends on the positions of the first three steps ρ1 , ρ2 and ρ3 through ρ˙ 1 =

1 ξ1 − ξ2 · ρ1 (1 − q) ln ρ1 − q 1 + ρ2 ρ1

1 ρ2

.

(C1)

ρ2 is always larger then ρ1 so the denominator is always negative. The direction of motion of the first step is thus given by the sign of the numerator 23

1 1 1 1 ρ3 − + 2g N = ξ1 − ξ2 = 3 − ρ1 ρ2 ρ2 + ρ3 (ρ3 − ρ2 )3 (ρ2 − ρ1 )

!

,

(C2)

which must be positive to avoid bounding of the first and second steps. Fixing ρ2 and ρ3 , the value of ρ∗1 which minimizes N is found by solving

∂N 6g 1 = − ∗2 = 0 . ∗ 4 ∂ρ1 ρ∗ (ρ2 − ρ1 ) ρ1

(C3)

1

The only solution between ρ2 and the origin is ρ∗1

= ρ2 +

s

3g − 2

s

q

ρ2 6g +

3g . 2

This is indeed a minimum since the second derivative of N with respect to ρ1 is always positive when 0 ≤ ρ1 ≤ ρ2 . Substituting ρ∗1 into N we find that N (ρ∗1 ) =

1 ρ2 +

q

3g 2

−

q

√

ρ2 6g +

3g 2

−

2gρ3 1 2g + q √ . q 3 − ρ2 3g (ρ2 + ρ3 )(ρ3 − ρ2 )3 − ρ2 6g + 3g 2 2

(C4)

The first three terms of N (ρ∗1 ) decay when ρ2 is large. The fourth term however will remain finite unless the difference ρ3 − ρ2 diverges with ρ2 . If the fourth term does remain finite, N (ρ∗1 ) becomes negative, the velocity of the first step becomes positive and the first step cannot pass ρ∗1 on its way to the origin. The above analysis indicates that the step density vanishes near the facet edge. We consider two scenarios for the collapse of the first step. In the first scenario the step starts the collapse from a position below ρ∗1 . In this case the √ separation between the first and second steps is initially divergent since ρ2 − ρ∗1 diverges as ρ2 . But the initial configuration in the collapse of one step is the final configuration of the former collapse period. Thus the former period ended with a divergent distance ρ3 − ρ2 . In the second scenario the first step starts to collapse from a position above ρ∗1 . If the first step collapses alone (i.e ρ2 is large throughout the collapse period), the distance ρ3 − ρ2 must grow to allow ρ1 to pass through ρ∗1 . Either way there must be some point in the collapse period where the separation ρ3 − ρ2 diverges. At this point the step density near the facet edge vanishes.

∗ ∗∗ 1

E-mail: [email protected] E-mail: [email protected] W. K. Burton, N. Cabrera and F. C. Frank, Philos. Trans. R. Soc. London, Ser. A 243, 299 (1951). 2 N. Israeli and D. Kandel, Phys. Rev. Lett. 80 3300 (1998). 3 S. Tanaka, C. C. Umbach, J. M. Blakely, R. M. Tromp and M. Mankos, J. Vac. Sci. Technol. A 15(3) 1345 (1997); Mat. Res. Soc. Symp. Proc. 440 25 (1997). 4 J. Blakely, C. Umbach and S. Tanaka, Dynamics of Crystal Surfaces and Interfaces 23 (1997). 5 C. C. Umbach, M. E. Keeffe and J. M. Blakely, J. Vac. Sci. Technol. A 9, 1014 (1991).

24

6

K. Yamashita, H. P. Bonzel and H. Ibach, Appl. Phys. 25, 231 (1981); E. S. Fu, M. D. Johnson, D. J. Liu, J. D. Weeks and E. D. Williams, Phys. Rev. Lett. 77, 1091 (1996). 8 K. Morgenstern, G. Rosenfeld and G. Comsa, Phys. Rev. Lett. 76, 2113 (1996). 9 W. W. Mullins, J. Appl. Phys. 28, 333 (1957); 30, 77 (1959). 10 H. P. Bonzel and W. W. Mullins, Surf. Sci. 350, 285 (1996). 11 H. P. Bonzel and E. Preuss, Surf. Sci. 336, 209 (1995); Appl. Phys. A 35 1 (1984). 12 J. Hager and H. Spohn, Surf. Sci. 324, 365 (1995). 13 F. Lancon and J. Villain, Phys. Rev. Lett. 64, 293 (1990). 14 M. Ozdemir and A. Zangwill, Phys. Rev. B 42, 5013 (1990). 15 A. Rettori and J. Villain, J. Phys. France 49, 257 (1988). 16 S. Tanaka, N. C. Bartelt, C. C. Umbach, R. M. Tromp and J. M. Blakely, Phys. Rev. Lett. 78, 3342 (1997). 17 C. Duport, A. Chame, W. W. Mullins and J. Villain, J. Phys. I France 6, 1095 (1996). 18 E. Adam, A. Chame, F. Lancon and J. Villain, J. Phys. I France 7, 1455 (1997). 19 M. Uwaha, J. Phys. Soc. Jpn. 57, 1681 (1988). 20 J. Villain, Europhys. Lett. 2, 531 (1986). 21 M. V. Ramana Murty and B. H. Cooper, Phys. Rev. B 54, 10377 (1996). These Monte Carlo simulations agree with15,14 . 22 M. A. Dubson and G. Jeffers, Phys. Rev. B 49, 8347 (1994). 23 Z. Jiang and C. Ebner, Phys. Rev. B 40, 316 (1989). 24 J. D. Erlebacher and M. J. Aziz, Surf. Sci. 374, 427 (1997). These Monte Carlo simulations agree with15 and disagree with14 . 25 W. Selke and P. M. Duxbury, Phys. Rev. B 52, 17468 (1995). 26 G. S. Bales and A. Zangwill, Phys. Rev. B 41, 5500 (1990). 27 P. Nozieres, J. Phys. I France 48, 1605 (1987). 28 V. I. Marchenko and A. Ya. Parshin, Sov. Phys. JETP 52(1), 129 (1980). 29 A. F. Andreev and Yu. A. Kosevich, Sov. Phys. JETP 54(4), 761 (1982). 7

25