Crystal growth from a supersaturated melt: relaxation of the solid-liquid dynamic stiffness
CCrystal growth from a supersaturated melt: relaxation of the solid-liquid dynamicstiffness
Francesco Turci and Tanja Schilling Theory of Soft Condensed Matter, Physics and Materials Science Research Unit,Universit´e du Luxembourg, L-1511 Luxembourg, Luxembourg (Dated: October 11, 2018)We discuss the growth process of a crystalline phase out of a metastable over-compressed liquidthat is brought into contact with a crystalline substrate. The process is modeled by means ofmolecular dynamics. The particles interact via the Lennard-Jones potential and their motion islocally thermalized by Langevin dynamics. We characterize the relaxation process of the solid-liquidinterface, showing that the growth speed is maximal for liquid densities above the solid coexistencedensity, and that the structural properties of the interface rapidly converge to equilibrium-likeproperties. In particular, we show that the off-equilibrium dynamic stiffness can be extracted usingcapillary wave theory arguments, even if the growth front moves fast compared to the typicaldiffusion time of the compressed liquid, and that the dynamic stiffness converges to the equilibriumstiffness in times much shorter than the diffusion time.
I. INTRODUCTION
When two phases of a material enter into contact, theyform an interface characterized by a typical free energycost. The interface free energy , surface tension or inter-facial tension is a global measure of such a cost. It en-codes information on the shape, the local curvature, andthe roughness of the region that distinguishes the twophases. Strictly speaking, an interfacial tension is prop-erly defined only if the two phases are coexisting and atthermodynamic equilibrium: in this case we can define acanonical ensemble of equilibrium interfaces from whichan interfacial free energy contribution can be computed.On the theoretical side, much work has been done onthe equilibrium properties of interfaces. While liquid-vapor systems have been largely discussed and under-stood [1–5], liquid-liquid and solid-liquid interfaces arestill under investigation. For equilibrium solid-liquid in-terfaces, even the very simplified cases of hard-spheresor Lennard-Jones particles are still subject of extensivestudies aimed at the correct evaluation of equilibriumproperties in connection with coarse-grained theories [6–9].Extending the concept of equilibrium interfacial ten-sion to non-equilibrium is not straightforward. In fact,while it is always possible to define an interface betweentwo distinct phases even during non-equilibrium dynam-ics (both steady states and relaxation processes), it isnot in general clear under which conditions the interfa-cial structural properties are actually related to a freeenergy cost or rather need to be seen from the purelykinetic point of view. Formally, it is always possible todefine an effective interfacial tension [10, 11], for exam-ple starting from the anisotropy of the pressure tensor,from the composition gradient across the interface or, aswe shall see in the following, from the fluctuations of theshape of the interface itself. There is simulation and ex-perimental evidence that the description of growth phe-nomena requires the introduction of an effective dynamic tension in order to properly stabilize the interface duringthe transient states for fluid-fluid interfaces [12–14], butthere still is a substantial lack of consensus on the deter-mination of the proper relaxation times related to it andtheir interpretation [15–18].In the present study we concentrate our attention onsolid-liquid interfaces and the relaxation process associ-ated with crystallization. As compared to relaxing liquid-liquid interfaces, crystal growth poses an additional chal-lenge due to the elastic and viscous properties of thegrowing, off-equilibrium solid phase. While scaling the-ory can be useful in order to provide general argumentsabout kinetic models underlying solid-on-solid growth[19–23], we devote our study to the analysis of the cap-illary effects related to the growth processes and theircounterparts at equilibrium, in order to characterize thecrossover from the dynamic to the equilibrium interfacialtension for solid-liquid interfaces.Our study is organized as follows: in section II we re-call some fundamental concepts related to capillary wavetheory (CWT) at equilibrium and out of equilibrium; insection III we present our model for heterogeneous crystalgrowth and provide an illustration of the growth kinet-ics; in section IV we proceed to the description of the dy-namic interfacial properties, presenting two examples ofrelaxation and two independent means of numerically de-riving estimates for the dynamic stiffness based on CWT;we conclude in section V with the discussion of the re-sults in comparison with the equilibrium features of thesolid-liquid interface. II. CAPILLARY WAVE THEORY ANDSTIFFNESSES
For liquid-liquid interfaces in thermodynamic equilib-rium, shape flucutations are described by capillary wavetheory (CWT)[1, 24–27]. CWT is derived by applyingthe equipartition theorem to the distribution of Fouriermodes of the interfacial height fluctuations [1, 24–27]. Al- a r X i v : . [ c ond - m a t . s o f t ] M a y beit thought and developed for gas-liquid or liquid-liquidinterfaces, this simple theory has also been successfullyapplied to the analysis of solid-liquid interfaces [28–30].The main prediction of the theory, directly related toequipartition of the thermal energy over the possible de-grees of freedom, is the scaling of the height-height cor-relation function h ( (cid:126)r ) with the (two dimensional) radialposition (cid:126)r = x ˆ e x + y ˆ e y . This correlation function de-scribes the local deviations of the interfacial height fromthe mean interface position. In general, for crystal-liquidinterfaces, the structure of the crystal is reflected in ananisotropic form of the correlation function. Yet, forcrystal orientations with x, y symmetry, in the limit ofshort wave-vectors (cid:126)q one can write the power spectrumof height fluctuations as: (cid:104)| h eq ( q ) | (cid:105) = k B TL ˜ γq (1)where q = | (cid:126)q | , ˜ γ is the stiffness of the interface, T isthe temperature, k B is Boltzmann’s constant, L is thelateral size of a system of square cross section L × L and (cid:104)·(cid:105) denotes the equilibrium ensemble average. It istherefore sufficient to extract the power spectrum of theheight fluctuations in order to estimate, through a bestfit at low wave-vectors, the equilibrium stiffness of theinterface.Note that the stiffness, which is orientation dependent,is related to the interfacial tension by the following ten-sorial relation:˜ γ xy ( n z ) = γ ( n z ) + ∂ γ ( n ) ∂ n x ∂ n y (cid:12)(cid:12) n = n z (2)where n x , n y , n z correspond to the (3d) principal axes of ∂ γ ( n ) /∂ n i ∂ n j | n = n z . Such a relation allows to comparestiffnesses and tensions, once several crystal orientationsare taken into account [7, 31].One can define the temporal dynamical correlationfunction expressed in terms of the spatial height-heightcorrelation function as g h (∆ t ) = (cid:104) h (cid:48) ( (cid:126)q, t ) h (cid:48) ( (cid:126)q, t + ∆ t ) (cid:105) (3)where h (cid:48) ( (cid:126)q, t ) = h ( (cid:126)q, t ) − (cid:104) h ( (cid:126)q, t ) (cid:105) t is the deviation fromthe time-averaged profile. This two-time correlationfunction is characterized by a typical decorrelation capil-lary time τ decorr , that is related to the relaxation of theinterface close to equilibrium and that determines thedynamics of the modes on the interface.Integrating the spectrum in eq. 1 between a minimumcutoff length l cut and the linear size of the cross section L ,and applying a convolution approximation[32], one findsthe scaling of the mean interfacial width: w = w + k B T γ ln( L/l cut ) (4)where w refers to an hypothetical intrinsic interfacialwidth, defined beyond any roughening due to capillary effects [24]. While the intrinsic interfacial width is hardlymeasurable, this expression allows for the computation ofthe interfacial stiffness ˜ γ by means of finite size scalinganalysis.CWT is an appropriate theory for rough liquid-gasand liquid-liquid interfaces: in the case of liquid-solidinterfaces, the applicability of the same theory is notas clear, since it makes the strong approximation thata solid, crystalline phase can be considered as anextremely viscous fluid phase, with consistent stiffnessesand capillary times [33]. Yet, it has been numericallyshown [34] that the logarithmic scaling of CWT can beeffectively used as a means to correctly estimate theinterfacial stiffness for solid-liquid mixtures [6, 35].When two phases are put into contact, a relaxationprocess of the interfacial properties is initiated, whicheventually converge into the properties of the equilib-rium interface as it is described by capillary wave theory.The shape of the interface during the relaxation processis partly determined by the nature of the driving forcethat creates the non-equilibrium state. The specific formof the power spectrum of the height-height fluctuationsobeys different scaling laws in the case of sharp interfaces,linear gradients of concentration or temperature, and in-terfaces moving at constant velocities [36–38]. In ref. [36]the spectrum of non-equilibrium stationary height fluc-tuations is derived for an interface that is moving at con-stant velocity v . It converges in the limit of small veloc-ities v →
0, to (cid:104)| h noneq ( q ) | (cid:105) trajectories = k B TL ˜ γ ( q + 4 π /l ) (5)where l is a crossover length due to the chemical potentialgradient and the latent heat of the liquid-to-solid transi-tion, that has a pinning effect on fluctuations occurringon scales larger than l . The consequence of the drivingis thus, in this limit, similar to the effect of gravity onliquid-liquid equilibrium interfaces.As suggested by equation 5, in the case of a station-ary state, a non equilibrium dynamic surface tension canbe derived directly via a fit of the spectrum of heightfluctuations. We use this procedure in the case of a re-laxation process driven by the chemical potential differ-ence; we follow the time dependence of the spectrum andextract from it an estimate for the effective dynamic ten-sion. This is an alternative to the computation of gra-dients and anisotropies across the interface, which alsodefine dynamic stiffnesses, such as in the case of compo-sition gradients in Kortweg’s theory of effective surfacetension [39], the pressure tensor anisotropy [18] or non-equilibrium linear response to excitations [40, 41].In the following, we will analyze different relaxationdynamics related to the crystallization of a Lennard-Jones liquid in contact with a crystalline substrate us-ing CWT and fluctuating hydrodynamics in order to ex-tract the dynamic behavior of the crystal-liquid interfa-cial stiffness that is related to the tension by eq. 2. III. GROWTH PROCESSA. A model for crystal growth on a substrate
We consider a simple model of out-of-equilibrium dy-namics modeling crystal growth and the motion of asolid-liquid interface and its capillary mode fluctuations.We perform Molecular Dynamics simulations of crystal-lization of an over-compressed liquid that is in contactwith a substrate, both composed of identical particles ofmass m interacting via a truncated and shifted Lennard-Jones potential u ( r ) = u LJ ( r ) − u LJ ( r cut ) for r < r cut and 0 otherwise, with u LJ ( r ) = 4 (cid:15) [( σ/r ) − ( σ/r ) ] and r cut = 4 σ .We first prepare separately an over-compressed,metastable liquid at high density and an immobile defect-free crystal of particles exposing the [100] orientation tothe liquid perpendicular to the z axis of the simulationbox, with unit cell length a = (cid:112) /ρ sol where ρ sol ( T ) isthe coexistence density for the solid phase at a given tem-perature (see Appendix A for the simulation details). Weconsider liquid densities ρ ranging from about 0 . ρ sol ( T )up to maximum 1 . ρ sol ( T ) avoiding configurations wherecrystallites are present due to homogeneous nucleation.The system undergoes isothermal and isochoricLangevin dynamics – isochoric, to avoid the additionalcomplication due to effects of a barostat on the dynam-ics. This implies that the density of the over compressedliquid does not remain constant during the growth of thecrystal: the interface progresses along the z-axis drivenby a time-dependent chemical potential difference thatgradually diminishes, eventually leading to a conditionfor which the liquid reaches the equilibrium coexistencedensity. B. Growth Profiles and Growth Speeds
In order to distinguish the crystal from the liquid phasewe use a set of observables that measure the degree oflocal order in the neighborhood of every particle. Weuse the so called locally averaged bond order param-eters ¯ q , ¯ q , inspired by Steinhardt’s order parameters[42] and defined in [43] (see the Appendix for more de-tails). These order parameters take into account the firstand the second shell of neighbours around a given parti-cle: this allows for the distinction between different crys-talline structures (face centered cubic, body centered cu-bic, hexagonal close packed). Using these observables,we follow the phase transition from the over-compressedliquid, tracking the gradual increase of crystalline order.As shown in fig. 1, the growth of the crystal involves theformation of a highly ordered region of fcc structures sep-arated from the remaining liquid particles by few layersof hcp-like particles. However, the growing crystal re-laxes only slowly towards its equilibrium configuration.At any time, there are extended clusters of impurities(green regions), a signature of the fact that the relax- t=30t=10 t=120t=150 FIG. 1. (color online) Four snapshots representing the time-evolution of the crystal growth on top of a perfectly orderedtemplate (gray particles) from a lateral side point of view.The density of the liquid is ρ/ρ sol = 1, the temperature is T = 0 . (cid:15) and the system has dimensions 25 × × a ,with time measured in units of t LJ := (cid:112) mσ /(cid:15) . The degreeof crystallinity, determined using the ¯ q , ¯ q order parameters,blends from essentially fcc-like (red) to liquid (blue) througha 2-3 layers thick interface of hcp-like particles. ation inside the crystal proceeds much more slowly thanthe formation of new crystalline layers.The local order observables permit to identify the in-terface and describe its dynamics. It is sufficient to trackthe coarse-grained laterally averaged (cid:104) ¯ q ( z, t ) (cid:105) x,y profilein order to determine the position and the width of theinterface (we denote with (cid:104)·(cid:105) x,y the average over lateralcross sections of the system). The profile is modeled witha hyperbolic tangent function (cid:104) ¯ q ( z, t ) (cid:105) x,y = C + C tanh (cid:18) z − z ( t ) w ( t ) (cid:19) (6)which provides the values of the position z ( t ) and thewidth w ( t ) via a four parameters best fit procedure (seefig. 2).In order to improve the statistics, we consider a sampleof 100 independent crystallization trajectories, where werestart the simulation from different equal-energy config-urations of the metastable liquid. Density profiles, localorder parameter profiles, interface trajectories and alikequantities result from a same-time average over such asample, with time measured in Lennard-Jones units.We use the crystallinity profiles defined by eq. 6 inorder to track the position and the speed of the growingcrystal front. In fig. 3 we show average interface trajecto-ries z ( t ). Depending on the initial difference in densitiesbetween the substrate and the over-compressed liquid,saturation of the interface position to its final value isreached at different times after an initial linear growthregime. We are interested in the characteristic initial ve-locity v , when the system is far from equilibrium anddriven by a large chemical potential difference. For thelargest over-compressions the initial linear growth regimelasts sufficiently long that the crystal completely fills thebox.The crystallinity profiles shown in fig. 2 can be com-pared with the laterally averaged density profiles ¯ ρ ( z ) = (cid:104) ρ ( x, y, z ) (cid:105) x,y . While the local bond order parametermonotonically connects the solid to the liquid region, theaveraged density profile is marked by the presence of adepletion zone in the vicinity of the interface (similarlyto what has been observed in equilibrium for both hardspheres [44] and Lennard-Jones particles [45]), which ispartly compensated by a higher density of the first crys-tal layers immediately in contact with the liquid. Thedepth of the depletion zone depends on time, while thelocal bond order profile simply translates along the z axis.The change in density spatially precedes the change in ¯ q ,effectively determining two interface positions that differby approximately one unit cell.By fitting the profiles we obtain the initial growthvelocity as a function of the initial liquid density (seefig. 4). The initial growth velocity has a maximum atliquid densities slightly higher than the coexistence soliddensity ρ sol , in contrast to what has been reported inthe case of hard-spheres [44]. For both temperaturesconsidered here, the maximum is reached at liquid den-sities such that ∆ ρ/ ∆ ρ coex ≈
2, with ∆ ρ/ ∆ ρ coex =( ρ − ρ liq ) / ( ρ sol − ρ liq ), where ρ liq , ρ sol indicate the equi-librium coexistence densities for the liquid and the solidphase respectively. At even higher densities the com-petition between the increasing driving force due to thechemical potential difference and the dramatic decreaseof the diffusion constant [46] leads to a decrease of thegrowth velocity. Yet, compared to the hard-spheres case,the decrease is much less pronounced. Notice that thegrowth velocity is generally very fast compared to thetypical diffusion time: for example, the peak velocity for T = 0 . (cid:15) is v peak ≈ . a/t LJ = 0 . σ/t LJ , to becompared with the diffusion time τ D = σ /D ≈ t LJ and velocity v D = σ/τ D ≈ . σ/t LJ . (As the homoge-neous nucleation rate rapidly increases with the degreeof over-compression, we end our analysis at a density forwhich the velocity is still about 60% of the peak velocity.)A simple model that reproduces some qualitative fea-tures of the growth velocities is the classical Wilson-Frenkel law [47–49]. It results from an estimated balancebetween the rate of activated attachment and detachmentof the liquid particles on the crystal surface v ( T ) = A kin (cid:20) − exp (cid:18) − g sol − g liq k B T (cid:19)(cid:21) (7)where A kin is a kinetic prefactor proportional to the dif-fusion coefficient D , and g sol , g liq are the solid and liq-uid free-energies respectively. Using accurate simulationdata for the Lennard-Jones equations of state and coexis-tence densities [50], we plug into eq. 7 the correspondingvalues for the free energy differences as a function of theliquid density and, using the diffusion coefficient com-puted by Molecular Dynamics simulation of large, liquid . . . ¯ ρ ( z ) , ¯ q ( z ) t = 20 t LJ t = 400 2 4 6 8 10 12 14 16 z/a . . . ¯ ρ ( z ) , ¯ q ( z ) t = 70 0 2 4 6 8 10 12 14 16 z/at = 160 FIG. 2. (color online) Time evolution of the coarse-grainedlaterally averaged density (red continuous line) and ¯ q (dashedblue) profiles, normalized in such a way that the perfect crys-tal phase has value 1 and the bulk liquid takes value 0. Theblack dotted line corresponds to the fit of the ¯ q profiles ac-cording to eq. 6 in order to determine the position of theliquid-solid interface and its width. Notice that the ¯ q profileshows two interfaces (wall-crystal and crystal-liquid respec-tively): we concentrate our study on the properties of themoving crystal-liquid interface. bulk systems, we obtain the red dots of figure 4 (rescaledsuch that the lowest density point coincides with thevelocity curves from simulation). We observe that theWilson-Frenkel law does not allow to reproduce the largeinitial increase in growth speed. In fact, the resulting es-timate for the speed is largely dominated by the decreaseof the diffusion coefficient at large densities, which is notcompensated by the chemical potential difference term.The failure of the Wilson-Frenkel approach is not sur-prising, since it is based on bulk, close to equilibriumproperties, and it has already been demonstrated [51]that the specific, non-diffusive dynamics at the interfacemay lead to understimations of the actual velocity of atleast one order of magnitude. More than that, one shouldalso consider that the equilibrium free energies can onlybe a rough estimate of the actual driving force acting be-tween the two phases of the present model, where somerelaxation mechanisms occur on time scales far beyondthe actual duration of the simulation. Finally, the for-mation of the initial layers is very fast, with barely nonucleation time related to it, thus the activation barrierseems to be suppressed [52]. ρ liq = 0.96 σ -3 z / σ t / t LJ FIG. 3. Average interface position (extracted from best fits ofeq. 6 to the ¯ q profiles) for a system at temperature T = 1 . (cid:15) and solid coexistence density ρ sol = 1 . σ − for differentinitial liquid densities ρ liq . Note that the speed of the interfaceis non-monotonic in T , i.e. the curve that corresponds to T = 1 . (cid:15) lies below the T = 1 . (cid:15) curve. T=0.7793 ε L=20 a t/t LJ z a
100 50 v t L J / a Δρ / Δρ coex FIG. 4. Velocity of the crystal-liquid interface as a functionof the ratio ∆ ρ/ ∆ ρ coex = ( ρ − ρ liq ) / ( ρ sol − ρ liq ) for differentdensities ρ of the metastable, over-compressed liquid. Thegray area corresponds to the range of coexistence densities,above which the crystalline phase is the only stable phase.Inset: evolution of the position of the interface as a functionof time for ρ/ρ sol = 1 , T = 0 . (cid:15) for different linear lengths L x = L y = 20 a, a, a, a : finite size effects play only amarginal role in the determination of the growth velocity. IV. NON-EQUILIBRIUM INTERFACIALSTIFFNESSA. Small chemical potential difference
If the chemical potential difference between the liquidand the crystal is small from the beginning, the growingcrystal rapidly depletes the liquid, and the motion of the interface slows down quickly. Under these conditions wecan observe the relaxation of the interfacial propertiesinto their equilbrium values on the scale of our simula-tion.We set T = 1 . (cid:15) , the initial liquid density ρ =1 . σ − , and ρ sol = 1 . σ − . Using the equation ofstate fitted in [50] we estimate the initial chemical po-tential difference between the solid and the liquid phaseto be µ sol − µ = − . (cid:15) per particle. Under these con-ditions, growth stops after approximately 16 crystallinelayers are formed. We simulated systems of squared crosssection L × L with L = 25 , , , , , , , a andlongitudinal length L z such that we always have 64 fcclayers between the two initial interfaces. We extractedthe time evolution of w ( t ) (eqn. 6) averaged over 100trajectories, each of which represents the growth of twoopposite crystal fronts that we assume to be independent(which is true if the liquid region between the interfacesremains much larger than the typical width of the inter-face and the correlation length in the liquid). Regardlessof the area of the cross section, after a rapid initial in-crease, the interfacial width saturates.The first remarkable feature is that the whole processis extremely fast: if we call τ D = σ /D the self-diffusiontime of a particle in the over-compressed liquid, the inter-face attains its final position z in a time of about 2 − τ D .Thus the growth process is not dominated by the diffu-sion of liquid material in the vicinity of the solid-liquidinterface, but by an adsorption mechanism.Before we analyse the evolution of the interfacialwidth, we briefly discuss the adsoprtion mechanism.Fig. 5 shows the joint probability distribution P (¯ q , ∆¯ q )for the difference ∆ q = ¯ q ( t + ∆ t ) − ¯ q ( t ) for differ-ent times t , with a time interval of ∆ t = 10 t LJ ≈ τ D / q ≈ .
16 while the solid region shows the fcc peak andthe signatures of the growing crystal around ¯ q = 0 . q = 0, because the majority of the particles are eitherin the liquid or in the crystal and do not change phase.Particles in the bulk, that have exceptionally low/highvalues of ¯ q are likely to move closer to the average bulkvalue of their phase. This is the origin of the oblique tiltin the distributions. However, particles that lie close tothe interface and are incorporated into the crystal show adifferent dynamics and are responsible for the main dif-ferences that one can observe between the four panels:At early times the distribution is strongly asymmetric,particles with 0 . < ¯ q < . t , while the reverse process is suppressed (blue dip inthe bottom-centre of the first panel). This asymmetryis significantly reduced already 2∆ t later at t = 30 t LJ and there is no noticeable difference between the last twopanels, showing that the local order around interface par-ticles fluctuates in a quasi-equilibrium manner alreadyvery early. FIG. 5. (color online) Time evolution of the joint probability distribution P (¯ q , ∆¯ q ) for a system of linear size L = 80 a atdensity ρ/ρ sol = 0 .
96 and temperature T = 1 . (cid:15) , computed for a time-interval of ∆ t = τ D / L=20 a
30 4080 T=1.967 ερ =1.08 σ -3 w ( t ) / σ t / t LJ
10 100
FIG. 6. (color online) Time evolution of the square of theinterface width determined from an hyperbolic tangent fit ofthe ¯ q profiles as a function of time for different representativesystem sizes at density ρ/ρ sol = 0 .
96 and temperature T =1 . (cid:15) . We conclude therefore that while at very early timesthe presence of the template accelerates the formation ofnew crystalline layers in such a way that typically liq-uid particles (¯ q ≈ .
16) are converted into solid parti-cles (∆¯ q ≈ = +0 .
2) but never converted back, the latergrowth dynamics (after about 1 τ D ) consists of a bal-anced, close to equilibrium exchange between liquid andsolid particles in the vicinity of the interface.Now we return to the analysis of the width. At equi-librium eq. 4 has to hold, which defines an equilibriumstiffness ˜ γ eq . This suggests to define a dynamic stiffness ˜ γ ( t ) extracted from the scaling of w L ( t ) with the laterallinear size L at a time t . In fig. 7 we see, for the rangeof lateral sizes considered here, a logarithmic scaling ofthe width at any time. Based on this we compute ˜ γ ( t ) asshown in fig. 8. Ignoring the very early times (for whichan interface width cannot be defined using eq. 6), we seethat the dynamic stiffness starts from very high valuesand rapidly relaxes to an average value that we interpretas the equilibrium stiffness ˜ γ eq .To test this interpretation, we compute the stiffness w ( L ) / σ L/a
25 50 75
FIG. 7. Scaling of the the interface width with the system size(logarithmic scale on the horizontal axis) at different timesgoing from t = 10 t LJ to t = 190 t LJ , bottom to top. For easeof visualization, the curves are shifted by a constant amountproportional to t . from the height-height correlations. We construct thediscretized height function h ij by mapping interface par-ticles onto a regular grid interpolating by Shepard’smethod [30]. In detail, we use the local bond orderparameters in order to distinguish solid and liquid par-ticles and select only those particles that satisfy ¯ q > .
025 and ¯ q < .
37 and have a number of solid neigh-bors 5 < n sol <
12 [53]. For the largest cross sec-tion L × L = 6400 a this procedure identifies roughly3000 particles at the interface between the solid andthe liquid (see fig. 9). Then, the ( x i , y i ) coordinatesof the particles are cast onto a regularly spaced gridof spacing ∆ x = ∆ y = 2 σ on top of which the z grid j coordinate is obtained as z grid j = (cid:80) i w ij z i / (cid:80) i /r ij ,where the summation runs over particles neighboringwithin a shell of radius 4∆ x from the grid point j and r ij = (cid:112) ( x i − x j ) + ( y i − y j ) . One obtains then thediscretized height h grid = z grid − (cid:104) z grid (cid:105) and by a Dis-crete Fourier Transform one derives the power spectrum (cid:104)| h ( q ) | (cid:105) trajectories shown in fig. 10. γ σ / ε v / v D t/t LJ z ( t ) / σ t / t LJ FIG. 8. (color online) Time evolution of the interface posi-tion (empty circles) and the interfacial stiffness (red circles)as computed from the width scaling with the lateral systemsize in the case of small chemical potential difference. Thehorizontal dotted line denotes the equilibrium value of thestiffness (within an interval of confidence of 3 standard devia-tions), computed from the low q limit of the heigh-height cor-relation function. In the inset, the interface speed (green dots)compared with the diffusion speed (dashed lines) v D = D/σ :the two cross at approximately t ≈ t LJ .FIG. 9. (color online) Interface particles and interpolatedinterface structure during the crystal growth using the Shep-ard’s method. On a cross section of L x × L y = 6400 a roughly3000 particles (red dots) are isolated according to their degreeof crystallinity (see text) and are the pivot for the computa-tion of the height-height correlation function (cid:104)| h ( (cid:126)q ) | (cid:105) . At early times, the spectrum shows a large plateau atlow q and at higher q it decreases as q − α with α ≈ . q the 1 /q decayis established; the amplitude of the fluctuations (relatedto the broadening of the interface width) increases andalready at t = 70 t LJ (ca. 1.5 τ D ) the spectrum is hardlydistinguishable from the equilibrium spectrum. Yet, theinterface is still moving with a speed that is fast com-pared to diffusion, and new particles are added to thecrystal. This appears to affect only the very low q re-gion, as a remnant of the effect of the pinning force andthe initially flat profile.Fitting the equilibrium profile with eq. 1 we obtain anestimate for the equilibrium stiffness ˜ γ eq = 1 . (cid:15)/σ that we plot as a reference value also in fig. 8.The scaling of the widths suggests that the interfacerapidly evolves into a quasi-equilibrium regime. We un-derpin this observation by computing the typical relax-ation time of the capillary fluctuations. We sample thefunction g t (∆ t ) = (cid:104) h ( (cid:126)q, t ) h ( (cid:126)q, t + ∆ t ) (cid:105) trajectories (8)parametrized by the starting time t , and where the aver-age is performed on equal-time configurations of distincttrajectories. We consider several initial spectra at dif-ferent times t and follow the subsequent evolution of thecorrelation function g t (∆ t ) in the later time steps. Infig. 11 we show the decay of these correlation functions,compared to the decay of equilibrium fluctuations definedin 3: in particular, in fig. 11(b) we show, by rescaling the g t (∆ t ) onto the equilibrium g eq (∆ t ), that, while the ini-tial relaxation process is marked by faster decays (shortercapillary times), long before the velocity of the interfacedecreases below the diffusion speed, the capillary timealready converges to its equilibrium value. B. Large chemical potential difference
If the chemical potential difference is large, growth isarrested by the contact between the two interfaces thatmove in opposite directions. Up to that point they movewith a large, almost steady velocity, and the growthmechanism is strong adsorption, see fig. 12. In the fol-lowing we explore the evolution of the interfacial height-height correlation for the fast, initial part of the growthprocess.Using the same techniques as applied in the case ofsmaller driving forces, we track the interface widths for asystem at temperature T = 0 . (cid:15) with ρ = ρ sol =0 . σ − (see fig. 13) corresponding to a theoreticalchemical potential difference µ sol − µ = − . (cid:15) per par-ticle. Exploiting the system size scaling, we compute thedynamic surface stiffness (fig. 14). Similarly to the pre-vious smaller ∆ µ case, we compare the resulting stiffnesswith the one obtained from the capillary spectrum of anequilibrated system at the same temperature. Again, af-ter an initial rapid drop, even at very early times the equilibriumt = 10 t LJ B T/L γ q 〈 | h ( q ) | 〉 / σ − − − − − q σ /2 π FIG. 10. (color online) Power spectrum of the height fluctu-ations for the model with small supersaturation. Notice therapid convergence to the equilibrium profile (continuous blackprofile). g ( Δ t ) / σ Δ t / t LJ equilibriumt=9 t LJ t=4t=2 A t g ( Δ t ) / σ Δ t / t LJ (a) (b) FIG. 11. (color online) Two-times correlation functions forthe small chemical potential regime: the equilibrium curve(black dots) is compared with the time dependent ones (a).In panel (b) a scale factor A t is used in order to compare thetrends. dynamic stiffness computed from the width scaling is inagreement with the equilibrium behavior.The signature of the non-equilibrium dynamics can befound only in the detailed shape of the height-heightspectrum (fig. 15): the plateau region at low q is alwayspresent, due to the non-negligible chemical potential dif-ference (corresponding to the high speed of the interface).The expected CWT scaling does not have sufficient timeto emerge. Thus we test how the stationary, low-velocitylimit of eq. 5 applies to the measured spectra. We per-form a single parameter fit using, as input, the stiffnessobtained from the L scaling of the widths extracted from¯ q ( z ) profiles. We get as a unique fitting parameter thecrossover length l the time evolution of which (far from FIG. 12. (color online) Joint probability distribution P (¯ q , ∆¯ q ) at initial times ( t = 10 t LJ ) and just before themerge of the two oppositely growing interfaces ( t = 70 t LJ )for the strong driving case. Notice that the distribution isstrongly asymmetric at both early and late times. Adsorpor-tion is therefore the main mechanism, rapidly transformingliquid particles into crystalline particles, with no backwardprocess. being stationary) is shown in the inset of fig. 14. Thiscrossover length brings the signature of the slowly re-ducing pinning force given by the chemical potential dif-ference, and it diverges in the limit of thermodynamicequilibrium when ∆ µ = 0 (cid:15) .Thus, the spectra (and the dynamic stiffness) revealthe non-stationary relaxation process of the growing crys-tal lattice. In particular, the fully non-equilibrium spec-tra are clearly different from the predictions of CWT orfluctuating hydrodynamics [36]. Yet, when dealing withglobally averaged quantities such as the interface width,it is still possible to observe the logarithmic scaling withthe system size, which permits to extract a dynamic stiff-ness still consistent with the equilibrium reference value. L=20 a
30 5080 ρ / ρ sol =1 T=0.7793 ε w ( t ) / σ t / t LJ
20 40 60
FIG. 13. Time evolution of the square of the interface widthdetermined from a hyperbolic tangent fit of the ¯ q profiles asa function of time for different representative system sizes atdensity ρ/ρ sol = 1 and temperature T = 0 . (cid:15) per particle.The maximum is reached when the interfaces begin to merge.The dashed lines are guides for the eye. t/t LJ l/ a γ σ / ε t / t LJ FIG. 14. Dynamic stiffness as a function of time for thestrongly driven system. The obtained values are in agree-ment with equilibrium values (gray-shaded area). In the inset:crossover scale l as a function of time in lattice constant units.Notice that throughout the relaxation, the crossover scale isfinite and below the system cross-section size L = 80 a . t = 10 t LJ 〈 | h ( q ) | 〉 / σ − − − − − q σ /2 π FIG. 15. (color online) Height-height correlation function (cid:104)| h ( q ) | (cid:105) for the growing crystal-liquid interface at ρ/ρ sol = 1 . and T = 0 . (cid:15) . The dotted lines represent the expressionin eq. 5 where we use the value of the interfacial stiffness ˜ γ independently obtained from the analysis of the widths of the¯ q ( z ) profiles. V. CONCLUSIONS
Simulating a simple model under isochoric, isother-mal conditions, we identify some characteristic featuresof crystallization from an undercooled liquid on a sub-strate: First, we have observed that the fastest initialgrowth speed is achieved when the liquid density is evenhigher than the crystalline template density. Then we have discussed two dynamical regimes at small and largedriving forces, fixed by the chemical potential differencebetween the substrate and the over-compressed liquid: atsmall driving forces we see that liquid particles get ini-tially adsorbed onto the growing crystal and, when theinterface velocity relaxes and becomes of the order of thediffusion speed, a quasi-equilibrium regime is attained.A dynamic stiffness can be defined using the scaling ofthe interface width. And the height-height fluctuationsrapidly converge to a CWT-like scaling regime. In thestrongly driven case, the spectra are very different fromthe CWT regime at any time but it is still possible tocompute a dynamic stiffness that is of the same order ofmagnitude of the equilibrium reference stiffness.The possibility of treating the driven dynamics withconcepts derived from CWT is in line with previous stud-ies of simplified models for liquid-liquid interfaces out ofequilibrium, such as under shear [54]. The fact that evenin a strongly driven dynamical regime the effective dy-namic stiffness appears to rapidly converge towards thecorresponding equilibrium value is particularly promis-ing for a coarse-grained description of the dynamic thatuses equilibrium properties in input, such as the DynamicDensity Functional Theory (DDFT) or, at a higher level,a Phase-field crystal approach[55]. DDFT has been suc-cessfully applied for nonequilibrium dynamics of hardspheres, for instance in the case of sedimentation [56],and the present results suggest that a proper choice ofthe functional form for the free energy functional couldreproduce the dynamics of the growing front, includingthe time evolution of the dynamic tension.
ACKNOWLEDGEMENTS
We thank Martin Oettel for insightful discussions.This project has been financially supported by the IN-TER/DFG/11/06 on Crystallization. Computer simula-tions presented in this paper were carried out using theHPC facility of the University of Luxembourg.
Appendix A: Simulation details
We simulated Lennard-Jones liquid particles at tem-perature T = 0 . , . , (cid:15) and variable densities ρ confined by two fcc walls of fixed particles. The wallshad a surface of A = 30 × a where a is the fcc latticespacing a = (cid:112) /ρ lattice , with ρ lattice = ρ sol ( T ) corre-sponding to the solid coexistence density at the giventemperature. The wall particles form three crystallinelayers per surface and interact with the liquid particleswith the identical interaction potential holding betweenliquid particles, with cutoff at r cut = 4 σ . This value isparticularly appropriate for the correct sampling of crys-talline structures. The distance along the z axis betweenthe two substrates is of at least 32 a , depending on the liq-uid density. We considered system sizes with lateral cross0sections raging from 25 a to 80 a , meaning approximately9 · to 9 · particles.We perform isochoric Langevin dynamics simulationswith a time step ∆ t = 0 . (cid:112) mσ /(cid:15) and friction coeffi-cient γ = 0 . t − and track the progression of crystal-lization using the averaged local bond order parameters¯ q , ¯ q proposed in [43]. Their definition requires the com-putation of the complex vector q l ( i ) q lm ( i ) = 1 N b ( i ) N b ( i ) (cid:88) j =1 Y lm ( r ij ) , (A1)where N b ( i ) corresponds to the number of nearest neigh- bors of particle i and Y lm ( r ij ) reads as the spherical har-monics. Averaging over the neighbors of particle i andparticle i itself¯ q lm ( i ) = 1˜ N b ( i ) ˜ N b ( i ) (cid:88) k =0 q lm ( k ) , (A2)and summing over all the harmonics we finally get¯ q l ( i ) = (cid:118)(cid:117)(cid:117)(cid:116) π l + 1 l (cid:88) m = − l | ¯ q lm ( i ) | . (A3) [1] R. Evans, Advances in Physics , 143 (1979)[2] M. Nijmeijer, A. Bakker, C. Bruin, and J. Sikkenk, TheJournal of chemical physics , 3789 (1988)[3] S. W. Sides, G. S. Grest, and M.-D. Lacasse, Phys. Rev.E , 6708 (1999)[4] I.-F. W. Kuo and C. J. Mundy, Science , 658 (2004)[5] M. Wertheim, The Journal of Chemical Physics , 2377(2008)[6] R. L. Davidchack, J. R. Morris, and B. B. Laird, Journalof Chemical Physics , 094710 (2006)[7] A. H¨artel, M. Oettel, R. E. Rozas, S. U. Egelhaaf, J. Hor-bach, and H. L¨owen, Phys. Rev. Lett. , 226101 (2012)[8] M. Oettel, S. Dorosz, M. Berghoff, B. Nestler, andT. Schilling, Phys. Rev. E , 021404 (2012)[9] X. Wang, J. Mi, and C. Zhong, Journal of ChemicalPhysics , 164704 (2013)[10] S. E. May and J. V. Maher, Phys. Rev. Lett. , 2013(Oct 1991)[11] W.-J. Ma, P. Keblinski, A. Maritan, J. Koplik, and J. R.Banavar, Phys. Rev. Lett. , 3465 (1993)[12] D. Truzzolillo, S. Mora, C. Dupas, and L. Cipelletti,Phys. Rev. Lett. , 128303 (2014)[13] P. Cicuta, A. Vailati, and M. Giglio, Applied optics ,4140 (2001)[14] C.-Y. Chen, C.-W. Huang, H. Gadˆelha, and J. A. Mi-randa, Phys. Rev. E , 016306 (Jul 2008)[15] J. R. Castrej´on-Pita, A. A. Castrej´on-Pita, E. J. Hinch,J. R. Lister, and I. M. Hutchings, Phys. Rev. E ,015301 (Jul 2012)[16] Y. D. Shikhmurzaev, IMA journal of applied mathemat-ics , 880 (2005)[17] J. R. Henderson, European Physical Journal-Special Top-ics , 61 (2011)[18] A. V. Lukyanov and A. E. Likhtman, Journal of ChemicalPhysics , 034712 (2013)[19] Z.-W. Lai and S. D. Sarma, Phys. Rev. Lett. , 2348(1991)[20] A.-L. Barab´asi, Fractal concepts in surface growth (Cam-bridge university press, 1995)[21] G. Gilmer and P. Bennema, Journal of Applied Physics , 1347 (2003)[22] J.-G. Yoon, H. K. Oh, and S. J. Lee, Phys. Rev. B ,2839 (1999) [23] S. Pal, D. Landau, and K. Binder, Phys. Rev. E ,021601 (2003)[24] F. P. Buff, R. A. Lovett, and F. H. Stillinger, Jr, Phys.Rev. Lett. , 621 (1965)[25] M. P. Fisher, D. S. Fisher, and J. D. Weeks, Phys. Rev.Lett. , 368 (1982)[26] D. Bedeaux and J. D. Weeks, Journal of ChemicalPhysics , 972 (1985)[27] K. Binder, M. Muller, F. Schmid, and A. Werner, Ad-vances in Colloid and Interface Science , 237 (2001)[28] J. Morris, Phys. Rev. B , 144104 (2002)[29] J. R. Morris and X. Song, Journal of Chemical Physics , 3920 (2003)[30] R. E. Rozas and J. Horbach, Europhysics Letters ,26006 (2011)[31] D. S. Fisher and J. D. Weeks, Phys. Rev. Lett. , 1077(1983)[32] K. Binder and M. M¨uller, International Journal of Mod-ern Physics C , 1093 (2000)[33] J. Hernandez-Guzman and E. Weeks, PNAS , 15198(2009)[34] T. Zykova-Timan, R. E. Rozas, J. Horbach, andK. Binder, J. Phys.: Condens. Matter , 464102 (2009)[35] J. Hoyt, M. Asta, and A. Karma, Phys. Rev. Lett. ,5530 (2001)[36] A. Karma, Phys. Rev. E , 3441 (1993)[37] A. Vailati and M. Giglio, Phys. Rev. E , 4361 (1998)[38] P. Cicuta, A. Vailati, and M. Giglio, Phys. Rev. E ,(2000)[39] K. D, Arch. Neerlandaises Sci. Exactes Naturelles , 1(1901)[40] G. Loglio, U. Tesei, N. Innocenti, R. Miller, and R. Cini,Colloids and Surfaces , 335 (1991), ISSN 0166-6622,a collection of papers presented at the XIth. EuropeanConference on the Chemistry of Interfaces[41] E. I. Franses, O. A. Basaran, and C.-H. Chang, CurrentOpinion in Colloid and Interface Science , 296 (1996),ISSN 1359-0294[42] P. J. Steinhardt, D. R. Nelson, and M. Ronchetti, Phys.Rev. B , 784 (1983)[43] W. Lechner and C. Dellago, Journal of Chemical Physics , 114707 (2008)[44] K. Sandomirski, E. Allahyarov, H. L¨owen, and S. U. Egel-haaf, Soft Matter , 8050 (2011) [45] H. E. A. Huitema, M. J. Vlot, and J. P. van der Eerden,Journal of Chemical Physics , 4714 (1999)[46] K. Meier, A. Laesecke, and S. Kabelac, Journal of Chem-ical Physics , 9526 (2004)[47] K. A. Jackson, Interface Science , 159 (2002)[48] H. A. Wilson, Phil. Mag. , 238 (1900)[49] J. Frenkel, Physik Z. der Sowjet Union , 498 (1932)[50] E. A. Mastny and J. J. de Pablo, Journal of ChemicalPhysics , 104504 (2007)[51] A. Kerrache, J. Horbach, and K. Binder, EurophysicsLetters , 58001 (2008) [52] F. Turci, T. Schilling, M. H. Yamani, and M. Oettel, Eu-ropean Physical Journal-Special Topics , 421 (2014)[53] A bond with a neighboring particle is named “solid” ifthe Steinheradt’s order parameter product (cid:126)q · (cid:126)q > . , 021126 (2010)[55] H. Emmerich, H. L¨owen, R. Wittkowski, T. Gruhn, G. I.T´oth, G. Tegze, and L. Gr´an´asy, Advances in Physics ,665 (2012)[56] C. P. Royall, J. Dzubiella, M. Schmidt, and A. vanBlaaderen, Phys. Rev. Lett.98