Melting and re-entrant melting of polydisperse hard disks

By using computer simulations, we systematically investigate the equilibrium phase behaviour of polydisperse hard disks in 2D. We find that with increasing the particle size polydispersity, the first-order hexatic-fluid transition becomes weaker and completely switches to be continuous at the polydispersity larger than $7\%$ following the famous Kosterlitz-Thouless-Halperin-Nelson-Young scenario. Simultaneously, the packing fraction range for stable hexatic phase increases dramatically by one order of magnitude compared with the narrow stable range in systems of monodisperse hard disks. More interestingly, in hard-disk systems of slightly higher polydispersity, re-entrant crystal-hexatic and hexatic-fluid transitions occur at high density. The re-entrant transitions were originally proposed in 3D polydisperse hard-sphere systems [Phys. Rev. Lett. 82, 1979 (1999)], but proven impossible due to the strong fractionation [Phys. Rev. Lett. 91, 068301 (2003)]. Therefore, our results suggest a new fundamental difference between the melting transition of polydisperse systems in 2D and 3D.

Melting of two-dimensional solids has been heavily discussed since it was proven that long-wavelength thermal fluctuations prevent the long-range positional order in 2D systems [1][2][3]. A significant attempt to reach the general description of 2D melting was the Kosterlitz-Thouless-Halperin-Nelson-Young (KTHNY) theory, that predicts a new phase of matter, i.e., the hexatic phase, interposed between the usual solid and fluid phases, having quasilong-range orientational order and short-range positional order. Because of the existence of hexatic phase, the crystal melting in 2D was predicted to undergo two continuous transitions from crystal to hexatic and hexatic to fluid, respectively [4][5][6][7]. However, the KTHNY theory does not rule out the first-order transition due to other effects [8], e.g., the grain-boundary induced melting [9,10]. For systems of monodisperse hard disks, arguably the simplest benchmarking model system for testing theories of 2D melting, after a long debate [11][12][13][14][15][16][17][18][19], it was recently settled that the melting of crystal undergoes a two-step process consisting of a continuous crystalhexatic transition followed by a first-order hexatic-fluid transition [20,21], and the shape and softness of particles also play important roles in the 2D melting [22][23][24].
It was found that pinning a small fraction of particles can change the melting transition in hard-disk systems significantly [25,26]. Moreover, simulations of binary hard-disk mixtures showed that the presence of tiny amounts of small particles can eliminate the stability of hexatic phase [27]. These highlight that the melting transition in 2D is subtle. A recent experiment with colloidal hard spheres in 2D [28] appears to support the two-stage melting found in simulation [20,21]. However, most of experimental systems always possess certain degree of (continuous) particle size polydispersity, and a specific question of both fundamental interest and practical relevance is "can the polydispersity change the nature of melting transition in 2D?" To this end, we investigate the equilibrium phase behaviour of a 2D polydisperse hard-disk system (PHDS) with Gaussianlike size polydispersity. We find that with increasing the polydispersity, the first-order hexatic-fluid transition becomes weaker, and completely switches to be continuous following the KTHNY scenario in PHDS with 7% size polydispersity. Concurrently, the packing fraction range for stable hexatic phase increases dramatically by one order of magnitude compared to that in monodisperse hard-disk systems. More surprisingly, in PHDS with slightly higher polydispersity, we observe re-entrant crystal-hexatic and hexatic-fluid transitions at high density, which were proven impossible in 3D polydisperse hard-sphere systems [29].
To model the effect of polydispersity, we consider the system of volume V containing N polydisperse hard disks based on the semigrand canonical ensemble, in which the chemical potential difference between particles of different size is fixed [30][31][32][33]. In this work, we use the following function for the distribution of chemical potential difference where σ is the particle diameter changing from 0 to ∞ with k B and T the Boltzmann constant and temperature of the system, respectively. ν is the polydispersity parameter, and in the ideal gas limit Eq. (1) gives a Gaussianlike particle size distribution centered around σ 0 with the standard deviation ν. This models PHDS in equilibrium with the dilute reservoir of Gaussian-like particle size distribution, e.g., the very top region in sedimentation experiments [28]. We perform Monte Carlo (MC) simulations in the semigrand canonical ensemble (N V T − ∆µ), in which each particle can randomly change its diameter σ under the control of Eq. (1). To accelerate the equilibration and sampling, we implement the rejection free event-chain MC algorithm [34]. By performing MC simulations, we calculate the equation of state (EOS) for a system of 256 2 = 65, 536 hard disks with various ν from 0.005 to 0.08σ 0 . As shown in Fig 1a, one can see that when the polydispersity parameter is small, i.e., ν/σ 0 ≤ 0.05, there is clearly a Mayer-Wood loop in the EOS [35] implying a firstorder transition from fluid with increasing the density of the system. Similar to the monodisperse hard-disk system, i.e., ν = 0, the coexisting phase with fluid in PHDS is a hexatic phase [20,21,33]. To characterize the structural difference between fluid and hexatic phases, we calculate the six-fold bond orientation order parameter Ψ 6 = 1 N N k=1 ψ 6 (r k ) with ψ 6 (r k ) = where θ kj is the angle between the vector connecting particle k with its neighbour j and a chosen fixed reference vector. N k is the number of first neighbours for particle k based on the Voronoi tessellation of the system. As shown in Fig. 1a, the density difference between the coexisting hexatic and fluid phases becomes smaller with increasing ν, and when ν/σ 0 ≥ 0.07, the density jump from fluid to hexatic phase disappears, while Ψ 6 increases dramatically at ρ hex (Fig. 1a inset). This suggests that the first-order fluid-hexatic transition becomes weaker with increasing ν, and changes to be continuous at high polydispersity, i.e., ν/σ 0 ≥ 0.07, following the celebrated KTHNY scenario [4][5][6][7]. It was shown that the system finite size effect is important in studying 2D melting [20,36]. Therefore, we perform MC simulations for PHDS of various numbers of particles from N = 64 2 to 512 2 at ν/σ 0 = 0.07, and the EOSs for different system sizes are shown in Fig. 1b. When the system size is small, i.e., N = 64 2 , there is a pronounced Mayer-Wood loop in EOS, signaturing a first-order fluid-hexatic transition consistent with Ref [33]. Surprisingly, with increas-ing the system size, the first-order fluid-hexatic transition becomes weaker, and changes to be continuous at N ≥ 256 2 . Further increasing N does not change EOS significantly, which suggests the finite size effect negligible in our system of N = 256 2 particles.
With increasing the density of the system, the hexatic phase crystallizes, and the hexatic-crystal transition point can be obtained by checking the pair correlation function g(x, 0) − 1 along the major axis of system switching from an exponential decay to a power law decay upon crystallization [20]. The resulting phase diagram is plotted in the representation of s/ σ vs ρσ 2 0 in Fig. 1c, where s = σ 2 − σ 2 is the actual particle size polydispersity in the system. One can see that with increasing the polydispersity of the system, the density range for stable hexatic phase increases substantially for s/ σ above 0.07. To understand the physics behind this enhanced stability of hexatic phase in PHDS, we calculate the fraction of topological defects in the system using the method in Ref [23]. Particle k is identified as a topological defect if its disclination charge q a = N k − 6 = 0. These topological defects form clusters which can be classified based on the total disclination charge and Burgers vector. The defect clusters having nonzero Burgers vectors and zero disclination charges are called dislocations, of which the existence with few defect clusters of nonzero disclination charges is evidence for the hexatic phase [23]. We ensure that the fraction of particles in the defect clusters with nonzero disclination charge in hexatic phase is always very small in the hexatic phase found in our simulations. The calculated fractions of topological defects and dislocations in our simulations with various ν/σ 0 are plotted as functions of (ρ − ρ hex )σ 2 0 in Fig 2a, where ρ hex is the lowest density of stable hexatic phase, and the vertical dashed lines locate hexatic-crystal transitions. One can i /V is the packing fraction of the system with σ i the diameter of particle i. The state points obtained from simulations at each ν/σ 0 from 0.005 to 0.08 are shown as the symbols, which are color coded with σ /σ 0 . see that typically when the fraction of dislocations decreases to below about 1 ∼ 1.5%, the system crystallizes. When the polydispersity is small, i.e., ν/σ 0 ≤ 0.05, the dependence of both fractions of defects and dislocations on density do not change significantly with increasing ν, and the density ranges for stable hexatic phase below the dashed lines in Fig. 2a are almost the same. However, at ν/σ 0 ≥ 0.07, the fraction of topological defects increases to around 10% at ρ hex , which also boosts the fraction of dislocations to about 4.5% at ρ hex . Along with increasing the density at ν/σ 0 = 0.07 and 0.08, the fraction of dislocations decreases to below the threshold for crystallization at much higher density compared with the system of ν/σ 0 ≤ 0.05. This suggests that the size polydispersity of hard disks creates more topological defects, which subsequently increases the fraction of dislocations in the system destroying the quasi-long-range positional correlation in the solid and driving the formation of hexatic phase.
As the chemical potential difference ∆µ(σ) is fixed in our simulations, the average particle size in the system decreases with increasing density. Therefore, the packing fraction of the system φ does not linearly depend on the density, while the actual particle size distribution remains very close to a Gaussian-like distribution centered around σ [37]. At very high density, when the polydispersity is high enough, e.g., ν/σ 0 0.08, the packing fraction of the system can even decrease with increasing the density. This decreasing packing fraction at high density signatures a re-entrant melting transition, which we explain later. In experiments, one of the most relevant parameters is the packing fraction of the system. Thus we plot the low pressure phase diagram of PHDS in the representation of φ vs s/ σ in Fig. 2b.
Here the low pressure means the range of pressure close to the fluid→hexatic and hexatic→crystal transitions below the re-entrant transitions. One can see that in systems of larger size polydispersity, the average particle size is smaller, and the actual packing fraction range for stable hexatic increases from 0.2 ∼ 0.3% in the monodisperse hard-disk system to about 2 ∼ 3% in the PHDS with s/ σ 0.07. Moreover, at fixed ν, increasing the density of the system along with the fluid-hexatic and hexatic-crystal transitions decreases the actual polydispersity s/ σ in the system, which is due to the formation of more ordered structures eliminating the particle size fluctuation in the system.
Another interesting feature of the phase diagram in Fig. 1c is that when the size polydispersity of hard disks is around 0.08 to 0.1, at very high density, increasing density triggers re-entrant melting transitions from crystal to hexatic and hexatic to fluid. As shown in Fig. 3a, for PHDS with ν/σ 0 = 0.08, with increasing the density from ρσ 2 0 = 1.6 to 1.6424, the system transforms from a hexatic phase with the short-range positional correlation, i.e., an exponential decay g(x, 0) − 1 ∼ exp(−x), to a crystal with quasi-long-range positional order, i.e., a power law decay g(x, 0) − 1 ∼ x −α with α ≤ 1/3, where the six fold orientation correlation function g 6 (r) = ψ * 6 (r + r)ψ 6 (r ) [38] remains almost the same (Fig. 3b). At much higher density, increasing the density from ρσ 2 0 = 3.64 to 3.67, g(x, 0) − 1 changes again from a power law decay to an exponential decay with almost the same g 6 (r) (Fig. 3b), which suggests a re-entrant transition from crystal to hexatic phase. For hard-disk systems with higher polydispersity, i.e., ν/σ 0 = 0.082 ∼ 0.0835, we plot Ψ 6 as functions of ρσ 2 0 in Fig. 3c, and one can see that with increasing ρ, Ψ 6 first increases from 0 to around 0.8 indicating the formation of an ordered phase, and further increasing ρ leads to the drop of Ψ 6 to 0 implying a re-entrant melting transition into a disordered phase. By checking the scaling of g(x, 0) − 1 in the system, we ensure that the ordered phase formed is a hexatic phase. To understand the mechanism behind the re-entrant melting of hexatic phase, we plot s/ σ as functions of ρσ 2 0 in the lower panel of Fig. 3c, from which one can find a clear correlation between the change of Ψ 6 and s/ σ . Especially, at high density, during the decrease of Ψ 6 , s/ σ increases dramatically to the value even higher than that in the low density fluid. This suggests that the entropy associated with particles of high polydispersity stabilizes the amorphous structures against ordered ones at high density. Moreover, along with the decrease of Ψ 6 with density at high pressure, we do not observe any Mayer-Wood loops in EOS [37], indicating that the re-entrant melting of hexatic phase at high density is continuous similar to that in soft particle systems [39]. Moreover, it was shown that the strong particle size fractionation in polydisperse hard spheres eliminates the possibility of re-entrant melting in 3D polydisperse hard-sphere crystals [29]. Therefore, in Fig. 3d, we plot the probability density distribution of particle size p(σ) for PHDS of various density with ν/σ 0 = 0.0835. One can see that within the whole density range, there is always a single peak in p(σ) at fixed ρ, indicating no particle size fractionation in the system. This also suggests that re-entrant melting can indeed exist in 2D polydisperse particle systems, if there is no particle size fractionation [40]. In addition, to check whether the amorphous phase formed at high density is a fluid or glass, we perform event driven molecular dynamics (EDMD) simulations starting from the equilibrated configurations obtained from our MC simulations at ν/σ 0 = 0.0835, and the calculated mean square displacement (MSD) for various densities are shown in Fig. 3e. One can see that not only the diffusion coefficient increases by nearly two orders of magnitude along with the re-entrant melting of hexatic phase (Fig. 3e inset), but also the plateau in MSD almost disappears. This suggests that the high density amorphous phase is a diffusive fluid, and in the equilibrium sedimentation of such PHDS, there should be a "floating hexatic phase" sandwiched in two fluids. For PHDS with ν/σ 0 ≥ 0.084, with increasing the density, we do not observe any ordered phase in our simulations.
In conclusion, by performing large scale computer simulations, we investigate the effect of particle size polydispersity on the phase behaviour of hard-disk systems. We find that with increasing the size polydispersity of hard disks from zero, the first-order melting transition from hexatic phase to fluid becomes weaker, and completely switches to be continuous at ν/σ 0 ≥ 0.07 following the celebrated KTHNY scenario. Simultaneously, the density range for stable hexatic phase gets substantially enlarged. Compared with monodisperse hard-disk systems, where hexatic phase is only stable in a very small range of packing fraction of 0.2 ∼ 0.3% [20], the packing fraction range for stable hexatic phase in PHDS with s/ σ 0.07 is about 2 ∼ 3%. This suggests new directions in searching for the hexatic phase in 2D polydisperse particle systems. More interestingly, in PHDS with even higher polydispersity, i.e., 0.08 ≤ ν/σ 0 0.0835, we find that at very high density, increasing density of the system can trigger re-entrant transitions from crystal to hexatic and hexatic to fluid phases depending on ν. In polydisperse systems, re-entrant transitions, especially the re-entrant melting, were originally proposed theoretically for 3D hard-sphere crystals [40], but proven impossible due to the strong fractionation in the 3D hard-sphere system that was not taken into account in the theory [29]. We find that the absence of strong fractionation in 2D polydisperse systems can indeed induce re-entrant transitions at high density, suggesting a new fundamental difference between the melting transition of polydisperse systems in 2D and 3D.