Soft-core particles freezing to form a quasicrystal and a crystal-liquid phase

Systems of soft-core particles interacting via a two-scale potential are studied. The potential is responsible for peaks in the structure factor of the liquid state at two different but comparable length scales and a similar bimodal structure is evident in the dispersion relation. Dynamical density functional theory in two dimensions is used to identify two unusual states of this system: a crystal-liquid state, in which the majority of the particles are located on lattice sites but a minority remains free and so behaves like a liquid, and a 12-fold quasicrystalline state. Both are present even for deeply quenched liquids and are found in a regime in which the liquid is unstable with respect to modulations on the smaller scale only. As a result, the system initially evolves towards a small-scale crystal state; this state is not a minimum of the free energy, however, and so the system subsequently attempts to reorganize to generate the lower-energy larger-scale crystals. This dynamical process generates a disordered state with quasicrystalline domains and takes place even when this large scale is linearly stable, i.e., it is a nonlinear process. With controlled initial conditions, a perfect quasicrystal can form. The results are corroborated using Brownian dynamics simulations.

Figure 1

Phase diagram in the $(a,{\rho }_0{R}^2)$ plane, where ${\rho }_0$ is the average density. The system exhibits three phases: a liquid and two crystalline phases. Crystal A is hexagonal, with a large lattice spacing, while crystal B, which is also hexagonal, has a smaller lattice spacing. The short-dashed blue line is the linear instability threshold (spinodal) for the uniform liquid. The long-dashed green horizontal line is where the two principal peaks in the static structure factor $S\left(k\right)$ have the same height, while the pink dotted line terminating in a circle is the locus where the two peaks in the dispersion relation $\omega \left(k\right)$ have equal height. The circle marks the point where the smaller $k$ peak disappears.

Figure 2

The top panel shows the density profile in the $(x/R,y/R)$ plane near an interface between coexisting crystal A and crystal B phases, for $a=0.75$. Because the lattice spacings of the two crystal structures are not commensurate, defects along the interface are necessarily present. The bottom panel shows the density profile between the uniform liquid and the [1,1] interface of the crystal B phase, for $a=0.5$.

Figure 3

Correlation functions characterizing the liquid phase for increasing density ${\rho }_0$ as indicated in the key, for fixed $a=0.8$. The top panel shows the radial distribution function $g\left(r\right)$, the middle panel the static structure factor $S\left(k\right)$, and the bottom panel the dispersion relation $\omega \left(k\right)$.

Figure 4

Correlation functions characterizing the liquid phase for a range of values of $a$ as indicated in the key, for fixed ${\rho }_0{R}^2=1.2$. The top panel shows the radial distribution function $g\left(r\right)$, the middle panel the static structure factor $S\left(k\right)$, and bottom panel the dispersion relation $\omega \left(k\right)$.

Figure 5

Dispersion relation $\omega \left(k\right)$ for fixed ${\rho }_0{R}^2=3.5$ and various values of $a$, as indicated in the key.

Figure 6

Dispersion relation $\omega \left(k\right)$ at a series of points along a diagonal path in the phase diagram passing through the point at ${\rho }_0{R}^2=2.95$ and $a=1.067$ at which the two modes $k_{1}R=3.12$ and $k_{2}R=6.03$ are simultaneously marginally unstable. Note that $k_2/k_1=1.932$.

Figure 7

Density profile in the top panel of Fig. 2 displayed in terms of the logarithm of the density $ln\left[\rho \left(\mathbf{r}\right)R^2\right]$ plotted in the $(x/R,y/R)$ plane. This representation allows one to see the fine structure of the density profile away from the principal peaks. Note in particular the honeycomb structure surrounding each of the peaks in the crystal A phase on the left of the interface.

Figure 8

The top panels show $ln\left[\rho \left(\mathbf{r}\right)R^2\right]$ in the $(x/R,y/R)$ plane for a system of $N=600$ particles with $(a,{\rho }_0{R}^2)=(0.8,6)$ confined in a square region of side $L=10R$ obtained from BD simulations (top left) and DFT (top right). The system forms crystal A with a density profile consisting of an array of peaks surrounded by a connected network within which the particles are free to move; this is the crystal-liquid state. The bottom panel is a snapshot from the BD simulation where each particle coordinate is plotted as an open circle.

Figure 9

(a)–(d) Plots of $ln\left[\rho \left(\mathbf{r}\right)R^2\right]$ in the crystal A phase in the $(x/R,y/R)$ plane for fixed $\beta \mu =39$ at the state points: (a) $(a,{\rho }_0{R}^2)=(0.75,4.1)$, (b) $(0.9,3.8)$, (c) $(1.05,3.5)$, and (d) $(1.3,3.1)$. The bottom figure shows a plot of the fraction of mobile particles that are in the liquid part of the density surrounding the density peaks. For $a<0.75$ crystal A is no longer the thermodynamic equilibrium crystal structure and is replaced by crystal B.

Figure 10

Snapshots of $ln\left[\rho \left(\mathbf{r}\right)R^2\right]$ in the $(x/R,y/R)$ plane obtained via Picard iteration for $a=0.8$ and ${\rho }_0{R}^2=3.5$, revealing the evolution towards the equilibrium state for the same state point as the results displayed in the upper panel of Fig. 13. The dispersion relation at this state point is displayed in Fig. 12. The panels along the top row, from left to right, correspond to times $t=30$, 32, and 35 and along the bottom row to $t=40$, 50, and 200. Note that the system first forms the small length scale crystal (at time $t\approx 30)$. It then tries to form the longer length scale crystal. However, due to the small length scale already imprinted on the system, it cannot form a perfect large length scale crystal and ends up forming a disordered system with domains of QC ordering. The Picard iteration used to generate these figures does not locally conserve particle number (although it does conserve the total density in the system; see Sec. 2), but is much faster than the full DDFT and gives qualitatively similar results; compare this figure with Fig. 11, which is calculated with DDFT.

Figure 11

Snapshots of $ln\left[\rho \left(\mathbf{r}\right)R^2\right]$ in the $(x/R,y/R)$ plane obtained from DDFT, for $a=1.067$ and ${\rho }_0{R}^2=3.5$. The dispersion relation at this state point is displayed in Fig. 12. The panels along the top row, from left to right, correspond to times $t/{\tau }_B\equiv t^{*}=1$, 2, and 5 and along the bottom row to $t^{*}=10$, 20, and 40, where ${\tau }_B\equiv \beta R^2/\mathrm{\Gamma }$ is the Brownian time scale. Note that the system first forms a small length scale crystal $\left(t^{*}=2\right)$. It then tries to form the longer length scale crystal, initiated from a grain boundary; see panels for $t^{*}=5$ and 10. However, because of the small length scale already imprinted on it, the system cannot form a perfect long length scale crystal and ends up forming a disordered system with domains of QC ordering; see the final stationary profile at $t^{*}=40$.

Figure 12

Dispersion relation at the state point $a=0.8$ and ${\rho }_0{R}^2=3.5$ (top, corresponding to Fig. 10) and $a=1.067$ and ${\rho }_0{R}^2=3.5$ (bottom, corresponding to Fig. 11). In both cases, QCs form at these state points. In the upper panel (Fig. 10) only one mode is unstable, corresponding to the smaller length scale crystal B. In the lower panel (Fig. 11) two modes are unstable, but the growth rate for the smaller length scale crystal B is much larger than that for the larger length scale crystal A.

Figure 13

The left panels show plots of $ln\left[\rho \left(\mathbf{r}\right)R^2\right]$ in the $(x/R,y/R)$ plane obtained from DFT for $(a,{\rho }_0{R}^2)=(0.8,3.5)$. The right panels are the corresponding Fourier transforms. The latter exhibit 12-fold symmetry, which is indicative of QC ordering. The top density profile is obtained from random initial conditions, while the lower profile was formed starting from an initial density profile having QC symmetry.

Figure 14

Grand potential density as a function of $a$ for fixed $\beta \mu =39$ for the two different crystal structures and also the QC solution displayed in Fig. 13. Near $a=0.75$ there is a point where all three have almost the same value of the grand potential, but the QC solution is never the global minimum (see the inset). The crystal A phase is of crystal-liquid type throughout the range of $a$ shown.

Figure 15

Plots of $ln\left[\rho \left(\mathbf{r}\right)R^2\right]$ in the $(x/R,y/R)$ plane for density profiles obtained when starting from the perfect QC displayed on the bottom left of Fig. 13, which has $R_s/R=1.855$, and then following the solution as $R_s$ is varied, for fixed $\mu$ and domain size. When $R_s$ is decreased, the QC remains stable until $R_s/R=1.77$, at which point the QC profile becomes linearly unstable and the Picard iteration then falls on to the crystal A profile (which contains defects), displayed above left. When $R_s$ is instead increased, the perfect QC solutions become unstable at $R_s/R=2.03$, where the iteration switches to a different branch of QC solutions (above middle), before this new QC state becomes unstable at $R_s/R=2.19$ and evolves to the crystal B branch of solutions displayed above right.