Solidification in soft-core fluids: Disordered solids from fast solidification fronts

Using dynamical density functional theory we calculate the speed of solidification fronts advancing into a quenched two-dimensional model fluid of soft-core particles. We find that solidification fronts can advance via two different mechanisms, depending on the depth of the quench. For shallow quenches, the front propagation is via a nonlinear mechanism. For deep quenches, front propagation is governed by a linear mechanism and in this regime we are able to determine the front speed via a marginal stability analysis. We find that the density modulations generated behind the advancing front have a characteristic scale that differs from the wavelength of the density modulation in thermodynamic equilibrium, i.e., the spacing between the crystal planes in an equilibrium crystal. This leads to the subsequent development of disorder in the solids that are formed. In a one-component fluid, the particles are able to rearrange to form a well-ordered crystal, with few defects. However, solidification fronts in a binary mixture exhibiting crystalline phases with square and hexagonal ordering generate solids that are unable to rearrange after the passage of the solidification front and a significant amount of disorder remains in the system.

Figure 1

The radial distribution function $g\left(r\right)$ for a GEM-4 fluid with bulk chemical potential $\mu =0$ obtained from the HNC closure to the OZ equation (dashed lines) and from the RPA DFT via the test particle method (solid lines) for several values of $\beta \epsilon$. For clarity, the results for $\beta \epsilon =1$ and 5 have been shifted vertically. The results correspond to the state points $(\beta \epsilon ,{\rho }^b{R}^2)=(1,0.36)$, $(5,0.14)$, and $(10,0.088)$. As $\beta \epsilon$ increases, the RPA approximation becomes increasingly poor; nevertheless, even for (fairly low density) state points such as $\beta \epsilon =10$ the agreement is surprisingly good; recall that the RPA approximation improves as the density is increased.

Figure 2

Phase diagrams of the one-component 2D GEM-4 and GEM-8 model fluids. The solid lines are the binodals, i.e., loci of coexisting liquid and solid phases. The dashed lines are the spinodal-like instability lines along which the metastable liquid phase becomes linearly unstable.

Figure 3

Equilibrium density profile at the free interface between coexisting liquid and solid phases in the GEM-4 model when $\beta \epsilon =1$ and $\beta \mu =17.0$.

Figure 4

The front speed (a) as a function of the density of the metastable liquid into which the front propagates and (b) as a function of the chemical potential for a GEM-4 fluid with temperature $k_{B}T/\epsilon =1$. The red solid line is the result from the marginal stability analysis and the black dashed line is the result from numerical computations from profiles such as that displayed in Fig. 5. The black circles denote (a) the densities ${\rho }_l$, ${\rho }_s$ at liquid-solid coexistence and (b) the coexistence value $\beta \mu \approx 17.0$.

Figure 5

Density profile across a solidification front advancing from left to right into an unstable GEM-4 liquid with bulk density $\rho R^2=8$ and temperature $k_{B}T/\epsilon =1$, calculated from DDFT. The top panel shows the full 2D density profile $\rho (x,y)$ while the panel below shows the 1D density profile $\rho \left(x\right)$ obtained by averaging over the $y$ direction, parallel to the front. The bottom panel shows $ln\left(\right|\rho \left(x\right)-{\rho }_0\left|R^2\right)$ in order to reveal the small-amplitude oscillations at the leading edge of the advancing front.

Figure 6

The wave number $k^{*}$ of the stripe state produced behind the front as a function of density for the GEM-4 fluid with $\beta \epsilon =1$, obtained from Eq. (29) together with the wave numbers $k_{\mathrm{r}}$ of the 1D oscillations at the leading edge of the front and $k_{\mathrm{eq}}\equiv 2\pi /\lambda$, where $\lambda$ is the distance between lattice planes in the equilibrium hexagonal state. This wavelength is very different from the wavelength of the oscillations produced by the advancing front, $2\pi /k^{*}$.

Figure 7

Density profiles obtained from DDFT for an unstable GEM-4 fluid with bulk density ${\rho }_0{R}^2=8$. To facilitate clear portrayal of the front structure we plot the quantity $ln\left(\right|\rho \left(\mathbf{r}\right)-{\rho }_0\left|R^2\right)$. Solidification is initiated along the vertical line $x=0$ at time $t^{*}=0$. This produces two solidification fronts, one moving to the left and the other to the right, moving away from the line $x=0$. The upper profile is for the time $t^{*}=1$ and the lower for $t^{*}=1.4$. We see significant disorder as the front creates density modulations that are not commensurate with the equilibrium crystal structure.

Figure 8

Top panel: The angle distribution $p\left(\theta \right)$ at times $t^{*}=2.2$, 3.2, and 4.4 after the initiation of a solidification front for a GEM-4 fluid with bulk density ${\rho }_0{R}^2=8$ (cf. Fig. 7) computed from the triangles of a Delauney triangulation on the density peaks of the profile from DDFT (middle panel: $t^{*}=2.2$; bottom panel: $t^{*}=4.4$).

Figure 9

The linear stability limit for a binary mixture of GEM-8 particles with $\beta \epsilon =1$ and $R_{22}/R_{11}=1.5$ and $R_{12}/R_{11}=1$, plotted in the total density $\rho \equiv {\rho }_1+{\rho }_2$ vs concentration $\varphi \equiv {\rho }_1/\rho$ plane. The black circles denote the state points corresponding to the density profiles displayed in Fig. 10.

Figure 10

Equilibrium crystal structures for a GEM-8 binary mixture with $\beta {\epsilon }_{ij}=\beta \epsilon =1$ for all $i,j=1,2$, $R_{22}/R_{11}=1.5$, $R_{12}/R_{11}=1$, with total average density $\rho R_{11}^2=4$ and concentrations (a) to (e) $\varphi =0$, 0.1, 0.25, 0.5, and 0.9. The structures are shown in terms of the quantity $[{\rho }_1\left(\mathbf{r}\right)-{\rho }_2\left(\mathbf{r}\right)]R_{11}^2$, with regions where ${\rho }_1\left(\mathbf{r}\right)>{\rho }_2\left(\mathbf{r}\right)$ colored black. All profiles correspond to local minima of the free energy, but we have not checked whether they correspond to global minima at the given state points. We observe a binary square lattice structure in (c), a binary hexagonal lattice structure in (b) and (d), and a simple hexagonal lattice in (a) and (e), where the minority species particles occupy the same lattice sites as the majority species particles, in contrast to the lattice structures in (b)–(d). The density profiles of species 1 or 2 in case (e) are very similar to the profile shown, the only difference being the height of the density peaks.

Figure 11

Snapshots of a solidification front in a GEM-8 mixture with $\beta {\epsilon }_{ij}=\beta \epsilon =1$ for all $i,j=1,2$, $R_{22}/R_{11}=1.5$, and $R_{12}/R_{11}=1$, advancing from left to right into an unstable fluid with $\rho R_{11}^2=8$ and $\varphi =0.5$, in terms of the quantity $[{\rho }_1\left(\mathbf{r}\right)-{\rho }_2\left(\mathbf{r}\right)]R_{11}^2$. Density peaks of species 1 are colored black while the peaks of species 2 are white. The front was initiated at time $t=0$ along the line $x=0$. The top profile corresponds to time $t^{*}=0.6$ while the lower profile corresponds to $t^{*}=3$.

Figure 12

Analysis of the density peaks in the density profile in a GEM-8 mixture with $\beta {\epsilon }_{ij}=\beta \epsilon =1$ for all $i,j=1,2$, $R_{22}/R_{11}=1.5$, and $R_{12}/R_{11}=1$ and average total density $\rho R_{11}^2=8$ and concentration $\varphi =0.5$, formed by a solidification front initiated along the line $x=25$ at time $t=0$. The diagrams along the top row correspond to time $t^{*}=2$, shortly after the solidification front has exited the domain and before the structure has had time to relax, while the diagrams along the bottom row correspond to time $t^{*}=400$, when the profiles no longer change in time—the system has reached a local minimum of the free energy. Left: Voronoi diagrams; the construction reveals the disorder created by the front. The hexagons and squares correspond to the two competing crystal structures. Middle: Delauney triangulation; domains of the hexagonal phase (equilateral triangles) are highlighted in red, while the remainder, including the right-angled triangles of the square phase, are shown in black. Right: The density maxima are color coded according to the triangle type they belong to as follows: right-angled triangles are black, equilateral triangles are gray (red online), and scalene triangles are open circles. Comparing the upper to the lower diagrams, we see that over time there is an increase in the size of the domains of the two different crystal structures.

Figure 13

Time evolution of the bond angle distribution function from Delauney triangulation, corresponding to the results in Fig. 12.

Figure 14

(a) The speeds $c_s$, $c_{hs}$, and $c_h$ defined in the text as a function of $\gamma$ computed from the model system (A8, A9) for $\lambda =1$ and $\lambda =2$. The results for $\lambda =1$ agree with those in Ref. [25]. The full range of values of $\gamma$ is shown including the Maxwell points ${\gamma }_M$, where $c=0$, and the location of the critical values ${\gamma }_1$ and ${\gamma }_2$, where $c_h=c_s$ and $c_{hs}=c_s$, respectively. (b) The location of the Maxwell point and the critical values ${\gamma }_1$ and ${\gamma }_2$ as a function of the nonlinear coupling coefficient $\lambda$. The dotted line shows $-2.5{\gamma }_1$ and indicates that, in the range considered, the ratio ${\gamma }_M/{\gamma }_1$ is nearly constant.