Solid-to-molecular-orientational-hexatic melting induced by local environment determined defect proliferations

Zhanglin Hou(侯章林) Jieli Wang(王杰利) Ying Zeng(曾颖) Zhiyuan Zhao(赵志远)Xing Huang(黄兴) Kun Zhao(赵坤) and Fangfu Ye(叶方富)

1Wenzhou Institute,University of Chinese Academy of Sciences,Wenzhou 325001,China

2Key Laboratory of Systems Bioengineering(Ministry of Education),School of Chemical Engineering and Technology,Tianjin University,Tianjin 300072,China

3School of Physical Science,University of Chinese Academy of Sciences,Beijing 100049,China

4Research Center of Computational Physics,School of Mathematics and Physics,Mianyang Teachers’College,Mianyang 621000,China

5Sichuan Provincial People’s Hospital,University of Electronic Science and Technology of China,Chengdu 610054,China

6Institute of Fundamental and Frontier Sciences,University of Electronic Science and Technology of China,Chengdu 610054,China

7Beijing National Laboratory for Condensed Matter Physics,Institute of Physics,Chinese Academy of Sciences,Beijing 100190,China

8Songshan Lake Materials Laboratory,Dongguan 523808,China

9Oujiang Laboratory(Zhejiang Laboratory for Regenerative Medicine,Vision and Brain Health),Wenzhou 325001,China

Keywords: two-dimensional melting, phase transition, molecular-orientational forming liquid crystal, local environment,melting mechanism

1. Introduction

Understanding two-dimensional (2D) melting is a longterm subject in condensed matter and statistical physics.[1,2]In late 1950s,Bernal introduced geometrical degrees of freedom to understand the entropy differences between solid and liquid phases, which has opened a door for an intuitive understanding of 2D melting.[3,4]Afterwards, motivated by the absence of true long-range order in two dimensions,[5]Berezinskii,[6,7]as well as Kosterlitz, Thouless, Halperin, Nelson and Young,developed a novel theory (BKTHNY theory) for 2D melting based on the evolution of topological defects in systems,[2,8–13]which states a two-step melting process:the unbinding of pairs of dislocations drives a continuous transition from solid to liquid crystal phase (x-atic), and subsequently, the dissociation of pairs of disclinations(i.e.single dislocations)drives the second continuous transition from liquid crystal phase to isotropic fluid. The BKTHNY theory strongly inspires researchers to investigate 2D melting behavior.[14–28]On the other hand,with the advances of detection and analysis techniques,other transition behaviors that the BKTHNY mechanism cannot cover have also been discovered,[19,22,24–26,28–34]suggesting that our understanding of 2D melting is still not complete.

The melting mechanism suggested by the BKTHNY theory have been tested in different particle systems,[2,8–10,24,28]and three types of 2D melting scenarios have been revealed up to date,[17,19,22,26,28,35–37]which are a first-order transition directly from solid to isotropic fluid, a two-step melting with a continuous solid-to-liquid crystal transition and a firstorder liquid crystal-isotropic fluid transition, and a two-step melting with continuous transitions for both solid-liquid crystal and liquid crystal-isotropic fluid. For example, both hard disks[19,22,35]and regular polygons with an edge numbern ≥7[28]show a two-step melting behavior including a continuous solid-hexatic transition and a weakly first-order hexatic-fluid transition, which is different from the BKTHNY theory. By contrast,melting of hard hexagons fully obeys the BKTHNY theory as investigated by Andersonet al.[28]These results imply that the melting transition is closely related to the symmetries of particle shape and their assembled structures.[28,31]However, the underlying mechanism is still unclear. Particularly,how the melting mechanism is affected by the symmetry of initial crystal configurations remains largely unknown.

In this work, by performing Monte Carlo (MC) simulations, we systematically studied the melting behavior of 60◦hard rhombs[see Fig.1(a)]in a square box. This system has two types of crystal structures including a complex hexagonal crystal [cHX, Fig. 1(b)] whose primitive cell contains three rhombs, and a simple rhombic crystal [RB, Fig. 1(c)] whose primitive cell contains just one rhomb. Both cHX and RB structures are able to fully cover the surface (i.e., the packing density is 1). In addition, because rhombs have internal angles of 60◦and 120◦,if only centers of particles are considered, the RB has the same symmetry as a hexagonal crystal,

whereas cHX has a symmetry of honey-comb lattice (i.e., a symmetry of inverse-hexagonal lattice). Thus,this 60◦rhomb system provides an ideal platform to investigate the dependence of the melting mechanism on initial crystal structures.Our results show that the local environment of particles is a key factor in 2D phase transitions and the melting behavior is initial-configuration-dependent. Our findings would benefit a better understanding of 2D melting.

This paper is organized as the follows: In Section 2, we introduce the simulation methods used in this work. In Section 3,the results are shown and discussed. Lastly,the conclusion is given in Section 4.

2. Simulation methods

3. Results and discussion

3.1. Simulation results: initial-configuration-dependent melting behavior

Fig.1.Schematic illustration of simulation models and typical configurations of each phase obtained from Monte Carlo(MC)NPT ensembles. (a)Sketch of a 60◦rhomb with an edge length L.(b)Illustration of a complex hexagonal crystal(cHX,green rhomb indicating a primitive cell of the lattice)that can be viewed as an assembled structure of a three-particle-based hexagonal unit(red hexagon). (c)Illustration of a simple rhombic crystal(RB,green rhomb indicating a primitive cell of the lattice) constructed by single particle unit(red rhomb). (d)–(g)Snapshots of I,Hmo, RB and cHX phases. I:isotropic fluid, Hmo: hexatic in molecular-orientational order. Insets at the lower left corners: FFTs of corresponding snapshots.

For both initial structures of cHX and RB crystals, systems melt first from a crystal to an intermediate liquid crystal phase Hmoand subsequently from the Hmoto an I phase.Transitions for both steps seem to be continuous.Thus,the melting transitions in both cHX and RB systems are consistent with the BKTHNY theory.[2,8,9,11,28]Although the phase sequence is similar during melting transition in the cHX and RB systems,they do show a different solid–Hmotransition point, which isφA=0.812 for the cHX system andφA=0.828 for the RB system. Apparently, such a difference in the transition point results from the different initial crystal structures.

Fig. 2. Monte Carlo simulation results of the cHX system using NPT ensembles. (a) Equation of state and (b) order parameters including 2-, 4-, and 6-fold molecular orientational orders Φ2, Φ4, and Φ6, 4- and 6-fold bond orientational orders Ψ4 and Ψ6, and translational order with cHX lattice ScHX. The vertical dashed lines indicate the cHX–Hmo and Hmo–I transition points. [(c),(d)]Spatial correlation functions of 6-fold molecular orientational orderand translational order gcHX. Green dashed lines in(c)and(d)decay as r−1/4 and r−1/3, which are the BKTHNY predictions on liquid crystal-isotropic fluid and solid-liquid crystal transition points for gmo6 and gs,respectively.

Fig. 3. Monte Carlo simulation results of the RB system using NPT ensembles. (a) Equation of state and (b) order parameters including 2-, 4-, and 6-fold molecular orientational orders Φ2, Φ4, and Φ6, 4- and 6-fold bond orientational orders Ψ4 and Ψ6, and translational order with RB lattice SRB. The vertical dashed lines indicate the RB–Hmo and Hmo–I transition points. [(c), (d)] Spatial correlation functions of 6-fold molecular orientational orderand translational order gRB. Green dashed lines in(c)and(d)decay as r−1/4 and r−1/3,which are the BKTHNY predictions on liquid crystal-isotropic fluid and solid-liquid crystal transition points for nd gs,respectively.

3.2. Definition of LPC-4NNs

To further understand the mechanism underlying the observed different melting behaviors, we employed an analysis based on local polymorphic configurations (LPCs). The LPC method has been used to analyze the glass transition in Penrose-kite systems[39]and to understand the emergent tetratic order in hard kite systems,[38]suggesting that it could be used to analyze the phase behavior of anisotropic particles.By applying the LPC method to the rhomb system, rhombic particles are first classified based on their local environment formed by their four nearest neighbors,i.e.,LPC-4NNs. Here,four nearest neighbors are chosen in view of the structural characteristics of cHX and RB lattices as well as the mechanical stability of edge-edge contact between neighboring particles. For each selected particle, its neighboring particle with a relative pointing angle between them being around 60◦/120◦or 0◦/180◦is categorized as a 60◦/120◦or 0◦/180◦neighbor,respectively; while a neighbor with a relative pointing angle being out of 60◦/120◦and 0◦/180◦is categorized as a fluctuating neighbor.Then,for the selected particle,based on the category information of its 4 nearest neighbors,it can be classified into one of six types with each one having a different LPC-4NN [Fig.4(a)]: Types 1–5 have zero to four 0◦/180◦neighbors,respectively. Type 6 has 2 or more fluctuating neighbors,and particles of type 6 are called disordered particles (DPs).For simplicity,for particles that have one fluctuating neighbor,the threshold value for judging 60◦/120◦or 0◦/180◦neighbors is increased from 15◦to 30◦,and the fluctuating neighbor under the threshold value of 15◦was re-categorized under the threshold value of 30◦, so that these particles can be classified into one of types 1–5 (see details in the supplementary information and Fig. S3). Figure 4(b) shows a representative configuration of the cHX system in which particles are colored according to their LPC-4NNs.

Fig. 4. The model of LPC-4NNs and a representative configuration with particles colored according to LPC-4NNs. (a) Schematic graphs showing six types of LPC-4NNs.Brown and dark violet dashed rectangle illustrate type 1(consistent with cHX structure)and type 5(consistent with RB structure), respectively. Green dashed rhombs in LPC-4NNs: indicating that its neighboring rhombs can slide along the edge of the target particle to the position of the green dashed rhomb without changing the type of LPC-4NN. (b) A representative configuration at φA=0.800 in the cHX system,in which particles are colored by their LPC-4NN types. Colors of the six LPC-4NN types are same as in(a).

3.3. Partial translational and molecular-orientational orders: effect of each type of LPC-4NN on melting

LPC-4NNs provide an intuitive view of local environment around target particles, which can be closely related to the global order of the system.[38,39]Based on LPC-4NNs,we measured the partial translational and molecular orientational orders for the particles that have the same type of LPC-4NNs.The results are shown in Fig. 5. The global orders obtained from all particles of the system (i.e., no matter what is type of particle’s LPC-4NN)are also shown,which can be used as references to judge the contribution of the particles with certain type of LPC-4NNs. We can see that in the cHX system[Fig. 5(a)], before melting, the globalScHXbegins to drop as soon as non-type-1 particles appear. Moreover, after type-2 particles appear, the partialScHXof them is roughly equal to the globalScHX. Then in the Hmoregion,as the global translational order has been destroyed,the measured globalScHXdecays quickly to below 0.4.Compared with the globalScHX,the partialScHX’s of type-1 and type-4 particles are higher whereas the partialScHX’s of type-3 and type-6 particles are lower,indicating that type-3 and type-6 particles are more detrimental than type-1 and type-4 particles in destroying the translational order of the system.Interestingly,for the whole crystal region,the partialScHXof type-1 particles seems to keep constant and exhibits a plateau. These results suggest that, when the system melts from cHX, the generation of other LPC-4NNs different from type 1 will reduce the globalScHX, though it has little effect on the remaining structures formed by type-1 particles. By contrast, for the RB system [Fig. 5(c)], before melting the globalSRBis roughly equal to the partialSRBof type-5 particles even when non-type-5 particles appear, indicating the dominant role of type-5 particles in the system.Near the RB–Hmotransition,type-4 and type-5 particles contribute positively while those of types 1–3 contribute negatively to the global order. In the crystal region,the overall trend of the partialSRBfrom type-5 particles shows a gradually decrease with decreasingφA. In both systems, near and above solid–Hmotransitions,the number of DPs is too few(less than 1%)to get statistically meaningful results.

Similar analysis for the 6-fold molecular orientational orderΦ6has also been performed[Figs.5(b)and 5(d)],and the results show that near Hmo–I transition,Φ6is mainly destroyed by DPs. This is same for the cHX and RB systems, as is expected.

Fig. 5. Translational order [(a), (c)] and 6-fold molecular orientational order [(b), (d)] measured for each type of particle at different φA.Cyan-filled region indicates the global order of the system. Results shown in the figure are averaged over more than 100 particles for each type of particle. The dashed lines indicate the transition points.

3.4. Proliferation of point/line defects induces the cHX–Hmo/RB–Hmo transitions

To further elucidate the roles of different LPC-4NNs in the melting transitions, we tracked the evolution of different LPC-4NNs during the melting of the two systems. The results are shown in Fig. 6. For the cHX–Hmotransition, at the beginning(e.g., atφA=0.844), particles mostly belong to type 1 [Fig. 6(b)]. Type-2 and type-3 particles can be observed sporadically and often appear in a form of point-like groups[e.g.,isolated star-like clusters in Figs.6(b)and 6(c)]consisting of three type-3 and six type-2 particles. AsφAdecreases toφA=0.825 [Fig. 6(c)], the amount of type-2 and type-3 particles increases, so more point-like groups are formed and some neighboring groups connect to each other to form bigger clusters. Then, atφA=0.805, the system melts into Hmophase. The numbers of both type-2 and type-3 particles reach a peak [Fig. 6(a)], and these particles connect to form a big cluster that can percolate the system [Fig. 6(d)]. Other types of particles also appear inside the big cluster. In contrast,the evolution of LPC-4NNs shows a different pattern during the RB–Hmotransition. AtφA=0.837[Fig.6(g)]where the system is still in RB phase, most particles belong to type 5, but particles of types 1–4 can emerge and often appear simultaneously. However, different from the point-like groups observed in the cHX system, these particles of types 1–4 form line-shaped clusters in the RB system. At a lowerφA=0.830[Fig.6(h)],with more non-type-5 particles appear,the number of line-shaped clusters increases and they begin to crosslink to form a big network, which can percolate the system even when the system is still in the RB phase. Then atφA=0.823[Fig.6(i)],the expanded network essentially covers the majority of particles and the system melts into Hmophase. For both systems,as the system continues to melt,more and more DPs appear [Figs. 6(e) and 6(j)], which destroy theΦ6order and drive the system from Hmoto I phase.

The observed different evolution patterns of LPC-4NNs seem to suggest that the generation mechanisms of the cHX melting and RB melting for the LPC-4NNs that are not consistent with the crystal structures (hence defects associated) are different. To better illustrate this point,we consider defects in the two systems. As shown in Fig.7(a),in the cHX structure,if a three-particle-based hexagonal unit(red dashed hexagon)rotates±π/3, it will generate six type-2 particles and three type 3-1 particles. Then, when two such rotated hexagonal units are neighbors,type 3-2 and type 3-3 particles are created.Similarly, type-4 and type-5 particles can be generated when three and four such rotated hexagonal units come to close with each other. Since all other types of LPC-4NN except type 1 are not compatible with the crystal structure, these generated particles of types 2–5 thus contribute to the defects of the system and drive the melting of cHX. Based on this picture, in the context of occurrence probability,the generation of different LPC-4NNs is sequential in the cHX system where type 2 and type 3-1 appear first,next followed in turn by type 3-2 and type 3-3,type 4 and type 5. This sequential generation is also consistent with the evolution curves of different LPC-4NNs shown in Fig.6(a).

Fig.6. Evolution of melting process characterized by types of particles according to LPC-4NNs. [(a),(f)]The fraction of of particles of types 1–6 as a function of packing fraction in the cHX system(a)and the RB system(f). The colored backgrounds represent different phase windows shown in Figs.2 and 3. Representative configurations of the cHX system(b)–(e)and the RB system(g)–(j)with particles colored by its type based on LPC-4NNs.

In contrast, in the RB structure, most LPC-4NNs that are not compatible with the RB structure are associated with twinning-like line defects. For example, Fig. 7(b) top shows some typical line defects of 60◦rhombs, in which particles of type 1, type 2, type 3-1, type 3-3 and type-4 are found.When more than 2 line defects meet, type 3-2 particles can be generated [Fig. 7(b) bottom]. Thus, during the melting of the RB structure,all types of LPC-4NNs that are not compatible with the crystal structure can appear simultaneously[consistent with the evolution curves shown in Fig. 6(f)], which can work synergistically in disrupting the crystal structures and drive the RB–Hmotransition. This is different from the sequential appearance of defect-associated LPC-4NNs in the cHX case. Such differences also result in different solid–Hmotransition points shown in Figs.2 and 3,where RB–Hmotransition happens earlier(i.e.,at higherφA)than cHX–Hmotransition. This is because the RB–Hmotransition is more sensitive to any type of defect-associated LPC-4NNs as these LPC-4NNs occur simultaneously, whereas for the cHX–Hmotransition, the sequential appearance of defect-associated LPC-4NNs leads to a gradual reduction of crystal orders and thus pushes the transition point to a lowerφA. We would like to point out that all these aforementioned analyses are based on defects that are loosely indicated by those LPC-4NNs that are not consistent with crystal structures. They are not equivalent to typically used dislocations and disclinations that are rigorously defined. Thus,to quantitatively characterize the defects identified in this work(e.g.,core energies of defects)and interactions between defects, new methods need to be developed,which are beyond the scope of this study.

Fig.7. Illustration of the generation mechanisms of LPC-4NNs in the cHX system (a) and RB system (b). (a) In the cHX system, rotation of a three-particle-based hexagonal unit (red dashed hexagon) excites LPC-4NNs of types 2 and 3-1(green rectangle),and several neighboring rotated units can excite LPC-4NNs of type 3-2, type 3-3, type 4 and type 5. (b)Illustration of twinning-like line defects including LPC-4NNs of type 1, type 2, type 3-1, type 3-3 and type 4 (top) in the RB system. When more than two line defects meet, type 3-2 LPC-4NN(bottom)can be generated.

4. Conclusion

In summary, by MC simulations, we have studied the melting behavior of 60◦rhombs from two different initial crystal structures, cHX and RB. Using the LPC-based analysis method, we found that the mechanisms of the meltings of two systems are different although both of them obey the BKTHNY theory. In the cHX melting, the proliferation of point defects formed by LPC-4NNs of types 2–6 that are not compatible with the cHX crystal structure drives the cHX–Hmotransition; whereas in the RB melting, line defects consisting of LPC-4NNs of types 1–4 and type 6 that are not compatible with the RB crystal structure are responsible for destroying the crystal order and lead to the RB–Hmotransition.Moreover,defect-associated LPC-4NNs are generated sequentially in the cHX system while they can be generated simultaneously in the RB system. These differences also correlate to the observed different solid–Hmotransition point which isφA=0.812 for the cHX–Hmotransition andφA=0.828 for the RB–Hmotransition. We note that here we mainly focus on the defects mediated solid–Hmotransitions. In principle,similar analysis can also be performed for the Hmo–I transition.However,such analysis will require a new way to detect Hmoassociated disclinations(e.g.,±1/6 disclinations). Moreover,there is no difference between the Hmo–I transition of the cHX and RB systems. Given such considerations,it is more appropriate to leave such defect analysis for the Hmo–I transition in future work.

The different melting manners in the cHX and RB systems reveal the important role of local environment formed by neighboring particles,which leads to the initial-configurationdependent behavior in our work. Thus,the phase behavior of 2D anisotropic particles is not only dependent on the details of particle-self,e.g.,shape and symmetry,but also on the initial configuration preset in the system. The results obtained in this study brings new challenges and opportunities for understanding 2D melting on the theoretical level and demonstrates the potential of the LPC-based analysis method as a powerful technique in understanding the phase behavior of anisotropic particles.Based on LPC,it may also help to control the assembled structures in a condensation process(which is the reverse of the expansion process used in MC simulations),by manipulating the evolution pathways of LPCs. With the advances in the fabrication of anisotropic colloids,for example,through lithographic techniques,[40]we would expect that controlled assembling through such LPC-based regulations can be tested experimentally in the future.

Acknowledgements

We thank the Post-Doctor Association of WIUCAS.This work was carried out at the National Supercomputer Center in Tianjin and the calculations were performed on TianHe-1(A).

Project supported by the National Natural Science Foundation of China(Grant Nos.11874277,21621004,12104453,and 12090054)and the Strategic Priority Research Program of Chinese Academy of Sciences(Grant No.XDB33030300).