Density affects the nature of the hexaticliquid transition in twodimensional melting of softcore systems
Abstract
We find that both continuous and discontinuous hexaticliquid transitions can happen in the melting of twodimensional solids of softcore disks. For three typical model systems, Hertzian, harmonic, and Gaussiancore models, we observe the same scenarios. These systems exhibit reentrant crystallization (melting) with a maximum melting temperature happening at a crossover density . The hexaticliquid transition at a density smaller than is discontinuous. Liquid and hexatic phases coexist in a density interval, which becomes narrower with increasing temperature and tends to vanish approximately at . Above , the transition is continuous, in agreement with the KosterlitzThoulessHalperinNelsonYoung theory. For these softcore systems, the nature of the hexaticliquid transition depends on density (pressure), with the melting at being a plausible transition point from discontinuous to continuous hexaticliquid transition.
pacs:
64.70.D, 82.70.Dd, 61.20.JaTwodimensional melting is one of the most fascinating and puzzling phase transitions strandburg (); dash (); gasser (). In contrast to the firstorder nature in three dimensions, the possible existence of an intermediate phase between liquid and solid, e.g., the hexatic phase, confuses the nature of twodimensional melting. According to the KosterlitzThoulessHalperinNelsonYoung (KTHNY) theory, the transitions from solid to hexatic and from hexatic to liquid are both continuous, accompanied by the disappearance of quasilongrange positional and orientational orders, respectively kt (); hn1 (); hn2 (); young (). Many experiments and simulations have confirmed the twostage melting proposed by the KTHNY theory murray (); zahn1 (); von (); keim (); lin (); lin1 (); lee (); qi (); prestipino1 (); shiba (), while there are still exceptions strandburg (); alba (); chui (); lansac (). The continuity of the hexaticliquid transition also remains a matter of debate daan (); marcus ().
Recent studies have suggested that the nature of the hexaticliquid transition is sensitive to the details of interparticle potential, including range, softness, length scale, and so on lee (); prestipino (); dudalov (); werner1 (). For instance, it has been confirmed that the hexaticliquid transition of hard disks is first order werner (); engel (); dijkstra (). In contrast, twodimensional melting of ultrasoft Gaussiancore particles was claimed to be consistent with the KTHNY theory prestipino (). By tuning the exponent of the inverse power law interparticle potential and hence the particle softness, Kapfer and Krauth observed the intriguing evolution of the hexaticliquid transition from discontinuous to continuous werner1 ().
Consider a widely studied model system with finite range, purely repulsive and softcore particle interaction
(1) 
where is the separation between particles and , is the particle diameter, is the Heaviside function, is the characteristic energy scale, and is a tunable parameter. At low temperatures and low densities, this system behaves as a hard sphere (disk) system xu (). Its melting temperature increases with density up to the maximum value at a crossover density . Above , the melting temperature instead decreases with increasing density, exhibiting reentrant crystallization (melting) prestipino (); jagla (); alan (); yu (). As shown in Fig. 1 of the phase diagram for Hertzian repulsion () in two dimensions, multiple reentrant crystallizations with different crystal structures occur successively with increasing density. Therefore, both the hard and ultrasoft particle limits can be achieved by the same model, just by varying the density. It is then interesting to know if both continuous and discontinuous hexaticliquid transitions can occur in the same system.
By systematically studying the twodimensional melting of Hertzian and harmonic () systems over a wide range of densities, we indeed observe both types of the hexaticliquid transition. Interestingly, the crossover density may act as the transition point between the two types. When , the transition is discontinuous, showing the coexistence of liquid and hexatic phases. The density region of the coexistence decreases with increasing temperature and tends to vanish at . When , the transition is continuous. We further verify that the same scenario exists for Gaussiancore model. Therefore, we propose that density affects the nature of the hexaticliquid transition for softcore particles exhibiting reentrant crystallization.
Our systems are rectangular boxes containing disks with diameter and mass . The systems have a side length ratio to accommodate the perfect triangular structure. Periodic boundary conditions are applied in both directions. We set the units of mass, energy, and length to be , , and . The time is thus in units of . The temperature is in units of , with being the Boltzmann constant. The density is calculated as .
The liquid, hexatic, and solid phases are identified from correlation functions of the bondorientational and positional order parameters according to the KTHNY theory binder (); daan (); lee (); prestipino (); buldyrev (); vilaseca ():
(2)  
(3) 
where is the separation between particles and located at and respectively, is the wave vector satisfying the periodic boundary conditions and at the first peak of the static structure factor, and denotes the average over configurations and particles. The local bondorientational order parameter for particle is defined as
(4) 
where the sum is over all nearest neighbors of particle determined by the Voronoi tessellation, and is the angle between and a reference direction.
For the liquid phase, both and show exponential decay corresponding to shortrange order. The hexatic phase has quasilongrange bondorientational order and shortrange positional order, resulting in a powerlaw decay of , with , and an exponential decay of . For the solid phase, with and shows almost no decay due to the quasilongrange positional order and longrange bondorientational order. In the Supplemental Material SM (), we show some examples of the correlation functions and also the subblock scaling analysis bagchi () to distinguish different phases.
We first study systems of Hertzian and harmonic repulsions. They have been widely employed in simulation and theoretical work and have been shown to approximate well interactions of various experimental systems such as polyNisopropylacrylamide colloids, granular materials, and foams zhang (); behringer (); eric (). Both repulsions are soft core with positive definite Fourier transform SM (), leading to reentrant crystallization likos (). Upon compression, there occurs a sequence of reentrant crystallizations with different solid structures william (). In this work, we concentrate only on the first one with the triangular structure.
Figure 1 is obtained by quenching hightemperature states with a slow rate using constanttemperature and constantpressure molecular dynamics simulations note1 (). We have verified that our quench rate is slow enough that even slower quench rates will not change the phase diagram significantly. The phase diagram shows approximate locations of the phase boundaries, which slightly vary with system size due to finite size effects. The maximum melting temperature for Hertzian (harmonic) repulsion estimated from the phase diagram is approximately () at a crossover density () or pressure () SM ().
The inset to Fig. 1 shows the isobaric equation of state across the phase boundaries on both sides of and approximately at . When , the density jumps up across the transitions from liquid to solid. When , the system exhibits a waterlike anomaly with the density of solid being lower than that of liquid. We find that the absolute value of the fast density change decreases when approaching from either side. The melting at may behave as a turning point with daan1 (). As shown in the inset to Fig. 1, there is almost no sign of a density discontinuity when note2 ().
The melting at looks special at least for the continuity in density. It is interesting to figure out what role it plays in the twodimensional melting of softcore systems. To probe the details of the melting, we simulate much larger systems up to using parallel LAMMPS package lammps () in an or ensemble and on both sides of .
We calculate the equilibrium isothermal equation of state in the ensemble across the transitions from solid to liquid. Figure 2(a) shows for Hertzian disks calculated at and . The curve displays a MayerWood loop mayer (), characterizing phase coexistence. The loop is due to interface free energy between coexistent phases in finite size systems furukawa (); schrader (). We fit the curve with a order polynomial, and determine the boundaries of coexistence by the Maxwell construction. Seen from Fig. 2(a), it is the coexistence of hexatic and liquid phases, because these two phases exist on both sides of the coexistence.
The interface free energy per particle is calculated as half of the area encircled by the polynomial curve and the horizontal line of the Maxwell construction. With increasing system size, the MayerWood loop flattens, so tends to decrease with increasing . Figure 2(b) shows that , further demonstrating the discontinuous nature of the hexaticliquid transition at werner (); jlee ().
Moreover, we find that the density interval of the phase coexistence decreases with increasing temperature approaching from the side. As shown in Fig. 2(c), can be fitted well with a powerlaw scaling relation: , where and are interaction dependent fitting parameters. The value of used in Fig. 2(c) is () for Hertzian (harmonic) repulsion, in good agreement with estimated from the phase diagram. It is thus plausible to conjecture that the hexaticliquid transition at becomes continuous.
What may happen for melting at ? In Fig. 2(d), we show at the same temperature as for Fig. 2(a), but on the higher density side of . Across the transitions, monotonically increases with note5 (). Therefore, the hexaticliquid transition is continuous and agrees with the KTHNY theory. We have also verified that the same phenomenon occurs at all other temperatures.
In Fig. 3, we further compare the system size dependence of the isobaric density , enthalpy and average bondorientational order calculated in the ensemble on both sides of note4 (), where denotes average over particles and configurations. When , all quantities apparently tend to be discontinuous with increasing system size, while they do not show such a tendency when .
Figures 2 and 3 provide robust evidence to suggest that the hexaticliquid transition undergoes a transition from discontinuous to continuous, with the melting at being a possible transition point. In Section IV of the Supplemental Material SM (), we provide another evidence by showing that the correlation length in the liquid phase tends to diverge approaching the maximum melting temperature from the side. Two different types of hexaticliquid transition can be achieved in the same system, just by tuning the density. Now there comes the question of whether the scenario is specific to systems described by Eq. (1) or exists in other softcore systems. Next, we will examine the widely studied Gaussiancore model and show that our observations are not unique to Hertzian and harmonic repulsions.
The potential between interacting particles and for the Gaussiancore model is , with all parameters having the same meanings as for Eq. (1). We set a potential cutoff at and shift the potential to make sure that both the potential and force vanish at . We also use the same set of units as for Hertzian and harmonic systems. The Gaussiancore model exhibits reentrant crystallization with maximum melting temperature happening at and estimated from the phase diagram of systems SM ().
Figure 4 compares isothermal for Gaussiancore model calculated in the ensemble on both sides of and at . Like Hertzian and harmonic repulsions, Fig. 4(a) shows that at has a clear MayerWood loop, so the hexaticliquid transition here is discontinuous. The inset to Fig. 4(a) shows that the coexistent region also decreases with increasing temperature and can be well fitted with , where agrees well with estimated from the phase diagram. Again, for Gaussiancore model, melting at is likely to become continuous. In contrast, the continuity of the transitions above is robust. The curve at shown in Fig. 4(b) is rather straight across the melting with an almost density independent compressibility.
By studying three representative softcore models exhibiting reentrant crystallization, we find that both continuous and discontinuous hexaticliquid transitions happen in the same system. The type of the transition is determined by density. Our data suggest that the melting point at the maximum melting temperature may be the demarcation between the two types of transitions. Note that Hertzian and harmonic models are quite different from Gaussiancore model SM (), but they still behave similarly in the hexaticliquid transition. Although it is impossible to check all models, based on our study, we are inclined to believe that our observations generalize to softcore systems with reentrant crystallization. Anyhow, our study reveals the unknown extraordinary features of twodimensional melting of softcore systems, which can be tested in experimental systems such as star polymers watz ().
In addition to the hexatic phase, the existence of the analogous tetratic phase upon the melting of solids with square lattice structure has been reported and discussed takamichi (); daan2 (); peng (). However, compared to the hexatic phase, the tetratic phase is much less studied. One possible reason is that the square lattice structure is more difficult to form than the triangular lattice. Hertzian and harmonic models exhibit multiple reentrant crystallizations with various solid structures, which are ideal to investigate the tetratic phase and other intermediate phases. It would be interesting to know next if we are able to observe different intermediate phases in these simple model systems and if the melting of various types of solids follows similar scenarios or not.
We are grateful to Werner Krauth and Peng Tan for helpful discussions. This work is supported by National Natural Science Foundation of China No. 21325418 and 11574278, National Basic Research Program of China (973 Program) No. 2012CB821500, and Fundamental Research Funds for the Central Universities No. 2030020028. We also thank the Supercomputing Center of University of Science and Technology of China for computer times.
References
 (1)

(2)
[]ningxu@ustc.edu.cn
 (3) K. L. Strandburg, Rev. Mod. Phys. 60, 161 (1988).
 (4) J. G. Dash, Rev. Mod. Phys. 71, 1737 (1999).
 (5) U. Gasser, J. Phys.: Condens. Matter 21, 203101 (2009).
 (6) J. M. Kosterlitz and D. J. Thouless, J. Phys. C: Solid State Phys. 6, 1181 (1973).
 (7) D. R. Nelson and B. I. Halperin, Phys. Rev. B 19, 2457 (1979).
 (8) B. I. Halperin and D. R. Nelson, Phys. Rev. Lett. 41, 121 (1978).
 (9) A. P. Young, Phys. Rev. B 19, 1855 (1979).
 (10) C. A. Murray and D. H. Van Winkle, Phys. Rev. Lett. 58, 1200 (1987).
 (11) K. Zahn, R. Lenke, and G. Maret, Phys. Rev. Lett. 82, 2721 (1999).
 (12) H. H. von Grünberg, P. Keim, K. Zahn, and G. Maret, Phys. Rev. Lett. 93, 255703 (2004).
 (13) P. Keim, G. Maret, and H. H. von Grünberg, Phys. Rev. E 75, 031402 (2007).
 (14) B. J. Lin and L. J. Chen, J. Chem. Phys. 126, 034706 (2007).
 (15) S. Z. Lin, B. Zheng, and S. Trimper, Phys. Rev. E 73, 066106 (2006).
 (16) S. I. Lee and S. J. Lee, Phys. Rev. E 78, 041504 (2008).
 (17) W. K. Qi, Z. Wang, Y. Han, and Y. Chen, J. Chem. Phys. 133, 234508 (2010).
 (18) S. Prestipino, F. Saija, and P. V. Giaquinta, J. Chem. Phys. 137, 104503 (2012).
 (19) H. Shiba, A. Onuki, and T. Araki, Europhys. Lett. 86, 66004 (2009).
 (20) C. AlbaSimionesco, B. Coasne, G. Dosseh, G. Dudziak, K. E. Gubbins, R. Radhakrishnan, and M. SliwinskaBartkowiak, J. Phys.: Condens. Matter. 18, R15 (2006).
 (21) S. T. Chui, Phys. Rev. Lett. 48, 933 (1982).
 (22) Y. Lansac, M. A. Glaser, and N. A. Clark, Phys. Rev. E 73, 041501 (2006).
 (23) P. Bladon and D. Frenkel, Phys. Rev. Lett. 74, 2519 (1995).
 (24) A. H. Marcus and S. A. Rice, Phys. Rev. Lett. 77, 2577 (1996).
 (25) S. Prestipino, F. Saija, and P. V. Giaquinta, Phys. Rev. Lett. 106, 235701 (2011).
 (26) D. E. Dudalov, Yu. D. Fomin, E. N. Tsiok, and V. N. Ryzhov, J. Phys.: Conf. Series 510, 012016 (2014).
 (27) S. C. Kapfer and W. Krauth, Phys. Rev. Lett. 114, 035702 (2015).
 (28) E. P. Bernard and W. Krauth, Phys. Rev. Lett. 107, 155704 (2011).
 (29) M. Engel, J. A. Anderson, S. C. Glotzer, M. Isobe, E. P. Bernard, and W. Krauth, Phys. Rev. E 87, 042134 (2013).
 (30) W. Qi, A. P. Gantapara, and M. Dijkstra, Soft Matter 10, 5449 (2014).
 (31) N. Xu, T. K. Haxton, A. J. Liu, and S. R. Nagel, Phys. Rev. Lett. 103, 245701 (2009); X. Wang, W. Zheng, L. Wang, and N. Xu, Phys. Rev. Lett. 114, 035502 (2015).
 (32) E. A. Jagla, J. Chem. Phys. 111, 8980 (1999).
 (33) A. B. de Oliveira, P. A. Netz, T. Colla, and M. C. Barbosa, J. Chem. Phys. 124, 084505 (2006).
 (34) Yu. D. Fomin, N. V. Gribova, V. N. Ryzhov, S. M. Stishov, and D. Frenkel, J. Chem. Phys. 129, 064512 (2008).
 (35) K. Binder, S. Sengupta, and P. Nielaba, J. Phys.: Condens. Matter 14, 2323 (2002).
 (36) S. V. Buldyrev, G. Malescio, C. A. Angell, N. Giovambattista, S. Prestipino, F. Sajia, H. E. Stanley, and L. Xu, J. Phys.: Condens. Matter 21, 504106 (2009).
 (37) P. Vilaseca and G. Franzese, J. NonCrystal. Solids 357, 419 (2011).
 (38) See Supplemental Material at http://link.aps.org/ supplemental/10.1103/PhysRevLett.117.085702 for additional information about the criterion of reentrant crystallization, some results for harmonic and Gaussiancore models, more supporting information for Hertzian model, and the identification and visualization of different phases.
 (39) K. Bagchi, H. C. Andersen, and W. Swope, Phys. Rev. Lett. 76, 255 (1996).
 (40) Z. Zhang, N. Xu, D. T. N. Chen, P. Yunker, A. M. Alsayed, K. B. Aptowicz, P. Habdas, A. J. Liu, S. R. Nagel, and A. G. Yodh, Nature (London) 459, 230 (2009).
 (41) T. S. Majmudar, M. Sperl, S. Luding, and R. P. Behringer, Phys. Rev. Lett. 98, 058001 (2007).
 (42) K. W. Desmond, P. J. Young, D. Chen, and E. R. Weeks, Soft Matter 9, 3424 (2013).
 (43) C. N. Likos, A. Lang, M. Watzlawek, and H. Löwen, Phys. Rev. E 63, 031206 (2001).
 (44) W. L. Miller and A. Cacciuto, Soft Matter 7, 7552 (2011).
 (45) Our phase diagram shows exactly the first two regions below in Fig. 2 of Ref. william (). Due to low data resolution, the maximum of the melting curve for square solid structure does not show up in Fig. 2 of Ref. william ().
 (46) P. Bolhuis, M. Hagen, and D. Frenkel, Phys. Rev. E 50, 4880 (1994).
 (47) As discussed in Section III of the Supplemental Material SM (), even though at , there is still a thin layer of hexatic phase between solid and liquid.
 (48) http://lammps.sandia.gov/.
 (49) J. E. Mayer and W. W. Wood, J. Chem. Phys. 42, 4268 (1965).
 (50) H. Furukawa and K. Binder, Phys. Rev. A 26, 556 (1982).
 (51) M. Schrader, P. Virnau, and K. Binder, Phys. Rev. E 79, 061104 (2009).
 (52) J. Lee and J. M. Kosterlitz, Phys. Rev. Lett. 65, 137 (1990).
 (53) As shown in Fig. S5(b) of the Supplemental Material SM (), the system size effects are already small for the largest systems studied, so it is plausible to expect that the MayerWood loop and phase coexistence as observed at do not occur.
 (54) Refer to Fig. S3 of the Supplemental Material SM () for a phase diagram in the plane.
 (55) M. Watzlawek, C. N. Likos, and H. Löwen, Phys. Rev. Lett. 82, 5289 (1999).
 (56) T. Terao, J. Chem. Phys. 139, 134501 (2013).
 (57) K. W. Wojciechowski and D. Frenkel, Comp. Met. Sci. Technol. 10, 235 (2004).
 (58) Y. Peng, Z. Wang, A. M. Alsayed, A. G. Yodh, and Y. Han, Phys. Rev. Lett. 104, 205703 (2010).
Appendix A Supplemental Material
a.1 I. Criterion of reentrant crystallization (melting) and phase diagrams
