Shape and symmetry determine two-dimensional melting transitions of hard regular polygons

The melting transition of two-dimensional (2D) systems is a fundamental problem in condensed matter and statistical physics that has advanced significantly through the application of computational resources and algorithms. 2D systems present the opportunity for novel phases and phase transition scenarios not observed in 3D systems, but these phases depend sensitively on the system and thus predicting how any given 2D system will behave remains a challenge. Here we report a comprehensive simulation study of the phase behavior near the melting transition of all hard regular polygons with $3\leq n\leq 14$ vertices using massively parallel Monte Carlo simulations of up to one million particles. By investigating this family of shapes, we show that the melting transition depends upon both particle shape and symmetry considerations, which together can predict which of three different melting scenarios will occur for a given $n$. We show that systems of polygons with as few as seven edges behave like hard disks; they melt continuously from a solid to a hexatic fluid and then undergo a first-order transition from the hexatic phase to the fluid phase. We show that this behavior, which holds for all $7\leq n\leq 14$, arises from weak entropic forces among the particles. Strong directional entropic forces align polygons with fewer than seven edges and impose local order in the fluid. These forces can enhance or suppress the discontinuous character of the transition depending on whether the local order in the fluid is compatible with the local order in the solid. As a result, systems of triangles, squares, and hexagons exhibit a KTHNY-type continuous transition between fluid and hexatic, tetratic, and hexatic phases, respectively, and a continuous transition from the appropriate"x"-atic to the solid. [abstract truncated due to arxiv length limitations].

The phase behavior of two-dimensional (2D) solids is a fundamental, long-standing problem in statistical mechanics. Whereas three-dimensional (3D) solids characteristically exhibit first order (or discontinuous) melting transitions, 2D solids can melt by either continuous or first order melting transitions and may exhibit an intermediate, so-called "x-atic" ordered phase that is somewhere between a fluid and a solid. Previous studies [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18] that examine two-dimensional melting find three distinct scenarios [19,20]. One, the system can exhibit a continuous fluid-to-x-atic-to-solid transition. The x-atic phase has quasi-long-range (power-law decay) correlations in the bond order but only short-range (exponential decay) correlations in positional order. The hexatic phase, with six-fold bond order, is the most well-known example. The existence of continuous fluid-to-solid transitions was predicted by the Kosterlitz-Thouless-Halperin-Nelson-Young (KTHNY) theory [21][22][23] and has been confirmed in experiments with electrons [1] and spherical colloids [2,4,8,9]. The KTHNY theory of two-step melting is based upon the behavior of topological defects in the form of dislocations and disclinations. The theory envisions that pairs of dislocations unbind to drive * sglotzer@umich.edu the transition from solid to hexatic, and then pairs of disclinations unbind to drive the transition from hexatic to fluid. Two, the system can exhibit a first-order fluidto-solid (or solid-to-fluid) transition, with no intervening phase. Both this and the first scenario were realized in a system of charged polystyrene microspheres, depending on the particle diameter, which was postulated to have an effect on defect core energies [7]. Three, the system can exhibit a first-order fluid-to-x-atic and a subsequent continuous x-atic-to-solid transition. This combination of transitions is intermediate to the one-step fluid-tosolid first order transition and the two-step continuous KTHNY behavior. It was first experimentally observed in neutral micron-scale colloidal spheres [5], and has been observed recently in simulations of hard disks in two dimensions [12,24] and under quasi-2D confinement of hard spheres where out-of-plane fluctuations are limited [14].
All three melting scenarios have been observed in experimental studies of different systems, with a variety of long and short range interactions. Recent simulation work [18] finds two of the three scenarios: point particles with hard core repulsion interactions follow the third scenario and softer potentials lead to continuous melting. In this paper, we report the occurrence of all three distinct melting scenarios in a single family of hard, regular polygons. Hard polygons have a rotational degree of freedom, which creates the possibility for more complex entropic forces and more diverse solid phases than observed for hard disks. By varying the number of polygon edges, we show that the melting transition scenario for a system of any given polygon is determined by the anisotropy of emergent entropic interactions along with the symmetry of the particles relative to that of the solid phase. In particular, we show that systems of hexagons are a perfect realization of the KTHNY theory, exhibiting melting from the solid to the hexatic phase with an increase in the dislocation density, then from the hexatic to the fluid with an increase in the disclination density. We find that systems of triangles and squares also show a continuous KTHNY-type melting transition, while systems of pentagons and 4-fold pentilles have a first order melting transition that occurs in a single step. Finally, we show that systems of regular polygons with n ≥ 7 behave like disks with a first order fluid-to-hexatic and continuous KTHNY-type hexatic-to-solid transition.
We focus our study on hard, convex, regular polygons because we aim to discover general unifying principles of 2D melting by filling in gaps in existing literature: regular triangles [13,16,25], squares [11,17,26], pentagons [3,10], and heptagons [10] have been previously studied by both experiment and simulation. These studies were instrumental in identifying and characterizing possible intermediate phases (triatic, tetratic, and hexatic). We present new results for regular hexagons, octagons, etc. up to 14-gons and clarify the results of previous simulations of hard polygons using very large simulations to conclusively determine the orders of the various melting transitions, where previous studies were too small to be conclusive. Section II includes detailed comparisons between our results and previous simulation and experimental results.
We demonstrate that changing only the number of edges on convex regular polygons is sufficient to generate a rich array of different melting behavior, including all three known 2D melting scenarios. We leave for future studies investigation of, e.g. rounding of polygons, where experiment [13] and simulation [27] have revealed additional phases.

I. METHODS
We investigate large systems of N identical polygons with n edges that interact solely through excluded volume interactions in a box of area A box . Particle a has position r a and orientation angle θ a . The circumcircle diameter of the polygons is denoted as σ. The majority of our work focuses on regular polygons (n-gons) with the area of a single particle A = σ 2 n/2 sin(2π/n). We also include in our study, the 4-fold pentille [28], which is the Voronoi cell of the Cairo pentagon pattern and thus tiles space (see figure S16 in the Supplementary Information for the tiling configuration). Figure 1 shows all thirteen polygons and summarizes their melting behavior.

A. Computational strategy
We executed Monte Carlo (MC) simulations as large as 1024 2 particles to obtain high precision equations of state and sample long range correlation functions to conclusively identify hexatic and solid phases, which recent hard disk simulations [12,24] found necessary. The MC simulations were performed in the isochoric (constant area) ensemble using the HPMC [29] module of HOOMDblue [30][31][32]. HPMC (for hard particle Monte Carlo) is a parallel simulation implementation on many CPUs or GPUs using MPI domain decomposition.
Each simulation begins with the particles placed in a square, periodic simulation box. We perform simulations for a long enough time to reach and statistically sample thermodynamic equilibrium, see Supplementary Information section II for a complete simulation protocol description, and figure S14 for an example. Typical simulation runs initially form many domains in the system, which coalesce together over a long equilibration period of several hundred million trial moves per particle.
Occasionally, two infinite, twinned domains form in the system. Such configurations are metastable, they are at a higher pressure than a corresponding single domain and the larger domain grows very slowly into the smaller one as the simulation progresses. We remove stuck simulations and rerun them with different random number seeds until we obtain a cleanly equilibrated single domain sample, except in a few cases where multiple attempts to do so failed (e.g. φ = 0.714 in figure S10).

B. Equation of state
We compute the equation of state P * (φ) using volume perturbation techniques [34,35] to measure pressure in isochoric simulations. We report pressure in reduced units P * = P σ 2 /k B T . In addition to the system density φ = N A/A box , we determine the averaged local density field on a grid, where we choose the cut-off r c = 20σ and H is the Heaviside step function.
With the fluid density φ f and the density of the solid (or hexatic) phase φ s the latent heat L = P ∆A box of a first-order phase transition is where P * c is the coexistence pressure estimated by Maxwell construction.

C. Order parameters
We use three order parameters that were previously effective at identifying the hexatic phase in the hard disk system [12,24]. Each order parameter is a complex number on the unit circle. We visualize order parameter fields directly in the x-y-plane of the system by mapping the complex values of an order parameter to a color wheel (see Figure 2). Short-range order shows up as colors rapidly cycling through the color wheel, quasi-long-range order appears as patches, long-range order appears as a single solid color across the box, and two-phase regions show two separate behaviors in a single simulation box. Independent simulation runs result in different system orientations in the box. We rotate the order parameters so that the average in a given frame is colored green so that images may be compared by eye.
The positional order parameter identifies how well the position r a of particle a fits on a perfect lattice with reciprocal lattice vector q 0 , as depicted in Figure 2. When all particles have the same phases in χ, they are in a perfect solid. Defects cause χ to rotate. We choose q 0 as the brightest peak in the structure factor [12] computed with the following procedure: (i) Initialize a density grid with roughly 8 × 8 pixels per particle. (ii) At the center of each particle, place a Gaussian on this grid with standard deviation 1 10 σ. (iii) Take the fast Fourier transform (FFT) of the density grid to get the discretized S( q). (iv) Smooth with a Guassian filter, standard deviation 2 pixels S( q) → S smooth ( q).
(v) Choose q 0 from the location of the brightest pixel in S smooth ( q).
The bond orientation order parameter for k-fold rotational symmetry We map order parameters to a color wheel for visualization. Each order parameter is rotated by its average in a given frame so that the average order parameter is colored green. The color wheel is the color part of the cubehelix [37] color map at constant apparent luminance (v = 0.5, γ = 1, s = 4.0, r = 1, h = 1).
identifies the orientation of the p nearest neighbors around particle a. Here α ab is the angle the separation vector r b − r a makes with respect to the positive x-axis, and NN p (a) is the set of p nearest neighbors of a, see Figure 2 for a graphical definition. We omit p in the subscript when it is equal to k and write ψ a k = ψ a k,k . Some authors suggest using a morphometric approach [36] to compute the bond orientation order parameter, which requires computing a Voronoi diagram of the system of particles. We do not adopt this scheme because we find the use of p fixed neighbors sufficient as it generates order parameter fields fully consistent with the defects present in the system.
The body orientation order parameter identifies the orientation of a particle accounting for sfold symmetry. θ a is the angle that rotates particle a from a reference frame into a global coordinate system (see Figure 2). It allows us to analyze the presence of rotator phases in which ξ a s decays to zero rapidly as a function of the separation distance.

D. Correlation functions
Correlation functions measure the behavior of the order parameters as the separation r ba = | r b − r a | between a pair of particles increases. The correlation functions for bond orientation order and body orientation order are implemented in our analysis code as where v = ψ k,p or ξ n and v * is the complex conjugate. The sampling S(m, N ) randomly selects m particles out of N without replacement.
One can compute a correlation function from the positional order parameter χ, but it is extremely sensitive to the choice of q 0 . We observe that peaks misidentified by only a single pixel result in an apparent lack of quasilong-range positional order. Instead, we follow Ref. [12] and compute positional correlation functions from which oscillates, so we show only identified peaks. The signal then decays to about 10 −3 at large r. The function g ψ ( r, α) is the discretized 2D pair correlation function obtained by correlated averaging over individual measurements with index i, where R(β) is the rotation matrix that rotates a vector by the angle β, ∆r is the bin size, and d κ is the coarsegrained delta function Since the system rotates from frame to frame, it must be aligned to the bond orientation order parameter before averaging. To do this, we compute g ψ ( r i , α) over many frames with large m. We align each frame using the bond orientation order parameter averaged over all particles in the frame, α = arg( ψ k,p ). Averaging the separately aligned g ψ ( r i , α) and then computing C g (r) significantly reduces noise [24].

E. Sub-block scaling analysis
The sub-block scaling of χ is a sensitive measure of the density of the hexatic-to-solid transition [6,15]. The first density that sits on or under a line of slope −1/3 is at the hexatic-to-solid transition. Similarly, sub-block scaling in ψ determines the density at which the hexatic melts into the fluid, with a line of slope −1/4.
We perform a sub-block scaling analysis on the positional order parameter χ and the bond orientational order parameter ψ a k,p . For this analysis the simulation box of length L is divided into squares of side length L b .
where v = ψ k,p or χ, N L b is the number of particles in a sub-block and denotes averaging over sub-blocks in the same snapshot. Standard errors are calculated by averaging the values of v(L b ) over independent frames.

F. Topological defect analysis via cell-edge counting
We generate statistics on topological defects using a Voronoi tessellation of the set of particle centers to count edges of Voronoi cells. Each particle a is assigned the number of adjacent Voronoi cells n a and the disclination charge q a = n a − 6. In the presence of wellseparated topological defects, disclinations can be found by identifying particles with non-zero disclination charge, and dislocations correspond to bounds pairs of a fivecoordinated particle (q = −1) and a seven-coordinated particle (q = +1). However, while disclination charges with |q| > 1 are very rare in our hard particle systems, particles with non-zero disclination charge are often not clearly separated or bound into pairs but instead agglomerate into larger clusters. This makes it ambiguous to identify the locations of individual disclinations and dislocations.
To overcome this ambiguity, we cluster defects if they are in adjacent Voronoi cells and calculate the total disclination charge q = q a of the cluster [38]. For disclination-neutral clusters (q = 0), the total Burgers vector b of the cluster can be determined from the disclination charges and their positions, whereẑ is the unit vector along the out-of-plane axis and the cross product is performed in 3D space. Because Burgers vectors are topologically restricted to be lattice vectors, the vector we obtain from this computation is then snapped to its closest lattice point in a hexagonal lattice whose lattice spacing is that of the ideal lattice at the density of the simulation frame. We count the number of particles in three types of defect clusters, overall neutral (q = 0, b = 0), Burgers-charged (q = 0, b = 0), and disclination-charged (q = 0) as a function of density.

G. Phase determination
We identify first-order transitions using the following criteria: (i) a two-phase region is evident in isochoric sim-  loop [12,24] which decreases in height as N increases. In contrast, continuous transitions have a monotonic equation of state and only a single phase at a single density is present in each frame across the entire transition. The signature of two phase regions includes a bimodal local density histogram along with two different correlation length behaviors seen in the ψ and χ order parameter fields coincident with the low and high local densities. In all cases where we find a Mayer-Wood loop, we observe two-phase regions, and in all cases where we find a continuous equation of state, we observe only single phases across the transition. Correlation lengths vary significantly in the simulations. In the fluid phase, positional order χ decorrelates instantly with correlation lengths of only 1σ. Bond orientational order ψ k,p persists a bit further with correlation lengths of tens of σ but still exhibits clear short-range order. Previous authors have used the dynamic Lindemann criterion to locate melting transition points [8,39]. We focus on the correlation functions here because of their utility in discriminating hexatic and solid phases, which the dynamic Lindemann criterion is unable to do. In the hexatic phase, χ oscillates visibly through the full color wheel along a given direction forming stripes (Fig. 3b), a behavior indicative of short-range order. Unbound dislocations exist at the end of each stripe. For continuous phase transitions, the oscillation period gets larger as density increases. Solid phases lead to patchy motifs in χ (quasi-long-range order). Two-phase regions show a combination of two of these motifs in a single system.
Our observations confirm in all cases that we have run simulations sufficiently large to allow the fluctuations in the order parameter fields of the fluids and x-atics to occur several times across the box. Only in the solid phase regions are the order parameter fields essentially constant and vary by less than one period.

II. RESULTS AND DISCUSSION
We performed identical analyses for each of the thirteen polygons investigated in this work. Within the main text we include representative plots and snapshots to illustrate examples and explanations of the three melting scenarios. The Supplementary Information contains detailed plots for all of the shapes in figures S1-S13, including the equation of state, local density histogram, subblock scaling analysis, correlation functions, snapshots colored by all of the order parameters and structure factors.

A. KTHNY behavior
Our data show that systems of triangles, squares, and hexagons have continuous fluid-to-x-atic and continuous x-atic-to-solid transitions. We refer to the monotonic equations of state in Figure S1, Figure S2, and Figure S5 for the triangle, square, and hexagon data, respectively. Triangles and hexagons have hexatic order with a signal in ψ 6 , and squares have tetratic order with a signal in ψ 4 . Figure 3a-c demonstrates the continuous transitions for hexagons, with a single density at each state point and gradually developing order in the order parameter fields.

Decay of correlation functions
The first state point showing quasi-long-range bond orientation order in the hexagon system is φ = 0.686. At this density, positional order is extremely short-ranged, persisting for only a few particle diameters. As density increases, the decay length of the positional order becomes longer, up to hundreds of particle diameters at φ = 0.702, and begins to diverge. At φ = 0.708, positional order switches to quasi-long-range and the system is in the solid phase as determined by sub-block scaling. KTHNY theory predicts a slope of −1/4 in ψ 6 scaling at the fluid-to-hexatic transition, a perfect match for the φ = 0.686 line in Figure 3c. The theory also predicts a slope of −1/3 for χ scaling, which similarly matches the scaling for φ = 0.708. At the same density, the directly computed positional correlation C g (r) length is beginning to diverge so it is difficult to determine the exact density of the hexatic-to-solid transition from correlation functions alone.

Topological defects and local order
The KTHNY theory envisions a picture of tightly bound defects transitioning to free dislocations at the solid-hexatic transition, and to free disclinations at the hexatic-fluid transition. Interestingly, both Bernard et al. [12] and this work show that defects are not free, but rather form large and complex clusters, often highly anisotropic, indicating medium-ranged effective interdefect interactions. In the extreme case that defects cluster into one-dimensional strings, such a scenario would lead to grain-boundary induced melting [40,41]. Figure 4c shows the density dependence of the number of particles comprising clusters of defects that are overall neutral, those with a Burgers charge, and those with disclination charge. A large number of free dislocation clusters with non-zero Burgers vectors simultaneous with few free disclinations is evidence for the hexatic phase. We observe that the count of particles in Burgerscharged clusters begins increasing just below φ = 0.710 at the solid-to-hexatic transition. The count of particles in disclination-charged clusters remains low, and starts ramping up slowly in the middle of the hexatic phase at φ = 0.696. Throughout the hexatic phase, we find many more particles in Burgers-charged clusters than in disclination-charged clusters. This is consistent with the two-step KTHNY melting scenario and the continuous phase transition we observe for hexagons.  The defect count algorithm is expensive, so we do not run it on the largest size simulations we performed. We also show the phase diagram overlaid on the counts. Error bars at one standard deviation are smaller than the symbol size. Figure 5a shows the structure of the defects in a system of hexagons at φ = 0.690, in the hexatic phase. The density of free dislocations (red and green pairs) is much higher than that of free disclinations (there is a single lone green particle in the image), but the structure is more complicated because the defects combine into large clusters. The left panel of Supplementary Movie 1 shows how defects migrate as the Monte Carlo simulation progresses via local moves for a simulation of hexagons at φ = 0.690. Dislocations are unstable and highly mobile, quickly hopping between sites and popping in and out of existence. We find that the count of nearest neighbors is very sensitive to small changes in particle coordinates.
We also compute the isoperimetric quotient IQ = 4πA/P 2 for each Voronoi cell, where A is the area and P is the perimeter of the cell. Regular polygons have the value IQ = π/ n tan π n . More elongated Voronoi cells have low IQ and indicate locations of large crystal deformity, which makes IQ an approximate indicator of defects. Figure 5b and the right panel of Supplementary Movie 1 color particles by IQ −IQ. Blue indicates large deviation from IQ and yellow indicates little deviation. If IQ > IQ , the particle is darkened, as these particles are closest to having regular crystal environments.
The coloring scheme based on the isoperimetric quotient IQ is a continuous quantity and, as seen in the movie, is less sensitive to small particle displacements. In contrast, the number of edges of a Voronoi cell changes discontinuously across Monte Carlo steps. This suggests IQ −IQ might be a more robust indicator of defect concentrations on short time scales, while cell-edge counting is simple and reliable for long time averages over uncorrelated frames. We also note that in systems of triangles and squares, defects are even more delocalized and Voronoi cell-edge counting is unable to identify defects, while IQ − IQ can still identify areas where defects are present.

Particle alignment
At all densities we simulate for triangles, squares, and hexagons, the body order parameter ξ s closely matches that of the bond order parameter ψ k,p . We quantify this by computing the cross-correlation between the two quantities | ψ k,p ξ * s | in the fluid, hexatic, and solid phases (Table I). We expect high bond-body cross-correlation at high density when particle rotation becomes locked in, as found in previous studies [10]. However, a strong signal is present even in the fluid for triangles, squares, and hexagons.
The cross-correlation results demonstrate alignment of the edges of neighboring particles. This can be interpreted as an effect of directional entropic forces [42]. We quantify the alignment with the potential of mean force and torque (PMFT) [43], computed from runs in the highest density pure fluid phase for each polygon. For each particle, we determine the relative position of its neighbors and construct histograms using many frames to obtain reliable averages. The free energy F of a given configuration is related to the negative logarithm of the histogram. We average an area in the PMFT 25σ from the center to set F = 0 as a baseline. Figure 6 shows the results of this calculation. All lowsymmetry regular polygons (n < 7) have distinctly separated wells for edge-to-edge contacts at F = −1.5k B T . This indicates preferential edge-to-edge alignment in the solid by strong directional entropic forces. But only when the symmetry of the particle shape is compatible with the symmetry of the bond order, as found in hexagons, triangles, and squares, does such edge-to-edge alignment promote local solid motifs in the fluid. Such a pre-ordered fluid allows the correlation length to increase smoothly from short-range to quasi-long-range and the transition to the solid to be continuous.
Triangles and squares have a distinct overall behavior from hexagons due to the delocalized defects. They have a smeared out phase transition between φ = 0.71 and φ ≈ 0.8 (Figure 7). Their equations of state are smooth and slowly increasing with barely detectable kinks (Figures S1 and S2). In comparison, the hexagon equation of state has sharp kinks and almost levels off through the transition ( Figure S5). Interestingly the melting transi- tion of hexagons becomes more evidently continuous as the system size increases and could, in fact, be mistaken for a first order transition in small systems of only a few thousand particles.
Our data agree qualitatively with previous studies. Simulations of equilateral triangles using an approximate event-driven molecular dynamics method showed a continuous fluid-to-liquid crystal-like phase transition [25]. Monte Carlo simulations revisiting that system reported a fluid-to-hexatic transition at φ = 0.70 and a hexatic-tosolid transition at φ = 0.87 [16]. Ref. [16] finds a chiral phase at φ ≥ 0.89, which we do not observe as we focus on the melting behavior at lower density φ < 0.80. Although no experiments have yet been reported on hard, mathematically regular triangles, rounded triangular colloidal platelets exhibit a fluid-to-hexatic transition, but the order of the transition was not determined [13]. Ref. [13] also finds a phase with local chiral symmetry breaking that we do not observe for hard regular triangles at φ < 0.85. Monte Carlo simulations of squares show a continuous fluid-to-tetratic transition at φ = 0.7 [26]. Experiments on vibrated granular squares (LEGOs) find tetratic orientational order in the range φ = 0.70 to 0.74 [17]. We are not aware of any studies on systems of hard hexagons to compare with.

B. First-order fluid-to-solid
Our data clearly show that systems of regular pentagons have a first-order transition directly from the fluid to the solid. This is seen in Figure 3d-f at φ = 0.688. There is a stripe of low density (local density Φ = 0.676) next to a stripe of solid (Φ = 0.694) with quasi-longrange positional order. The two phase region starts at φ = 0.680 and the pure solid phase starts at φ = 0.694. The two phase region is coincident with a Mayer-Wood loop in the equation of state ( Figure S3). The symmetry of the pentagon is not compatible with hexagonal ordering in the solid, so even though it has strong entropic edge-edge bonds, the body-bond cross-correlation (Table I) is zero and the phase transition becomes firstorder.
The one-step melting process of pentagons also shows up in the defect counts Figure 4a,b. We observe sharp kinks just below the pure solid phase at φ = 0.692 for counts in clusters with Burgers and disclination charges. One-step melting is supported by this coincident increase in the number of free dislocations and free disclinations.
Our results are consistent with previous studies of pentagons. Monte Carlo simulations on small systems of pentagons showed a transition from a fluid to a hexagonal rotator crystal at φ = 0.68 [10]. The same fluid to rotator crystal transition was also reported with rounded colloidal pentagons at φ = 0.66 [44]. The data in that work Disk results (n → ∞) are from [12,24]. The label 5 * refers to the 4-fold pentille. The n = 3 solid is a honeycomb lattice with alternating triangle orientations, the n = 4 the solid is a square lattice, and the n ≥ 5 solids are all hexagonal.
suggests a possible hexatic phase, but is inconclusive due to small system sizes within the field of view of the camera. Our simulations indicate there should be no hexatic phase with zero rounding. Pentagon phase behavior has also been studied in systems of 5-fold symmetric molecules [45] and with vibrating shaking tables [3,46]. None of these previous studies, however, conclusively demonstrated the first order, one step nature of the melting transition in pentagons.
Of all the regular shapes we studied, only triangles, squares, and hexagons fill space. These are also the only three KTHNY-type shapes. To test if another space filling polygon behaves similarly, we conducted simulations of the 4-fold pentille [28], an irregular pentagon with two different edge lengths that tiles space with 4-fold symmetry. We find that the 4-fold pentille behaves like the regular pentagon, with a first-order fluid-to-solid transition and no intermediate hexatic, though the transition occurs at a higher pressure ( Figure S4). At high fluid densities, directional entropic forces ( Figure 6) are blurred by the edge lengths and the resulting ten-fold particle body order is not compatible with either the tiling or the hexagonal solid motifs, so the transition is first-order. This suggests the space-filling property itself is not the factor triggering KTHNY melting but instead the similarity between the local order in the dense fluid and the local order in the solid. C. Disk-like behavior Figure 7 summarizes the phases found for all regular polygons studied in this work and includes the hard disk phase diagram (n → ∞) for comparison. Disks have a first-order fluid-to-hexatic transition, a narrow region of stability for the hexatic phase, and a continuous hexaticto-solid transition. A similar smooth behavior with n is found in the latent heat of the first-order transitions (Figure 8). The coexistence pressure decreases monotonically as n increases, approaching that of disks. The 4-fold pentille is an outlier compared to the regular polygons with almost a factor of two higher coexistence pressure.
We executed simulations up to n = 14 in this work and find that all regular polygons with n ≥ 7 have first-order fluid-to-hexatic and continuous hexatic-to-solid phase transitions. Figure 3g-i illustrates this for a system of octagons. In the two-phase region we see a mixture of fluid and hexatic phases. This is the same behavior as seen for disks, the only difference is that the transition shifts to lower packing fraction as n decreases. This shift is expected from the increasing anisotropy and relative strength of directional entropic forces with decreasing n. At n = 7 the start of the phase transition is at φ = 0.680, increases to φ = 0.692 by n = 12, 14 then again increases to φ = 0.702 for disks. We observe a small jump in critical density from n = 14 to disks, despite the very close coexistence pressures. We surmise that all regular polygons with n ≥ 7 have a first-order fluid-to-hexatic transition followed by a continuous hexatic-to-solid transition, with the transition range shifting higher in density with larger n.
Of the studied polygons with n ≥ 7, only the dodecagon with n = 12 has a particle symmetry compatible with the bond order in the solid. Although for this polygon strong directional entropic forces have the potential to drive a continuous transition from the fluid to the hexatic and solid phase, the simulations show a clear first-order transition for dodecagons. There is no bodybond cross-correlation in the dodecagon fluid (Table I), but it does appear weakly in the hexatic phase. The relatively short edge lengths in the dodecagon lower the strength of edge-edge alignment to the point where it is not strong enough to encourage local hexagonal motifs in the fluid, and as a result the transition from bond orientation disorder to order becomes sharp. All regular polygons with n ≥ 7 have smeared out F = −1k B T contact wells ( Figure 6) that encircle the polygon, indicating weaker edge-edge alignment than triangles, squares, and hexagons, allowing a much more isotropic behavior. This is why all regular polygons n ≥ 7 share the same fluidsolid transition properties as disks.
Regular polygons with n ≥ 7 exhibit an intermediate hexatic phase with a very narrow region of stability and a first-order transition to the fluid. As a result, the defect counts (Figure 4d-k) for these polygons is not as clear. The Burgers-charged and disclination-charged counts for these polygons have sharper kinks than the hexagons, but they do not parallel each other as closely as in the pentagons. Despite this, the hexatic-to-solid transition follows KTHNY predictions. In particular, the sub-block scaling analysis for χ predicts the hexatic-to-solid transition density with high accuracy, as shown in Figure 3i where the scaling line for φ = 0.704 falls almost exactly on the dotted line predicted by theory.

III. CONCLUSION
In 1988, Strandburg wrote [19]: "For a decade now the nature of the two-dimensional melting transition has remained controversial." At that time the three fluidto-solid transition scenarios discussed here were considered, but in many cases the available evidence was inconclusive as to which system shows which scenario and whether all scenarios indeed occur. 2D simulations have come a long way since 1988. With the advancement of high-performance computing we can now obtain correlation functions with high enough precision and rely on additional analysis techniques like order parameter fields, sub-block scaling, and cell-edge counting. Together, these tools reliably identify x-atic phases and resolve the nature of the 2D melting transition. Today, we can confidently state that the three scenarios discussed already 30 years ago indeed occur and, in fact, can be observed in a single system, namely the hard polygon family studied in the present work.
As we have demonstrated, the polygon's shape symmetry with respect to the lattice of the solid phase together determine the melting scenario. Systems of triangles, squares, and hexagons follow the KTHNY scenario well. They promote strong directional entropic forces that preorder the fluid into symmetries that are compat-ible with the solid. Pentagons and plane-filling 4-fold pentilles have similar strong anistropic aligning forces in the fluid, but their symmetry is incompatible with the hexagonal order of their solid phase. As a result they display a clear one-step first-order melting of the solid to the fluid with no intermediate phase. Regular polygons with sufficiently many (n ≥ 7) edges have no preferential alignment. They show a close resemblance to hard disk behavior and exhibit a first order fluid-to-hexatic phase transition and a continuous hexatic-to-solid transition, which is the intermediate scenario between KTHNY and one-step first-order melting.
Our results show that some 2D particles exhibit a KTHNY melting scenario when local coupling is strong and a hard-disk melting scenario when local coupling is weak. Our findings agree with a recent study [18] of disks interacting via soft r −m potentials; for m > 6, disks are hard enough to exhibit hard-disk melting, whereas for m ≤ 6, they exhibit KTHNY melting. In our system, we correlate "softness" achieved with n ≤ 6 with strong anisotropy in the entropic force field. For polygons with n ≤ 6, the anisometry is sufficient to allow other particles to approach at distances well inside the corresponding circumcircle of the polygon. This "overlap" creates the equivalent of a soft effective interaction. Rounding polygon edges limits how close neighboring particles may approach, and e.g. systems of sufficiently rounded hexagons should demonstrate disk-like behavior -we leave such a study for future work. More generally, the melting transition depends strongly on the symmetry of local interactions, not just the strength, when one has anisotropic interactions. In the present study, this anisotropic coupling is provided by emergent entropic forces, but we predict that anisotropic coupling of any origin is sufficient to produce the variety of melting scenarios we have observed in polygons. All data analysis in this work was performed using Freud, an in-house Python-driven high-performance toolkit developed by the Glotzer Group. Matthew Spellings implemented the complex correlation function routine and the cubeellipse color map. Eric S. Harper implemented the PMFT analysis code and parallelized all modules. We thank Wenbo Shen for suggesting that we examine correlations in the body orientation order parameter, and Greg van Anders for a critical read of the manuscript.