Emergent elastic fields induced by topological phase transitions: Impact of molecular chirality and steric anisotropy.

SignificanceChirality, the property of an object that cannot be superimposed on its mirror image, plays an essential role in condensed matter, such as magnetic, electronic, and liquid crystal systems. Topological phases emerge in such chiral materials, wherein helical and vortex-like structures-called skyrmions-are observed. However, the role of elastic fields in these topological phases remains unexplored. Here, we construct a molecular model of two-dimensional crystals incorporating steric anisotropy and chiral interactions to elucidate this problem. The coupling between the elastic fields and phase transitions between uniform, helical, and half-skyrmion phases can be utilized to switch these topological phases by external forces. Our results provide a fundamental physical principle for designing topological materials using chiral molecular and colloidal crystals.

Understanding and controlling phase transitions into spatially modulated structures are of importance in technological applications, such as spatial light modulation by cholesteric liquid crystals with helical structures (1). Their microscopic controllability can be utilized to control materials properties, including mechanical rigidity (2). From the viewpoint of fundamental physics, on the other hand, the nature of formation processes of spatially modulated phases remains elusive compared with phase transitions between uniform states, though there have been extensive studies on the role of critical fluctuations (3–5). Skyrmions and half-skyrmions (merons) are representative objects of the modulated structures. They are of interest in magnetic systems (6–8), Bose–Einstein condensates (9, 10), quantum Hall systems (11), dielectrics (12, 13), liquid crystals (14, 15), and active matter (16, 17). These objects are manipulated in designing magnetoelectric, optical, and mechanoresponsive materials (18–20). In liquid crystals, in particular, molecular chirality plays an essential role in the formation of helical and half-skyrmion phases (1, 19, 21, 22).

In theoretical description of chiral liquid crystals, a spontaneous twisting term reflecting chirality is incorporated to examine the formation of mesoscale cooperative structures (1, 14, 21, 23). This term resembles Dzyaloshinskii–Moriya interaction in chiral magnets (6, 7), incorporating liquid crystalline symmetry of the order parameter. Thus, the appearance of helical and half-skyrmion phases in cholesteric liquid crystals is explained (14). Recently, furthermore, the half-skyrmion phase in metal–organic frameworks (MOFs) was theoretically proposed using a similar lattice model (24). Therefore, it is reasonable to expect the existence of similar phase transitions in other chiral molecular systems, such as organic crystals (25, 26), colloidal crystals (27, 28), and biological systems (29), by utilizing the same theoretical framework, though the existence of half-skyrmion structures in these substances has not yet been investigated.

However, there is a crucial difference between liquid crystals and molecular crystals. In molecular crystals and MOFs, the change of molecular configurations associated with phase transitions induces lattice distortion, resulting in the emergence of the elastic field, though the role of the elastic field remains elusive. We call this elastic field an “emergent elastic field,” which is absent in cholesteric liquid crystals due to its liquid nature in molecular translation. The emergent elastic field is an analog of the emergent electromagnetic field in condensed matter, which is known to play an essential role in manipulating skyrmions (30, 31). Therefore, it is of crucial importance to understand how the emergent elastic field is linked to topological phase transitions in molecular crystals.

Here, we reveal an essential role of the emergent elastic field associated with topological phase transitions into helical and half-skyrmion phases in a crystal. We develop a molecular dynamics model incorporating intermolecular chiral and spheroidal steric interactions. These interactions model pair interaction between anisotropic colloidal molecules with chirality. In spheroidal molecules, molecular rotation induces lattice distortion. By controlling the steric anisotropy and chiral twisting interaction between adjacent molecules, phase transitions into helical and half-skyrmion phases in a crystal are controlled. The long-range nature of elastic correlation is essential to control phase behavior and domain formations. As a consequence of this elastic coupling, these phases can be switched by varying temperatures and by applying external electromagnetic fields and anisotropic stresses. The comparison between magnetic solids, liquid crystals, and our systems is summarized in Table 1. Our work reveals a link between topological phase transitions and the elastic fields, which provides a control system for designing molecular and colloidal crystals with tunable electro- and magneto-mechanical properties.

Table 1.

Characteristic features of our system

	Chiral magnets	Chiral liquid crystals (cholesterics)	Monolayer chiral crystals (our system)
Interaction	Dyzaloshinskii–Moriya (6, 7)	Liquid crystal elasticity (1)	Intermolecular twisting, steric repulsion
Characteristic phases	Skyrmion crystal, etc. (7).	Blue phases, half-skyrmions, etc. (14, 15, 19, 21)*	Half-skyrmions
Size of (half-)skyrmions	5 to 100 nm (7, 18)	100 nm to μm (15)	Several times particle diameter†
Emergent field	EM fields (7, 30, 31)	Nematic director fields (14, 15, 19)	Elastic fields‡
Phase controllability	EM fields, heat, strain (7, 32)	EM fields, heat, bending, fluid flow (1)	EM fields, heat, anisotropic stress
Melting of (half-)skyrmions	KTHNY (33, 34)	Not available	KTHNY

Here, EM field denotes the electromagnetic field.

*Full-skyrmions and quarter-skyrmions have also been observed (19).

†∼ 10 nm for molecular crystals and m for colloidal crystals.

‡Produced by topological defects and phase ordering.

Results

Phase Diagram with Respect to Twist Parameters.

We construct a molecular dynamics model by incorporating twisting interactions between adjacent molecules, which is schematically described in Fig. 1. The alignment of neighboring molecules is stable when they exhibit a certain twist angle determined by q0, and the rigidity of this twist is given by K2 (see Materials and Methods for the definitions of these parameters). These molecules are confined to realize a monolayer geometry, as shown in Fig. 1, forming two-dimensional hexagonal crystals (see for larger thickness geometries). There are two-step phase transitions between liquids, orientationally disordered crystals, and orientationally ordered crystals shown in Fig. 1. The melting temperature is much higher than that of the orientational phase transition. In the present study, we focus on low-temperature range, wherein the two-dimensional hexagonal crystal is stable (see Materials and Methods for details). Using this model system, we demonstrate phase controllability via the manipulation of q0 and K2 at low temperatures, as displayed in Fig. 1. For small K2 and q0 values, almost all the molecules are oriented normal to the layer to form a uniform phase, as displayed in Fig. 1. In this state, the structure factor (see Materials and Methods for this definition) exhibits only the Bragg peaks representing the nearest-neighbor particles, and no mesoscopic structure forms. As the twist rigidity K2 is increased, the number of molecules that orient tangential to the monolayer increases. For large q0 values, long-range orders of helical and half-skyrmion states are formed, as displayed in Fig. 1 , respectively. Conversely, for small q0 values (), particle orientation does not vary continuously in space, and orientational defects are dispersed heterogeneously, exhibiting a halo structure at low wavenumbers, as shown in Fig. 1. As K2 is increased, the number of defects increases and eventually results in the division of ordered particles into compartments (Fig. 1). The transitions between these structures are gradual, except for the sharp transition between the helical and the half-skyrmion phases (). Therefore, we successfully control the topological structure of this system by controlling the twisting interactions. Here, it should be noted that the mother phase (high-temperature orientationally disordered crystal phase) always forms a hexagonal lattice. The relationship between underlying crystal structures and topological structures is beyond the scope of this work.

Fig. 1.

Phase behavior in a monolayer geometry. (A) Molecular twisting in our model. The favored twist angle is given by q0, and the rigidity of the twist is determined by K2 (see Materials and Methods for this definition). (B) Monolayer geometry of this study (see Materials and Methods for details). (C) Phase diagram with respect to q0 and K2 at low temperatures. The phase boundaries are shown for reference. The star symbol denotes the phase point examined in Fig. 3. (D) Color notation for the snapshots in E–I and the following figures. (E–I) Diagonal view of the real spatial structure and corresponding structure factor (see Materials and Methods for this definition) at each phase. The temperature is T = 0.05 in this figure.

Fig. 3.

Thermal stability of the half-skyrmion phase. (A) Hexatic correlation function divided by the radial distribution function (see Materials and Methods for this definition) of the half-skyrmion cores. The hexatic correlation increases as the temperature decreases via the characteristic power-law decay. (B) Pair correlation function of the half-skyrmion cores. It exhibits a power-law decay in the solid phase, whereas it decays exponentially in the hexatic and liquid phases. (C–E) Real spatial structures at T = 0.42 (C; half-skyrmion crystal), T = 0.44 (D; half-skyrmion hexatic), and T = 0.46 (E; half-skyrmion liquid). We only display particles with . (F–H) Corresponding Delaunay triangulation. The red (yellow) circles denote seven (eight)-member coordinated vortices, whereas the blue (green) crosses denote five (four)-member coordinated vortices. and in this figure, which are indicated by the star symbol in Fig. 1. (Scale bars [C–H]: .).

Emergence of Elastic Fields in Topological Phases.

The characteristic feature of our system is the emergence of strain and stress fields. In this study, a two-dimensional crystalline state is realized at a high temperature because the average density of this system is high, and orientational phase ordering proceeds without particle migration. Because each molecule exhibits steric anisotropy, the tangential orientation generates strain and stress around the molecule, which is similar to the generation of an elastic field with a long-range spatial correlation around an inclusion, a precipitate, and a defect (35–38). As mentioned in the introduction, we call this elastic field an emergent elastic field. The concept of the emergent elastic field can explain the domain formation shown in Fig. 1, wherein the helical state is divided into domains with different wave vectors, whose angles are approximately . From the real space structure depicted in Fig. 2, we calculate the emergent strain fields and stress fields in Fig. 2 , respectively (see Materials and Methods for the definition). As shown in Fig. 2 , a large anisotropic strain is induced inhomogeneously because molecular orientations in different helical domains produce different strain components. If a uniform helical structure is formed, it produces a large uniform strain that generates large anisotropic pressure and thus becomes mechanically unstable. By forming domain structures, the anisotropic stress is macroscopically reduced and localized, as shown in Fig. 2 . Such domain formation does not occur for the half-skyrmion state shown in Fig. 1 because half-skyrmions form an isotropic hexagonal solid at low temperatures (see for the elastic fields of uniform and half-skyrmion phases, respectively). The emergent elastic field is not relevant in cholesteric liquid crystals due to its liquid nature in molecular translation, but it can be crucial in chiral smectics (1) and elastomers (39) because molecular rotations induce layer compression and network deformation, respectively. The investigation of (half-)skyrmion structures in these substances, molecular crystals, and colloidal crystals should be conducted to reveal the impact of the emergent elastic fields.

Fig. 2.

Emergent elastic field of a helical state. (A) Particle configuration corresponding to Fig. 1, wherein the system is divided into two major domains with different particle orientations (indicated in blue and green). (B–D) Strain field of this configuration (see Materials and Methods for this definition). The volumetric strain represents the volume dilation and compression and is small. The emergence of a large macroscopic uniaxial strain in C and a large shear strain in D results from the formation of a helical state, wherein tangentially aligned particles produce large strains along their orientations. (E–G) Stress field of this configuration (see Materials and Methods for this definition). (E) Isotropic stress corresponding to the volumetric strain via the bulk modulus (volume compressibility). (F and G) Uniaxial stress (F) and shear stress (G), which exhibit the same spatial pattern as the particle configuration, unlike the strain correspondence. (Scale bars: .).

In this study, we consider monolayer systems without deformation toward the third dimension (z-direction). To examine the stability of the planar monolayer without bending, we calculate the displacement and the force field along the z-direction in . Both the displacement and force fields do not exhibit long-wavelength fluctuation, while their heterogeneous distribution is confirmed in short wavelength. This result implies that a macroscopic deformation of the substances would not occur.

Thermal Phase Transitions of the Helical and the Half-Skyrmion Phases.

The phase diagram shown in Fig. 1 indicates the stable structure at a low temperature. Previous studies often examined the phase transitions between topological phases by varying the temperature and external perturbations, such as the magnetic field. In this section, we investigate the thermal stability of the helical and the half-skyrmion phases. For the former phase, we demonstrate the thermal hysteresis of the helical phase in . As the temperature increases, the helical state transforms into the half-skyrmion state, which is indicated by the abrupt change in the system’s potential energy. This phase transition is reversible and exhibits a hysteresis loop, implying its first-order nature. For the latter phase, the vortex structure in the half-skyrmion phase is stable over a wide temperature range. However, another phase transition with respect to the long-range ordering of half-skyrmions emerges. The transitions between the liquid, hexatic, and crystalline half-skyrmions occur, as shown in Fig. 3. The positional and bond-orientational order of the half-skyrmions exhibits a characteristic two-dimensional melting behavior based on the Kosterlitz–Thouless–Halperin–Nelson–Young (KTHNY) scenario (33, 34, 40–43). This is confirmed by the correlation functions on the center of mass positions of the half-skyrmions, as displayed in Fig. 3 . The hexatic correlation function in Fig. 3 (see Materials and Methods for this definition) decays exponentially for (characteristic of the liquid phase), algebraically for with an exponent close to –1/4 (characteristic of the hexatic phase), and approaches a constant value for (characteristic of the solid phase). The hexatic and solid phases are distinguished by the decay of the pair correlation function along with the direction of the bond-orientation angle (see Materials and Methods for further details), which is displayed in Fig. 3. The pair correlation function decays exponentially at T = 0.44 (hexatic phase), whereas it exhibits a power-law decay at T = 0.42 with the exponent close to –1/3, which is characteristic of the stability limit of the solid phase in the two-dimensional melting theory. The real-space topological structures are presented in Fig. 3 , wherein the filtered particle configurations and corresponding Delaunay triangulation results are displayed in the upper (Fig. 3 ) and lower (Fig. 3 ) rows, respectively. In KTHNY two-dimensional melting theory (40–42), only dislocation pairs with a zero Burgers vector are excited in the solid phase. They separate into single dislocations (disclination pairs) in the hexatic phase, and the liquid phase is characterized by the disclination-unbound phase, in which isolated point defects disperse. This characteristic feature is captured in the Delaunay triangulation. Our results based on the KTHNY scenario are consistent with those found recently in magnetic skyrmion systems (33, 34).

Phase Switching and Re-Entrant Phase Transition by an External Field.

Another method for controlling topological structures is to apply an external magnetic/electric field to substances, as is often examined in magnetic skyrmion systems (7) and liquid crystals (1). In this paper, we examine magnetic field effects for our convenience, but the same conclusion holds for electric field effects because static field effects without impurities are considered. In rod-like uniaxial molecular systems, the magnetic susceptibility exhibits both (positive anisotropy) and (negative anisotropy), where () is the susceptibility of a molecule parallel (perpendicular) to the molecular long axis, respectively (1). For positive (negative) anisotropy, the orientation of the molecular long axis becomes parallel (perpendicular) to the external field. This feature is reminiscent of the perpendicular magnetic anisotropy systems in transition-metal oxides (44), in which the magnetic easy axis is tangential to the thin films. In this study, we neglect the dipole–dipole interaction between the molecules. Then, the molecular response to the external field becomes second-order, often assumed in liquid crystals (see Materials and Methods for further details) (1). We present the field response of a half-skyrmion phase in Fig. 4, where H denotes the magnitude of the field applied normal to the monolayer (see for other phases). In this work, we decrease temperature at a fixed magnetic field to examine equilibrium structures without hysteresis, though cyclic field effect, including hysteresis loop, is often studied experimentally. Starting from the half-skyrmion state shown in Fig. 4, molecules with positive (negative) anisotropy begin to align parallel (perpendicular) to the external field (see Materials and Methods, Eq. , and the following sentences). For the positive anisotropy case shown in Fig. 4, molecules at the edges of the half-skyrmions change their orientation to align parallel to the external field, resulting in the coalescence of half-skyrmions to form helical structures (H = 2). By further increasing this field, an increasing amount of particles align parallel to the external field, and a uniformly aligned state is eventually realized at H = 10. The transition from the helical state to the uniform state is similar to the field responses of cholesteric liquid crystals, wherein uniformly aligned domains are separated by sharply twisting walls during the intermediate stage (ref. [1]; see their figure 6.13). Therefore, the transition pathway for positive anisotropy demonstrated by our system is reasonable, while it sharply contrasts with magnetic skyrmion systems, in which the helical state transforms into the skyrmion state under the application of an external magnetic field.

Fig. 4.

Response to external fields under the field-cooling condition. (A–C) The external field effects on the half-skyrmion phase (A) with positive (B) and negative (C) magnetic anisotropy are shown, where the magnetic anisotropy is defined by , and H denotes the magnitude of the external magnetic field applied normal to the monolayer (Scale bars [A-C]: 10σ.). (D) The structure factors obtained at H = 0 and H = 4 with negative anisotropy imply that the topology of these structures is the same. The red circles and blue hexagons represent the Bragg spots corresponding to the nearest molecules and half-skyrmions, respectively, indicating that the distance between the half-skyrmion cores becomes smaller. , and T = 0.05 in this figure.

For negative anisotropy systems, conversely, a curious phase transformation emerges, as displayed in Fig. 4. First, the initial half-skyrmion state transforms into a helical-like striped structure with a small pitch under a weak external field (H = 2). Subsequently, another half-skyrmion state appears at H = 4, which finally transforms into a two-dimensional zig-zag structure at H = 6 (the morphology of the two-dimensional structure depends on the twist rigidity K2, as shown in ). Although the real-space structures of the half-skyrmion states in Fig. 4 look different, their topological structures exhibit similar symmetries. This is confirmed by the structure factor displayed in Fig. 4. Both states have sixfold symmetrical peaks, as highlighted by the blue hexagons, which represent the spatial correlation of the vortices. The difference between these peaks is represented by the size of the hexagons, indicating that the vortex size is reduced under the external field. This re-entrant phase transition is a unique feature of negative anisotropy systems, implying that interactions between the steric and magnetic/dielectric anisotropy result in rich phase behavior in molecular solids. As described in Fig. 2, an elastic field emerges via phase transformations, which implies that our system exhibits electro- and magneto-striction as a cross-coupling effect, which demonstrates the great potential of this system for the design of functional materials.

Phase Switching by Local Anisotropic Stress.

Finally, we present a crucial role of the emergent elastic fields on phase controllability. As revealed in Fig. 2, the helical domain and the emergent strain field exhibit the same spatial patterns. This suggests that domain orientation can be controlled by external anisotropic strain and stress. To see this, we display the mechanical response of a helical state in Fig. 5. To realize uniaxial stretching in a limited region, we apply and for the particles with , where x and y are relative displacement from the center of the system, and F0 specifies the magnitude of external force. We choose in this figure. When external mechanical stress is applied locally (the circled region in Fig. 5), both the helical pattern inside and outside the stressed region transform considerably, as displayed in Fig. 5 (the white broken circle denotes the stressed region). For the former, the helical pitch becomes parallel to the compression direction, which indicates that molecules inside the stressed region orient to the elongation direction, reducing the uniaxial stress inside the circular region. For the latter, helices change their pitch direction parallel to the dipole field. This is because both the strain and the stress field are induced outside the stressed region due to the long-range quadrupolar nature of elastic correlation, as displayed in Fig. 5 . Both the strain field and the stress field exhibit angle-dependent long-range correlation under the external stress, where they become positive along and and negative along . After the external stress is removed, the helical state inside the circular region transforms into a half-skyrmion structure, whereas the helical domains outside the stressed region remain almost unchanged, as displayed in Fig. 5. The remnant strain in relaxed states also exhibits the same spatial heterogeneity as molecular orientation (Fig. 5), whereas the stress is localized (Fig. 5). Thus, we successfully control topological phases by applying mechanical stress.

Fig. 5.

Response to local anisotropic stress. (A) The helical state for initial unstressed condition. The circle at the center and the arrows represent the area and the orientation of applied anisotropic stress in B. (B) Molecular configuration under the stress, where the molecular configuration is perturbed globally. (C) The molecular configuration after stress removal. The half-skyrmion structure is maintained at the central area after the stress is relaxed. (D–F) Uniaxial component of the emergent strain field for initial state (D), stressed state (E), and relaxed state (F). (G–I) Uniaxial component of the emergent stress field for initial state (G), stressed state (H), and relaxed state (I). , and T = 0.2 in this figure (Scale bars: 10σ.).

Discussion and Summary

In this study, we succeeded in controlling the topological phase transition between the helical and the half-skyrmion phases. This transition is reversible without plastic deformation. This reversibility is attributed to the fact that the effective aspect ratio of our molecule is close to unity: Lattice distortion induced by molecular rotation is not large, so that crystal defects do not form. When particles with a large aspect ratio, such as fd-viruses (45) and cellulose nanocrystals (46), are utilized in a solid phase, particle rotation induces plastic deformation, which results in irreversible phase transformation. Therefore, a particle should have a shape deformed from the isotropic spherical shape only slightly, and indeed it can be synthesized experimentally (47–53). A densely packed solid thin film of these particles will be a candidate to exhibit phase transition between the helical and the half-skyrmion phase. Furthermore, it is possible to measure elastic heterogeneity by atomic force microscopy for molecular and nanoparticle crystals (54) and by confocal microscopy for crystals composed of micrometer-scale colloidal particles (55). These measurements can examine the relationship between the formation of topological phases and the emergent elastic field.

In summary, we presented a physical principle for topological phase control using material parameters and external electromagnetic fields in a model molecular solid. The most important aspect of our model is the emergent elastic field produced by the interactions between the molecular steric anisotropy and twisting interactions. We also presented a method for controlling the phases of this system using external fields, which is achieved due to the competition between the emergent elastic field and the paramagnetic/paraelectric response. We identified a topologically re-entrant phase transition induced by an external field in negative anisotropy systems. These results are summarized in Table 1. In this paper, we focus our attention on a quasi-two-dimensional monolayer for simplicity. Thickness dependence of the stability of the helical and half-skyrmion structures should be studied in the future (see also ). We also considered the electromagnetic effects in a simple manner by neglecting the electromagnetic interactions between the induced dipoles. The inclusion of dipole–dipole interaction results in controllable polar orders with large mechanical responses (56). This cross-coupling may also be a key to control emergent elastic fields associated with topological phase transitions induced by anisotropic mechanical stresses (32, 57).

Materials and Methods

Molecular Dynamics Model.

We construct a simple molecular model exhibiting phase transitions into helical and half-skyrmion phases by applying conventional knowledge in liquid crystals (1). The potential energy of our system is given bywhere ϵ and σ denote the characteristic energy and length in our model; and are the absolute value and unit vector of the intermolecular displacement, respectively; and denotes the molecular orientation of the uniaxial molecules. A represents the symmetric steric repulsion used to mimic spheroidal molecules in a condensed phase (58, 59), in which η represents the magnitude of the anisotropy. For small η values, a particle can be regarded to have an ellipsoid-like shape, with the aspect ratio being . We choose η  =  2 in this work so that . B represents twisting interactions arising from molecular chirality that adjacent molecules favor to align with a twist, where the favored twist angle and twist rigidity are given by q0 and K2, respectively (60, 61). This is the discretized form of the twist Frank energy in liquid crystal theory (1), wherein the bilinear term is the same as the Dzyaloshinskii–Moriya interaction in magnetic systems ( with ) (7). In this paper, we assume that the intermolecular interaction has only the short-range repulsive steric term in order to examine the role of the twist interaction and the steric anisotropy in the simplest manner. We note that inclusion of van der Waals attractive interaction does not change the qualitative features of the formation of the half-skyrmion and helical phases.

represents the external field effects of the molecular orientation, which is defined aswhere is the external field, and denotes the anisotropic susceptibility. Here, () is the susceptibility of a molecule that is parallel (perpendicular) to the molecular long axis. A molecule is oriented parallel (perpendicular) to the external field for positive (negative) values. This form is often adopted for liquid crystals (1). A small number of molecules are known to exhibit negative magnetic anisotropy (62), whereas various liquid crystalline molecules that exhibit negative dielectric anisotropy are found (1). In numerical simulations, we normalize by redefining as .

arises from the confinement due to the walls of the system. In this study, we assume a monolayer of spheroid-like molecules without surface anchoring to examine the impact of molecular twisting on the molecular configuration of the system in a simple manner. We therefore assumewhere z is the z-coordinate of the i-th molecule, and L is the distance between the two walls. By utilizing , it can be easily confirmed that the molecules form a monolayer without undulation along the z axis (see for larger L values). in this monolayer geometry, and hence A  =  0 when two adjacent molecules align perpendicular to the intermolecular displacement , while B  =  0 when two molecules exhibit the twisting angle .

The equations of motion with respect to the molecular position and orientation are given bywhere m is the molecular mass, I is the moment of inertia with respect to the molecular long axis, and is the unit tensor. By using , Eq. may be rewritten as

The molecular torque is given by . From Eq. , the ground state of molecular orientation is given by , where U is the first term in Eq. , representing intermolecular interaction. This equation indicates the balance between the molecular field (the left-hand side) and the torque arising from the external field (the right-hand side), where the stable root of the latter is () for (), consistent with the continuum description in liquid crystals (1). Thus, in Fig. 4 , molecular orientation aligns parallel and perpendicular to strong external fields, respectively.

Eqs. and are integrated during the time evolution under the canonical ensemble using the Nosé–Hoover thermostat (63). In this study, we assume densely packed systems. The volume fraction is set to 0.3, where N is the number of the particles, and V is the system volume. Temperature T is noted in the unit of , where is the Boltzmann constant. The melting temperature is much higher than the orientational transition temperature in this condition. First, we prepare the liquid state at T = 4 and anneal slowly () to realize two-dimensional hexagonal crystals at T = 2 (the melting temperature is ). At this temperature, the lattice constant becomes and is common for all the simulation parameters examined in Fig. 1. Then, we lower the temperature to T = 0.01 with the cooling rate of to examine the ground-state structure in Fig. 1. The stress tensor in Figs. 2 and 5 is noted in the unit of . We adopt a periodic boundary condition in the x and y directions. For most numerical simulations, we choose N = 4, 000, except for Fig. 3 (N = 64, 000 to examine long-range correlations) and (N = 16, 000) and (N = 64, 000).

Structure Factor.

As displayed in Figs. 1, 3, 4, and 5 in the main text, the topology of a structure is characterized by using the arrangement of the molecules that are aligned perpendicular to the z axis. Therefore, we define a filtered density distribution to determine the topological structure, where for and otherwise. Accordingly, we calculate the corresponding structure factor , where , and i is the imaginary unit.

Emergent Elastic Fields.

In Figs. 2 and 5 and , we present in-plane local strain and stress fields. The strain tensor for each particle is defined bywhere is the coordination number of the i-th particle, in which the cutoff length is the first minimum of the radial distribution function. This sum is obtained over the coordinated particles, and is the first maximum of the radial distribution function required to normalize the strain tensor such that at the ground state. The volumetric, uniaxial, and shear strains are defined as , , and , respectively. The strain field displayed in Figs. 2 , 5 , and is obtained by coarse-graining the strain tensor as

where is the weight function.

The local stress field is calculated using the Irving–Kirkwood formula as (64)where the first and second terms denote the kinetic and interaction (configuration) terms, respectively. is the translational velocity of the particle, and represents interparticle forces arising from the pair interaction term (the first term in Eq. . In Figs. 2 , 5 , and , we also apply coarse-graining to the stress field using the weight function .

Identification of Half-Skyrmions.

In order to carry out the structural analysis among the half-skyrmions, we identify the centers of the half-skyrmions. We consider a molecule as a member of a half-skyrmion when . Among such molecules, two molecules are defined to belong to the same cluster if their distance is smaller than 1.5 (nearest-neighbor particles). Then, the position of a half-skyrmion is defined by , where is the number of particles in the same cluster.

Correlation Functions and Delaunay Triangulation on Half-Skyrmions.

The hexatic correlation function (the angular average is used in this study) on half-skyrmions are calculated from , which is defined aswhere the bracket denotes the angular, space, and sample averaging processes; and , i is the imaginary unit, and θ is the bond-orientational angle between the displacement vector and x axis. Using the radial distribution function on half-skyrmions, the scaled hexatic correlation function is displayed in Fig. 3. In Fig. 3, we show the pair correlation function along with the direction of the bond-orientational angle. We then perform Delaunay triangulation to obtain Fig. 3 .