Two-Step Melting in Two Dimensions: First-Order Liquid-Hexatic Transition

Melting in two spatial dimensions, as realized in thin films or at interfaces, represents one of the most fascinating phase transitions in nature, but it remains poorly understood. Even for the fundamental hard-disk model, the melting mechanism has not been agreed upon after 50 years of studies. A recent Monte Carlo algorithm allows us to thermalize systems large enough to access the thermodynamic regime. We show that melting in hard disks proceeds in two steps with a liquid phase, a hexatic phase, and a solid. The hexatic-solid transition is continuous while, surprisingly, the liquid-hexatic transition is of first order. This melting scenario solves one of the fundamental statistical-physics models, which is at the root of a large body of theoretical, computational, and experimental research.

Figure 1

Phase coexistence for ${1024}^2$ thermalized hard disks at density $\eta =0.708$. (a) Color-coded local orientations ${\Psi }_k$ showing long orientational correlations [blue region, see (b), (c)] coexisting with short-range correlations [see (d)]. (e) Local densities (averaged over a radius of $50\sigma$), demonstrating the connection between density and local orientation (see 22). In (b), (c), and (d), disks with five (seven) neighbors are colored in gray (black).

Figure 2

Equilibrium equation of state for hard disks. The pressure is plotted vs volume per particle [$v=V/N$) (lower scale) and density $\eta$ (upper scale)]. In the coexistence region, the strong system-size dependence stems from the interface free energy. The Maxwell constructions (horizontal lines) suppress the interface effects (with a convex free energy) for each $N$. Stripe [(c), for $N={1024}^2$] and bubble configurations (b), (d) are shown in the coexistence region, together with two single-phase configurations (a), (e). The interface free energy per disk $\beta \Delta f$ (hatched area) scales as $1/\sqrt{N}$ (f).

Figure 3

Configuration-averaged two-dimensional pair correlation. $g(\Delta \mathbf{r})$ is obtained from the two-dimensional histogram of periodic distances $\Delta {\mathbf{r}}_{ij}={\mathbf{r}}_i-{\mathbf{r}}_j$. (a) Pair correlation $g(\Delta \mathbf{r})$ at density $\eta =0.718$ for small $\Delta \mathbf{r}=(\Delta x,\Delta y)$. Each disk configuration is oriented with respect to $\Psi$. The excellent contrast between the peak and the bottom values of $g(\Delta \mathbf{r})$ at $|\Delta \mathbf{r}|\gtrsim 2\sigma$, of about $(16:0.2)$, provides evidence for the single-phase nature of the system. (b) Cut of the sample-averaged $g(\Delta \mathbf{r})-1$ for $\Delta \mathbf{r}=(\Delta x,0)$. Decay is exponential for $\eta =0.718$ and algebraic for $\eta =0.720$. (See [22] for positional and orientational correlation functions.)

Figure 4

Approach to thermal equilibrium from different initial conditions. (a), (b): ${1024}^2$ hard disks at density $\eta =0.708$, after a quench from (a) a high-density crystal and from (b) a low-density liquid, showing coarsening leading to phase separation (color code for ${\Psi }_k$ as in Fig. 1b, see also 22). Each of the runs takes about one week of CPU time. (c) Absolute value of the sample orientation for the simulations in (a) and (b), compared to runs with the local Monte Carlo algorithm from the same initial conditions (time in attempted displacements per disk). The correlation time of the event-chain algorithm, on the order of ${10}^6$ displacements per disk, estimated from (c), agrees with the correlation time estimated in our production runs with $6\times {10}^7$ total displacements per disk.