key: cord-0820305-bgmsccxl
authors: Li, Yan-Wei; Ciamarra, Massimo Pica
title: Attraction tames two-dimensional melting: from continuous to discontinuous transitions
date: 2020-05-28
journal: Physical review letters
DOI: 10.1103/physrevlett.124.218002
sha: b11d9f57e0f4b9606460f2499b8e5fb95db1986e
doc_id: 820305
cord_uid: bgmsccxl

Two-dimensional systems may admit a hexatic phase and hexatic-liquid transitions of different natures. The determination of their phase diagrams proved challenging, and indeed those of hard-disks, hard regular polygons, and inverse power-law potentials, have been only recently clarified. In this context, the role of attractive forces is currently speculative, despite their prevalence at both the molecular and colloidal scale. Here we demonstrate, via numerical simulations, that attraction promotes a discontinuous melting scenario with no hexatic phase. At high-temperature, Lennard-Jones particles and attractive polygons follow the shape-dominated melting scenario observed in hard-disks and hard polygons, respectively. Conversely, all systems melt via a first-order transition with no hexatic phase at low temperature, where attractive forces dominate. The intermediate temperature melting scenario is shape-dependent. Our results suggest that, in colloidal experiments, the tunability of the strength of the attractive forces allows for the observation of different melting scenario in the same system.

Two-dimensional (2D) systems with short-range interactions melt either via a first-order solid/liquid transformation or via a two-step process with subsequent solid/hexatic and hexatic/liquid transitions. The twostep scenario may further follow the Kosterlitz-Thouless-Halperin-Nelson-Young (KTNHY) paradigm [1] [2] [3] , with continuous solid-hexatic and hexatic-liquid transitions, or the mixed one [4] , where a discontinuous hexatic-liquid transition follows a continuous solid-hexatic one. The possibly enormous value of the hexatic correlation length makes it difficult to ascertain which of the above melting scenarios a system follows. However, the increase in computational power and the development of novel algorithms, and careful experiments, allowed to make progresses in recent years. For instance, it is now ascertained [4, 5] that, in hard disks, a discontinuous hexaticliquid transition follows a continuous solid-hexatic one. In hard regular polygons, the melting transition depends on the number of edges, e.g., hexagons and squares following the KTHNY melting scenario and pentagons the first-order one [6] . The melting scenario of 2D systems interacting via power-law potentials [7, 8] has been demonstrated to depend on the stiffness of the interaction [9] . In these recently settled cases, density drives the melting transition, and temperature plays no role as the interaction potentials lack an energy scale.

At both the molecular and colloidal scale, attractive forces are prevalent, and the phase behavior is both temperature and density-dependent. The effect of attractive forces on 2D melting remains, however, controversial. Indeed, Nelson noticed that attraction may lead to a variety of phase diagrams, illustrating possible scenarios with the hexatic phase occurring in an intermediate temperature range [10, 11] , for Lennard-Jones (LJ) particles. This would imply that a weak attraction promotes the hexatic phase, while a strong one suppresses it. In attractive systems the existence of the hexatic phase is controversial, as this phase has been observed in some studies [12, 13] , but not in others [14, 15] . The complete mapping of the phase diagram of LJ particles is a recent, but still debated, achievement [16] ; indeed, at high-temperature, where attractive forces are negligible, LJ particles have been found not to follow the melting scenario of 1/r 12 [9, 16] ones.

Here we demonstrate, via the numerical determination of the temperature-density phase diagram of attractive hexagons, pentagons, squares, and LJ point particles, that attraction universally influences the melting scenario by suppressing the hexatic phase and promoting discontinuous transitions.

We simulate attractive hexagons (N = 48071), pentagons (N = 20449) and squares (N = 20521), as well as Lennard-Jones point particles (N = 318 2 ), under periodic boundary conditions, in the canonical ensemble using the GPU-accelerated GALAMOST package [17] . We construct the extended polygons by lumping together LJ point particles equally spaced along the perimeter, as shown in Fig. S1 [18] . The resulting short-ranged attractive interaction, detailed in the Supplemental Material [18] , allows estimating the size of the polygonal particles and their interaction energy scale, we adopt as our units of length and energy, respectively. For the considered state points and interaction, the values of N we consider are large enough for finite-size effects to be negligible, as we prove in Fig. S3 in Supplemental Material [18] . We verify thermal equilibration by ascertaining that the same final state is reached in simulations starting from a liquid-like configuration and an ordered one, as illustrated in Fig. S2 [18] .

Attractive hexagons at high temperatures. -We be- gin by reporting results for the determination of the phase diagram of attractive hexagons. At high temperature, the pressure of attractive hexagons increases monotonically with the density, as illustrated in Fig. 1(a) . We associate to each particle a local density defined as

where H is the Heaviside step function, r c = 50. Results are robust with respect to choice of r c , unless it becomes very small, or of the order of the system size. Figure 1 (b) shows the distribution of the local density, which is always unimodal. The density dependence of the pressure and of the local density distribution exclude the presence of a discontinuous transition with a coexistence phase. We identify the different pure phases investigating the spatial decay of the correlation function of the translational, c(r), and of the bond-orientational order, g 6 (r). The translational correlation function is c(r = | r i − r j |) = e i G·( ri− rj ) , where G is one of the first Bragg peaks, identified by the static structure factor [4, 19] . The bond-orientational correlation function is

Here, n is the number of nearest neighbors of the particle and θ i m is the angle between ( r m − r i ) and a fixed arbitrary axis. The value of k reflects the rotational symmetry of the crystal structure: k = 4 for squares, k = 6 for the other particles.

At high density, the system is in the solid phase. Consistently, we observe the translational correlation function to decay as c(r) ∝ r −η with η ≤ 1/3, a consequence of the Mermin-Wagner theorem [20] , and the bond-orientational correlation function to reach a con- stant, as illustrated in Figs. 1(c) and 1(d) for state point 1 . At lower density, c(r) decays exponentially, while g 6 (r) has a power-law decay, g 6 (r) ∝ r −η6 with η 6 < 1/4. This occurs, for instance, at state points 2 -4 , and indicates that the system is in the hexatic phase. Further lowering the density, the system enters the liquid phase, where both correlation functions decay exponentially.

These findings demonstrate that, at high temperature, LJ hexagons follow the KTHNY scenario [1] [2] [3] . This result agrees with a previous investigation of the melting transition of hard hexagons [6] , the role of attractive forces being negligible at high temperatures.

Attractive hexagons at intermediate temperature.

-As the temperature decreases the equation of state of attractive hexagons flattens, in a range of densities, and develops a Mayer-Wood [21] loop for T 0.53, as illustrated in Fig. 2 (a). Since pressure loops are induced by the interfacial free energy of coexisting phases [4, 22] , this indicates the presence of a first-order transition. Within the coexisting region, the distribution of the local density becomes extremely broad and well described by the superposition of two Gaussian functions, as shown in Fig. 2(b) . Besides, the distribution becomes system-size dependent, with a bimodal character more apparent in larger systems, as we show in Fig. S3 [18] . These findings further support the presence of coexisting phases.

We determine the coexistence boundaries via the Maxwell construction with the pressure curve fitted by either a fifth-or a tenth-order polynomial. Outside of the coexistence region, we identify the pure phases investigating the translational and the bond-orientational correlation functions, as summarized in Figs. 2(c) and 2(d). We observe the solid phase (e.g., 1 ), where c(r) decays algebraically and g 6 (r) is extended, the hexatic phase (e.g., 2 ), where c(r) decays exponentially and g 6 (r) is extended, and the liquid phase (e.g., 7 ) where both correlation functions decay exponentially. These results indicate that, at intermediate temperatures, attractive hexagons follows the mixed melting scenario with a continuous solid-hexatic transition anticipating a discontinuous hexatic-liquid transition.

We visualize the different phases by colour-coding each particle according to the angle ∆α i k between the global Ψ k =

In the solid and hexatic phase, the long-range or quasi-long-range nature of the bondorientational order leads to snapshots with a uniform colour, as in Fig. 2 (e) 1 and 2 . In the liquid phase, Fig. 2 (e) 7 , the snapshot appears almost randomly coloured, due to the short-range of the bond-orientational order. In the coexistence phase, Fig. 2 (e) 5 , the coexistence of hexatic and liquid phases lead to that of regions of uniform colour and regions randomly coloured.

Attractive hexagons at low temperature. -As the tem- perature decreases, the coexistence region widens, and the hexatic phase shrinks, as apparent in Fig. 2(a) . At low enough temperatures, therefore, melting occurs via a first-order liquid-solid transition with no hexatic. The pressure loop and the bimodal character of the local density distribution within the coexistence region, which we illustrate in Fig. 3 , confirms that the system undergoes a discontinuous transition at low temperature. Furthermore, snapshots of the system indicate that the coexisting phases are of solid and liquid type, as in Fig. 3(c) . We can, therefore, exclude an intermediate hexatic phase, further supporting a first-order melting scenario at low temperature.

Shape and temperature dependence of the melting scenario. -The phase diagram of Fig. 4(a) summarizes the results we have obtained so far: for attractive hexagons, melting is of KTHNY type at high temperature and becomes first of mixed type and then first-order as the temperature decreases. The system has no liquid-gas transition, nor a solid-solid transition, as the attraction range, we determine in Table I in the Supplemental Material [18] , is small but not much smaller than the particle size. This suppresses the liquid-gas critical point [23] without promoting a solid-solid transition [23] [24] [25] . In the solid, the body-orientation of the hexagonal particles is long-ranged, thus excluding the presence of a plasticcrystal phase.

We investigate the universality of the role of attraction in the 2D melting by determining the phase diagram for different particle shapes: squares, pentagons and LJ point particles (disks). Squares crystallize in the square lattice, all other shapes in the hexagonal one. For pentagons, shape-frustration is not able to inhibit crystallization at the lowest temperature (T = 0.19) and the highest density (φ = 0.854) we studied. Details on the phase determination are in Figs. S4 and S5 [18] , for squares and pentagons, and elsewhere for disks [26] . The resulting phase diagrams are in Figs. 4(b)-(d) .

At high temperature, the polygonal particles follow the melting scenario previously reported for hard particles [6] , KTHNY for squares and first-order for pentagons. LJ disks follow the mixed scenario as r −12 particles [9] . As in hexagons, in both squares and pentagons, the liquid-gas transition is suppressed, and no plasticcrystal phase occurs. In LJ disks, the liquid-gas critical point occurs at low-temperature and low-density, well outside the parameter space we have investigated.

Regardless of the high-temperature behaviour, melting always occurs via a first-order solid-liquid transition at low temperature, and the coexistence region broadens as the temperature decreases. These findings imply that particles' shapes fix the melting scenario at high temperature, attractive forces at low temperatures. At intermediate temperatures, conversely, both shape and attraction may be relevant. In disks, the hexatic phase disappears as the temperature decreases, making the melting transition first-order. In squares, the equation of state within the tetratic region becomes flat as the temperature decreases without the coexistence region shrinking, as we illustrate in Fig. S4 [18] . Hence, no intermediate mixed scenario separates the high-temperature continuous melting and the low-temperature discontinuous one. In this respect, squares differ from hexagons and disks. In pentagons, the transition is always first-order.

These results consistently demonstrate that attraction influences the melting scenario of 2D systems by promoting the emergence of a coexistence region, if this is not already present in the high-temperature limit, as well as widening it. The widening of the coexisting region leads to the disappearance of the hexatic/tetratic phase, and hence to a first-order melting transition.

Conclusions.

-Different system properties affect the melting scenario in 2D [6, 9, [27] [28] [29] . In this context, the influence of attractive forces was unclear, despite their prevalence at both the molecular and colloidal scale. We have found that attractive forces induce a discontinuous transition and widen the coexistence region at the expense of the hexatic phase, making the low-temperature melting transition first-order. Hence, attractive forces never induces the hexatic phase or widens the hexatic region, at variance with previous speculations [10, 11] . We suggests that attractive forces always promote the discontinuous transition, as we demonstrated this to occur in systems which melt according to different scenarios at high temperature.

Theoretically, our results suggest that the dislocation core energy, E c , is suppressed at low temperatures, in the presence of attractive forces. Conversely, at hightemperature E c /k b T ≫ 1, and a continuous two-step melting scenario may occur, according to the KNTHY theory. We are looking forward to the experimental investigation of our predictions in colloidal systems, where the tuning of the strength of the attractive forces, e.g. via the depletion interaction, should allow for the observation of different melting scenario in the same system, e.g. colloidal hexagonal-shaped particles.

The interparticle interaction of our polygonal particles inherits their discrete rotational symmetry, as the attraction range is small compared to the particle size, as we detail in the Supplemental Material [18] . As the attraction range increases, this discrete rotational symmetry vanishes and the interaction becomes more rotationally symmetric. Hence, while we have not explicitly investigated the role of the attraction range on the phase behaviour, we anticipate that on increasing the attraction range the phase diagrams of the polygonal particles evolve towards that of the LJ point particles.

We acknowledge support from the Singapore Ministry of Education through the Academic Research Fund MOE2017-T2-1-066 (S), and are grateful to the National Supercomputing Centre (NSCC) of Singapore for providing computational resources. We thank Joyjit Chattoraj for helpful discussions

We construct the polygons by lumping together N d point particles equally spaced along the perimeter, as illustrated in Fig. 5(a) . Point particles of different polygons interact via a LJ potential, V LJ (r) = 4ǫ σ r 12 − σ r 6 + c for r ≤ r c = 2.5σ, V LJ (r) = 0 otherwise, where the constant c enforces V LJ (r c ) = 0. The radius of the circumcircle is 5 √ 2/2σ for squares, pentagons and hexagons. The interaction energy U n (r, θ) of two polygonal particles with n sides is the sum of the interactions of their constituent point particles. This interaction depends on the separation between the polygons, r, and their relative orientation, θ. The angular dependence inherits the rotationally symmetry of the particles, U n (r, θ) = U n (r, θ + 2π n ). For the different particles, we illustrate U n (r, θ), with the energy expressed in units of ǫ, and the distances in units of σ, in Fig. 5(b) . For the number of point particles N d we use to construct the polygons, the interaction energy is de-facto linear in N d . The discrete representation of the polygonal particles is thus not affecting our results.

For θ = 0, U n (r, θ = 0) is well described by a LJ potential for extended particles,

(1) I : For the polygonal particles, we provide the values of the parameters ǫn and dn describing the interparticle interaction for θ = 0, given the number of beads N d we use to describe the particles, as well as the particle size r0 and the ratio between the attraction range and the particle size, λn.

as we illustrate in Fig. 5(c) . Hence, for θ = 0, the interaction potential is zero at r 0 = d n + σ, and reaches its minimum value ǫ n at r min = d n + 2 1/6 σ. The parameter λ n = 2(rmin−r0) r0 is a measure of the width of the attractive range with respect to the particle size. We summarized the values of these parameters in Table I .

For LJ disks, we use standard LJ units: ǫ, σ and m are our units of energy, distances and masses, where m is the mass of each particle.

For polygonal particles, we use units facilitating comparisons with colloidal-scale experiments, where the depletion interaction controls the strength and the range of the interparticle interaction, and the area fraction measures the density. The minimum of the interaction energy, ǫ n , is the unit of energy, and the distance at which the interaction potential between aligned particles is zero, r 0 , is the unit of length. The unit of mass is the mass of the particle, mN d , with m mass of the point particles used to construct the extended polygons, and N d their number. The values of these quantities in standard LJ units are in Table I . Interpreting our unit of length as the diameter of the circle inscribed in the polygonal particles, we evaluate their area A n . We measure the degree of crowding via the area fraction φ = ρA n , with ρ number density.

We ensure that our simulations reach thermal equilibrium by checking for the convergence of simulations starting from a liquid-like and a crystalline configuration, at long times. The liquid configuration is prepared by quickly compressing a dilute configuration to the target density. We illustrate this analysis for the hexagonal particles, focusing on three values of the control parameters in the hexatic phase ( 1 ), the hexatic-liquid ( 2 ) and the solidliquid ( 3 ) coexistence regions (see Fig. 6(a) ).

Regardless of the initial configuration, the global bond-orientational parameter reaches the same asymptotic value at long-times, for state points in the hexatic phase and the hexatic-liquid coexistence regions, as we illustrate in Figs. 6(b) and 6(c). We do not always find this convergence for state points within the solid-liquid coexistence region, as exemplified in Fig. 6(d) ; nevertheless, the system always reaches the solid-liquid coexistence phase, as demonstrated by the snapshots of the final configurations we provide in Fig. 6(e) . Hence, the different long-time values of the global bond-orientational parameters reflect the different geometrical shapes of the coexisting phases. As the system coarsen, on a time-scale not accessible in our simulations, the observed difference should disappear.

Equilibration runs of at least t = 10 5 , for point particles, and t = 10 4 , for polygonal systems, have been carried out before collecting data. We have verified that this equilibration time ensured that the system reaches thermal equilibrium at the considered state points. 

We investigate the size dependence of the melting scenario of hexagonal particles at three different temperatures, T = 1.40, 0.53, and 0.35. At these temperatures, melting follows the KTHNY, the mixed and the first-order scenario, respectively. We vary the number of polygons from N = 10661 to 48071.

The equation of state may exhibit, in principle, Mayer-Wood loops across both continuous and discontinuous transition. Across continuous transition, the amplitude of the loop scales as N −1 and the size dependence is lost even in moderate system sizes. Indeed, we do not observe any size dependence in the equation of state at high temperature, as we illustrate in Fig. 7(a) , the system melting via two continuous transitions.

Across discontinuous transitions, a size dependence commonly occurs as the amplitude of the Mayer-Wood loops scales as N −1/2 . We do not observe any size dependence in hexagonal systems even in the presence of discontinuous transition, as illustrated in Figs. 7(b) and 7(c) , the equation of state being flat in the coexistence region, even for the smallest systems we have investigated. This peculiarity makes easy the identification of the coexisting densities.

Consistently with these results, the local density distribution results size-independent at high temperature. At intermediate temperatures, the system melts via the mixed scenario, and the local density distribution is size-dependent in the hexatic-liquid coexistence region. In particular, the distribution broadens with the system size increases, suggesting a bimodal behaviour in larger systems. The size dependence is also weak at T = 0.35, where the system melts via a first-order transition, and within the coexistence region, the local density distribution is bimodal.

Squares and pentagons behave similarly to the hexagons, as concern the size dependence. Conversely, we observe the N −1/2 size dependence across the liquid-hexatic transition of LJ point particles. 

At high-temperature squares melt via continuous transitions as no pressure loop occurs in the equation of state, as in Fig. 8(a) , T = 0.274. As the temperature decreases the pressure flattens, the transition becoming discontinuous below T 0.256. Consistently, at high-temperature, the distribution of the local density is unimodal at all densities, while at low temperature, the distribution is bimodal within the coexistence region, as in Figs. 8(b) and 8(c). We identify the pure phases studying the decay of the correlation functions of the translational and of the bond-orientational order parameters. In squares, the bond-orientational correlation function is g 4 (r = | r i − r j |) = ψ 4 ( r i )ψ * 4 ( r j ) , given the expected square symmetry of the crystal, and the intermediate phase between the crystal and the liquid phase is tetratic, rather than hexatic.

We illustrate the decay of the correlation functions, for the state points 1 -5 identified in Fig. 8(a) , in Figs. 8(d) and 8(e). We observe the solid phase (e.g., 1 ), where c(r) decays algebraically and g 4 (r) is extended, the tetratic (e.g., 2 and 3 ) phase, where c(r) decays exponentially and g 4 (r) decays algebraically (or extended in the observation window), and the liquid one (e.g., 5 ) where both correlation functions decay exponentially. Snapshots of the system with the colour of the squares reflecting the angle between their local bond-orientational order parameter and the global one, we illustrate in Fig. 8(f) , support the above identification of the phases. In the solid and tetratic phase, the uniform colour reflects the long-range or quasi-long-range nature of the bondorientational order. Random colours characterize the liquid phase, as for state point 8 , while in the coexistence phase homogeneously coloured patches coexist with more random ones.

We remark that the temperature dependence of the equation of state in hexagons and square is qualitatively different, although both systems follow the KTNHY scenario at high-temperatures and the first-order one at low temperatures. In hexagons, the hexatic phase gradually shrinks while the coexistence region widens, as the temperature decreases; in squares, conversely, the coexistence phase abruptly replaces the tetratic one. A mixed melting scenario thus occurs in hexagons, but not in squares, at intermediate temperatures. 

In the main text, we have presented the phase diagram in the T vs φ − φ l plane for polygonal systems and in the T vs ρ − ρ l plane for disks, where φ l (ρ l ) is the highest area fraction (number density) of the pure liquid phase, at each temperature. Translating the density axis by φ l facilitates comparing different systems.

We present the same phase diagrams in the T vs φ (or ρ) plane, to ease experimental comparisons, in Fig. 10 . The highest area fraction (or number density) of the liquid phase, and the lowest one of the solid one, have a system-specific temperature dependence. In disks, both densities decrease with the temperature, while in pentagons the opposite occurs. In hexagons and squares, the maximum liquid area fraction decreases while the minimum solid area fraction increases. 

Defects and Geometry in Condensed Matter Physics

for additional information about the interparticle interaction, unit, thermal equilibration, finit-size effects, phase determination for squares and pentagons, and the phase diagram in T -φ plane

and 6 ) where both correlation functions decay exponentially. Pentagons, therefore, melt via the first-order scenario at all temperatures. Snapshots of the system

Equation of state, at selected values of the temperature. Different symbols correspond to different phases, as in the legend. Black lines are from polynomial fits. For selected state points, we illustrate the distribution of the local density in (b) and (c), and the decay of the translational and bond-orientational correlation functions in (d) and (e)