Spiral dislocation dynamics in crystal evaporation

SURFACE

SCIENCE 41 (1974) 77-101 0 North-Holland

SPIRAL

DISLOCATION CRYSTAL

Publishing Co.

DYNAMICS

IN

EVAPORATION T. SUREK*

Department McGill

of Mining and Metallurgical Engineering, University, Montreal, Quebec, Canada

G. M. POUND Department of Materials Science and Engineering, Stanford University, Stanford, California, U.S.A.

J. P. HIRTH Metallurgical

Engineering

Department,

Ohio State

University,

Columbus,

Ohio, U.S.A.

Received 12 June 1973 A model of concentric, circular disc-shaped holes is used to describe the ledge-wise evaporation of low index crystal surfaces. Within the assumptions of the model, an exact solution is obtained to the problem of repeated nucleation of such discs on preferred surface sites. The model also provides an approximate solution to the transient development of the spiral shape, a problem heretofore unsolved. The results are compared with earlier theories of spiral evaporation. Required conditions are proposed for observing low values of the evaporation coefficient.

1. Introduction A most significant breakthrough in the modern theory of crystal growth and evaporation has been the work of Frank’), who identified the role played by dislocations. Accordingly, a screw dislocation emergent at a point on the low index surface provides a ledge on the surface with a height equal to the component of the Burgers vector normal to the surface. Crystal evaporation of vicinal surfaces then proceeds by the terrace-ledge-kink (TLK) mechanismz-J), whereby atoms dissociate from kink positions in ledges, diffuse along the low index terraces and finally desorb into the vapor phase. At high temperatures** (T/r,,, 20.5 for molecular or covalent crystals, and * Present address: Gordon McKay Laboratory, Harvard University, Cambridge, Massachusetts 02138, U.S.A. ** The presence of straight ledges at lower temperatures implies that nucleation of a double kink is slow compared to lateral kink motion, hence the former should be rate controlling. The present analysis is concerned only with the high temperature region of behavior. 77

78

T. SUREK,

0. M. POUND

AND J. P. HIRTH

T/T,X0.05for monatomic metals, where r,, is the melting temperature), the TLK mechanism leads to a smoothly curved, spiral-shaped Iedge front which rotates at a constant angular velocity if steady state is reached. The shape of the curved spiral is determined by kinetic considerations for the movement of the curved ledge fronts and, in part, by capillarity considerations for the limiting radius of curvature at the spiral center. The model completely describes the evaporation kinetics when local equilibrium of adatoms is maintained adjacent to the ledges and where all ledges are associated with a spiral source, i.e., where ledge nucleation is excluded. In cases where there is some constraint to dissociation at the ledges?, the present treatment still describes the ledge dynamics and gives the surface diffusion contribution to the evaporation coefficiente). Numerous experimental observations5-7) of growth and evaporation spirals provide qualitative and, in some cases, direct quantitative support 8, Q) for the spiral (and the TLK) mechanism. The complete quantitative understanding of crystal growth (and evaporation) has been limited, however, by the mathematical difficulties which are encountered in the theoretical analyses. Two of the problems which have thus far escaped analytical treatment are the transient development of the spiral shape, and the exact mathematical description of the steady-state spiral shape under a range of experimental conditions. The consideration of these problems constituted the main thrust of the present work. The earlier spiral treatments were centered on the prediction of the steadystate growth rates and spiral morphologies. The most quoted result relating to the spiral mechanism is the solution by Cabrera and Levineis). They considered the case where the ledge motions are governed by the surface diffusion of adatoms and where the overlap of the adatom diffusion fields from adjacent spiral turns is negligible. Under these conditions, and with the assumption that the strain energy associated with the spiral dislocation is negligible, the steady-state spacing between successive spiral turns is given by% it,, =

19Y ______ 19'*=- (kT/V)In(p/p,)’

(1)

and the steady-state angular velocity is wcL = 21rv,/l,,

= 0.33 0,/r*.

(2)

t The transfer of atoms from kink to adsorbed state can be affected, for example, by entropy constraints caused by restricted rotational degrees of freedom of the molecules in the condensed phase, or by the dissociation of polymeric species during various stages of the stepwise kinetics. t The subscript CL is used to denote the variables in the Cabrera-Levine analysis.

SPIRAL

DISLOCATION

Here y is the specific solid-vapor

DYNAMICS

IN CRYSTAL

interfacial

EVAPORATION

free energy

surface, V is the atomic volume, T is the temperature constant. The quantity p/p, < 1 is the undersaturation p/p,> 1 and the sign in eq. (1) is reversed]; therefore,

79

of the spiral ledge

and k is Boltzmann’s ratio [for growth, r* is the critical disc

size which is in unstable equilibrium with the undersaturated vapor phase. The equilibrium vapor pressure is pe, the actual pressure is p, and the maximum velocity of a straight monatomic ledge is a,. The steady-state growth or evaporation rate is given by J=wc,d/2rr, where d is the density of atoms in the close-packed surface. Eqs. (1) and (2) are likely to be satisfied near equilibrium (i.e., p/p, N 1), and lead to Jcc (1 -p/pe)‘, which is the parabolic evaporation (growth) law at low undersaturations (supersaturations). Most evaporation experiments of interest, however, are Langmuirrr) or free evaporation experiments where the gross rate of vaporization is measured in high vacuum (p/p,=O). In this case, lzcL+ 0, and therefore the assumptions leading to eqs. (1) and (2) are no longer satisfied. The effect of overlapping diffusion fields from adjacent spiral turns at large undersaturations was considered by Burton, Cabrera and Franks) and by Cabrera and Colemanra). These solutions do not deal with the transient dynamics of the ledges, however, and they are only rough approximations to the steady-state behavior. Hirth and Pound13) replaced the spiral by concentric circular ledge fronts in order to evaluate the effect of overlapping diffusion fields. However, they treated only an asymptotic steady-state behaviour which perforce excludes the important transient dynamics. The basic problem in the TLK model of evaporation lies in the timedependent dynamics of ledge nucleation and subsequent ledge motions and interactions. A number of researchers 14-1s) have recognized that this problem is amenable to simulation on the digital computer. Recently, we described a computer simulation treatment of crystal evaporation; that treatment considered the kinetics of straight, parallel ledges which emanate from crystal edge sourcesrs). In the present paper, we describe an analogous treatment which is intended to represent spiral dislocation sources of ledges. Although this analysis is concerned mainly with crystal evaporation or dissolution, with a view to comparing the results with the predictions of earlier spiral treatments lop r2pis), many of the quantitative results are directly applicable to the reciprocal case of crystal growth from the vapor or dilute solution by the spiral mechanism, provided that the assumptions enumerated herein are satisfied. The results of the present calculations turn out to qualitatively resemble those of earlier theoretical models of spiral evaporation, but differ quantitatively. In particular, it is shown that the Hirth-Pound asymptotic steady-state condition13) is only a very special case of the dynamic spiral system. Also,

80

the concept

T. SUREK,

of a steady-state

G. M. POUND

AND .I. P. HfRTH

spiral shape is shown to be meaningful

even at

very large undersaturations. The present results can be used to delineate the conditions for observing low values of the evaporation coefficient when the evaporation kinetics is controlled by spiral dislocation sources of ledges.

2. Model of concentric circular discs A schematic top view fig. 1. An instant later because of the diffusion being the rotation of the

of a portion of an evaporation spiral is shown in in time, all the ledge fronts will have advanced of adatoms away from the ledges, the net result spiral as a whole (shape-preserving at steady state).

Fig. 1. Schematic top view of a ledge spiral associated with a screw dislocation-surface intersection in crystal evaporation.

The problem of the transient development of the spiral shape involves the solution of a time-dependent diffusion equation with moving boundary conditionsr7). It is beheved that the latter problem cannot be solved by analytical methods, whereas a numerical simulation of the transient shape of the spiral, analogous to the computer modeling of the dendrite shapel*), is deemed to be a formidabfe, but not necessarily impossible, task. Hence an approximate model is adopted here in order to yield a tractable solution. At any instant of time, some cross-section of the spiral (say at OA in fig. 1) will resemble the train of ledges emanating from a crystal edge sourceis).

SPIRAL

DISLOCATION

DYNAMICS

IN CRYSTAL

EVAPORATION

81

Two approximations are made: (i) the radius of curvature of the spiral front at a point B along OA is assumed to be equal to the radial coordinate ra of the point; and (ii) the centers of curvature of all points B along OA are assumed to coincide. Both of these approximations are very good, as can be shown by replacing the spiral with successive Archimedian segments of one turn each. In this case, the radius of curvature at r2, for example, is found to be 0.99r,; clearly, further out from the spiral center, both approximations get better. DIRECTION

OF

LEDGE

MOTION

-R Sink iR=&‘Z)

Source (R=Ol

Fig. 2. Cross section of crystal surface evaporating by the motion of concentric, discshaped holes. Source is a spiral dislocation emergent at the surface, or a preferred site for disc nucleation. There is an equivalent source a distance I away.

In essence, in the above approximations, the spirdl is replaced by concentric, disc-shaped holes (see fig. 2). By analogy to the straight, parallel ledgesis), the kinetic model in fig. 2 is used in the following to approximate the transient to steady state dynamics of ledges emanating from spiral dislocation sources. Of course, the following solution applies exactly to the problem of evaporation proceeding by the successive nucleation of concentric, discshaped holes at a preferred surface site, e.g., a dislocation-surface intersection. The consequences of this model with regard to the adatom concentration profiles (and hence the ledge velocities) and the problem of disc “nucleation” at the center of the spiral are the subjects of the next two sections of this paper. 3. Surface diffusion problem The first step in the solution of the surface diffusion problem is to compute the adatom concentration profiles in a direction normal to the spiral ledge fronts. By equating the divergence of the surface diffusion flux to the net

82

T.

SUREK,

G. M. POUND

AND J.P.

HIRTH

vaporization Aux, we find that the adatom concentration between two concentric discs of radii Ri and Ri_l (see fig. 2) is given by 1s) ?Z(R) = Aila (QR) + BiKe (QR) + no.

(3)

Here IO and KO are the modified Bessel functions of zero order, no is the adatom concentration in equilibrium with the vapor phase, and the quantity Q is given by Q=2/Z, where x is the mean free path for two-dimensional random walk of an adatom on the low index surface. The constants A i and Bj are considered further below. In the derivation of eq. (3), the following assumptions are made: (i) the effect of the strain energy associated with the spiral dislocation-surface intersection is negligible; (ii) the angular or e-dependence of the adatom concentration profile can be neglected; (iii) the concentration of adatoms at the ledge fronts is given by the Gibbs-Thomson condition n (P) = n, exp (- yVl&T) f

(4)

where p is the radius of curvature of the ledge front, and n, is the concentration of adatoms in equilibrium with the crystal; and (iv) all time dependences of the adatom concentration, including the effects of moving ledge boundaries, are neglected. For the evaporation of metal crystals, assumption fi) is usually valid. However, the effect of the strain energy can be included, by proper parametrization, into the solution of the ledge nucleation problemls), where its greatest influence is expected. Assumption (ii) merely reflects the concentric, circular disc model described earlier. Inherent in assumption (iii) is the condition that the influence of the crystallographic orientation of the ledges can be neglected. Assumption (iv) can be justified by considerations similar to the case of straight, parallel ledges sa,al). In the case of an actual steady-state spiral, the effect of the moving ledge fronts on the adatom concentration turns out to be negligible 17) when (R2w/D,)(dn/83)+ 1, where w is the angular velocity of the spiral and D, is the surface self-diffusion coefficient: to the extent that assumption (ii) holds, the same condition applies to the discs. This condition is invariably valid for real systems. Eq. (4) provides the boundary conditions for the evaluation of the constants in eq. (3); thus we find

and

Ei 10(QRi-1)- Et- 110 (QRi) Ko (QRl)'~o(QRi-I)-Ko (QRi-1) 10(QRi),'

Bi=

Vb)

SPIRAL

DISLOCATION

DYNAMICS

IN CRYSTAL

EVAPORATION

83

where

YVIRikT) - P/Al-

4 = % Cexp (-

(6)

The velocity of the ledge front at Ri has an “inward” and an “outward” component due to the diffusive fluxes of the adatoms onto the terraces Ai+l and &, respectively (see fig. 2). From the mass conservation condition for the diffusing adatoms at Ri, we have vi = dR,ldt

= vi, in + vi, OUt,

(7)

where and vi, out =

-

(Dsa~Q/ny)

CA

1,

(Q&I

-

Bi Kl

(QRi)I .

(9)

Here I1 and Ki are the modified Bessel functions of order one, a, is the width of one row of atoms and n,, is the number of atoms per unit length of ledge. The constants Ai+ 1 and Bi+ 1 are given by the general formalism in eq. (5). Eq. (8) applies to the ledge fronts in fig. 2 with L< i< (IV- I), while eq. (9) applies to all (L + 1)
Therefore

(R)

=

ne[exp

(QR) + y1 0. Io(QR,)

(-- yI’/R&T)

- P/P.] 1,

(10)

eqs. (7) and (8) can still be used if one sets A N+l=ENIIO

(Q&jand&+l=O.

(11)

The outward component of the velocity of the disc at R, can be calculated by assuming that the adatom concentration profile is symmetric about R=1/2 in fig. 2. In this case, the formalism of eqs. (7) and (9) turns out to apply with ELK, (QW) AL=K,(QWIo(QR,) +MQWKo(Q&) and (12) BL

= K,

(Q&')Io(QRL)

+ 11 (QWKo(QR,)'

where EL is given by eq. (6). For large values of QRi (say QRiS80), it is more convenient to use, instead of eqs. (8) and (9), approximate velocity expressions l?) which can be readily derived from the asymptotic expansions for the modified Bessel functions.

84

T, SUREK,

G. M. POUND

AND J. P. HIRTH

4. Disc nucleation at the source The final requisite to a solution for the spiral problem is a description of the source of ledges in this model. By analogy to the straight ledge sourcesl6), a “disc creation criterion” is sought, i.e., some critical value for the radius of the Nth disc (RN in fig. 2> before the (N-f- 1)th disc, with some critical radius R N+l, can be nucleated. In terms of the spiral parameters in fig. 1, we seek a relationship between the radius of curvature at the tip of the spiral (r,) and the variable that describes the proximity of the first spiral turn, the radial coordinate r2. The usual assignment of ri =r* and r2 z 19r* [cf. eq. (l)] is clearly not valid at large undersaturations, where r * -+ 0. Here the effect of overlapping diffusion fields from adjacent spiral turns is to cause a reduction in the local driving force or undersaturation, and thus to result in the “back-force effect”ls) at the center of the spiral. This leads to the result that the limiting radius of curvature at the spiral center must be greater than the critical disc size r* given by eq. (1). For the model of concentric, circular discs, the analysis of the adatom diffusion profiles leads to a functional relationship r1 =f (r2, p/p,, y VQ/kT) between the size of the disc-shaped hole (rl) which is in unstable equilibrium within the larger disc (r2) which last nucleated at the source. The quantities p/p, and y VQ/kT are system or material constants. The exact relationship is given byi9) [cf. eq. (I)]

where n(rl) is given by eq. (10) with R=r, and R,=r,. Eq. (13) constitutes a necessary condition between rl and r,; a second relationship is obtained by invoking a geometrical constraint between r, and r2 for dislocation spiral sources, or a nucleation rate criterion for the true problem of repeated nucleation of concentric discslg). If, for example, the spiral in fig. 1 can be approximated by an Archimedian spiral from 0 to r2, the second reIationship turns out to be Pi = r,/4n.

(14)

Fig. 3 schematically illustrates the simultaneous solutions of eqs. (13) and (14), denoted by r: and rl in the figure. The effects of varying the ratio rz/rl from 4n are readily examined in this model, as described in more detail below. The above considerations lead to the following “disc creation criterion” for dislocation spiral sources in the present model: when the preceding disc has grown to a size r:, a new disc of radius r: is nucleated [in fig. 2, when

SPIRAL DISLOCATION

DYNAMICS

IN CRYSTAL

85

EVAPORATION

12 Fig. 3.

Schematic illustration of the evaluation of the quantities rl* and ra* used in the “disc creation criterion” for the dynamic system in fig. 2.

TABLE

Summary of computer

YVQlkT

PI&Z

1

experiments for typical values of yVQ/kT, (all quantities are in dimensionless units)

Qr’

Qn *

pipe and ra*/rl*

=4x

QrP

Ii*

Azx,ss @p.)

t BS (exp.)

ass (exp.)

- 0.6 - 0.6 ~0.8 1.1 2.9

0.96 0.96 0.94 0.906 0.620

0.001 0.001 0.001 0.001 0.001

0.001 0.1 0.5 0.9 0.99

0.000145 0.000434 0.00144 0.00949 0.0995

0.02872 0.02976 0.03641 0.06419 0.16549

0.36088 0.37395 0.45753 0.80662 2.07957

0.33216 0.34419 0.42112 0.74243 1.91408

1.019 1.133 2.066 11.03 162.4

0.01 0.01 0.01 0.01 0.01

0.001 0.1 0.5 0.9 0.99

0.00145 0.00434 0.0144 0.0949 0.995

0.06247 0.06495 0.08125 0.16317 0.99503

0.78504 0.81617 1.02104 2.05043 12.5039

0.72257 0.75122 0.93979 1.88726 11.5089

1.101 1.229 2.330 16.13 1406.5

1.10 1.2 1.42 2.85 28.0

0.91 0.90 0.86 0.624 0.071

0.1 0.1 0.1 0.1

0.001 0.1 0.5 0.9

0.0145 0.0434 0.144 0.949

0.14017 0.14850 0.21257 0.94919

1.76148 1.86616 2.67125 11.9278

1.62131 1.71766 2.45868 10.9786

1.514 1.738 4.129 135.7

2.58 2.72 3.91 27.0

0.661 0.640 0.488 0.074

RN=r:, the (N+l)th disc forms with RN+$ =r:]. The parameters r: and rt are, of course, functions of the system constants p/p, and y VQ/kT. Values of rf, r; and r * (in units of I/Q) are given in table 1 for r:/r: =4d. The relative values of r: and r* are particularly noteworthy, since the ratio r:/r * > 1 indicates the extent of diffusion overlap from the first spiral turn at the spiral center, i.e., the back-force effect described above. t Complete solutions of eq. (I 3) for a range of p/pe and y VQ/kT values are given in ref. 19.

86

T. SUREK,

G. M. POUND

AND J. P. HIRTH

5. Computer simulation model of spiral evaporation In this section, we present the explicit dynamic model for dislocation spiral sources on the basis of the preceding theory. The problem of evaporation proceeding by the motion of concentric, disc-shaped holes is depicted in fig. 2. As in the mode1 for straight ledge sources16), the subscript L, denotes the leading disc in the system, while N denotes the last disc nucleated at the source. At any instant of time, therefore, there are (N-L + 1) discs on the surface (between source and sink), and (L- 1) discs have already been annihilated at the sink by discs from an equivalent source a distance 1 away (initially, L=N= 1). As in the case of the straight ledge sources, it is convenient to introduce dimensionless time and distance variables, as follows : QR, is replaced by Ri,

Q,$ is replaced by ii,

Kt is replaced by t ,

where K=

RaOQ2ne= a0 - -- pe ?v

a, (2nmkT)*’

(15)

Here Iziis the width of the ith terrace spacing in fig. 2 and m is the molecular mass. The dimensionless source-to-sink distance is Q1/2 and the critical discsize is Qr*, while the dimensionless nucleation parameters are denoted by Qr: and Qr,*. The dynamics of the discs in fig. 2 are described by a set of coupled ordinary differential equations [cf. eq. (7)] for the dimensionless velocities vi = dR,/dt for L < i
SPIRAL

DISLOCATION

DYNAMICS

IN CRYSTAL

87

EVAPORATION

vaporization are monitored during the simulation experiments as functions of the dimensionless time variable. The net vapor flux from the surface is given by J = 01”(& - P)/(2xmW

7

(16)

where c(, is the evaporation coefficient. For the kinetic model in fig. 2, we can replace CC,by an average evaporation coefficient

= T ai (areah[T

(area).]

_

I,

where (area)! is the area of the ith annular region in fig. 2, i.e., (area),= = rr (RF- 1- Rf), and Ci(area)i = n (@/a)‘. The quantity C(~can be obtained [using eq. (3)] from the average adatom concentration over the ith terrace in fig. 2; to a very good approximation, Eli =

(2/Ai) tanh (&/2) {[n(Ri) + n (Ri- 1)]/2ne} 7

(17)

where the term in the braces represents the lowering in the adatom concentration at the ledges because of the curvature effect [cf. eq. (4)J. 6. Results of computer simulation experiments Table 1 summarizes the results of computer experiments for typical values of the system constants pip, and y VQ/kT, and r,*/r: = 4~. The parameters selected were such as to give a wide range of values for the initial ledge spacing near the source, li, = Qr: - Qrf . The values of Q1/2 were between 10 (for small ,I,,) and 150 (for large 12J, but subject to certain conditions, as stated below, the value of QZ/2 had no effect on the eventual steady-state dynamics of the ledges. The annihilation parameter was il, = 0.0033 for all the simulation runs (for a typical value of Q- 105cm-1, lZAw1 atomic spacing). The Iast three columns in table I show the results of the computer experiments. The quantity Ahz,,ss is the dimensionless time interval between nucleation events at steady state, while i,, and OZ$, denote the steady-state ledge spacing and evaporation coefficient [cf. eq. (17)], respectively, approached in those regions of the surface where the effects of curvature [eq. (4)] are negligible (i.e., regions where the dynamic behaviour of concentric discs can be approximated by that of straight, parallel ledges). The uncertainties in the experimental data are _I one unit in the last significant digit unless otherwise indicated. The time development of the surface profile during a typical simulation run to steady-state conditions is shown in fig. 4. Here yVQ/kT=O.Ol and p/p,=OS, hence the nucleation parameters are Qr~=O.O8125 (i.e., approxi-

88

T. SURE& G. M. POUND AND J. P. HIRTH

mately 5.6 times the thermodynamic critical radius Qr*) and Qr,* = 1.02104. The surface profile is shown at times Zi when the ith disc nucleates at the source. The times ri and the number of discs on the surface at that instant (Nledge) are given in the figure. It is seen that a steady-state ledge spacing I,,= 1.42kO.01 develops by the time the 50th ledge nucleates at the source (z~,,). There are 15 discs (approximately +Ql/A,,) on the surface at steady state.

Ass

=

1.42 ? .Ol

4t

i

2 x

3-

0

4

8 DISTANCE,

Fig. 4.

12

16

20

R -

Surface profile at times ~6 when the ith disc nucleates at the source. First disc (denoted by x ) impinges at sink at t = 21.5. Steady state is reached by ~50.

Fig. 5 shows the average evaporation coefficient (a) (or equivalently, the rate of vaporization in eq. (16)) as a function of time for p/p,=O.5 and yVQ/KT=O.OOl, 0.01 and 0.1. Line b in the figure corresponds to the case examined in fig. 4. Typically, (CC) increases approximately linearly with time until the leading disc impinges at the sink, after which (u) approaches its steady-state value at an ever decreasing rate (unless steady state is reached by the time the first disc impinges, as in line c in fig. 5). Variations in (a) about the steady-state value rx,, occur because of transient variations in the ledge spacings near the source and the sink, respectively. The transient and steady-state behavior of the concentric-disc system correlates, in general, to that of the straight, parallel ledges, as can be seen

SPIRAL

DISLOCATION

DYNAMICS

IN CRYSTAL

EVAPORATION

89

t

&

TIME

-

Fig. 5. Average evaporation coefficient versus time for p/p, = 0.5 and yVQ/kT = 0.001, 0.01 and 0.1. Here, x denotes impingement of first disc at sink. Scatter bars denote the extent of the transient variation in (a).

by a comparison of our fig. 4 and fig. 2 of ref. 16+. In order to point out the salient similarities between the two models, whenever appropriate, we now enumerate the main results of the present calculations. (a) The source is uniquely characterized by the system constants yVQ/kT and p/p,, and the ratio r,*/:r, i.e., by the set of variables (yVQ/kT, p/p,, r: and rz). The nucleation parameter I, played a similar role in the case of the straight ledge sources. (b) Steady-state surface profile (L,,) and evaporation rate (LX,,)obtain for a given source independent of the values of Q1/2 and AA, as long as Q1/2X6 and Q1/22 IO&,, analogous to the straight ledge case. In addition a third condition, Q1/2x7VQ/O.O1 kT ensures that there are regions of the surface where capillarity [cf. eq. (4)] is negligible. To some extent, the value of ,&, determines the eventual ASSand IX,, values. This can be seen from table 1, where similar values of Izi, result in approximately the same steady-state ledge spacings, independent of the other system constants. This behavior is analogous to the straight ledge case, where the value of A,, was only a function of the initial ledge spacing At. (c) As in the case of the straight ledges, the sink merely provides the necessary external boundary condition which allows steady state to be reached t It is interesting to note that the initial ledge spacings are approximately figures: rZrn= 0.93979 in fig. 4 and II= 1.O in fig. 2. The steady-state smaller, however, in the case of the concentric discs.

equal in the two ledge spacing is

90

T. SUREK,

G.M.

POUND

AND

J.P. HIRTH

over the entire surface, but does not affect the dynamics of the discs otherwise. The trivial condition 1, < QZ/2 assures the validity of these results. (d) Steady-state conditions obtain near the source first, then they “propagate” toward the sink. Consequently, the steady-state nucleation rate at the source is reached in a time which is much less than that required to obtain steady state over the entire surface. These results are demonstrated by the surface profile plots and the zi values in fig. 4. (e) Based on the result in (d), we can predict the steady-state ledge spacing and evaporation coefficient in those regions of the surface where capillarity is negligible. By equating the steady-state nucleation rate (a,,) to the rate at which ledges pass a given point of the surface, we obtain

%s= The steady-state given by

llbz,,,,

evaporation

(18)

= [2 (1 - p/p,) tan11 (L/2)3/L. coefficient

a,, = (2/L)

in these regions

of the surface

tanh (L/2).

is

(19)

Values of ,I,,, IX,, and o,,, based on the values of ArN,ss obtained in the computer experiments, are given in table 2. The predicted steady-state values are seen to be in good agreement with the experimental values in table 1. TABLE

2

Steady-state predictions for dislocation spiral sources for the model of concentric, circular discs Y VQIkT

PIP,

I

C&S

[es. ;;“s)l

AIt

[eq.(Wl

0.001 0.001 0.001 0.001 0.001

0.001 0.1 0.5 0.9 0.99

0.47 0.49 0.63 1.123 2.914

0.981 0.883 0.484 0.0906 0.00616

0.982 0.980 0.968 0.906 0.616

0.138 0.146 0.209 0.381 1.00

0.01 0.01 0.01 0.01 0.01

0.001 0.1 0.5 0.9 0.99

1.106 1.14 1.43 2.884 28.13

0.908 0.814 0.429 0.0620 0.000711

0.909 0.904 0.858 0.620 0.0711

0.383 0.389 0.49 0.997 16.62

0.1 0.1 0.1 0.1

0.001 0.1 0.5 0.9

2.611 2.753 3.977 27.14

0.6606 0.5755 0.2422 0.00737

0.661 0.639 0.4844 0.0737

0.99 1.035 1.518 16.16

t The quantity A1 = Iss - lr,, represents the diffusion and capillarity-induced of the ledge fronts, as discussed in the text.

accelerations

SPIRAL

DISLOCATION

(f) The time required the following

empirical

DYNAMICS

to reach steady-state relationship

IN CRYSTAL

91

EVAPORATION

conditions

can be estimated

from

:

where the factor n is between 1 (for large A,,) and 10 (for small A,,). Therefore, the transient stage of the ledge dynamics is increased for a low density of and/or large ledge sources (i.e., large Q//2), and for small undersaturations values of y VQ/kT [i.e., conditions which constitute a small back-force effect at the source lg)]. More importantly, the real time-scale of the transient stage is increased for small values of the system constant K [see eq. (15)J which was used to make the time variable dimensionless. Finally, we consider the effects of choosing the ratio rz/rl different from 47c, which was based on an Archimedian-shape approximation for the first spiral turn [cf. eq. (14)]. No such preconception of the spiral shape is necessary in our model, since any numerical ratio r2/r1, together with the functional relationship between the critical disc sizes [eq. (13)], will yield a set of nucleation parameters r: and r: (cf. fig. 3). Since the steady-state dynamics of the discs were found to be closely related to the initial ledge spacing near the source [result (b) above], it is sufficient to consider the effect of r2/r1 on the value of &,. It was found that for large undersaturations (which are the cases of experimental interest) and for small values of y VQ/kT, the effect of the ratio rz/rl was generally small. For example, for p/p, = 0.001, the value of &,, is approximately doubled for an increase in rf/r,* from - 2n: to - 8~ for the range of y VQ/kT values considered in this analysis. For yVQ/kT GO.01, this variation in pi” leads to a change in CI,, (and thus in the rate of vaporization) of less than 5%.

7. Transient and steady-state One of the primary

objectives

shapes of evaporation spirals

of this work was to solve for the steady-state

spiral shape at all undersaturations. In the model of concentric, circular discs, this can be achieved by associating the dynamics of the discs in fig. 2 between two successive nucleation events with the rotation of the cross-section OA in fig. 1 from 8 = 0 to 8 = 27~.A steady-state shape is indicated by the fact that each ledge front in fig. 2 moves from Ri to R,_l during one complete cycle. Fig. 6 shows the first five turns of the spiral that corresponds to the steadystate surface profile in fig. 4t. The spiral was constructed from the kinetic information for the discs for z 54 < t < rs5 ; the angular velocity of the spiral is t As a useful simplification, the center of the spiral is placed at (QYI*, 0) instead of (0,O) as in fig. 1. In this way, the central turn of the spiral is obtained in the same manner as all the other turns.

92

T. SUREK,

-rPe

= 0.5

YVQ kT

= 0.01

sg.

:

20

G.M.

POUND

1800

AND J.P.

HIRTH

X,” =

0.93979

A,,

=

I42

t

=

-&to

co,,=

r

01 T55

0.4292

0”

Fig. 6.

Steady-state

spiral shape for p/p== 0.5 and yVQ/kT= 0.01. Data taken from 554 =Z t < 555 in fig. 4.

given by mang= ~RW,, = 2.697 (in radians per dimensionless unit time). The ledge spacing between successive turns is seen to increase rapidly from li, =0.93979 for 8=0 to its steady-state value A,,= 1.42kO.01. Since y VQ/O.Ol kT= 1.O for this case, the effect of capillarity is negligible beyond the first turn of the spiral. There is, however, overlap of adatom diffusion fields from adjacent spiral turns, since the A-values are less than the spacing I, = 6 (in dimensionless units) which is required for no diffusion overlap43 13). Fig. 6 is typical of the cases examined in table 1: in each case, there is a rapid increase of the ledge spacing from li, to A,,, with A<&, in general, for R
SPIRAL

L 9 yvc) kf

DISLOCATION

DYNAMICS

= 0.5 ii 0.0,

IN CRYSTAL

EVAPORATION

‘in 180’ /

t Cd,,

93

= 0.93979 = T4 tar, = 0.4424

Fig. 7. The shape of the spiral in fig. 6 during the transient stage of evaporation, i.e., s4 < t =G 56.

the spiral during the transient stage of evaporation. In a manner analogous to fig. 6, the shape of the spiral for z,
94

T. SUREK,

C.M.

POUND

AND J.P. HIRTH

8. Discussion

The main purpose of these calculations was to evaluate the role of spiral dislocations in the kinetics of evaporation of otherwise perfect single crystals. Because the model of circular discs is analogous to the treatment of the crystal edge sourceis), we have continually emphasized, in the presentation of the results, the similarities between the dynamic behaviour of the discs and that of the straight, parallel ledges. In this section, we compare our results with the predictions of earlier theories of evaporation spirals, and examine our own predictions regarding the quantitative aspects of the transient and steady-state kinetics of dislocation spirals in crystal evaporation. 8.1. COMPARISON

WITH

THE CABRERA-LEVINEAND

CABRERA-COLEMAN

THEORIES Figs. 8, 9 and 10, respectively, show the dependence on the degree of undersaturation (1 -p/pJ of the steady-state spiral ledge spacing, nucleation frequency and evaporation coefficient in the various spiral models. The subscript ss denotes the values obtained in the present analysis, as given in table 2. The variables in the Cabrera-Levinelo) and the Cabrera-Colemanlz) theories are denoted by the subscripts CL and CC, respectively. The CCvalues are obtained from [cf. eqs. (1) (18) and (19)] Act = 19 Qr;,

(21)

and @CC

=

(1

-

P/P,)

ECC

=

iI2

(f -

~/P~)

tan11

&GQII/&c~

(22)

where r,*= - (y VWllnCn (O)/nJ, and it represents the radius of a discshaped hole which is in unstable equilibrium with the local undersaturation existing at the center of a disc of radius 19 r,* (cf. modification of the CabreraColeman theory in ref. 19). In the above expression, n(O) is given by eq. (IO) with R=O and lt,=l9 r,*. The value of & is obtained from eq. (1); the other CL-values are calculated from a modification of eq. (2), by replacing u, in order to account for the overlap of diffusion with v,tanh(l9Qr*/2) fields from adjacent spiral turns t. In our dimensionless units, we find

Since the various models involve somewhat “arbitrary” numerical constants (i.e., 19 for both CC and CL and 47~for the present work), it is clearly more appropriate to consider the functional dependence of the variables on t The need for the above modifi~tion is apparent (2)-x23 at large undersaturations where r* +O.

when one considers that WCL in eq.

SPIRAL

DISLOCATION

SPIRAL

Fig. 8.

Steady-state

spiral

DYNAMICS

LEDGE

IN CRYSTAL

SPACING,

i

EVAPORATION

-

ledge spacing versus the degree the various models.

of undersaturation

IO0

t

3 5 - lo-’ 2 0 E 2

0

lo-*

2

W

12

2 z

/ IO-? ,

1;

IO-'-3 I Id

I

I

I

lO-2

10-l

IO0

(I-P/Pe)-

Fig. 9.

Steady-state nucleation frequency as a function of the degree of undersaturation in the various spiral models (ss - present work, CL - Cabrera and Levine).

in

96

T. SURE&

G. M. POUND

EVAPORATION

Fig. 10.

AND J. P. HIRTH

COEFFICIENT,

a! -

Evaporation coefficient at steady state versus the degree of undersaturation spiral dislocation sources (ss - present work. CL - Cabrera and Levine).

for

(1 ---p/p,), rather than a direct numerical comparison of the results. Pertinent to our considerations is the value of the back-force effectrg), i.e., the reduction in the local undersaturation at the spiral center because of the overlap of the diffusion field from the first spiral turn. When the back-force effect is small (i.e., small undersaturations and large values of y VQ/LT), we see from fig. 8 that the spiral ledge spacing IX (1 -p/p,)-“, with m= I in all three models. The value of the exponent m decreases as the back-force effect increases; for example, at large undersaturations, A.,,N &cc (1 -p/p,)- %for large rVQ/kT, and &s~ACCcc (1 -p/p,)-* for small yVQ/kT. The value of JCL is meaningless at large undersaturations; the result of the back-force considerations is that 1 converges to some finite value as p/p, -+ 0 (see I,, and d,$ in fig. 8). The dependence of the steady-state nucleation frequency (and therefore the angular velocity or the rate of evaporation) on the degree of undersaturation is shown in fig. 9. The CC-values [eq. (22)] almost coincide with those in the present work, and are therefore not shown in the figure. Two types of functional dependences are apparent in fig. 9 : firstly, WCK(1 -p/p,) when the

SPIRAL

DISLOCATION

DYNAMICS

IN CRYSTAL

EVAPORATION

97

back-force effect is large; and secondly, woe (1 --p/p,)” when the back-force effect is small. These cases correspond to the familiar “linear” and “parabolic” growth (or evaporation) rates, respectively, in the earlier theories. The close agreement between w,, and wcL resuits from our modification of the CL-theory [cf. eq. (23)]. In general, we thus find that: (i) there is agreement between our results and CC at all undersaturations and for all values of the pertinent material parameters (i.e., ~~~/~~); and (ii) there is agreement among al1 the theories at very small undersaturations only. The significant result is that, in general, it is incorrect to apply the relationship A= 19r” to dislocation spiral sources. The results of the back-force calculations (and the quantitative extension of the CC-theory)lQ) show that the radius of curvature at the spiral center (rt) can be several orders of magnitude greater than r*, particularly for small values of yVQ/kT and large undersaturations. In addition, we believe that there is no apparent significance in the number 19 when the back-force effect is appreciable; the ratios &,/Qr: in tables 1 and 2 vary between 16 and 29. Part or all of this variation may be attributed, however, to the arbitrary choice of t-,*/r F = 4n in our model. 8.2. T~ESHAFESOFEVAPORATIO~AN~GROWTHSPIRA~ There is no exact theory to compare to the transient development of the spiral shape in the present calculations. As one would expect, the time dependence of the surface profile in fig. 4 is consistent with the qualitative predictions of the continuum kinematic theoryQQ-24). The theory in its present form, however, is not applicable to the ledge fronts near the source where capiliarity is not negligible. This conclusion was demonstrated experimentally in the dissolution of lithium fluoride by Ives and HirthQs). The steady-state spiral shape is consistent with the physical expectation that dA/dR should be greater than 0 near the spiral center when the ledge kinetics are adatom di~usion-controlled~6). The spreading of the spiral ledge fronts takes place partly because of the diffusion induced acceleration of the ledges and partly because of the effect of decreasing capillarity away from the spiral center. The quantity AL?LIz=Iz,, --&, in table 2 is a measure of the total ledge acceleration; the relationship between Ah and /lin is not as straightforward, however, as in the analogous case of the straight, parallel ledges (cf. fig. 3 in ref. 16). Finally, in a recent paper, Keller2’) presented experimental evidence for the transient shapes of monatomic spirals during the evaporation of NaCl single crystals. He determined the &distributions of the spirals after various evaporation times, and presented the results in the form of our fig. 4. These experiments not only give a unique verification of our computer simulation

98

T. SUREK,

G. M. POUND

models, but also provide quantitative crystal growth and evaporation. 8.3. COMPARISONWITHTHE

AND J. P. HIRTH

support for the TLK mechanism of

HIRTH-POUNDTHEORY

For much the same reasons as in the case of the straight ledges, the present calculations show that the Hirth-Pound asymptotic steady-state behaviorla) is not exhibited by ledges originating from spiral dislocation sources. For li, < &, where i, = 6 represents the ledge spacing for which diffusion overlap is negligible4* Is), the initial tendency of the leading ledge spacing is indeed to increase toward Jo (see disc denoted by x in fig. 4). This is the result of the diffusion and capillarity-induced accelerations, as discussed above. The physical requirement for a finite system (i.e., some sink-condition for the discs), but more importantly, the propagation of steady-state conditions away from the source, precludes the situation in which an appreciable portion of the surface is covered by ledge spacings -Jo. The trivial exception to this is some particular choice for the source parameters which would lead to &-,I, (from table 1, it appears that il,, N 6 for R,, N 4). For li” > Jo, Hirth and Pound predict that the evaporation kinetics will be dominated by ledges from crystal edge sources. The computer simulations show that sources with the highest nucleation frequency w,, will control the eventual steady-state dynamics 1s). Hirth and Pound are correct in their assessment of the region affected by capillarity: for our dynamic disc system at steady state, R.<;1,, for R 5 5 y VQlO.01kT. Also, their formulation of the spiral kinetics is essentially correct; indeed, we have adopted their model of concentric, circular discs in the present paper. The source of error in their application of their infinite surface treatment to a finite crystal lies in the concept of a “perturbed” region in which h+l,. This asymptotic steady-state assumption leads to erroneous conclusions regarding the theoretical value of the evaporation coefficient CI,in eq. (16) for a finite system. As in the case of the straight ledges, the results of the present calculations show that the Hirth-Pound value 4.13) of 01,= 3 has no general significance as a limiting law in the evaporation of real crystals. As shown in table 2, any value of 0 < c(,< 1 is possible under steady-state conditions in the present model. In addition, we differ from Hirth and Pound’s order-of-magnitude estimaters) that no steady-state evaporation rate obtains when the ledge spacings are < 2, near the spiral center. The present calculations show that steady-state spiral shapes exist even at very large undersaturations; in fact, the spiral shape is convergent when p/p, + 0 (cf. ASScurves in fig. 8). Also, it is apparent that the transient information for the p/p, = 0.001 cases in table 1 can be interpreted as adequately representing the kinetics for pJp,=O. The

SPIRAL DISLOCATION

concept

of a steady-state

DYNAMICS

IN CRYSTAL

EVAPORATION

spiral shape may break down, however,

99

for cases

when the strain energy of the spiral dislocation is large. Under these conditions, there is no barrier to the successive nucleation of discs in our model, and a hollow, macroscopic pit will resultrg). The present calculations can be used to delineate the conditions for observing low values of a, when the evaporation kinetics are controlled by spiral dislocation sources of ledges. As seen from fig. 10, evaporation coefficients which are appreciably less than unity can be expected for spiral sources with large values of the parameter yVQ/kT and also for small undersaturations. In addition, a low value of c(, can be measured during the transient stage of evaporation (cf. fig. 5). As shown in section 6, the real time-scale of the transient stage is increased for small values of the system constant K in eq. (15), i.e., typically for materials with low values of the equilibrium vapor pressure pe, Finally, the present results, which are approximate for spiral dislocation sources, are exact for the problem of evaporation proceeding by the repeated nucleation of disc-shaped holes on some preferred surface site, such as a dislocation (edge or screw)-surface intersection. The present calculations are therefore pertinent to the understanding of the formation and observability of etch pits on such sites. The steady-state surface slopes in our model are given by Q/z/&,, where h is the height of the spiral steps. Typically, none of the cases examined in table 1 correspond to the steep-sided etch pits which are observed. The inclusion of the dislocation strain energy in our model [parametrically represented by rz/rl <4x in eq. (14)], however, leads to a decrease in the value of Ain (and therefore n,,) and hence to steeper surface profiles. Higher index planes can also result if the spacings of the spiral dislocation sources were such as to violate any of the conditions on the parameter QZ/2 [cf. result (b) in section 61, in which case the steady-state ledge spacing would not be reached over the entire surface. For very steep pits, of course, back-flux effects in the vapor phase would also have to be considered 2s). 9. Summary A computer simulation technique, involving a model of concentric, discshaped holes, was used to approximate the transient-to-steady state kinetics of ledges emanating from spiral dislocation sources in crystal evaporation. This approximation yielded a tractable solution to the otherwise very complicated problem of dislocation spiral kinetics. The exact analysis of the true spiral shape at all undersaturations must eventually involve a numerical solution of the time-dependent diffusion equations with moving boundary conditions.

loo

T. SVREK,

G.M.

POUND

AND

J.P.

HiRTH

In terms of a comparison of our results with the predictions of earlier theories of spiral evaporation, the significant findings of the present work are : (i) the prediction of a steady-state spiral shape at all undersaturations; (ii) the conclusion that the relationship a= 19r * [cf. eq. (l)] cannot be generally applied in problems involving dislocation spiral sources; and (iii) the demonstration that the evaporation coefficient can have any value between 0 and I (rather than the limiting value of a”=$ found in the Hirth-Pound analysis) for the case when spiral dislocations control the evaporation kinetics. Inasmuch as our results are in dimensionless, reduced units, they apply to any actual material for which the parameters in the system constants yVQ/kT, p/pe and K are known. Although we were mainly concerned with surface diffusion kinetics as the rate-controlling mechanism in ledge-wise evaporation, the present results and models can be readily applied to give the ledge dynamics and surface diffusion contribution to the evaporation coefficient in cases where the ledge dynamics are affected by other processes as well. The usefulness of these computer simulations to study the formation and stability of macroscopic steps comprised of many interacting monatomic ledges is also apparent; recently, such models were used to study the modes of oscillation which can occur during the perturbation of monatomic step trains in crystal growth2Q). Acknowledgements

This work was performed at Stanford University. Financial support from the National Science Foundation and the Stanford Center for Materials Research (T.S., G.M.P.), the Office of Naval Research (J.P.H.) and the National Research Council of Canada (T.S.) is gratefully acknowledged. The authors are grateful to the IBM Scientific Research Center in Palo Alto for the use of the Center’s computer facilities, and to Dr. F. F. Abraham for illuminating discussions. References 1) F. C. Frank, DiscussionsFaraday Sot. 5 (1949)48, 67. 2) W. Kossel, Nachr. Ges. Wiss. Gijttingen (1927) 135. 3) W. K. Burton, N. Cabrera and F. C. Frank, Phil. Trans. Roy. Sot. London A 243 (1950) 299. 4) J. P. Hirth and G. M. Pound, J. Chem. Phys. 26 (1957) 1216. 5) E. D. Dukova, Dokl. Akad. Nauk SSSR 121 (1958) 288; Soviet Phys. Dokl. 3 (1958) 703. 6) A. R. Verma, Crysiaf Growrh and ~~sloca~io~s (Butterworths, London, 1953). 7) A. J. Forty, Advan. Phys. 3 (1954) 1. 8) H. Bethge, in: Molecular Processes on Solid Surfaces, Eds. E. Drauglis, R. D. Gretz and R. I. Jaffee (McGraw-Hill, New York, 1968) p. 569. 9) K. W. Keller, Phys. Status Solidi 36 (1969) 557.

SPIRAL

DISLOCATION

DYNAMICS

IN CRYSTAL

EVAPORATION

101

10) N. Cabrera and M. M. Levine, Phif. Mag. 1 (1956) 450. 11) I. Langmuir, Phys. Rev. 2 (1913) 329. 12) N. Cabrera and R. V. Coleman, in: The Art and Science of Growing Crystals, Ed. J. J. Gilman (Wiley, New York, 1963) p. 3. 13) J. P. Hirth and G. M. Pound, Acta Met. 5 (1957) 649. 14) L. D. Hulfett, Jr. and F. W. Young, Jr., J. Electrochem. Sot. 113 (1966) 410. 15) U. Bertocci, Surface Sci. 15 (1969) 286. 16) T. Surek, G. M. Pound and J. P. Hirth, J. Chem. Phys. 55 (1971) 51.57. 17) T. Surek, Ph. D. Thesis, Stanford University, 1972. 18) W. Oldfield, G. T. Geering and W. A. Tiller, Mater. Sci. Engng. 2 (1967) 91. 19) T. Surek, J. P. Hirth and G. M. Pound, J. Crystal Growth 18 (1973) 20. 20) W. W. Mullins and J. P. Hirth, J. Phys. Chem. Solids 24 (1963) 1391. 21) P. Bennema and G. H. Gilmer, in: Theory of Crysral Growth, Proc. Summer School on Crystal Growth, Noordwijkerhout, June 1971 (North HoUand, Amsterdam, 1973). 22) F. C. Frank, in: Growth and Perfection of Crystah, Eds. R. H. Doremus, B. W. Roberts and D. Turnbuli (Wiley, New York, 1958) p. 41 I. 23) N. Cabrera and D. A. Vermilyea, ref. 22, p, 393. 24) A. A. Chernov, Dokl. Akad. Nauk SSSR 117 (1957) 983. 25) M. B. Ives and J. P. Hirth, J. Chem. Phys. 33 (1960) 517. 26) T. Surek, J. P. Hirth and G, M. Pound, Rost Kristallov 11/12, Proc. of the P. K. IV Conf. on Crystal Growth, Tsakhkadzor, Armenian SSR, September 1972, in press. 27) K. W. Keller, Rost Kristallov 11/12, in press. 28) J. P. Hirth and G. M. Pound, Progr. Mater. Sci. 11 (1963) 1. 29) P. Bennema and R. van Rosmalen, Rost Kristallov 11/12, in press. 30) F. B. Hildebrand, Introduction to Numerical Analysis (McGraw-Hill, New York, 1965)