Melting and re-entrant melting of polydisperse hard disks

Because of long-wavelength fluctuations, the nature of solids and phase transitions in 2D are different from those in 3D systems, and have been heavily debated in past decades, in which the focus was on the existence of hexatic phase. Here, by using large scale computer simulations, we investigate the melting transition in 2D systems of polydisperse hard disks. We find that, with increasing the particle size polydispersity, the melting transition can be qualitatively changed from the recently proposed two-stage process to the Kosterlitz-Thouless-Halperin-Nelson-Young scenario with significantly enlarged stability range for hexatic phase. Moreover, re-entrant melting transitions are found in high density systems of polydisperse hard disks, which were proven impossible in 3D polydisperse hard-sphere systems. These suggest a new fundamental difference between phase transitions in polydisperse systems in 2D and 3D. The hard disk model is generally applied to study melting in two dimensional colloidal solids, which for idealized 2D systems proceeds through a solid to hexatic - hexatic to fluid process, but impurities perturb the hexatic phase for real systems. The paper reports Monte Carlo and molecular dynamics simulations of a 2D system of polydisperse hard disks, finding that increasing polydispersity decreases stability of hexatic phase, and that even for polydisperse systems there are re-entrant transitions at high density, which is not observed for 3D systems.

M elting 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, which has quasi-long-range orientational order and short-range positional order. The hexatic phase is predicted to be interposed between the usual solid and fluid phases. Therefore, the melting in 2D was predicted to undergo two continuous transitions from solid to hexatic and hexatic to fluid, respectively [4][5][6][7][8] . However, the KTHNY theory does not rule out the first-order transition due to other effects 9 , e.g., the grain-boundary induced melting 10,11 . In arguably the simplest benchmarking model system in 2D, i.e., the system of monodisperse hard disks, the physics of melting transition has been debated for long [12][13][14][15][16][17][18][19][20] . It was recently settled that the melting of solids in the system of monodisperse hard disk undergoes a two-stage process consisting of a continuous solid-hexatic transition followed by a first-order hexatic-fluid transition 21,22 , and the shape and softness of particles also play important roles in the 2D melting [23][24][25] .
It was found that pinning a small fraction of particles can change the melting transition in hard-disk systems significantly 26,27 . Moreover, simulations of binary hard-disk mixtures showed that the presence of tiny amounts of small particles can eliminate the stability of hexatic phase 28 . These highlight that the melting transition in 2D is subtle. A recent experiment with colloidal hard spheres in 2D 29 appears to support the two-stage melting found in simulation 21,22 . However, most of experimental systems possess certain degree of (continuous) particle size polydispersity, and polydisperse hard disks have also been widely employed as a model system to investigate the glass transition [30][31][32] . Yet the effect of polydispersity on the nature of phase transitions in 2D remains unknown. To this end, we investigate the equilibrium phase behavior of a 2D polydisperse hard-disk system (PHDS) with Gaussian-like particle 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 around 7% size polydispersity. Concurrently, the packing fraction range for stable hexatic phase increases significantly 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 solid-hexatic and hexatic-fluid transitions at high density, which were proven impossible in 3D polydisperse hard-sphere systems 33 .

Model.
To model the effect of polydispersity, we consider a 2D 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 [34][35][36][37] . 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 Gaussian-like 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 29 .
Phase diagram. By performing Monte Carlo (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 38 implying a first-order 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 21,22,37 . To characterize the structural difference between fluid and hexatic phases, we calculate the sixfold bond orientation order parameter hΨ 6 i ¼ 1 where θ kj is the angle between the vector connecting particle k with its neighbor 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 significantly at ρ hex (Fig. 1b). Here we note that |〈Ψ 6 〉| in an infinitely large system, should be zero for hexatic phase, while in our large finite system it can be a positive number indicating the existence of some orientational order in the system. However, we do not use |〈Ψ 6 〉| to identify the existence of hexatic phase, for which we always check the change of correlation functions to confirm. This suggests that the first-order fluidhexatic 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] . This, to some extent, is similar to the melting of soft spheres 23 . Because with increasing the particle size polydispersity, the distribution of distance between the nearest neighbors becomes wider, which is similar to introducing a "soft" repulsion between particles, and it has been also found in the melting of soft spheres in 2D, that the transition type can switch from a first-order transition to a continuous KTHNY scenario with increasing the "softness" of the potential 23 .
Finite size effect on the melting in polydisperse hard disks. It was shown that the system finite size effect is important in studying 2D melting 21,39 . 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. 1c. When the system size is small, i.e., N = 64 2 = 4096, there is a pronounced Mayer-Wood loop in EOS, signaturing a firstorder fluid-hexatic transition consistent with ref. 37 . Interestingly, with increasing 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. The change from a first-order like transition to a continuous transition with increasing system size can be interpreted as following. A finite 2D system with periodic boundary conditions can be seen as a system wrapped on a torus, and this increases the cooperation between neighboring particles as a result of extra "communication" via paths encircling the torus, which was shown being able to change the transition type from continuous to first order like in small systems 40 . By contrast, this effect does not appear in systems with open free boundaries, while in real experiments, the open boundary walls can induce order in the fluid. Moreover, for ν/σ 0 < 0.07, we ensure that most of the coexisting densities do not change significantly with further increasing the system size to N = 512 2 , while the only exception is for systems at ν/σ 0 = 0.05 (see Supplementary Fig. 1), of which the exact boundary is out of the reach with our present computation capability.
Stabilization of the hexatic phase. With increasing the density of the system, the hexatic phase solidifies, and the hexatic-solid 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 solidification 21 . The resulting phase diagram is plotted in the representation of s/ 〈σ〉 vs ρσ 2 0 in Fig. 1d, where s ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi hσ 2 i À hσi 2 q 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. 24 . 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 24 . Defect clusters having nonzero Burgers vectors and zero disclination charges are called dislocations, and the dislocations as well as the defect clusters with nonzero disclination charges can destroy the bond orientational order in 2D solids. Although the total amount of topological defects is not directly related to the stability of hexatic phase, the dislocations are evidence for the hexatic phase along with few defect clusters with nonzero disclination charges 24 . The fraction of particles in defect clusters with nonzero disclination charge is much smaller, i.e., more than one order of magnitude smaller, than that of dislocations in the hexatic phase found in our simulations, although the changing trends of these two types of defects are very similar (see Supplementary Fig. 2).
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-solid transitions. One can see that typically when the fraction of dislocations decreases to below about 1~1.5%, the system solidifies. 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 , and also the fraction of dislocations increases 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 solidification at much higher density compared with the system of ν/σ 0 ≤ 0.05. The fraction of defect clusters with nonzero disclination charges also changes similarly in systems of different polydispersity ( Supplementary Fig. 2). Therefore, this suggests that the size polydispersity of hard disks creates more topological defects, which could subsequently increase the fraction of dislocations as well as the defect clusters with nonzero disclination charges in the system. This destroys the quasi-long-range positional correlation in the solid 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 ϕ ¼ π 4 P i σ 2 i =V does not linearly increase with the density, while the actual particle size distribution remains very close to a Gaussian-like distribution centered around 〈σ〉 (see Supplementary Fig. 3). 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 in μVT ensemble, which is different from the NVT ensemble 41 . 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 → solid transitions below the re-entrant transitions. One can see that in systems of larger particle 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-solid 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.
Re-entrant meltings. Another interesting feature of the phase diagram in Fig. 1d 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 solid 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 solid with quasi-long-range positional order, i.e., a power law decay g(x, 0) − 1~x −α with α ≤ 1/3, where the sixfold orientation correlation function g 6 ðrÞ ¼ hψ Ã 6 ðr′ þ rÞψ 6 ðr′Þi 42 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 solid 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/〈σ〉. Here we ensure that our simulations at high density have been well equilibrated and independent with the initial configurations ( Supplementary Fig. 4). Actually, it has been shown that the equilibration in systems of continuously polydisperse particles is much faster than that in monodisperse systems [43][44][45] . At high density, during the decrease of |〈Ψ 6 〉|, s/〈σ〉 increases significantly to the value even higher than that in the low density fluid. This suggests that by increasing the density of the system, the standard deviation of particle size increases, and the combinatorial entropy in the system associated with the variation of particle size increases, which stabilizes the amorphous structures against the ordered ones in the system. Moreover, along with the decrease of |〈Ψ 6 〉| with density at high pressure, we do not observe any Mayer-Wood loops in EOS ( Supplementary Fig. 5). This implies that the re-entrant melting of hexatic phase at high density is continuous, which is similar to that in soft particle systems 46,47 . Furthermore, 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 solids 33 . 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 48 . Here we note that although all our simulations starting with different particle size distributions converge to the same result at fixed polydispersity and density, there still exist a possibility that the system may posses an equilibrium multimodal particle size distribution which is not accessible with direct simulations 48 .
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 Here ρ and ρ hex are the density of the system and the lowest density of stable hexatic phase, respectively. The vertical dashed lines indicate the hexatic-solid transition points. b Low pressure phase diagram of polydisperse hard disks in the representation of ϕ vs s/〈σ〉, where ϕ 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 . The error bars are smaller than the symbols configurations obtained from our MC simulations at ν/σ 0 = 0.0835. In the EDMD simulation, the system evolves via a time-ordered sequence of elastic collision events, which are described by Newton's equations of motion. The spheres move at constant velocities between collisions, and the velocities of the respective particles are updated when a collision occurs. All collisions are elastic and preserve energy and momentum 49 . 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 of Gaussian-like particle size distribution, 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.

Discussion
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% 21 , 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 solid to hexatic and hexatic to fluid phases depending on ν. In polydisperse systems, re-entrant transitions, especially the re-entrant melting, were originally predicted theoretically for 3D hard-sphere crystals 48 , but proven impossible due to the strong fractionation that was not taken into account in the theory 33 . We find that the absence of strong fractionation in 2D polydisperse systems can indeed induce re-entrant transitions at high density. These show a new difference between phase transitions in polydisperse hard-sphere systems in 2D and 3D, and suggest a new direction in investigating phase transitions in 2D systems by considering the particle size polydispersity, which was overlooked in the past.

Methods
We perform MC simulations in the semigrand canonical ensemble (NVT − Δμ) using a square simulation box with periodic boundary conditions in both directions, 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, with which the pressure P in the system can be calculated using the mean excess chain displacement 50 where x final and x initial are the initial and final positions of each event chain along the direction of the chain, respectively. L c and ρ are the chain length and number density of the particles, respectively, with 〈⋅〉 chains calculating the average over all event chains. For each simulation, we perform about 10 8 -10 9 MC sweeps, and each MC sweep on average consists of 1 event chain with the chain length of L c = 100σ 0 and 500 trials of randomly changing the diameter of a randomly selected particle.

Data availability
All data that support the findings of this study are available from the corresponding author on reasonable request.

Code availability
Codes used for the numeric simulations are available on request from R Ni at Nanyang Technological University, Singapore.