Eur. Phys. J. Special Topics 177, 5–21 (2009) c EDP Sciences, Springer-Verlag 2009
DOI: 10.1140/epjst/e2009-01165-8

THE EUROPEAN PHYSICAL JOURNAL SPECIAL TOPICS

Regular Article

Homogenization and two-scale models for liquid phase epitaxy Ch. Eck1,a and H. Emmerich2 1 2

IANS, University of Stuttgart, Germany CME/GHI, RWTH Aachen University, Germany

Abstract. We consider models for liquid phase epitaxy without and with elasticity. The models are based on continuum models for fluid flow and transport of adatoms in the liquid solution and a BCF–model for the growth of the solid phase. Using homogenization by formal asymptotic expansion, we obtain two–scale models that are appropriate to describe the evolution of microstructures in the solid phase for processes of technically relevant macroscopic length scales. The two–scale models consist of macroscopic equations for fluid flow and solute transport in the liquid and microscopic cell problems for the growth and elastic deformation of the solid. For the case without elasticity and a phase field approximation of the BCF–model, an estimate of the model error is presented.

1 Introduction Liquid phase epitaxy is a technique for the growth of thin semiconductor compounds that combines solution growth and epitaxy and thereby takes advantage of the potentialities of both processes. It proved tremendous superiority to the early attempts to fabricate such thin films by means of e.g. polishing bulk crystals. Despite this technological importance there is still a lack of well-founded models that take into account all relevant processes. Available models typically only consider a part of the full process, as e.g. dimensionally reduced problems, models with quiescent liquids, models without elasticity, etc. In this article we consider more sophisticated models that combine a full continuum mechanical description of fluid flow and diffusion of adatoms in the liquid phase with a description of the attachment of atoms on the solid surface, the surface diffusion and the growth of monoatomistic steps in the solid described by a BCF model. In a first model we neglect the possible elastic deformation of the growing solid phase. In a second model we also include elastic deformation that may be caused by a misfit of lattice parameters between substrate and epitaxial layer. In many applications of epitaxy one observes a specific microstructure of the liquid–solid phase interface consisting e.g. of spirals, pyramids, or steps. Due to the small scale of these structures it is rather complicated to simulate their evolution within one model for an application problem of a technically relevant length scale. For this purpose we develop homogenized models that combine macroscopic equations for fluid flow and diffusion of adatoms in the liquid and microscopic equations for the evolution of the microstructure and the elastic deformation. This is done by combining a traditional homogenization approach for the elastic solid and the BCF model with a boundary layer expansion for the equations in the liquid. We discuss the approximation of the BCF free interface problem by a correspoding phase field model. Such a model may be useful because it simplifies both the mathematical analysis and numerical simulation. a

e-mail: [email protected]

6

The European Physical Journal Special Topics

Liquid solution

QL

S0

Epitaxial layer QE

S

Fig. 1. Liquid phase epitaxy.

2 The model We consider the growth of a thin epitaxial layer at the bottom of a domain Q ⊂ R3 filled with a liquid solution that contains the adatoms as depicted in Fig. 1. The adatoms in the liquid solution are transported to the solid layer by diffusion and convection, there they may deposit to the surface of the solid layer. On the surface the atoms diffuse until they either reach a monoatomistic step or they desorb to the solution. The process is described by a free boundary problem with two free boundaries: the interface S between liquid solution and epitaxial layer and the interfaces Γ of the monoatomistic steps on the epitaxial surface. For a time t ∈ I with time interval I = (0, T ) the domain Q is partitioned into a domain QL (t) ⊂ Q filled by the liquid solution and a domain QE (t) filled by the epitaxial solid phase. The corresponding time-space domains are ZL := {(t, x) | t ∈ I, x ∈ QL (t)} and ZE := {(t, x) | t ∈ I, x ∈ QE (t)}. The interface S = S(t) between QL and QE is represented by the graph of a function over the bottom surface S0 of Q, S(t) = {x + h(t, x)e3 | x ∈ S0 } , (1) where e3 is the 3rd unit vector. The liquid and solid domains are therefore given by QL (t) = {x ∈ Q | x3 > h(t, x)}

and

QS (t) = {x ∈ Q | x3 < h(t, x)} . The process of liquid phase epitaxy is described by the following equations: – a Navier–Stokes system for the fluid flow, ∇·v =0 ∂t v + (v · ∇)v − η∆v + ∇p = 0

(2)

with fluid velocity v : ZL → R3 , pressure p : ZL → R and given (kinetic) viscosity η of the liquid, – a convection–diffusion equation ∂t cV + v · ∇cV − DV ∆cV = 0

(3)

for the mass specific concentration cV : ZL → R in the liquid solution, where DV is a diffusion constant, – a Burton–Cabrera–Frank (BCF) model for the diffusion of adatoms on the epitaxial surface and the growth of the monoatomistic steps, ∂t cS = DS ∆cS +

1 V 1 c − cS τV τS

for t ∈ I, x ∈ S0 \ Γ(t),

cS = ceq (1 + κΩγ/(kB T ))
vΓ = DS Ω

∂cS ∂n

and

(4) (5)

 for t ∈ I,

x ∈ Γ(t).

(6)

Advances in the Multi-Scale Computational Design of Condensed Matter Interfaces

A B

−vS n

h

S

7

n

Fig. 2. Conservation of matter at the moving interface.

Here cS is the surface density of adatoms on the surface S0 , measured by the number of atoms per unit area, DS is a surface diffusion coefficient, τV is related to the mean time for the deposition of adatoms from the solution to the surface, τS is the mean time for the desorption of adatoms from the surface to the solution, ceq is the equilibrium surface concentration of atoms at the monoatomistic steps, κ is the curvature of the steps, Ω is the area of a single atom, γ is the step stiffness, kB T describes the thermal energy for a fixed  S denotes the difference temperature T , and vΓ is the velocity of the steps. The bracket ∂c ∂n of the normal derivatives on both sides of the interface;
S ∂c = ∇cS+ · n+ + ∇cS− · n− , ∂n where cS± (x) = lim cS (x + rn± ) for x ∈ Γ with the normal vector n+ = n of Γ and n− = −n. r→0

The volume equations (2), (3) and the surface equation (4) are coupled by boundary conditions that describe the conservation of total mass and mass of the adatoms. Let us consider a small box B at the interface with height h = −vnS ∆t and surface area A, as depicted in Fig. 2. Here, ∆t is a small time increment and vnS the velocity of the interface, measured in the direction of the normal vector of S that points to the exterior of QL (t). The box B is filled by liquid material at time t and by solid material at time t + ∆t. Let V and E denote the densities of the liquid and solid adatom phase. At time t the box B contains the mass V |B|, at time t + ∆t the mass is E |B| with volume |B| of the box. The total mass that flows into the box is given by vn V ∆t A. Together with the relation |B| = −AvnS ∆t the conservation of mass is V vn = (V − E )vnS .

(7)

The conservation of mass for the adatoms only is derived in an analogous way: The mass of adatoms in the box is |B|V cV at time t and |B|E at time t+∆t. During the time increment ∆t,  V  the mass flux of atoms into the box is A ∆t cV V vn − DV ∂c ∂n . Hence the interface condition is   ∂cV DV = cV V vn − vnS + vnS E . (8) ∂n The rate of solidified mass in the epitaxial layer is −E vnS . By comparison with the transfer V S term τcV − τcS in the BCF-model (4) we get −E vnS

mA = JS



cV cS − τV τS

 ,

(9)

where mA is the mass of a single atom and JS is the density of the surface measure for the surface S, parametrized over S0 as described in (1). Obviously JS = 1 + |∇h|2 . This factor arises, because −E vnS describes the rate of attached solid mass as a surface density with respect mA S A V to S, while m τV c − τS c corresponds to the surface density related to S0 . As a consequence relations (7) and (8) can be modified to   V  1 c mA 1 cS vn = − − and (10) JS V E τV τS

8

The European Physical Journal Special Topics

 mA  ∂cV 1 − cV = DV ∂n JS



cS cV − τS τV

 .

(11)

The tangential fluid velocity at the interface S vanishes, vt (t, x) = 0 for t ∈ I, x ∈ S(t). The movement of the interface S is given by the condition   mA cV cS ∂t h = − . E τV τS

(12)

The described equations are completed by initial conditions v(0, x) = v0 (x) andcV (0, x) = cV0 (x) for x ∈ QL (0), cS (0, x) = cS0 (x) for x ∈ S0 , by an initial partition of Q into a liquid domain QL (0) and a solid domain QE (0) with interface S(0) that may be defined by an initial condition for h(0, x), by initial steps Γ(0) on the epitaxial layer S0 , and by boundary conditions for the volume equations (2), (3) and the surface equation (4). We may prescribe velocities and forces on different parts of the boundary, v(t, x) = b(t, x) for t ∈ I and x ∈ Γ1 (t)

η    Dv + (Dv) n − pn (t, x) = g(t, x) for t ∈ I and x ∈ Γ2 (t) 2 with Γ1 (t) ∪ Γ2 (t) = ∂QL (t) \ S(t), and homogeneous Neumann conditions for the diffusioin fields ∂cV DV (t, x) = 0 for t ∈ I and x ∈ ∂QL (t) \ S(t), ∂n DS

∂cS (t, x) = 0 for t ∈ I and x ∈ ∂S0 . ∂nS

3 Homogenization and two-scale model The described process of liquid phase epitaxy involves various different length scales: the height of the epitaxial layer is measured in atom diameters, while the continuum models for diffusion and fluid flow, and the model for the adatom diffusion on the epitaxial surface, are valid for a much larger scale. Moreover, the epitaxial surface exhibits a specific microstructure, whose scale is larger than that of an atom diameter but typically much smaller than that of the diffusion and fluid flow processes. In this section we use homogenization techniques to derive a model that is suitable for such situations. It is necessary to choose suitable scaling properties of some coefficients in dependence of a scale parameter ε that represents the size of the epitaxial microstructure. The size of the adatoms is scaled proportional to ε, this leads to Ω = ε2 Ω0 ,

hA = εh0A and mA = ε3 m0A .

(13)

The epitaxial microstructure of scale ε requires the relations DS = ε2 DS0 ,

τV = ε3 τV0 ,

ceq = ε−2 c0eq and γ = ε−1 γ0 .

(14)

In the limit ε → 0 two phenomena must be described: oscillations of smaller and smaller scale in the epitaxial microstructure, and boundary layers for the solutions of (2), (3) at the interface, where these oscillatory data on the epitaxial surface are transported into the domain. These

Advances in the Multi-Scale Computational Design of Condensed Matter Interfaces

9

phenomena require a combination of the standard homogenization technique for the oscillations and the technique of matched asymptotic expansions for the computation of the boundary layer. The homogenization is done for a simple periodic setting with periodicity cell Y ⊂ R2 of area 1. A periodicity cell is a bounded Lipschitz domain with the property that R2 can be represented as union of shifted copies of Y with empty intersection,   R2 = z + Y , z1 + Y ∩ z2 + Y = ∅ for z1 = z2 ∈ IY , z∈IY

where IY ⊂ R2 is a countable set of shifts. For the simplest example with the unit square Y = (0, 1)2 the set IY is Z2 with the set Z of integers. Such a periodic situation is generated by periodic initial conditions of the type S ε (0, x)

= ε−2 cS0 (x, x/ε) forx ∈ S0 and

cVε (0, x) = ε cV 0 (x) for x ∈ QL (0).

(15)

The function cS0 (x, y) is periodic in its second variable y with a periodicity cell Y ; cS0 (x, y) = cS0 (x, y + z) for y ∈ R2 , z ∈ IY . The index ε here and in the sequel indicates the problem of scale ε. The scaling of cS by ε−2 makes sense, because cS measures the number of atoms per unit area, and the area of single atoms is proportional to ε2 . Hence the scaling by ε−2 just keeps the area covered by adatoms constant. The scaling of cV is motivated by the exchange of matter between the solid and the liquid phases. The thickness of the solid phase is proportional to ε, therefore the total mass of adatoms that desorb to the surface is also proportional to ε. The scaling by ε accounts for this rate of mass exchange; it is reasonable for a dilute solution. The outer expansions for the volume fields v, p and cV are valid far away from the surface Sε (t). They are given by the power series vε (t, x) = V 0 (t, x) + ε V 1 (t, x) + · · · , pε (t, x) = P 0 (t, x) + ε P 1 (t, x) + · · · , cVε (t, x) = ε C0V (t, x) + ε2 C1V (t, x) + · · · . The notation by capital letters is done in order to distinguish the outer from the inner expansion presented below. It turns out that only the coefficients of lowest order ε0 are needed; they solve the system ∇ · V0 = 0, ∂t V0 + (V0 · ∇)V0 − η∆V0 + ∇P0 = 0, ∂t C0V + V0 · ∇C0V − DV ∆C0V = 0 for t ∈ I, x ∈ Q. The inner expansion for v, p, cV is valid close to Sε (t). Since this surface is a priori unknown, the inner expansion is done in terms of the modified variable x

= x − hε (t, x)e3 . Let us, for a moment, omit the index ε there. To (t, x) → u(t, x) we denote by (t, x

) → u

(t, x

) the transformed function, defined via u(t, x) = u

(t, x

(t, x)). The gradient and the Laplacian transform according to ∇u = ∇ u − ∇h ∂3 u

and

− 2∇h · ∇∂3 u

+ |∇h|2 ∂32 u

. ∆u = ∆ u − ∆h ∂3 u The time derivative transforms as

− ∂3 u

∂t h. ∂t u = ∂ t u

10

The European Physical Journal Special Topics

We now omit the tilde at the transformed variable x

and the transformed functions. The transformed volume equations (2), (3) are ∇ · v − ∇h · ∂3 v = 0,

(16)

∂t v − ∂t h ∂3 v + (v · ∇)v − (v · ∇h)∂3 v   − η ∆v − 2(∇h · ∇)∂3 v + |∇h|2 ∂32 v − ∆h ∂3 v + ∇p − ∇h ∂3 p = 0

(17)

∂t cV − ∂t h ∂3 cV + v · ∇cV − v · ∇h ∂3 cV   +DV ∆cV − 2∇h · ∇∂3 cV + |∇h|2 ∂32 cV − ∆h ∂3 cV = 0.

(18)

The interface conditions (10) and (11) transform according to   V  1 c mA 1 cS v= − − n, JS V E τV τS       DV ∇h · ∇cV − 1 + |∇h|2 ∂3 cV = 1 − cV 

This follows with the normal n = √

1 1+|∇h|2



(19)

m A S mA V c − V c τS τ

 .

(20)

 ∇h and JS = 1 + |∇h|2 . The functions v, p, cV −1

are expanded in power series vε (t, x) = v0 (t, x, x/ε) + ε v1 (t, x, x/ε) + · · · , pε (t, x) = p0 (t, x, x/ε) + ε p1 (t, x, x/ε) + · · · , cVε (t, x) = ε cV0 (t, x, x/ε) + ε2 cV1 (t, x, x/ε) + · · · . This expansion combines both the homogenization asymptotics and the inner expansion for the boundary layer. The  former is realized by Y –periodic oscillations of the components (y1 , y2 ) in uV (t, x, y) for u ∈ v, p,cV , the latter via the component y3 . As a consequence, the functions y → uj (t, x, y) for u ∈ v, p, cV }, j = 0, 1, 2, . . ., are defined in the domain Y × (0, +∞), with periodic boundary conditions at the lateral surface (y1 , y2 ) ∈ ∂Y . The condition for y3 → +∞ is derived by matching the inner and outer expansion and the condition for y3 = 0 is given by (19), (20). The height function h = hε is also expanded into a series, hε (t, x) = εh0 (t, x, x/ε) + ε2 h1 (t, x, x/ε) + · · ·

for x ∈ S0 .

(21)

For the adatom density cS we assume a power series of two–scale type, cSε (t, x) = ε−2 cS0 (t, x, x/ε) + ε−1 cS1 (t, x, x/ε) + · · · . The surface Γ = Γε (t) is represented as Γε (t) = S0 ∩



  ε z + Γε (t, εz)

z∈IY

with Γε (t, εz) ⊂ Y . It is assumed that the surfaces Γε (t, x) depend in a suitable continuous sense on t, x and ε, and that they converge to some surface Γ0 (t, x) as ε → 0. For the velocity vΓ ε and the curvature κε asymptotic representations of the following type are assumed vΓε = ε vΓ0 + ε2 vΓ1 + · · · and κε = ε−1 κ0 + ε0 κ1 + · · · .

Advances in the Multi-Scale Computational Design of Condensed Matter Interfaces

11

In order to be precise, the surface Γε (t, x) is represented by the image of Γ0 (t, x) for a suitable function π(t, x, ·) : Γ0 (t, x) → Y ;   Γε (t, x) = π(t, x, y) | y ∈ Γ0 (t, x) . Then vΓε and κε are evaluated at (t, x, π(t, x, y)) while the coefficients vΓj , κj are evaluated at (t, x, y) for y ∈ Γ0 (t, x). A suitable choice for π(t, x, y) could be the (closest) intersection of the normal to Γ0 (t, x) at y with Γε (t, x). These formal expansions are used in the relations (4), (5), (6), (12), (16)–(20), and the different orders of ε are collected. We first consider the equations on the surface S0 ⊂ R2 . The O(ε−2 ) terms of (4) and (5) as well as the O(ε)-terms in (6) give a microscopic BCF-model, ∂t cS0 − DS0 ∆y cS0 =

1 V 1 c − cS0 on Y, τV0 0 τS

cS0 = c0eq (1 + κ0 Ω0 γ0 /(kB T ))  S ∂c vΓ0 = DS0 Ω0 ∂n0

(22)

 on Γ0 (t, x).

(23)

This problem is given in the periodicity cell, y ∈ Y ; there is one such problem for every point x ∈ S0 . The function cV0 here is the trace of the leading order term from the inner expansion at y3 = 0. Obviously the leading order terms here already deliver a full BCF-model on Y , with cV0 interpreted as part of the given data. It therefore suffices to consider (16)–(18) and the coupling conditions (19), (20). The terms of leading order in (16)–(18) yield ∇y · v0 − ∇y h0 · ∂y3 v0 = 0,

(24)

∆y v0 − 2(∇y h0 · ∇y )∂y3 v0 + |∇y h0 |2 ∂y23 v0 − ∆y h0 ∂y3 v0 = 0 and

(25)

L0 cV0 = 0

(26)

in Y × (0, +∞) with the leading order differential operator of the convection–diffusion equation L0 = ∆y − 2∇y h0 · ∇y ∂y3 + |∇y h0 |2 ∂y23 − ∆y h0 ∂y3 . These equations are supplemented by periodic boundary conditions at y ∈ ∂Y × (0, +∞), and by boundary conditions for y3 = 0 and y3 → +∞. The conditions at y3 = 0 follow from the O(ε0 )-terms in (20), (19): v0 = 0, B0 cV0 = 0 with leading order boundary operator of the convection–diffusion equation     B0 = DV ∇y h0 · ∇y − 1 + |∇y h0 |2 ∂y3 . The conditions for y3 → +∞ are lim v0 (t, x, y) = V0 (t, x) and

y3 →+∞

lim cV0 (t, x, y) = C0V (t, x).

y3 →+∞

Since (25) and (26) are elliptic equations of second order, it is easy to conclude that v0 = 0 and that cV0 is independent of y. This gives the required boundary condition V0 = 0 on S0

12

The European Physical Journal Special Topics

for the outer expansion of the velocity field and also shows that cV0 = cV0 (t, x) is a purely “macroscopic” variable. It remains to consider the convection–diffusion equation (18). The O(ε0 )-terms in (18) and the O(ε)-terms in (20) give L0 cV1 − ∆y h0 ∂x3 cV0 = 0 in Y × (0, +∞) B0 cV1 + B1 cV0 =

and

(27)

m0A S m0A V c − 0 c0 on Y × {0} τS 0 τV

(28)

with

    B1 = DV ∇y h0 · ∇x − 1 + |∇y h0 |2 ∂x3 . The boundary conditions at the remaining parts of the boundary are periodic boundary conditions for (y1 , y2 ) ∈ ∂Y as well as the condition lim cV1 (t, x, y) = C1V (t, x) which results y3 →+∞

from the matching to the outer expansion. From this problem we derive the boundary condition for the leading order term C0V of the outer expansion. This is done by combining the above problem with the leading order matching condition for the fluxes, DV ∂3 C0V = DV ∂x3 cV0 + DV Since C0V and cV0 are independent of y, the limit

lim ∂y3 cV1 .

y3 →+∞

lim ∂y3 cV1 is also independent of (y1 , y2 ). In

y3 →+∞

order to derive the desired condition, we integrate relation (27). This integration is done with respect to the variable y = y + h0 (y1 , y2 )e3 . The domain of integration is    Y := y ∈ R3 | ( G y1 , y2 ) ∈ Y, h0 ( y ) ≤ y3 ≤ yM in terms of the variables y and   GY := y ∈ R3 | (y1 , y2 ) ∈ Y, 0 ≤ y3 ≤ yM − h0 (y) in terms of y. In the following the parameter yM will take the limit yM → +∞. Thus we obtain       L0 cV1 − ∆y h0 ∂x3 cV0 dy = ∆y
 y c1V − ∆y
h0 ∂x3  c0V d 0=
Y G

GY

 =

Y ×{yM }

∂y
3  c1V dsy






+ Y ×{0}

   ∇y cV1 · ∇y h0 − ∂y3 cV1 1 + |∇y h0 |2 − |∇y h0 |2 ∂x3 cV0 dsy .

Here the last integral is written in terms of the variables y. Combined with the boundary condition (28) the previous relation leads to  ∂y
3  c1V dsy
Y ×{yM }



=−




Y ×{0}

DV−1



m0A S m0A V c − 0 c0 τS 0 τV



− ∇y h0 · ∇x cV0 + ∂x3 cV0

 dsy .

 Here the term Y ×{0} ∇y h0 · ∇x cV0 dsy cancels due to the periodic boundary condition for h0 and the fact that cV0 is independent of y. From the matching condition we now conclude  V V ∂y
3  c1V dsy
DV ∂x3 C0 |x3 =0 = DV ∂x3 c0 + DV lim yM →+∞

 =− with the cell average cS0 =

 Y

cS0 dsy .

m0A S c τS 0

−

m0A V c τV0 0



Y ×{yM }

Advances in the Multi-Scale Computational Design of Condensed Matter Interfaces

13

Let us sum up the result. We obtain the following two–scale model: – Macroscopic equations for fluid flow and solute diffusion ∂t cV0 + v0 · ∇cV0 − DV ∆cV0 = 0,

(29)

∇ · v0 = 0,

(30)

∂t v0 + (v0 · ∇)v0 − η∆v0 + ∇p0 = 0

(31)

to be solved on I × Q. These equations follow from the outer expansion, if the capital letters C0V , P0 and V0 there are replaced by lower case letters cV0 , p, v. – A microscopic BCF-model 1 V 1 c − cS0 for y ∈ Y, τV0 0 τS  cS0 = c0eq (1 +
κ0SΩ0 γ0 /(kB T ))  ∂c0 for y ∈ Γ0 (t, x), vΓ0 = DS0 Ω0  ∂ny ∂t cS0 − DS0 ∆y cS0 =

to be solved for every x ∈ S0 on I × Y . – Coupling conditions  S  ∂cV0 c0 cV0 0 DV = mA − 0 , ∂n τS τV

(32)

(33)

(34)

v=0 (35)  S = Y c0 (t, x, y) dy. These equations serve as boundwith the microscopic mean value ary conditions for equations (29) and (31) on the part S0 of the boundary. cS0 (t, x)

It is easily observed that the equations for fluid flow decouple from the other equations. The fluid velocity can be therefore computed in a preprocessing step, and the main computation then consists in solving the remaining equations with given fluid velocity. The two-scale model describes the limit ε → 0 of the model with scale ε. In reality, ε will never be zero. The two-scale model is an approximation for the problem of scale ε, this approximation becomes more accurate as ε becomes smaller.

4 Phase field model In the phase field version of the original BCF-model (4)–(6) we introduce a continuous function Φ that approximates the number of monoatomistic steps. In a static setting, the phase field minimizes the free energy    
ξ 1 S 2 S |∇Φ| + f (Φ) − kB T (c − ceq )Φ dx. F(c , T, u, Φ) = ceq γΩδ 2 ξ S0 The function f here is a multi–well potential with minima at integer values, for example f (Φ) = − cos(2πφ); the parameter ξ controls the thickness of a diffusive transition layer between the steps. The parameter ∆ is computed via the solution ϕ of the one–dimensional equation −ϕ (x) + f  (ϕ(x)) = 0 with boundary conditions ϕ(−∞) = 0, ϕ(+∞) = 1, and ϕ(0) = 1/2 by  +∞    (ϕ (x))2 + f (ϕ(x)) dx. ∆= −∞

14

The European Physical Journal Special Topics

The solution ϕ here determines the shape of a diffuse transition for a monoatomistic step.

 ξ 2 The quantity S0 ceq γΩ∆ 2 |∇Φ| + 1ξ f (Φ) dx defines the free energy of the interface curve  c γΩ dsx . Dynamic relaxation leads to the standard phase field equation Γ eq τ ∂t Φ = −DΦ F(cS , T, u, Φ)   1 = ceq γΩ∆ ξ∆Φ − p (Φ) + kB T (cS − ceq ) ξ or, after a suitable rescaling αξ 2 ∂t Φ − ξ 2 ∆Φ + f  (Φ) + g(cS , Φ) = 0

(36)

with a modified time relaxation parameter α and g(cS , Φ) =

 ξ kB T  S c − ceq . δ ceq Ωγ

(37)

This equation is supplemented by a surface diffusion equation ∂t cS − DS ∆cS + Ω−1 ∂t Φ =

cV cS − . τV τS

(38)

In the phase field version of the two–scale model we replace relation (32)–(33) by the microscopic phase field equations cV0 cS0 ∂t cS0 − DS0 ∆y cS0 + Ω−1 , (39) 0 ∂t Φ0 = 0 − τV τS α0 ξ02 ∂t Φ0 − ξ02 ∆y Φ0 + f  (Φ0 ) + g0 (cS0 , Φ0 ) = 0. cS0

(40) 2 S

corresponds to the limit of ε c and cV0 As in the corresponding BCF-version of the model, −1 V to the limit of ε c . The relation between the original parameters α, ξ and their two-scale counterparts is given by the scaling condition ξ = εξ0 , α = ε−2 α0 . The function g = gε in (36) is assumed to depend on ε via gε (cS , Φ) = g0 (ε2 cS , Φ). For the example (37) we have g0 (cS0 , Φ0 ) =

 ξ0 kB T  S c − c0eq . ∆c0eq Ω0 γ0 0

The phase field version of the two–scale model is therefore given by: – a macroscopic diffusion–convection equation and Navier–Stokes system ∂t cV0 + v0 · ∇cV0 − DV ∆cV0 = 0,

(41)

∇ · v0 = 0,

(42)

∂t v0 + (v0 · ∇)v0 − η ∆v0 + ∇p0 = 0

(43)

to be solved in the liquid domain Q. – the phase field version of a microscopic Burton–Cabrera–Frank model ∂t cS0 − DS0 ∆y cS0 + Ω−1 0 ∂t Φ0 =

cV0 cS0 − , τV0 τS

α0 ξ02 ∂t Φ0 − ξ02 ∆y Φ0 + f  (Φ0 ) + g0 (cS0 , Φ0 ) = 0 n

(44) (45)

to be solved for every x ∈ S0 on a unit cell Y ⊂ R with periodic boundary conditions,

Advances in the Multi-Scale Computational Design of Condensed Matter Interfaces

– and the coupling condition ∂cV DV 0 = m0A ∂n



cS0 cV − 00 τS τV

15

 ,

v0 = 0

(46) (47)

on S0 . It is possible to prove the existence and uniqueness of solutions to this model under appropriate assumption with appropriate boundary- and initial conditions, for details of the mathematical analysis see [12].

5 Model error The models presented above are derived by formal asymptotic expansions. In the derivation we have assumed that expansions in power series of the described form actually exist. It is, however, not guaranteed that this assumption is true. As a consequence, there is some need of further justification of the model. A good method to justify our result is the derivation of an estimate for the model error. Such an estimate is indeed available for the present model. Since the technicalities are rather cumbersome and involve sophisticated mathematical techniques, we report here the results only; the details can be found in [12]. We compare the solutions of the two–scale model (43)–(46) with the solutions of the original problem (2), (3), (4), (5), (6), (10)–(12) with the appropriate scaling (13), (14) of the physical parameters. The latter problem is therefore given by the Navier–Stokes system and the convection–diffusion equation ∇ · vε = 0, ∂t vε + (vε · ∇)vε − η∆vε + ∇pε = 0,

(48)

∂t cVε + vε · ∇cVε − DV ∆cVε = 0

(49)

in the liquid domain Qε = {(t, x) ∈ Q | x3 > hε (t, x)}, the phase field version cVε cSε − , τV0 τS

(50)

α0 ξ02 ∂t Φε − ε2 ξ02 ∆Φε + f  (Φε ) + g(cSε , Φε ) = 0,

(51)

−3 ∂t cSε − ε2 DS0 ∆cSε + ε−2 Ω−1 0 ∂t Φε = ε

of the microscopic BCF-model and the coupling conditions   V  S 1 cε m0A 1 3 cε vεn = − −ε and JSε V E τV0 τS    3 cSε m0A  ∂cVε cVε V = 1 − cε − 0 DV ε ∂n JSε τS τV

(52)

(53)

on S. Since the functions cS0 and Φ0 depend on both space scales x and y ∼ x/ε, we cannot compare them directly. The solution of the two scale model should give an approximation for the solution of the original model with given scale ε, if the microscopic variable y is replaced by y/ε. We therefore define reconstructions S ε ε cS,ε 0 (t, x) = c0 (t, x, x/ε) and Φ0 (t, x) = Φ0 (t, x, x/ε).

Due to the scaling properties (15) we introduce rescaled functions

cSε (t, x) = ε2 cSε (t, x),

cVε (t, x) = ε−1 cVε (t, x).

(54)

16

The European Physical Journal Special Topics

The estimate for the model error is done in terms of the norms 1/2  2 |u(t, x)| dx , uL∞ (I;L2 (S0 )) = supt∈I   uH 1 (I;L2 (S0 )) =

I

 u(t)L2 (Qε (t)) =

S0

S0

1/2 |∂t u(t, x)|2 dx dt ,

Qε (t)

1/2 |u(t, x)|2 dx

.

The estimate is valid under appropriate assumptions on the regularity of the solutions that are described in [12]. It is given by     S   cε − cSε + Φε − Φε0 H 1 (I;L (S )) 0 L (I;L (S )) ∞

2

0

2

   2 cεV − cV0 (t)L + sup  t∈I

2 (Qε

0

 ≤ Cε1/2 . + (v − v )(t) ε 0 L (Q (t)) 2 ε (t))

(55)

This estimate shows an error of order ε1/2 , instead of ε1 as could be expected from the form of the asymptotic expansions. This deterioration has two reasons. First, the macroscopic fluid flow and convection–diffusion processes happen in the whole domain Q in the case of the two–scale model and in the liquid domain Qε (t) for the original model. These domains differ by a volume of size proportional to ε, which leads to an error of order ε1/2 in L2 –type estimates. Second, for a general domain S0 and a given scale parameter ε it must be expected that there are “incomplete” unit cells at the boundary that do not coincide with shifted copies of εY . In this case the √ microscopic problems are not correct at the boundary; and this leads to an additional error of ε in the L2 (S0 )-norm. The advantage of the new two-scale model is the decoupling of scales: the mesh-size of a computational grid for a numerical discretization of the two-scale model can be chosen independently of the scale of the microstructure. The price to pay is an increase of the dimension: the original model consists of partial differential equations on a three-dimensional (spatial) domain, coupled to differential equations on the two-dimensional surface. In the two-scale model, there are still macroscopic differential equations on the three-dimensional fluid volume and the microscopic problems have to be solved on a two-dimensional unit volume, but there is a microscopic problem for every point of the two-dimensional epitaxial surface S0 . This adds up to four dimensions, compared to three dimensions for the original model. Hence it is expected that the two-scale model is superior to the original model for a sufficiently small scale ε of the microstructure only. If the mesh size of a discretization for the original problem is larger than ε, then the original model does not resolve the microstructure and for such a situation the two-scale model is definitely superior. It can be also expected that the two-scale model is superior, if the mesh size for the original model is aε with a not too small number a; then the microstructure is resolved by the original model only very coarsely.

6 Liquid phase epitaxy with elasticity In this section we add elastic deformations of the growing solid layer to the model. Elastic deformations typically occur as a consequence of misfit strains between the solid phase and the substrate. We add elasticity via the equations of linear elasticity, −∇ · σ(u) = 0

(56)

with stress tensor σ(u). The stress tensor σ(u) = (σij (u))3i,j=1 is given by Hooke’s law for linear 3 elasticity σij (u) = k,=1 aijk ek (u) in terms of the linearized strain tensor e(u) = (eij (u))3i,j=1

Advances in the Multi-Scale Computational Design of Condensed Matter Interfaces

17

Γ(S E)

S0\SE SE

Fig. 3. Geometry of a step.

with components eij (u) = 12 (∂i uj + ∂j ui ). The elastic deformation of the solid is generated by the misfit strain between the atoms in the epitaxial solid layer and the atoms of the substrate on which it grows. A simple model for the misfit strain is the boundary condition σ(u)n = b on the bottom of the domain, where b is the misfit stress. To the Burton Cabrera Frank model we add an additional term coming from the elastic energy in the interface condition. The BCF model with elasticity, ∂t cS = DS ∆cS +

1 V 1 c − cS τV τS

cS = ceq (1 + κγΩ/(kB T )) +


∂cS vΓ = DS Ω ∂n

for x ∈ S0 \ Γ(t),

hA σ(u) : e(u) 2kB T

and

(57)

(58)

 for x ∈ Γ(t).

(59)

3 The symbol : denotes the scalar product of two tensors, σ(u) : e(u) = i,j=1 σij (u)eij (u). The condition at the free boundary is obtained by minimization of a free energy functional   hA σ(u) : ε(u) dx kB T (cS − eeq ) dx + F(cS , T, u, SE ) = S0 \SE SE 2  + γΩceq dsx Γ(SE )

with respect to variations of the surface part SE covered by the solid phase, see also Fig. 3. This functional is composed of the free energies of the surface diffusion process (the first integral), the elastic deformation (the second integral) and the free boundary (the third integral). The Eqs. (2), (3) in the liquid, the coupling conditions (10) and (11), and the condition (12) for the movement of the free boundary remain unchanged. This implies that the elastic deformation does not influence the liquid domain. This is an approximation that is justified, since the elastic deformation is small compared to the overall growth of the solid layer. The boundary condition for the elasticity equations follows from the equilibrium of forces between the fluid and the elastic solid,  1  σ(u)nS + η ∇v + (∇v) nL − p nL = 0. 2

(60)

The index L or S at the normal vector n indicates its orientation: nL points to the exterior of QL while nS points to the exterior of QS . In order to perform the homogenization limit for the case with elasticity, we need suitable scaling properties for the elastic field. The displacement field is some fraction of an atom diameter and therefore scales as u ∼ ε. The other quantities scale as described in (13), (14). The outer expansions for the volume fields v, p and cV are done as described in Section 3; the result is (29)–(31). The expansion for the BCF model is also done as in Section 3, the only

18

The European Physical Journal Special Topics

modification is the additional elastic energy on the right hand side of (58). The leading order terms of (57), (58) and (59) yield the microscopic cell problem ∂t cS0 − DS0 ∆y cS0 =

cV0 cS − 0 0 τV τS

on Y,

   γ0 Ω0 h0A  σy (u0 ) : ey (u0 )  cS0 = c0eq 1 + κ0 +  kB T 2kB T  ∂cS   0  vΓ0 = DS0 Ω0  ∂ny

(61)

on Γ0 (t, x).

(62)

The elasticity equations are defined on an ε-dependent domain QSε (t) = {x ∈ Q | x3 < hε (t, x)}. The height function hε is expanded into a two–scale series as in (21). The expansion of the 3 (hε (t, x) − ε h0 (t, x, x/ε))e3 , elasticity equations is done in terms of the variables x

= x − hεx(t,x) 3

Sε (t) = {x ∈ Q | 0 < where e3 is the 3rd unit vector of R . This transform maps QSε (t) to Q x3 < ε h0 (t, x, x/ε)} and therefore reduces the ε-dependence to a simple scaling relation. As a consequence, the problem obtained from the expansion of the elasticity equations will be defined on an ε-independent domain. The transformed function u

is defined by u(t, x) = u

(t, x

(t, x)). The space derivatives transform according to ∂xi u =

3 

∂xj u

∂xi x

j .

j=1 3 Due to the two–scale expansion (21) for hε we find x

= x − hεx(t,x) (ε2 h2 (t, x, x/ε) + · · · ); this implies ∂xi x

j (t, x) = ∆ij + O(ε) for x3 < hε (t, x). As a consequence the leading order term in the asymptotic expansion of the elasticity equations satisfies

−∇y · σy (u0 ) = 0, with σy;ijk (u) = aijk ey;k (u) and ey;ij (u) = 12 (∂yi uj + ∂yj ui ). This equation must be solved for (y1 , y2 ) ∈ Y , 0 ≤ y3 ≤ h0 (t, x, y) with periodic boundary conditions for (y1 , y2 ) ∈ ∂Y . The boundary traction σy (u0 )n will be coupled to the inner expansion of the Navier–Stokes system, on the lower part of the boundary we have the prescribed misfit stress σy (u0 ) · ny = b0 . The higher order terms in the expansion of the displacement field will not be needed. The inner expansion for v, p, cV is again done as described in Section 3. For the convection– diffusion equation we obtain the same result as for the case without elasticity, because elasticity does not couple directly to the convection–diffusion process. For the Navier–Stokes system, the only difference is the boundary condition (60) that describes the equilibrium of forces. The leadin order terms in the Navier–Stokes system and the boundary conditions (10) and (11) lead to (34), (35). The problem of first order in the Navier–Stokes system yields the Stokes system ∇y · v1 = 0, −η∆y v1 + ∇y p0 = 0. This system is supplemented by periodic boundary conditions for (y1 , y2 ) ∈ ∂Y and by the boundary condition   V  1 c0 1 cS0 v1 = JS−1 − − ny for y ∈ SY . 0 V E τV0 τS The limit condition for y3 → +∞ again follows from the matching of the inner to the outer expansion,   ∇V0 = lim ∇x v0 + ∇y v1 and P0 = lim p0 , y3 →+∞

y3 →+∞

Advances in the Multi-Scale Computational Design of Condensed Matter Interfaces

19

where V0 and P0 denote the leading order terms of the outer expansion. This leads to the condition   η η ∇y y1 + (∇y y1 ) · e3 − p0 e3 = ∇V0 + (∇V0 ) · e3 − P0 e3 . (63) lim y3 →+∞ 2 2 The solution of this Stokes system is required to formulate the boundary condition for the elasticity equations  η σy (u0 ) · nyS + ∇y v1 + (∇y v1 ) · nyL − p0 nyL = 0. (64) 2 Summing up all derived relations we obtain the following two-scale model: – Macroscopic equations for fluid flow and solute diffusion in the bulk, ∂t cV0 + v0 · ∇cV0 − DV ∆cV0 = 0,

(65)

∇ · v0 = 0,

(66)

∂t v0 + (v0 · ∇)v0 − η∆v0 + ∇p0 = 0

(67)

to be solved on Q. – Microscopic equations for elastic deformation and fluid flow close to the interface. These problems are given on the domains QY S (t, x) : = {y ∈ Y × (0, +∞) | y3 < h0 (t, x, y)}

and

QY L (t, x) : = {y ∈ Y × (0, +∞) | y3 > h0 (t, x, y)} by −∇y · σy (u0 ) = 0 on QY S (t, x), ∇y · v1 = 0

and

−η ∆y v1 + ∇y p1 = 0 on QY L (t, x), 

σy (u0 )ny = b 



on Y × {0},

∇y v1 + (∇y v1 ) · nyL − p0 nyL = 0 on SY (t, x),

 V  c0 cS 1 1 − τ0S on SY (t, x) and v1 = JS−1 V − E τV0 0     η  · e3 − p1 e3 = η2 ∇v0 + (∇v0 ) · e3 − p0 e3 2 ∇y y1 + (∇y y1 )

σy (u0 ) · nyS +

limy3 →+∞

η 2

with periodic boundary conditions for (y1 , y2 ) ∈ ∂Y . The notation of the pressure here is changed to p1 , to avoid any confusion with the pressure of the macroscopic fluid flow equation. – A microscopic BCF model ∂t cS0 − DS0 ∆y cS0 =

cV0 cS0 − for y ∈ Y, τV0 τS

(68)



cS0

=

c0eq

 γ0 Ω0 h0A σy (u0 ) : ey (u0 ) and 1 + κ0 + kB T 2kB T
S ∂c0 vΓ0 = DS0 Ω0 for y ∈ Γ0 (t, x) ∂ny

to be solved on Y for every x ∈ S0 .

(69)

20

The European Physical Journal Special Topics

– Coupling conditions DV

∂cV0 cS cV = 0 − 0 , ∂n τS τV

(70)

v=0 (71)  S with the microscopic mean value = Y c0 (t, x, y) dy. This equation serves as boundary condition for equation (65) on the part S0 of the boundary. cS0 (t, x)

The macroscopic Navier–Stokes system and convection diffusion equation here decouples again from the other equations. The fluid velocity can be computed in a preprocessing step, and the main computation then consists in solving the remaining equations with given macroscopic fluid velocity.

7 Phase field model The phase field model for epitaxial growth is adapted to the case with elastic deformations by adding elastic energy to the free energy functional. This leads to     ξ 1 S |∇Φ|2 + f (Φ) − kB T (cS − ceq )Φ ceq γΩ∆ F(c , T, u, Φ) = 2 ξ S0  hA + σ(u) : e(u)Φ dx 2 with multi–well potential f that has minima at integer values, e.g. f (Φ) = − cos(2πΦ) and parameter ∆ to be chosen as described in Section 4. The usual time relaxation and rescaling leads to αξ 2 ∂t Φ − ξ 2 ∆Φ + f  (Φ) + g(T, cS , u, Φ) = 0 (72) with g(T, cS , u, Φ) =

ξ ∆



 kB T hA (ceq − cS ) + σ(u) : e(u) . ceq γΩ 2ceq γΩ

The BCF model can be replaced by the phase field equation (72) and a slightly modified surface diffusion equation cV cS ∂t cS + Ω−1 ∂t Φ − DS ∆cS = − , τV τS with additional term Ω−1 ∂t Φ that describes the conservation of adatoms.

8 Summary We presented models for liquid phase epitaxy without and with elasticity. The models are based on continuum models for the flow of the liquid solution, the adatoms in the liquid phase and the elastic deformation, and a BCF model for the growth of the expitaxial solid that is adapted to the presence of elastic deformations. The homogenized models are derived by formal asymptotic expansions, they consist of macroscopic equations for fluid flow and transport of adatoms in the liquid and microscopic models for the growth of the solid. In the model with elastic deformations, the microscopic models also include the elasticity equations and a Stokes system in a semi–infinite domain that describes the elasticity and coupling of microscopic elastic forces to the macroscopic flow field. For both models we presented a version with a phase field model for the growth of the solid phase. For the model without elasticity, a rigorous justification of the homogenization result is available, for the model with elasticity this is still an open question. The homogenized models can be used as an approximation for processes with very small spatial scale of microstructures that arise in the growing solid. The accuracy of this approximation increases with decreasing microscale.

Advances in the Multi-Scale Computational Design of Condensed Matter Interfaces

21

References 1. R. Asaro, W.A. Tiller, Metall. Trans. 2, 1789 (1972) 2. E. Bauser, Atomic mechanisms in semiconductor liquid phase epitaxy, edited by D.T.J. Hurle, Handbook of Crystal Growth, Vol. 3 (North–Holland, Amsterdam, 1994) 3. T. Blesgen, U. Weikard, Electr. J. Differ. Eq. 89, 1 (2005) 4. E. Bonnetier, A. Chambolle, SIAM J. Appl. Math 62, 1093 (2002) 5. W.K. Burton, N. Cabrera, F.C. Frank, Phil. Trans. Roy. Soc. 243, 299 (1951) 6. G. Caginalp, Arch. Ration. Mech. Anal. 92, 205 (1986) 7. M. Carrive, A. Miranville, A. P´etrus, Adv. Math. Sci. Appl. 10, 539 (2000) 8. V. Chalupecki, Ch. Eck, H. Emmerich, A two–scale model for liquid phase epitaxy, Preprint 196, DFG SPP 1095 “Mehrskalenprobleme” (2006) 9. A.A. Chernov, T. Nishinaga, edited by I. Sungawa, Growth shapes and stability in Morphology of Crystals (Terra Scientific Publ. Co., 1987), p. 270 10. W. Dorsch, S. Christiansen, M. Albrecht, P.O. Hansson, E. Bauser, H.P. Strunk, Surf. Sci. 896, 331 (1994) 11. V. Chalupecki, Ch. Eck, H. Emmerich, A two–scale model for liquid phase epitaxy, Preprint 196, DFG SPP 1095 “Mehrskalenprobleme” (2004) 12. Ch. Eck, H. Emmerich, Math. Methods Appl. Sci. 32, 12 (2009) 13. H. Emmerich, Ch. Eck, Cont. Mech. Thermodynamics 17, 373 (2006) 14. H. Emmerich, Cont. Mech. Thermodyn. 15, 197 (2003); Cont. Mech. Thermodyn. 17, 373 (2006) 15. E. Fried, M.E. Gurtin, J. Mech. Phys. Solids 51, 487 (2003) 16. H. Garcke, Ann. Inst. H. Poincar´e, AN 22, 165 (2005) 17. H. Garcke, Proc. R. Soc. Edinb., Sect. A, Math. 133, 307 (2003) 18. H. Garcke, M. Rumpf, U. Weikard, Interfaces Free Bound. 3, 101 (2001) 19. R. Ghez, Int. J. Heat Mass Transfer 23, 425 (1980) 20. M.A. Grinfeld, Doklady Akademii Nauk SSSR 290, 1358 (1986) 21. M.E. Gurtin, Physica D 92, 178 (1996) 22. M.E. Gurtin, M.E. Jabbour, Arch. Ration. Mech. Anal. 163, 171 (2002) 23. L.D. Khutoryanskii, P.P. Petrov, Sov. Phys. Crystallogr. 23, 571 (1978) 24. F.P.J. Kuijpers, G.F.M. Beenker, J. Cryst. Growth 48, 411 (1979) 25. M.B. Small, E. Ghez, E. Giess, Liquid Phase Epitaxy, edited by D.T.J. Hurle, Handbook of Crystal Growth, Vol. 3 (North–Holland, Amsterdam, 1994) 26. V.V Jikov, S.M. Kozlov, O.A Oleinik, Homogenization of Differential Operators and Integral Functionals, (Springer, Berlin – Heidelberg, 1994) 27. W. Lu, Z. Suo, J. Mech. Phys. Solids 49, 1937 (2001) 28. H. M¨ uller–Krumbhaar, J. Chem. Phys. 63, 5131 (1975) 29. W.W. Mullins, R.F. Sekerka, J. Appl. Phys. 35, 444 (1964) 30. E.M. Sparrow, J.L. Gregg, Trans. ASME J. Heat Transfer 82C, 294 (1960) 31. B.J. Spencer, J. Tersoff, Equilibrium shapes of islands in epitaxially strained solid films, edited by K. Golden et al., Mathematics of Multiscale Materials, Proceedings of the IMA Workshops, Minneapolis, 1995–1996, IMA Vol. Math. Appl. 99 (Springer, New York, 1998), p. 255 32. N. Tokuda, J. Cryst. Growth 67, 358 (1984) 33. W.R. Wilcox, J. Cryst. Growth 12, 93 (1972) 34. L.O. Wilson, N.L. Schryer, J. Fluid Mech. 85, 479 (1978) 35. Y. Xiang, SIAM J. Appl. Math. 63, 241 (2002)