Hydrodynamic theory of freezing: Nucleation and polycrystalline growth

TL;DR

This paper presents a hydrodynamic theory combining fluctuating hydrodynamics and phase-field crystal models to explore crystal nucleation and polycrystalline growth, revealing mechanisms for grain formation and spherulite-like structures.

Contribution

It introduces a novel hydrodynamic approach to simulate various nucleation processes and growth front phenomena in undercooled liquids, including grain orientation and defect formation.

Findings

Growth front nucleation leads to new grain orientations.

Dislocation formation occurs at the solid-liquid interface.

Density wave interference causes spherulite-like nanostructures.

Abstract

Structural aspects of crystal nucleation in undercooled liquids are explored using a nonlinear hydrodynamic theory of crystallization proposed recently [G. I. Toth et al., J. Phys.: Condens. Matter 26, 055001 (2014)], which is based on combining fluctuating hydrodynamics with the phase-field crystal theory. We show that in this hydrodynamic approach not only homogeneous and heterogeneous nucleation processes are accessible, but also growth front nucleation, which leads to the formation of new (differently oriented) grains at the solid-liquid front in highly undercooled systems. Formation of dislocations at the solid-liquid interface and interference of density waves ahead of the crystallization front are responsible for the appearance of the new orientations at the growth front that lead to spherulite-like nanostructures.

Equations

$$\frac{\partial\mathbf{p}}{\partial t}+\nabla\cdot(\mathbf{p}\otimes\mathbf{v})=\nabla\cdot\bigg{[}\mathbf{R}(\rho)+\mathbf{D}(\mathbf{v})+\mathbf{S}\bigg{]},\\$$

$$\frac{\partial \rho }{\partial t}+\nabla \cdot \mathbf{p}=0.$$

$$⟨S_{ij}^{\mathbf{r},t}{S}_{kl}^{{\mathbf{r}}^{\prime },t^{\prime }}⟩$$

$$\frac{F}{n_0{k}_{B}T}=\int dV\left\{f(n)-\frac{C_2}{2}(\nabla n{)}^2-\frac{C_4}{2}({\nabla }^{2}n{)}^2\right\},$$

$$f(n)=(1-C_0)\frac{n^2}{2}-\left(\frac{a}{2}\right)\frac{n^3}{3}+\left(\frac{b}{3}\right)\frac{n^4}{4}.$$

Full text

Crystal nucleation, i.e., stochastic formation of crystal grains via fluctuations that are able to grow plays an essential role in the making of polycrystalline and nanostructured matter ref1 ; ref2 , including the formation of polycrystalline growth forms such as disordered dendrites ref3 , and spherulites ref4 . The latter structures appear via the formation of new grains at the solidification front, a process termed “growth front nucleation” (GFN ref5 ). This phenomenon has been successfully addressed by conventional phase-field methods relying on coarse grained fields ref5 . Unfortunately, these works cannot provide details on the micromechanism of GFN. Although molecular simulations and theory recently shed light on many details of nucleation and liquid ordering preceding nucleation ref6 , GFN seems to be out of scope for such studies. Continuum theories working on the molecular scale offer a complementary approach to molecular simulations, and may deliver additional information on the relevant nanoscale processes, such as the mechanism of polycrystalline growth.

A fairly successful continuum approach, termed the Phase-Field Crystal (PFC) model, was developed recently to address the microscopic aspects of crystallization ref7 ; ref8 . The majority of the PFC studies were done assuming diffusive dynamics, which approximates reasonably crystalline aggregation in suspensions of micron size colloidal particles. In a recent work, we have proposed a hydrodynamic theory of freezing (HPFC) ref9 that applies for solidification in normal liquids. Our approach relies on fluctuating nonlinear hydrodynamicsref10 , and employs the free energy functional of the PFC model in determining the reversible stress tensor. This model recovers the proper dispersion relation for long wavelength acoustic phonons, a steady state front velocity, which is inversely proportional to the viscosity (as opposed to the time dependent front velocity observed in the case of diffusive dynamics), and describes the stress relaxation reasonably ref9 . It is thus expected to be able to capture defect formation, and therefore polycrystalline growth on the nanoscale.

Herein, we first employ the HPFC model for describing homogeneous and heterogeneous crystal nucleation, and then to the formation of new grains at the solidification front (GFN). While there were other PFC-based hydrodynamic models put forward recently ref11 , it is only the HPFC for which steady state growth, $v\propto\mu_{S}^{-1}$, and proper capillary wave spectrum were demonstrated.

In the HPFC model, we start from momentum transport and continuity equations used in fluctuating nonlinear hydrodynamics ref10 :

[TABLE]

[TABLE]

Here $\mathbf{p}(\mathbf{r},t)$ is the momentum, $\rho(\mathbf{r},t)$ the mass density, $\mathbf{v}=\mathbf{p}/\rho$ the velocity, $\nabla\cdot\mathbf{R}=-\rho\nabla\frac{\delta F[\rho]}{\delta\rho}\approx-\rho_{0}\nabla\frac{\delta F[\rho]}{\delta\rho}$ the divergence of the reversible stress tensor, $\frac{\delta F[\rho]}{\delta\rho}$ the functional derivative of the free energy with respect to density, $\rho_{0}$ a reference density, and $\mathbf{D}=\mu_{S}\{(\nabla\otimes\mathbf{p})+(\nabla\otimes\mathbf{p})^{T}\}+[\mu_{B}-\frac{2}{3}\mu_{S}]\mathbf{I}(\nabla\cdot\mathbf{v})$ the dissipative stress tensor, $\mu_{S}$ and $\mu_{B}$ are the shear and bulk viscosities, respectively, while the fluctuation-dissipation theorem yields the following covariance tensor for the momentum noise $\mathbf{S}$:

[TABLE]

Here $k_{B}$ is Boltzmann’s constant and $T$ the temperature.

We adopt the free energy functional of the PFC model with parameters used in Ref. ref9 :

[TABLE]

where the local free energy density is

[TABLE]

The order parameter is the normalized mass density: $n=\tilde{\rho}-1$, where $\tilde{\rho}=\rho/\rho_{0}$ ($\rho_{0}=m_{0}n_{0}$ is the reference density, where $m_{0}$ is the atomic mass and $n_{0}$ the number density of the reference liquid). $C_{0}$ can be related to the bulk modulus of the reference liquid via $K_{0}=(1-C_{0})(n_{0}k_{B}T)$, while $C_{2}$ and $C_{4}$ are responsible for elasticity ref9 . This formalism transforms into the usual Swift-Hohenberg formalism as follows ref8 : the reduced temperature and reduced density are expressed as $\epsilon=-(r-t^{2}/3)$ and $\psi=\phi+t/3$, where $r=4(1-C_{0})Z-1$, $\phi=n/X$, $t=-a\sqrt{3Z/b}$, $X=\sqrt{3/4bZ}$, and $Z=|C_{4}|/C_{2}^{2}$. The HPFC model differs from the original PFC model in only the dynamic equations.

The combination of fluctuating hydrodynamics with atomic scale theory is supported by recent results, which indicate that fluctuating hydrodynamics remains valid down to the nanoscale ref12 . To avoid interatomic flows in the crystal owing to the large density gradients, we employ coarse-grained momentum and density fields when computing the velocity field: $\mathbf{v}=\hat{\mathbf{p}}/\hat{\rho}$, an approximation used in the advective and viscous dissipation terms ref9 . While the HPFC model is not restricted to two dimensions (2D), an extensive testing was performed so far in 2D (see Ref. ref9 ).

Following Ref. ref9 , the properties of liquid iron were used to fix the model parameters for the 2D simulations. The kinetic equations were solved in 2D, using a pseudo-spectral scheme with a second order Runge–Kutta time stepping, while employing periodic boundary conditions on square grids of sizes ranging from $2048^{2}$ to $8196^{2}$. Accordingly, the presented results refer to a hypothetical 2D iron, which can in principle be realized in molecular dynamics simulations, and are expected to be relevant to crystallization in thin metal films. The melting point corresponds to $\epsilon_{L}=0.0923$, whereas the scaled liquid density is $\psi_{L}=-0.1982$. The linear stability limits taken at constant density or at constant temperature are $\epsilon_{c}=0.1178$ and $\psi_{c}=-0.1754$. The reduced temperature, and thus the undercooling, was tuned by varying $C_{0}$.

The structure of the solid phase is characterized by the number density map, whose peak positions are analyzed in terms of Voronoi polygons, and the bond-order parameter, $g_{6}=\sum_{j}exp\{i6\theta_{j}\}$, where $\theta_{j}$ is the angle corresponding to the $j$-th neighbor in the laboratory frame. The crystal grains in polycrystalline cases were identified on the basis of the orientation map obtained as the phase angle of the complex hexatic bond-order parameter $g_{6}$. The Voronoi polygons were colored gray, blue, yellow and red, when having 4, 5, 6, and 7 neighbors, respectively. $|g_{6}|$ characterizes the degree of disorder, and its phase specifies the local crystallographic orientation. Examples of these fields are displayed in Fig. 1.

The section of the phase diagram, where HPFC simulations were performed is shown in Fig. 2. First, we studied growth initiated by a small (atom size) potential well. Close to the liquidus circular crystals grow, whereas approaching the stability limit, the growth form becomes hexagonal first with rounded corners, evolving into hexagons with pointed corners and concave edges. This behavior is attributed to a change of the interface structure: The respective density distributions are displayed in Fig. 3. Panels (a) and (b) show a closeup of the growth fronts, whereas panels (c) and (d) present the respective density profiles, and the time dependencies of the average crystal radius. Near equilibrium, the solid-liquid interface extends to 6-7 interatomic distances, becoming considerably sharper, when approaching the stability limit. The respective change in anisotropy is expected to be responsible for the different crystal shapes. Apart from an initial transient, the crystals display essentially linear growth as predicted in ref9 .

Crystallization kinetics: The momentum noise in Eq. 1 gives rise to density fluctuations, which together with molecular scale density waves from the reversible stress tensor, lead to the formation of crystal-like fluctuations (homogeneous nucleation), followed by crystal growth, yielding to polycrystalline freezing. We investigated this in the metastable liquid domain slightly backward from liquid instability: $\epsilon=0.1158$ and $\psi_{L}=-0.1982$. The results are summarized in Fig. 4. The Johnson-Mehl-Avrami-Kolmogorov expression, $Y(t)=1-exp\{-[(t-t_{0})/\tau]^{p}\}$ ref2 was fitted to the temperature dependent crystalline fraction evaluated from the number of density peaks. Here $t_{0}$ is the incubation time, $\tau$ a characteristic time related to the nucleation and growth rates, and $p$ is the Avrami-Kolmogorov exponent indicative to the mechanism of crystallization ref2 . In our case, from the average slope of the $ln\{-ln(1-Y)\}$ vs. $ln(t-t_{0})$ plot, we obtained a $p=3.31\pm 0.03$ [Fig. 4] implying linear growth with slightly increasing nucleation rate. This result differs considerably from the $p$ strongly decreasing with increasing $Y$ reported for diffusive dynamics ref13 .

Particle induced ’nucleation’: To model a foreign particle (substrate), an extra term, $V(\mathbf{r})\psi(\mathbf{r},t)$, was added to the free energy density, where $V(\mathbf{r})=0$ outside the substrate, whereas it is a periodic potential inside that determines the crystal lattice of the substrate. To test whether the free growth limited model of Greer and coworkers ref14 remains valid in the presence of hydrodynamics, we employ a square shape substrate of made of a square lattice, whose lattice constant coincides with that of the forming triangular crystal. This approximates the ideal wetting Greer and coworkers assumed between the substrate and the crystal. The undercooling has been increased by multiplying $C_{0}$ by factors $\xi=0.999,0.998,0.995,$ and $0.990$. In Fig. 5, we present long-time solutions as a function of $\xi C_{0}$. In the HPFC model (as well as in the PFC model with diffusive dynamics ref15 ), there exists a critical undercooling for a given particle size, beyond which free growth takes place.

Growth Front Nucleation (GFN): Formation of new grains at the solidification front has been identified as the mechanism by which complex polycrystalline growth forms appear ref3 ; ref4 ; ref5 . This phenomenon has been successfully modeled by phase-field methods employing orientation fields to monitor the local crystallographic orientation. In these models, new grains form either by quenching orientational defects (bundles of dislocations) into the crystal (at large undercoolings), or via branching in directions of low grain-boundary energies (at small undercoolings). The orientation field approach, became fairly successful in capturing complex growth structures ref5 . Yet it is desirable to clarify the microscopic background of GFN. It is suspected that in substances of different molecular structures/interactions different GFN mechanisms may occur. Herein, we employ the HPFC model to study the formation of new grains at the solid-liquid interface of 2D hexagonal crystals at high undercoolings.

We made our first attempts to model GFN years ago ref16 using the original PFC model. At supersaturations beyond the liquid stability limit $[\psi>\psi_{c}=-(\epsilon/3)^{1/2}]$, we observed that without noise crystal seeds evolved into ordered polycrystalline structures. The growing crystal was surrounded by concentric density waves, which initiated crystallization accordingly: in six directions these waves helped the growth of the original crystal, whereas in other directions a large number of defects formed and new orientations appeared that fitted better to the local direction of the density waves (Fig. 6). Later works performed without noise in the unstable liquid regime indicate the transition of a flat single crystal front into polycrystalline growth and eventually into glassy freezing with increasing driving force ref17 . These findings could be regarded as signs of elementary processes of GFN. However, when a noise term obeying the fluctuation-dissipation theorem is added to the equations of motion in the unstable regime, these solutions are suppressed by explosive nucleation. Inside the metastable regime, we were unable to observe GFN, probably owing to the lack of an extended layer of molecular scale density waves ahead of the front.

In the HPFC model, we were able to observe polycrystalline growth forms in the metastable liquid (see Fig. 7). The kinetic equations were solved on a $2048^{2}$ rectangular grid. A shallow molecular size potential well was used to initiate freezing. It induced concentric density waves (akin to the ’onion structures’ predicted in ref18 , then a small hexagonal single crystal formed, but as it grew new orientations appeared gradually via two mechanisms of GFN: (i) Dislocations entered in the growing crystallites at cusps centers. These appear due to the interaction of the stress field of the growing nanoscale crystallite with density fluctuations. (ii) Small crystallites nucleate in the close vicinity of the solid-liquid interface, which apparently originate from the interference of the density waves ahead of the the rough solid-liquid interface (see Fig. 8). The two mechanisms are clearly visible in the snapshots of the density, orientation and Voronoi maps (Figs. 7 and 8).

We found that mechanism (i) occurs in a relatively broad range of undercoolings or densities, however, the formation rate of dislocations decreases with decreasing driving force, a finding that might be associated with a decreasing growth velocity (see the discussion below). In contrast, mechanism (ii) can only be observed in the close vicinity of the linear stability limit of the liquid. For example, at the reduced temperature $\epsilon=0.0923$, mechanism (ii) appears in the scaled density range $-0.1778<\psi_{0}<\psi_{c}=-0.1754$; whereas at fixed scaled density of $\psi_{0}=-0.1982$, nucleation ahead of the front occurs in the reduced temperature range $0.1158<\epsilon<\epsilon_{c}=0.1178$. This observation correlates with the finding that in the HPFC model, the thickness of the liquid layer, in which liquid ordering in the form of density waves takes place, increases towards the stability limit.

Although our simulations are on the nanoscale, we observe spherulite-like structures (Fig. 9) similar to those seen on larger scales. We note that polycrystalline growth structures can form spontaneously on this scale as shown by molecular dynamics simulations ref19 and experiments on organic and metallic nanospherulites ref20 . Our work might be directly relevant to them and to carbon nanostructures.

An intriguing question is, why the HPFC model recovers GFN in the metastable liquid domain, whereas the diffusive PFC model does not. This might be related to either the differences in the interface structure or defect dynamics, or in both. For example, in the case of diffusive PFC, a fast diffusionless growth mode characterized by a broad interface occurs at large driving forces (still in the metastable liquid) ref21 , in which healing of the defects can be relatively easy, avoiding the formation of dislocations at the perimeter of the growing crystal, preventing thus GFN. In the HPFC model, the relatively sharp inter face may be unable to do this. We discuss these issues below.

The original PFC model incorporates elastic interactions as well as crystal plasticity combined with diffusive dynamics, however without hydrodynamic modes of the evolution of the density field. Defect dynamics, including vacancy motion, dislocation glide, climb, and annihilation, and grain boundary melting were studied extensively on the diffusive time scale ref7 ; ref8 ; ref22 . A few studies went beyond this approximation by incorporating faster processes enabled by the inclusion of a second order time derivative into the equation of motion (see the MPFC model in Ref. ref23 ), or via linearizing the Navier-Stokes equations ref24 . In both of these quasi-hydrodynamic models the atomic positions are relaxed rapidly at early times in a manner consistent with elasticity theory, while the late time defect motion, including vacancy diffusion, grain boundary kinetics, and dislocation climb, is governed by diffusive dynamics. This suggests that defect dynamics alone may not be responsible for the observed differences in GFN.

A work by Majaniemi et al. ref24 explored differences between dynamics of mass distribution in the MPFC model (they termed it Type-2 model) and a more sophisticated linearized hydrodynamic model (termed Type-3 model) in the non-equilibrium case of a crack relaxation under stress. They found that the dynamic mass distribution the two models predict can be quite different, which they attributed to the transport differences of the quasi-phonons appearing in the MPFC model and the acoustic phonons of the Type-3 model. Remarkably, in the Type-3 linearized hydrodynamics model, which is close to the present full hydrodynamic model, the system relaxed to a multigrain structure, as opposed to the single crystal solution from the MPFC model. The behavior predicted by the Type-3 model appears to be similar to what we observed in our full hydrodynamic HPFC simulations, in which, after engulfing low density fluctuations into the solid, dislocations and new grains were formed at the interface. This phenomenon is more frequent in the cusps of the interface, and the dislocations may be of misfit origin.

Apparently, we have stress at the growth front already in the early single crystal state, as reflected by the non-uniform orientation field that indicates slight misorientations (relative rotation of the Voronoi cells) between the two sides of the corners of the crystallite and at the central region of its sides. To see this clearly, we performed simulations without momentum fluctuations for the highly nonequilibrium liquid, starting crystallization with a small cluster that evolves due to a weak potential of 2D hexagonal structure. The potential contains a central well and 6 neighboring wells of equal depth given by the single mode solution for the particle density for this structure.

The orientation field of a crystallite, grown from such a cluster, is shown in Fig. 10(a). The deformation and the related stress that causes the inhomogeneity of the orientation field appears to originate from a radial density change shown in Fig. 10(b). The latter panel displays the coarse grained particle density obtained by FIR (Finite Impulse Response) filtering ref25 from the particle density. The spatial modulation of coarse grained density at the interface might indicate a weakly oscillating growth velocity. Despite the presence of the hydrodynamic mode of density relaxation, density depletion is observed ahead of the solidification front, which is probably responsible for the curving of the sides of the crystal, yielding orientation and stress fields antisymmetric to the lines across the tips and centers of the opposite sides. An essentially similar, though more noisy coarse grained distribution is obtained in the presence of momentum fluctuations [see Fig. 11(a)].

We observe the formation of dislocations at the interface both with [Fig. 11(b)] and without momentum fluctuations, however at a much later stage in the latter case: In the simulations without noise, we do not observe dislocations below an equivalent radius ($R_{c}=2\sqrt{2\pi N/3}\approx 450$, where $N$ is the number of particles in the crystallite, see Figs. 10). Beyond this size the stress at the interface is large enough to initiate the formation of cusps, and the nucleation of dislocations in them. A possible mechanism for this can be the Asaro-Tiller-Grinfeld instability ref26 ; ref27 . We find that these misfit dislocations appear at a much smaller size in the presence of momentum noise ($R_{c}\approx 150$). These findings indicate that the fluctuations play a key role in the formation of defects. This conclusion is further supported by a recent PFC study, which has shown that the formation of misfit dislocations is helped by an increasing strength of the density fluctuations ref27 . We note furthermore, that the noise influences the appearance of dislocations, which leads to symmetry breaking ref27 . Apparently, in the presence of density fluctuations the dislocations appear fairly randomly, although normally in cusps forming at the interface.

To test the role of fluctuations further, we have determined the wavenumber spectrum of the density fluctuations in the HPFC model emerging from the momentum fluctuations, and included the same noise in the equation of motion of the diffusive PFC model under the same $\epsilon$ and $\psi_{0}$ values (under these conditions the diffusionless ’fast growth mode’, takes place in the diffusive PFC model ref21 ). Despite these, no GFN was observed in the diffusive PFC model. This observation suggests that under equivalent conditions the diffusive PFC model is less susceptible to defect formation than the HPFC model.

It appears that a combination of faceted growth, a weak depletion at the interface, and the presence of density fluctuations is needed for initiating GFN. Work is underway to investigate further the microscopic aspects of these phenomena.

Summarizing, we applied the HPFC model to solidification problems in 2D. We demonstrated that

(i) radial growth happens at a steady state rate,

(ii) crystallization takes place via homogeneous nucleation and steady state growth, and the respective Avrami-Kolmogorov exponent is $p=3.31\pm 0.03$,

(iii) the free growth limited mode of particle induced crystallization ref14 remains valid down to the nanoscale for simple liquids, and

(iv) the model predicts two modes for growth front nucleation.

While (i) to (iii) are important demonstrations of the capabilities of the model, the most significant result is (iv): The HPFC appears to be the first atomic scale model that yields growth front nucleation in the metastable liquid regime. The identification of the basic mechanisms of growth front nucleation in simple liquids might help to control microstructure evolution in such systems.