An integral equation method for epitaxial step-flow
growth simulations
Jingfang Huang
a,*, Ming-Chih Lai
b, Yang Xiang
ca
Department of Mathematics, University of North Carolina, CB#3250, Phillips Hall, Chapel Hill, NC 27599, United States b
Department of Applied Mathematics, National Chiao Tung University, Hsinchu 30050, Taiwan c
Department of Mathematics, Hong Kong University of Science and Technology, Clear Water Bay, Kowloon, Hong Kong Received 14 June 2005; received in revised form 12 December 2005; accepted 3 January 2006
Available online 14 February 2006
Abstract
In this paper, we describe an integral equation approach for simulating diffusion problems with moving interfaces. The solutions are represented as moving layer potentials where the unknowns are only defined on the interfaces. The resulting integro-differential equation (IDE) system is solved using spectral deferred correction (SDC) techniques developed for gen-eral differential algebraic equations (DAEs), and the time dependent potentials are evaluated efficiently using fast convo-lution algorithms. The numerical solver is applied to the BCF model for the epitaxial step-flow growth of crystals, for which the solutions are calculated accurately instead of using quasi-static approximations. Numerical results in 1 + 1 dimensions are compared with available results in the literature.
2006 Elsevier Inc. All rights reserved.
Keywords: Integral equation; Potential theory; Diffusion equation; Jump conditions; Epitaxial step-flow; Spectral deferred correction; Fast algorithms
1. Introduction
Many problems in science and engineering require the numerical solution of diffusion or diffusion type equations with moving interfaces. Examples include the study of crystal growth in semiconductor manufac-turing and simulations of diffusion–reaction process in animals and plants. For problems of this kind, most existing schemes are based on finite difference or finite element methods. In this paper, as an alternative approach, we present an efficient and accurate scheme using integral equation ideas.
Integral equation methods (IEM) have traditionally been neglected as a general numerical scheme. In these methods, the solution is represented as convolutions of the underlying Green’s function (fundamental solu-tion) with given or unknown density functions, which are defined on the interface (layer potentials) or inside the computational domain (volume potentials). Direct evaluations of these volume and layer potentials require
0021-9991/$ - see front matter 2006 Elsevier Inc. All rights reserved. doi:10.1016/j.jcp.2006.01.006
* Corresponding author. Fax: +1 919 962 9345.
E-mail addresses:[email protected](J. Huang),[email protected](Y. Xiang).
a prohibitive amount of work: when n discretization points are used in a numerical quadrature, usually O(n2) operations are required. However, the situation is changing rapidly due to the introduction of fast convolution algorithms. The first such algorithm is the fast Fourier transform (FFT) introduced by Cooley and Tukey in
1965[10]. For uniform grid points, the total amount of work is reduced to O(n log n). Starting from 1980 or so,
several new techniques have been developed for ‘‘arbitrary’’ grid distributions, including the O(n log n) FFT based methods (e.g., particle mesh Ewald (PME), particle–particle particle–mesh (P3M), and precorrected
FFT (pFFT)[13,38,54]) and multipole expansion based methods (e.g., O(n log n) tree code[2,5]and the
asymp-totically optimal O(n) fast multipole method (FMM) [24,27]). With the fast algorithm accelerations, integral
equation approach has been proven extremely powerful for large scale scientific problems, especially in the fields of electromagnetics and fluid dynamics. Compared with finite difference and finite element methods, integral equation methods have several advantages: they are insensitive to the complexity of geometry; with proper formulation, they are unconditionally stable; and they easily allow for adaptive mesh refinement and parallel computing.
Enormous progress has also been made for the diffusion equation. In[29], an algorithm for solving the pure
initial value problem in d + 1 dimensions
Utðx; tÞ ¼ DU ðx; tÞ for x 2 Rd; t >0;
Uðx; 0Þ ¼ f ðxÞ for x2 Rd
(
was presented by Greengard and Strain. To evaluate the ‘‘initial potential’’ Uðx; tÞ ¼ ð4ptÞd=2
Z Rd
ejxyj2=4tfðyÞ dv
y; ð1Þ
their ‘‘fast Gauss transform’’ (FGT) requires O(N + M) work to determine U(x, t) at N points given the initial
data f(y) at M points, while direct method requires O(NM). In[28], an algorithm for the rapid evaluation of
the history dependent heat potentials was developed. The potentials were divided into a ‘‘history’’ part (rep-resenting the influence of the density at distant times) and a ‘‘local’’ part (reflecting the influence of the density over the most recent time). The history part is updated using Fourier series, and the local part is evaluated using high order product rules. More recently, a new version of the fast Gauss transform was introduced
by Greengard and Sun in [30].
In this paper, we apply integral equation approach to diffusion problems with moving interfaces, in
partic-ular, the epitaxial step-flow growth of crystals described by the BCF model[9,58,70]. In this model, the crystal
surface consists of terraces with different heights separated by atomic-height steps. The adatoms are deposited onto the crystal surface, and then diffuse on the terraces until they meet and get incorporated into the steps.
The diffusion process is modeled by the diffusion equation /t= DD/ + F, where /(x, t) is the adatom density,
x2 Rd, d = 1 or 2, t2 R, D is the diffusion constant, and F(x, t) is the deposition flux which is usually assumed
to be a constant. The incorporation of adatoms to the steps is described by the interface conditions defined only on the steps. Traditionally, this BCF model was solved by the quasi-static approximation in the regime where the deposition rate of the adatoms onto the surface is much less than the adatom hopping rate on the
terraces[9,58,70], or by the finite difference method for problems with equidistant steps when the model can be
transformed into a diffusion process on a terrace with fixed boundaries[48,53,70]. Recently, several phase field
models were proposed in[41,46,50,55,63], simulation methods based on the level set framework[49]were
pro-posed in[11,12], and a finite element method was developed in[4]. However, to our knowledge, very limited
integral equation solvers have been proposed in the literature. The Green’s function representation was used
by Liu and Metiu[46]for linear instability analysis for uniform step trains, and an integral equation method
coupled with level set framework was developed by Sethian and Strain [65] to simulate solidification-type
problems.
In our integral equation approach, by introducing the moving layer potentials, the original BCF model is transformed into an integro-differential equation (IDE) system where the unknown density functions are only defined on the steps between different terraces, hence dramatically reduce the total number of unknowns. The layer potentials are efficiently evaluated using fast convolution techniques, including those introduced by
Greengard and Strain in[28]. For the IDE system, most existing solvers are based on either the backward
have also been tested. In[15], spectral integration is coupled with Gaussian quadratures and deferred/defect corrections. The resulting spectral deferred correction (SDC) methods turn out to be very competitive with best existing ODE initial value problem solvers. SDC methods were also applied to partial differential
equa-tions (PDEs) and differential algebraic equaequa-tions (DAEs) in[7,43,47,61,62]. In this paper, we generalize SDC
techniques to IDE systems. Preliminary numerical results show that higher order accuracy can be obtained for problems with moving interfaces.
This paper is organized as follows. In Section2, we present the mathematical model for epitaxial step-flow
growth, which is a representing instance of general diffusion problems with moving interfaces. In Section3,
starting from jump conditions for general moving layer diffusion potentials, we formulate the IDE system
for the BCF model. In Section4, we present a spectral deferred correction scheme for integrating general
DAEs and show how it can be generalized to IDEs. The efficient evaluation of layer potentials is also dis-cussed. One difficulty in the numerical simulation is that when two steps are close to collision, adaptive
tem-poral discretization is required to correctly capture the local dynamics. In Section5, by studying solutions near
collision, we present ‘‘reduced order’’ local collision model to overcome this difficulty. Numerical simulation
results in 1 + 1 dimensions are presented in Section6.
2. BCF model for epitaxial step-flow growth of crystals
Diffusion type equations with moving interfaces have been used to model different biological, chemical, and physical processes including drug delivery, tumor growth, fluid–solid interactions, and manufacturing process of micro-electro-mechanical systems (MEMS). A specific example is the epitaxial growth of crystals.
Epitaxial growth is the growth of crystalline film on a crystalline substrate following the same structure as the substrate. The crystal surface consists of terraces separated by atomic-height steps. Adatoms are deposited onto the surface, and diffuse on the terraces until they meet and get incorporated into the steps.
As a result, the steps move forward and the crystal grows, as shown inFig. 1. This growth mode is referred
to as the step-flow growth, which can be described by the BCF model proposed by Burton et al.[9,58,70].
The best way to grow a crystal is to grow it on an infinite monotonic surface with parallel, equidistant steps
[58,70].
Detailed description of the BCF model can be found in[58,70]. In this section, for simplicity of
discus-sion, we focus on the model in 1 + 1 dimensions, in which the steps are straight and parallel, and diffusion on the terrace is uniform along the direction of steps. Without loss of generality, in the remainder of this section, the surface is assumed to be infinite, and the height of each step is assumed to be +1, where the height of a step is defined as the height of the terrace on the left minus the height on the right (in atomic-height units). The atomic-height is positive if the left is higher than the right, and negative otherwise. Letting
deposition
diffusion attachment
terrace
step
{wj(t)} with < wj1(t) < wj(t) < wj+1(t) < be the locations of steps at time t, the equations describing the adatom diffusion and step motion can be written as
/t¼ D/xxþ F ; wjðtÞ < x < wjþ1ðtÞ; D/x kþ/¼ 0; x¼ wþjðtÞ; D/xþ k/¼ 0; x¼ w jþ1ðtÞ; w0jðtÞ ¼ a2 kþ/j x¼wþj þ k /j x¼wj ; 8 > > > > > < > > > > > : ð2Þ
where /(x, t) is the adatom density on the terraces, F is the deposition flux which is assumed to be a constant,
D is the diffusion constant on the terraces, a is the lattice constant, k+and kare the hopping rates of an
ada-tom to the upward step and downward step, respectively.
The first equation in system (2)describes the deposition and diffusion processes of adatoms on a terrace
between the jth and (j + 1)st steps. The second and third equations describe the incorporating process of adatoms into the steps, which serve as boundary conditions of the diffusion problem. Note that there is
a jump in the adatom density /(x, t) at a step wj(t). The fourth equation gives the velocity of each step.
Usually the adatoms prefer to go to the upward step rather than to the downward step, since an energy
barrier exists near the downward step. This is called the Schwoebel barrier [64]. The Schwoebel effect
stip-ulates that
kþ P k. ð3Þ
In this paper, we present an integral equation approach for solving the BCF model(2)as an alternative to the
quasi-static approximation and other methods mentioned in the introduction section. We shall compare our simulation results with available analytical and numerical results, which are reviewed briefly below.
In the literature, the quasi-static approximation /t 0 was used when the deposition process is much
slower than the diffusion process, i.e., when the Peclet number
Pea
2l2F
D 1; ð4Þ
where l is the average spacing between steps[58,70]. Under the quasi-static approximation, the original BCF
model (2)is simplified and analytical solution can be obtained. The step velocity in this regime is
w0jðtÞ ¼ a2Fljþ lj1 2 þ a2Fl j 2 1 kk1þ 1 kþþ 1 kþ lj D a2Fl j1 2 1 kk1þ 1 kþþ 1 kþ lj1 D ; ð5Þ
where lj= wj+1 wjis the spacing between the jth and (j + 1)st steps. The readers are referred to [58,70]for
more results in this regime. Some recent results can be found in, e.g.,[14,17,42,56,59,60,69,72,73].
More physics can be incorporated in the BCF model(2), such as the nucleation of new steps[17,60], and the
step–step and adatom–step elastic interactions[14,69,72,73]. These can also be handled within the framework
of our integral equation method. We shall present an example which includes the deterministic nucleation
model of Elkinani and Villain, and Politi and Villain[17,60], in which a new terrace with width 2a is nucleated
at the center of a top terrace when its width is greater than a critical value lc. In the regime Pe 1, when
k+= +1 (i.e., no barrier for the adatoms to get incorporated into an upward step), lc satisfies
l3cðlcþ 6lsÞ ¼ 12D=F , where ls= D/k a is a measure of the strength of the Schwoebel effect [17,60].
In the regime Pe 1, the deposition process and the diffusion process on the terraces are comparable, and
the steps move fast[22,23,37,46,48,53,70,71]. In this regime, the BCF model(2) predicted oscillations in the
step velocity[37,53,70]. These oscillations were suggested as a possible mechanism for the intensity oscillations
observed in reflecting high-energy electron-diffraction (RHEED) studies of growing surfaces[53], and
condi-tions under which these oscillacondi-tions occur were investigated[37]. However, nucleation of new terraces is
impor-tant in this regime due to the relative high adatom density compared with that in the regime Pe 1. Research
results on both the step velocity oscillations and the nucleation indicated that nucleation is more likely to be
the primary mechanism for this RHEED intensity oscillation [37,48]. The nucleation process depends on
longer valid in this regime. Traditionally, problems in this regime were studied numerically for equidistant step arrays, so that the problems can be transformed into a diffusion process on a terrace with fixed boundaries
[48,53,70], and available nucleation models in this regime used the rate equations for the total number of
terraces or islands coupled with the adatom density[11,12,48]. In this paper, we shall focus on the integral
equation method for solving the full BCF system(2), and examine the method by comparing the simulation
results with available results on the step velocity oscillations without nucleation[37,53,70]. Since the adatom
density is solved accurately and efficiently using our integral equation method, nucleation model based on adatom density in this regime within the framework of the integral equation method can be developed without essential mathematical difficulties. This involves more physics and detailed discussions will be presented in the future.
3. Integral equation formulations
In this section, we study analytical properties of general moving layer diffusion potentials, the results in
1 + 1 dimensions are shown in Theorems 1 and 2. Using moving layer potentials, the original BCF model
can be reduced to an IDE system where the unknowns are only defined on the steps. 3.1. Jump conditions of moving layer potentials
Traditional potential theory solves the diffusion equation with stationary boundary conditions in d spatial
dimensions using a combination of the initial potential to satisfy the initial condition (see Eq.(1)), the volume
potential for the source term, and the layer potentials for different boundary conditions[32,44]. The purpose
of this section is to study the moving layer potentials and their jump properties.
Assuming C(t) is the time dependent surface of a domain in Rd, for a given density function q(y(t), t),
y(t)2 C(t), the single layer potential is defined for any x 2 Rdas
Sqðx; tÞ ¼ Z t
0 Z
CðsÞ
Gðx yðsÞ; t sÞqðyðsÞ; sÞ dsyðsÞds; ð6Þ
where sy(s) denotes the surface integral over C(s) and G(x, t) is the free space Green’s function for the heat
equation given by
Gðx; tÞ ¼ ð4ptÞd=2ekxk24t . ð7Þ
In many computations, one may prefer the periodic Green’s function instead of the free space one. The peri-odic Green’s function can be represented either in the Fourier domain using Fourier series or in the physical domain by the method of images. When d = 1, for the interval [p, p], the Fourier domain representation is given by Kðx; tÞ ¼ 1 2p X1 k¼1 ek2t eikx; ð8Þ
while in physical domain, it is the lattice sum of G(x, t) in Eq.(7)given by
Kðx; tÞ ¼ 1ffiffiffiffiffiffiffi 4pt p X 1 k¼1 eðx2pkÞ2=4t. ð9Þ
The equivalence of these two representations is a particular instance of the Poisson summation formula[16].
Similarly, the double layer potential with density l(y(t), t) (y(t)2 C(t)) is given by
Dlðx; tÞ ¼ Z t 0 Z CðsÞ o omyðsÞ
Gðx yðsÞ; t sÞlðyðsÞ; sÞ dsyðsÞds; ð10Þ
For completeness of concepts, for a given density function f(x, t), the volume potential is defined by Vfðx; tÞ ¼ Z t 0 Z XðsÞ Gðx y; t sÞf ðy; sÞ dvyds; ð11Þ
where dvyrepresents the volume integration over the time dependent physical domain X(t) in Rd. The volume
potential satisfies the equation ðVfðx; tÞÞt¼ DVfðx; tÞ þ f ðx; tÞ.
When d = 1, the following jump conditions for single and double layer potentials defined on a moving inter-face can be derived.
Theorem 1. Assume the interface at time t locates at y(t), the single-layer potential Sqðx; tÞ ¼R
t
0Gðx yðsÞ;
t sÞqðyðsÞ; sÞ ds is continuous "x 2 R and t P 0, and satisfies
(1) o
otSqðx; tÞ ¼ o2
ox2Sqðx; tÞ 8x 6¼ yðtÞ, t > 0, (2) Sq(x, 0) = 0"x 6¼ y(0),
(3) limx!yðtÞoxoSqðx; tÞ limx!yðtÞþ o
oxSqðx; tÞ ¼ qðyðtÞ; tÞ.
In the theorem, the + sign gives the limit as x approaches y(t) from the direction x > y(t), and the sign is
the limit from the other direction x < y(t). For double layer potentials, we have
Theorem 2. Assume the interface at time t locates at y(t), the double-layer potential Dlðx; tÞ ¼
Rt 0 o
oyGðx yðsÞ;
t sÞlðyðsÞ; sÞ ds is defined "x 2 R and t P 0, and satisfies
(1) o
otDlðx; tÞ ¼ o2
ox2Dlðx; tÞ "x 6¼ y(t), t > 0. (2) Dl(x, 0) = 0"x 6¼ y(0),
(3) limx!y(t)Dl(x, t) limx!y(t)+Dl(x, t) =l(y(t), t), (4) limx!yðtÞ oDlðx;tÞ ox limx!yðtÞþ oDlðx;tÞ ox ¼ lðyðtÞ;tÞ 2 y 0ðtÞ.
The proof of above theorems uses the simple but tedious local Taylor expansions and hence we neglect the details.
For standard boundary value problems where C is independent of time, above theorems are reduced to the
well-known jump conditions for stationary (y0(t) = 0) layer potentials as in the following corollary, which can
be found in standard textbooks including[32]:
Corollary 1. In one dimensional space, when x approaches a point x0on the interface, we have
(1) [Sq(x, t)] = 0, (2) o omxSqðx; tÞ h i ¼ qðx0; tÞ, (3) [Dl(x, t)] =l(x0, t), (4) o omxDlðx; tÞ h i ¼ 0.
Here [Æ] represents the jump across the boundary from x < x0to x > x0.
Jump conditions for higher order derivatives and higher dimensions for both single and double layer poten-tials can be derived similarly. The detailed results with applications will be reported at a later time.
3.2. Integro-differential equation system
In this section, we focus on the BCF model for epitaxial crystal growth. Using both single and double layer potentials, we show how the original partial differential equations can be reduced to an IDE system.
For simplicity, we first non-dimensionalize model system(2)by rescaling length by L0(the length scale we
are interested in and contains many steps) and time by T ¼ L2
0=D. Let /0= a2/be the coverage of a lattice site
(average number of atoms per atom site), F0¼ a2L20F =D be the dimensionless deposition flux, and
k0 ¼ kL0=Dbe the dimensionless hopping rates, the model system then becomes (still using x, t, /, wj, F,
and k±for the dimensionless variables and quantities)
/t¼ /xxþ F ; wjðtÞ < x < wjþ1ðtÞ; /x kþ/¼ 0; x¼ w þ jðtÞ; /xþ k/¼ 0; x¼ wjþ1ðtÞ; w0jðtÞ ¼ kþ/j x¼wþ j þ k /j x¼w j. 8 > > > > < > > > > : ð12Þ
This system contains three dimensionless parameters F, k+and k. Assuming periodic boundary conditions
for an interval of length L containing N steps, the Peclet number (see Eq.(4)) in the dimensionless system
be-comes Pe¼ l2F, where l¼ L=N is the average spacing between steps.
The first step to solve this system is to represent the solution /(x, t) as /ðx; tÞ ¼ Ft þ Z L 0 Kðx y; tÞ/ðy; 0Þ dy þX N j¼1 Z t 0 Kðx wjðsÞ; t sÞqjðsÞ ds þ Z t 0 oKðx wjðsÞ; t sÞ owj ljðsÞ ds " # . ð13Þ
Here, K(x, t) is the corresponding periodic Green’s function, qj(t) = q(wj(t), t) and lj(t) = l(wj(t), t) are the un-known density functions at time t. For non-periodic boundary conditions, additional layer potentials are
added for the left and right boundaries, respectively. Note that using Theorems 1 and 2and properties of
the initial potential, this representation automatically satisfies the first equation in(12). Applying jump
con-ditions to Eq.(13), the system to be solved becomes
w0j¼ kþ/þþ k/ ¼ /þ x / x ¼ w0j 2 ljþ qj " # ; ð14Þ
with boundary conditions for the jth step lim x!w j ð/xþ k /Þ ¼ 0; ð15Þ lim x!wþ j ð/x kþ/Þ ¼ 0. ð16Þ
We refer to this formulation (Eqs.(14)–(16)) as an IDE system. The unknowns are the density functions qj(t)
and lj(t), and the step locations wj(t). One obvious advantage of this formulation is that the unknowns are
only defined on the steps, not the interior of the terraces. Such formulations are usually referred to as surface integral formulations in the integral equation methods literature. We postpone more explicit descriptions of
these equations to Section4, where spectral deferred correction techniques are introduced to derive higher
or-der in temporal direction. 4. Numerical techniques
In this section, we discuss the efficient solution of the IDE system(14)–(16). We focus on two ideas: (a) how
to derive higher order accuracy in temporal direction by iteratively refining the approximate solution using a first order method to the ‘‘error equation’’, and (b) how different potentials can be efficiently evaluated. 4.1. Spectral deferred correction methods
The basic idea of deferred and defect correction methods is to use a low order method and iteratively
quadrature with Picard integral equation, Dutt et. al. proposed the spectral deferred correction methods (SDC) for ODE initial value problems. Preliminary results show that the resulting numerical methods are comparable to best existing schemes in efficiency but with better accuracy and stability properties. In this sec-tion, we generalize this technique to differential algebraic equations (DAEs).
For simplicity, we restrict our discussion to a general DAE initial value problem HðyðtÞ; y0ðtÞ; tÞ ¼ 0;
where y(0) = y0and y0ð0Þ ¼ y00 are specified initially and satisfy Hðy0; y00;0Þ ¼ 0. Suppose we want to march
from time t0= 0 to Dt, and p Gaussian nodes are given by t1, t2, . . . , tpin [0, Dt]. The spectral deferred
correc-tion scheme works by searching for a polynomial approximacorrec-tion y0(t) P(t) such that
H y0þ Z ti 0 PðsÞ ds; P ðtiÞ; ti ¼ 0 for i¼ 1; 2; . . . ; p. ð17Þ
This formula is often referred to as the collocation or pseudo-spectral formulation. Because of the use of Gaussian nodes, P(t) can be constructed as a Legendre polynomial expansion where the coefficients are
com-puted using Gaussian quadrature. Also,Rti
0 PðsÞ is derived using spectral integration[25]hence the numerical
unstable differentiation operator is avoided and it is possible to derive high order methods.
The deferred correction methods first assume an approximate solution exists and is denoted by P[0](t).
Introducing the error dðtÞ ¼ P ðtÞ P½0ðtÞ; the error equation is defined as
H y0þ Z ti 0 P½0ðsÞ ds þ Z ti 0 dðsÞ ds; P ðtiÞ þ dðtiÞ; ti ¼ 0.
A low order scheme is then applied to derive di= d(ti) at Gaussian nodes, i.e., using implicit Euler method and
solving equations H y0þ Z ti 0 P½0ðsÞ ds þX i j¼1 ðtj tj1Þdj; PðtiÞ þ di; ti ! ¼ 0 ð18Þ
for i = 1, . . . , p. The approximate solution is then updated using P½0ðtÞ ¼ P½0ðtÞ þ dðtÞ;
where d(t) is the interpolating polynomial of (ti, di) at Gaussian nodes. This procedure is repeated until the
error is below a prescribed tolerance.
Pseudo-Code: Spectral Deferred Correction Method Comment [Initial approximation]
Initially, y0(t) is approximated by a constant function P½0ðtÞ ¼ y0
0. Comment [Compute successive corrections.]
while id(t)i > tol, do
(1) Use a low order method and solve the discretized error Eq.(18).
(2) Use Gaussian quadrature and compute the interpolating polynomial d(t).
(3) Update the approximate solution P[0](t) = P[0](t) + d(t).
end do
4.2. SDC Methods for IDE Systems
Generalization of SDC methods to Eqs. (14)–(16)is straightforward. In this section, we discuss the error
In a time marching scheme, assume all density functions and locations of steps are derived at previous times
up to t Dt, we can rewrite the solution as
/ðx; tÞ ¼ t þ Z L 0 Kðx y; tÞ/ðx; 0Þ dy þX N j¼1 Z tDt 0 Kðx wjðsÞ; t sÞqjðsÞ ds þ Z tDt 0 oKðx wjðsÞ; t sÞ owj ljðsÞ ds " # þX N j¼1 Z t tDt Kðx wjðsÞ; t sÞqjðsÞ ds þ Z t tDt oKðx wjðsÞ; t sÞ owj ljðsÞ ds " # . ð19Þ
We further represent each step location wj(s) as a function of its velocity w0jðsÞ using Picard integral equation
wjðsÞ ¼ wjðt DtÞ þ
Z s
tDt
w0jðaÞ da. ð20Þ
The unknowns become the densities {qj(s)}, {lj(s)} and the step velocities w0jðsÞ for s from t Dt to t and
j = 1, . . . , N. We assume that when s2 [t Dt, t], approximate solutions to the unknowns are given by
q½0j ðsÞ, l ½0 j ðsÞ and ðw ½0 j ðsÞÞ 0
, the errors are then introduced as dqjðsÞ ¼ qjðsÞ q ½0 j ðsÞ; dljðsÞ ¼ ljðsÞ l ½0 j ðsÞ; dw0jðsÞ ¼ w0jðsÞ ðw ½0 j ðsÞÞ 0 . 8 > > < > > : ð21Þ
Plug these errors back to the original IDE system and the error equations follow.
For the low order method, consider one step in the marching scheme from ti1to tifor i = 1, . . . , p, where
t0= t Dt and t1, . . . , tpare the p Gaussian points in [t Dt, t]. Take the case i = 1 as an example (and for the
ease of notations), linear approximations to the error functions give dqjðsÞ ¼ dqjðt0Þ þ ðdqjðt1Þ dqjðt0ÞÞ st0 t1t0 ; dljðsÞ ¼ dljðt0Þ þ ðdljðt1Þ dljðt0ÞÞ st0 t1t0 ; dw0 jðsÞ ¼ dw0jðt0Þ þ ðdwj0ðt1Þ dw0jðt0ÞÞ st0 t1t0 . 8 > > > > < > > > > :
The unknowns are then reduced to dqjðt1Þ, dljðt1Þ, and dw0jðt1Þ. These are the error function values at t1which
can be derived using the following procedures:
The first step is to derive the locations of steps using an explicit method. Plug the approximation dw0jðgÞ dw0jðt0Þ into the Picard integral equation wjðsÞ ¼ w
½0 j ðsÞ þ Rs t0dw 0 jðgÞ dg, we have an approximation ~ wjðsÞ ¼ w½0j ðsÞ þ ðs t0Þdw0jðt0Þ ð22Þ
for s2 [t0, t1]. We want to mention again that w½0j ðsÞ is a polynomial of degree p 1. Next, as the locations of
steps are known, we solve the boundary conditions in Eqs.(15) and (16). For the unknown errors, letting x
approaching the jth step and assuming that t1 t0is small hence higher order terms of t1 t0are neglected,
these equations approximately take the form
1 2 w0j 2 dljðt1Þ þ dqjðt1Þ ! 1 2k þd ljþ o/ ox x¼~wjðt1Þ þ kþ/j x¼~wjðt1Þ¼ g1;jðt1Þ; ð23Þ 1 2 w0j 2 dljðt1Þ þ dqjðt1Þ ! 1 2k dljþ o/ ox x¼~wjðt1Þ k/jx¼~wjðt1Þ¼ g2;jðt1Þ; ð24Þ
where w0j¼ ðw½0j ðt0ÞÞ0þ dw0jðt0Þ and all available quantities are collected in g1,j(t1) and g2,j(t1). Derivation of these equations uses the tedious local Taylor expansion analysis and is performed using mathematica (same
procedure as in the proof of jump conditions for moving layer potentials). The first and third terms in Eq.(23)
come from the right limit of /xin Eq.(16)when x approaches ~wjðt1Þ. The first term is due to the
discontinu-ities in the derivatives of the single and double layer potentials when calculating the right limit of /x. The
sec-ond and fourth terms in Eq.(23)come from the right limit ofk+/in Eq.(16)when x approaches ~wjðt1Þ. The
second term is due to the discontinuity in the double layer potential when computing the right limit of /. Eq.
(24)is obtained similarly from Eq.(15).
In the final step, the velocity dw0j is updated according to
dw0jðt1Þ ¼ ðw½0j ðt1ÞÞ0 2ðq½0ðt 1Þ þ dqjðt1ÞÞ 2þ l½0ðt 1Þ þ dljðt1Þ ; ð25Þ
which is the error equation form of Eq.(14).
Recently in[39], it was shown that for linear problems, the original SDC method is equivalent to solving a
preconditioned linear system using Neumann series expansion. This observation was then coupled with New-ton–Krylov methods, and the performance of the accelerated scheme was dramatically improved, especially for stiff systems. Incorporating this new technique into our step-flow simulator is being proposed.
4.3. Efficient evaluation of diffusion potentials
In this section, we describe a fast convolution technique for the efficient evaluation of the history dependent layer potentials required by the algorithm presented in previous section. This technique was first introduced by
Greengard and Strain in[28]. The readers are referred to [26,28,29,44,67,68]for recent advancements. Since
this technique is not yet widely known, in the following, we describe the main ideas using the single layer potential as an example. It is given by
Z t
0 Z
CðsÞ
Kðx yðsÞ; t sÞqðyðsÞ; sÞ dsyds;
where K(x, t) is the periodic Green’s function, q(y(t), t) is a given density defined on the moving interface C(t). Because of the history dependency, if N points are used in discretizing C(t), J time steps are marched, and the potential is evaluated at M points for each step, then direct numerical evaluation of the potential requires
O(J2MN) operations. The Greengard–Strain technique, on the other hand, can reduce the work to the
asymp-totically optimal O(J(N + M)).
In [28], analytical results are presented in d dimensions. For simplicity of notations, in the following, we
focus on d = 1. In this case, the periodic Green’s function is given by Eq. (8) or (9), C(t) becomes the set
of N moving steps each denoted by wj(t), and the single layer potential becomes
Z t
0
XN
j¼1
Kðx wjðsÞ; t sÞqðwjðsÞ; sÞ ds.
By introducing a small constant D, usually chosen the same as Dt, the first step of the technique rewrites the layer potential as the sum of
Z tD 0 XN j¼1 Kðx wjðsÞ; t sÞqðwjðsÞ; sÞ ds and Z t tD XN j¼1 Kðx wjðsÞ; t sÞqðwjðsÞ; sÞ ds.
Following the terminology in [28], the first part is referred to as the history part, and the second the local
part. Notice that in the history part, as t s > D, the Fourier domain representation of K(x wj(s), t s)
decays rapidly as a function of jkj. Thus the history part can be represented by a rapidly decaying Fourier
Z tD 0 XN j¼1 Kðx wjðsÞ; t sÞqðwjðsÞ; sÞ ds Xq k¼q akðtÞeikx;
where q is a parameter depending only on the required accuracy and D, and the coefficients ak(t) are given by
akðtÞ ¼ Z tD 0 XN j¼1 1 2pe k2ðtsÞ eikwjðsÞqðw jðsÞ; sÞ ds; ð26Þ
which can be updated recursively in a marching scheme.
Assume that at time t the Fourier series expansionPqk¼qakðtÞeikxfor the history part is given. At time t + D,
Z tþD 0 XN j¼1 Kðx wjðsÞ; t þ D sÞqðwjðsÞ; sÞ ds ¼ Z t 0 XN j¼1 þ Z tþD t XN j¼1 X q k¼q akðt þ DÞeikxþ Z tþD t XN j¼1 Kðx wjðsÞ; t þ D sÞqðwjðsÞ; sÞ ds.
Therefore evaluation of the potential contains two parts: updating the Fourier coefficients ak(t + D) from ak(t),
and evaluating the local integral
Z tþD
t
XN
j¼1
Kðx wjðsÞ; t þ D sÞqðwjðsÞ; sÞ ds. ð27Þ
4.3.1. Updating the Fourier coefficients ak(t + D)
To calculate the Fourier coefficients ak(t + D), notice that
Xq k¼q akðt þ DÞeikx Z t 0 XN j¼1 Kðx wjðsÞ; t þ D sÞqðwjðsÞ; sÞ ds.
Divide this integral into two partsR0tPNj¼1¼R0tDPNj¼1þRtDt PNj¼1and use the Fourier domain representation
for K(x, t) in Eq.(8), the first part becomes
Z tD 0 XN j¼1 Kðx wjðsÞ; t þ D sÞqðwjðsÞ; sÞ ds ¼ Z tD 0 XN j¼1 X1 k¼1 1 2pe k2ðtþDsÞ eikðxwjðsÞÞqðw jðsÞ; sÞ ds;
which, by comparing with ak(t) in (26), is approximately
Xq
k¼q ek2Da
kðtÞeikx.
For the second part we have Z t tD XN j¼1 Kðx wjðsÞ; t þ D sÞqðwjðsÞ; sÞ ds Xq k¼q bkðt þ DÞeikx; where bkðt þ DÞ ¼ Z t tD XN j¼1 1 2pe k2ðtþDsÞ eikwjðsÞqðw jðsÞ; sÞ ds ð28Þ
can be calculated using numerical quadratures. Notice that when we have N steps (interfaces), the total
num-ber of operations for calculating bk(t + D) for k =q, . . . , q is O(qN), hence the total amount of work to derive
the coefficients {ak(t + D)} is O(qN) using
akðt þ DÞ ¼ ek
2D
The evaluation of Fourier series at M points requires O(qM) operations, therefore, evaluating the history part
for a total of J time steps requires O(qJ(N + M)) operations instead of O(J2NM).
4.3.2. Evaluating the local part
For the local part in(27), traditional trapezoidal rule or Gaussian quadratures in time could not provide an
accurate solution due to the singularity of the Green’s function. In our solver, we apply a slightly modified product rule, where the singularity of the Green’s function is extracted, and the remaining smooth part in the integrand is approximated by a polynomial. This approximation is then analytically evaluated, as explained by the following procedure.
Consider the simplified moving layer potential in one dimensional space,
Z t 0 eðxyðsÞÞ2s ffiffiffi s p qðyðsÞ; sÞ ds; ð30Þ
notice that y(s) = y(0) + sh(s) where h(s) is a smooth function, the singular part of the integrand can be
rep-resented as eðxyð0ÞÞ2s =pffiffiffis, and the remaining part becomes smooth and can be approximated by a Legendre
polynomial expansion e2ðxyð0ÞÞhðsÞsh2ðsÞqðyðsÞ; sÞ X P k¼0 ckLkðsÞ. Using formula Z t 0 ex2 s ffiffiffi s p spds¼ x2pþ1Gamma 1 2 p; x2 t ;
where Gamma(a, z) is the incomplete Gamma function defined by
Z 1
z
ta1et;
the local part can then be evaluated to higher order.
In the current implementation, the local part is evaluated directly using above techniques. When the num-ber of steps N is fewer than several thousands, we found the computation efficiency is very acceptable. How-ever, when N and M are large, evaluating the local part at each time step requires approximately O(NM) work, which will be the dominating cost in our solver. Reducing this cost is currently being studied and it was found
that techniques similar to the new version of fast Gauss transform[30]can be used to reduce the amount of
work to O(N + M). Discussion of this accelerated method will be reported at a later time. 5. Local collision model
In the numerical simulations, we noticed that for certain initial and parameter settings, there is a possibility that two steps collide and hence form a ‘‘singularity’’. Such collisions have been observed for steps with the
same as well as opposite height signs as shown inFig. 2. In the left plot, the convergence test is carried out for
collisions of two steps with heights +1 and1, respectively. In the right plot, both of the colliding steps have
height +1. Convergence tests show that the collisions are very unlikely the results of numerical truncation error. Analytical understanding of these ‘‘singularities’’ and conditions under which they may happen are still research topics being investigated.
For most front tracking methods, when the steps are close to collision, numerical simulations may require extremely fine mesh in both temporal and spatial directions to be employed near the singularity. This kind of singularity or near singularity formations have been one of the major numerical difficulties, especially for inte-gral equation based schemes.
In our algorithm, instead of adaptive mesh refinement, we introduce the ‘‘local collision model’’: when the steps are ‘‘close’’ to each other, we assume they march at constant velocity; and once they collide, we adopt a
simple physical model where two steps either cancel each other when they have opposite signs (see the first
image inFig. 2), or stick to each other and move as a new step when they have the same sign (see the second
image inFig. 2). The ‘‘constant speed’’ assumption is mostly based on our refined time-step numerical
simu-lations. InFig. 2, it can be seen that when two steps are close to collision, the leading order approximation of
the step speed is given by a constant. The cancelation of two steps with opposite signs after the collision is physical, it means that the gap between two terraces has been filled by adatoms and the two terraces are
con-nected to form a new one (see the first image inFig. 2). It is also physical that two steps with same sign move
together after they get close, because there is a strong short-range repulsive interaction between such two steps
[14,69,58,72,73], the step above cannot surpass the one below.
When p (P 2) steps move together, we adopt the method used by Elkinani and Villain[17]. The speed of a
step bunch consisting of p steps wj, wj+1, . . . , wj+p1, when all steps have height +1, is v¼ ðk/jx¼w
jþ
kþ/jx¼wþ
jþp1Þ=p. The step bunch may dissociate when step wj+p1 tends to move faster than the other steps,
i.e., when k/jx¼w
j=ðp 1Þ < k
þ/j x¼wþ
jþp1. Also, if more accurate local properties need to be resolved, the full
short-range repulsive interactions can be added to our model. 6. Numerical examples
In this section, we present several simulation results for epitaxial step-flow growth in 1 + 1 dimensions. We compare our numerical results with available analytical and numerical solutions for motions of steps in slow
deposition regime (Pe 1) and motion of a uniform step train in fast deposition regime (Pe 1).
When a uniform step train moves in steady state, there is an analytical solution to the BCF model(12). In
our first numerical experiment, we consider a uniform step train and compare our simulation results with ana-lytical solutions. We use n = 4 steps and periodic boundary conditions in [p, p]. The length of each terrace is
uniformly lj¼ l ¼2pn 1:57. We set k+= k= 1, F ¼ð1el2Þð2þlÞ. The Peclet number is hence approximately
1.7. Assume initially /(x, 0) is the steady state solution to system(12), then it is straight forward to verify that
the solution w0j¼ 1; /ðx; tÞ ¼ 2 1el F þ exþt el1 F ðx tÞ ( ð31Þ
solves the BCF model. InFig. 3, we study numerical error in the location of first step w1(t) as a function of
time. The analytical solution is given by a linear function. In the plot, q is the number of terms in the Fourier expansion for the history part. It can be seen that approximately 110 Fourier terms can guarantee seven digits
for Dt = 103. Same conclusion applies to other variables as well.
3.196 3.1965 3.197 3.1975 3.198 0.1 0.11 0.12 0.13 0.14 0.15 0.16 step locations time
collision convergence test case I
dt=0.05 dt=0.025 dt=0.0125 dt=0.00625 5.20031 5.2004 5.2005 5.2006 5.2007 1.05 1.1 1.15 1.2 1.25 1.3 1.35 1.4 1.45 step locations time case II dt=0.005 dt=0.00025 dt=0.000125
It is well known that the Schwoebel barrier (k+> k) stabilizes uniform step trains[64]. Our next example
considers this stabilizing effect. We consider a step train with n = 50 steps, each with height1. The average
terrace length is given by l¼ 2p=n 0:13. For the Schwoebel barrier, we set k+
= 109and k= 1. We set
F = 0.006 and hence the Peclet number is approximately 104. The initial location of each step is set to
wj¼ 2pj=n þ Rjl where Rjis a random variable uniformly distributed in [0.005, 0.005]. Due to the Schwoebel
effect, the initially perturbed step train quickly converges to a uniform step train, and the velocity of each step
converges to w0j¼ F l 0:00075. Our numerical simulation results agree with these analytical results, as shown
in Fig. 4. In the left, we plot the variance of step locations defined as 1 n 1 Xn j¼1 ðwjþ1 wj lÞ 2 ;
where wn+1is defined as w1+ 2p. In the right, we plot the velocity of the first step as a function of time.
In our third example, we numerically verify the step velocity oscillations in a fast step-flow observed in
pre-vious simulations[53,70]. The result is shown inFig. 5. As in[70], we start the simulation from a uniform step
train, with parameters k+= 109, k= 102, F = 0.6, and l¼ 2p, hence the Peclet number is approximately 24.
The initial / is set to zero, representing the growth from a clean surface. It can be seen that in the simulation, oscillations occur in step velocity while the solution converges to the steady state. This simulation result agrees
with those in[53,70](see, for example, Fig. 6.25 in[70]). See Section2for a brief review and discussion on this
step velocity oscillation.
Applications of our method to problems with stationary interfaces are straight forward. To further test the correctness of our code, we study a problem with fixed interfaces
0 0.2 0.4 0.6 0.8 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 q=20 time error 0 0.2 0.4 0.6 0.8 0 0.5 1 1.5 2 2.5 3x 10 –3 q=50 time error 0 0.2 0.4 0.6 0.8 0 0.5 1 1.5 2 2.5 3 3.5 x 10 –5 q=80 time error 0 0.2 0.4 0.6 0.8 –2 0 2 4 6 8 10 12 x 10 –8 q=110 time error
ut¼ uxxþ 1; ux u ¼ t 2et; when x¼ 0; uxþ u ¼ t; when x¼ L; w0¼ 0. 8 > > > < > > > :
Note that this example is used only to test the correctness of our code for stationary interfaces. The system is
slightly different from BCF system(12)and not necessarily a model for step-flow growth. The analytical
solu-tion to this problem is given by u(x, y) = ex+t+ t and the double layer density is given by l(t) = et(1 eL),
where L is the length of each periodic interval and we assume the solution is periodically expanded. InFig. 6,
we plot the error of the computed double layer density function. In the calculation, we set L = 2p/4 and
Dt= 103. For the first 1000 time steps, the numerical results are accurate in the first 6 digits.
To test the order of SDC accelerated solver, in the next example, we present convergence results for our method with different number of SDC corrections for a simple problem with 4 steps. The initial locations
of steps are not uniform, hence analytical solution is not available. InTable 1, we present convergence study
of our method using two Gaussian–Radau points (t2is the right end point). M1 is the solution derived using
the scheme introduced in Section4.1, but with only one correction, and M2 is the method with two
correc-tions. The convergence rate is determined by the ratio xDtxDt=2
xDt=2xDt=4, where Dt is the step size in the marching
0 0.5 1 1.5 2 2.5 3 3.5 0 2 4 6 8 10 12 14 16 18 20 Velocity Oscillation time v elocity
Fig. 5. Oscillatory flow of an equi-spaced array of ‘‘fast’’ steps.
0 0.01 0.02 0.03 0.04 0.05 0 1 2 3 4 5 6 7 8 x 10 -4 time v elocity
Velocity of first step
0 500 1000 1500 2000 0 1 2 3 4 5 6 x 10 -4 time v a riance Location variance
scheme, and xDtis the location of the third step at time t = 1. It can be shown from the numerical results that the method with two corrections is second order. However, we want to mention that for the same time step size, efficiency of the second order solver slows down significantly. This is due to the iterative refining proce-dure and evaluation of special functions. Currently, we are trying to improve the efficiency by coupling the
GMRES accelerated SDC methods introduced in[39]with precomputed tables, and by choosing the optimal
step size and linear equation solver for a prescribed accuracy requirement.
Finally, we perform simulations for two systems of non-uniform arrays of steps in the regimes Pe = O(1)
and Pe 1, respectively. InFig. 7, when the Pe number is approximately 0.4, we consider 10 steps which are
non-uniformly distributed initially and the heights of steps are set to +1. We use k+= 107and k= 10. In the
left, the locations of steps are plotted as a function of time, and in the right, the relative locations (with respect to the first step) are shown. From the plots, oscillations of the step locations can be observed, which we refer to as the ‘‘water wave’’ phenomenon.
InFig. 8, we consider a surface which initially has 100 steps and the heights are randomly chosen as +1 or1.
The Peclet number is approximately 0.004 for this problem and we choose k+= 109and k= 467. In the local
collision model, we assume that two steps either cancel each other if they have opposite heights, or they move together as a new step when the sum of heights is nonzero. The cutoff distance is the lattice constant a = 0.0012.
Also, nucleation of new steps on local maximal terraces is allowed and the critical nucleation length is lc= 0.063.
The Schwoebel effect is weak, indicated by a small ls: ls= D/k a = 0.0009 < a. We plot several snapshots of
the surface inFig. 8for t = 0, 2000, 4000, 6000, and 8000, respectively. In this simulation, both collisions and
nucleations are frequently encountered. The number of steps increases from 100 to approximately 200 during the simulation process (hence the Peclet number becomes smaller) and the ‘‘coarsening process’’ (big mounds are formed by combining smaller ones) can be observed.
0 0.5 1 1.5 2 2.5 3 0 1 2 3 4
5x 10 Double Layer Density Error
time
Error
Fig. 6. Double layer density function error.
Table 1
A second order method based on SDC
Dt D= 0.01 D/2 D/4 D/8 D/16 D/32
xDt(M1) 5.667252 5.668801 5.669575 5.669962 5.670156 5.670252
Ratio 2.00 2.00 1.99 2.02
xDt(M2) 5.670149 5.670299 5.670336 5.670346 5.670348 5.670349
Also in this simulation, when t = 2000, the surface profile shows a ‘‘selected slope’’ (see left ofFig. 9), which is a result of the deterministic nucleation on the top terraces and the weak Schwoebel effect. This can be ana-lyzed using the quasi-static approximation as follows: Consider a top terrace with a uniform step train on its
right, see for example, the left peak in the left ofFig. 9. Assume the distance between steps in the uniform step
train is l, and the width of the top terrace is l0. Since the Schwoebel effect is weak, it can be neglected when
calculating the velocity of the downward step wjof the top terrace right above the uniform step train, which
from Eq.(5) is given by
0 10 20 30 0 5 10 15 20 25 30 35 40 x(j,t) time Step Locations 0 2 4 6 8 0 5 10 15 20 25 30 35 40 time
Relative Step Locations
Fig. 7. Step location oscillations.
0 1 2 3 4 5 6 7
step locations
step heights
Fig. 8. Surface coarsening phenomenon.
2.8 2.9 3 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8
t=2000
3.8 4 4.2 4.4 4.6 4.8 5 5.2 5.4
t=8000
dwj=dt¼ 1 2dl0=dt¼ 1 2a 2Fl 0þ 1 2a 2 Fl.
In a steady state, the width of the top terrace l0increases approximately from 0 to lc between two
nucle-ation events in time t = 1/(a2F). If we assume that the uniform step train has a fixed slope, i.e., l is fixed,
we can solve the ODE and get l = lc/(e 1) = 0.58lc, and the selected slope is a/l = 0.033. The role of the
Schwoebel effect here is to stabilize the uniform step train. The left image in Fig. 9 shows the comparison
of the uniform step train at t = 2000 and the continuum surfaces with the selected slope. They agree quite well.
The surfaces with selected slope are not stable according to our simulations: new patterns form after long time. There are some continuum approximations for the stationary states (derived from quasi-static approx-imation) which model the competition between the Schwoebel and up–down asymmetry effects of the stepped
surface [60,72,73]. The analytical expression for the stationary profile (right branch) is given by
hðxÞ ¼ a 2 3ls ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ls a 1 mð0Þ 2 lsx 6a2 s þls a log ls a ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ls a 1 mð0Þ 2 lsx 6a2 s ls aþ 1 mð0Þþ ls a logjmð0Þj 0 @ 1 A þ hð0Þ;
where the peak is at x = 0, h(0) is the height of the peak, and m(0) is the slope on one side of the peak. From
the right ofFig. 9, some portions of our simulation result at t = 8000 do show this stationary profile, in which
m(0) is taken to be a/l1, where l1is the width of the terrace right below the top one. As discussed in[17,60], this
stationary profile does not apply to the valleys due to the deterministic nucleation and the weak Schwoebel effect: deep valleys tend to be healed as time increases based on the quasi-static approximation of the step dynamics; while the stationary profile gives an infinitely deep crack. This is also seen in our simulation results. The standard deviation of the surface width measures the roughness of the surface, and usually increases with
time according to a power law tb[58]. In this simulation, we found that b = 0.43. The exponent b 0.5 has been
seen in some experiments [18,19], kinetic Monte Carlo simulations [1,66,74], and analysis using continuum
models without slope selection[31,45].
7. Conclusion and future work
In this paper, we discuss an integral equation approach for solving diffusion problems with moving inter-faces. Jump conditions of moving layer potentials are studied and the results are applied to the BCF model which describes the epitaxial step-flow growth of crystals. The partial differential equation is reduced to an integro-differential equation system where the unknowns are only defined on the steps. For this system, we introduce two recently developed numerical techniques: higher order can be derived using spectral deferred correction ideas in the temporal direction; and moving layer potentials can be efficiently evaluated using fast convolution methods designed for dense convolution type matrices. The resulting simulation toolbox has been tested and applied to the BCF model in 1 + 1 dimensions with different parameter settings. Numerical results are compared with available analytical and simulation results.
Recently, several new techniques have been proposed to further improve the accuracy and efficiency of the current solver. These include the acceleration of the SDC techniques using Krylov subspace methods and the efficient evaluation of the ‘‘local part’’ in diffusion potentials of both volume and layer types.
We have focused on the integral equation method in this paper. Applications of the method in 1 + 1 dimensions include the study and validation of the quasi-static approximation under various circumstances, and when it is no longer valid, the properties of the adatom density and related nucleation models. Real
crystal surfaces are 2 + 1 dimensional, and straight steps are unstable under the Schwoebel effect[3,56–59].
Generalization of the solver to 2 + 1 dimensions is being considered. Except for several technical details including the discretization of steps and local collision and nucleation models, such generalization is
straightforward. More physics such as the step–step and adatom–step elastic interactions [14,69,72,73]
and more realistic nucleation models can be incorporated within the framework of our integral equation method. These generalizations and applications are being proposed and results along these directions will be reported in the future.
Acknowledgements
We would like to express our gratitude to Prof. Weinan E of Princeton University for helpful discussions concerning the epitaxial step-flow growth models. We would also like to thank Prof. Leslie Greengard of Cou-rant Institute at New York University for many interesting discussions concerning fast solvers for the heat equation. The work of JH was supported by a UNC Junior Faculty Development Award, and by the National Center of Theoretical Sciences of Taiwan during his visit. The work of YX was supported by Hong Kong RGC Competitive Earmarked Research Grant 604604, and the work of ML was supported by NSC of Taiwan under Research Grant NSC-93-2113-M-009-008.
References
[1] J.G. Amar, F. Family, Effects of crystalline microstructure on epitaxial growth, Phys. Rev. B 54 (1996) 14742–14753. [2] A.W. Appel, An efficient program for many-body simulations, SIAM J. Sci. Stat. Comput. 6 (1985) 85–103.
[3] G.S. Bales, A. Zangwill, Morphological instability of a terrace edge during step-flow growth, Phys. Rev. B 41 (1990) 5500–5508. [4] E. Bansch, F. Hauser, O. Lakkis, B. Li, A. Voigt, Finite element method for epitaxial growth with attachment-detachment kinetics, J.
Comput. Phys. 194 (2004) 409–434.
[5] J. Barnes, P. Hut, A hierarchical O(n log n) force calculation algorithm, Nature 324 (1986) 446–449.
[6] K. Bohmer, P. Hemker, H.J. Stetter, The defect correction approach, in: K. Bohmer, H.J. Stetter (Eds.), Defect Correction Methods. Theory and Applications, Springer-Verlag, 1984, pp. 1–32.
[7] A. Bourlioux, A.T. Layton, M.L. Minion, High-order multi-implicit spectral deferred correction methods for problems of reactive flow, J. Comput. Phys. 189 (2003) 351–376.
[8] K.E. Brenan, S.L. Campbell, L.R. Petzold, Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations, SIAM, Philadelphia, 1995.
[9] W.K. Burton, N. Cabrera, F.C. Frank, The growth of crystals and the equilibrium structure of their surfaces, Philos. Trans. Roy. Soc. London 243A (1951) 299–358.
[10] J.W. Cooley, J.W. Tukey, An algorithm for the machine calculation of complex Fourier series, Math. Comput. 19 (1965) 297–301. [11] R.E. Caflisch, M.F. Gyure, B. Merriman, S. Osher, C. Ratsch, D.D. Vvedensky, J.J. Zinck, Island dynamics and the level set method
for epitaxial growth, Appl. Math. Lett. 12 (1999) 13–22.
[12] S. Chen, B. Merriman, M. Kang, R.E. Caflisch, C. Ratsch, L.T. Cheng, M. Gyure, R.P. Fedkiw, C. Anderson, S. Osher, A level set method for thin film epitaxial growth, J. Comput. Phys. 167 (2001) 475–500.
[13] T. Darden, D. York, L. Pedersen, Particle mesh Ewald: an N log(N) method for Ewald sums in large systems, J. Chem. Phys. 98 (1993) 10089–10092.
[14] C. Duport, P. Politi, J. Villain, Growth instabilities induced by elasticity in a vicinal surface, J. Phys. I 5 (1995) 1317–1350. [15] A. Dutt, L. Greengard, V. Rokhlin, Spectral deferred correction methods for ordinary differential equations, BIT 40 (2) (2000) 241–
266.
[16] H.P. Dym, H.P. McKean, Fourier Series and Integrals, Academic Press, San Diego, 1972.
[17] I. Elkinani, J. Villain, Growth roughness and instabilities due to the Schwoebel effect: a one-dimensional model, J. Phys. I 4 (1994) 949–973.
[18] W.C. Elliott, P.F. Miceli, T. Tse, P.W. Stephens, Temperature and orientation dependence of kinetic roughening during homoepitaxy: a quantitative X-ray-scattering study of Ag, Phys. Rev. B 54 (1996) 17938–17942.
[19] H.J. Ernst, F. Fabre, R. Folkerts, J. Lapujoulade, Observation of a growth instability during low temperature molecular beam epitaxy, Phys. Rev. Lett. 72 (1994) 112–115.
[20] L. Fox, Some improvements in the use of relaxation methods for the solution of ordinary and partial differential equations, Proc. Roy. Soc. London A 190 (1020) (1947) 31–59.
[21] R. Frank, C.W. Ueberhuber, Iterated defect correction for the efficient solution of stiff systems of ordinary differential equations, BIT 17 (1977) 146–159.
[22] R. Ghez, S.S. Iyer, The kinetics of fast steps on crystal surfaces and its application to the molecular beam epitaxy of silicon, IBM J. Res. Develop. 32 (1988) 804–818.
[23] R. Ghez, H.G. Cohen, J.B. Keller, The stability of growing or evaporating crystals, J. Appl. Phys. 73 (1993) 3685–3693. [24] L. Greengard, The Rapid Evaluation of Potential Fields in Particle Systems, MIT Press, Cambridge, MA, 1988. [25] L. Greengard, Spectral integration and two-point boundary value problems, SIAM J. Numer. Anal. 28 (1991) 1071–1080. [26] L. Greengard, P. Lin, Spectral approximation of the free-space heat kernel, Appl. Comput. Harmonic Anal. 9 (1) (2000) 83–97. [27] L. Greengard, V. Rokhlin, A fast algorithm for particle simulations, J. Comput. Phys. 73 (1987) 325–348.
[28] L. Greengard, J. Strain, A fast algorithm for the evaluation of heat potentials, Commun. Pure Appl. Math. 43 (1990) 949–963. [29] L. Greengard, J. Strain, The fast Gauss transform, SIAM J. Sci. Stat. Comput. 12 (1991) 79–94.
[30] L. Greengard, X. Sun, A New Version of the Fast Gauss Transform, Documenta Mathematica, Extra Volume ICM, III, 1998, pp. 575–584.
[31] L. Golubovic, Interfacial coarsening in epitaxial growth models without slope selection, Phys. Rev. Lett. 78 (1997) 90–93. [32] R.B. Guenther, J.W. Lee, Partial Differential Equations of Mathematical Physics and Integral Equations, Prentice-Hall, 1988.
[33] B. Gustafsson, W. Kress, Deferred correction methods for initial value problems, BIT 41 (2001) 986–995.
[34] B. Gustafsson, L. Hemmingsson-Franden, Deferred correction in space and time, J. Sci. Comput. 17 (1–4) (2002) 541–550. [35] E. Hairer, S.P. Norsett, G. Wanner, Solving Ordinary Differential Equations I, Non-Stiff Problems, Springer-Verlag, Berlin, 1987. [36] E. Hairer, G. Wanner, Solving Ordinary Differential Equations II, Springer, 1996.
[37] S. Harris, Onset of fast step-velocity oscillations during growth by molecular-beam epitaxy, Phys. Rev. B 51 (1995) 4415–4417. [38] R. Hockney, J. Eastwood, Computer Simulation Using Particles, McGraw-Hill, New York, 1981.
[39] J. Huang, J. Jia, M. Minion, Accelerating the convergence of spectral deferred correction methods, J. Comput. Phys., in press (doi:10.1016/j.jcp.2005.10.004).
[40] W. Kress, B. Gustafsson, Deferred correction methods for initial boundary value problems, J. Sci. Comput. 17 (1–4) (2002) 241–251. [41] A. Karma, M. Plapp, Spiral surface growth without desorption, Phys. Rev. Lett. 81 (1998) 4444–4447.
[42] J. Krug, On the shape of wedding cakes, J. Stat. Phys. 87 (1997) 505–518.
[43] A.T. Layton, M.L. Minion, Conservative multi-implicit spectral deferred correction methods for reacting gas dynamics, J. Comput. Phys 194 (2) (2004) 697–714.
[44] P. Lin, On the numerical solution of the heat equation in unbounded domains, Ph.D. Thesis, New York University, 1993. [45] B. Li, J.G. Liu, Epitaxial growth without slope selection: energetics, coarsening, and dynamics scaling, J. Nonlinear Sci. 14 (2004)
429–451.
[46] F. Liu, H. Metiu, Stability and kinetics of step motion on crystal surfaces, Phys. Rev. E 49 (1994) 2601–2616.
[47] M.L. Minion, Semi-implicit spectral deferred correction methods for ordinary differential equations, Commun. Math. Sci. 1 (3) (2003) 471–500.
[48] A.K. Myers-Beaghton, D.D. Vvedensky, Nonlinear model for temporal evolution of stepped surfaces during molecular-beam epitaxy, Phys. Rev. B 42 (1990) 9720–9723.
[49] S. Osher, J.A. Sethian, Fronts propagating with curvature-dependent speed: algorithms based on Hamilton–Jacobi formulations, J. Comput. Phys. 79 (1) (1988) 12–49.
[50] F. Otto, P. Penzler, A. Ratz, T. Rump, A. Voigt, A diffuse-interface approximation for step flow in epitaxial growth, Nonlinearity 17 (2004) 477–491.
[51] V. Pereyra, On improving an approximate solution of a functional equation by deferred corrections, Numer. Math. 8 (1966) 376–391. [52] V. Pereyra, Iterated deferred corrections for nonlinear operator equations, Numer. Math. 10 (1966) 316–323.
[53] G.S. Petrich, P.R. Pukite, A.M. Wowchak, G.J. Whaley, P.I. Cohen, A.S. Arrott, On the origin of RHEED intensity oscillations, J. Cryst. Growth 95 (1989) 23–27.
[54] J.R. Phillips, J. White, A precorrected-FFT method for capacitance extraction of complicated 3-D structures, in: Proceedings of ICCAD-94, 1994, pp. 268–271.
[55] O. Pierre-Louis, Phase field models for step flow, Phys. Rev. B 68 (2003) 021604. [56] O. Pierre-Louis, Dynamics of crystal steps, C. R. Phys. 6 (2005) 11–21.
[57] O. Pierre-Louis, C. Misbah, Y. Saito, J. Krug, P. Politi, New nonlinear evolution equation for steps during molecular beam epitaxy on vicinal surfaces, Phys. Rev. Lett. 80 (1998) 4221–4224.
[58] A. Pimpinelli, J. Villain, Physics of Crystal Growth, Cambridge University Press, 1998.
[59] P. Politi, G. Grenet, A. Marty, A. Ponchet, J. Villain, Instabilities in crystal growth by atomic or molecular beams, Phys. Rep. 324 (5– 6) (2000) 271–404.
[60] P. Politi, J. Villain, Ehrlich-Schwoebel instability in molecular-beam epitaxy: a minimal model, Phys. Rev. B 54 (1996) 5114–5129. [61] A. Rangan, Adaptive solvers for partial differential and differential-algebraic equations, Ph.D. Thesis, University of California at
Berkeley, 2003.
[62] A. Rangan, Deferred correction methods for low index differential algebraic equations, Preprint.
[63] A. Ratz, A. Voigt, Various phase-field approximations for epitaxial growth, J. Cryst. Growth 266 (2004) 278–282. [64] R.L. Schwoebel, Step motion on crystal surfaces, J. Appl. Phys. 40 (1969) 614–619.
[65] J.A. Sethian, J. Strain, Crystal growth and dendritic solidification, J. Comput. Phys. 98 (1992) 231–253.
[66] P. Smilauer, D.D. Vvedensky, Coarsening and slope evolution during unstable epitaxial growth, Phys. Rev. B 52 (1995) 14263–14272. [67] J. Strain, Fast potential theory II. Layer potentials and discrete sums, J. Comput. Phys. 99 (1992) 251–270.
[68] J. Strain, Fast adaptive methods for the free-space heat equation, SIAM J. Sci. Comput. 15 (1992) 185–206.
[69] J. Tersoff, Y.H. Phang, Z. Zhang, M.G. Lagally, Step-bunching instability of vicinal surfaces under stress, Phys. Rev. Lett. 75 (1995) 2730–2733.
[70] J.Y. Tsao, Fundamentals of Molecular Beam Epitaxy, Academic, San Diego, 1993.
[71] K. Voigtlaender, H. Risken, E. Kasper, Modified growth theory for high supersaturation, Appl. Phys. A 39 (1986) 31–36. [72] Y. Xiang, Derivation of a continuum model for epitaxial growth with elasticity, SIAM J. Appl. Math. 63 (1) (2002) 241–258. [73] Y. Xiang, W. E, Misfit elastic energy and a continuum model for epitaxial growth with elasticity, Phys. Rev. B 69 (2004) 035409. [74] Z. Zhang, J. Detch, H. Metiu, Surface roughness in thin-film growth: the effect of mass transport between layers, Phys. Rev. B 48
