Mesoscale computational study of the nano-crystallization of amorphous Ge via a self-consistent atomistic - phase field coupling
MMesoscale computational study of the nano-crystallization ofamorphous Ge via a self-consistent atomistic - phase field coupling
C. Reina , , L. Sandoval , and J. Marian ∗ Science and Technology Principal Directorate,Lawrence Livermore National Laboratory, Livermore, CA 94551, USA Department of Mechanical Engineering and Applied Mechanics,University of Pennsylvania, Philadelphia, PA 19104-6315, USA and Theoretical Division T-1, Los Alamos National Laboratory, Los Alamos, NM 87545, USA (Dated: September 14, 2018)Germanium is the base element in many phase-change materials, i.e. systems that can undergoreversible transformations between their crystalline and amorphous phases. They are widely usedin current digital electronics and hold great promise for the next generation of non-volatile mem-ory devices. However, the ultra fast phase transformations required for these applications can beexceedingly complex even for single component systems, and a full physical understanding of thesephenomena is still lacking. In this paper we study nucleation and growth of crystalline Ge fromamorphous thin films at high temperature using phase field models informed by atomistic calcula-tions of fundamental material properties. The atomistic calculations capture the full anisotropy ofthe Ge crystal lattice, which results in orientation dependences for interfacial energies and mobilities.These orientation relations are then exactly recovered by the phase field model at finite thicknessvia a novel parametrization strategy based on invariance solutions of the Allen-Cahn equations. Bymeans of this multiscale approach, we study the interplay between nucleation and growth and findthat the relation between the mean radius of the crystallized Ge grains and the nucleation ratefollows simple Avrami-type scaling laws. We argue that these can be used to cover a wide region ofthe nucleation rate space, hence facilitating comparison with experiments.
PACS numbers: 46.35.+zKeywords: phase transition, Germanium, phase field, multiscale
I. INTRODUCTION
The last two decades have seen an increas-ing interest in phase-changing materials (PCM),such as Ge-Sb-Te (GST), as promising candi-dates for next-generation nonvolatile memorydevices . These materials display ultrafast andreversible phase transitions between their crys-talline and amorphous states at high temper-ature, while displaying remarkable stability atlow temperatures. Their potential originatesfrom the striking differences in electrical and op-tical properties exhibited by each phase, whichcan reach several orders of magnitude in someinstances . These two properties of GST ma-terials —reversible kinetics and high contrastin phase properties— makes them ideal can-didates for memory applications. Generally, aphase transformation is induced by some con-trolled external heat source that drives the sys- tem above the transition threshold.However, a complete understanding of thesefast phase transformations at ultra high heat-ing rates and their remarkable properties is stilllacking. Even for single component systemssuch as Ge, the crystalline structures induced bylaser-heating of amorphous thin films are verycomplex and depend on many factors such asthe deposition history and morphology of theamorphous structure, the intensity and dura-tion of the heating pulse, and the substratestructure, among others . The interplay be-tween all these processes governs the kinetic andthermodynamic evolution of PCM, which mustbe understood in terms of the microscopic prop-erties of each phase. Modeling and simulationcan provide physical insight into the governingphenomena, as well as a path to understandingthe physical evolution of PCM under the condi-tions just discussed. a r X i v : . [ c ond - m a t . m t r l - s c i ] D ec The goal of this paper is two-fold: (i) topresent a novel multiscale model for phasetransformations in a heat bath at constant tem-perature; and (ii) to obtain scaling laws for Ge-crystal mean radius evolution as a function ofthe nucleation rate during the recrystallizationprocess. Our approach consists of a thermody-namically consistent phase field model (TCPF)that reproduces exactly the interface energeticsand kinetics of atomistically computed crystal-lization fronts. As an additional feature, theinterface thickness may be chosen arbitrarilylarge while preserving this exact atomistic-to-continuum coupling, thus delivering a highly ef-ficient multiscale computational model.The paper is organized as follows. We beginin Section II by discussing the physical processescaptured by the model (nucleation, growth, roleof anisotropy, etc.) and the boundary condi-tions employed. In Section III we describe themathematical formulation of the TCPF model,and prove that the energetics and kinetics ofplanar fronts at constant temperature obey uni-versal relations when expressed in appropriatenon-dimensional units. These relations leadto a direct correspondence between the phasefield parameters and physical quantities thatcan be measured via atomistic simulations, andthus provide with an accurate coupling strat-egy between the two scales. Additionally, inthe same section, we present detailed informa-tion on the evaluation of continuum data fromatomistic calculations, namely, the free energiesfor each of the faces, the interface velocities andthe transformation pathway between the amor-phous and crystalline phase. Next, in SectionIV we show the results of the atomistic calcu-lations as a function of temperature and inter-face orientation, compute the invariance solu-tions of the Allen-Cahn equations and obtainthe related phase field parameters. Followingthese results, we use the developed multiscalemodel to analyze, in Section V, the interplaybetween nucleation and growth during crystal-lization kinetics under isotropic and anisotropicgrowth conditions. We conclude the paper inSection VI, with a summary of the results.
II. PHYSICAL MODEL
When an amorphous Ge thin film is heatedabove the glass transition temperature ( T g ) butbelow the melting temperature ( T m ), it sponta-neously crystallizes. The difference in free en-ergies between the two phases constitutes thedriving force for crystallization, which decreaseswith temperature. At the same time, the mobil-ities and nucleation rates rapidly increase as thetemperature rises. For most PCMs, there is anintermediate temperature at which an optimalcompromise is achieved, and the crystallizationrate is maximized. For Ge, such temperatureis around 1100 K, cf. Sec. IV A, which can beachieved, for instance, by high-intensity pulsedlaser heating . Under those conditions, crys-tals nucleate and grow until full impingement isachieved. Complete crystallization occurs overtens of nanoseconds, and the resulting crystalshave a mean radius of tens of nanometers .Experiments suggest that nucleation occurshomogeneously in space throughout the thinfilm with no preferential texture, which typi-cally indicates heterogeneous nucleation at uni-formly distributed imperfections and/or irregu-larities in the amorphous microstructure. Froma theoretical point of view, nucleation is ill-defined despite substantial efforts , and herethe nucleation rate is simply taken taken as aparameter of the model, which can be changedto understand its effect on the overall crystal-lization process.Once crystals have nucleated, they grow byfront propagation in the surrounding amor-phous material, or, if in contact with anothercrystal, by grain boundary motion. In this pa-per, we focus solely on the former, and leavethe study of coarsening of the recrystallizednanograined material for future studies.The mesoscopic time and length scales in-volved in the crystallization process (micronsand tens of nanoseconds) are prohibitive forfully atomistic simulations. Hence, in this pa-per we adopt a computational approach basedon the phase field equations parameterized withatomistic data. The details of the method aredescribed in the following section. In the re-mainder of the paper it will be assumed thatthe film is thin enough for critical nuclei tospan the full thickness and for the problem tobe treated as two-dimensional. Furthermorewe will consider that the amorphous thin filmis free-standing and capable of fully releasingstresses via bending. III. METHODSA. The phase field model
The phase-field method is a computa-tional approach widely used in modelingmorphological and microstructural evolu-tion in materials , including solidification ,fracture , dislocation dynamics , crystalnucleation and growth , and diffusionprocesses , among others. The method uses aso-called phase-field variable, φ , to distinguishamong the different material phases and im-plicitly track interface motion. For example,in this work we take φ = 1 to represent thecrystalline phase, and φ = 0 for the amorphousregion. The interface between the two phasesis then described as a smooth transition of thevariable φ through a small thickness.In phase field models, the equilibrium con-figuration of a material point in homogeneousstate is assumed to be characterized by a Gibbsfree energy g ( φ, T, p, . . . ) that depends on thephase field variable and a set of thermodynamicstate variables such as temperature and pres-sure. The total free energy of the system isthen described via a Landau-Ginzburg free en-ergy functional of the form G = (cid:90) Ω (cid:104) g ( φ, T, p, . . . )+ m ( θ, T, p, . . . )2 ( ∇ φ ) (cid:105) d Ω , (1)where the second term accounts for the ad-ditional contribution due to the presence ofthe interface, which is characterized by a non-vanishing value of ∇ φ . The variable m in Eq. 1sets the interfacial energy and it is a functionof the state variables and the misorientation θ ,defined as the relative angle between the crystal lattice θ C and the normal to the interface, i.e. θ = θ C − arccos( ∇ φ |∇ φ | ). The existence of a totalfree energy of the form just described is a stan-dard postulate in transport processes in contin-uum mechanics . For the system of interest inthis paper, we will assume in the following thatthe only thermodynamic variable dependence iswith the temperature, which is assumed to becontrolled by thermostatting to a heat bath.The system’s time evolution is then modeledusing Allen-Cahn’s equation ∂φ∂t = − M ( T, θ ) δ G δφ = M ( T, θ ) (cid:34) m ( T, θ ) ∇ φ − (cid:18) ∂g∂φ (cid:19) T (cid:35) , (2)where M > . This simpleform of the kinetic relation at constant temper-ature guarantees a decrease of total free energyof the system G [ T, φ ] along the solution path, inaccordance to the second law of thermodynam-ics.We further assume that the local free energyfor intermediate values of φ can be computedfrom the amorphous free energy g A ( T ) and thecrystalline one g C ( T ) as g ( T, φ ) = g A ( T ) − q ( φ )∆ g ( T ) + B ( T ) g dw ( T, φ ) , (3)with q ( φ ) = φ (10 − φ + 6 φ ) ,g dw ( T, φ ) = 16 φ (1 − φ ) , ∆ g ( T ) = g A ( T ) − g C ( T ) . (4)∆ g ( T ) > B ( T ) represents the en-ergy barrier along the transformation path (cid:16) B ( T ) = g ( T, φ = 1 / − g A ( T )+ g C ( T )2 (cid:17) . Thisfunctional form guarantees that the free energyhas two local minima at all temperatures at φ = 0 and φ = 1, which correspond respectivelyto the pure amorphous and crystalline phases.The kinetic equation (2) can then be written as1 M ( T, θ ) ∂φ∂t = m ( T, θ ) ∇ φ − B ( T ) g (cid:48) dw ( φ )+ ∆ g ( T ) q (cid:48) ( φ ) , (5)or equivalently, as α ( T, θ ) (cid:15) ( T, θ ) ∂φ∂t = (cid:15) ( T, θ ) ∇ φ − g (cid:48) dw ( φ )+ ∆¯ g ( T ) q (cid:48) ( φ ) , (6)where (cid:15) ( T, θ ) := m ( T, θ ) B ( T ) , (7) α ( T, θ ) := 1 M ( T, θ ) m ( T, θ ) , (8)∆¯ g ( T ) := ∆ g ( T ) B ( T ) , (9)and (cid:15) has units of length, α has units of aninverse diffusivity, and ∆¯ g is non-dimensional.Equation (6) determines the growth of existingnuclei. Nucleation, for its part, is treated explic-itly in this work by introducing cylindrical crit-ical nuclei in the simulation with random orien-tation and with a probability that is consistentwith a certain nucleation rate. It is assumedthat the nuclei span the entire thickness of thefilm, so that their critical radius R c ( T ) in 2Dcan be estimated via classical nucleation theoryas R c ( T ) = 2 g int ( T, θ )∆ g ( T ) , (10)where g int is the interfacial free energy. B. Invariance relations of the phase fieldequation
1. Steady state velocity
The evolution of a flat crystalline-amorphousinterface at constant temperature obeys an in-variance relation that we proceed to character-ize. To that end, consider a set of parameters (cid:15) and α and the corresponding ones of a referenceproblem ( (cid:15) , α ). The solutions characterizingthe motion of the front for both sets of param-eters are denoted by φ and φ respectively, andare given by Eq. (6) α(cid:15) ∂φ∂t = (cid:15) ∇ φ − g (cid:48) dw ( φ ) + ∇ ¯ g q (cid:48) ( φ ) ,α (cid:15) ∂φ ∂t = (cid:15) ∇ φ − g (cid:48) dw ( φ ) + ∇ ¯ g q (cid:48) ( φ ) . (11)These equations lead to a planar front mov-ing at constant velocity in steady state, that wecall v and v respectively. With the change ofvariables φ ( x, t ) = φ ( x − vt ) = φ ∗ ( x ∗ ) ,φ ( x, t ) = φ ( x − v t ) = φ ∗ ( x ∗ ) , (12)Eqs. (11) transform into the following ordinarydifferential equations − α(cid:15) v ∇ ∗ φ ∗ = (cid:15) ∇ ∗ φ ∗ − g (cid:48) dw ( φ ∗ ) + ∇ ¯ g q (cid:48) ( φ ∗ ) , − α (cid:15) v ∇ ∗ φ ∗ = (cid:15) ∇ ∗ φ ∗ − g (cid:48) dw ( φ ∗ ) + ∇ ¯ g q (cid:48) ( φ ∗ ) , (13)where ∇ ∗ represents the gradient with respectto coordinates x ∗ . If we now define y ∗ = (cid:15) (cid:15) x ∗ (14)and denote φ ∗ ( x ∗ ) = φ ∗ (cid:16) (cid:15)(cid:15) y ∗ (cid:17) = φ ∗∗ ( y ∗ ), thenthe two equations become − α(cid:15)(cid:15) v ∇ ∗∗ φ ∗∗ = (cid:15) ∇ ∗∗ φ ∗∗ − g (cid:48) dw ( φ ∗∗ )+ ∇ ¯ g q (cid:48) ( φ ∗∗ ) − α (cid:15) v ∇ ∗ φ ∗ = (cid:15) ∇ ∗ φ ∗ − g (cid:48) dw ( φ ∗ )+ ∇ ¯ g q (cid:48) ( φ ∗ ) , (15)which are identical for α (cid:15) v = α(cid:15)v . In otherwords, the non-dimensional steady state veloc-ity ( α(cid:15)v ) satisfies an invariance relation of theform α ( T, θ ) (cid:15) ( T, θ ) v ( T, θ ) = I v (∆¯ g ( T, θ )) (16)
2. Interfacial energy
Similarly, the non-dimensional interface en-ergy (cid:16) g int ( T,θ ) B ( T ) (cid:15) ( T,θ ) (cid:17) can be shown to satisfy g int ( T, θ ) B ( T, θ ) (cid:15) ( T, θ ) = I g (∆¯ g ( T, θ )) . (17)This relation can be obtained by defining theinterface energy as the difference between thesum of energies of the bulk systems and the en-ergy of the two-phase configuration separatedby a sharp interface. Choosing a reference framethat moves with the interface steady state ve-locity, the surface energy reads g int ( T, θ ) = (cid:90) (cid:96) (cid:104) m ( T, θ )2 ( ∇ ∗ φ ∗ ) + B ( T ) g dw ( φ ∗ ) − ∆ g ( T ) q ( φ ∗ ) + g A ( T ) (cid:105) dx ∗ − h (cid:96) g A ( T ) + g C ( T )) , (18)where (cid:96) is the range of the coordinate orthogo-nal to the interface and h (cid:96) = | (cid:96) | represents thesystem length in such direction. Without loss ofgenerality, we consider that the interface is atthe center of the domain under consideration.Then, at constant temperature, g int ( T, θ ) = m ( T, θ ) (cid:34) (cid:90) (cid:96) ( ∇ ∗ φ ∗ ) dx ∗ (cid:35) + B ( T ) (cid:34) (cid:90) (cid:96) g dw ( φ ∗ ) dx ∗ (cid:35) + ∆ g ( T ) (cid:34) h (cid:96) − (cid:90) (cid:96) q ( φ ∗ ) dx ∗ (cid:35) . (19)By recourse to the changes of variableexpressed by Eq. (14) and the equivalence φ ∗∗ ( y ∗ ) = φ ∗ ( x ∗ ), this equation can be rewrit- ten as g int ( T, θ ) B ( T ) (cid:15) ( T, θ ) = (cid:15) (cid:34) (cid:90) (cid:96) ( ∇ ∗ φ ∗ ) dx ∗ (cid:35) + 1 (cid:15) (cid:34) (cid:90) (cid:96) g dw ( φ ∗ ) dx ∗ (cid:35) + ∆¯ g ( T, θ ) (cid:15) (cid:34) h (cid:96) − (cid:90) (cid:96) q ( φ ∗ ) dx ∗ (cid:35) = I g (∆¯ g ( T, θ )) (20)and the sought-after result is obtained.
C. Coupling strategy
The phase field model at constant tempera-ture described by Eq. (6) requires that the fol-lowing magnitudes be defined: g ( T, φ ), M ( T, θ ), m ( T, θ ) and B ( T ). In light of the invariance re-lations previously derived, these parameters canbe obtained directly via atomistic simulations ofthe following kind:(a) Free energy calculations using thermo-dynamic integration of bulk crystallineand amorphous phases at zero pressureand constant temperatures. These yield g C ( T ) and g A ( T ).(b) Transition energy barrier between theamorphous and crystalline states at zeropressure and constant temperature. Toobtain B ( T ) we equate the transition pathto an activated process described by ageneral configurational nonlinear many-body reaction coordinate .(c) Free energy calculations at zero pres-sure and constant temperature to com-pute g int ( T, θ ). As well, molecular dy-namics simulations of moving fronts insteady state to extract interface velocities v ( T, θ ).Equations (16) and (17) allow to explicitlyrelate the atomistic data with the phase fieldparameters M ( T, θ ) and m ( T, θ ) as m ( T, θ ) = g int ( T, θ ) (cid:112) B ( T ) I g (∆¯ g ) , (21) τ ( T, θ ) = 1 M ( T, θ ) = g int ( T, θ ) v ( T, θ ) I v (∆¯ g ) I g (∆¯ g ) . (22)The details of such computations for the differ-ent subsystems are described in the followingsection. D. Free energy calculations from atomisticsimulations
To compute the free energies of bulk phaseswe use thermodynamic integration. TheHelmholtz free energy is defined as F = E − T S, where E and S are the internal energy and theentropy. The approach employed here to com-pute F for solid systems involves two steps.First, a free energy F is obtained at a tem-perature T in each case. Second, the internalenergy is computed for the desired range of tem-peratures as E = (cid:104) V + K (cid:105) , where V and K are the potential and kineticenergies, and, under the usual assumption ofergodicity, (cid:104)·(cid:105) denotes time average. The freeenergy at a given temperature T can then beobtained by solving the Gibbs-Helmholtz inte-gral FT = F T − (cid:90) T E ( τ ) τ dτ, (23)where E ( T ) is made to be a smooth integrablefunction that expresses the temperature depen-dence of the internal energy.To obtain F , we employ the λ -integrationtechnique. The use of this technique is well doc-umented in the literature and here we onlyprovide a brief overview. It is assumed that the system at T is represented by the following en-ergy function˜ U ( λ ) = (1 − λ ) k U (cid:48) + λ k U, (24)where k ≥ U ( r ) and U (cid:48) ( r )are, respectively, the potential energy functionsof the system under study and a reference sys-tem whose free energy is known. For practicalreasons, U (cid:48) is typically taken to be the potentialenergy function of an analytical or a numericalsystem of coupled harmonic oscillators. Fromhere, F is straightforwardly obtained by inte-grating along a λ trajectory between 0 and 1 F ( U ) = F ( U (cid:48) ) ++ (cid:90) k (cid:2) λ k − (cid:104) U (cid:105) − (1 − λ ) k − (cid:104) U (cid:48) (cid:105) (cid:3) dλ. (25)In this work, the calculation hypotheses war-ranted the use of linear thermodynamic inte-gration ( k = 1). The numerical behavior of Eq.(25) under other choices of k has been reviewedin the literarure , and we refer the reader tothose works for more details.When internal transport processes derivedfrom the existence of viscous flow, e.g. as in liq-uids, are important, Eq. (25) can no longer beused in conjunction with a reference harmoniccrystal to compute F . In such cases, the freeenergy is obtained as F = F g ( ρ , T )++ (cid:90) ρ dρ (cid:20) p ( T , ρ ) − ρk B T ρ (cid:21) , (26)where p is the pressure, ρ is the system’s den-sity at a temperature T above the supercriticaltemperature (here we have taken T = 1500 K)and F g ( ρ, T ) is the free energy of an ideal gas atdensity and temperature ρ and T . The inte-grand in the second term of the r.h.s. of Eq. (26)represents an isothermal expansion from ρ tozero density ( i.e. infinite volume), where thesystem effectively behaves as an ideal gas. Thisexpansion should be reversible, which meansthat no first-order transition —e.g. liquid-gas—should be traversed. For its part, the equationof state p ( T , ρ ) is obtained from a set of canon-ical ensemble calculations of liquid Ge at T .This is shown in Fig. 1. The resulting dataare fitted to a 3 rd degree polynomial and theintegral in Eq. (26) is solved to yield the freeenergies. −3.5−3.0−2.5−2.0−1.5−1.0−0.50.00.54.35 4.40 4.45 4.50 4.55 4.60 4.65 4.70 4.75 4.80 P r e ss u r e [ G P a ] r a [ · m −3 ]3 rd −degree polynomial fit FIG. 1. Equation of state for liquid Ge at 1500 K. D [ m / s ] T [K] FIG. 2. Self-diffusivity of amorphous Ge as a func-tion of temperature.
To establish the limits of applicability ofEqs. (25) and (26), we compute the self-diffusivity of amorphous Ge as a function oftemperature. Results are shown in Figure 2.Our calculations clearly demonstrate that in-ternal transport is important only above 900 K,at which point Eq. (26) is the pertinent one forcomputing free energies.
E. Calculations of interfacial free energies
For ideal (perfectly-flat, stress-free) inter-faces, the interfacial free energy γ can be definedas the change in free energy due to an increasein interface surface area Aγ = (cid:18) ∂F∂A (cid:19) N,p,T . From a practical standpoint, interfacial free en-ergies are computed using the Dupr´e equation γ ij = γ i + γ j − W ij , (27)where γ i is the surface free energy of solid i ina vacuum and W ij > adhesion . The above formula results froma two-step construction in which (i) slabs withexposed area A are created in a vacuum frombulk systems i and j , and (ii) these new surfacesare brought into contact, releasing an amountof work W ij . Using the bulk free energy densi-ties, f i , computed in Section III D by thermo-dynamic integration, we can express Eq. (27)as γ ij = F ij (Ω ij ) − ( f i Ω i + f j Ω j )2 A , where F ij is the free energy of a system withtotal volume Ω ij = Ω i + Ω j of a simulationcell containing a volume Ω i of phase i and avolume Ω j of phase j separated by the inter-face. As for the bulk phases, F ij is computedusing λ -integration. The calculation is simpli-fied due to the fact that both phases correspondto monoatomic solid systems and thus the samereference thermodynamic system can be usedfor both. Because cohesion in crystalline solidsis anisotropic –while in amorphous solids or liq-uids is isotropic–, we have repeated this proce-dure for a number of different surface orienta-tions to assess the orientation dependence of theinterfacial free energy.As a final point, we note that, in real-ity, atomistic interfaces are intrinsically roughand accurate methods for computing interfacialfree energies that account for such atomic-scaleroughness have been developed . F. Molecular dynamics calculations ofinterface mobility
To calculate interface mobilities, large-scaleMD simulation of crystal-amorphous mixing areperformed. For this, a bicrystal with the desiredinterface orientation is constructed and the sys-tem is left to evolve at a given temperature.The temperature must be kept constant on theboundaries of this box using some suitable the-mostat. It is important to ensure that no tem-perature control is applied in the interface re-gion, as this could result in thermodynamic ar-tifacts. Following Hoyt and Song we measurethe rate of change in potential energy as thefront advances, dU/dt , and obtain the velocityas v = − A Ω h l dUdt , where A is again the interface surface area, Ωis the volume per atom in the crystalline phaseat the temperature of interest, and h l is the la-tent heat, which in this case is simply equal tothe internal energy difference between the twophases. As above, the factor of two in the aboveexpression reflects the use of periodic boundaryconditions, which introduces two interfaces inthe simulation volume. G. Interatomic potential and boundaryconditions
Next, we discuss the most salient features ofthe potential energy function U ( r ) employedhere. We use a Stillinger-Weber potential pa-rameterized by Posselt and Gabriel for Ge .The potential builds on the parameterization byDing and Anderson , which yields the correctvalues for the cohesive energy and lattice con-stant. The elastic constants, melting point, aswell as the formation and migration energies ofpoint defects for diamond-structure Ge are re-produced within certain limits. The energeticsof other crystalline phases and the structure ofthe liquid are described reasonably well. All the simulations were run in the isobaric-isothermal ensamble N pT –where N is the num-ber of particles, p is the pressure,and T theabolute temperature– using periodic boundaryconditions in three dimensions. The pressureand temperature were controlled using a Nos´e-Hoover thermostat with a damping constant of0.1 ps in both cases. IV. ATOMISTIC RESULTS ANDPARAMETRIZATION OF THE PHASEFIELD MODELA. Molecular dynamic results
First, we calculate the free energy of the crys-talline and amorphous phases as a function oftemperature. Figure 3(a) shows the results foreach case, with the crystalline phase being themore stable one up to a temperature of 1350K. This value can be considered as the melt-ing point for the numerical force field U ( r ) em-ployed here, about 10% larger than the experi-mental value of 1210 K . Figure 3(b) shows thevariation of the atomic density ρ a and the ther-mal expansion coefficient α t with of tempera-ture. As the figure shows, the amorphous struc-ture is always approximately 10% denser thanthe crystalline phase, and is seen to undergoa smooth transition at ≈
800 K. This tempera-ture correlates well with the point at which theself diffusivity is markedly nonzero (cf. Fig. 2),and can be assigned to the reversible liquid-glasstransition. The temperature at which we haverecorded this transition is in excellent agree-ment with experimental measurements of theglass transition temperature for Ge, T g ≈ . This glass-liquid transition manifests itselfas a marked dip in the value of α t . The glasstransition is not itself a phase transition, buta laboratory phenomenon, and depends on thecooling rate and the boundary conditions. Ourcalculations were obtained for systems contain-ing 13824 atoms at a cooling rate of 10 K perps from the liquid melt with 100 ps annealingsevery 100 K.Next, we construct a bilayer system consist- −4.8−4.6−4.4−4.2−4.0−3.8−3.6 0 200 400 600 800 1000 1200 1400 1600 F [ e V pe r a t o m ] T [K] crystallineamorphous (a) r a [ · m − ] a [ · − K − ] Temperature [K]crystallineamorphous (b)
FIG. 3. Results of single phase atomistic simu-lations at zero pressure: (a) free energy and (b)atomic density as a function of temperature. ing of crystalline and amorphous half crystals tomeasure interface properties. We calculate in-terfacial free energies and interface velocities asa function of temperature. To ascertain the ori-entation dependence of the crystal-amorphousinterface, three different crystallographic direc-tions have been considered, namely [100], [110]and the [111]. Figure 4(a) reveals the temper-ature evolution of the interfacial free energy,which is self-consistent in that it vanishes atthe melting point. The [100] orientation is the most stable one at low temperatures, followedby the [111] and the [110] directions. At theglass transition temperature, however, all in-terfacial free energies converge, likely due to a‘wetting’ effect of the crystalline surface in con-tact with a liquid-like system. For their part,the interface velocities also display a behaviorqualitatively similar to the self-diffusion coeffi-cient for intermediate temperatures (cf. Fig. 2).This amounts to practically zero mobility below T g and a monotonic increase above it. For verylarge temperatures though, the driving force de-creases till vanishing at the melting tempera-ture, and induces a decrease in the velocity. Thecompromise between driving force and mobilityleads to a peak in the interface velocity around1100 K.The last parameter to be computed fromatomistic simulations is the energy barrier alongthe crystal-amorphous transformation path, B ( T ), introduced in Section III A. This calcu-lation involves defining a suitable reaction co-ordinate between the initial –amorphous– andfinal –crystalline– states. A natural choice ofinitial and final states is to start from a crys-talline specimen and melt it to the liquid state,followed by a sensible cooling into the amor-phous structure. However, in most constrainedminimization techniques used for this purpose–e.g. nudged elastic band or drag methods–,the reaction coordinate vector is constructedas a linear interpolation of all degrees of free-dom within each end state. Due to the sizableamount of mixing and thermal transport oc-curring during the melting and cooling phases,however, we have seen that this constructionleads to large variations in the reaction coordi-nate that make convergence difficult (and thephysical meaning of the procedure question-able). Instead, we proceed to use the bilayergeometry described above and let the inter-face propagate spontaneously turning the amor-phous half system into a crystalline one. Oncethe conversion is complete, we extract a volu-metric sample from within the converted regionand use it as one of the end states, while us-ing the corresponding atoms contained withinthat volume in the amorphous region as the sec-0 I n t e r f a c e f r ee ene r g y [ J / m ] I n t e r f a c e f r ee ene r g y [ e V / A ] T [K] [100][110][111] (a) I n t e r f a c e v e l o c i t y [ m / s ] T [K][100] [110] [111] (b) FIG. 4. Calculations of a planar interface at zeropressure: (a) surface free energy and (b) front ve-locity as a function of temperature for three surfacenormals. ond end state. Replicas for constrained mini-mization are generated by linear interpolationbetween both end states. This procedure es-tablishes a direct link between both structureswhile minimizing the amount of internal trans-port, providing a more physical transition path-way.Figure IV A shows the relaxed minimum en-ergy path at 0 K for a system containing 13824atoms and 16 replicas. The static energy bar-rier is approximately B (0) ≈ .
49 eV peratom. Finding the temperature dependence in-trinsic in B ( T ) is not trivial and requires us-ing elaborate thermodynamic integration tech-niques ( e.g. refs. ). Here, for lack of a better estimate, we assume that B ( T ) = B (0). −0.2−0.10.00.10.20.30.40.50.0 0.2 0.4 0.6 0.8 1.0 E ne r g y [ e V pe r a t o m ] f FIG. 5. Minimum energy path at zero tempera-ture and pressure for conversion of an amorphoussystem into an ordered crystal. A barrier of 0.49eV per atom is seen to separate both phases, whilethe energy difference between the stable amorphousand crystalline is 0.14 eV per atom, consistent withFig. 3(a) at 0 K. The roughness in the energy pathis a consequence of the multidimensional nature ofthe transition, with many atoms undergoing trans-formations at once.
B. Numerical analyses of the phase fieldmodel
The evolution equation for the phase fieldvariable, c.f. Eq. (6), is numerically solved inthis work by means of a semi-implicit, backwardEuler discretization scheme in time α(cid:15) φ j +1 − φ j ∆ t = (cid:15) ∇ φ j +1 − g (cid:48) dw ( φ j )+∆ g q (cid:48) ( φ j ) , (28)where ∆ t is the time step. The resulting or-dinary differential equation is then solved inspace via a first order finite difference schemefollowed by a discrete Fourier transform ( fftw package ). The grid size is denoted by h .In this section we illustrate the convergence ofthe numerical scheme to known analytical solu-tions and demonstrate good agreement with the1theoretical spatial and temporal convergencerates.
1. Static one-dimensional solution at constanttemperature
We first consider a planar front separatinga two phase region at constant temperaturewith an energy landscape described by a perfectdouble-well potential, i.e. with ∆ g = 0. Theequilibrium solution of this system is character-ized by Eq. (6), that is, α(cid:15) ∂φ∂t = (cid:15) ∇ φ − φ (1 − φ )(1 − φ ) . (29)This equation is analytically soluble and givesa static interface with the following profile φ ( x ) = 12 −
12 tanh (cid:16) x(cid:15) √ (cid:17) , (30)where the x axis is centered at the interface andorthogonal to it. For this solution, (cid:15) exactly cor-responds to the interface thickness, measured asthe support of 0 . < φ < . analytic − Energy
PhaseField
Energy analytic analytic = (cid:90) x : φ =0 . x : φ =0 . φ (1 − φ ) dx. (31)We observe a quadratic convergence of theenergy as it is expected from the first order finitedifference scheme used in these simulations.
2. Steady state one-dimensional solution atconstant temperature
A difference in free energy between the twophases (∆ g >
0) will induce a planar interfaceto move until a constant velocity is reached, as E rr o r ene r g y [ % ] e / h Phase field e =0.2Phase field e =0.3Error = 18.17 ( e / h ) −1.94 FIG. 6. Energy error as a function of the normalizedmesh size. Parameters of the simulation: (cid:15) = 0 . . α = 1, ∆ t = 0 . discussed in Sec. III. Figure 7 shows the value ofthe numerical steady state velocity as a functionof the time step for a constant (small) mesh size.As expected from a first order time scheme,the average velocity converges linearly with thetime step. However, if the mesh size used forthe spacial discretization is too large, the nu-merical solution does not lead to a steady statevelocity. Rather, the underlying mesh acts asa periodic potential inducing oscillations in thesteady state numerical solution, c.f. Fig. 8(a).Additionally, the mean velocity decreases withthe mesh size for a given constant time step,and can even vanish for very large mesh sizes,c.f. Fig. 8(b). This behavior, previously re-ported by other authors , is characteristic ofthe scheme used for solving the ordinary differ-ential equation at a given constant time step.Equivalent simulations using finite elements in-dicate that the average steady state velocitymay increase with the mesh size for sufficientlylarge driving forces. C. Parameters of the phase field model
The functions I v (∆¯ g ) and I g (∆¯ g ) describedin Eqs. (16) and (17), can now be computedwith a controlled numerical error lower than 1%.In Figures 9(a) and 9(b) we plot each invariant2 A v e r age s t ead y s t a t e f r on t v e l o c i t y D t Phase field v =−752 D t +5.24 FIG. 7. Average steady state velocity with respectto time step. (cid:15) = 0 . α = 1, ∆¯ g = 1, h = 0 . for different values of the phase field parame-ters. For small values of ∆¯ g , I v (∆¯ g ) and I g (∆¯ g )can be approximated by α(cid:15)v (cid:39) ∆ gB , (32) g int (cid:39) (cid:15)B. (33)The linear dependence given in Eq. (32) be-tween the interface velocity and the free energydifference is often assumed in the literature .For its part, Eq. (33) is the result of the equipar-tition of the interface energy at zero velocity, asshown in Appendix A.Equations (21) and (22) provide a directrecipe for computing m ( T, θ ) and M ( T, θ ) fromthe functions I v (∆¯ g ) and I g (∆¯ g ). These areshown in Figs. 10 and 11 for different valuesof B . Equivalently, the values of (cid:15) ( T, θ ) and α ( T, θ ), defined in Eqs. (7) and (8), are shownin Figs. 12 and 13 respectively. As the figuresand equations show, the phase field parameterscan be adjusted to deliver the appropriate in-terface energy and front velocity at constanttemperature for an arbitrary value of the freeenergy barrier B ( T ). In view of this obser-vation, and the computational cost associatedwith the computation of B ( T ), the value ob-tained in Sec. IV A is used. Other strategies forselecting B based solely on computational cost,i.e. fixing the interface thickness, may also be A v e r age s t ead y s t a t e f r on t v e l o c i t y h (a)Average velocity of the steady state numericalsolution with respect to the mesh size for a constanttime step. F r on t v e l o c i t y time h =1/5 h =1/10 h =1/15 h =1/20 h =1/25 (b)Evolution of the front velocity with time forseveral mesh sizes. FIG. 8. Numerical behavior of a planar amorphous-crystalline front with increasing mesh size. Parame-ters of the simulation: (cid:15) = 0 . α = 1, ∆ t = 0 . g = 2. used. In such case, the velocity and free energyof the interface would be exactly captured bythe phase field model with a finite and poten-tially thick interface. V. SIMULATIONS OFGE-NANOCRYSTALLIZATION
Next, the multiscale model developed in theprevious sections is used to examine the inter-play between the nucleation and growth mech-3 a e v D - g e =0.2, a =1 e =0.5, a =1 e =0.2, a =2 (a) −1.0−0.50.00.51.00.0 1.0 2.0 3.0 4.0 5.0 g i n t B e D - g e =0.2, a =1 e =0.5, a =1 e =0.2, a =2 (b) FIG. 9. Invariance relations for the steady statesolution of a planar crystalline-amorphous inter-face: (a) non-dimensional interface velocity and (b)free energy versus non-dimensional driving force.Parameters of the simulations: ∆ t = 5 × − , h = 0 . anisms in the crystallization process. For verifi-cation purposes, first we examine the ideal caseof isotropic growth and demonstrate its consis-tency with the Avrami relations. To that end,we consider a temperature of 1100 K, nucleationrates that span several orders of magnitude, andgrowth conditions dictated by the phase fieldparameters previously computed for the threesurface orientations. In these isotropic calcula-tions, the parameter θ is simply a differentiatingproperty for each growth condition.We begin by rendering Eq. (6) dimension-less with the following characteristic length and m [ Å ‘‘(cid:214) e V ] T [K] B=0.1 eVB=0.2 eVB=0.3 eVB=0.4 eVB=0.49 eV (a) m [ Å ‘‘(cid:214) e V ] T [K] B=0.1 eVB=0.2 eVB=0.3 eVB=0.4 eVB=0.49 eV (b) m [ Å ‘‘(cid:214) e V ] T [K] B=0.1 eVB=0.2 eVB=0.3 eVB=0.4 eVB=0.49 eV (c)
FIG. 10. Phase field parameter m as a functionof temperature for several values of the free energybarrier B : (a) (111) surface, (b) (110) surface and(c) (100) surface. timescale l c = R c ,t c = αR c , (34)where R c is the radius of the critical nuclei,given by Eq. (10). This specific choice of lengthand time scales has the property of deliveringnon-dimensional parameters (cid:15) and α that are4 M [ p s − e V − ] T [K] B=0.1 eVB=0.2 eVB=0.3 eVB=0.4 eVB=0.49 eV (a) M [ p s − e V − ] T [K] B=0.1 eVB=0.2 eVB=0.3 eVB=0.4 eVB=0.49 eV (b) M [ p s − e V − ] T [K] B=0.1 eVB=0.2 eVB=0.3 eVB=0.4 eVB=0.49 eV (c)
FIG. 11. Interface mobility as a function of temper-ature for several values of the free energy barrier B :(a) (111) surface, (b) (110) surface and (c) (100)surface. independent of θ . In particular, using Eq. (21),¯ (cid:15) := (cid:15)R c = m ( T, θ )∆ g ( T ) √ B g int ( T, θ ) = ∆¯ g ( T )2 I g (∆¯ g ( T )) , ¯ α := αR c t c = 1 . (35)Consequently, the resulting non-dimensional e [ Å ] T [K] B=0.1 eVB=0.2 eVB=0.3 eVB=0.4 eVB=0.49 eV (a) e [ Å ] T [K] B=0.1 eVB=0.2 eVB=0.3 eVB=0.4 eVB=0.49 eV (b) e [ Å ] T [K] B=0.1 eVB=0.2 eVB=0.3 eVB=0.4 eVB=0.49 eV (c)
FIG. 12. Phase field parameter (cid:15) as a function oftemperature for several values of the free energybarrier B : (a) (111) surface, (b) (110) surface and(c) (100) surface. phase field equation¯ (cid:15) ∂ ¯ φ∂ ¯ t = ¯ (cid:15) ¯ ∇ ¯ φ − g (cid:48) dw ( ¯ φ ) − ∆¯ g q (cid:48) ( ¯ φ ) (36)allows to examine all growth conditions witha single equation, which we proceed to couplewith an explicit nucleation strategy, with non-dimensional nucleation rates ranging between10 − and 10 − . The numerical results atten-5 a [ p s Å − ] T [K] B=0.1 eVB=0.2 eVB=0.3 eVB=0.4 eVB=0.49 eV (a) a [ p s Å − ] T [K] B=0.1 eVB=0.2 eVB=0.3 eVB=0.4 eVB=0.49 eV (b) a [ p s Å − ] T [K] B=0.1 eVB=0.2 eVB=0.3 eVB=0.4 eVB=0.49 eV (c)
FIG. 13. Phase field parameter α as a function oftemperature for several values of the free energybarrier B : (a) (111) surface, (b) (110) surface and(c) (100) surface dant to this coupled nucleation and growth pro-cess are shown in Fig. 14(a), where the time evo-lution of the mean radius during the nanocrys-tallization process is represented. These curvesmay be normalized with the time and mean ra-dius at 95% crystallization, delivering a univer-sal crystallizatian curve, as shown in Fig. 14(b).For their part, the time required to complete95% of the crystallization, and the mean ra- dius at that stage of crystallization, satisfy verysimple relations with the nucleation rate as de-picted in Figs. 15(a) and 15(b). The slope of thefit to the numerical points in both curves coin-cide with the exponents of Avrami’s relations,as detailed in Appendix B. m ean r ad i u s time J =10 −5 J =5 · −5 J =10 −4 J =5 · −4 J =10 −3 (a) no r m a li z ed m ean r ad i u s time / time for 95% crystallization J =10 −5 J =5 · −5 J =10 −4 J =5 · −4 J =10 −3 (b) FIG. 14. Mean radius evolution with time for sev-eral nucleation rates. Isotropic growth conditionsand spatially homogenous nucleation. (a) Datanon-dimensionalized using l c = R c and t c = αR c aslength and time scale respectively, (b) Mean radiusand time normalized with the corresponding valueat 95% crystallization. Simulation box: 400 × (cid:15) = 0 . α = 1, h = 0 .
04, ∆ t = 0 . Next, we proceed to examine the case ofanisotropic growth at T = 1100 K using the sur-face orientation dependent data for the phasefield parameters. In particular, we use a linearinterpolation scheme between the [111], [100]and [110] data previously computed, to define α ( θ ) and (cid:15) ( θ ) for an arbitrary value of the mis-orientation angle θ ∈ [0 , π ), c.f. Fig. 16. Thespecific spacing between the three data points6 −7 −6 −5 −4 −3 −2 −1 t i m e f o r % c r ys t a lli z a t i on nucleation rate t = 1.06 J −0.33 Phase field simulations (a) t i m e f o r % c r ys t a lli z a t i on mean radius 95% crystallization R = 0.58 t + 1.27Phase field simulations (b) FIG. 15. Scaling laws between the nucleationrate and the time and mean radius to achieve95% crystallization. Spatially homogeneous nucle-ation and isotropic growth conditions. Data non-dimensionalized with l c = R c and t c = αR c aslength and time scale respectively. has been chosen proportional to the number ofdistinct (cid:104) (cid:105) , (cid:104) (cid:105) and (cid:104) (cid:105) orientations in adiamond cell structure (4, 3 and 6 respectively).The resulting curves for the mean radius evo-lution as a function of time for different non-dimensional nucleation rates follow in a verysimilar manner to the isotropic case, c.f. Fig. 17.Despite the strong growth anisotropy exhibitedby Ge, c.f. Sec. IV A, the time for crystalliza-tion and the mean radius of the final microstruc-ture also follow fairly closely Avrami’s relations,c.f. Fig. 18(a) and 18(b). The physical rea-soning behind this result lies on the fact thatthe final microstructure does not exhibit grainswith a large aspect ratio, c.f. Fig. 19.. The ho-mogenous nucleation creates constraints to thegrowth of preexisting nuclei that are, on aver-age, equi-spaced, leading to anisotropic ratios ✓ ✓ = 3 .
513 2 ⇥✓ = 813 2 ⇡ [100][110] [111] FIG. 16. The values of (cid:15) and α for a given misorien-tation θ ∈ [0 , π ) will be obtained by linear interpo-lation of the values for orientations [111], [100] and[110] spaced as shown in the figure. no r m a li z ed m ean r ad i u s time / time for 95% crystallization J =9.00 · s −1 m −2 J =45.02 · s −1 m −2 J =90.03 · s −1 m −2 J =450.16 · s −1 m −2 J =900.32 · s −1 m −2 FIG. 17. Normalized mean radius evolution ver-sus normalized time for several nucleation rates.Spatially homogeneous nucleation and anisotropicgrowth conditions. of the final microstructure similar to the sim-ulations of isotropic growth, c.f. Fig. 21. Themeasure of anisotropy used in these calculationsis AR = (cid:114) λ λ , (37)where λ > λ are the eigenvalues of the gyra-tion tensor of each grain, see Fig. 20 for severalexample of AR.7 t i m e f o r % c r ys t a lli z a t i on [ n s ] nucleation rate [s −1 m −2 ] t = 20.7 · J −0.34 Phase field simulations (a) m ean r ad i u s % c r ys t a lli z a t i on [ n m ] time for 95% crystallization [ns] R = 0.16 R + 0.79Phase field simulations (b) FIG. 18. Scaling laws between the nucleation rateand the time and mean radius to achieve 95% crys-tallization. Spatially homogeneous nucleation andanisotropic growth conditions.
VI. SUMMARY AND CONCLUDINGREMARKS
The results shown in this paper indicate thatthe Ge nanocrystalline structure resulting fromhighly anisotropic growth and homogenous nu-cleation is fairly isotropic in the range of nu-cleation rates examined. As a result, the rela-tions between the crystallization time, final av-erage mean radius and nucleation rate followfairly closely Avrami’s relations. These rela-tions could be potentially used to extract nu-cleation rates from experimental measurementsof time-to-crystallization and grain size or viceversa. However, it must be kept in mind thatthe nucleation rates used here for numericalconvenience may differ from those found in lab-oratory experiments depending on the experi-mental conditions , and any exercise of ex-trapolation for the anisotropic curves is not backed by our model.The numerical strategy chosen to obtainthese results consists on phase field simulationsthat exactly deliver the interface kinetics andenergetics obtained via atomistic calculations.This exact mapping between the atomistic andcontinuum solutions is based on two invariancerelations of the Allen-Cahn equations at con-stant temperature. More precisely, we showin this paper that the non-dimensional steadystate velocity and interfacial free energy of pla-nar crystallization fronts satisfy universal rela-tions with the non-dimensional free energy dif-ference between the crystalline and amorphousphase. These two relations allow for an ex-plicit parametrization of the phase field equa-tions with the atomistic data calculated frommolecular dynamic simulations in suitably con-structed subsystems using state-of-the-art inter-atomic potentials.
Appendix A: Free energy of planar stillinterface
Consider a flat crystalline-amorphous inter-face with an energy landscape described by aperfect double-well potential, i.e. with ∆ g = 0.The equilibrium solution consists of a static in-terface that satisfies0 = m ∇ φ − B g (cid:48) dw ( φ ) . (A1)For this specific case, the interface profile isalso the minimizer of the interface free energyfunctional, c.f. Eq. (18), g int = m (cid:34) (cid:90) (cid:96) ( ∇ φ ) dx (cid:35) + B (cid:34) (cid:90) (cid:96) g dw ( φ ) dx (cid:35) , (A2)where x to be the coordinate orthogonal to theinterface. Then, the integrand of this functionalcan be seen as a Lagrangian, that exclusivelydepends on φ and dφdx and does not depend ex-plicitly on the variable x , L (cid:18) φ, dφdx (cid:19) = m ∇ φ ) + B g dw ( φ ) . (A3)8As a result, by Noether’s theorems, the cor-responding Hamiltonian (Legendre transform ofthe Lagrangian) is independent of the variable x . More precisely, L − ∂L∂ ( dφ/dx ) dφdx = constant , (A4)or equivalently, − m ∇ φ ) + B g dw ( φ ) = constant . (A5)The constant may be evaluated far away fromthe interface, where φ is either 0 or 1, deliveringa zero value. As a result, the interface energyis partitioned equally between the two terms inEq. (A2), g int m (cid:34) (cid:90) (cid:96) ( ∇ φ ) dx (cid:35) = B (cid:34) (cid:90) (cid:96) g dw ( φ ) dx (cid:35) . (A6)The value of either of these two integrals maybe evaluated for the equilibrium solution, c.f.Eq. (30), giving g int = √ (cid:15)B (cid:39) (cid:15)B. (A7)We further remark that the value of the in-tegrals in Eq. (A6) are not very sensitive tothe specific interface profile. A simple piece-wise linear profile with interface thickness (cid:15) , i.e. φ ( x ) = x(cid:15) for 0 < x < (cid:15) , would deliver g int = 1615 (cid:15)B (cid:39) (cid:15)B. (A8) Appendix B: Crystallization time and meanradius at crystallization as a function of thenucleation rate
Avrami’s equation describes the kinetics ofcrystallization under the assumption of homoge-neous nucleation, constant nucleation rate andconstant growth rate for the radius of each nu- clei independently of its size. For a two dimen-sional system, it readslog ( − log(1 − f C )) = log (cid:32) π ˙ R J (cid:33) + 3 log t, (B1)where f C is the volume fraction of the crys-talline phase and ˙ R is the growth rate. Thisequation immediatly delivers the relation be-tween the time required to acheive 95% crys-tallization ( f C = 0 .
95) and the nucleation rate.Namely, it is of the form t . = CJ − / , (B2)where C = (cid:16) π ˙ R (cid:17) − / ( − log 0 . / .On its side, the mean radius at 95% crystal-lization ( ¯ R . ) satisfies0 . VN . = π ¯ R . , (B3)where V is the total volume and N . is thenumber of nuclei composing the nanocrystallinestructure. The later, can be computed as N . = (cid:90) t . ˙ N dt = (cid:90) t . JV (1 − f C ( t )) dt = JV (cid:90) CJ − / exp (cid:32) − π ˙ R t J (cid:33) dt. (B4)By use of the change of variables τ = tJ / ,Eq. B4 simplifies to N . = J / V (cid:90) C exp (cid:32) − π ˙ R τ (cid:33) dτ = C J / , (B5)with C solely dependent on ˙ R . Then, the meanradius at crystallization can be expressed as¯ R . = (cid:18) . VπC (cid:19) / J − / (B6)or ¯ R . = (cid:18) . VπC C (cid:19) / t . . (B7)9 ACKNOWLEDGEMENTS
This work performed under the auspices ofthe U.S. Department of Energy by LawrenceLivermore National Laboratory under Contract DE AC52-07NA27344. C. Reina acknowledgessupport from the Lawrence Fellowship programat Lawrence Livermore National Laboratory.J. Marian acknowledges support from DOE’sEarly Career Research Program. ∗ [email protected] Matthias Wuttig. Phase-change materials: To-wards a universal memory?
Nature Materials ,4:265 – 266, 2005. Matthias Wuttig and Noboru Yamada. Phase-change materials for rewriteable data storage.
Nature Materials , 6:824 – 832, 2007. Liliya Nikolova, Thomas LaGrange, Mark J.Stern, Jennifer M. MacLeod, Bryan W. Reed,Heide Ibrahim, Geoffrey H. Campbell, FedericoRosei, and Bradley J. Siwick. Complex crystal-lization dynamics in amorphous germanium ob-served with dynamic transmission electron mi-croscopy.
Phys. Rev. B , 87:064105, Feb 2013. Chubing Peng, Lu Cheng, and M. Mansuripur.Experimental and theoretical investigations oflaser-induced crystallization and amorphizationin phase-change optical recording media.
Jour-nal of Applied Physics , 82(9), 1997. Zhimei Sun, Jian Zhou, and Rajeev Ahuja.Unique melting behavior in phase-change ma-terials for rewritable data storage.
Phys. Rev.Lett. , 98:055505, Feb 2007. L. Nikolova, T. LaGrange, B.W. Reed, M. J.Stern, N. D. Browning, G.H. Campbell, J. C Ki-effer, B.J. Siwick, and F. Rosei. Nanocrystal-lization of amorphous germanium films observedwith nanosecond temporal resolution.
AppliedPhysics Letters , 97(20):203102–203102–3, 2010. S. R. Stiffler, Michael O. Thompson, and P. S.Peercy. Nucleation of amorphous germaniumfrom supercooled melts.
Applied Physics Letters ,56(11), 1990. K. Lu. Nanocrystalline metals crystallized fromamorphous solids: nanocrystallization, struc-ture, and properties.
Materials Science and En-gineering: R: Reports , 16(4):161 – 221, 1996. D. Schwen, E. Martinez, and A. Caro. On the an-alytic calculation of critical size for alpha primeprecipitation in fecr.
Journal of Nuclear Materi-als , 439(1–3):180 – 184, 2013. L. Q. Chen. Phase-field models for microstruc-ture evolution.
Annual Review of Materials Sci- ence , 32:113–140, 2002. WJ Boettinger, JA Warren, C Beckermann, andA Karma. Phase-field simulation of solidifica-tion.
ANNUAL REVIEW OF MATERIALSRESEARCH , 32:163–194, 2002. Christian Miehe, Martina Hofacker, and FabianWelschinger. A phase field model for rate-independent crack propagation: Robust algo-rithmic implementation based on operator splits.
Computer Methods in Applied Mechanics andEngineering , 199(45–48):2765 – 2778, 2010. M. Koslowski, A.M. Cuiti˜no, and M. Ortiz. Aphase-field theory of dislocation dynamics, strainhardening and hysteresis in ductile single crys-tals.
Journal of the Mechanics and Physics ofSolids , 50(12):2597 – 2635, 2002. Laszlo Granasy, Tamas Pusztai, David Saylor,and James A. Warren. Phase field theory of het-erogeneous crystal nucleation.
PHYSICAL RE-VIEW LETTERS , 98(3), JAN 19 2007. CE Krill and LQ Chen. Computer simulationof 3-D grain growth using a phase-field model.
ACTA MATERIALIA , 50(12):3057–3073, JUL17 2002. J.W. Cahn, P. Fife, and O. Penrose. A phase-field model for diffusion-induced grain-boundarymotion.
Acta Materialia , 45(10):4397 – 4413,1997. Oliver Penrose and Paul C. Fife. Thermody-namically consistent models of phase-field typefor the kinetic of phase transitions.
Physica D:Nonlinear Phenomena , 43(1):44 – 62, 1990. J. W. Cahn and J. E. Hilliard. Free energy ofa nonuniform system. i. interfacial free energy.
Journal of Chemical Physics , 28(2):258 – 267,1958. Samuel M. Allen and John W. Cahn. A micro-scopic theory for antiphase boundary motion andits application to antiphase domain coarsening.
Acta Metallurgica , 27(6):1085 – 1095, 1979. In the derivation of this equation it has beenneglected the dependence of θ on ∇ φ . E.A. Carter, Giovanni Ciccotti, James T. Hynes,and Raymond Kapral. Constrained reaction co-ordinate dynamics for the simulation of rareevents.
Chemical Physics Letters , 156(5):472 –477, 1989. In our problem stresses are assumed to be zero,and therefore the Helmholtz free energy equalsthe Gibbs free energy at all temperatures. Defined as those where self-diffusion and/or self-relaxation processes occur on time scales beyondthe thermal equilibration time of the system. Mathias Rousset Tony Lelievre, Gabriel Stoltz.
Free Energy Computations: A MathematicalPerspective . Imperial College Press, 2010. J. M. Rickman and R. LeSar. Free-energy calcu-lations in materials research.
Annual Review ofMaterials Research , 32(1):195–217, 2002. Haluk Resat and Mihaly Mezei. Studies on freeenergy calculations. i. thermodynamic integra-tion using a polynomial path.
The Journal ofChemical Physics , 99(8), 1993. James N. Glosli and Francis H. Ree. The melt-ing line of diamond determined via atomisticcomputer simulations.
The Journal of ChemicalPhysics , 110(1), 1999. D. Frenkel and B. Smit.
Understanding Molec-ular Simulation: From Algorithms to Applica-tions . Academic Press, Orlando, FL, 2002. Jeremy Q. Broughton, James V. Lill, and J. KarlJohnson. c sphase diagram: A full free-energyanalysis. Phys. Rev. B , 55:2808–2817, Feb 1997. D. Y. Sun, M. Asta, and J. J. Hoyt. Ki-netic coefficient of ni solid-liquid interfaces frommolecular-dynamics simulations.
Phys. Rev. B ,69:024108, Jan 2004. H. Song and J.J. Hoyt. A molecular dynamicssimulation study of the velocities, mobility andactivation energy of an austenite–ferrite interfacein pure fe.
Acta Materialia , 60(10):4328 – 4335, 2012. M. Posselt and A. Gabriel. Atomistic simulationof amorphous germanium and its solid phase epi-taxial recrystallization.
Phys. Rev. B , 80:045202,Jul 2009. Kejian Ding and Hans C. Andersen. Molecular-dynamics simulation of amorphous germanium.
Phys. Rev. B , 34:6987–6991, Nov 1986. H. Tracy Hall. The melting point of germaniumas a function of pressure to 180,000 atmospheres.
The Journal of Physical Chemistry , 59(11):1144–1146, 1955. C. A. Angell. The amorphous state equivalentof crystallization: new glass types by first or-der transition from liquids, crystals, and biopoly-mers.
Solid State Sciences , 2:791?805, 2000. M. R. Gilbert, P. Schuck, B. Sadigh, and J. Mar-ian. Free energy generalization of the peierls po-tential in iron.
Phys. Rev. Lett. , 111:095502, Aug2013. Matteo Frigo and Steven G. Johnson. The designand implementation of FFTW3.
Proceedings ofthe IEEE , 93(2):216–231, 2005. Special issue on“Program Generation, Optimization, and Plat-form Adaptation”. Alain Karma and Wouter-Jan Rappel. Quanti-tative phase-field modeling of dendritic growthin two and three dimensions.
Phys. Rev. E ,57:4323–4349, Apr 1998. P. Germain, K. Zellama, S. Squelard, J. C. Bour-goin, and A. Gheorghiu. Crystallization in amor-phous germanium.
Journal of Applied Physics ,50(11), 1979. M. Avrami. Kinetics of Phase Change. II Trans-formationTime Relations for Random Distribu-tion of Nuclei.
The Journal of Chemical Physics ,8:212–224, 1940. (a)(b)(c) FIG. 19. Final nano-crystalline structure in-duced by spatially homogeneous nucleation andanisotropic growth conditions: (a) J = 9 × ,(b) J = 90 × and (c) J = 900 × s − m − .Each color corresponds to a given crystallographicorientation, which is randomly assigned at nucle-ation. L L AR = 1 AR = L /L AR = 1 AR = 1
FIG. 20. Anisotropy ratio (AR) for several illustra-tive shapes. a v e r age A R time / time for 95% crystallization J =10 −5 J =5 · −5 J =10 −4 J =5 · −4 J =10 −3 (a) a v e r age A R time / time for 95% crystallization J =9.00 · s −1 m −2 J =45.02 · s −1 m −2 J =90.03 · s −1 m −2 J =450.16 · s −1 m −2 J =900.32 · s −1 m −2 (b)(b)