Interface kinetics of freezing and melting of Si and Na

6 May 1996

PHYSICS

LETTERS

A

Physics Letters A 214 (1996) 71-75

ELSEVIER

Interface kinetics of freezing and melting of Si and Na M. Iwamatsu I, K. Horii Deprtment

of Computer Engineering, Hiroshima City University, 151-S Ozuka, Numata-cho, Asaminami-ku, Hiroshima 731-31, Jupun Received 15 January 1996; accepted for publication 1 March 1996 Communicated by J. Flouquet

Abstract We examine the interface kinetics of freezing and melting of semiconducting Si with a diamond structure and metallic Na with a bee structure using the time-dependent Ginzburg-Landau or Cahn-Allen equation. The various parameters of this free energy are related to the structural data at the liquid phase of those materials from density functional theory. The asymmetry of the interface velocity between superheating and undercooling is partly accounted for by the asymmetry of the functional form of the free energy. PACS: 64.70.D~; 8 I. 10 Kqworrls: Freezing; Crystal; Melt; Density functional

1. Introduction

The microscopic steady state interface kinetics of freezing and melting has been the subject of continuous interest for nearly a century [ I 1. Traditionally theoretical approaches fall into two categories: the one based on the phenomenological transition state theory originally proposed by Wilson and Frenkel [I] and its various modifications [2,3] and the more microscopic theory based on the so-called time-dependent Ginzburg-Landau [ 3,4] or the equivalent Calm-Allen [ 5,6] equations. The latter, though more microscopic, assumes order parameters whose physical meaning is not necessarily clear. Furthermore the theory usually assumes a very simple model free energy which contains parameters whose actual magnitudes are also not clear.

However, the recent advent of the density functional theory of freezing [ 7,8] enables one to deduce the model parameters of this simple free energy from first principles or from experimental data. In this short communication, we are going to examine the applicability of the latter formalism to the freezing and melting of silicon (Si) and sodium (Na). The solid-liquid transition can be studied in a unified formalism based on the density functional method [ 7,8]. In particular, various microscopic parameters which appear in the Ginzburg-Landau or the Calm-Allen equations can be deduced from the scattering experiments on the liquid phase for a bee element [9] and the semiconductors Si and Ge with a diamond structure [ IO]. Recently we have shown [ 111 that using the Ginzburg-Landau equation with parameters deduced from scattering experiments, it is possible to reproduce the microscopic calculation [ 121 of the free energy of the formation of homogeneous nucleation of Na. In this note, we will show that the observed asymmetry of the freezing

’ E-mail: [email protected] 0375-960

l/96/$12.00 Copyright 0 1996 Published by Elsevier Science B.V. All rights reserved PII SO375-9601(96)00163-6

and melting velocities is partly accounted for by the asymmetry of the functional form of the free energy.

2. Steady-state

where C” denotes the second derivatived2C( K) /dK’. At the melting point T = 7’,,,,from the condition of the crystal-melt coexistence of the uniform solid and liquid phase, we obtain [ 91

inter-facial velocity @ = 2( U&,UJ) ‘I2

The order parameters of the crystal-melt transition arc the Fourier components IRK of the local density u(r).

and the equilibrium u(r) of the crystal

(6) (uniform)

order parameter

(7)

140 = L42n,/U&

where K denotes reciprocal lattice vectors of the crystal phase, II~ and II, denote the average densities of the liquid and solid. respectively. In bee elements such as Na and the semiconductors Si and Ge with a diamond structure, only one Fourier component M= UK, which is sometimes called “crystallinity”, is suflicient to retain in ( I ). Relevant reciprocal lattice vectors are the K = ( I IO) lattice for the bee structure [9], and the K = (220) lattice for the diamond structure [ IO]. Then, a simple Landau type expansion of the free energy density of an inhomogeneous solid relative to that of the liquid AF in terms of this order parameter is given by [9]

where ~0 = II, is the number density of the crystal, kT is the temperature, ~2, ~3, ~14and 0 are the expansion coefficients. We have neglected the difference of the densities between the solid phase 11% and the liquid phase /I/, which amounts at most to a few percent and can be included by renormalizing the expansion coefticicnt (14 [ 91. The expansion coefficient a2 is related to the structure factor S(K) and the corresponding direct correlation function C(K) of the liquid phase 17.8 I S(K)

= [ I -C(K)]-’

(3)

ua =

where Q,, is the u2 at T = T,,. The equilibrium order parameter of the melt is 14 = 0. Stroud and coworkers [ 9, IO] have used this formalism to calculate the crystal-melt interfacial tension at the crystal-melt phases coexistence. In order to study the crystal growth at undercooling and the melting at superheating, we will assume furthermore that the above free energy (2) depends on temperature through the height of the main peak of structure factor S(K). Therefore only the coefficient (12(T) is a function of temperatureT [ 11,121. We will further expand the structure factor near the melting point as follows,

Inserting (8) into (2), we obtain a desired model free energy functional which can describe the crystalmelt interface of off-coexistence (superheated or undercooled) conditions. Since the order parameter is nonconserved, the equation of motion for the order parameter u(T) is given by [4,6] &I at’

---

r

6AF (9)

6u ’

n&T,

where the free energy A F is given by (2). This equation gives the time-dependent Ginzburg-Landau equation. A steady state solution of (9) of the form

through u(Z)

(I? = 12/S(K),

(4)

= u(z

~ rt)

(10)

is given by [4]

and the square gradient coefficient 0 is given by UO I4(Z ~ IV) =

h = ~--C”(K).

(5)

(11)

1 +exp[ua@(i

P~lt)l

and the steady state velocity analytic formula

~1 is also given by an 0

4

0.2

X c/m-

I).

(12)

f($> o -0.2

The order of magnitude of da*,/dT is N IO-‘/K and 6 is N 10-s m. The free energy of the form (2) has a strong asymmetry in the stability of the liquid and solid phases. The stability of these phases is seen from the local part of the free energy of (2)

(13)

f(4)

=e&+&(I

is the nondimensional

-&2

(14)

free energy density and

I$ = U/110

(15)

is the nondimensional order parameter, where we have introduced the parameter E defined by

( 16) which corresponds to the temperature. The nondimensional I’rcc energy density f( 4) has three extrema at 4 +. 4, given by (17)

(18) and c,b = 0. For superheating

below E, = $ (0 < t < $) the solid phase 4+ is metastable and the liquid phase 4 = 0 is stable. Above this temperature E, the solid phase is no longer stable but becomes completcly unstable. Therefore E, is the classical spinodal temperature (Fig. I ) . On the other hand, for undercooling above up = -I ( -I < c < 0) the liquid phase 4 = 0 is metastable and the solid phase 4+ is stable. Below EJ, the liquid phase becomes completely unstable but a new

-0.4 0

0.5

I

0 Fig. 1.Free energy density f(d) at various tempemtures l . c > 0: above the melting point: c < 0: below the melting point.

metastable state at 4 = $_ appears (Fig. I ), which always survives at lower temperatures below E! = - I. Therefore the liquid phase remains metastable until fairly deep undercooling (E > -1). Again, E( is the classical spinodal temperature (Fig. I ). Below this temperature a new metastable state at 4 = & might correspond to the amorphous or glassy state. Clearly the unstable (spinodal) state is closer to coexistence for a superheated solid than for an undercooled liquid. This strong asymmetry of the free energy for superheating and undercooling comes from the fact that temperature influences the quadratic term 11’ of the order parameter. This is a nice mathematical realization of the physical fact that the liquid to solid transition involves symmetry breaking but the solid to liquid transition does not. Therefore the latter is more easily realized. This is in great contrast with the liquidvapor nucleation [ 131 where the control parameter of the phase transition, the chemical potential, influences the linear term of the order parameter. Recently we showed that this strong asymmetry has a direct influence on the nucleation behavior [ I 11 using this density functional model. This strong asymmetry directly reflects on the behavior of the interface velocity of the freezing and melting front given by ( 12). In fact at the spinodal Ts (E = $), where the solid phase becomes unstable, we find the melting velocity

u(T,)

= -ii-

( >&g, f

( 19)

m

and at the spinodal Tp (E = -I ), where the liquid phase becomes unstable, we find the freezing velocity

Table I Used parameter

Fig. 2. The interfacial velocity for Si. The melting temperature is 16.50 K. The dots are the result of molecular dynamics simulation [ I7 I. The full line is the theoretical curve ( 12). The parameters which were used to draw the theoretical curve are summarized in Table I.

for interfacial

element

$1’6

Si Na

80 450

-g2-

c*(T) = +-

f (

&= m>

-2&)(T).

(20)

Therefore, for the same amounts of superheating and undercooling, the interface velocity of melting is generally faster than that of freezing. Recently several authors observed such an asymmetry in pulsed laser annealing experiments on silicon (Si) [ 14,151 and germanium (Ge) [ 161. A similar asymmetry has been found by molecular dynamics experiments on Si [ 17 1, Na [ 181 and a Lennard-Jones liquid [ 191. In Figs. 2 and 3 we show the theoretical formula ( 12) compared with the molecular dynamics simulations of Si ]I71 and Na [ 181 using the parameters given in Table I. The experimental results for Si [ 141

&+

x IO-’

(K-‘)

2.5 1.0

d=C N

d(gK)2

Fig. 3. The same as Fig. 2. but for Na. The dots are the result of molecular dynamics simulation I I8 I. The melting temperature is 390 K.

(m/s)

( 12)

are known to be well reproduced by the molecular dynamics simulation [ 171. We see that the general trend of the asymmetry is fairly well reproduced by this simple model. Although we did not attempt to make a nonlinear optimization of those parameters, the used values of the parameters are also reasonable compared with those estimated by others [ 12,9,10]. Finally we note that for bee elements at the melting point T,, when we know S(K) - 2.8 [ 121, then it is reasonable to assume a2111 N 4.3. Furthermore, that the structure factor S(K) has a universal form given by the hard-sphere liquid structure factor with the hard-sphere diameter cr. Then 0 =

T (K)

velocity

I .6rp)

(21)

where we have used the facts that d2C/d(gK)2 N - 1.8 from the data of sodium and potassium [ 121 and that the packing fraction 77 = moc3/6 N 0.45. If we assume that these two values uzrn and b are universal, then we can estimate the magnitude of the velocity from a universal formula,

(22) and ( 12). Therefore the steady state velocity is roughly proportional to the lattice constant. The crystal is denser when the velocity is lower if the microscopic mobility r is constant.

3. Conclusion In this note a mathematically simple density functional model of a crystal-melt phase transition was used to examine the steady-state velocity of the crystal-melt transition using the time-dependent Ginzburg-Landau equation. This model shows a remarkable asymmetry between freezing (melt to crystal transition) and melting (crystal to melt transi-

M. Iwnmtr~srr. K. Horri / Physics Letters A 2 I4

tion), which reflects the same asymmetry seen in the velocity of the interface. Furthermore, we found that: ( 1) The spinodal for superheated solids should be closer to the undercooled liquid, and the maximum superheating is i of the maximum undercooling. ( 2) The above asymmetry reflects the interface velocity, and the melting interface is faster than the freezing interface. Although this density functional model can qualitatively explain the observed interfacial velocity from computer simulation and experiments, it certainly cannot include the complex dynamic process of the atomic Icvel. In particular the atomic arrangement near the surface, which is shown to be responsible for the crystal-fact dependence of the velocity [ 16,191, cannot be taken into account properly. A complementarity analysis using the molecular dynamics or the Monte Carlo simulations is certainly necessary to look into the detailed atomic process of these phenomena.

12 I J.Q. Broughton,

1 I 1 H.A.

Philos.

Mag. SO (1900)

Phys. Z. Sowjetunion

I

151 J.W. Cahn and J.E. Hilliard, 161 S.M. Allen

498

Phys. Rev.

J. Chem. Phys. 96 ( 1992) 5755.

J. Chem. Phys. 28 (1958)

and J.W. Cahn, Acta. Metal.

17 I T.V. Ramakrishnan

and M. Yussouff,

27 (1979)

258. 1085.

Phys. Rev. B 23 ( 1979)

5871. 18 I D.W. Oxtoby

and A.D.J. Haymet,

J. Chem. Phys. 76 ( 1982)

6262.

[ 91

W.H. Shih, Z.Q. Wang, X.C. Zeng and D. Stroud, Phys. Rev. A 39 (1987)

261 I.

I01 X.C. Zeng and D. Stroud, J. Phys. Condens. Matter

I ( 1989)

1779.

11I

M. lwamatsu

and K. Horii,

121 P Harrowell

unpublished

and D.W. Oxtoby,

( 199.5).

J. Chem.

Phys. 80 (1984)

1639. [ 131 M. Iwamatsu,

[ I41

C.J. Calvin, 46 (1985)

I I.51

I I ]

Matter

5 (1993)

and P.S. Peercy, Appl.

7537. Phys. Lett.

644. J.Z. Tischler

and D.M.

Mills.

J. Mater.

Res. I

144.

J.Z. Tischler, M.D.

J. Phys. Condens. J.W. Mayer

B.C. Larson, (1986)

I17 118 II9

238;

(1932)

and P.R. Harrowell,

14) S.K. Ghan, J. Chem. Phys. 67 ( 1977)

52 (1988)

Wilson,

J. Frenkel.

Oxtoby

and K.A. Jackson,

1496.

3834.

II6 I References

G.H. Gilmer

Lett. 49 ( 1982)

( 3 1 D.W.

15

CI9961 71-75

B.C. Larson and D.M. Mills,

Appl.

Phys. Lctt.

1785.

Kluge and J.R. Ray Phys. Rev. B 39 (1989)

C.J. Tymczak

E. Burke, J.Q. Broughton 89 (1988)

1738.

and J.R. Ray, J. Chem. Phys. 92 ( 1990) 7.520.

1030.

and G.H. Gilmer,

J. Chem. Phys.