Simulation of Kinetic Model for Island Dynamics

1
Simulation of Kinetic Model for Island Dynamics

• Lovell David Shao
• Materials group meeting - Oct. 27th, 2004

2
Outline

• I. Kinetic model for a step edge
• II. Implementation issues
• III. Numerical schemes
• IV. Some results
• V. Work in progress

3
I. Kinetic model for a step edge

• Based on the paper by
• R.E. Caflisch, W. E, M.F. Gyure, B. Merriman,
  and C. Ratsch, Kinetic model for a step edge in
  epitaxial growth, Phys. Rev. E, 59 (1999), pp.
  6879-6887.
• This model determines the velocity of a step edge
  based strictly on kinetics with no
  near-equilibrium assumptions. It agrees with the
  BCF model when the system is in equilibrium.
• This model differs from our current island
  dynamics model in that it includes edge adatoms
  and kinks, in addition to adatoms on the terrace.

4

• ? density of adatoms on the terrace
• ? density of edge adatoms along the step edge
  (have one nearest neighbor)
• k density of kinks along the step edge (have
  two nearest neighbors) k kr kl where kr is
  the kink that moves to the right kl is the one
  that moves to the left.

5

• Modeling equations
• Adatom diffusion on the terrace
• BCs
• F deposition flux
• M nucleation rate
• desorption rate, adatom
  diffusion coeff on the terrace
• v normal velocity
• n unit normal vector
• f net flux of adatoms to the step edge from
  upper terrace
• f- net flux of adatoms to the step edge from
  lower terrace

6

• (2) Diffusion of edge adatoms along the step edge
• edge diffusion coefficient
• net flux of edge adatoms to kinks
• s arclength variable along the step edge
• (3) Convection of kinks along the step edge
• w kink velocity
• g net gain in kink pairs due to
  nucleation/breakup
• h net loss in kink pairs due to
  creation/annilihation

7

• Constitutive laws (a lattice constant)
• (4)
• (5) (? angle between crystallographic
  axis n)
• (6) (derived from conservation of mass
  ? curvature)
• Normal velocity
• (7)
• are formulated in terms of diffusion
• coefficients and coordination numbers using
  kinetic mean
• -field approximation. They describe the dynamics
  at the
• microscopic level and are used as parameters in
  the
• macroscopic model.

8

• Validation of kinetic theory
• (for the SOS model)
• At equilibrium, the simulation of the kinetic
  model using KMC is consistent with BCF theory.
• At steady state, the difference between the
  theoretical and computational results (KMC) for
  the kink density k is very small.

9
II. Implementation issues

• For a circular island on a cartesian grid, the
  PDEs to be solved are
• (1) in ? (? ? U ?-)
• on ? (from the ? side)
• on ?- (from the ?- side)
• Periodic BCs on the boundaries of the cartesian
  grid
• (2) ( ) (
  ) on ?
• (3) on ?,
• (4) , vn kwcos(?)
• R.E. Caflisch and B. Li, Analysis of island
  dynamics in epitaxial growth of thin films, SIAM
  Mult. Model. Simul., 1 (2003), pp. 150-171.

10

• Solving the diffusion equation with mixed
  boundary condition (requires the normal
  derivative) on an irregular domain is difficult.
• In addition, since Dt is of the order
  to , we have to use implicit schemes for
  time stepping, or else dt is too small (requires
  too much work).
• To solve the edge diffusion and kink convection
  along ? on a cartesian grid using only the
  levelset function is not so straightforward. And
  we need to use implicit time stepping again b/c
  De is also large.

11

• Its not so clear what will happen when ? p/4
  (? is the angle between the unit normal vector of
  the island boundary and one of the principal
  crystallographic axes).
• The kinetic model for a step edge might not
  handle this case so well we might need to add
  more physics to the current kinetic model so that
  it incorporates this process on the corners.

12
III. Numerical Schemes

• For adatom diffusion on the terrace, we
  use
• 2nd order centered differencing for the Laplacian
  operator ? (interior points)
• Backward Euler (implicit method) for the
  timestepping
• note CN isnt used b/c oscillation occurs if
  timestep is large.
• So for the interior points, we have
• But what about the boundary points (points that
  are close to the boundary
• the boundary is within a distance of h)? How do
  we solve for there?

13

• Its hard to incorporate the mixed BC into the
  5-pt stencil for diffusion we cant use the GFM
  (Ghost Fluid Method) because it assumes that we
  have Dirichlet BC on the boundary.
• There are some works on solving the Poisson
  equation on irregular domains with jump condition
  using similar ideas as GFM (by X.D. Liu, R.
  Fedkiw, and M. Kang), but for our problem, we
  dont have a jump condition
• . Our problem is basically solving the
  diffusion problem in
• separate domains (? and ? -) with
  different mixed BCs.

14

• There is also IIM (Immersed Interface Methods)
  developed by LeVeque and Li which can be applied
  to mixed BCs (robin) on an irregular domain (A.
  Fogelson and J. Keener) with 2nd order accuracy.
  However, the scheme is very complicated and hard
  to implement because it requires the careful
  selection of the stencil used and the
  coefficients used in the stencil.
• Since, for a start, we can live with a 1st order
  accurate method, what we did
• to solve for the boundary points is
• This is like perturbing the boundary by O(h)
  since we basically replaced the diffusion
  equation at the boundary points with the equation
  of the mixed BC. And because of this, our method
  solving the diffusion equation can only be at
  most 1st order accurate.

15

• Its ok (for now) that we dont need 2nd order
  accuracy for ? because the normal velocity is not
  depended on the normal derivative of ? but on
  k,w,?.
• We have a non-symmetric matrix A? b to solve,
  so GMRES is used.
• The rate of convergence is okay, like 28
  iterations for grid size of
• 128 x 128, but can be improved (using
  preconditioned GMRES or may be some other
  iterative methods for non-symmetric matrix).

16

• For the diffusion of edge adatoms along the
  boundary ?, we use a semi-
• implicit scheme using the levelset function
  developed by J. Xu and H.
• Zhao. The edge diffusion equation we need to
  solve is
• Note that we need to use Backward Euler for the
  timestep because De is large
• timestep will be very small if explicit
  scheme is used!
• Its difficult to implement the surface diffusion
  operator on a cartesian grid, so what Xu
  Zhao did is to decompose
• and only let be applied at time n1 (thats
  why its only semi-implicit).

17

• Because we have the level set function u, so we
  can determine the unit normal vector n and the
  curvature ? on the island boundary.
• Since we are on a cartesian grid, but we want to
  diffuse the edge adatom density ? along the
  island boundary, what we need to do is to extend
  ? constant in the normal direction to the island
  boundary. Then we can solve for ? on the whole
  cartesian grid.
• So the final equation we have to solve is
• The matrix (discrete operator) on the left is
  symmetric, so Conjugate
• Gradient (CG) is used to get rapid convergence.
  This method is 2nd order
• in space, 1st order in time.

18

• For the convection of kinks along the boundary ?,
• We rewrite the equation as
• where since we have the level
  set function u.
• Now we dont have kr and kl, but using the
  constitutive laws
• we get kr (k tan(?))/2
• kl (k tan(?))/2

19

• Since the methods to solve the convection
  equation is well established,
• we will use
• Upwind scheme for the convective term (1st order
  one-sided diff. or WENO3)
• TVD-RK3 for the explicit time stepping.
• The source/sink terms g,h are depended on ?, ?
  -, and ? at time n so they
• can be evaluated without any problem.
• And again, since the kink density k is only on
  the island boundary, we
• need to extend k constant in the normal direction
  from the boundary so that
• we can compute on the cartesian grid.

20

• Numerical Algorithm
• at time n,
• Extend ?,k constant in the normal direction
• Solve the edge diffusion equation for ?
• Solve the edge convection equation for k
• Update the level set function u to time n1
• Using the new level set function u at time n1,
  solve the diffusion equation for ? with the mixed
  BC at time n1 (which depends on ? at time n1)
• Update the normal velocity vn and kink velocity w
  for the next iteration using ?, ?, k at time n1.
• Repeat steps 1) 6)

21
IV. Some Results

• Set up for all the results,
• Our system is 0,100 x 0,100
• Deposition flux F 0 or 1.0
• Diffusion coeff Dt 1E5.0
• Edge diff coeff De 1E3.0
• Intial edge adatom density ?(x,0) 0
• The adatom density ? on the terrace and the kink
  density k
• are initialized according to the initial island
  shape (a circular
• island, a star-shaped island, etc).

22

• Result 1 rate of convergence
• Initial island shape a circular island centered
  at x50,y50
• with radius of 23.1
• initial adatom density ?(x,0) 0.04
• initial kink density k(x,0) tan(?)
• a) Grid size 32x32, 64x64, 128x128
• dt 1E-5
• Number of iterations 10
• Final stopping time 0.0001
• Order of Accuracy (Max Norm) 1.18563
• Order of Accuracy (L1 Norm) 1.3811
• Order of Accuracy (L2 Norm) 1.3847

23

• b) Grid size 64x64, 128x128, 256x256
• dt 1E-5
• Number of iterations 10
• Final stopping time 0.0001
• Order of Accuracy (Max Norm) 0.843912
• Order of Accuracy (L1 Norm) 0.721405
• Order of Accuracy (L2 Norm) 0.828958
• c) Grid size 64x64, 128x128, 256x256
• dt 1E-5
• Number of iterations 20
• Final stopping time 0.0002
• Order of Accuracy (Max Norm) 0.966984
• Order of Accuracy (L1 Norm) 0.891336
• Order of Accuracy (L2 Norm) 0.950801
• d) Grid size 128x128, 256x256, 512x512
• dt 1E-6
• Number of iterations 20

24

• Graphs
• For N50, dx2, dt0.00001, T100 (final time
  0.001)
• initial zeroth levelset

25

• Graphs
• adatom density ?(x,y) at t 0.001

26

• Graphs
• edge adatom density ? at t 0.001

27

• Graphs
• initial kink density k at t 0 kink density
  k at t 0.001

28

• Graphs
• initial angle ? at t 0 angle ? at t 0.001

29

• Result 2 A circular island grow into a square
• Initial island shape a circular island centered
  at x50,y50 with radius of 13.7
• Adatom density ?(x,t) 0.4 is kept constant for
  all time
• Initial kink density k(x,0) tan(?)
• Deposition flux F 0
• Since the kink density is highest at the angles ?
  ?/4 for a
• circular island, and if the flux to the island
  boundary is kept
• constant i.e. ? constant, then we should expect
  that the
• circular island should eventually grow into a
  square island.
• This is because the normal velocity is the
  greatest at ? ?/4
• and smallest at ? 0. And this is what we
  observe in our
• simulation.

30

• For N50, dx2, dt0.00001, T10000 (final time
  0.1)
• island boundary at t0 (red) and at
  t0.1 (blue)

31
V. Work in progress

• Test numerical schemes with different island
  shapes (a star-shaped island, a perturbed
  circular island, a square island), and see how
  corners/kinks behave in the growing process.
• Check order of accuracy for edge adatom density ?
  and kink density k.
• Check rate of convergence for a step edge with
  the analytic solution.
• Test simulation of islands merging and watch out
  for any numerical issues.
• Add more physics to the kinetic model to improve
  the description of how edge adatoms and kinks
  move at the corners.

32

• THE END.