- Email: [email protected]
A kinetic model for the premelting of a crystalline structure
ARTICLE IN PRESS
Physica A 387 (2008) 134–144 www.elsevier.com/locate/physa
A kinetic model for the premelting of a crystalline structure Yuri Djikaev, Eli Ruckenstein Department of Chemical and Biological Engineering, SUNY at Buffalo, Buffalo, NY 14260, USA Received 18 May 2007 Available online 19 August 2007
Abstract An analytical kinetic approach to examine the premelting phenomenon is suggested by using a first passage time analysis. Premelting is considered to occur when the time of formation of a Frenkel type defect in the surface monolayer becomes sufficiently small. The mean time of defect formation on the surface lattice, i.e., the mean time necessary for a selected (surface-located) molecule to leave its lattice site and form a Frenkel defect, is calculated by using a first passage time analysis. The model is illustrated by numerical calculations for a crystalline structure composed of molecules interacting via the Lennard-Jones (LJ) potential. The lattice vectors in the plane parallel to the free surface of the crystal were assumed to be equal (to the lattice parameter) and the angle between them was varied. The model predictions of the Tammann temperature (of premelting) are very sensitive to the parameters of the LJ potential. In all the cases considered, the temperature dependence of the mean first passage time has two clearly distinct regimes: at low temperatures the dependence is sharp and at high temperatures it is weak. r 2007 Elsevier B.V. All rights reserved. PACS: 05.10.Gg; 05.70.Fh; 05.70.Np Keywords: Premelting; Tammann temperature; First passage time analysis
1. Introduction Although some theories of melting consider the process as a breakdown of the crystal lattice which occurs uniformly and homogeneously throughout the solid at the melting point [1], there is evidences that melting is a heterogeneous process involving the nucleation of a liquid phase at some preferred sites of the solid (the free surface, grain boundaries, large dislocations and disclinations, etc.) and subsequent growth of the liquid phase. The free surface of a crystalline structure constitutes a natural ‘‘heterogeneous’’ site for the onset of the crystal-to-liquid phase transition (as it does for the liquid-to-crystal transition [2]). Tammann [3] was the first to point out that the surface premelts before melting occurs. This idea is consistent with the empirical criterion proposed by Lindemann [1] for melting in the bulk according to which melting might be expected to occur when the root mean amplitude of thermal vibrations of an atom exceeds a certain threshold value of more or about 10% of the distance to the nearest neighbor in the crystalline Corresponding author. Tel.: +1 716 645 2911x2214; fax: +1 716 645 3822.
E-mail addresses: [email protected] (Y. Djikaev), [email protected] (E. Ruckenstein). 0378-4371/$ - see front matter r 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.physa.2007.08.022
ARTICLE IN PRESS Y. Djikaev, E. Ruckenstein / Physica A 387 (2008) 134–144
135
structure. Tammann suggested [3] that the outermost layer of the crystal should acquire some disorder far below the bulk melting point. This behavior is a result of the higher freedom of motion for the surface-located molecules. The surface melting (often referred to as premelting) involves the formation of a thin layer which acquires some disorder at a temperature significantly below the melting one. The surface premelting can provide a plausible explanation to why a liquid can be cooled below the crystallization point (i.e., supercooled), whereas a solid cannot be heated above the melting point (i.e., superheated). Such a hysteresis suggests that there is no free-energetic barrier for the melting phase transition which occurs as spinodal decomposition (a barrier exists for homogeneous crystallization which thus always occurs via nucleation). As a possible cause, one can suggest that the surface of a solid, which premelts at temperatures far below the melting point, constitutes a large natural ‘‘heterogeneous’’ center ‘‘whereupon’’ the melting phase transition occurs in a barrierless way. Experimentally the premelting phenomenon was apparently first detected by Lyon and Somorjai [4] who studied the structures of the clean (1 1 1), (1 1 0), and (1 0 0) crystal faces of platinum as a function of temperature by means of low-energy electron diffraction and observed the formation of disordered surface structures at temperatures far below the melting temperature of T m ¼ 2043 K. Direct experiments on the surface initiated melting were also carried out by Frenken et al. [5] using Rutherford backscattering in conjunction with ion-shadowing and blocking. That experiment revealed a reversible order–disorder transition on the (1 1 0) surface of a lead crystal well below its melting point of T m ¼ 600:7 K. Since then, other techniques have been employed such as calorimetry, electron, neutron, and X-ray diffraction, microscopy, ellipsometry, and helium scattering. Although most experiments were carried out under equilibrium conditions, melting tended to be initiated at the surface even when the crystalline solid was heated very quickly so that equilibrium conditions were not established [6]. With the advent of superpowerful computers, molecular dynamics (MD) simulation methods have become an invaluable and major tool in studying the premelting phenomenon at a microscopic level [7–9]. Among recent applications of these methods one can mention the simulation of premelting in AgBr [10], premelting in Cr2O3 [11], the premelting of a clean Al(1 1 0) surface, [12], etc. Still, theoretical models of premelting have proven to be useful in elucidating various aspects of this phenomenon (some examples are given in Ref. [13]). In the thermodynamic approach premelting is regarded as a particular case of wetting, namely the wetting of a solid by its own melt [14]. Another phenomenological model of premelting was developed [15–17] in the framework of the Landau–Ginsburg theory wherein the free energy of the system is expanded in a power series with respect to an order parameter whereof the temperature dependence is then analyzed. A microscopic theory of melting of two-dimensional solids on the basis of the dislocations pairs model was developed in Ref. [18] and discussed in Ref. [13]. Another microscopic theory of surface melting [19] was based on a lattice model of a solid and explored the fact that by using a discrete reference lattice and a mean-field approach, the partition function of a system of particles interacting via a pairwise potential could be calculated. An efficient theoretical approach [20] to study premelting was developed on the basis of the density functional theory. In that approach the grand canonical free energy of an inhomogeneous system (for given chemical potential(s), volume, and temperature) is regarded as a functional of the local density, the equilibrium spatial profile of which minimizes the free energy functional. As an attempt to clarify some kinetic aspects of premelting, in this paper we present a microscopic model of the process on the basis of a first passage time analysis. The model assumes that the thermal motion of a molecule located in the surface layer of a finite crystalline system is determined only by its interactions with other surface-located (in the lattice sites) molecules. The main idea of the model is to determine the mean time of a defect formation on the lattice surface, i.e., the mean time necessary for a selected (surface-located) molecule to leave its lattice site and go wandering on the surface. The temperature dependence of this time exhibits two clearly distinct regimes, the high- and low-temperature ones, with a relatively non-sharp transition between them. The onset temperature of the high-temperature regime can be roughly associated with the Tammann temperature of premelting. It should be noted that, unlike the boiling or melting transitions whereof the temperatures can be unambiguously defined, there exists no strict definition of the Tammann temperature. The transition of the surface layer of a crystalline body from an orderer state into a
ARTICLE IN PRESS Y. Djikaev, E. Ruckenstein / Physica A 387 (2008) 134–144
136
premelted one is apparently a gradual process and one can associate the Tammann temperature to that which corresponds to one or another degree of ‘‘surface’’ disordering. The paper is structured as follows. In Section 2 we present a model for a semi-infinite crystalline structure with a surface layer. The basics of the mean first passage time analysis and its application to premelting are given in Section 3, while Section 4 contains numerical calculations based on our model. The results are discussed and the conclusions are summarized in Section 5. 2. Semi-infinite crystalline system of pairwise interacting molecules Consider a semi-infinite crystalline system such that in the Cartesian coordinates it will occupy the space 1oxo1; 1oyo1; 1ozp0 so that the plane z ¼ 0 constitutes a free surface of the crystalline body. This crystalline structure is defined by three lattice vectors a, b, and c with a and b in the plane parallel to the plane x–y of the Cartesian coordinate system. The angle between a and b will be denoted Y. The crystal lattice sites are occupied by units hereafter assumed to interact via the Lennard-Jones (LJ) pair potential and henceforth termed ‘‘molecules’’. In reality, they can be atoms, actual molecules, groups of atoms or molecules, colloidal particles, etc. Important is the pairwise and central character of their interactions. Consider a molecule at a lattice site with the coordinates x ¼ 0; y ¼ 0; z ¼ 0 (i.e., in the origin of the coordinate system) and denote it by s (this is one of the molecules belonging to the surface layer of the crystalline structure). Clearly, the molecule s is not frozen at its lattice site but performs constant thermal motion of (nonlinear) vibrational type. The x–y motion of the molecule s is completely determined by the x and y components of the forces acting thereupon. Assuming a continuum distribution of molecules below the surface layer, it is clear that F sxy , the xy component of the total force acting on the molecule s, is determined only by its interactions with the other surface-located molecules. Indeed, due to the continuum distribution of molecules in the interior of the crystal, for any elementary volume thereof there exists a counterpart located symmetrically (with respect to the axis z) such that the xy components of the forces exerted on the molecule s by these two elements are compensated. As a result, the force F sxy can be represented as F sxy ¼
N 1 X
f is ,
(1)
i¼1
where f is is a force exerted by the molecule associated with the lattice site i upon the selected molecule s and N is the total number of surface located molecules (including the selected one). Except for the molecule s, all other surface located molecules are assumed to be ‘‘frozen’’ at the corresponding lattice sites, so that for a given crystalline structure F sxy is a function of the x and y coordinates of molecule s only: F sxy F sxy ðxs ; ys Þ. It should be noted that, in reality, molecule s performsPmotion not only along the x- and y-axis, but also along the z-axis. Therefore, strictly speaking, the force N1 and its x and y i¼1 f is has also a z component P components depend on the coordinate z. However, it is clear that the z component of N1 i¼1 f is does not affect the x–y motion of molecule s. Besides, the temperature involved in premelting is much lower than the sublimation temperature. For this reason the motion in the z direction will be neglected. Let us assume that the interaction potential is pairwise and denote the potential between molecules i and s by fis ðx; yÞ.PThe overall potential of interactions of molecule s with all other surface located molecules is then js ðx; yÞ ¼ N1 i¼1 fis ðx; yÞ. Fig. 1 presents the function js ¼ js ðx; yÞ for the potential fis ðx; yÞ of the LJ type " 6 # Z 12 Z fis ðx; yÞ ¼ 4e , (2) ris ris where e and Z are the standard parameters of the LJ potential. Fig. 1a shows the contour plot of a threedimensional landscape determined by js ¼ js ðx; yÞ, whereas in Figs. 1b–d various sections of this landscape are plotted, corresponding to x ¼ 0 (similar to y ¼ 0), y ¼ x (similar to y ¼ x), and x ¼ 0:5 (similar to y ¼ 0:5), respectively (all the distances are given in units of the lattice parameter L, the potential is given in units of kB T, and the LJ parameters are e=kB T ¼ 0:5, Z=L ¼ 0:5). Clearly, the force F sxy and the potential
ARTICLE IN PRESS Y. Djikaev, E. Ruckenstein / Physica A 387 (2008) 134–144
137
Fig. 1. (a) Contour plot and various profiles of the potential field for a selected molecule (molecule s) associated with the lattice site ð0; 0Þ (for a lattice with a ¼ bð¼ LÞ; Z=L ¼ 0:5; Y ¼ p=2, and e=kB T ¼ 0:5). The darker regions correspond to a lower potential C and the lattice sites are indicated as points. (b)–(d) The profiles of the landscape (determined by Cðx; yÞ) along its various sections (indicated in the figure panels).
js ðx; yÞ are related through F sxy ¼ =js ðx; yÞ.
(3)
The potential field wherein the molecule s performs its thermal motion is a complicated array of potential wells with infinite barriers between them in some directions and finite ones in others. Two important features of this potential field should be noted. First, the site s of the lattice is surrounded by a series of local wells in each of which the molecule s can be located with equal probability. The selected molecule s, considered to be ‘‘at site s’’, would actually be located in one of these wells or would be in the process of transition from one of them into another, so that at any given moment the molecule s can actually be found (with various probabilities) in any place in the vicinity of the lattice site to which it is attributed. Second, there are eight pathways whereby the selected molecule s can leave the vicinity of its native lattice site, cross over a finiteheight barrier, and fall into another local well in the potential field. This would correspond to the formation of a surface version of a Frenkel type defect of the crystalline structure [21]. The probability (or the frequency) of such transitions depends on Z; e; L, and temperature T (hereafter expressed in K). Once these transitions
ARTICLE IN PRESS 138
Y. Djikaev, E. Ruckenstein / Physica A 387 (2008) 134–144
become very frequent, the surface layer of the lattice can be considered to have lost its ordered structure. Thus, for a given crystalline system (i.e., given Z; L and e) one can expect to detect premelting with increasing temperature T by examining the change of the mean first passage time, necessary for the molecule s to leave the vicinity of its native lattice site, with variations in the temperature T. In the following section, we outline the basics of a mean first passage time analysis and its applications to the problem of premelting. 3. A mean first passage time analysis Upon close examination of the landscape in Figs. 1a–d, one can conclude that the molecule s can leave the vicinity of its original site by passing over a barrier having the shape of a hyperbolic paraboloid (‘‘saddle’’) with the path of steepest descent parallel either to the x- or y-axis (Fig. 2a). Although the barrier is two-dimensional, it is very narrow over the saddle point hence one can assume that the molecule crosses the barrier by means of one-dimensional motion either along the x- or y-axis. One can therefore formulate a mean first passage time analysis for the chaotic motion of a molecule in a one-dimensional potential well. Let us consider the thermal chaotic motion of a molecule in a one-dimensional potential field CðxÞ (assumed to be expressed in kB T units) as shown in Fig. 2b. The field has the shape of a well with one of the boundaries (at x ¼ 0 in Fig. 2) impenetrable to the molecule and the other penetrable (at x ¼ l) so that the molecule can be regarded to have escaped from the well once it reaches this boundary. The mean first passage time of a molecule escaping from the well is calculated on the basis of a kinetic equation governing the chaotic motion of the molecule in this potential well. In general, the chaotic motion of a molecule subjected to an external force is governed by the Fokker–Planck equation for the single-molecule distribution function with respect to its coordinates and momenta, i.e., in the phase space [22–24]. Assuming that the relaxation time for the velocity distribution function is small compared to the characteristic time scale of the escape process [22], the Fokker–Planck equation reduces to the Smoluchowski equation [22–24] for the diffusion of a particle in an external force field which can be written in the form (see, e.g., Eq. (312) in Ref. [22]) qpðr; tjr0 Þ F ¼ = D =pðr; tjr0 Þ pðr; tjr0 Þ , (4) qt kB T where pðr; tjr0 Þ is the probability of observing a molecule in the volume element dr at r at time t given that initially it was at r0 , D is the diffusion coefficient, and F FðrÞ is the force exerted on the molecule by an
Fig. 2. (a) The vicinity of the point ð1; 0:5Þ for the potential field shown in Fig. 1a. In this vicinity the potential has a shape of a hyperbolic paraboloid with the path of steepest descent parallel to the x-axis. (b) One-dimensional potential well CðxÞ wherein a molecule is assumed to perform thermal chaotic motion before leaving the vicinity of its original lattice site. The well boundary at x ¼ 0 is assumed to be of infinite height, whereas at x ¼ l it is finite.
ARTICLE IN PRESS Y. Djikaev, E. Ruckenstein / Physica A 387 (2008) 134–144
139
external field. If the latter can be derived from some potential fðrÞ as F ¼ =fðrÞ, Eq. (4) can be rewritten as qpðr; tjr0 Þ ¼ =½DefðrÞ=kB T =ðefðrÞ=kB T pðr; tjr0 ÞÞ. (5) qt For the one-dimensional potential field CðxÞ ¼ fðxÞ=kB T and constant diffusion coefficient D, the Smoluchowski equation (5) reduces to qpðx; tjx0 Þ q CðxÞ q CðxÞ ¼D e ðe pðx; tjx0 ÞÞ . (6) qt qx qx The mean minimum time, necessary for a molecule in the well to cross over the external boundary of the well and leave it, depends on the initial position x0 of the molecule. It is convenient to use the backward Smoluchowski equation [23,25–28] which expresses the dependence of the transition probability pðx; tjx0 Þ on x0 : qpðx; tjx0 Þ q q ¼ DeCðx0 Þ eCðx0 Þ pðx; tjx0 Þ . (7) qt qx0 qx0 The probability that a molecule, initially at a distance x0 from the inner boundary of the well, will remain in the well after time t is given by the so-called survival probability, [25–28] Z l Qðtjx0 Þ ¼ dxpðx; tjx0 Þ ð0ox0 olÞ, (8) 0
where l is the width of the potential well, i.e., the distance between the inner and outer boundaries of the well (see Fig. 2b). The probability for the dissociation time to be between 0 and t is equal to 1 Qðtjx0 Þ, and the probability density for the dissociation time to be between 0 and t is given by qQ=qt. The first passage time is provided by [25–28] Z 1 Z 1 qQðtjx0 Þ dt ¼ tðx0 Þ ¼ t Qðtjx0 Þ dt, (9) qt 0 0 where the rightmost equality is obtained by integrating the integral on its LHS by parts. The equation for the first passage time is obtained by integrating the backward Smoluchowski equation (7) with respect to x and t over the entire range and using the boundary conditions Qð0jx0 Þ ¼ 1 and Qðtjx0 Þ ! 0 as t ! 1 for any x0 (the latter means that if the molecule is initially in the well, sooner or later it will leave the well, whereas the former trivially says that the molecule, initially in the well, cannot jump out of it instantaneously). Thus, one obtains Cðx0 Þ d Cðx0 Þ d De e tðx0 Þ ¼ 1. (10) dx0 dx0 One can solve Eq. (10) by assuming a reflecting boundary condition dt=dx0 ¼ 0 at x0 ¼ 0 and a complete absorption boundary condition tðx0 Þ ¼ 0 at x0 ¼ l. Consequently, tðx0 Þ is given by Z Z y 1 l dy eCðyÞ dx eCðxÞ . (11) tðx0 Þ ¼ D x0 R The average dissociation time, t¯ (often referred to as the mean first passage time), is obtained by averaging tðx0 Þ with the Boltzmann factor over all possible initial positions x0 : Z 1 l t¯ ¼ dx0 eCðx0 Þ tðx0 Þ, Z 0 Z Z¼
l
dx0 eCðx0 Þ .
(12)
0
Note that the mean first passage is inversely proportional to the diffusion coefficient D of a molecule in the potential well. Although this quantity may not be readily available for calculating the exact value of t¯ , even the
ARTICLE IN PRESS 140
Y. Djikaev, E. Ruckenstein / Physica A 387 (2008) 134–144
dimensionless quantity t¯ D=Z2 as a function of kB T=e can provide some information on the premelting behavior of the system because it relates the average time of formation of a surface defect of Frenkel type to the temperature. A drastic change in this behavior would indicate a qualitative change in the state of the surface layer of a crystalline structure with changing temperature. 4. Numerical evaluations As an illustration of the above model (described in Section 2), we carried out numerical calculations for crystalline lattices consisting of molecules interacting via the LJ potential, Eq. (2). The lattice vectors in the plane x–y (parallel to the surface) were assumed to be of equal length, a ¼ bð¼ LÞ, and lattices with various Z=L and Y (angle between the vectors a and b) were considered. At fixed Z=L and Y the main quantity of interest is the dimensionless quantity t¯ D=Z2 (representative of the mean first passage time necessary for a selected molecule s to leave the vicinity of its original lattice site by overcoming a barrier) as a function of the dimensionless temperature kB T=e. The results of numerical calculations are shown in Figs. 1, 3–5. In Fig. 1a (representing the case a ¼ bð¼ LÞ; Z=L ¼ 0:5; Y ¼ p=2, and kB T=e ¼ 2) the darker regions correspond to a lower potential C and the lattice sites, indicated as points, lie in the centers of the white regions. As clear from this figure, the selected molecule s can leave the vicinity of its original lattice site by crossing the barrier between one of the eight local potential wells and another of the wells not belonging to the site s. Although both potential wells are two-dimensional and the barrier between them has a ‘‘saddle’’-like shape, the bottleneck of the barrier is very narrow hence the mean first passage time for this transition can be estimated by considering the transition to occur via one-dimensional motion of the molecule along the path of the steepest descent from one well into the other. Eq. (12) can then be applied to evaluate the dimensionless ratio t¯ D=L2 . Figs. 1b–d present the profiles of the landscape (determined by Cðx; yÞ) along its various sections. They show that molecule s cannot leave the vicinity of its original lattice site (with coordinates x ¼ 0; y ¼ 0) along x ¼ 0, y ¼ 0, y ¼ x, and y ¼ x, but only along x ¼ 0:5; x ¼ 0:5; y ¼ 0:5, and y ¼ 0:5. Fig. 3a shows the dependence of the dimensionless mean first passage time t¯ D=L2 on the dimensionless temperature kB T=e for the case a ¼ bð¼ LÞ; Y ¼ p=2. Each curve corresponds to a fixed ratio Z=L (as specified in the caption of Fig. 3a). The examination of the slopes of these curves indicates that there are two distinct regimes for the variation of the mean first passage time with changes in temperature (e assumed fixed). At high temperatures (large kB T=e) the mean first passage time is hardly sensitive thereto, whereas at low temperatures (small kB T=e) the mean first passage time sharply changes with changing temperature. One can conjecture that these two regimes correspond to two different states of the surface layer of the crystalline structure, associating the low-temperature state with an ordered structure and the high-temperature state with a disordered structure. The transition between these two states is not very sharp but occurs rather smoothly with changing
Fig. 3. (a) Dimensionless mean first passage time t¯ D=L2 vs the dimensionless temperature kB T=e for a ¼ bð¼ LÞ; Y ¼ p=2. Each curve corresponds to a fixed ratio Z=L (0.53, 0.52, 0.51, 0.5, 0.49, 048, 0.47 from right to left). (b) The dimensionless Tammann temperature kB T m2 =e as a function of Z=L for Y ¼ p=2.
ARTICLE IN PRESS Y. Djikaev, E. Ruckenstein / Physica A 387 (2008) 134–144
141
Fig. 4. (a) Dimensionless mean first passage time t¯ D=L2 as a function of the dimensionless temperature kB T=e. Each curve corresponds to a fixed Y (0:99p=2; 0:98p=2; 0:97p=2; 0:96p=2; 0:95p=2; 0:94p=2; 0:93p=2 from right to left). (b) The dependence of the (dimensionless) Tammann temperature kB T m2 =e on Y for Z=L ¼ 0:5.
temperature. The onset of the high-temperature regime can probably be regarded as a rough estimate of the Tammann temperature at which the premelting of a crystalline system is said to occur. The uncertainty of a thus determined Tammann temperature is consistent with the somewhat arbitrary and non-rigorous character of the original definition of the Tammann temperature. In Fig. 3b the dimensionless Tammann temperature kB T m2 =e (determined as the temperature at which t¯ is by an order of magnitude higher than its value in a hightemperature regime) is plotted as a function of Z=L. One can note a decrease in kB T m2 =e with decreasing Z=L. The LJ potential in the soft sticky dipole (SSD) model of water (in which a water molecule is considered to have only a single, center of mass, interaction site with a tetrahedrally coordinated sticky potential that regulates the tetrahedral coordination of neighboring molecules) are e ¼ 15:3 kJ=mol and Z ¼ 3:01 A˚ [29]. For the cubic ice L ¼ 6:36 A˚ [30]. Thus, in this model one can estimate the Tammann temperature for the cubic ice as T m2 ’ 185 K which is in agreement with the experimentally observed temperature of about 200 K for the premelting of the (0 0 0 1) face of hexagonal ice [31]. To examine the sensitivity of the above model to the angle Y, we calculated the dimensionless mean first passage time t¯ D=L2 as a function of temperature for several values of Y (with a ¼ bð¼ LÞ and Z=L ¼ 0:5). The results are presented in Fig. 4. In Fig. 4a each curve corresponds to a Y specified in the figure caption. One can conclude that the dimensionless mean first passage time t¯ D=L2 is sensitive to Y, and so is the Tammann temperature. This can be explained by the drastic deformation that the landscape of the potential field Cðx; yÞ suffers with even a small change in Y. The smaller is Y at a given Z=L, the higher is the surface number density of molecules (lattice sites on a free surface of the crystal structure). Hence, the premelting of the surface layer occurs easier (the Tammann temperature decreases) with increasing surface density of lattice sites (i.e., with decreasing Y). This is clearly seen in Fig. 4b which presents the dependence of the (dimensionless) Tammann temperature kB T m2 =e on Y for Z=L ¼ 0:5. The above results (Figs. 1,3,4) describe the process of initial formation of a Frenkel defect in a defectless surface layer of the lattice. In order to somewhat clarify how the presence of a defect (i.e., a vacancy) in the structure of the surface layer affects the possible pathways of the evolution of its neighbors, we considered also a lattice with a pre-existing surface Frenkel defect (Fig. 5), where a molecule s0 has left the vicinity of its original lattice site ðx; yÞ ¼ ð0; 0Þ and now occupies the position ðx; yÞ ¼ ð0:65; 1:35Þ (all distances are expressed in units of L) marked as a point in a circle. In this case, the potential field for the molecule s1 associated with the lattice site ð0; 1Þ (one of the nearest neighbors of the vacancy at ð0; 0Þ) is shown in Figs. 5a (contour plot of the landscape with the darker regions representing lower values of C and the lattice sites indicated as points in the centers of white regions) and 5b–h. As expected, the potential field for molecule s1 has become very asymmetric. As shown in Figs. 5d–g, six pathways (whereby molecule s1 can leave the vicinity of its original lattice site to become a Frenkel type defect) are not affected by the existence of the defect, whereas two pathways become more probable because of decreases in the potential barriers.
ARTICLE IN PRESS 142
Y. Djikaev, E. Ruckenstein / Physica A 387 (2008) 134–144
Fig. 5. (a) Contour plot and various profiles of the potential field for a selected molecule (molecule s1 ) associated with the lattice site ð0; 1Þ in a lattice with a pre-existing surface Frenkel defect, where a molecule s0 has left the vicinity of its original lattice site ðx; yÞ ¼ ð0; 0Þ and now occupies the position ðx; yÞ ¼ ð0:65; 1:35Þ marked as a point (for a lattice with a ¼ bð¼ LÞ; Z=L ¼ 0:5; Y ¼ p=2, and e=kB T ¼ 0:5). The darker regions correspond to a lower potential C and the lattice sites are indicated as points. (b)–(h) The profiles of the landscape (determined by Cðx; yÞ) along its various sections as indicated in the figure panels.
ARTICLE IN PRESS Y. Djikaev, E. Ruckenstein / Physica A 387 (2008) 134–144
143
5. Concluding remarks It was demonstrated both experimentally and by simulations that the surface melting (often referred to as premelting) involves the formation of a thin disordered layer at a temperature significantly below the melting one. Several theoretical models of premelting were developed to elucidate various aspects thereof. All previous models have focused on obtaining the temperature of the premelting-like transition at the crystal free surface (i.e., the Tammann temperature) by using ‘‘static’’ formalisms (classical thermodynamic, Landau–Ginsburg, lattice, or density functional theories). In this paper, we have proposed an analytical kinetic approach to premelting on the basis of a first passage time analysis. The model assumes that the thermal motion of a molecule located in the surface layer of a finite crystalline system is determined only by its interactions with other surface-located (in the lattice sites) molecules. The main idea of the model is to determine the mean time of the defect formation on the surface lattice, i.e., the mean time necessary for a selected (surface-located) molecule to leave its lattice site and form a Frenkel defect. Thus, the model is aimed at predicting conditions under which the structure of the surface layer of a crystalline system may acquire some disorder due to a sharp decrease in the time of formation of Frenkel type defects. As an illustration, we carried out numerical calculations for a crystalline structure composed of molecules (base units) interacting via the LJ potential. The lattice vectors in the plane parallel to the free surface of the crystal were assumed to be equal (to the lattice parameter) and the angle between them was varied. In all cases considered, the temperature dependence of the mean first passage time has two clearly distinct regimes. At low temperatures it is strong, whereas at high temperatures it is weak. The onset of the high-temperature regime provides an estimate of the Tammann temperature. It should be noted that an important approximation has been made in evaluating the mean first passage time necessary for a surface molecule to leave the vicinity of its original lattice site to generate a Frenkel defect. In this vicinity there are eight local potential wells wherein the molecule is located most of the time. These potential wells are two-dimensional. The molecule leaves its original lattice site vicinity by crossing a barrier (between this vicinity and the rest of the lattice) having the shape of a hyperbolic paraboloid. To calculate the mean first passage time of this process, we assumed the motion of the molecule in the local well to be onedimensional and to occur along the path of the steepest descent. This approximation can be avoided by starting with a two-dimensional Smoluchowski equation (describing the chaotic motion of a molecule in a two-dimensional force field) and deriving the expression for the mean first passage time. References [1] F. Lindemann, Z. Phys. 11 (1910) 609. [2] Y.S. Djikaev, et al., J. Phys. Chem. A 106 (2002) 10247–10253; Y.S. Djikaev, et al., J. Chem. Phys. 118 (2003) 6572–6581. [3] Tammann, Z. Phys. Chem. Stoechiom. Verwandtschalft 68 (1910) 205. [4] H.B. Lyon, G.A. Somorjai, J. Chem. Phys. 46 (1967) 2539. [5] J.W.M. Frenken, P.M.J. Mare´e, J.F. van der Veen, Phys. Rev. B 34 (1986) 7506. [6] H. Hakkinen, U. Landmann, Phys. Rev. Lett. 71 (1993) 1023. [7] O. Tomagnini, F. Ercolessi, S. Iarlori, F.D. Di Tolla, E. Tosatti, Phys. Rev. Lett. 76 (1996) 1118. [8] F.D. Tolla, F. Ercolessi, S. Iarlori, Surf. Sci. 211/212 (1989) 55. [9] S. Iarlori, P. Carnevali, F. Ercolessi, E. Tosatti, Surf. Sci. 211/212 (1989) 75. [10] A.K. Ivanov-Schitza, G.N. Mazob, E.S. Povolotskaya, S.N. Savvin, Solid State Ionics 173 (2004) 103–105. [11] M.A.S.M. Barrera, J.F. Sanz, L.J. A´lvarez, J.A. Odriozola, Phys. Rev. B 58 (1998) 6057–6062. [12] R. Zivieri, G. Santoro, V. Bortolani, Phys. Rev. B 62 (2000) 9985–9988. [13] E. Ruckenstein, in: S.A. Stevenson, J.A. Dumesic, R.T.K. Baker, E. Ruckenstein (Eds.), Metal-support Interactions in Catalysis, Sintering, and Redispersion, Van Nostrand Reinhold, New York, 1987; E. Ruckenstein, Mater. Sci. Res. 16 (1984) 199, in: G.C. Kuczynski, A.E. Miller, G.A. Sargent (Eds.), Sintering and Heterogeneous Catalysis, Plenum Press, New York, 1984. [14] J.F. van der Veen, in: H. Traub (Ed.), Phase Transitions in Surface Films, vol. 2, Plenum Press, New York, 1991. [15] W. Lipowsky, Phys. Rev. Lett. 49 (1982) 1575. [16] W. Lipowsky, W. Speth, Phys. Rev. B 28 (1983) 3983. [17] W. Lipowsky, U. Breuer, K.C. Prince, H.P. Bonzel, Phys. Rev. Lett. 62 (1989) 913.
ARTICLE IN PRESS 144
Y. Djikaev, E. Ruckenstein / Physica A 387 (2008) 134–144
[18] J.M. Kosterlitz, D.J. Thouless, J. Phys. C: Solid State Phys. 5 (1972) L124; J.M. Kosterlitz, D.J. Thouless, J. Phys. C: Solid State Phys. 6 (1972) L124. [19] A. Traynov, E. Tossati, Phys. Rev. B 38 (1989) 6961. [20] R. Ohnesorge, H. Lwen, H. Wagner, Phys. Rev. E 50 (1994) 4801. [21] C. Kittel, Introduction to Solid State Physics, Wiley, New York, 1968. [22] S. Chandrasekhar, Rev. Mod. Phys. 15 (1949) 1. [23] C.W. Gardiner, Handbook of Stochastic Methods, Springer, New York, Berlin, 1983. [24] N. Agmon, J. Chem. Phys. 81 (1984) 3644. [25] G. Narsimhan, E. Ruckenstein, J. Colloid Interface Sci. 128 (1989) 549. [26] E. Ruckenstein, B. Nowakowski, J. Colloid Interface Sci. 137 (1990) 583. [27] B. Nowakowski, E. Ruckenstein, J. Colloid Interface Sci. 139 (1990) 500. [28] S. Park, M.K. Sener, D. Lu, K. Schulten, J. Chem. Phys. 119 (2003) 1313. [29] M.-L. Tan, J.T. Fischer, A. Chandra, B.R. Brooks, T. Ichiye, Chem. Phys. Lett. 376 (2003) 646–652. [30] L.G. Dowell, A.P. Rinfret, Nature 188 (1960) 1144–1148. [31] X. Wei, P.B. Miranda, Y.R. Shen, Phys. Rev. Lett. 86 (2001) 1554–1557.
Physica A 387 (2008) 134–144 www.elsevier.com/locate/physa
A kinetic model for the premelting of a crystalline structure Yuri Djikaev, Eli Ruckenstein Department of Chemical and Biological Engineering, SUNY at Buffalo, Buffalo, NY 14260, USA Received 18 May 2007 Available online 19 August 2007
Abstract An analytical kinetic approach to examine the premelting phenomenon is suggested by using a first passage time analysis. Premelting is considered to occur when the time of formation of a Frenkel type defect in the surface monolayer becomes sufficiently small. The mean time of defect formation on the surface lattice, i.e., the mean time necessary for a selected (surface-located) molecule to leave its lattice site and form a Frenkel defect, is calculated by using a first passage time analysis. The model is illustrated by numerical calculations for a crystalline structure composed of molecules interacting via the Lennard-Jones (LJ) potential. The lattice vectors in the plane parallel to the free surface of the crystal were assumed to be equal (to the lattice parameter) and the angle between them was varied. The model predictions of the Tammann temperature (of premelting) are very sensitive to the parameters of the LJ potential. In all the cases considered, the temperature dependence of the mean first passage time has two clearly distinct regimes: at low temperatures the dependence is sharp and at high temperatures it is weak. r 2007 Elsevier B.V. All rights reserved. PACS: 05.10.Gg; 05.70.Fh; 05.70.Np Keywords: Premelting; Tammann temperature; First passage time analysis
1. Introduction Although some theories of melting consider the process as a breakdown of the crystal lattice which occurs uniformly and homogeneously throughout the solid at the melting point [1], there is evidences that melting is a heterogeneous process involving the nucleation of a liquid phase at some preferred sites of the solid (the free surface, grain boundaries, large dislocations and disclinations, etc.) and subsequent growth of the liquid phase. The free surface of a crystalline structure constitutes a natural ‘‘heterogeneous’’ site for the onset of the crystal-to-liquid phase transition (as it does for the liquid-to-crystal transition [2]). Tammann [3] was the first to point out that the surface premelts before melting occurs. This idea is consistent with the empirical criterion proposed by Lindemann [1] for melting in the bulk according to which melting might be expected to occur when the root mean amplitude of thermal vibrations of an atom exceeds a certain threshold value of more or about 10% of the distance to the nearest neighbor in the crystalline Corresponding author. Tel.: +1 716 645 2911x2214; fax: +1 716 645 3822.
E-mail addresses: [email protected] (Y. Djikaev), [email protected] (E. Ruckenstein). 0378-4371/$ - see front matter r 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.physa.2007.08.022
ARTICLE IN PRESS Y. Djikaev, E. Ruckenstein / Physica A 387 (2008) 134–144
135
structure. Tammann suggested [3] that the outermost layer of the crystal should acquire some disorder far below the bulk melting point. This behavior is a result of the higher freedom of motion for the surface-located molecules. The surface melting (often referred to as premelting) involves the formation of a thin layer which acquires some disorder at a temperature significantly below the melting one. The surface premelting can provide a plausible explanation to why a liquid can be cooled below the crystallization point (i.e., supercooled), whereas a solid cannot be heated above the melting point (i.e., superheated). Such a hysteresis suggests that there is no free-energetic barrier for the melting phase transition which occurs as spinodal decomposition (a barrier exists for homogeneous crystallization which thus always occurs via nucleation). As a possible cause, one can suggest that the surface of a solid, which premelts at temperatures far below the melting point, constitutes a large natural ‘‘heterogeneous’’ center ‘‘whereupon’’ the melting phase transition occurs in a barrierless way. Experimentally the premelting phenomenon was apparently first detected by Lyon and Somorjai [4] who studied the structures of the clean (1 1 1), (1 1 0), and (1 0 0) crystal faces of platinum as a function of temperature by means of low-energy electron diffraction and observed the formation of disordered surface structures at temperatures far below the melting temperature of T m ¼ 2043 K. Direct experiments on the surface initiated melting were also carried out by Frenken et al. [5] using Rutherford backscattering in conjunction with ion-shadowing and blocking. That experiment revealed a reversible order–disorder transition on the (1 1 0) surface of a lead crystal well below its melting point of T m ¼ 600:7 K. Since then, other techniques have been employed such as calorimetry, electron, neutron, and X-ray diffraction, microscopy, ellipsometry, and helium scattering. Although most experiments were carried out under equilibrium conditions, melting tended to be initiated at the surface even when the crystalline solid was heated very quickly so that equilibrium conditions were not established [6]. With the advent of superpowerful computers, molecular dynamics (MD) simulation methods have become an invaluable and major tool in studying the premelting phenomenon at a microscopic level [7–9]. Among recent applications of these methods one can mention the simulation of premelting in AgBr [10], premelting in Cr2O3 [11], the premelting of a clean Al(1 1 0) surface, [12], etc. Still, theoretical models of premelting have proven to be useful in elucidating various aspects of this phenomenon (some examples are given in Ref. [13]). In the thermodynamic approach premelting is regarded as a particular case of wetting, namely the wetting of a solid by its own melt [14]. Another phenomenological model of premelting was developed [15–17] in the framework of the Landau–Ginsburg theory wherein the free energy of the system is expanded in a power series with respect to an order parameter whereof the temperature dependence is then analyzed. A microscopic theory of melting of two-dimensional solids on the basis of the dislocations pairs model was developed in Ref. [18] and discussed in Ref. [13]. Another microscopic theory of surface melting [19] was based on a lattice model of a solid and explored the fact that by using a discrete reference lattice and a mean-field approach, the partition function of a system of particles interacting via a pairwise potential could be calculated. An efficient theoretical approach [20] to study premelting was developed on the basis of the density functional theory. In that approach the grand canonical free energy of an inhomogeneous system (for given chemical potential(s), volume, and temperature) is regarded as a functional of the local density, the equilibrium spatial profile of which minimizes the free energy functional. As an attempt to clarify some kinetic aspects of premelting, in this paper we present a microscopic model of the process on the basis of a first passage time analysis. The model assumes that the thermal motion of a molecule located in the surface layer of a finite crystalline system is determined only by its interactions with other surface-located (in the lattice sites) molecules. The main idea of the model is to determine the mean time of a defect formation on the lattice surface, i.e., the mean time necessary for a selected (surface-located) molecule to leave its lattice site and go wandering on the surface. The temperature dependence of this time exhibits two clearly distinct regimes, the high- and low-temperature ones, with a relatively non-sharp transition between them. The onset temperature of the high-temperature regime can be roughly associated with the Tammann temperature of premelting. It should be noted that, unlike the boiling or melting transitions whereof the temperatures can be unambiguously defined, there exists no strict definition of the Tammann temperature. The transition of the surface layer of a crystalline body from an orderer state into a
ARTICLE IN PRESS Y. Djikaev, E. Ruckenstein / Physica A 387 (2008) 134–144
136
premelted one is apparently a gradual process and one can associate the Tammann temperature to that which corresponds to one or another degree of ‘‘surface’’ disordering. The paper is structured as follows. In Section 2 we present a model for a semi-infinite crystalline structure with a surface layer. The basics of the mean first passage time analysis and its application to premelting are given in Section 3, while Section 4 contains numerical calculations based on our model. The results are discussed and the conclusions are summarized in Section 5. 2. Semi-infinite crystalline system of pairwise interacting molecules Consider a semi-infinite crystalline system such that in the Cartesian coordinates it will occupy the space 1oxo1; 1oyo1; 1ozp0 so that the plane z ¼ 0 constitutes a free surface of the crystalline body. This crystalline structure is defined by three lattice vectors a, b, and c with a and b in the plane parallel to the plane x–y of the Cartesian coordinate system. The angle between a and b will be denoted Y. The crystal lattice sites are occupied by units hereafter assumed to interact via the Lennard-Jones (LJ) pair potential and henceforth termed ‘‘molecules’’. In reality, they can be atoms, actual molecules, groups of atoms or molecules, colloidal particles, etc. Important is the pairwise and central character of their interactions. Consider a molecule at a lattice site with the coordinates x ¼ 0; y ¼ 0; z ¼ 0 (i.e., in the origin of the coordinate system) and denote it by s (this is one of the molecules belonging to the surface layer of the crystalline structure). Clearly, the molecule s is not frozen at its lattice site but performs constant thermal motion of (nonlinear) vibrational type. The x–y motion of the molecule s is completely determined by the x and y components of the forces acting thereupon. Assuming a continuum distribution of molecules below the surface layer, it is clear that F sxy , the xy component of the total force acting on the molecule s, is determined only by its interactions with the other surface-located molecules. Indeed, due to the continuum distribution of molecules in the interior of the crystal, for any elementary volume thereof there exists a counterpart located symmetrically (with respect to the axis z) such that the xy components of the forces exerted on the molecule s by these two elements are compensated. As a result, the force F sxy can be represented as F sxy ¼
N 1 X
f is ,
(1)
i¼1
where f is is a force exerted by the molecule associated with the lattice site i upon the selected molecule s and N is the total number of surface located molecules (including the selected one). Except for the molecule s, all other surface located molecules are assumed to be ‘‘frozen’’ at the corresponding lattice sites, so that for a given crystalline structure F sxy is a function of the x and y coordinates of molecule s only: F sxy F sxy ðxs ; ys Þ. It should be noted that, in reality, molecule s performsPmotion not only along the x- and y-axis, but also along the z-axis. Therefore, strictly speaking, the force N1 and its x and y i¼1 f is has also a z component P components depend on the coordinate z. However, it is clear that the z component of N1 i¼1 f is does not affect the x–y motion of molecule s. Besides, the temperature involved in premelting is much lower than the sublimation temperature. For this reason the motion in the z direction will be neglected. Let us assume that the interaction potential is pairwise and denote the potential between molecules i and s by fis ðx; yÞ.PThe overall potential of interactions of molecule s with all other surface located molecules is then js ðx; yÞ ¼ N1 i¼1 fis ðx; yÞ. Fig. 1 presents the function js ¼ js ðx; yÞ for the potential fis ðx; yÞ of the LJ type " 6 # Z 12 Z fis ðx; yÞ ¼ 4e , (2) ris ris where e and Z are the standard parameters of the LJ potential. Fig. 1a shows the contour plot of a threedimensional landscape determined by js ¼ js ðx; yÞ, whereas in Figs. 1b–d various sections of this landscape are plotted, corresponding to x ¼ 0 (similar to y ¼ 0), y ¼ x (similar to y ¼ x), and x ¼ 0:5 (similar to y ¼ 0:5), respectively (all the distances are given in units of the lattice parameter L, the potential is given in units of kB T, and the LJ parameters are e=kB T ¼ 0:5, Z=L ¼ 0:5). Clearly, the force F sxy and the potential
ARTICLE IN PRESS Y. Djikaev, E. Ruckenstein / Physica A 387 (2008) 134–144
137
Fig. 1. (a) Contour plot and various profiles of the potential field for a selected molecule (molecule s) associated with the lattice site ð0; 0Þ (for a lattice with a ¼ bð¼ LÞ; Z=L ¼ 0:5; Y ¼ p=2, and e=kB T ¼ 0:5). The darker regions correspond to a lower potential C and the lattice sites are indicated as points. (b)–(d) The profiles of the landscape (determined by Cðx; yÞ) along its various sections (indicated in the figure panels).
js ðx; yÞ are related through F sxy ¼ =js ðx; yÞ.
(3)
The potential field wherein the molecule s performs its thermal motion is a complicated array of potential wells with infinite barriers between them in some directions and finite ones in others. Two important features of this potential field should be noted. First, the site s of the lattice is surrounded by a series of local wells in each of which the molecule s can be located with equal probability. The selected molecule s, considered to be ‘‘at site s’’, would actually be located in one of these wells or would be in the process of transition from one of them into another, so that at any given moment the molecule s can actually be found (with various probabilities) in any place in the vicinity of the lattice site to which it is attributed. Second, there are eight pathways whereby the selected molecule s can leave the vicinity of its native lattice site, cross over a finiteheight barrier, and fall into another local well in the potential field. This would correspond to the formation of a surface version of a Frenkel type defect of the crystalline structure [21]. The probability (or the frequency) of such transitions depends on Z; e; L, and temperature T (hereafter expressed in K). Once these transitions
ARTICLE IN PRESS 138
Y. Djikaev, E. Ruckenstein / Physica A 387 (2008) 134–144
become very frequent, the surface layer of the lattice can be considered to have lost its ordered structure. Thus, for a given crystalline system (i.e., given Z; L and e) one can expect to detect premelting with increasing temperature T by examining the change of the mean first passage time, necessary for the molecule s to leave the vicinity of its native lattice site, with variations in the temperature T. In the following section, we outline the basics of a mean first passage time analysis and its applications to the problem of premelting. 3. A mean first passage time analysis Upon close examination of the landscape in Figs. 1a–d, one can conclude that the molecule s can leave the vicinity of its original site by passing over a barrier having the shape of a hyperbolic paraboloid (‘‘saddle’’) with the path of steepest descent parallel either to the x- or y-axis (Fig. 2a). Although the barrier is two-dimensional, it is very narrow over the saddle point hence one can assume that the molecule crosses the barrier by means of one-dimensional motion either along the x- or y-axis. One can therefore formulate a mean first passage time analysis for the chaotic motion of a molecule in a one-dimensional potential well. Let us consider the thermal chaotic motion of a molecule in a one-dimensional potential field CðxÞ (assumed to be expressed in kB T units) as shown in Fig. 2b. The field has the shape of a well with one of the boundaries (at x ¼ 0 in Fig. 2) impenetrable to the molecule and the other penetrable (at x ¼ l) so that the molecule can be regarded to have escaped from the well once it reaches this boundary. The mean first passage time of a molecule escaping from the well is calculated on the basis of a kinetic equation governing the chaotic motion of the molecule in this potential well. In general, the chaotic motion of a molecule subjected to an external force is governed by the Fokker–Planck equation for the single-molecule distribution function with respect to its coordinates and momenta, i.e., in the phase space [22–24]. Assuming that the relaxation time for the velocity distribution function is small compared to the characteristic time scale of the escape process [22], the Fokker–Planck equation reduces to the Smoluchowski equation [22–24] for the diffusion of a particle in an external force field which can be written in the form (see, e.g., Eq. (312) in Ref. [22]) qpðr; tjr0 Þ F ¼ = D =pðr; tjr0 Þ pðr; tjr0 Þ , (4) qt kB T where pðr; tjr0 Þ is the probability of observing a molecule in the volume element dr at r at time t given that initially it was at r0 , D is the diffusion coefficient, and F FðrÞ is the force exerted on the molecule by an
Fig. 2. (a) The vicinity of the point ð1; 0:5Þ for the potential field shown in Fig. 1a. In this vicinity the potential has a shape of a hyperbolic paraboloid with the path of steepest descent parallel to the x-axis. (b) One-dimensional potential well CðxÞ wherein a molecule is assumed to perform thermal chaotic motion before leaving the vicinity of its original lattice site. The well boundary at x ¼ 0 is assumed to be of infinite height, whereas at x ¼ l it is finite.
ARTICLE IN PRESS Y. Djikaev, E. Ruckenstein / Physica A 387 (2008) 134–144
139
external field. If the latter can be derived from some potential fðrÞ as F ¼ =fðrÞ, Eq. (4) can be rewritten as qpðr; tjr0 Þ ¼ =½DefðrÞ=kB T =ðefðrÞ=kB T pðr; tjr0 ÞÞ. (5) qt For the one-dimensional potential field CðxÞ ¼ fðxÞ=kB T and constant diffusion coefficient D, the Smoluchowski equation (5) reduces to qpðx; tjx0 Þ q CðxÞ q CðxÞ ¼D e ðe pðx; tjx0 ÞÞ . (6) qt qx qx The mean minimum time, necessary for a molecule in the well to cross over the external boundary of the well and leave it, depends on the initial position x0 of the molecule. It is convenient to use the backward Smoluchowski equation [23,25–28] which expresses the dependence of the transition probability pðx; tjx0 Þ on x0 : qpðx; tjx0 Þ q q ¼ DeCðx0 Þ eCðx0 Þ pðx; tjx0 Þ . (7) qt qx0 qx0 The probability that a molecule, initially at a distance x0 from the inner boundary of the well, will remain in the well after time t is given by the so-called survival probability, [25–28] Z l Qðtjx0 Þ ¼ dxpðx; tjx0 Þ ð0ox0 olÞ, (8) 0
where l is the width of the potential well, i.e., the distance between the inner and outer boundaries of the well (see Fig. 2b). The probability for the dissociation time to be between 0 and t is equal to 1 Qðtjx0 Þ, and the probability density for the dissociation time to be between 0 and t is given by qQ=qt. The first passage time is provided by [25–28] Z 1 Z 1 qQðtjx0 Þ dt ¼ tðx0 Þ ¼ t Qðtjx0 Þ dt, (9) qt 0 0 where the rightmost equality is obtained by integrating the integral on its LHS by parts. The equation for the first passage time is obtained by integrating the backward Smoluchowski equation (7) with respect to x and t over the entire range and using the boundary conditions Qð0jx0 Þ ¼ 1 and Qðtjx0 Þ ! 0 as t ! 1 for any x0 (the latter means that if the molecule is initially in the well, sooner or later it will leave the well, whereas the former trivially says that the molecule, initially in the well, cannot jump out of it instantaneously). Thus, one obtains Cðx0 Þ d Cðx0 Þ d De e tðx0 Þ ¼ 1. (10) dx0 dx0 One can solve Eq. (10) by assuming a reflecting boundary condition dt=dx0 ¼ 0 at x0 ¼ 0 and a complete absorption boundary condition tðx0 Þ ¼ 0 at x0 ¼ l. Consequently, tðx0 Þ is given by Z Z y 1 l dy eCðyÞ dx eCðxÞ . (11) tðx0 Þ ¼ D x0 R The average dissociation time, t¯ (often referred to as the mean first passage time), is obtained by averaging tðx0 Þ with the Boltzmann factor over all possible initial positions x0 : Z 1 l t¯ ¼ dx0 eCðx0 Þ tðx0 Þ, Z 0 Z Z¼
l
dx0 eCðx0 Þ .
(12)
0
Note that the mean first passage is inversely proportional to the diffusion coefficient D of a molecule in the potential well. Although this quantity may not be readily available for calculating the exact value of t¯ , even the
ARTICLE IN PRESS 140
Y. Djikaev, E. Ruckenstein / Physica A 387 (2008) 134–144
dimensionless quantity t¯ D=Z2 as a function of kB T=e can provide some information on the premelting behavior of the system because it relates the average time of formation of a surface defect of Frenkel type to the temperature. A drastic change in this behavior would indicate a qualitative change in the state of the surface layer of a crystalline structure with changing temperature. 4. Numerical evaluations As an illustration of the above model (described in Section 2), we carried out numerical calculations for crystalline lattices consisting of molecules interacting via the LJ potential, Eq. (2). The lattice vectors in the plane x–y (parallel to the surface) were assumed to be of equal length, a ¼ bð¼ LÞ, and lattices with various Z=L and Y (angle between the vectors a and b) were considered. At fixed Z=L and Y the main quantity of interest is the dimensionless quantity t¯ D=Z2 (representative of the mean first passage time necessary for a selected molecule s to leave the vicinity of its original lattice site by overcoming a barrier) as a function of the dimensionless temperature kB T=e. The results of numerical calculations are shown in Figs. 1, 3–5. In Fig. 1a (representing the case a ¼ bð¼ LÞ; Z=L ¼ 0:5; Y ¼ p=2, and kB T=e ¼ 2) the darker regions correspond to a lower potential C and the lattice sites, indicated as points, lie in the centers of the white regions. As clear from this figure, the selected molecule s can leave the vicinity of its original lattice site by crossing the barrier between one of the eight local potential wells and another of the wells not belonging to the site s. Although both potential wells are two-dimensional and the barrier between them has a ‘‘saddle’’-like shape, the bottleneck of the barrier is very narrow hence the mean first passage time for this transition can be estimated by considering the transition to occur via one-dimensional motion of the molecule along the path of the steepest descent from one well into the other. Eq. (12) can then be applied to evaluate the dimensionless ratio t¯ D=L2 . Figs. 1b–d present the profiles of the landscape (determined by Cðx; yÞ) along its various sections. They show that molecule s cannot leave the vicinity of its original lattice site (with coordinates x ¼ 0; y ¼ 0) along x ¼ 0, y ¼ 0, y ¼ x, and y ¼ x, but only along x ¼ 0:5; x ¼ 0:5; y ¼ 0:5, and y ¼ 0:5. Fig. 3a shows the dependence of the dimensionless mean first passage time t¯ D=L2 on the dimensionless temperature kB T=e for the case a ¼ bð¼ LÞ; Y ¼ p=2. Each curve corresponds to a fixed ratio Z=L (as specified in the caption of Fig. 3a). The examination of the slopes of these curves indicates that there are two distinct regimes for the variation of the mean first passage time with changes in temperature (e assumed fixed). At high temperatures (large kB T=e) the mean first passage time is hardly sensitive thereto, whereas at low temperatures (small kB T=e) the mean first passage time sharply changes with changing temperature. One can conjecture that these two regimes correspond to two different states of the surface layer of the crystalline structure, associating the low-temperature state with an ordered structure and the high-temperature state with a disordered structure. The transition between these two states is not very sharp but occurs rather smoothly with changing
Fig. 3. (a) Dimensionless mean first passage time t¯ D=L2 vs the dimensionless temperature kB T=e for a ¼ bð¼ LÞ; Y ¼ p=2. Each curve corresponds to a fixed ratio Z=L (0.53, 0.52, 0.51, 0.5, 0.49, 048, 0.47 from right to left). (b) The dimensionless Tammann temperature kB T m2 =e as a function of Z=L for Y ¼ p=2.
ARTICLE IN PRESS Y. Djikaev, E. Ruckenstein / Physica A 387 (2008) 134–144
141
Fig. 4. (a) Dimensionless mean first passage time t¯ D=L2 as a function of the dimensionless temperature kB T=e. Each curve corresponds to a fixed Y (0:99p=2; 0:98p=2; 0:97p=2; 0:96p=2; 0:95p=2; 0:94p=2; 0:93p=2 from right to left). (b) The dependence of the (dimensionless) Tammann temperature kB T m2 =e on Y for Z=L ¼ 0:5.
temperature. The onset of the high-temperature regime can probably be regarded as a rough estimate of the Tammann temperature at which the premelting of a crystalline system is said to occur. The uncertainty of a thus determined Tammann temperature is consistent with the somewhat arbitrary and non-rigorous character of the original definition of the Tammann temperature. In Fig. 3b the dimensionless Tammann temperature kB T m2 =e (determined as the temperature at which t¯ is by an order of magnitude higher than its value in a hightemperature regime) is plotted as a function of Z=L. One can note a decrease in kB T m2 =e with decreasing Z=L. The LJ potential in the soft sticky dipole (SSD) model of water (in which a water molecule is considered to have only a single, center of mass, interaction site with a tetrahedrally coordinated sticky potential that regulates the tetrahedral coordination of neighboring molecules) are e ¼ 15:3 kJ=mol and Z ¼ 3:01 A˚ [29]. For the cubic ice L ¼ 6:36 A˚ [30]. Thus, in this model one can estimate the Tammann temperature for the cubic ice as T m2 ’ 185 K which is in agreement with the experimentally observed temperature of about 200 K for the premelting of the (0 0 0 1) face of hexagonal ice [31]. To examine the sensitivity of the above model to the angle Y, we calculated the dimensionless mean first passage time t¯ D=L2 as a function of temperature for several values of Y (with a ¼ bð¼ LÞ and Z=L ¼ 0:5). The results are presented in Fig. 4. In Fig. 4a each curve corresponds to a Y specified in the figure caption. One can conclude that the dimensionless mean first passage time t¯ D=L2 is sensitive to Y, and so is the Tammann temperature. This can be explained by the drastic deformation that the landscape of the potential field Cðx; yÞ suffers with even a small change in Y. The smaller is Y at a given Z=L, the higher is the surface number density of molecules (lattice sites on a free surface of the crystal structure). Hence, the premelting of the surface layer occurs easier (the Tammann temperature decreases) with increasing surface density of lattice sites (i.e., with decreasing Y). This is clearly seen in Fig. 4b which presents the dependence of the (dimensionless) Tammann temperature kB T m2 =e on Y for Z=L ¼ 0:5. The above results (Figs. 1,3,4) describe the process of initial formation of a Frenkel defect in a defectless surface layer of the lattice. In order to somewhat clarify how the presence of a defect (i.e., a vacancy) in the structure of the surface layer affects the possible pathways of the evolution of its neighbors, we considered also a lattice with a pre-existing surface Frenkel defect (Fig. 5), where a molecule s0 has left the vicinity of its original lattice site ðx; yÞ ¼ ð0; 0Þ and now occupies the position ðx; yÞ ¼ ð0:65; 1:35Þ (all distances are expressed in units of L) marked as a point in a circle. In this case, the potential field for the molecule s1 associated with the lattice site ð0; 1Þ (one of the nearest neighbors of the vacancy at ð0; 0Þ) is shown in Figs. 5a (contour plot of the landscape with the darker regions representing lower values of C and the lattice sites indicated as points in the centers of white regions) and 5b–h. As expected, the potential field for molecule s1 has become very asymmetric. As shown in Figs. 5d–g, six pathways (whereby molecule s1 can leave the vicinity of its original lattice site to become a Frenkel type defect) are not affected by the existence of the defect, whereas two pathways become more probable because of decreases in the potential barriers.
ARTICLE IN PRESS 142
Y. Djikaev, E. Ruckenstein / Physica A 387 (2008) 134–144
Fig. 5. (a) Contour plot and various profiles of the potential field for a selected molecule (molecule s1 ) associated with the lattice site ð0; 1Þ in a lattice with a pre-existing surface Frenkel defect, where a molecule s0 has left the vicinity of its original lattice site ðx; yÞ ¼ ð0; 0Þ and now occupies the position ðx; yÞ ¼ ð0:65; 1:35Þ marked as a point (for a lattice with a ¼ bð¼ LÞ; Z=L ¼ 0:5; Y ¼ p=2, and e=kB T ¼ 0:5). The darker regions correspond to a lower potential C and the lattice sites are indicated as points. (b)–(h) The profiles of the landscape (determined by Cðx; yÞ) along its various sections as indicated in the figure panels.
ARTICLE IN PRESS Y. Djikaev, E. Ruckenstein / Physica A 387 (2008) 134–144
143
5. Concluding remarks It was demonstrated both experimentally and by simulations that the surface melting (often referred to as premelting) involves the formation of a thin disordered layer at a temperature significantly below the melting one. Several theoretical models of premelting were developed to elucidate various aspects thereof. All previous models have focused on obtaining the temperature of the premelting-like transition at the crystal free surface (i.e., the Tammann temperature) by using ‘‘static’’ formalisms (classical thermodynamic, Landau–Ginsburg, lattice, or density functional theories). In this paper, we have proposed an analytical kinetic approach to premelting on the basis of a first passage time analysis. The model assumes that the thermal motion of a molecule located in the surface layer of a finite crystalline system is determined only by its interactions with other surface-located (in the lattice sites) molecules. The main idea of the model is to determine the mean time of the defect formation on the surface lattice, i.e., the mean time necessary for a selected (surface-located) molecule to leave its lattice site and form a Frenkel defect. Thus, the model is aimed at predicting conditions under which the structure of the surface layer of a crystalline system may acquire some disorder due to a sharp decrease in the time of formation of Frenkel type defects. As an illustration, we carried out numerical calculations for a crystalline structure composed of molecules (base units) interacting via the LJ potential. The lattice vectors in the plane parallel to the free surface of the crystal were assumed to be equal (to the lattice parameter) and the angle between them was varied. In all cases considered, the temperature dependence of the mean first passage time has two clearly distinct regimes. At low temperatures it is strong, whereas at high temperatures it is weak. The onset of the high-temperature regime provides an estimate of the Tammann temperature. It should be noted that an important approximation has been made in evaluating the mean first passage time necessary for a surface molecule to leave the vicinity of its original lattice site to generate a Frenkel defect. In this vicinity there are eight local potential wells wherein the molecule is located most of the time. These potential wells are two-dimensional. The molecule leaves its original lattice site vicinity by crossing a barrier (between this vicinity and the rest of the lattice) having the shape of a hyperbolic paraboloid. To calculate the mean first passage time of this process, we assumed the motion of the molecule in the local well to be onedimensional and to occur along the path of the steepest descent. This approximation can be avoided by starting with a two-dimensional Smoluchowski equation (describing the chaotic motion of a molecule in a two-dimensional force field) and deriving the expression for the mean first passage time. References [1] F. Lindemann, Z. Phys. 11 (1910) 609. [2] Y.S. Djikaev, et al., J. Phys. Chem. A 106 (2002) 10247–10253; Y.S. Djikaev, et al., J. Chem. Phys. 118 (2003) 6572–6581. [3] Tammann, Z. Phys. Chem. Stoechiom. Verwandtschalft 68 (1910) 205. [4] H.B. Lyon, G.A. Somorjai, J. Chem. Phys. 46 (1967) 2539. [5] J.W.M. Frenken, P.M.J. Mare´e, J.F. van der Veen, Phys. Rev. B 34 (1986) 7506. [6] H. Hakkinen, U. Landmann, Phys. Rev. Lett. 71 (1993) 1023. [7] O. Tomagnini, F. Ercolessi, S. Iarlori, F.D. Di Tolla, E. Tosatti, Phys. Rev. Lett. 76 (1996) 1118. [8] F.D. Tolla, F. Ercolessi, S. Iarlori, Surf. Sci. 211/212 (1989) 55. [9] S. Iarlori, P. Carnevali, F. Ercolessi, E. Tosatti, Surf. Sci. 211/212 (1989) 75. [10] A.K. Ivanov-Schitza, G.N. Mazob, E.S. Povolotskaya, S.N. Savvin, Solid State Ionics 173 (2004) 103–105. [11] M.A.S.M. Barrera, J.F. Sanz, L.J. A´lvarez, J.A. Odriozola, Phys. Rev. B 58 (1998) 6057–6062. [12] R. Zivieri, G. Santoro, V. Bortolani, Phys. Rev. B 62 (2000) 9985–9988. [13] E. Ruckenstein, in: S.A. Stevenson, J.A. Dumesic, R.T.K. Baker, E. Ruckenstein (Eds.), Metal-support Interactions in Catalysis, Sintering, and Redispersion, Van Nostrand Reinhold, New York, 1987; E. Ruckenstein, Mater. Sci. Res. 16 (1984) 199, in: G.C. Kuczynski, A.E. Miller, G.A. Sargent (Eds.), Sintering and Heterogeneous Catalysis, Plenum Press, New York, 1984. [14] J.F. van der Veen, in: H. Traub (Ed.), Phase Transitions in Surface Films, vol. 2, Plenum Press, New York, 1991. [15] W. Lipowsky, Phys. Rev. Lett. 49 (1982) 1575. [16] W. Lipowsky, W. Speth, Phys. Rev. B 28 (1983) 3983. [17] W. Lipowsky, U. Breuer, K.C. Prince, H.P. Bonzel, Phys. Rev. Lett. 62 (1989) 913.
ARTICLE IN PRESS 144
Y. Djikaev, E. Ruckenstein / Physica A 387 (2008) 134–144
[18] J.M. Kosterlitz, D.J. Thouless, J. Phys. C: Solid State Phys. 5 (1972) L124; J.M. Kosterlitz, D.J. Thouless, J. Phys. C: Solid State Phys. 6 (1972) L124. [19] A. Traynov, E. Tossati, Phys. Rev. B 38 (1989) 6961. [20] R. Ohnesorge, H. Lwen, H. Wagner, Phys. Rev. E 50 (1994) 4801. [21] C. Kittel, Introduction to Solid State Physics, Wiley, New York, 1968. [22] S. Chandrasekhar, Rev. Mod. Phys. 15 (1949) 1. [23] C.W. Gardiner, Handbook of Stochastic Methods, Springer, New York, Berlin, 1983. [24] N. Agmon, J. Chem. Phys. 81 (1984) 3644. [25] G. Narsimhan, E. Ruckenstein, J. Colloid Interface Sci. 128 (1989) 549. [26] E. Ruckenstein, B. Nowakowski, J. Colloid Interface Sci. 137 (1990) 583. [27] B. Nowakowski, E. Ruckenstein, J. Colloid Interface Sci. 139 (1990) 500. [28] S. Park, M.K. Sener, D. Lu, K. Schulten, J. Chem. Phys. 119 (2003) 1313. [29] M.-L. Tan, J.T. Fischer, A. Chandra, B.R. Brooks, T. Ichiye, Chem. Phys. Lett. 376 (2003) 646–652. [30] L.G. Dowell, A.P. Rinfret, Nature 188 (1960) 1144–1148. [31] X. Wei, P.B. Miranda, Y.R. Shen, Phys. Rev. Lett. 86 (2001) 1554–1557.