Phase-Field Simulations at the Atomic Scale in Comparison to Molecular Dynamics

Early solidification is investigated using two different simulation techniques: the molecular dynamics (MD) and the phase-field (PF) methods. While the first describes the evolution of a system on the basis of motion equations of particles, the second grounds on the evolution of continuous local order parameter field. The aim of this study is to probe the ability of the mesoscopic phase-field method to make predictions of growth velocity at the nanoscopic length scale. For this purpose the isothermal growth of a spherical crystalline cluster embedded in a melt is considered. The system in study is Ni modeled with the embedded atom method (EAM). The bulk and interfacial properties required in the PF method are obtained from MD simulations. Also the initial configuration obtained from MD data is used in the PF as input. Results for the evolution of the cluster volume at high and moderate undercooling are presented.


Introduction
The phase-field is a powerful method to describe solidification phenomena [1,2] on the mesoscopic length scale. It has been used to model homogeneous and heterogeneous nucleation, microstructure formation in solids, and motion of grain boundaries. The phase-field models includes formulations for pure substances [3], for multicomponent systems [4], and for polycrystalline structure [5] and solidification in eutectic [6,7], peritectic [7], and monotectic [8] systems. Other simulation techniques for dendrite growth are cellular automata [9] and hybrid methods such as the multiscale diffusion Monte Carlo (DMC) [10,11].
The phase-field method requires previous knowledge of the material properties of the system in study. The input includes bulk properties such as density, heat capacity and latent heat, and others such as interfacial and kinetic growth coefficients, being the latter properties which are hardly accessible in experiments. Here molecular simulations play a fundamental role, since they provide a link between an interaction potential and all the required properties. Kinetic coefficients, interfacial free energies, and their dependence on the crystal orientation, of pure metals and alloys, can be directly obtained from simulations of inhomogeneous liquidcrystal systems [12][13][14][15].
There are previous phase-field studies in which the material parameters have been obtained from molecular simulations to model dendritic growth in pure Ni [16,17], CO 2 hydrates [18], and binary alloys. In this work we compare MD and PF simulation results for the isothermal growth velocity of a nanoscopic spherical crystalline cluster of Ni embedded in the melt. For this purpose we use precise thermodynamic and kinetic data obtained from molecular dynamics and use equivalent initial configuration in the phase-field simulations.

Properties of Ni (EAM)
The properties of Ni obtained from molecular dynamic simulations are given in this section. The embedded atom potential of Foiles [19] is used. The thermal dependence of bulk properties is determined from independent simulations of equilibrated homogeneous systems at pressure = 0. Liquid phase is prepared by melting of an initial crystalline phase, by setting the temperature of the thermostat well above the melting temperature. Crystalline phase is obtained by 2 The Scientific World Journal relaxation of the perfect crystalline phase. The samples are prepared in NpT simulation runs and the production in NVE ensemble.
Densities of the bulk phases are fitted to the function ( ) = + ⋅ + ⋅ 2 , for the fcc solid phase ( = s) = 8901.6 kg/m 3 , = −0.20379 kg/m 3 K, and = −0.0000614202 kg/m 3 K 2 while for the liquid phase ( = ℓ) = 8992.26 kg/m 3 , = −0.667037 kg/m 3 K, and = 0.0000331612 kg/m 3 K 2 . The parameters were obtained in the temperature intervals 1000 K to 3000 K for the liquid and 300 K to 1900 K for the solid. In both cases parts of the metastable phases are included in the fit procedure.
The latent heat is obtained from the difference between the enthalpy of the liquid and the solid and fitted tõ0( ) = + ⋅ + ⋅ 2 with = −15980.2 J/kg, = 324.745 J/kg K, and = −0.0810984 J/kg K 2 . These values are determined in the range from 1000 K to 1900 K. The latent heat expressed in units of energy per volume of the solid phase is given by The specific heat capacity at constant volume is given bỹV ( ) = + ⋅ where = 419.452 J/kg K and = 0.020388 J/kg K 2 for the solid and = 563.024 J/kg K and = −0.06952 J/kg K 2 for the liquid. The thermal diffusivity for the model is 2.1⋅10 −7 m 2 /s [20] while the experimental values range from 120 ⋅ 10 −7 m 2 /s [21] to the recent value 8.7⋅10 −7 m 2 /s [22] (obtained in microgravity conditions). The value obtained from simulations is at least 4 times lower than that in experiments, and the difference is due to the electronic contribution to thermal transport in metals which are not accounted in classical molecular dynamics.
The melting temperature of the model is precisely estimated from simulations of inhomogeneous solid-liquid system at the point where the velocity of growth is zero. For the model = 1748 K. The kinetic coefficient is estimated from linear relationship between planar growth velocity and undercooling close to coexistence. Growth velocities are obtained from independent simulations for the crystal orientations [100], [110], and [111]. The kinetic growth coefficient depends also on the orientation of the crystal. Its orientational dependence is represented here by the expansion with = 4 1 + 4 2 + 4 3 , = 6 1 + 6 2 + 6 3 , and = 2 The interfacial stiffness and interfacial free energy were obtained from the analysis of capillary waves spectrum of large solid-liquid interfaces at coexistence (i.e., at = for = 0). The orientational dependence of the interfacial free energy is described by a cubic harmonic expansion with 0 being an orientational averaged value of the interfacial free energy, while the coefficients describe the strength of the anisotropy. Interfacial stiffness is defined as̃= + 2 /̂̂, where and indicate unit vectors tangent to the interface plane [23]. The parameters of (2) are obtained from the system of equations for the stiffnesses̃and their numerical values are obtained from the analysis of the capillary waves spectrum [14].

Phase-Field Model for Pure Material
Isothermal crystal growth for pure material is modeled by the variables of the internal energy density and of two phasefields 0 ≤ s ≤ 1 (solid) and 0 ≤ ℓ ≤ 1 (liquid). The phase-field variables fulfill the constraint s + ℓ = 1, so a single phase-field variable = s is sufficient to describe the evolution of the phase boundaries in the system. Note that ℓ = 1 − . The variable ( ⃗ , ) denotes the local fraction at ⃗ of the considered phase at time .
The phase-field model is based on an entropy functional to ensure consistency with classical irreversible thermodynamics The bulk entropy density depends on the internal energy density and the phase-field variable ( ⃗ , ). The functions ( , ∇ ) and ( ) reflect the thermodynamics of the interfaces and = 4 is a small length scale parameter related to the thickness of the diffuse interface.
The Scientific World Journal 3 From (3) one derives the equations for the nonconserved phase-field variable by taking the functional derivatives in the form where = ( , ∇ ) is an anisotropic relaxation parameter dependant on temperature, (̂) = ( )/ (̂) with (1), and it is related to the kinetic coefficient 0 .
The evolution of the phase-field variables is described by The thermodynamic relation = + gives us ( , )/ = −(1/ )( ( , )/ ), detailed in [24]. Here is the bulk free energy density where ℎ( ) is a monotonous function on is the melting temperature. For the isothermal case we can define two constants ( , ) =: ∑ ∈{s,ℓ} ℎ( ) and consider here that the free energy contribution to the equation of motion (6) Equation (6) is solved by a time dependent forward euler scheme with a second order spatial discretisation on a regular Cartesian grid.

Preparation of the Sample.
Growth simulations begin from a configuration made of a stabilized crystalline spherical cluster embedded in the melt. The preparation of the initial configuration involves a sequence of NpT simulation runs. First the crystalline solid phase is equilibrated at = 0 starting from atoms arranged in a perfect fcc structure in a cubic box with an initial density close to the experimental value at room temperature ( = 8.9 g/cm 3 ). The temperature in this step is chosen equal to the one used later in the simulation of growth. In the second step two regions are defined in the simulation box: the melting zone and the atoms in the crystalline phase (which is chosen here as a sphere of radius = 25Å). The region outside the crystalline zone is melted by increasing the temperature well above the melting point while the atoms in the solid region are constrained to remain at fixed positions. The temperature is linearly increased while an isotropic barostat is applied. In the next step, the stable melt at high temperature is cooled down to the temperature of the crystal. The melt reaches a metastable (undercooled) state. At the end of this run, a crystalline cluster embedded in the liquid phase at the same temperature is obtained. In the last step, the relaxation of the atoms at the interface is made by an NpT simulation of some picoseconds in which all particles are allowed to move.
In order to obtain an equivalent initial configuration for the PF method, the local order parameter values given at atom positions are mapped on a regular grid with cell-width of 1Å. The order parameter in a given grid point is given by the average of all particles contained in cell u = ( , , ) and its first two neighbor cells k = ( , , ) with |u − k| < 2.
The order parameter values obtained from MD simulation were scaled so that the volume of the crystal cluster is equal in the coarse-grained and in the original initial configurations. The distribution of the order parameter in the bulk phases is Gaussian centered at with standard deviation . The definition of order parameter in the phase-field method is as follows: = 0 if 6 6 < := ℓ − 2 ℓ , and = 1 if

Planar
Growth. An analytic formulation of the phasefield method in the sharp interface limit assumes a temperature gradient at the interface. The formulation reads The Scientific World Journal here Δ denotes the undercooling in the interface, V is the velocity, and is a kinetic coefficient. The curvature of a two dimensional interface is the sum of the two principal curvatures 1/ 1 , 1/ 2 ; so = 2/ if 1 = 2 = . The Gibbs-Thomson-coefficient Γ = / with surface tension and being the latent heat.
In equilibrium, that is, for V = 0, (11) reads as Δ = Γ = Γ/ (in 2D). It follows that the critical radius = Γ/Δ . A small spherical cluster with radius is in equilibrium. For example, at = 1550 K the critical radius = 9.7Å. A cluster smaller than melts. The phase-field simulations corroborate this result: for = 1550 K the nucleus grows for = 10Å and melts for = 9.5Å.
For a planar interface (11) is simplified to Δ = V, taking account of the absence of curvature. The kinetic coefficient adapted from MD is = 1/ . So the phase-field relaxation parameter is Simulations of a planar front in comparison to the analytic solution and molecular dynamics simulations are shown in Figure 1. Additionally we list the kinetic coefficient in Table 1. Both methods converge at low undercooling to the solution given by (11). At higher undercooling velocities obtained from PF simulations still exhibit a linear dependence with undercooling Δ ; this is given in the range from 1400 K to ; the fit of the kinetic coefficient is made in this range and agrees with the expected kinetics from (1).
The MD simulations indicate that this linearity is not more valid; the velocity is lower than that expected from extrapolation of the linear relation V = Δ . This effect can also be recognized in the results of Hoyt et al. [15]. Many effects influence the growth velocity of MD at higher undercooling, with the increase of order in the liquid bulk being the most important which leads to the formation of small crystalline clusters in the bulk liquid. The phase-field shows this nonlinearity only for very high undercoolings Δ > 500 K.

Spherical Cluster Growth
Growth simulations starting from the samples prepared as described in Section 4.1 are performed using both MD and PF methods. At microscopic level NpT simulations are performed at = 0 using the Andersen thermostat and barostat. The number of particles is = 256000 which corresponds to a cubic box of length of about 145Å. The initial radius of the spherical crystalline cluster is 25Å which corresponds to about 5.5 ⋅ 10 3 particles (see Figure 2(a)); that is much larger than the size of the critical cluster at = 1500 K which is made of about 10 3 particles [27]. For ≥ 1600 K the initial size of the cluster is lower than the critical value and it melts as expected. PF simulations are performed imposing the isothermal condition starting with the same configuration as in MD but converted as described in Section 4.2. Results for the evolution of the cluster volume are shown in Figure 3; snapshots from the MD and PF simulations are shown in Figure 2. At moderate undercooling ( > 1550 K) the MD and PF show comparable growth rates. For higher undercooling the growth in the phase-field model is notably faster than in the MD method. The discrepancy becomes more pronounced at lower temperature and higher simulation time. Assuming that the cluster conserves its spherical shape during the first stage of growth we define the radius as = (3 /4 ) 1/3 and the radial velocity as V = / . Figure 4 shows the evolution of radial velocity in the MD simulation; three regimes can be identified; initially the velocity tends to decrease, then exhibits a linear increase, and finally decreases again. The first decrease is due to the preparation method of the sample; the system relaxes while the interface is formed. Growth and relaxation cannot be decoupled; both occur simultaneously in this initial regime. In the second regime the growth velocity increases linearly in time. The decrease of velocity in the last regime is due to the finite size of the sample; when the cluster is large enough it interacts with its periodic images through the periodic boundaries of the simulation box and loses its spherical shape. In the PF simulation no extra relaxation regime is observed; the dynamic is determined by the initial definition after the interface is relaxated. Moreover, in the PF the cluster is not affected by finite size effects as in the MD simulation since no periodic boundaries are required. In comparison the PF exhibits an apparent linear increase from the beginning. This means that the starting radius of the nucleus is not the same for both simulation methods, so a comparison for the same radius instead of the time is more adequate. Figure 5 shows the radial velocity versus cluster radius at different temperatures. Results obtained from PF simulations are well described by Gibbs-Thomson equation (11). The analytical velocity is V = avg (Δ − 2Γ/ ). So we use ( ) = + / to fit the data with the parameters and for each temperature . We applied the fit in the region starting with = 30Å, where the PF has relaxated the interface. These fits are also applied for the MD results in the second regime of growth, that is, for crystalline clusters of size of about 40Å ≾ ≾ 60Å (the points in Figure 5). According to the Gibbs-Thomson equation the radial velocity reaches an asymptotic value (V ,∞ = V , → ∞ ), an ideal situation where the cluster is large enough and conserves its shape. This limit cannot be directly obtained from simulations first because the simulations are performed in a finite volume and also because the crystalline cluster would adopt a nonspherical shape due to the anisotropy of the interfacial energy and growth coefficient. We estimate this hypothetical limit as a basis of comparison between the PF and MD simulations. For this purpose we extrapolate the fit and obtain V ,∞ = ; one value is obtained at each temperature. The asymptotic The Scientific World Journal radial velocities obtained from the PF and MD simulations are shown in Figure 6. The results of both methods lie between the linear relations V 100 = 100 Δ and V 111 = 111 Δ . For the MD the results show a similar behavior as observed in the planar growth case; the deviation from the linear relation at high undercooling observed in the PF simulations is due to the correction introduced in the kinetic coefficient. As shown in Section 4.3 the crystal grows in different crystal orientations with a different velocity. In the PF the kinetic anisotropy is responsible for this. The planar front simulations from Figure 1 show that the kinetics are very well reproducible from the PF simulations. If we assume that we have a perfect sphere the mean kinetic can be calculated as a surface area integral of the unit sphere over the kinetic coefficient; we get avg = ∫ 2 (̂) / ∫ 2 1 ≈ 0.175. So V(Δ ) = avg Δ is the excepted growth velocity for a sphere; it is shown as the solid line in Figure 6.
At this early stage of growing the crystal is very similar to a sphere. The preferred growing in [100] directions and the surface anisotropy formed the crystal to its equilibrium shape; for Ni this is similar to an octahedron; in the limit this shape has only [111] directions. So the area of the [100]   area measurement method assumes a smooth nondeformed interface; this is not given for PF interfaces if there are driving forces like curvatures and undercoolings. We use a factor estimated for a sphere of radius 25Å to correct this. The surface area is multiplied by this factor 1.089. Figure 7 shows the sphericity over the radius; the crystal becomes non spherical for longer growing and faster for higher undercoolings. This means that the velocity of growing is below avg especially for higher undercoolings. The nonspherical shape is shown in Figure 2(h).

Conclusion
Many have shown the linking between MD and PF [16,17]. So it is a common way to obtain MD parameters for PF simulations. The PF simulations on the natural PF scale, for example, mesoscopic scale, produce correct results compared 8 The Scientific World Journal with experiments [29,30]. We have confirmed that a PF with MD parameters is correct even on the atomistic scale. This is not surprising, since we are very close to the sharp interface limits. The mesoscopic view at the atomistic scale provides the same growth rates as an atomic simulation method, although the atoms are not resolved in the PF. The fluctuations of the interface on the atomistic length scale will be investigated in a consecutive work. However, you will not simulate too high undercooled melts, where the mesoscopic view of the PF will not consider all the nonlinearity from the MD interface. We have used the sphericity to measure the shape and explain the discrepancy to the analytic linear velocity.
While it is known that PF works correctly on mesoscopic scale and with this study it is shown that PF also works on the atomistic scale, one is able to do multiscale simulations starting on the atomistic scale with one simulation method.