A critical edge number revealed for phase stabilities of two-dimensional ball-stick polygons

Two-dimensional (2D) systems are one of the most attractive topics in condensed matter physics because of their fundamental and practical importance. (Quasi-) 2D materials such as graphene^1,2 and Moiré superlattice³ exhibit novel physical properties originated from the massless Dirac cones^4,5,6,7 and the corresponding topological effects^8,9, including superconductivity, semi-conductivity, and insulation. The study on confined 2D ice is heuristic for understanding a wide range of natural phenomena¹⁰, and the ice surface has already been used for controlling self-assembly processes¹¹ and promoting chemical reactions¹². To gain insights from numerous structures and physical properties, a central focus is to understand the phase behaviors of 2D molecular systems, including phase stabilities, crystalline morphologies, and associated phase transitions.

The phase behaviors in 2D systems are more complicated than in three-dimensional systems¹³. First proposed by Landau and Peierls and then proved by Mermin and Wagner^14,15, a 2D crystal has a quasi-long-range translational order as well as a long-range bond-orientational order at a finite temperature. Unlike the 3D case, these two orders do not have to vary synchronously in 2D, leading to the possible intermediate phase of hexatic (with a six-fold symmetry)¹⁶ or tetratic (with a four-fold symmetry)¹⁷, characterized by a quasi-long-range bond-orientational order along with a short-range translational order¹⁶, which has no analog in 3D. Correspondingly, there are three different melting scenarios in 2D: (1) the KTHNY theory^16,18,19,20 when crystal melts into hexatic and then into liquid via two continuous phase transitions; (2) the hard-disk-like behavior²¹ when crystal melts into hexatic via a continuous phase transition and then into liquid via a discontinuous phase transition; (3) the regular solid-liquid phase transition²² when crystal melts into liquid directly without the appearance of the hexatic phase. Most recently, the simulation on truncated rhombs suggested the possibility of the fourth melting scenario when crystal melts into hexatic via a discontinuous phase transition and then into liquid continuously²³.

Realistic systems contain complex intermolecular interactions and molecular shapes, which are known to have dramatic influence on the 2D melting scenario, exceeding the capability of analytical calculations²⁴. Thanks to the fast development of simulation tools and data-analysis methods, complex model systems are possible to be investigated extensively. Roughly speaking, interaction leads to the enthalpy effect while shape is directly related to the entropy effect. Therefore, complex phase behaviors in 2D material systems may be analyzed from the viewpoint of the competition between entropy and enthalpy.

The entropy-limited case was investigated by studying hard polygons with only the volume-repulsive interaction (see e.g., Refs. ^{17,23,25,26,27,28}). A comprehensive simulation study by Anderson et al.²⁵ provides a complete picture of the melting scenarios of hard regular polygons, suggesting that triangle, square, and hexagon follow the KTHNY theory, pentagon follows a simple solid-liquid transition, and polygons with 7 or more edges follow the hard-disk-like behavior. They have also simulated a space-filling pentagon to clarify that it is the orientational entropy and the confliction of geometrical symmetry between liquid and crystal, rather than the space-filling, that determine the melting scenario. The colloidal experiments on hard polygons have observed more solid structures, e.g., the rhombic crystal for squares²⁹, the frustrated glassy state for pentagons³⁰, and the rotator crystal phase for hexagons³¹, some of which can be observed in the simulations of corresponding round-corner polygons^31,32.

The enthalpy-limited case was studied by simulating various interactive point systems^{33,34,35,36,37,38} and performing experiments with confined monolayer colloidal spheres³⁹. It has been shown that the melting scenario can be tuned by varying the interaction parameters. One of the general conclusions is that there can exist different melting scenarios in the same system at different densities, which has been found in different interactive point-like particles and polygons composed of attractive Lennard-Jones (L-J) beads^36,38,40.

Only in recent years did researchers pay attention to more realistic interactive polygons considering both interaction and shape, which leads to qualitatively different phase behaviors. The simulations on L-J beads polygons⁴⁰ (square, pentagon, and hexagon) reveal that a wide solid-liquid coexistence region similar to the L-J point system emerges at low temperatures, while the melting scenario of corresponding hard polygon is followed at high temperatures. Different melting scenarios under various thermal conditions provide a basic understanding on the thermodynamic competition between two factors: Enthalpy plays a major role at low temperatures/pressures while entropy is dominant at high temperatures/pressures. Experiments on millimeter-sized polygons (triangle and square) with isotropic magnetic repulsion⁴¹ exhibit crystalline structures incompatible with the symmetry of monomer shape, e.g., the appearance of the triangular lattice structure for the square system, which results from the confliction between the isotropic interaction and the anisotropic shape.

Since these models are not suitable for modeling real molecular systems, the above results are still inadequate for understanding most 2D materials, in which novel phase behaviors may emerge due to the complex competition between entropy and enthalpy. The ball-stick polygons with an L-J ball on each vertex and adjacent balls connected by valence bonds are ideal for extending previous theoretical studies towards the understanding on real 2D materials: On the one hand, the ball-stick models mimic the molecular shape along with the intermolecular interaction; on the other hand, they can be directly compared with the well-studied model systems we described above. Because the L-J balls are on the vertices, although the anisotropies induced by shape and interaction of each molecule belong to the same symmetry group, they spatially conflict with each other, indicating a more complex competition between entropy and enthalpy, which is closer to realistic systems. Since the confliction is weaker for polygons with more edges where both the shape and interaction are more isotropic, the phase behaviors of polygons with fewer edges are expected to deviate more from the corresponding hard polygons. Moreover, no previous studies have compared phase stabilities between different polygons, so some fundamental questions, such as which physical factor is dominant for phase stabilities under different thermal conditions and how it relates to monomer properties, remain unanswered. Here and below, both monomer and molecule refer to one ball-stick polygon.

In this work, we investigate phase stabilities of 2D regular ball-stick polygons (triangle, square, pentagon, hexagon, and octagon) by molecular dynamics (MD) simulation supplemented with some numerical calculation methods. As summarized in Fig. 1, with increasing edge number, we discover a critical edge number \({n}_{{\mathrm{c}}}=4\) for the phase stabilities in this polygonal family below melting points: The triangular system is stabilized in a spin-ice-like glassy state rather than a crystalline structure at any finite temperatures before melting into liquid via a discontinuous phase transition; the square system exhibits a distorted square lattice structure, featured by a self-similar translational symmetry with a high concentration of topological defects; the pentagon, hexagon, and octagon systems are stabilized with the normal crystalline structures analogous to corresponding hard polygons. Moreover, for the crystalline state, \({T}_{{{{\rm{m}}}}}\) is higher for polygons with more edges at a higher pressure but exhibits a crossover for hexagon and octagon at a lower pressure. Detailed analysis suggests that the perfect match of interaction and shape in the hexagon system induces an ultra-stability at low pressures, revealing the entropy-dominant mechanism at a higher pressure and the enthalpy-dominant mechanism at a lower pressure.

As n increases, there exists a critical edge number \({n}_{{\mathrm{c}}}=4\) for the stability of the crystalline morphology.

Full size image

Triangle: ground-state degeneracy

Theoretical determination of ground states

The presence of the short-range L-J interaction results in two isotopic unit-cell dimer structures with exactly the same potential energy, splitting from the original one for the dense packing scheme of the hard regular triangle, as illustrated in Fig. 2a. Either type of the isotopes can be spatially duplicated to form the crystalline morphology in Fig. 2b, in which each unit-cell is composed of two neighboring monomers colored differently according to their body-orientations, as shown in the upper part of Fig. 2c. The spatial arrangement of unit-cells is shown in the lower part of Fig. 2c, which is a rhombic crystal phase^29,32 (\({\alpha }_{1}\ne {\alpha }_{2}\), \({\alpha }_{2}\) is smaller than \(\pi /2\) but larger than \(\pi /3\)) different from the crystalline structure of the hard triangular system²⁵ due to the existence of the L-J balls.

a Two isotopic structures of the unit-cell, splitting from the original space-filling scheme for the hard triangle. b Crystalline morphology. c Spatial arrangement of unit-cells. One of \({\alpha }_{1}\) and \({\alpha }_{2}\), which are the angles formed by two unit-cell primitive vectors, is larger than \(\pi /2\), corresponding to a rhombic crystal structure. Each point represents a unit-cell composed of two monomers with opposite body-orientations, as shown at the top of this picture. d Three stable local structures composed of the two isotopes with the local clusters in the three configurations all composed of six atoms, which are circled by red, purple, and green frames, respectively. e Three typical small hand-made structures constructed based on the three local structures in d. The top one is a part of crystal, noted as ‘crystal’, the middle one is a mixture of the first two local configurations in d, noted as ‘perfect-mixing’, and the bottom one is formed by mixing all the three local configurations in d, noted as ‘tower-like’. The local clusters are framed with the same color as in d. f Potential energies (\({E}_{{{{\rm{p}}}}}\)) of different structures under various cutoff radii (\({r}_{{{{\rm{c}}}}}\)) with the solid lines representing \({E}_{{{{\rm{p}}}},{{{\rm{lowest}}}}}\) and the dashed lines representing \({E}_{{{{\rm{p}}}},{{{\rm{avg}}}}}\). Physical quantities are expressed in L-J units.

Full size image

These two isotopes can also mix with each other arbitrarily, forming an amorphous state without any geometrical conflicts. In order to have a comprehensive understanding on all the global structures, it is necessary to enumerate the possible local configurations at first. As shown in Fig. 2d, the stacking of the two types of isotopes can form three non-equivalent local clusters. The exhaustive geometrical enumeration (details in SI) suggests that the first two local configurations in Fig. 2d can fill in the whole space by their own to construct a crystalline structure or by mixing with the other one to form an amorphous state, while the third one cannot fill in the whole space by itself and must mix with the first two to avoid the creation of geometrical defects.

If we count in only the L-J interaction from nearest-neighboring atoms, all the global structures without geometrical conflicts, including both crystalline and amorphous ones, have exactly the same potential energy when the system is so large that the energy cost on the boundary can be ignored, attributed to the fact that the atom clusters in all three of the above stable local structures, which are framed respectively by red, purple, and green in Fig. 2d, have 9 atom pairs. Therefore, a huge number of possible global structures with identical potential energy are all the ground states at \(T=0\), i.e., the triangular system has a ground-state degeneracy, which is very similar to the ground-state degeneracy in the spin-ice system⁴², where the oxygen atoms adopt a fixed lattice arrangement while hydrogen atoms can be either ordered or disordered without altering the potential energy.

When the interaction exceeds nearest neighbors, the potential energy is difficult to be analytically calculated, so we instead calculate it numerically. We first construct three typical structures with 30 monomers, as shown in Fig. 2e, where the top one composed of solely the first local structure in Fig. 2d, the middle one containing both the first and the second structures, and the bottom one formed by mixing all the three local structures together. For comparison, we calculate two quantities related to the potential energy \({E}_{{{{\rm{p}}}}}\) for the three structures in Fig. 2e: (1) \({E}_{{{{\rm{p}}}},{{{\rm{avg}}}}}\), the potential energy averaged over nine clusters, as circled in each structure; (2) \({E}_{{{{\rm{p}}}},{{{\rm{lowest}}}}}\), the lowest potential energy of the clusters. The first one quantifies the global property while the second one provides the information on the local stability. By combining these two quantities, we can have a good estimation on the ground state(s). These two quantities are calculated with the cutoff radius (\({r}_{{{{\rm{c}}}}}\)) ranging from \(3\sigma\) which is approximately the size of the atom cluster, to \(7\sigma\) at which the L-J interaction has a very small value about \({10}^{-5}\), where \(\sigma\) is the van der Waals (vdW) diameter, and the results are shown in Fig. 2f. These three structures have nearly the same potential energy, as the difference of the potential energy varies with \({r}_{{{{\rm{c}}}}}\) but remains always on the order of \({10}^{-4}\). Therefore, we affirm that the ground state in the triangular system is degenerate, with two crystalline morphologies purely formed by one of the isotopes and a huge amount of amorphous configurations formed by mixing different local structures in Fig. 2d. At any finite temperatures, the vast amount of possible ways to mix the isotopes results in a very large configurational entropy, which allows the amorphous state to have a lower free energy than the crystal state.

Confirmation of the ground-state degeneracy by simulation

To verify the above structural origin of the crystalline instability, we performed a replica exchange MD (REMD) simulation on 625 ball-stick triangles with 26 temperature replicas ranging from 2.55 to 6.1. An energy minimization on the final structure of the replica at the lowest simulated temperature has been performed after the REMD simulation, which is theoretically the most stable configuration at low temperatures⁴³. All the three local structures in Fig. 2d appear in the optimized structure, whose representative snapshot is shown in Fig. 3a. Moreover, benefiting from the fact that the bond length (1.5\(\sigma\)) is just a little larger than the equilibrium distance of \(\root 6\of{2}\sigma\), at finite temperatures, there exist some ‘meta-space-filling’ structures (colored in blue), where at least two atoms in the same cluster are connected by bond directly, leading to small geometrical defects and a little higher potential energy than the space-filling ones at zero temperature. As a result, they should vanish at \(T=0\) but could appear at a slightly higher temperature for they further enlarge the configurational entropy.

a A small part of the representative snapshot for the most stable structure determined by the replica exchange molecular dynamics simulation. The green and red frames mark different local configurations in Fig. 2d, and the blue frames illustrate the existence of meta-space-filling structures at finite temperatures. Red and green have the same meaning as in Fig. 2d, while we do not distinguish red and purple here. b Amorphization process characterized by the decay of bond-orientational order parameter \(|{\varPsi }_{4}|\) at three temperatures below the melting point. The inset is a small part of the snapshot at T = 2.86 after amorphization. c Caloric curve (potential energy \({E}_{{{{\rm{p}}}}}\) versus temperature T) and density curve (density \(\rho\) versus T), both with the discontinuous changes located at T = 2.88‒2.89. The error-bars represent standard deviations, resulting from the thermal fluctuations in one simulation trajectory. d Time evolution of the unit-cell order parameter d in the simulation starting from a locally disordered crystalline structure at T = 2.83, 2.86, and 2.88. Each curve plotted in b and d is based on the data from a single simulation. Due to the randomness of the dynamical process, the starting point and the rate of amorphization may take different values in different trajectories. Physical quantities are expressed in L-J units.

Full size image

As the amorphous state is more stable than crystal at finite temperatures, we expect that, driven by entropy, the crystal state will spontaneously destabilize into the amorphous state. To check this, MD simulations starting from crystalline morphology are performed at pressure \(P=0\) and various temperatures. As quantified by the bond-orientational order parameter \({\varPsi }_{4}\) (see Methods for the definition) plotted in Fig. 3b, which initially keeps a high value (larger than 0.8) corresponding to the rhombic crystal state^29,32 and then decays after certain simulation steps until a very low value due to its amorphous nature, the amorphization takes place at \(T\ge 2.85\), lower than the melting point \(T=2.88 \sim 2.89\) manifested by the sharp changes in the caloric and density curves plotted in Fig. 3c. The non-equilibrium amorphous state whose potential energy and density are close to those of crystal is featured by the mixing of two isotopes, as shown in the inset of Fig. 3b, in agreement with our theoretical prediction described above. At lower temperatures, the fact that amorphization does not occur within the finite simulation time should be attributed to the existence of relatively high energy barrier(s) between the crystalline and amorphous states. To verify this, we manually introduce a very small disordered region (about 30 monomers) into the crystalline configuration (see Methods for technical details) and run NPT simulations at \(T=2.83\), 2.86, and 2.88, respectively. The evolution of the two monomers in each unit-cell is traced by the unit-cell order parameter defined as \(d(t)=\frac{1}{n} < {\sum}_{{{{\rm{unit}}}}-{{{\rm{cell}}}}}\varTheta ({l}_{{{{\rm{c}}}}}-l(t)) > \), where the number of unit-cells n equals to half of the number of monomers, l(t) is the distance between the center-of-masses (COMs) of the two monomers originally in the same unit-cell at time t, \({l}_{{\mathrm{c}}}\) =5.77 is the cut-off value of l, and \(\varTheta\) is the Heaviside step function. This order parameter quantifies the fraction of remaining unit-cells at time t, taking a value close to 1 in crystal and 0 in liquid. The time evolution of d in \(6.4\times {10}^{7}\) steps after a relaxation of \(4\times {10}^{6}\) steps is shown in Fig. 3d, demonstrating that, once the energy barrier is overcome, the system will rearrange into the amorphous state continuously. We expect that d would finally drop to approximately zero after sufficiently long time, as a result of reorganization of the whole system. Although regular MD simulations at much lower temperatures would suffer from very slow dynamics, we have managed to verify that the amorphous state is still more stable than the crystalline state even at a temperature as low as T = 1.5 by running an additional REMD simulation focusing on the low-temperature region (see Supplementary Fig. S4 in SI). Moreover, we never observe an inverse pathway from the amorphous state back to the crystal, as appeared in the simulation for the hard-triangle system²⁵, in agreement with the fact that the amorphous state is more stable than the crystal benefiting from the large configurational entropy.

Glassy nature of the spin-ice-like amorphous state

Detailed analysis on dynamics provides a more comprehensive understanding on this amorphous state. As shown in Fig. 4a, in the amorphous state at \(T=2.83\), 2.86, and 2.88, the calculated mean-squared displacements (MSDs) \({\varDelta }^{2}\) of monomers grow in a strongly frustrated way, clearly distinguished from the one for crystal or liquid shown in Supplementary Fig. S5. We further plot the non-Gaussian parameter for MSD in Fig. 4b, defined as \(\alpha=\frac{3 < {(r(t)-r(0))}^{4} > }{5 < {(r(t)-r(0))}^{2}{ > }^{2}}-1\), which is expected to be 0 for the Brownian motion and larger than 0 when some monomers move faster than others. The calculated \(\alpha\) for this amorphous state is higher than in the crystal or liquid state, indicating a strong dynamical heterogeneity and the fact that the amorphous state is a glassy state from the dynamical perspective⁴⁴.

a Mean-squared displacements (MSDs) (Δ²) in the glassy state exhibiting frustrated diffusive behavior. b Non-Gaussian parameters for MSD \(\alpha\) vs. temperature T. The two values at T = 2.83 correspond to the crystalline state (smaller one) and the glassy state (larger one). c Schematic illustration of the spin-ice-like state defined on a topologically triangular lattice with the ‘4-in-2-out’ ice-rule. The left panel shows the six neighboring atom-clusters (circled in blue) of the central cluster (circled in red) while the right one illustrates that the six atoms in the cluster exhibit the ‘4-in-2-out’ feature along the direction indicated by the orange arrows. d Bond-orientational order parameter \({\varPsi }_{6,{\mathrm{atom}}}\) vs. T, whose discontinuous change suggests a first-order nature of the melting of the spin-ice-like state. The inset is a small portion of the snapshot exhibiting the local packing scheme, where the purple lines and circles highlight the six neighboring atoms of the central atom shadowed in purple arranged in a nearly dense-packing way. The error-bars in this figure represent standard deviations, resulting from the thermal fluctuations in one simulation trajectory. Physical quantities are expressed in L-J units.

Full size image

Overall, considering both the structural and dynamical characteristics, the 2D ball-stick triangular system is in a spin-ice-like glassy state at any finite temperatures below melting point⁴⁵, as a result of the strong confliction between interaction and shape. As shown in Fig. 4c, each stable local atom cluster has six neighboring clusters and contains six atoms with four close to each other while the other two away from them, so our system shares some common features with a spin-ice system defined on a triangular lattice with the ‘4-in-2-out’ ice rule⁴⁶. However, the ball-stick triangular system is more complicated than the spin-ice system due to the following two features: (1) Not only the pair interaction from nearest neighbors is considered; (2) The triangular lattice is defined topologically without well-defined lattice structure (no on-site oxygen atoms), so the spatial structure is totally disordered.

The detailed examination on Fig. 3a, b suggests that no matter the isotopes mix or not, each atom has six atomic nearest neighbors arranging in a nearly dense-packing way. This local structure is shown clearer in the inset of Fig. 4d, where the atom shadowed in purple has six neighbors circled by purple frames. Inspired by the spin-glass system^47,48 whose magnetization order parameter characterizes the order-disorder transition while the Edwards-Anderson order parameter measures the glass-liquid transition, here it is also possible to find an order parameter other than \({\varPsi }_{4}\) to distinguish the glassy state from the liquid one. We therefore employ the bond-orientational order parameter \({\varPsi }_{6,{{{\rm{atom}}}}}\) associated with all atoms indistinguishably of the presence of bonds instead of molecular COMs. The discontinuous drop of \({\varPsi }_{6,{{{\rm{atom}}}}}\) at T = 2.88~2.89 in Fig. 4d demonstrates that the spin-ice-like state melts into liquid via a first-order phase transition, in agreement with the caloric and density curves in Fig. 3c.

Furthermore, the MD simulation at P = 10 (details in SI) manifests that the crystalline morphology rearranges into a spin-ice-like state at finite temperatures and the melting transition from spin-ice-like to liquid is still a first-order phase transition, indicating that the revealed phase behavior of the 2D ball-stick triangular system is robust with respect to different thermal conditions.

Square: distorted square lattice

Based on the result from REMD simulation and the consideration on symmetry, the initial state was established as a square lattice shown in Fig. 5a, where the body-orientation of each monomer is not perpendicular to the primary vector of lattice to ensure a higher packing fraction. MD simulations at various temperatures were performed at P = 0. The apparent discontinuous changes in the caloric curve, density curve, and \({\varPsi }_{4}-T\) plotted in Fig. 5b manifest a first-order phase transition from solid to liquid with the melting point located between 2.92 and 2.93.

a Initial crystalline morphology. The two orange arrows are the primitive lattice vectors. b Caloric curve (potential energy \({E}_{{{{\rm{p}}}}}\) versus temperature T) and density curve (density \(\rho\) versus T). The inset is the bond-orientational order parameter \({\varPsi }_{4}\) at different temperatures. All the discontinuous changes of these three curves locate at T = 2.92~2.93. The error-bars represent standard deviations, resulting from the thermal fluctuations in one simulation trajectory. c The spatial arrangement of center-of-masses (the black dots) at T = 2.92. The polygon corresponding to each point is decided by the Delauney tessellation considering the periodic boundary condition and is colored according to its number of edges. d Bond-orientational correlation function (correlation \({{{{\rm{C}}}}}_{\psi }\) as a function of distance r) at the melting point for N = 2496 monomers and N = 22,500 monomers, indicating no decay compared to the brown line corresponding to an algebraical decaying behavior. e Translational correlation functions (correlation \({{{{\rm{C}}}}}_{g}\) as a function of r) at melting points for N = 22,500 monomers with different block sizes n (number of monomers in each block), decaying faster than algebraically (the black line) locally but slower on the long range. The distance is rescaled by the characteristic length of the coarse-graining a. Physical quantities are expressed in L-J units.

Full size image

The visual examination of Fig. 5c tells us that, in the solid state, although each monomer has four neighbors aligned nearly orthogonally, also reflected by the diffraction pattern shown in Supplementary Fig. S7a, it apparently distorts from the normal square lattice. This is evidenced by the fact that the Voronoi cells are close to squares but mostly have 6 edges. Moreover, the concentration of topological defects, defined as the fraction of Voronoi cells with more or less than 6 edges, is as high as about 0.42 (data in Supplementary Fig. S8b), significantly larger than the typical value of crystal (about 0.05) and even hexatic (about 0.123)^25,49,50. The quantitative characterization of this state exhibits two seemingly contradict facts: (1) The long-range bond-orientational correlation can be retained until melting, indicated by the correlation function in Fig. 5d, as well as the heatmap in Supplementary Fig. S7c colored according to the argument of \({\psi }_{4,i}\), where the whole picture holds almost the same color, in agreement with its very narrow distribution; (2) The translational correlation function exhibits a decaying behavior slightly faster than algebraical, as shown in Supplementary Fig. S8e.

To clarify this paradox, we perform simulations on a larger system with N = 22,500 monomers, whose results are consistent with the above (details in SI). With the larger simulation scale, we identify that, the translational correlation decays relatively fast locally but exhibits a typical algebraical decaying behavior in the long range, as illustrated by the red line in Fig. 5e, indicating that this solid state is still a crystal. Combined with the above features, we regard it as a distorted square lattice. The translational symmetry is further investigated by a ‘renormalization’ process, in which we divide the monomers into blocks, with n (=1, 4, 9, 25, 100, 225) monomers in each block, and then calculate the translational correlation of the block COMs. After rescaling the distance by the characteristic length of the block, which is just the distance of the first peak in the correlation function, different curves collapse together, as shown in Fig. 5e, exhibiting the so-called self-similarity, which guarantees the existence of the translational symmetry. Furthermore, considering the locally disordered nature, this self-similarity strongly implies that it is a critical point for the ordered and disordered solid states^51,52, neither totally disordered as the triangular system nor perfectly ordered as a standard crystalline structure for pentagon, hexagon, or octagon. Since this holds at different temperatures (e.g., also T = 2.8 in Supplementary Fig. S8f), the critical point is defined in the shape space in terms of a critical edge number, rather than in the thermal-parameter space.

When cooling down the distorted square lattice, the Voronoi cells have less dispersive areas but surprisingly higher concentrations of topological defects (Supplementary Fig. S8b, c), even for the case as low as T = 0.01. This supports the mechanism that the local distortions are originated from the intrinsic properties of the system, in terms of interaction and shape rather than thermal fluctuations, emphasizing its critical nature in the shape space.

We also perform simulations at P = 5 and P = 10. Under different pressures, the square system still follows a solid-liquid transition without the appearance of the tetratic phase and the solid structures also share similar features (see SI).

Pentagon, hexagon, and octagon: entropy-dominant versus enthalpy-dominant

Phase diagrams

For these three polygons, the ground-state crystalline structures: striped phase for pentagon, triangular lattice for both hexagon and octagon, can be stabilized until their melting temperatures, attributed to the even weaker confliction between interaction and shape. The snapshots of these crystalline configurations are shown in Fig. 6a, b, and g, respectively. Below we describe the melting scenarios of each polygon under various pressures, and the readers are referred to the Methods section and SI for more details on the calculation procedure and results.

a Striped phase of the pentagon system, in which monomers have two body-orientations comprising alternating lines. b Rotator crystal phase of the pentagon system, in which monomers have random body-orientations. c Body-orientation distributions for the two crystal phases of the pentagon system, where \(\theta\) represents the orientation of monomer with respect to the x-axis of the simulation box, ranging from 0 to \(2\pi /n\) for n-gons, and \(P(\theta )\) is the corresponding probability density. d Triangular lattice crystal phase of the hexagon system. e Rotator crystal phase of the hexagon system, in which monomers have two preferred body-orientations. f Body-orientation distributions for the two crystal phases of the hexagon system. g Triangular lattice crystal phase of the octagon system. h Rotator crystal phase of the octagon system. i Body-orientation distribution for the two crystal phases of the octagon system, where the y axes are scaled differently to exhibit the trimodal distribution in the rotator crystal phase.

Full size image

Our MD simulations indicate that the striped phase is the most stable phase at low temperatures for the ball-stick pentagon system, which is also theoretically the closest packing mode for hard pentagon⁵³. With increasing temperature, the striped phase first experiences a first-order solid-solid phase transition into the rotator crystal phase and then melts into the liquid phase via another first-order phase transition, the same as the melting scenario of hard pentagon²⁸. In the striped phase, monomers are aligned with two body-orientations to form alternative lines, as shown in Fig. 6a, rendered as two peaks in the body-orientation distribution (the blue line in Fig. 6c). In the rotator crystal phase shown in Fig. 6b, monomers have a totally random body-orientation, demonstrated by the flat distribution of body-orientation (the orange line in Fig. 6c), indicating that this rotator crystal phase is a continuous one⁵⁴. The phase behavior of the pentagon system is summarized in the phase diagram in Fig. 7a, in which the striped phase, rotator crystal phase, and liquid phase appear in sequence as temperature increases at any pressures.

a Pentagon. b Hexagon. c Octagon. The striped phase (S) and the triangular lattice crystal phase (TX) are colored in orange, the rotator crystal phase (RX) is colored in light blue, the hexatic phase (H) is colored in black (very narrow for hexagon and octagon), and the fluid phase (F) (liquid or gas) is colored in pink. Physical quantities are expressed in L-J units.

Full size image

The hexagon system exhibits a rich phase behavior with increasing pressure. At P = 0, the crystalline structure sublimates without the appearance of the liquid state, evidenced by a nearly-zero potential energy and an RDF curve with only one peak shown in Supplementary Fig. S11 in SI. At P = 0.5 and 1.0, the crystalline structure melts into the liquid state via a first-order phase transition. At P = 1.5, 2.5, and 5, a rotator crystal phase (Fig. 6e) appears between the crystal state and the liquid state, and both of the two phase transitions are discontinuous. The monomer COMs in the rotator crystal phase form a triangular lattice, but monomers have two body-orientations identified by the two broad peaks in Fig. 6f, so this rotator crystal phase is regarded as a discontinuous one⁵⁴. When the pressure increases to P = 7.5, the hexatic phase appears between the rotator crystal phase and the liquid phase. In this case, both transitions from the triangular lattice crystal to the rotator crystal and from the hexatic phase to the liquid phase are discontinuous, while the one from the rotator crystal phase to the hexatic phase is continuous, basically following the hard-disk-like behavior²¹. Even at P = 10, the melting scenario still follows a hard-disk-like behavior, instead of the KTHNY theory followed by the hard hexagon or the L-J beads hexagon^25,40, attributed to the appearance of the rotator crystal phase. The phase diagram is shown in Fig. 7b, in which the rotator crystal phase only appears at relatively high pressures and the hexatic phase emerges at even higher pressures.

The octagon system stabilizes at the triangular lattice up to a finite temperature when it transforms into a rotator crystal phase discontinuously. The rotator crystal of octagon has roughly three preferred body-orientations of monomers, as shown in Fig. 6g, h, similar to the case for hard octagon⁵⁴. However, the body-orientation distribution is very rough and the ranges between peaks and valleys are much smaller than that of the discrete rotator crystal phase for hexagon. The melting scenario follows a regular solid-liquid one at relatively low pressures (P = 0–5). The hexatic phase appears at higher pressures of P = 5–8 and leads to a hard-disk-like behavior, qualitatively the same as the L-J disk⁴⁰. The phase diagram of the octagon system is drawn in Fig. 7c, showing that there are always two solid states, and the hexatic region is very narrow at high pressures and disappears at low pressures.

Phase stabilities

To gain a deeper understanding on the mechanism of phase stabilities, we plot the melting points of different crystalline structures under different pressures in Fig. 8. The melting point increases with pressure and is higher for polygons with more edges at relatively high pressures, which is regarded as the entropy-dominant regime. Meanwhile, the melting-temperature curves for hexagon and octagon have a crossover point located between P = 1 and P = 1.5. At relatively low pressures (\(P\le 1\)), hexagon stays in the triangular-lattice crystal phase before melting into liquid, and has a higher melting point than octagon. By contrast, at relatively high pressures (\(P\ge 1.5\)), hexagon is in the rotator crystal phase before melting into liquid, and has a melting point lower than octagon. Therefore, the crossover of the melting point locates in the pressure interval where the rotator crystal phase for hexagon appears. Detailed examination on the local structure of triangular-lattice crystal morphology suggests that, as shown in Fig. 6d and much clearer in Fig. 9d, each atom has two neighboring atoms belonging to two different molecules, and these three atoms form a regular triangular local structure, which is the most stable structure of an L-J cluster composed of three atoms. In other words, the crystal phase for hexagon has the intermolecular interaction perfectly matches the molecular shape, leading to an ultra-stability at low pressures when the phase stability is dominated by enthalpy.

Lines in different colors represent different polygons, and different symbols represent different solid states: square for close packing (CP) and triangle for rotator crystal (RX). Physical quantities are expressed in L-J units.

Full size image

a Triangle, in which two monomers colored differently belong to the same unit-cell with different orientations. b Square. c Pentagon. d Hexagon. For each polygon, there is a non-zero torsional angle between the connecting line of neighboring center-of-masses determined by shape only (i.e., dense-packing pattern, represented by the dark purple line) and the one in the ground state of ball-stick polygon (represented by the orange line).

Full size image

The above phenomenon helps us to deepen our understanding on the competition between the entropy and enthalpy. At relatively low pressures, the interaction plays a leading role in phase stability, i.e., whether the local structure is interaction-favored determines the phase stability, which is regarded as ‘enthalpy-dominant’. Conversely, at relatively high pressures, the molecular shape of a polygon determines its phase stability, corresponding to an ‘entropy-dominant’ mechanism. The striped phase for pentagon and the triangular lattice crystal for octagon do not show the same ultra-stability as hexagon because the interaction and shape there do not perfectly match each other.

In this work, we have investigated phase stabilities of the 2D ball-stick polygonal family (triangle, square, pentagon, hexagon, and octagon) by MD simulation and numerical calculation. We have found a critical edge number of \({n}_{{{{\rm{c}}}}}=4\) for the crystalline stability of 2D ball-stick polygons. For \(n\; < \;4\), the ball-stick triangular system has no stable crystalline structure at any finite temperatures, instead it forms a spin-ice-like state characterized by amorphous structures and glassy-like dynamics due to the ground-state degeneracy originated from the mixing of two types of isotopic dimer structures. For \(n=4\), the ball-stick square system forms a distorted square lattice with a high concentration of topological defects but self-similar translational symmetry below melting. For \(n\; > \;4\), the ground-state crystal structures: striped phase for pentagon, triangular lattice for hexagon and octagon, can be stabilized up to the melting temperature. We have also studied the phase stabilities of the crystal states of different polygons by investigating the P‒T_m curves, and found a melting-point crossover for hexagon and octagon.

All the results in this work may be understood under the theoretical framework based on the competition between entropy and enthalpy. The 2D ball-stick polygons in this study exhibit spatially conflicted anisotropies induced by the molecular shape (entropy) and intermolecular interaction (enthalpy), resulting in a deviation of the stable structures from the space-filling pattern of the corresponding hard polygons, as shown in Fig. 9, in which the dark purple lines represent the connection of COMs of neighboring polygons determined solely by geometrical consideration, while orange lines represent the case in the ball-stick polygon systems. The critical edge number for the anisotropic confliction is 4: When n < 4, the confliction is strong enough to break the crystalline morphology and leads to a spin-ice-like glassy state; when n > 4, the confliction is weakened as a perturbation because the shape and interaction of monomers become more isotropic, resulting in phase behaviors close to the corresponding hard polygons; when n = 4, the effects of interaction and shape are subtly balanced with each other, served as a critical point in the shape space. For triangle, due to the vanish of ordered structures, no hexatic phase appeared in hard triangular system is observed²⁵. For square, the distorted lattice directly melts into liquid without going through the tetratic phase appeared in the hard square^17,25 or L-J beads square⁴⁰ system, attributed to the sufficiently large concentration of topological defects in the solid state. For polygons with more than 4 edges, the appearance (or disappearance) of a hexatic phase is basically the same as corresponding hard polygons²⁵. Moreover, under various thermal conditions, two different mechanisms for phase stabilities compete with each other: The entropy-dominant mechanism at relatively high pressures results in a higher melting point for polygons with more edges due to the smaller orientational entropy, while the enthalpy-dominant mechanism at relatively low pressures takes interaction details into account, inducing an ultra-stability for hexagon whose crystalline structure has an interaction-favored local structure.

“More is different”⁵⁵. People have accumulated extensive knowledge on both the entropy-limited case and the enthalpy-limited case for 2D systems, as well as the thermodynamic and spatial competitions between these two factors. Real 2D systems, however, stand in between, involving comparable entropy and enthalpy effects, whose phase behaviors are more abundant than either of the two above cases. This work reveals some interesting and unique physical features in 2D molecular systems, as well as establishes a theoretical framework which can be generalized to understand phase behaviors in more complex 2D systems. Besides its fundamental importance, our research is also valuable for the design of 2D materials at least in two aspects: (1) When there is a strong anisotropy induced by molecular shape or intermolecular interaction, the details of short-range interactions are crucial for designing various structures, such as amorphous materials used for frustrated magnetic materials⁴⁵, and distorted lattice which is generally existed in some compound systems, e.g., high-entropy alloys⁵⁶, widely used for tuning the mechanical properties⁵⁷ or the electron band structures⁵⁸; (2) Phase stability at high pressures is mainly determined by shape, but at low pressures it can be tamed by the match of interaction and shape, which can facilitate the design of materials with desired physical features, such as ultra-stability and heat-sensitivity.

In our simulation, the deformation or elongation of the polygon can be ignored due to the large strengths of the bond and angle interactions (see Methods for more details). Since the mechanism we have provided above does not explicitly depend on bond length, we argue that our conclusion is still effective for the normal case of molecules. However, many studies have shown exotic melting scenarios and ground-state configurations when the polygon is allowed to deform^59,60,61. Moreover, the melting scenarios of polygon mixtures may differ from identical polygon systems attributed to the much more complicated competition among various ingredients⁵⁰. Considerations on these two factors may lead to a comprehensive understanding in the shape space and reveal the critical behavior in the vicinity of \({n}_{{{{\rm{c}}}}}\), but are beyond the scope of our study. In real applications, 2D materials have to be supported on a substrate. Although many substrates are known to have little influences on crystalline structures^16,20, some well-designed substrates may be used for taming melting scenarios⁶². The consideration on the effect of substrates might result in a more complex competition between entropy and enthalpy, leading to the emergence of new structures, and providing more guidance on experiments and industrial applications of 2D materials.

Preparation of initial configurations

Molecules are modeled by regular polygons with a 12-6 L-J particle on each vertex, and each intra-molecular pair of adjacent vertices are connected by a covalent bond. Each n-gon monomer is composed of n atoms and n covalent bonds. We use L-J units in this work, where physical quantities are expressed in a dimensionless manner, and set L-J parameters m = 1, \(\sigma=1\), and \(\varepsilon=5\). Both the covalent bonds and the valence angles are represented by a strong harmonic potential with k = 1800\(\varepsilon /{\sigma }^{2}\). The equilibrium bond length in each polygon is set as short as 1.5\(\sigma\) to prevent the overlap of two monomers or two bonded atoms (A typical elongation of \(\varDelta l=0.1\sigma\) would correspond to a characteristic temperature equal to \(T=18\varepsilon=90\), significantly larger than the highest temperature we have used in all the simulations about T = 20). We set the cutoff radius of the L-J interaction to \({r}_{{{{\rm{c}}}}}=7\sigma\) in the simulation, at which the L-J potential has decayed to a very low value (less than \({10}^{-4}\)), ensuring a high accuracy on the estimation of potential energy. Moreover, we have performed several simulations with various values of \({r}_{{{{\rm{c}}}}}\), from 3\(\sigma\) to 10\(\sigma\). A \({r}_{{{{\rm{c}}}}}\) value out of this range is unnecessary to be tested because a smaller attraction range would result in a large error while a larger one should have no influences since the potential energy at \({r}_{{{{\rm{c}}}}}=10\sigma\) is as low as about \({10}^{-6}\). Packing modes of polygons with various \({r}_{{{{\rm{c}}}}}\) values exhibit no qualitative differences, as shown in Supplementary Fig. S16.

All the simulations started from the crystalline state unless otherwise specified and were performed at each temperature independently. To construct a crystalline structure, the unit-cell structure was first determined and then duplicated along two primary lattice vectors. In our polygonal systems, we determined the unit-cell structure by careful examination on the final structure from the REMD simulations described below. For triangle, the unit-cell contains two monomers with opposite body-orientations, as shown in Fig. 2c; for the other four polygons, the unit-cell is just a monomer. The exact values of the two primary lattice vectors were calculated by geometrical construction, taking the effective radius of each L-J particle into account. For the simulations of the triangular system starting from a crystalline structure with a locally disordered region, the initial configuration is prepared by the following procedure: At a relatively high temperature of T = 2.88, the triangular system becomes amorphous spontaneously. When this happens, we select a snapshot with a small disordered region (about 30 monomers). This snapshot is then used as the initial configuration for related simulations.

MD simulations

All the simulations in this work were performed with LAMMPS⁶³. After the hand-made initial crystalline structures were optimized, they were tailored to fit in a simulation box with appropriate side-length ratio and tilting angle under periodic boundary condition. All the results reported in this paper were obtained from the simulations with 2496 polygonal molecules. Additional simulations on larger scales, 22,500 monomers for the square system and about 5000 monomers for other polygons, found no qualitative differences in phase behaviors.

In NPT simulations, the shape of the simulation box controlled by the tilting angle and box sides ratio was invariant since the algorithm allows only isotropic shrinkage or stretching of the box⁶³. Although the simulation box is perfectly commensurate with the ground-state crystal, there is no guarantee that it is still appropriate for crystalline morphology at relatively high temperatures. Therefore, other than NPT simulations for the triangular system, simulations in the isotension-isothermal ensemble⁴³ were performed for the other four systems. In each simulation, a very short NVT relaxation had been carried out for \(6\times {10}^{3}\) steps before an NPT or isotension-isothermal equilibration run for \(8.3\times {10}^{6}\) steps was performed. Since the typical correlation time was determined to be less than \(1.6\times {10}^{3}\) steps, \(2\times {10}^{3}\) configurations were then sampled with the interval of \(2\times {10}^{3}\) steps. Because a longer simulation time might be required in each equilibration run, \(4\times {10}^{6}\) more steps were carried out before sampling near the melting point. In very rare cases, the crystal melts at T but stabilizes at a higher temperature \(T^{\prime}\). For these cases, we performed simulations at \(T^{\prime}\) initiating from the final configuration at T or reseeded the random numbers used in the simulations. In this way, the system always melted/stabilized at the liquid state, which ensured that the previous simulation at \(T^{\prime}\) was stuck in the superheated state. This was also confirmed by the fact that the simulations at other temperatures slightly higher than T resulted in melting as well. The Nosé-Hoover barostat and thermostat^64,65,66,67 were used to keep the system pressure and temperature constants. The simulation timestep was 0.005, much smaller than the characteristic time \(\tau=1/\sqrt{k}\), where k is the strength of the valence bond, which is the fastest degree of freedom in the system.

For each polygon, the lowest pressure we study is P = 0, and the highest is P = 10 (P = 8 for octagon). An even higher pressure may cause numerical instability.

REMD simulations

First, 625 monomers were put on the square lattice (25\(\times\)25) with the same body-orientation, different from the equilibrium structure of any polygons to avoid the influence from the initial configuration. REMD simulations were then carried out in the NVT ensemble with 26–30 replicas. We have used uneven temperature intervals for different polygons to ensure that the exchange probabilities are all around 20–35%, and the temperature distribution is shown in Supplementary Fig. S1. The simulation temperature was kept constant by the Langevin dynamics. All the simulations lasted \(7.2\times {10}^{7}\) steps and exchange attempts took place every \(3\times {10}^{5}\) steps. For each polygon, after the simulation was finished, an energy minimization was performed on the final configuration at the lowest temperature to quench thermal fluctuations, and the optimized structures are shown in Supplementary Fig. S2. Because the purpose of performing the REMD simulation was to investigate local structural features, the REMD simulation employed a reflect-wall boundary condition to avoid the influence of the box shape or the periodic boundary condition, which may prevent the formation of a dense but amorphous state like the spin-ice-like state that cannot fit into a regular box with periodic boundary condition seamlessly.

Determinations of phases and phase transitions

The bond-orientational order parameter is defined as \({\varPsi }_{n}=| < \frac{1}{N}{\sum}_{i}{\psi }_{n,i} > |\) and \({\psi }_{n,i}=\frac{1}{n}{\sum}_{j}{e}^{in{\theta }_{ij}}\), where N is the total number of unit-cell, i runs over all the representative points of the unit-cells while j runs over the n nearest neighbors of point i, < > means taking the ensemble average. For triangle, the representative points of the unit-cells are chosen as the COMs of molecules with the same initial body-orientation, and since we study the non-equilibrium dynamical process, \({\varPsi }_{4}\) is calculated for each snapshot instead of taking the ensemble average; for the other four polygons, they are chosen as the COMs of all molecules. The characterizations of the spin-ice-like state of the triangular system and the distorted lattice of the square system have been described in the main text, and it is worth emphasizing that the Delaunay tessellation/Voronoi diagram in Fig. 5c is calculated by the freud library⁶⁸ with periodic boundary condition. Below we focus on the identifications of the crystal, hexatic, and liquid states. The liquid state has an approximately zero bond-orientational order parameter \({\varPsi }_{n}\), the hexatic state has an exponentially decaying translational correlation function as well as a non-zero bond-orientational order parameter, and the crystal state has an algebraically decaying translational correlation function as well as a non-zero bond-orientational order parameter. The bond-orientational order parameter was calculated by the analysis tool in LAMMPS⁶³.

The calculation of the translational correlation function is a key point for the phase determination of 2D melting⁶⁹. In 2011, the pioneer work by Bernard and Krauth provided a method to evaluate the correlation function²¹, which was later modified by Li and Ciamarra for the L-J systems⁷⁰. The correlation function is defined as \({C}_{g}(r,\theta )= < \delta (r-({x}_{i}-{x}_{j}))\delta (r\,\tan \theta -({y}_{i}-{y}_{j})) > \), where θ is the orientation of the crystal axis and should be determined by scanning \({C}_{g}(r,\theta )\) at a fixed r (the first peak value corresponding to the orientation of the crystal axis).

Basically, a first-order phase transition can be identified by the discontinuous change in the caloric and density curves. We also confirm the discontinuous transitions by the corresponding order-parameters: body-orientational order parameter for solid-solid transition and bond-orientational order parameter for solid-liquid or hexatic-liquid transition.

The continuous transition for 2D crystal is theoretically a Kosterlitz–Thouless one instead of a strict second-order phase transition that can be evidenced by the divergence of heat capacity, so it is difficult to be identified by the potential energy or its derivative. Since we do not observe the KTHNY melting scenario in all our simulations, the only continuous phase transition we observe is the solid-hexatic transition. Therefore, we determine the phase transition point by the decaying behavior of the translational correlation function: At the phase transition point, it suddenly changes from algebraical decay to exponential decay.

All the results for caloric curves as well as the calculations of order parameters and correlation functions can be found in SI (texts from Supplementary S4–S6 and figures from Supplementary Figs. S10–S15).