Two-dimensional melting in simple atomic systems : continuous vs. discontinous melting

Here, we first deal with the case of a moderately soft (and longer ranged) potential with $\alpha=3.5$. Figure 2 shows the isothermal equation of state (at $T=0.7$) on the plane of pressure vs. density. This was obtained by the NH-MD simulations by decoupling the PR (isobaric) part from NH-PR MD equations of motions by taking $M_{v}=\infty$ which reduces to $NVT$ condition. We here define the density $\rho$ as $\rho\equiv N\sigma^{2}/V$ where $N$ is the total number of particles and $V$ the total volume (area in two dimensions) of the system. Densities were chosen from the range $\rho=1.56\sim 1.62$, with the density increment of $\Delta\rho=0.005$ and the pressure was evaluated by means of the virial expression (with $k_{B}=1$). This range of the density corresponds to the region of transition from liquid to solid. For each density, $10^{6}\sim 3\times 10^{6}$ steps of integration were carried out for equilibration beginning with a configuration of triangular lattice and, after equilibration, $10^{7}$ steps of integration were performed for thermodynamic calculations.

The number of particles employed was $N=3600$. In order to reduce the finite size boundary effect, we used a rhombic box (with the smaller side angle of $60$ degrees) for the shape of the system with periodic boundary conditions. However, independent results of ours from square box showed no significant difference (from those of rhombic box) with respect to the quantities of our interest.

The isothermal curve increases monotonically near the transition region satisfying the condition of mechanical stability (unlike the discontinuity of density in a first order transition) that the isothermal compressibility should be positive $K_{T}=(1/\rho)(\partial\rho/\partial P)_{T}>0$. We may identify the boundary of stable hexatic phase as the values of the density where an abrupt change in the isothermal compressibility occurs. In this way, we estimate the density of solid-hexatic transition as $\rho_{s-h}\simeq 1.6$.

Although the change in isothermal compressibility is less conspicuous near the hexatic-liquid boundary, we see that, near the density $1.58\leq\rho\leq 1.585$, there exists a crossover in the slope of the isothermal compressibility. Below, we give an estimation of the density of hexatic-liquid transition by applying a theoretical expectation from KTHNY theory on the decay exponent of the spatial correlation of the orientational order parameter (see below).

The fact that the pressure within the hexatic phase is monotonically increasing as the density increases (with the resulting isothermal compressibility kept positive) appears to be a very compelling evidence for a stable hexatic phase in thermal equilibrium.

In order to distinguish the orientational order of the phases, we have computed the bond-orientational correlation function $G_{6}(r)$ defined above.

Figure. 3 is the bond-orientational correlation function for the range of the density ($1.56\leq\rho\leq 1.615$). We find that, for the density range of $1.585<\rho\leq 1.595$, the averaged correlation functions exhibit algebraic decays with the decay exponent $\eta<1/4$ while, at lower densities, they exhibit decays with faster than power law behavior in the long distance limit. At $\rho=1.585$, the orientational correlation exhibit a slope of approximately $\eta=0.25$. Here, the crossover from a power law decay to exponential takes places with exponent approximately equal to $1/4$, and also this value of the crossover density agrees almost precisely with the density region exhibiting an abrupt change of the slope i.e., of the compressiblity in the equation of state (Figure 2).

Now, the fourth order Binder cumulant for the global bond-orientationl order parameter for sub-block systems of (linear) size $L$ is shown in Fig. 4. We can see that the density at the crossing point is around $\rho\simeq 1.585$. This is compatible with the boundary density ($1.585$) between the liquid and the hexatic phase which was shown above in isothermal curve for $NVT$ ensemble, and also compatible with the the boundary density between liquid and hexatic-like phase in terms of the power law decay of the bond-orientational order near $\rho=1.585$

It may be expected theoretically that the binder cumulants of local orientational order in hexatic phase collapse to a line because of the critical charateristic of the phase. However, we may not consider non-collapse of the Binder cumulants to a line as an evidence for non-existence of hexatic phases. This is because, even for the case of XY model, complete collapse was not found in the region where orientational order decay algebraically, but rather it exhibits a crossing point at the transition temperature oliveira .

Now, we turn to the linear susceptibility for the global bond-orientational order parameters near the melting transition. Shown in Fig. 5 is the sub-block susceptibility obtained from the fluctuation of the bond-orientational order parameter for sub-blocks of the system with linear sizes $L$ with $L=N/M_{b}$ where $M_{b}$ ranges from $11$ to $20$. We see that the suceptibility shows a broad peak region near the density $1.580<\rho<1.585$ which borders the liquid-hexatic phase boundary region. We also see that the suceptibility exhibits broader shape (showing weaker dependence on the density) in the liquid region compared with other cases that will be shown below.

Figure 6 shows a snapshot of the particle configuration for density $\rho=1.585$ within the hexatic phase region but close to the transition (to liquid phase) at the temperature $T=0.7$, which shows free dislocations (i.e., bounded pair disclinations). This shows rather clearly the fundamental role of defects leading to the power law decay of orientational correlations.

In order to further understand the nature of the hexatic phase we obtained the histogram distributionLee90 of the density order parameter for five different system sizes ($N=900,1600,3600,10000$) under constant external pressure and temperature of $T=0.7$, and $P=13.5$, where a hexatic phase is expected to occur from our measurement of the orientational correlations. In Fig. 7 we see that all the histograms exhibit single peaks, which implies that there exist a unique phase with minimum free energy. It is also observed that, as the number of particles increases, the peak height becomes higher and, at the same time, the width of the peak decreases. Also the position of the peak tends to shift to the lower density (within the hexatic regime) probably due to the development of long range fluctuation. This indicates that this region corresponds to a single phase region (unlike solid-liquid mixture) consistent with the absence of van der Waals loop in the pressure vs. density curve.

All of these observations lead us to the conclusion that there exist a thermodynamcally stable hexatic phase consistent with the KTHNY melting scenario for the case of $\alpha=3.5$.

Next, we examine an intermediate-range potential with $\alpha=6$, which corresponds approximately to the famous LJ potentialtclim2003 . In Fig. 8, the equation of state exhibits a weak van der Waals-like loop in the pressure vs. density, which indicates a first order transition. The unstable region ranges from $\rho\simeq 1.04$ to $\rho\simeq 1.065$. This is confirmed more rigorously by the histogram distributions of the density in $NPT$ ensemble simulations for different system sizes. In Fig. 9, the histogram distributions of the density for systems with $N=900$, $3600$, and $10000$ are shown from $NVT$ ensemble for $T=0.57$ and $P=1.85$. For the cases of $N=900$ and $3600$, we observe transitions between two peaks (through the valley of finite height between the peaks) with the resulting double peaks in the histograms. For the $N=10000$ systems, however, we can no longer observe crossing between the the coexisting phases. Instead, we could observe two different (separate) histograms that are determined by the initial states, depending whether the initial state is in the ordered solid phase or in the disordered liquid phase. Evidently, the free energy barrier increases with increasing system size, which indicates clearly that the transition is of first orderLee90 . Nevertheless, the system configurations in the coexistence region resembles those of the hexatic phases (Fig. 10), showing algebraically decaying orientational order (Fig 11). We see that the boundary between liquid and hexatic-like phase in terms of the orientational order is located around $\rho=1.05\sim 1.055$.

Furthermore, we also observe a hexatic-like feature in the $NPT$ ensemble, where the system goes through temporarily a hexatic-like phase before transiting into the other phase. This kind of characteristics in $NPT$ ensemble seems to have been already reported as ‘metastable hexatic phase’ in LJ system by Chen et alChen95 ; Somer97 ; Somer98 . Although they argued that large system size ($N\simeq 40000$) is necessary to observe this kind of features, we could observe such metastable hexatic phases even for systems with smaller sizes of $N=1600$. We think that this hexatic-like feature in a first order transition can be attributed to weakly first order nature of the transition. As shown below, the relative free energy barrier in this case is considerably lower than that in the case of shorter-ranged potential of $\alpha=12$, from which we presume that the metastable or unstable hexatic-like phase comes from defect proliferation by thermal fuctuations, but not from some mechanism leading to a true phase transition.

Now, the fourth order Binder cumulant for the global bond-orientationl order parameter for sub-block systems of (linear) size $L$ is shown in Fig. 12. We can see that the density at the crossing point is located near $\rho\sim 1.06$. This corresponds to a point in the middle of the unstable part of the van der Waals curve in the isothermal equation of state.

Next, we deal with the linear susceptibility for the bond-orientational order parameters. Shown in Fig. 13 is the sub-block susceptibility obtained from the fluctuation of the bond-orientational order parameter for sub-blocks of the system with linear sizes $L$ with $L=N/M_{b}$ where $M_{b}$ ranges from $10$ to $20$. We see that the suceptibility shows a sharper peak (as compared with the case of $\alpha=3.5$) at the density $\rho\simeq 1.045$ which is a little bit below the density ($\rho=1.05\sim 1.055$) where the orientational correlation exhibits a spatial decay exponent of $0.25$. This might be interpreted as a small evidence that the nature of the melting transition of this system is inconsistent with the expectation of the KTHNY theory.

Finally, we investigate the case of a short-ranged potential with $\alpha=12$. Figure. 14 shows the equation of state in the density region of $0.85\leq\rho\leq 1.075$ obtained from $NVT$ ensemble simulations at $T=0.57$. We can see that the equation of state exhibits a van der Waals-like region in the density with the unstable region ranging from $\rho\simeq 0.96$ to $\rho\simeq 1.04$ clearly indicating a first order melting transition.

The fourth order Binder cumulant for the global bond-orientationl order parameter for sub-block system size $L$ is shown in Fig. 15. We can see that the density at the crossing point is located near $\rho\sim 1.045$ which is located near the lower density limit of the metastable solid (spinodal) as shown in the $NVT$ isothermal equation of state.

First order nature of the melting transition is confirmed further with the double peak nature of the histogram distributions of the density from $NPT$ ensemble simulations for system sizes of $N=400,900,3600,10000$ as shown in Fig. 16. In this case, the free energy barrier is presumably much higher than the case of $\alpha=6$, and none of the systems with $\alpha=12$ exhibit any tunneling transitions between liquid-like states to solid-like states during $10^{8}$ MD steps of simulations with $\Delta T=0.002$. Therefore, double peak histogram distributions for each of the system sizes are actually obtained by combining two separate histograms, one with ordered initial states (higher density) and the other with disordered initial states (lower density), respectively.

Also, typical system configuration for $\alpha=12$ is shown in Fig. 17 for $\rho=1.0$ and $T=0.57$ corresponding to the coexisting region, where we can see that hexatic-like feature disappears, and that liquid phase region consisting of defect clusters coexists with solid region. Also in $NPT$ ensemble simulations of the melting process, we can no longer observe metastable or unstable hexatic phase, but observe a discontinuous abrupt change in density. From these observations we thus conclude that first order nature of the melting gets stronger as the potential range decreases.

Next, we look into the linear susceptibility for the bond-orientational order parameters. Shown in Fig. 18 is the sub-block susceptibility obtained from the fluctuation of the bond-orientational order parameter for sub-blocks of the system with linear sizes $L$ with $L=N/M_{b}$ where $M_{b}$ ranges from $12$ to $20$. We see that the suceptibility shows a peak at the density $\rho\simeq 0.96$ which is located near the limit of the metastable liquid (spinodal) as shown in the curve of the isothermal equation of state from $NVT$ ensemble.