Das Sh.P. Statistical Physics of Liquids at Freezing and Beyond

Подождите немного. Документ загружается.

3.3 The density-functional approach 141

The corresponding choice for the set of order parameters represents a saddle point for

[n] in the function space of density. The ordinary differential equations obtained for

the order parameters are given by

∂

∂γ

D[γ ]+

αβ





∂

∂γ

αβ

[γ ]{∇

(x)∇

(x)}

− #

αβ

[γ ]∇

∇

(x) −

∇

αβ

[γ ]

∇

(x)

= 0 (3.3.15)

for the index i =1,...,m running over the set. The above equations represent the

spatial variation of the local amplitudes γ

across the crystal–liquid interface in the

critical nucleus.

From what we have described above, the ﬁrst step in the construction of the density-

functional model of the nucleus in the melt therefore involves constructing (a) the density

function parametrized in terms of a set of local order parameters that are chosen so as to

represent a spherical nucleus in the melt and (b) the proper grand-canonical potential 

as a functional of the inhomogeneous density function. In the simplest approximation the

ﬁrst term in (3.3.13) is computed in terms of a low-order functional Taylor-series expan-

sion of the Ramakrishnan and Yussouff (1979) model described in Chapter 2. However,

this would imply an ad-hoc truncation of the density expansion for the crystalline state in

which the density ﬂuctuations are not small. Alternatively, 

is computed using the more

sophisticated weighted-density-functional methods described in Chapter 2. We consider

both cases here.

3.3.2 The critical nucleus

In the simplest approximation 

is computed (Harrowell and Oxtoby, 1984) in terms of

the expansion (3.3.10) around the uniform liquid state and keeping terms of up to second

order in the density ﬂuctuations around the uniform liquid state. The functional  is

to be optimized with respect to a suitably chosen density function, which represents the

crystal nucleus. We consider for the density n

(r) representing the spherical nucleus the

functional form (3.3.3) with the local amplitudes γ(r) ≡{η(r), A

(r)} being functions

solely of the radial distance r. This implies spherical symmetry for the order parameters

but not for the density function itself. Using the parametrization (3.3.3), we obtain for

the density ﬂuctuation δn

(x;γ) = n

(x;γ)− n

of the inhomogeneous state around the

uniform liquid state density n

δn

;γ(x

)) = n



η(x

) +



· x



. (3.3.16)

142 Crystal nucleation

On substituting (3.3.17) into (3.3.10), the grand potential function 

is obtained as



= n



˜η(x) +



(x)e

· x

⎡

⎣

˜η(x) +



(x)e

· x

⎤

⎦

− n



η(x) +



(x)e

· x

−





c(x, x



;γ)

η(x) +



(x)e

· x

⎧

⎨

⎩

η(x) +



(x)e

· x



⎫

⎬

⎭

, (3.3.17)

where ˜η =1 + η. Next, with the present choice of the density function, we obtain the

following values for the partial derivatives:

∂n

(r)

∂η

= n

, and

∂n

(r)

∂A

= n

· r

. (3.3.18)

Using the results (3.3.18) in the expression (3.3.12) the kernel #

αβ

is obtained in a purely

diagonal form,

αβ

⎧

⎨

⎩



αβ

, for γ

,γ

≡ η,



αβ

, for γ

≡ A

,γ

≡ A

(3.3.19)

We represent here c(k) as the spatial Fourier transform of the direct correlation c(r), which

is assumed to be translationally invariant for the uniform liquid state. We have used here

the notation c

= n

c(K

) for the liquid density n

and the double-primed quantity, i.e.,



, is the corresponding second derivative of c

with respect to the wave vector K

.Itis

useful to note the following relations in this respect:



dr c(r)rr =−n



c(K )



K =0

=−c



I, (3.3.20)



dr c(r)e

· r

rr =−n



c(K )



=−



, (3.3.21)

where

is the unit vector along K

. The subscript 0 in c

denotes the same quantity

for zero wave vector. Note that with the present choice (3.3.3) for the density function the

kernel #

αβ

is not only diagonal but also determined solely in terms of the static correlation

of the uniform liquid state.

The total difference of the grand potentials between the inhomogeneous and homoge-

neous states is obtained by evaluating the RHS of (3.3.13) in this case as

β= β







|∇η(r)|





{

· ∇A

(r)}



. (3.3.22)

3.3 The density-functional approach 143

Using the above expression for  in (3.3.15), a set of ordinary differential equations for

the order parameters {η, A

} is obtained. Differentiating (3.3.22) with respect to η and A

respectively, gives

−1





1 + η +



· x



= c

η −βU

, (3.3.23)

−1



· x



η +



· x



= c

− βU

. (3.3.24)

Note that the LHS of each of the above equations follows from the corresponding results

for the uniform system with local order parameters {η, A

}. We have used for simpliﬁcation

in the RHS of (3.3.23) the following deﬁnitions for U

and U

βU

(x) =



∇

η(x), βU

(x) =



(

· ∇)

(x). (3.3.25)

and U

are the Fourier components of an external ﬁeld U (x) that maintains the density

proﬁle in a uniform system with the space-dependent order-parameter ﬁelds {η(x), A

βU (x) = βU



βU

· x

(3.3.26)



∇

η(x) +



· x



(

· ∇)

(x). (3.3.27)

For the special case of the spherical nucleus considered here the order parameters η and

are merely radial functions and the deﬁning relations for U

and U

in (3.3.25) are

expressed with the following differential equations:

βU

(r) =







η(r ), (3.3.28)

βU

(r) =







(r). (3.3.29)

In (3.3.28) we have averaged over the angular dependence present in (3.3.26) due to

the directional operator (

· ∇)

. These two equations are coupled to the two equations

(3.3.23) and (3.3.24) obtained above. In order to solve for the optimum density distribution

in the critical nucleus corresponding to the saddle point in the space of the density function,

the last two equations are expressed as a single relation,

1 + η +



· x

= c

η −βU



· x

− βU

}. (3.3.30)

Equivalently, we have,

1 + η +



· x

= e

η−βU

exp





− βU

· x



. (3.3.31)

144 Crystal nucleation

Hence the order parameters η and A

are obtained in terms of the self-consistent equations

1 + η =

η−βU



exp





− βU

· x



(3.3.32)

and

η−βU



· x

exp





− βU

· x



. (3.3.33)

The coupled set of equations (3.3.26), (3.3.27), (3.3.32), and (3.3.33) constitutes the den-

sity distribution for the critical nucleus. To solve it the required input is the structure factor

for the uniform liquid.

The excess free energy per particle of the critical nucleus is then obtained from (3.3.1)

by evaluating  for the optimum density distribution. We show in Appendix A3.2 that

 is obtained as



= 4π



∞



[η, A

]−





dη



−











dr.

(3.3.34)

[η(r ), A

(r)] represents the local contribution to the excess free energy. The latter is

obtained from the corresponding functional for a uniform system characterized by {η, A

ThelasttwotermsontheRHSof(3.3.34) represent the nonlocal contribution due to the

variation of the order-parameter ﬁeld.

In the ﬁrst application of the density-functional method to the nucleation problem

Harrowell and Oxtoby (1984) solved the model equations for a simpliﬁed case by set-

ting c



= 0 and considered only a single RLV, K

, corresponding to the peak of the liquid

structure factor. It follows from eqns. (3.3.28) and (3.3.29), respectively, that βU

= 0 and

βU

(r) =







(r). (3.3.35)

η is now an algebraic function of A

as follows from (3.3.32). D

exhibits two minima

as a function of the order parameter A

, corresponding to the crystalline and uniform

liquid states, respectively. The liquid minimum exists down to the lowest supercooling

temperature, indicating the absence of a spinodal point in the present model. From the r

dependences on the order parameters η and A

for the critical nucleus corresponding to

various temperatures below the freezing point it follows that the crystal–liquid interface

is a few atomic diameters thick and that this is essentially independent of supercooling.

With lowering of the temperature the order parameter A

remains essentially constant near

r =0, being equal to the corresponding value for the solid state. This indicates that the

center of the nucleus has properties characteristic of the bulk solid. On the other hand,

the order parameter η in the same region steadily decreases with supercooling. Thus for

the critical nucleus the structure and average density, characterized by A

and η, respec-

tively, do not change at the same point in the interface. From the liquid side the crystalline

3.3 The density-functional approach 145

Fig. 3.6 The free-energy barrier  A

∗

(in units of k

T ) to nucleation vs. the extent of undercooling

T , from the DFT (solid) and CNT (dashed) with the capillary approximation. From Harrowell and

Oxtoby (1984).

 American Institute of Physics.

structure appears ﬁrst in the interface region and the density changes more inside the

nucleus. It is energetically less favorable to compress a disordered liquid and once the

order appears it is easier to compress. The free energy of the critical nucleus is computed

from the formula (3.3.34) in the present calculation. The result is shown in Fig. 3.6.The

corresponding result according to the CNT is estimated from the formula (3.1.17).Here

G

(noted as A in Fig. 3.6) was computed using the DFT model. For the surface free

energy γ

the DFT results for the planar interface were used. Figure 3.6 shows that the

dependence of the free energy on the supercooling as predicted in the CNT is qualitatively

the same as that in the present DFT model, though the former somewhat underestimates

the free energy of the critical nucleus. In the CNT the radius of the corresponding critical

nucleus is estimated from (3.1.18). The DFT approach by construction does not assume

any sharp interface. However, a radius for the critical nucleus can be deﬁned (Harrowell

and Oxtoby, 1984) as the distance at which A

falls to half of its value at the center of the

nucleus. Such an estimate indicates that the CNT somewhat underestimates the radius of

the critical nucleus in comparison with that predicted from the DFT.

3.3.3 The weighted-density-functional approach

The DFT model for the critical nucleus using a simple Ramakrishnan–Yussouff-type per-

turbation expansion for the free energy was subsequently improved using more advanced

techniques like the modiﬁed weighted-density-functional approximation (MWDA) by Shen

and Oxtoby (1996a, 1996b). We have already discussed the MWDA calculation of the

bulk free energy in the previous chapter. Here we will consider its speciﬁc application for

146 Crystal nucleation

identifying the critical nucleus. This involves the identiﬁcation of a set of order-parameter

ﬁelds {γ

(x)} characterizing the density function for the nucleus forming in the melt.

Let us ﬁrst identify a suitable set of parameter ﬁelds in deﬁning the density function. The

inhomogeneous density is expressed in terms of the usual RLV expansion in the following

form:

(x) =¯n

+ n



· x

, (3.3.36)

where ¯n

denotes the average density, n

is the average density of the solid and the μ

denote the Fourier coefﬁcients corresponding to the RLV K

. Note that the amplitude

is dependent only on the magnitude of K

. The average density ¯n

changes from the

density of the liquid (n

) to that of the solid (n

). The above expansion (3.3.36) in the RLVs

can also be identiﬁed with a simple parametrization of the inhomogeneous density as a

sum of the Gaussian density proﬁles centered at the lattice sites R

(x) =

¯n

)

3/2



}

−α(x−R

)

, (3.3.37)

where the inverse of α gives the width of the density proﬁle. a

is the lattice constant

with the conventional unit cell being in the shape of a cube of side 2a

and having four

particles in each cell for the f.c.c. lattice structure. 2a

is the volume of the primitive

unit cell such that the factor 2a

¯n

on the RHS of (3.3.37) represents the dependence

of the density function on the average density ¯n

. The two parametrizations (3.3.36) and

(3.3.37) are equivalent to each other (see the discussion accompanying eqns. (2.2.15) and

(2.2.16) in Chapter 2) with the choice μ

= exp

−K

/(4α)

for the Fourier amplitudes.

Therefore deﬁning the density function in terms of Gaussian density proﬁles is equivalent

to parametrizing all the RLV amplitudes μ

(m > 1) in (3.3.36) in terms of the amplitude

corresponding to the smallest nonzero K

= (μ

)

. (3.3.38)

Therefore the two order parameters ¯n

and μ

completely deﬁne the inhomogeneous den-

sity function. The grand potential D for the bulk system is a functional of the inhomoge-

neous density function characterized by a pair of values for the parameters {¯n

,μ

}.For

a Lennard-Jones system free energy is obtained as a sum of two different contributions.

The Lennard-Jones interaction potential (see eqn. (1.2.117) for its deﬁnition) is divided

into a short-range repulsive part and a long-range attractive part. The contribution to the

free energy from the repulsive part of the interaction is computed in terms of an equiva-

lent hard-sphere system using the MWDA as described in Chapter 2. (See eqns. (2.2.13)

and (2.2.20).) The contribution to the free energy from the attractive part of the interac-

tion is obtained in a mean-ﬁeld approximation involving the pair correlation function for

the uniform liquid state. For example, Shen and Oxtoby (1996a) obtained for a Lennard-

Jones system at T = 0.6 the stable bulk solid phase for ¯n

= 1.04 with the parameter

3.3 The density-functional approach 147

= 0.85 and the metastable bulk-liquid minimum for ¯n

= 0.89 and μ

= 0. All densi-

ties and temperatures for the Lennard-Jones system were expressed in dimensionless units

(see Chapter 1).

Now let us consider the technique being applied to the critical nucleus. The spherically

shaped critical nucleus nucleating from the melt is described in terms of the two order

parameters {¯n

,μ

}, which are treated as functions of the radial distance r from the cen-

ter. This form is similar to (2.4.32) used in the previous chapter for studying interfaces

using the weighted-density-functional model. The spatial dependence of the order param-

eters ﬁelds {¯n

,μ

} is controlled by the corresponding Euler–Lagrange equations (3.3.14)

obtained from the minimization of the grand potential. The optimum order-parameter

functions therefore follow from the solution of the simple differential equations given by

(3.3.15) as discussed above. The boundary conditions in this case are chosen in accord

with the requirement that at the center of the nucleus (r =0) and at large distance (r →∞)

the order parameters are equal to their characteristic values for the stable (crystalline)

and metastable (liquid) phases, respectively. The order-parameter proﬁles for the critical

nucleus obtained in the present case (Shen and Oxtoby, 1996a) are qualitatively sim-

ilar to those obtained in the model described above, with the Ramakrishnan–Yussouff

(Ramakrishnan and Yussouff, 1979) functional for the free energy. The results for the

optimum density proﬁle for the critical nucleus obtained with the MWDA are shown in

Fig. 3.7. The density proﬁle (¯n

) decays more rapidly across the interface than the structural

order parameter (μ

). However, the results obtained from the Ramakrishnan–Yussouf (RY)

free-energy functional and the MWDA approach follow opposite trends when compared

Fig. 3.7 The density (ρ) and structural order parameter (m), denoted by n

and μ

, respectively,

in the text, for the crystal nucleus at temperatures T =75 K (solid), T =70 K (long-dashed line),

T =60 K (short-dashed line), and T =50 K (dotted line). Reproduced from Shen and Oxtoby

(1996b).

 American Institute of Physics.

148 Crystal nucleation

Fig. 3.8 An illustration of Bain’s distortion. Two f.c.c. cells (each of size a) outlined with thin black

lines are shown. v, t,andu form an orthogonal triad of vectors. Along the t axis, the atoms are shown

with open circles on the t = 0 plane, dot-ﬁlled circles on the t = a/2 plane, and solid circles on the

t = a plane. With Bain’s distortion along the v axis, the cube outlined with thick dark lines becomes

b.c.c. From Shen and Oxtoby (1996a).

 American Physical Society.

with the corresponding CNT predictions. The free energy of the critical nucleus and its

size obtained with the CNT are less than the corresponding results from the DFT formu-

lated with the RY expansion. On the other hand, the CNT results are greater than the

corresponding quantities obtained in the MWDA model.

The nature of the critical nucleus is further analyzed by using the MWDA model (Shen

and Oxtoby, 1996b) to study the free-energy landscape for the system corresponding to

different possible crystalline structures. In this extension of the theory the parameter space

for the inhomogeneous density function includes a new order parameter, χ, which refers to

the symmetry of the crystal lattice. χ changes in a continuous manner monitoring Bain’s

distortion (Bain, 1924)(seeFig. 3.8) in the lattice structure, which is deﬁned in terms of

the following set of real-space vectors:



ˆκ

√

ˆκ



, a



ˆκ

√

ˆκ



, a

( ˆκ

+ˆκ

), (3.3.39)

where ˆκ

, ˆκ

, and ˆκ

represent an orthogonal triad of basic vectors (indicated in Fig. 3.8

as v, t, and u, respectively). The lattice symmetry corresponds to b.c.c. and f.c.c. for χ

taking the values 1 and

√

2, respectively. The RLVs corresponding to (3.3.39) are obtained

as linear combinations of the following set of basis vectors:

2π



−ˆκ

+ˆκ

√

ˆκ



, b

2π



ˆκ

−ˆκ

√

ˆκ



2π



ˆκ

+ˆκ

−

√

ˆκ



. (3.3.40)

3.3 The density-functional approach 149

Fig. 3.9 The landscape of the local grand-canonical free-energy density (ω) at coexistence tempera-

ture T =83.1Kand ¯n

= 0.89. The symbol m denotes the parameter μ

used in the text. The three

minima shown correspond to (a) the stable liquid with m =0, (b) the stable f.c.c. solid with χ =1.414

and m =0.818, and (c) the metastable b.c.c. solid with χ =1.0andm =0.750. The metastable b.c.c.

solid creates a saddle point near χ =1.0. From Shen and Oxtoby (1996a).

 American Physical

Society.

At ﬁxed values of the density ¯n

and temperature T , the landscape for the local grand-

canonical free-energy density F in the parameter space of χ and μ

exhibits (Shen and

Oxtoby, 1996b) a minimum for a stable f.c.c. structure



χ =

√



as well as one for a

metastable b.c.c. structure (χ =1). This landscape is displayed in Fig. 3.9 for ¯n

=0.89

and coexistence temperature T =83.1 K. In this extended description for studying the

critical nucleus, the order parameters which characterize the nucleus are given by the

set {¯n

,χ,μ

}. As the liquid is undercooled below the coexistence temperature, study of

the MWDA model obtains the proﬁle for χ(r ) corresponding to the critical nucleus. The

spatial dependence of χ represents a substantial b.c.c. character (χ =1) of the nucleus

near the interface changing over to f.c.c. structure



χ =

√



towards the center. With

increased supercooling the b.c.c. character spreads further inside, until at temperature 50 K

the nucleus is essentially all b.c.c. We note at this point that the computer studies of nucle-

ation of the crystal from a melt (ten Wolde et al., 1996) discussed next also suggest that in

the critical nucleus the interfacial region displays strong signatures of the b.c.c. symmetry

while near the center there is f.c.c. symmetry. This changes the surface free energy of the

nucleus. For a small-sized nucleus even the nature of the crystal symmetry near the center

is changed and is more b.c.c. like. Such a metastable b.c.c. phase has also been seen (Löser

et al., 1994) in experiments on rapidly cooled liquid metals whose stable crystalline phase

is f.c.c. An earlier theoretical study by Klein and Leyvraz (1986) using a ﬁeld-theoretic

model of the phase transition also suggested the possibility of a metastable b.c.c. phase in

the ﬁrst-order phase transition.

150 Crystal nucleation

3.4 Computer-simulation studies

The phenomenon of nucleation of a crystal from the melt has been studied using real exper-

iments as well as computer simulations. There are some obvious problems in making an

accurate study of the evolution of the nucleation process experimentally. It is difﬁcult to

study in an experiment the formation of the critical nucleus since the crystallite spends

only a microscopic time in this stage of its evolution as the top of the potential barrier

is approached. Computer simulation, which is a useful tool for studying the microscopic

description of the ﬂuid, has also been used for studying the nucleation process. In the

present section we discuss some of the recent developments in the study of the nucleation

process using such numerical methods. The straightforward approach of simulation for

this problem is to supercool the liquid below the freezing point and watch the formation

of the crystal nucleus. The free-energy barrier to the nucleation process falls as the inverse

square of the supercooling (see eqn. (3.1.17) for the nucleation barrier in the CNT) so that

at large supercooling the barrier is low. Mandell et al. (1976) simulated a small Lennard-

Jones system of 128 particles to study the nucleation process by monitoring the structure

factor. Honeycutt and Andersen (1984) investigated the effect of system size on nucleation.

Assuming the nucleation rate to be proportional to the volume, the observed nucleation rate

in the simulation is expected to grow with the number of particles. The observed nucleation

rate, however, decreased with the number of particles N, possibly indicating that the phe-

nomenon being observed was not homogeneous nucleation. The simulations were strongly

inﬂuenced by the periodic boundary conditions used in standard molecular-dynamic meth-

ods. Subsequent work by Swope and Andersen (1990) considered nucleation with one

million Lennard-Jones particles. The study shows that, although b.c.c. and f.c.c. crystalline

phases are formed in the early stage of nucleation, only the f.c.c. nuclei grow beyond the

post-critical stage. In all these cases the nucleation is studied at deep supercooling (more

than 40% below melting).

The free-energy barrier to critical-nucleus formation is not small enough unless the sys-

tem is deeply supercooled and as a consequence the rate of nucleation is too low to be

observed. For example, at 20% supercooling in a typical simulation of 10

or so particles

one has to wait up to 10

simulation steps to observe one nucleation event (ten Wolde

et al., 1996). Thus without deep supercooling the “brute-force” method for studying crys-

tallization with computer simulations fails. At least 50% supercooling is needed in order

to observe crystal formation over the time scale of a simulation. Studying the nucleation

of the crystal from a melt using computer simulations at low undercooling (close to the

freezing point) when the free-energy barrier is large requires the development of special

techniques. In this regard it is useful to identify suitable reaction coordinates characteriz-

ing the crystallinity of the system. The free energy for the nucleus is obtained from the

probability of occurrence of the corresponding value of the reaction coordinate, and this

probability is computed through efﬁcient sampling of the conﬁguration space. By locating

the maximum in the free-energy barrier against the reaction coordinate the critical nucleus

is identiﬁed.