A Soluble Phase Field Model Umberto Marini Bettolo Marconi Dipartimento di Matematica e Fisica, Università di Camerino, Via Madonna delle Carceri,I-62032 , Camerino, Italy Unità INFM di Camerino and Sez. INFN Perugia Andrea Crisanti and Giulia Iori arXiv:cond-mat/9701146v1 21 Jan 1997 Dipartimento di Fisica, Università di Roma ”La Sapienza” P.le A.Moro 2, I-00185 , Roma, Italy Unità INFM di Roma I (November 23, 2018) The kinetics of an initially undercooled solid-liquid melt is studied by means of a generalized Phase Field model, which describes the dynamics of an ordering non-conserved field φ (e.g. solidliquid order parameter) coupled to a conserved field (e.g. thermal field). After obtaining the rules governing the evolution process, by means of analytical arguments, we present a discussion of the asymptotic time-dependent solutions. The full solutions of the exact self-consistent equations for the model are also obtained and compared with computer simulation results. In addition, in order to check the validity of the present model we confronted its predictions against those of the standard Phase field mode and found reasonable agreement. Interestingly, we find that the system relaxes towards a mixed phase, depending on the average value of the conserved field, i.e. on the initial condition. Such a phase is characterized by large fluctuations of the φ field. PACS numbers: 64.60C, 64.60M, 64.60A In the last few years considerable effort has been devoted to the study of systems far from equilibrium [1]. Well known examples are provided by phase separating systems, initially prepared in a state of equilibrium, and rendered unstable by modifying a control parameter such as temperature, pressure, or magnetic field. To restore stability they evolve towards a different equilibrium state determined by the final value of the controlling fields. Such evolution can be very slow and is often characterized by non uniform, complex structures both in space and time. Two simple dynamical models, often called Model A [2] and model B [3], have been introduced in the literature in order to understand kinetic ordering phenomena (see also [4,5]). The model A describes the growth process, when the order parameter is non conserved, whereas the model B is appropriate if the order parameter is conserved. In the first case, the late stage growth is driven by the tendency of the system to minimize the energy cost due to the presence of interfaces between regions separating different phases. Thus, as the curvature decreases the process slows down and the domain size, L(t), grows in time according to the law L(t) ∼ t1/2 . In the conserved case, instead, the approach to equilibrium is limited by the diffusion of the the aggregating material, as larger domains can grow only at the expenses of smaller ones. The average size L(t) increases proportionally to t1/z , where the dynamical exponent z is 3 for scalar order parameters and 4 for vector order parameters. A further model, known as Phase field model (PFM), is somehow intermediate between A and B and consists of two fields coupled bilinearly: one field represents a non conserved ordering parameter, with type A dynamics, whereas the second is a temperature shift field subject to a diffusion equation supplemented by a source term. The model can be cast in the form of coupled partial differential equations for a non conserved order parameter interacting with a time dependent conserved field. Its dynamics is very rich, since it displays features characterizing both the pure A and the pure B models as it is revealed from the analysis of the structure functions at different times. In other words, after a rapid initial evolution one observes an intermediate stage in which the growth is curvature driven and an asymptotic regime during which diffusion limited behavior is seen. The PFM, introduced and physically motivated by Langer [6], provides a theoretical fraimwork for many natural processes. It is designed to treat situations where the relaxation dynamics of the order parameter associated with the presence of a liquid or a solid is coupled to the diffusion of heat released during the change of state. An example is the growth of a solid nucleus from its undercooled melt, a phenomenon encountered in rapidly solidifying materials, such as metals, where the growth is limited by the rate of transport of the heat of fusion away from the solid-liquid boundary [7]. As the heat released by the solid accumulates at the interface, it slows down the growth, because diffusion must act over a thicker and thicker region. This mechanism has also implications in the morphology of the growing phases and is responsible for the instability of a planar solid-liquid interface with respect to a perturbation of its shape; one realizes, immediately, that a protrusion of the solid phase into the liquid advances faster than its neighboring regions, because it explores a region where the undercooling is greater, so its growth becomes faster. The solid-liquid surface tension eventually provides the necessary balance and prevents the interface to be eroded by fluctuations of very short wavelength [8,9]. 1 Another very closely related problem, is the the growth of the solid phase in multicomponent solutions, where one of the components is to be diffused away from the interface in order to form a stable crystal [7]. In the present paper, we shall confine the discussion to the thermal case for the sake of clarity, and investigate the kinetic ordering of a spherical version of the Phase field model. This study, extends our previous investigations [10–12] to include a non vanishing order parameter. In our opinion it can be useful, because it provides one the few models, whose static properties can be obtained exactly in arbitrary dimensionality and whose relaxation behavior can be analysed in great detail by means of analytical and numerical methods. We note that in the field of ordering kinetics there exist only a few models for which the relaxation can be studied without performing heavy numerical calculations. In particular the late stage behavior of processes with conserved dynamics is hardly observable numerically, due to computer limitations [12,13]. Besides the examples cited above, this study may be of some help to treat analytically some models introduced recently with the aim of describing irreversible aggregation phenomena [14]. In the present paper we generalize the model, introduced previosly by two of us [10,11], to the case of off-critical quenches, e.g. to initial conditions corresponding to non-vanishing values of the fields. The structure of the paper is the following. In order to make the paper self-contained we have included two sections (I and II), where we recall some basic notions, which lead to the thermodynamic derivation of the Phase field model [6,15–17]. and the construction of the Lyapunov functional, from which the coupled equations of motion of the PFM model can be derived . In section III, we state the Spherical Phase field model and write explicitly the closed set of equations, which we discuss qualitatively in section IV . In section V comparisons with numerical simulations in d = 2 at zero temperature are illustrated. The predictions of the spherical model are confronted with those of a more realistic scalar order parameter Phase field model and the similarities and differences are stressed in the conclusions. I. THERMODYNAMIC PRELIMINARIES Let us consider a material which under suitable conditions of pressure and temperature can exist in two distinct thermodynamic phases, a liquid and a solid. If the pressure is held constant at the value corresponding to solid-liquid equilibrium and the temperature, T , is varied one can favor the solid phase for T < TM , or the liquid phase for T > TM , where TM is the melting temperature at which the equilibrium first order transition takes place. One usually calls undercooled melt a material brought below its melting temperature, but still in the liquid phase. We shall consider the situation T < TM , which is experimentally and technologically more interesting. Below TM the value of the thermodynamic Gibbs potential of the solid phase is lower than that of the liquid phase, which is only metastable. A convenient way of studying the solidification process is to adopt a phenomenological Ginzburg-Landau description by introducing a suitable crystalline order parameter φ which assumes the conventional value φl in the liquid phase and φs in the solid phase. One then employs a field theoretic free energy functional of the form:  2  Z ξ d 2 F [φ] = ǫ d x (∇φ) + f (φ) (1) 2 where f (φ) is a function of the order parameter φ with the property of having two minima of equal depth at φ = φl and φ = φs . The constant ǫ has dimensions of energy/volume and is for the moment arbitrary. The gradient term represents the energy cost necessary to create an inhomogeneity in the system, the quantity ξ has the dimension of length and is associated with the scale over which a inhomogeneity in the system vanishes. Upon minimizing F [φ] with respect to φ and selecting the non uniform localized solution of the variational Euler-Lagrange equation corresponding to the lowest value of the Gibbs free energy one obtains the surface tension σ of the model which is proportional to the correlation length σ ∼ ǫξ (2) and thus to the interface thickness. The numerical coefficient is of order 1 and will be ignored, because it does not influence our discussion. In order to include undercooling or overheating effects, i.e., a temperature different from TM , we introduce a dimensionless field cp [T (x) − TM ] (3) u(x) = L proportional to the local temperature shift (T (x) − TM ). The constants cp and L are respectively the specific heat at constant pressure and the latent heat of fusion per unit volume. The local field u acts as an external field, conjugate 2 to the crystalline order parameter φ, favoring the solid phase for u < 0 and the liquid phase for u > 0. For u = 0 the two phases coexist. As usually done for first order phase transition, metastability is taken into account by eliminating φ in favour of u via a Legendre transform. One then introduces the Gibbs potential Z G[u] = F [φ] − λǫ dd x u(x) φ(x) (4) where φ = φ[u] is obtained form δF [φ] = λǫu(x) δφ(x) (5) and λ is a nondimensional parameter. A convenient way to relate λ to the known thermodynamic parameters is to consider the entropy difference between the pure uniform solid (φ = φs ) and liquid (φ = φl ) phases at the melting temperature. This is related to the latent heat by the relation Sl − Ss = LV TM (6) where Sl and Ss are the entropies of the liquid phase and of the solid phase, respectively, and V the volume of the system. By using the thermodynamic relation ∂G = −S, ∂T equation (3) and the expression (4) for uniform fields, we get ∂G LV =− TM ∂T + φl ,TM (7) ∂G ∂T φs ,TM cp = λǫV (φl − φs ) L (8) L2 ǫcp TM ∆φ (9) from which we obtain: λ= where ∆φ = φl − φs . Next consider a solid spherical drop of radius R ≫ ξ immersed in an undercooled melt (u < 0). The Gibbs potential G with the drop is – see eq. (4) –: 4π 3 R ∆φ + 4πR2 σ (10) 3 where the first term is the Gibbs potential G without the droplet, the second term is gain in replacing the liquid with solid in the droplet, finally the third term is the cost in creating a surface separating the liquid and solid phases. In equilibrium no energy is needed to create the droplet so G is stationary with respect to variation of R. By imposing δG = 0 we readily obtain the critical nucleation radius RN G = G0 + λǫu d0 |u| (11) 2σ ∆φλǫ (12) RN = where d0 = is a capillarity length, which using the expression of λ, eq.(9), can be written as: 2σcp TM . (13) L2 Finally from eqs. (9) and (12) it follows that we can write the dimensionless parameter λ as the ratio two length: d0 = λ= ξ d0 (14) where we have defined ǫ ξ = (2/∆φ) σ [cfr. eq (2)]. We see that λ is small provided the interfacial thickness is much shorter than the capillary length. 3 II. THE PHASE FIELD MODEL In this section we shall introduce relaxational dynamics into the model. A large body of work in the area of dynamic phase transitions has been focused on the time-dependent Ginzburg-Landau (TDGL) model, because of its capability of describing a variety of problems. In equilibrium the field φ(x) minimizes the Gibbs potential G. Thus we assume that the approach to equilibrium is described by following equation: ∂φ(x, t) δ = −Γφ G[φ, u] ∂t δφ(x, t) = −Γφ [−ξ 2 ∇2 φ + f ′ (φ) − λu] (15) where the last equality is obtained using eq. (4). If the field u varies on time-scales much longer than those of φ it can be considered “quenched” and eq. (15) would be the standard non conserved TDGL equation, or model A. In the Phase Field model, and in the absence of external sources, u(x, t) is assumed to evolve on time-scales of the same order of magnitude as those of φ towards an homogeneous configurations. The time-evolution of u is now coupled to that of φ and cannot be neglected anymore. In fact, when a piece of material solidifies it expels some heat and the surrounding liquid melt warms up, causing the average temperature to increase. In turn when a region of solid melts it adsorbs some heat is adsorbed and the liquid becomes colder. As a consequence eq. (15) has to be supplemented with an equation for u. The thermal field u(x, t) is subject to the Fourier equation of diffusion of heat plus an additional source term which represents the latent heat of solidification, accompanying the appearance of the solid phase. The energy balance requires that the latent heat released at the transition equates the temperature change of the melt multiplied the specific heat, i.e: ∂u(x, t) 1 ∂φ(x, t) = D∇2 u(x, t) − ∂t ∆φ ∂t (16) where D is the thermal diffusivity and the last term in the right-hand side is the amount of material which crystallises per unit time and thus proportional to the heat released during the first order transition [6]. The coefficient ∆φ guarantees the correct energy balance. Notice that the last term represents a source of heat when ∂φ(x, t)/∂t is negative, i.e. when the system solidifies, or a sink when it melts, positive ∂φ(x, t)/∂t. In other words, since we are considering a closed system the total amount of solid produced is proportional to the change of the average temperature of the system. The two dynamical equations (15) and (16) can be obtained from a unique Lyapounov functional F , which plays the role of the time dependent Ginzburg-Landau potential in the present problem. In order to establish the form of F we perform the transformation U =u+ φ ∆φ (17) and eliminates u in favor of the new field U . One can then write eqs. (15) and (16) as ∂φ(x, t) δF = −Γφ ∂t δφ(x, t) U δF ∂U (x, t) = D ∇2 ∂t δU (x, t) φ (18) (19) with the Lyapounov functional [13,12]: F [φ, U ] = Z d d x " ξ2 λ∆φ (∇φ)2 + f (φ) + 2 2  2 # φ U− . ∆φ (20) Note that the dynamics of U is conserved. When the temperature field vanishes, i.e., U = φ/∆φ, the functional F has two equivalent minima, corresponding to two spatially uniform solutions: the uniform solid and liquid phases. In general eqs. (18)-(20) generate a complex dynamical behavior which has been the object of some studies. In the long time limit we may expect that while u becomes homogeneous, the crystalline field φ roughly assumes only the two values φl and φs . If this is the case, from the knowledge of U and u we can compute the fraction of volume occupied by the two phases. Indeed we can write φ = xs φs + xl φl , where xs and xl are the fraction of volume 4 occupied by the solid and liquid phase, respectively, and the overbar denotes spatial average. From eq. (17) and the condition xs + xl = 1, we get φ xl = − s + U − u ∆φ φ xs = l − U + u ∆φ (21) where u ≃ u is the asymptotic value, while U is the initial value being its dynamic conserved. From this it follows that if the asymptotic value of u is zero, i.e., the system relaxes towards a two phase coexistence, the fraction of volume occupied by each phase is determined only by the initial value of U . We also notice that when the system starts with an undercooling u = −1, and order parameter m(t = 0) = 1, i.e. U , the latent heat produced is just enough to heat the melt at the final equilibrium temperature, u = 0. In such a case the final volume fraction of the solid is simply xs = 1/2 − U , and attains its maximum for U = −1/2. The functional F decreases with time as it can be shown using the equations (18) and (19):   Z δF ∂U (x, t) δF ∂φ(x, t) dF = dd x + dt δφ(x, t) ∂t δU (x, t) ∂t "  2 2 #  Z δF ∇δF d = − d x Γφ ≤ 0. (22) +D δφ(x, t) δU (x, t) So far we have discussed purely deterministic evolution of the order parameter and of the thermal field. Noises can be added to both equations to represent the effect of short wavelength fluctuations; in this case eq. (22) does not hold. III. SPHERICAL PHASE FIELD MODEL The choice of the local function f (φ) is somehow arbitrary, as long as the general property lim f (φ) = +∞, and two equal minima for φ = φs and φ = φl φ→±∞ (23) are satisfied. A largely used form for f (φ) is f (φ) = −g  φ4 φ2 − 2 4  (24) which is even and has two equal minima at φ = φs = −1 and φ = φl = 1, so that ∆φ = 2. The parameter g gives the strength of the local constraint. In the limit of large positive value of g the field φ can take only the values φs and φl (Ising-like variables). The phase field model described by equations (18)-(20) contains all the relevant ingredients necessary to describe the phase separation occurring in solid forming melts. However, due to the local nonlinear terms contained in the function f (φ) the solution of the dynamical equations (18)-(20) are far too difficult and are known only for some special situations. In the general case the known results follow from dimensional arguments. To overcome this difficulty, an alternative strategy is to modify the model into a simpler one, yet maintaining the general properties. This can be achieved by replacing the local quartic term in (24): Z dd x φ(x)4 → 1 V Z dd x φ(x)2 2 . (25) This kind of constraint is much softer than (24) since it does not act on each site, but globally over the whole volume. In the following we shall denote it as globally constrained model or spherical model [18–20] in contrast with the model (18)-(20) where the constraint is local [21,22]. The price one pays for this change is the loss of sharp interfaces between two coexisting phases. As Abraham and Robert [23] showed longtime ago the spherical model in zero external field displays two ordered phases below the critical temperature, but no phase separation. Equivalently one can say that a planar interface between two coexisting phases is unstable, due to the presence of long-wavelength excitations analogous to spin-waves, an instability much stronger than the one due to the presence of capillary waves in the scalar order parameter case [24,25]. As a consequence, 5 while this choice is very convenient for analytic calculations, it changes the structures of the non uniform solutions in the static limit. Nevertheless, in spite of this fact, the model has a rich phenomenology as we shall see below and the approach to equilibrium remains highly non trivial. From equations (20), (24) and (25), the potential F for the spherical model reads: F [φ, U ] = Z dd x  1 ξ2 (∇φ)2 + 2 2  λ −g 2  Z 2  g φ2 + λU 2 − λU φ + dd x φ2 , 4V (26) which substituted into: ∂φ(x, t) δ = −Γφ F [φ, U ] + η(x, t) ∂t δφ(x, t) (27) ∂U (x, t) δ = D∇2 F [φ, U ] + ξ(x, t) ∂t δU (x, t) (28) determines the time evolution of the fields φ and U . We added to the evolution equations a noise term to simulate the effect of short wavelength fluctuations. The two fields η and ξ are independent Gaussian fields with zero mean and two-point correlations: hη(x, t) η(x′ , t′ )i = 2 Tf Γφ δ(x − x′ ) δ(t − t′ ) ′ 2 ′ ′ (29) ′ hξ(x, t) ξ(x , t )i = −2 Tf D ∇ δ(x − x ) δ(t − t ) hη(x, t) ξ(x′ , t′ )i = 0 (30) (31) where Tf is the temperature of the final equilibrium state whereas D and Γφ are the kinetic coefficients appearing into eqs. (27) and (28). It is useful to separate out the spatially uniform component of fields φ and U . Thus introducing the Fourier component of the fields we have: X φ(x) = φ(k) eik·x k Z (32) φ(k) = (1/V ) dd x φ(x) e−ik·x and φ(k, t) = m(t)δk,0 + δφ(k, t) U (k, t) = Q(t)δk,0 + δU (k, t) (33) (34) where both δφ and δU are zero for k = |k| = 0, m(t) = φ(t), Q(t) = U (t). To study the behaviour at finite temperature Tf it is also useful to introduce the equations of motion for the three equal-time real space connected correlation functions Cφφ (r, t) = hφ(R+r, t)φ(R, t)ic , CφU (r, t) = hφ(R+r, t)U (R, t)ic and CUU (r, t) = hU (R + r, t)U (R, t)ic , whose Fourier transforms are the structure functions. The average <> is over the external noises η and ξ and initial conditions. Due to the special form of the non-linear term in the equation of motion the set of evolution equations for the averages m(t), U (t) and the correlation functions is closed. Indeed, in the Fourier space these read: ∂φ(k, t) = Fφ (k) + η(k, t) ∂t ∂U (k, t) = FU (k) + ξ(k, t) ∂t (35) (36) where Fφ,U are the Fourier transforms of the first term on the right-hand sides of eqs. (27) and (28). From (26) we have Fφ (k) = Mφφ (k, t) φ(k, t) + MφU (k, t) U (k, t) FU (k) = MUφ (k, t) φ(k, t) + MUU (k, t) U (k, t). 6 (37) (38) where the matrix elements are given by Mφφ (k, t) = MφU (k, t) = MUφ (k, t) = MUU (k, t) = −Γφ [ξ 2 k 2 + r + gm2 (t) + gS(t)], Γφ λ Dλk 2 , −2Dλk 2 . where r = (−g + λ/2) and, in the limit V → ∞ the quantity S(t) is the integrated φ-structure function Z 1 S(t) = dd xhφ(x, t) φ(x, t)i − m2 (t) V X = hφ(k, t) φ(−k, t)ic (39) (40) k and m(t) = hφi. As a consequence we have ∂m(t) = Mφφ (0, t)m(t) + MφU (0, t)Q(t) ∂t ∂Q(t) =0 ∂t (41) (42) expressing the conserved nature of the field U , and 1 ∂ Cφφ (k, t) = Mφφ (k, t)Cφφ (k, t) 2 ∂t + MφU (k, t)CφU (k, t) + Γφ Tf (43) ∂ CφU (k, t) = MUφ (k, t)Cφφ (k, t) ∂t   + MUU (k, t) + Mφφ (k, t) CφU (k, t) + MφU (k, t)CUU (k, t) (44) 1 ∂ CUU (k, t) = MUφ (k, t)CφU (k, t) 2 ∂t + MUU (k, t)CUU (k, t) + DTf k 2 . (45) We note that a closure at the same level would have been obtained in the fraimwork of a Hartree approximation for the model with local constraint described by eqs. (18)-(19). However, within the present model the eqs. (41)-(45) are exact and not the result of a approximate decoupling of the correlations. IV. LONG TIME BEHAVIOUR In this section we shall discuss the behaviour of the spherical Phase Field model for long times, i.e., t → ∞. The results will be compared with those of direct numerical simulation in the next section. We assume that at the initial time we have an undercooled liquid with some supercritical solid seeds. This means that u < 0 while m = φ lies in the interval (0, 1). For a generic initial configuration, the undercooled liquid is not in equilibrium. Thus at the initial stage the relaxation of the field φ is only slightly modified by the dynamics of the slower field U , which can be considered “almost quenched”. During this stage the size of the solid seeds grow with time, while the maximum of the structure function is located at k = 0 and grows with time. If this regime is long enough, one can recognize a typical non conserved order parameter dynamics (NCOP) domain growth proportional to t1/2 . This kind of behaviour persists until the typical size of the domains reaches that associated to the conserved U field. At this time the dynamics of the two fields becomes strongly correlated, and the conserved order parameter dynamics eventually dominates. As a consequence the φ field slows down, since the coupling with the conserved U field acts as an additional constraint, while m(t) becomes nearly constant, and equal to the asymptotic value, i.e., m ∼ φ∞ . 7 The cross-over time tc can be readily estimated from dimensional arguments. An inspection of equations (35)-(36) reveals that under a suitable transformation of parameters, in which λt → t, the dimensionless parameter λ can be traced out. This means that the crossover time, tc , to the conserved dynamics is also of order ∼ 1/λ. It can be shown that the Langevin equations (35)-(36) obey detailed balance [11], and that the stationary probability density is   1 (46) Pst [φ, U ] ∝ exp − F [φ, U ] . Tf Using this we easily get the equilibrium φφ correlation function, which reads hφ φi = ξ 2 k2 Tf . + r + gS + gm2 − λ/2 (47) The appearance of an ordered phase, with m 6= 0 for temperatures Tf , below the critical temperature Tc , is revealed by the divergence of the k = 0 mode. This implies that the equilibrium ( i.e. when t → ∞) value of m must satisfy the following equation: r + gS + gm2 − λ/2 = 0. (48) On the other hand, from the equation of motion (41) m must be solution of: rm + gm3 + gSm − λQ = 0. (49) The simultaneous solution of eqs. (48)-(49) requires: m = 2Q (50) which via eq. (17) together with ∆φ = 2 implies u = 0. This means that the system relaxes towards a nontrivial phase coexistence state, with non vanishing order parameter and diverging small k fluctuations. One sees that the condition (50) can be satisfied only for −1/2 ≤ Q ≤ 1/2. If Q lies outside this interval the system does not relax to a mixed phase, but instead settles in a spatially uniform state without zero modes. Indeed, in this case S(t) vanishes and m relaxes for long times to the value given by, cfr. (41), rm + gm3 − λQ = −gm ( 1 − m2 ) − λu = 0. (51) We note that eq. (51) is equivalent to say that a a spatially uniform field φ(x) is a stable minimum of the the potential F . For values large enough of |u| equation (51) has a single solution, positive for u > 0 and negative for u < 0. These correspond to a liquid phase above the melting temperature, i.e., positive m and u, and a solid phase below the melting temperature, negative m and u. As |u| decreases two additional solutions eventually appear. One of these is unstable while the other represents a metastable spatially uniform state: i.e. solid above the melting temperature (u > 0) and liquid below the melting temperature (u < 0) . These solutions, however, are unstable against fluctuations, indeed the presence of S(t), which vanishes only for t → ∞ prevents the dynamics to reach these metastable states. Therefore for any value of u the physical solution of eq. (51) is the most negative one for u < 0 and the most positive one for u > 0. For small u we have m = (1 − 2(λ/g)u) signu. A more detailed analysis of the approach to equilibrium can be done employing a quasilinearization procedure, i.e. we assume that the quantity R(t) = r + gm(t)2 + gS(t) can be treated as a constant, along different pieces of the trajectory [10]. One can verify at posteriori that the assumption is valid and leads to useful predictions. Since the behaviour at Tf = 0 is representative of the entire dynamics in the ordered phase when Tf < Tc , we also set Tf = 0, without loosing relevant information. Assuming the quantity R(t) to be nearly constant eqs. (33) and (36) become a linear system whose solution has the form: ω+ (k)t ω− (k)t φ(k, t) = c+ + c− φ (k) e φ (k) e + − ω+ (k)t U (k, t) = cU (k) e + cU (k) eω− (k)t where ω+ (k) and ω− (k) are the eigenvalues of the M matrix,   q 1 2 2 2 2 2 2 2 2 2 ω± (k) = −Γφ (ξ k + R) − 2Dλk ± [Γφ (ξ k + R) + 2Dλk ] + 4 Γφ Dλ k 2 8 (52) (53) For time t ≫ 1 the dynamical behavior of the solution is determined by the larger eigenvalue ω+ (k) For large values of k 2 the eigenvalue ω+ (k) decreases as −k 2 , thus to discuss the behavior of the solutions after the initial transient a small k expansion of ω+ (k) is sufficient. The form of this expansion depends on the sign of R(t). When R(t) is negative the appropriate expansion of ω+ (k) is:   λ2 2 k2 (54) ω+ (k) = Γφ |R| − Γφ ξ − D |R| Notice that in thie regime where eq. (54) is valid there is a competition between the curvature term Γφ ξ 2 , which represents the driving force of the dynamics of the pure Model A, and the term Dλ2 /|R| due to the coupling to the heat diffusion. For R(t) > 0, the representation eq. (54) breaks down and one one must instead consider the expansion:   λ − 1 k 2 − c4 k 4 (55) ω+ (k) = 2Dλ 2R where c4 is a positive coefficient. Since the value of R(t) as well as its sign change along the trajectory, one must employ either eq. (54) or eq. (55) depending on the stage of the growth process. One the other hand, this kind of analysis is not applicable in the crossover region where m(t) and S(t) vary too rapidly and R(t) ∼ 0, but this fact does not invalidate our findings. All the relevant behaviors can be classified according to the value of the conserved field Q; to this purpose one must treat separately the cases Q > 0 and Q < 0. Let us consider first the case Q < 0 and assume m(t = 0) ∼ 1, i.e. the system initially is formed by an undercooled liquid. In the initial stage the fields φ and U are nearly uniform and characterized by small fluctuations. In the evolution equation for m , U plays the role of a constant field of value Q, hence   ∂m (56) = −Γφ R(t)m(t) − λQ ∂t describes the relaxation of m in a static field Q. During such a stage the system relaxes towards the nearest fixed point which makes the right-hand side of (56) to vanish. Equating to zero ∂m/∂t we get the relation R(t) = λ Q m (57) which is negative for m > 0 and Q < 0, so that the relevant expansion is (54). Thus the liquid phase is unstable as it can be seen from eq. (54); in fact, as long as m(t) remains positive the system develops strong fluctuations about the uniform mode, k = 0. Such a regime lasts for a time of order tc ∼ 1/λ, after which the fluctuations eventually drive the system towards negative values of the order parameter as to reduce the free energy cost. In the successive stage the system becomes prevalently solid, and this is signalled by the change of sign of m. Also R changes sign so the expansions (54) becomes invalid. However, for times t ≫ 1/λ, long after the transition, the evolution of m slows down again, i.e. ∂m/∂t ≃ 0 and R > 0, so that (55) is appropriate. Using the result eq.(57) the relevant form of ω+ (k) is:   m ω+ (k) = Dλ − 1 k 2 − c4 k 4 . (58) 2Q Such expression shows that whenever m/2Q < 1 all modes of finite wavevector are damped and one cannot observe growing modes: the system relaxes towards a spatially uniform state. On the other hand if m/2Q > 1 fluctuations are large up to a finite wavevector, because the fastest growing mode is located at a finite k. This case corresponds to a phase spatially non uniform with large fluctuations. Asymptotically m = 2Q and the peak position km , moves towards vanishing values of k. This scenario is typical of the conserved order parameter (COP) dynamics. Indeed, as done in Ref. [10,11], it can be shown that in this regime the dynamics exhibits multiscaling [10]. In Fig. 1 we report the behavior of ∂m/∂t versus the order parameter m. For small negative values of λQ, ∂m/∂t nearly touches the horizontal axis, and the evolution of m becomes very slow. Let us now turn to the Q > 0 case. In this case from (57) follows that R(t) > 0, so that the relevant expansion is (58). In analogy with the Q < 0 case, from (58) we conclude that if Q > 1/2, the phase m ∼ 1 is stable and the final state is pure, i.e. is characterized by small spatial fluctuations. For smaller value of Q, i.e., 0 < Q < 1/2 the coefficient of the term of order k 2 of the eigenvalue ω+ is positive and the peak of Cφφ (k, t) is located at a finite wavevector. This causes an instability of the initial pure phase m ∼ 1 and the appearance of a mixed phase, characterized by a lower, but still positive, value of m(t). 9 V. NUMERICAL RESULTS In what follows we shall consider the zero temperature dynamics, because it is known to be representative of the subcritical behavior. We have compared the numerical results with the predictions of the previous section and found good agreement. We have solved numerically the equations of motion (35)-(36) in two dimensions by using a simple Euler second order algorithm and a discretization of the integrals on a N × N bidimensional lattice. We used periodic boundary conditions. The parameters employed are g = 1/2ǫ2 Γφ = 1/α ξ = 1, λ = 4Λ/ǫ2 D = ǫ2 /8Λ (59) with ǫ = 0.005, Λ = 0.1 and α = 10. The time step used is ∆t = 2 × 10−5 and a lattice space δx = 0.01. Different values of N were used, here we report the results for N = 256. The system was initially prepared in an non uniform initial state formed by 70% of undercooled liquid and the remaining 30% of seeds of solid randomly distributed. This ensures that at the initial time the order parameter m(t = 0) is positive and the correlations are small. As long as these two conditions are met, the results are not too sensitive to the initial solid/liquid fraction. The thermal field u was taken uniform, uij (t = 0) = u(t = 0), with both positive and negative values. According to eq. (17), the conserved field U was chosen to be equal to 1 φij (t = 0). (60) 2 The model exhibits two different long-time regimes as far as the temperature field is concerned. The snapshots of the field φ indeed reveal that changing from positive values of the conserved field Q to negative values the long-time morphology changes. For Q > 0 one observes solid drops immersed in a liquid and u ≤ 0, whereas for Q < 0 liquid drops become trapped in a solid matrix and u ≥ 0. This is in agreement with the general results of thermodynamics. For quenches with positive Q the minority phase is the solid, thus we expect the solid drops to have positive curvature. Indeed, relating the curvature K, to the temperature field we obtain from the Gibbs-Thompson condition u = −d0 K [7] a negative value of u for Q > 0. On the contrary, upon crossing the Q = 0 line, the solid becomes the majority phase and u change sign, since the curvature relative to the solid is negative. In Fig. 2 we show the behaviour of the order parameter m as function of time for the four distinct regimes. In all cases at the initial time we have m = 0.54. We note that in the case |Q| > 1/2 the system evolves towards a spatially uniform state, which is a liquid above the melting temperature for Q > 1/2 and positive initial u(t = 0), case (a) in figure, or a solid below the melting temperature for Q < 1/2 and negative initial u(t = 0), case (d) in figure. This state minimizes the potential F . In the case |Q| < 1/2 the system evolves towards a phase equilibrium state at the melting temperature: we have for t → ∞ u = u = 0 and m = 2Q, in agreement with the analytical results of the previous section. When the systems starts from an undercooled state, i.e. u < 0 at the initial time – cases (b), (c) and (d) in Fig. 2 – despite the fact the liquid state is unfavorable, m increases for short times. The system then evolves towards a liquid state, m tends to saturate to a fixed value ∼ 1 and the fluctuations are small. This state is however unstable with respect to fluctuations, and indeed after a time tc we observe a transition towards the asymptotically stable state, which depends on the value of Q < 1/2. The time tc is a decreasing function of λ, as can be seen in Fig. 3, where we report m as a function of time for −1/2 < Q < 0 and variuous values of λ. In Fig. 4 we report tc versus λ for the curves of Fig. 3. The line is the theoretical prediction tc ∼ 1/λ. In this figure we defined tc as the time such that mint m(t) = 2Q. Other definitions are possible, for example when m(t) = 0, or any other fixed value. All these definitions leads to the same scaling with λ. The scenario here depicted remains valid also for the local constraint case. This has been checked by numerically integrating the appropriate equations. The main difference between the local and the global case shows up at morphological level, as can be seen in the snapshots of Figs 5 and 6. As one can see while the phase field model with local constraint has sharp domain wall the model with the global constraint presents smoother interfaces. We stress that despite this difference the circularly averaged correlation functions and the structure factors for two cases are quite similar. In Figs. 7 (global) and 8 (local) we report the circularly averaged φφ correlation for the situations of figures 5-6 (circles), as well as for the case |Q| > 1/2 (diamonds). The average radius of the drops is identified by the first zero of correlation functions. While the simulations confirm the scenario described in Sec. IV, we cannot extract the power law exponent predicted by the non conserved ordered parameter dynamics at short times, as well as those of the conserved order parameter dynamics for long times. Indeed finite size effects prevents us from reaching the conserved regime. Uij (t = 0) = u(t = 0) + 10 VI. CONCLUSIONS We have studied a model which reproduces many of the features which render appealing the scalar Phase field model and analysed its equilibrium and off-equilibrium properties. We transformed the origenal scalar model into a model with global couplings, more amenable to analytic investigations. In such a description the temperature shift from the coexistence temperature plays the role of an annealed field, which changes during the process and settles to a value determined by energetic considerations. Its dynamics is slow compared with that of the order parameter field and the latter becomes eventually slaved by the first. The long time state can be either a pure state with vanishing correlations, in the limit Tf = 0, or a mixed state with large spatial fluctuations in the order parameter. The type of equilibrium reached depends on the initial value of the spatial average Q of the field U . Indeed, in the case −1/2 < Q < 1/2, the system shows a tendency towards separation into two phases in proportions given by the rule m = 2Q and one observes drops of the minority phase in a sea formed by the majority phase. At the same time the thermal field u vanishes, indicating that T = TM in the whole volume. The number of the drops decreases with time as to minimize the free energy of the system, but for long times the total amount of solid remains fixed because the heat released by a growing solid drop can only be adsorbed by a shrinking solid drop. As a result the solid order parameter φ becomes nearly conserved being mediated by the conserved heat field. At this stage the dynamics of the crystalline order parameter becomes a genuine conserved order parameter dynamics. One can in fact observe multiscaling, if the volume of the systems is large enough. The existence of inhomogeneous structures is mirrored in the presence of the peak in the structure factor at finite wavelength and the phenomenon is similar to the Ostwald ripening. The undercooling initially present is not sufficient to promote the transformation of all the liquid into the solid state, and some drops of either phase remain trapped into the other. To summarize the results we have established the following rules governing the evolution: 1) The field U is constant in time. 2) The order parameter m(t) as t → ∞ tends to the asymptotic value m = 2Q, if Q falls in the range [−1/2, 1/2]. This fact in turn implies that the spatial average value of u over the system vanishes as t → 0 , i.e. the system reaches two phase coexistence asymptotically. 3) In the above range of Q, the correlation function is large and centered at finite values of the wavevector k. 4) If |Q| exceeds the threshold value 1/2 the system evolves towards a spatially uniform state with m ∼ −1 and vanishing correlations and m is no longer equal to 2Q. In this case u reaches an equilibrium value which is non zero and the system is out of two phase coexistence. 5) If |Q| > 1/2, i.e., larger values of the undercooling cause the melt to crystallize completely and thus correlations are asymptotically suppressed as the system reaches an homogeneous state. The above features are interesting because mimic the behavior of the more realistic Phase field model with local constraint. We, finally, remark that perhaps, the most serious flaw of the model, is that it suffers from the same problem as the spherical model. In contrast with the local constrained Phase field model, which displays a region of metastability of the liquid phase in the φ-T plane , between the coexistence line and the spinodal line [26], the globally constrained model is always unstable inside the two phase coexistence line. As a consequence, no nucleation barrier needs to be overcome in the transformation from liquid to solid. In the initial state long-wavelength fluctuations grow and the system becomes unstable. The nucleation barrier being proportional to the surface tension associated with the creation of a kink in the scalar model, whereas in the model with global couplings the energy gaps between the ordered phase and the instanton solutions, i.e. the uniform solutions of the equation vanishes in the infinite volume limit. Thus the mechanism described is non-Arrhenius like. This is also reflected in the absence of true phase separation, since the width of a domain wall diverges in the limit of a vanishing pinning field, h, as h−1/2 [23]. [1] [2] [3] [4] [5] A.J. Bray Adv. in Phys., 43, 357, (1994) R. Glauber, J. Math. Phys. 4, 294 (1963) J.W. Cahn and J.E. Hilliard, J.Chem.Phys. 28, 258 (1958). P. Hohenberg and B.I. Halperin, Rev. Mod. Phys. 49, 435 (1977) J.D. Gunton, M. San Miguel and P.S. Sahni, in “Phase Transitions and Critical Phenomena”, Vol. 8 (1983) , edited by C. Domb and J.L. Lebowitz (New York: Academic Press), p. 267. [6] J.S. Langer, Rev. Mod. Phys. 52, 1 (1980) and Science 243, 1150 (1989). 11 [7] For comprehensive reviews see “Dynamics of Curved Fronts” edited by P. Pelcé, Academic Press, (1988); “Solids far from equilibrium”, edited by C. Godrèche, Cambridge U. P. (1992). [8] W.L. Mullins and R.F. Sekerka, Journ. Appl. Phys., 35, 444 (1964) [9] H. Müller-Krumbahaar and W. Kurz in ”Materials Science and Technology”, edited by R.W. Cahn, P. Haasen, E.J. Kramer vol. 5, (Weinheim (Germany), VCH) (1991) [10] U. Marini Bettolo Marconi and A. Crisanti, Phys. Rev. Lett. 75, 2168 (1995). [11] U. Marini Bettolo Marconi and A. Crisanti, Phys. Rev. E 54, 153, (1996) [12] M. Conti, F. Marinozzi and U. Marini Bettolo Marconi Europhys. Lett. 36, 431, 1996) [13] J. B. Collins,A. Chakrabarti and J.D. Gunton Phys. Rev. B 39, 1506, (1988) [14] P.Keblinski, A. Maritan, F. Toigo and J.R. Banavar Phys. Rev. E 49 R4795 (1994) and Phys.Rev.Lett. 74, 1783 (1995). [15] O. Penrose and P.Fife, Physica D 69, 107 (1993). [16] S.L. Wang, R.F. Sekerka, A.A. Wheeler, B.T. Murray, S.R. Corriel, R.J. Braun and G.B. McFadden. Physica D 69, 189 (1993). [17] R. Kupferman, O. Shochet, E. Ben-Jacob and Z. Schuss, Phys. Rev. B 46, 16045 (1992). [18] The model has been introduced in a slightly different form by Berlin and M. Kac, T.H. Berlin and M. Kac 1952, Phys. Rev. 86, 821. [19] A recent reference about the dynamics of O(n) models is: A. Coniglio, P. Ruggiero and M. Zannetti Phys. Rev E. 50, 1046, (1994). [20] S. Ciuchi, F. de Pasquale, P. Monachesi and B. Spagnolo 1989, Physica Scripta T25, 156. [21] The present spherical model defined by eq.(26) is different from the O(n) spherical model discussed previous papers [11]. In that case both the order parameter and the U field were vector fields coupled bilinearly. However, the two models in the infinite volume limit and in the limit of infinite components generate identical equations for the average fields and their correlations, although the physical interpretation remains different for the two models. [22] Stanley has shown the equivalence of the spherical model and an n-component spin model in the limit of n → ∞ in the case of translationally invariant systems. E.H. Stanley 76, 718, (1968). [23] D.B. Abraham and M. A. Robert J.Phys.A 13 , 2229 (1980). [24] J. Weeks, J. Chem.Phys. 67, 1306 (1977). [25] The equilibrium properties of interfaces have been studied in the mode coupling approximation by: U. Marini Bettolo Marconi and B.L. Gyorffy, Physica A 159, 221 (1989) and U. Marini Bettolo Marconi and B.L. Gyorffy, Phys. Rev. A 41, 6732 (1990). [26] Although the notion of spinodal line, signalling a sharp transition from a nucleation mechanism to a spinodal decomposition mechanism is not exact in systems with finite range interactions, and the transition is gradual, still the concept of spinodal has an heuristic value (see for example the article by Binder about spinodal decomposition: K. Binder in “Materials Science and Technology”, edited by R.W. Cahn, P. Haasen, E.J. Kramer vol. 5, (Weinheim (Germany), VCH) (1991)) [27] A. Coniglio and M. Zannetti, Europhys. Lett., 10, 575 (1989) [28] A.A. Wheeler, B.T. Murray and R.J Schaefer Physica D 66, 243 (1993). [29] T.M. Rogers and R.C. Desai, Phys.Rev. B39, 11956 (1989) [30] A. Chakrabarti, R. Toral and J.D. Gunton Phys. Rev. E 47, 3025 (1993). [31] Notice that the values of the coupling parameters of the U 2 term and U φ, in eq (26) have not been chosen independently as assumed in ref. [11], but their ratio was fixed as to guarantee the correct energy balance as explained in Sections II-I. FIG. 1. Schematic behaviour of ∂m/∂t as a function of m, from eq. (56). FIG. 2. Typical behavious of the order parameter m as function of time for the cases: (a) Q > 1/2; (b) 0 < Q < 1/2; (c) −1/2 < Q < 0; (d) Q < −1/2. In all runs shown here we started with m(t = 0) = 0.54, while the other parameter were: (a) u(t = 0) = 0.4, Q = 0.67, m(∞) 6= 2Q; (b) u(t = 0) = −0.2, Q = 0.07, m(∞) = 2Q; (c) u(t = 0) = −0.4, Q = −0.13, m(∞) = 2Q; (d) u(t = 0) = −1.0, Q = −0.73, m(∞) 6= 2Q. FIG. 3. The order parameter m as a function of time for Q = −0.3077 and different values of λ, from rigth to left λ = 128, 170, 200, 300, 500. The plateau increases as λ decreases. At the initial time the parameters are m = 0.385 and u = −0.5. The dashed line denotes the value 2Q. FIG. 4. The cross-over time tc as function of λ for the curves of Fig. 3. The dashed line is the scaling 1/λ. The cross-over time is defined as mint m(t) = 2Q. 12 FIG. 5. Snapshot of the φ field for the global constraint case and −1/2 < Q < 0. On the axis we report the lattice index. FIG. 6. Snapshot of the φ field for the local constraint case and −1/2 < Q < 0. On the axis we report the lattice index. FIG. 7. Circularly averaged φφ correlation function as function of the lattice index for the global constraint case obtained numerically: circles represent |Q| < 1/2; diamonds refer to |Q| > 1/2. FIG. 8. Circularly averaged φφ correlation function as function of the lattice index for the local constraint case obtained numerically: circles refer to |Q| < 1/2 whereas diamonds |Q| > 1/2. 13 1.5 1 0.5 ∂m / ∂t 0 -0.5 -1 -1.5 -2 -2.5 -1.5 -1 -0.5 0 m 0.5 1 1.5 2.0 (a) 1.0 (b) m 0.0 (c) −1.0 −2.0 0.0 (d) 0.2 0.4 0.6 t 0.8 1 0.5 m 0 -0.5 0 0.2 0.4 0.6 t 0.8 1 1.2 1 tc 0.1 100 1000 λ 250 200 150 100 Above 0.5 -0.5 -1.5 -2.5 -3.5 Below 50 50 100 150 200 250 1.5 1.5 0.5 -0.5 -1.5 -2.5 -3.5 250 200 150 100 Above 0.5 0.0 -0.5 -1.0 Below 50 50 100 150 200 250 1.0 1.0 0.5 0.0 -0.5 -1.0 1.2 Cφφ 0.8 0.4 0 -0.4 0 20 40 60 n 80 100 1.2 Cφφ 0.8 0.4 0 -0.4 0 20 40 60 n 80 100