Theory of hybrid defects, with coupled orientational order parameters, on flat and curved surfaces

Abstract

Many physical systems involve two types of orientational order, which are coupled together. For example, ferroelectric nematic liquid crystals have coupled polar and nematic order, and tilted hexatic phases have coupled polar and hexatic order. In these systems, defect structures can be quite complex. Here, we investigate phases with two types of two-dimensional orientational order, m-atic and n-atic, where m and n are two distinct integers. We simulate these phases in a flat disk with strong radial anchoring, and on a spherical surface, because both of these geometries require the presence of defects. If the coupling between the two types of order is weak, then the defects are connected by a network of diffuse walls, and the system forms a stable domain structure. As the coupling increases, the domain walls become sharper and shorter. For very strong coupling, the higher-order defects merge into the lower-order defects, forming stretched defect cores.

---

1 Introduction

Topological defects are fundamental features of ordered phases in condensed matter physics, as discussed in several review articles and textbooks.^1–6 These defects play essential roles in the structure, dynamics, and statistical mechanics of hard and soft matter. Furthermore, observation of defects is often a useful method for recognizing ordered phases.

Many physical systems have different types of orientational order. This order may be polar, with one-fold rotational symmetry, as in ferromagnets. It may be nematic, with two-fold rotational symmetry, as in the most common phases of liquid crystals. It may also have a higher-order rotational symmetry, such as a tetratic phase with four-fold symmetry, or a hexatic phase with six-fold symmetry.⁷ In general, researchers use the term n-atic to refer to orientational order with n-fold rotational symmetry. In any phase with n-atic order, topological defects are points in two dimensions (2D). Around a point defect, the orientational order rotates through an angle of 2π/n radians.

In certain cases, a physical system might have more than one type of orientational order. For example, a liquid crystal might have strong nematic order (all molecules aligned up or down along a certain axis), along with weak polar order (51% of molecules pointing up and 49% of molecules pointing down along that axis). Alternatively, a thin film of liquid crystal might have hexatic order (six-fold symmetry in the orientations of intermolecular bonds), along with polar order (one-fold symmetry in the directions of molecular tilt).^8–10 In general, there might be a combination of m-atic and n-atic order, where m and n are two distinct integers.

In a phase with two types of orientational order, what happens to the topological defects? This question was addressed in two important papers in the 1980's. In a theoretical paper, Lee and Grinstein considered the general concept of a phase with 2D polar and nematic order.¹¹ They showed that this system forms “string defects,” each consisting of two ±1/2 point defects in the nematic order, connected by a domain wall in the polar order. Around the same time, Dierker, Pindak, and Meyer investigated liquid-crystal films with 2D polar (tilt) and hexatic order.¹² They found experimentally and explained theoretically that these films form “star defects,” each with five 1/6 point defects in the hexatic order at the tips of a star and one more at the center, connected by tilt domain walls. A further theoretical study considered how the star defects should evolve through transitions between phases with different couplings between tilt and hexatic order.¹⁰ The general concept of string defects and star defects is reviewed in Section 7.2 of the textbook.⁶

In this paper, we further explore the structure of topological defects in phases with different types of orientational order that are coupled together. We will refer to this general class of defects as “hybrid defects.”¹³ There are several motivations for this study:

First, the ferroelectric nematic phase has recently been discovered in newly synthesized liquid crystals.^14–18 Unlike the conventional nematic phase, the ferroelectric nematic phase has polar order of the molecules, leading to an electrostatic polarization. Hence, these materials have a combination of polar and nematic order, similar to the combination studied theoretically by Lee and Grinstein,¹¹ although the electrostatic interaction is more complex than the energy function in their model. Experiments have observed topological defects in ferroelectric nematic liquid crystals, and they have interesting extended structures, involving both polar and nematic order.^13,19,20

Second, simulations by Glotzer and collaborators have studied how particles pack to form a colloidal crystal on the surface of a sphere.²¹ The spherical geometry is important because it requires the presence of topological defects with a total topological charge of +2. In these simulations, the defect structure depends on the shape of the particles. If the particles are rounded cubes, they pack in a lattice with four-fold symmetry. The tetratic orientational order leads to eight defects, each with topological charge +1/4, around the sphere. By contrast, if the particles are rounded tetrahedra, they form a more complex woven motif, which we might interpret as a combination of nematic and tetratic orientational order. In that case, the simulations show a network of extended defects and domain walls around the sphere.

Third, models of active matter often exhibit a combination of polar and nematic order. In some cases, this combination leads to string defects, which have two +1/2 defects in the nematic order connected by a domain wall in the polar order.^22,23 The presence of these defects controls the dynamic evolution of the active material.

Fourth, researchers have recently analyzed the arrangement of biological cells in epithelial cell monolayers, which can exhibit active flow.^24–28 These analyses show that the arrangement has a combination of different types of n-atic orientational order, particularly of nematic and hexatic order. This combination has at least the potential to generate hybrid defects.

Fifth, Vafa and Doostmohammadi have recently published a theory for string defects in “nematopolar” phases with combined polar and nematic order.²⁹ They argue that these defects lead to a new universality class for the transition between nematic and nematopolar phases.

Here, we consider the general phenomenon of hybrid defects involving two distinct types of orientational order, m-atic and n-atic. We explore four combinations that commonly occur: (a) polar (m = 1) and nematic (n = 2), (b) nematic (m = 2) and tetratic (n = 4), (c) polar (m = 1) and hexatic (n = 6), and (d) tetratic (m = 4) and hexatic (n = 6). Note that n is a multiple of m in the first three cases, but not in the fourth case. We use a lattice model to simulate the combined orientational order in two geometries that require defects: a disk with strong radial anchoring (which requires a total topological charge of +1), and the surface of a sphere (which requires a total topological charge of +2). Several previous studies have investigated a single type of orientational order on a sphere, and have found the configuration of topological defects;^30–33 now we consider the corresponding problem with two types of orientational order.

The main results of this study are as follows. If the coupling between m-atic and n-atic order is weak, then the system forms a network of point defects connected by domain walls. The relationship between m-atic and n-atic order is fixed within each domain, and it shifts across each wall. As the coupling increases, the domain walls become sharper and shorter. The behavior for larger coupling depends on whether n is a multiple of m. If so, the system forms clusters with one defect in the m-atic order and n/m defects in the n-atic order; each cluster has a topological charge of 1/m. Eventually, for very strong coupling, the n-atic defects merge into the m-atic defect, so that the cluster becomes a single defect with a stretched core. By comparison, if n is not a multiple of m, then the system remains as a global network of domains and domain walls. We recognize that n is a multiple of m in all current experiments. It might be possible to make a system where n is not a multiple of m by designing colloidal particles with strong m-fold and n-fold-symmetric features in their shape, similar to the simulations of Glotzer and collaborators.²¹

2 Phases with a single type of orientational order

To develop the model, we first simulate phases with a single type of m-atic orientational order. We consider two geometries that require the presence of defects. The first geometry is a flat disk with strong radial anchoring on the edge. This anchoring condition requires the total topological charge of +1 inside the disk. The second geometry is a spherical surface, with orientational order in the local tangent plane. On a sphere, topology requires a total topological charge of +2.

For the flat disk, we set up a lattice model by constructing a finite-element mesh, using the mesh algorithm pygmsh, a Python interface to Gmsh. The circular disk geometry is defined parametrically by the radius R = 1 of the disk, and a triangular discretization is generated using Gmsh's built-in meshing algorithms. The mesh size parameter a = 0.1 controls the element size of the mesh. On each site i of the mesh, we define an orientation by the unit vector ŝ_i = (cos θ_i, sin θ_i). For a single type of m-atic orientational order, we use the Hamiltonian

In this expression, the sum is over nearest-neighbor sites i and j in the mesh, and the coefficient J_m > 0 represents the strength of the m-atic interaction. On edge sites, the orientation is constrained to point radially outward, with θ_i = tan⁻¹(y_i/x_i). We perform Monte Carlo simulations of this lattice model. In this procedure, we begin at a high temperature T, and gradually reduce T to zero. At each step, we make a trial change of θ_i at one site, and accept or reject this change following the Metropolis algorithm. The trial change may be either a small rotation or a large rotation of 2π/m (except for polar order m = 1, where the large rotation is π), so that the orientation is free to explore all of the minima of the interaction potential.

Fig. 1 shows results for polar (m = 1), nematic (m = 2), tetratic (m = 4), and hexatic (m = 6) interactions. In this figure, we visualize the results by drawing m arrows at each mesh site, corresponding to the orientations θ_i, θ_i + 2π/m, …, θ_i + 2π(m − 1)/m. This visualization is appropriate because all m of those arrows represent equivalent directions, all with the same energy, showing the m-fold symmetry. We also color each triangular cell of the mesh to indicate the local energy density.

Fig. 1 Disk with a single m-atic type of orientational order. (a) Polar. (b) Nematic. (c) Tetratic. (d) Hexatic. At each site of the mesh, the m arrows indicate the local orientational order with m-fold symmetry. On each triangular cell, the color (yellow to red) indicates the local energy density (average of the three bond energies around the cell).

In each case, the simulated annealing algorithm reaches an energy minimum, in which the system has well-defined m-atic order everywhere except at point defects. As expected, the disk has m defects, each with topological charge +1/m, for a total topological charge of +1. Around each defect, the orientation of the arrows rotates through an angle of 2π/m. The defect cores have high energy density, as indicated by the red color in the figure.

For simulations on a sphere, we follow a similar procedure, but we must be careful with the coordinate system. The 2D spherical coordinate system has singularities at the north and south poles, and these singularities might bias the defect structure on the sphere. To avoid any such bias, we explicitly work in 3D Cartesian coordinates, generalizing the approach from previous studies of polar and nematic order on a sphere.^32,33 We construct a finite-element mesh on a spherical surface with radius 1, again using the mesh algorithm in the Python library, trimesh package. The function trimesh.creation.icosphere() generates a sphere using the icosphere method, which starts from an icosahedron, a regular polyhedron with 20 faces, and iteratively subdivides the triangles to create a more spherical shape. On each site i of the mesh, with position r̂_i = (x_i,y_i,z_i), we define an orientation by the unit vector ŝ_i = (s_ix,s_iy,s_iz). We require that all of the vectors are normalized to |ŝ_i| = 1, and maintain this normalization at each step of the simulation. To force these vectors into the local tangent plane, we include a term in the Hamiltonian of the form

with a large coefficient C > 0. To express the m-atic interaction in terms of the 3D vectors, we define θ_ij as the angle between ŝ_i and ŝ_j, so that cos θ_ij = ŝ_i·ŝ_j. The m-atic interaction then becomesUsing trigonometric identities, it can be transformed intoWe perform Monte Carlo simulations of the total lattice Hamiltonian , beginning at high temperature T and gradually reducing T to zero. At each step, we make a trial change of ŝ_i at one site, and accept or reject this change following the Metropolis algorithm. The trial change may be either a small rotation about a random axis or a large rotation of 2π/m about r̂_i (except for polar order m = 1, where the large rotation is π), so that the orientation is free to explore all of the minima of Hamiltonian.

Fig. 2(a) shows results for polar order, visualized on the sphere. On each mesh site, the arrow represents the unit vector ŝ_i. On each triangular cell of the mesh, the color indicates the local energy density. The polar order is well-defined everywhere except at two point defects, each with topological charge +1, for a total topological charge of +2. One of the defects is clearly visible in the figure, as shown by the red color for high energy density. The other defect cannot be seen because it is on the back of the sphere.

Fig. 2 Sphere with a single m-atic type of orientational order. (a) Polar, shown on a sphere, only one defect visible. (b) Polar, projected onto a cylinder, only one defect visible. (c) Polar, Mercator projection onto flat plane, both defects visible. (d) Nematic, Mercator projection, all four defects visible. (e) Tetratic, Mercator projection, all eight defects visible. (f) Hexatic, Mercator projection, all twelve defects visible. As in Fig. 1, the arrows indicate the local orientational order with m-fold symmetry, and the color indicates the local energy density.

For a more global view of the entire defect configuration, we must modify our visualization. First, we project the sphere onto a cylinder, as shown in Fig. 2(b). Next, we unroll the cylinder into the flat plane, as in Fig. 2(c). This construction is the Mercator projection, which is well-known in map-making. The same construction was recently used to visualize colloidal crystals on a sphere.²¹ For our current study, the main advantage of the Mercator projection is that it shows the orientational order on the entire sphere. In the system with polar order, it shows both of the defects. For that reason, we will use the Mercator projection of the sphere for the rest of this paper. Of course, we must recognize that the Mercator projection distorts the structure near the north and south poles.

Fig. 2(d)–(f) show similar results for nematic (m = 2), tetratic (m = 4), and hexatic (m = 6) orientational order. In each case, we draw m arrows at each mesh site, corresponding to the orientation ŝ_i rotated through angles 0, 2π/m, …, 2π(m − 1)/m about the local normal vector r̂_i. This visualization is appropriate because all m of those arrows represent equivalent directions in the local tangent plane. We can see that the sphere has well-defined m-atic order everywhere except at point defects. There are 2m defects, each with topological charge +1/m, for a total topological charge of +2.

3 Phases with two types of orientational order

We now generalize the model to simulate phases with two types of orientational order, m-atic and n-atic, which interact with each other. In this generalization, we combine the Hamiltonian functions from the previous section, corresponding to m-atic and n-atic order. For a flat disk, we add the functions H_m and H_n from eqn (1), to obtainSimilarly, for a sphere, we add the functions of eqn (2)–(4) to obtain . One might naively think that the Hamiltonian of eqn (5) is just two independent terms that are added together. However, it implicitly creates a coupling between m-atic and n-atic order, because the two terms involve the same angles θ_i and θ_j on neighboring sites. This coupling can be understood both physically and mathematically:

Physically, the coupling arises from the shape of the potential energy, which is plotted in Fig. 3. First, consider the case of polar (m = 1) and nematic (n = 2) order, as in part (a). This potential energy has its deepest minimum at θ_i − θ_j = 0 (mod 2π), and its second-deepest minimum at θ_i − θ_j = π (mod 2π). If all nearest-neighbor pairs are near the deepest or second-deepest minimum, the system has nematic order but not polar order. If all pairs are near the deepest minimum only, the system has both nematic and polar order, and these two types of order are aligned with each other. In the ground state, nematic order is along some axis ±n̂, and polar order is along +n̂ only. Hence, these two types of order are coupled together.

Fig. 3 Plots of the potential energy H_m,n as a function of the relative orientation of neighboring sites, for the four examples considered in this paper. (a) Polar and nematic. (b) Nematic and tetratic. (c) Polar and hexatic. (d) Tetratic and hexatic. In each case, the coefficients are J_m = J_n = 1.

The same consideration applies in part (b), for nematic (m = 2) and tetratic (n = 4) order, and in part (c), for polar (m = 1) and hexatic (n = 6) order. In these cases where n is a multiple of m, the system can have n-atic order only, or it can have both n-atic and m-atic order, with favored alignment between the two types of order.

Part (d), for tetratic (m = 4) and hexatic (n = 6) order, is more complex because n is not a multiple of m. In that case, the values of J₄ and J₆ determine which minima are deepest, and which are second-deepest, and hence which phases can occur. However, whenever the system has both tetratic and hexatic order, the potential energy favors alignment between these two types of order. We will discuss that case further in Section 3.4 below.

Mathematically, we can develop a mean-field theory for the Hamiltonian of eqn (5), This theory is expressed in terms of complex order parameters Ψ_m = ψ_me^imϕ_m = 〈e^imθ〉 and Ψ_n = ψ_ne^inϕ_n = 〈e^inθ〉. Under rotation through an angle γ, these order parameters transform as Ψ_m → e^imγΨ_m, , and likewise for Ψ_n. In the free energy, all terms must be invariant under rotation. Hence, we expect coupling terms of the form

where a, b, c, and d are nonnegative integers with

m(a − b) + n(c − d) = 0. (7)

Indeed, we can construct the mean-field free energy F = E − TS = 〈H_m,n〉 + k_BT〈log ρ〉, where ρ(θ) is the probability distribution function for the angle θ. When we expand it as a power series in the order parameters, we obtain terms of the form in eqn (6). For polar and nematic order, the lowest-order coupling term has the form ψ₁²ψ₂ cos[2(ϕ₁ − ϕ₂)]. As a result, the strength of the angular coupling depends on the order parameter magnitudes as ψ₁²ψ₂. For nematic and tetratic order, the same argument gives a coupling of the form ψ₂²ψ₄ cos[4(ϕ₂ − ϕ₄)]. For polar and hexatic order, the analogous term is ψ₁⁶ψ₆ cos[6(ϕ₁ − ϕ₆)]. For tetratic and hexatic order, it is ψ₄³ψ₆² cos[12(ϕ₄ − ϕ₆)].

This mean-field theory is a generalization of the theory for nematic and tetratic order in a previous paper by our group.³⁴ We will not use in the rest of this article. We only mention it here in order to show mathematically that the Hamiltonian of eqn (5) creates a coupling between m-atic and n-atic order.

3.1 Polar and nematic

As a first example, we begin with the combination of polar and nematic order. This is the simplest case, which was originally considered by Lee and Grinstein.¹¹ It is also the case that is most relevant for ferroelectric nematic liquid crystals.

Fig. 4(a)–(d) shows a series of simulations of the Hamiltonian H_1,2 on a disk with radial boundary conditions. The nematic interaction strength is fixed at J₂ = 2, and the polar interaction strength increases over the range J₁ = 0.2, 0.5, 0.8, and 1.5. The arrows indicate the local orientational order. In most of the disk, away from any defects, the phase has polar order, as shown by alignment of neighboring arrows. It also has nematic order, shown by the alignment of neighboring arrow shafts, neglecting the arrow heads. In those respects, it is similar to the polar phase shown in Fig. 1(a). The difference occurs in the defects.

Fig. 4 Phase with combined polar and nematic order. (a)–(d) Simulation on disk with J₂ = 2 and J₁ = 0.2, 0.5, 0.8, and 1.5, respectively. (e)–(h) Mercator projection of sphere with C = 5, J₂ = 1, and J₁ = 0.2, 0.5, 0.8, and 1.5.

In Fig. 4(a), where the polar coefficient J₁ is smallest, the phase has two point defects with high energy density, indicated by the red color. These two points are topological defects in the nematic order, with topological charge of +1/2 each. They are connected by a domain wall with moderate energy density, indicated by yellow. Across the domain wall, the polar order reverses, while the nematic order remains unchanged. In Fig. 4(b)–(d), as the polar coefficient J₁ increases, the domain wall becomes shorter, and the two point defects become closer together. In all cases, the entire structure (two point defects connected by a domain wall) can be regarded as a single “string defect,” with topological charge of +1, as required by the boundary conditions.

Fig. 4(e)–(h) shows analogous simulations on a sphere, illustrated in the Mercator projection. Here, the normal interaction is fixed at C = 5, the nematic interaction at J₂ = 1, and the polar interaction increases over the range J₁ = 0.2, 0.5, 0.8, and 1.5. The results are generally similar to the simulations on the disk. Once again, the system forms string defects of topological charge +1, each consisting of two +1/2 nematic defects connected by a domain wall. The sphere has two string defects, for a total topological charge of +2, as required by topology. As the polar coefficient J₁ increases, the domain walls become shorter and the pairs of nematic point defects become closer together. In the limit of strong polar interaction, each pair of nematic defects effectively merges into a single polar defect with an extended core.

These simulation results are consistent with the theory of string defects, as reviewed in Section 7.2 of the textbook.⁶ The free energy of a string defect has two parts. First, there is the repulsion between the two nematic defects, which is 2πKq₁q₂ log(R/L), where K is the orientational elastic constant, q₁ and q₂ are the separate topological charges (q₁ = q₂ = +1/2), R is the system size, and L is the distance between the defects. Second, there is the free energy of the domain wall itself, which is ε_WL, where ε_W is the line tension (energy per length). The total is

Minimizing the total free energy over L gives the equilibrium lengthHere, the effective line tension ε_W is the difference between the energies for antiparallel and parallel alignment of neighboring arrows, normalized by the typical spacing a between mesh points. We estimate ε_W ≈ 2J₁/a. The effective elastic constant K is the second derivative of the potential H_1,2 at the minimum. We estimate K ≈ (J₁ + 4J₂). Hence, the equilibrium length of a string defect should beFor J₁ = 0, this predicted length is infinite, meaning that the two nematic defects are free to move apart, not connected by a domain wall. As J₁ increases, this predicted length decreases, so that the nematic defects are more tightly bound by the domain wall, and the string becomes shorter. This trend agrees with the simulations.

3.2 Nematic and tetratic

Next, we consider the combination of nematic and tetratic order. Mathematically, this problem is equivalent to the combination of polar and nematic order, but with the angle 2θ substituted in place of θ. For that reason, we expect to see string defects, similar to the string defects seen in the previous case. Hence, the main purpose of the new simulations is visualize defects in nematic and tetratic order, so that they can be compared with other research.

Fig. 5 presents simulations of a phase with nematic and tetratic order. Parts (a)–(d) show the phase on a disk with radial anchoring, and parts (e)–(h) show the phase on a sphere (in the Mercator projection). At each mesh site i, the orientational order is represented by a + symbol, decorated in black and blue. The colorless + symbol shows the tetratic order, and the black/blue decoration shows the nematic order. On the disk, each + symbol is constructed by drawing two black arrows at the angles θ_i and θ_i + π, and two blue arrows at the angles θ_i ± π/2. Likewise, on the sphere, the black arrows are drawn along the vector ŝ_i, rotated by 0 and π about the local normal. The blue arrows are drawn along ŝ_i, rotated by ±π/2 about the local normal. The tetratic interaction favors alignment of neighboring + symbols, regardless of color, as if they were four-fold symmetric objects. The nematic interaction favors alignment of neighboring black with black, and neighboring blue with blue.

Fig. 5 Phase with combined nematic and tetratic order. (a)–(d) Simulation on disk with J₄ = 2 and J₂ = 0.2, 0.4, 0.8, and 1, respectively. (e)–(h) Mercator projection of sphere with C = 5, J₄ = 0.25, and J₂ = 0.05, 0.1, 0.2, and 0.5.

As in the previous case of polar and nematic order, the simulations of nematic and tetratic order show point defects with high energy density (red color), as well as domain walls with intermediate energy density (yellow color). The point defects are topological defects in the tetratic order, with topological charge of +1/4 each. The domain walls are lines where the nematic order reverses, so that the black and blue colors are exchanged, while the tetratic order remains the same.

On the disk, in Fig. 5(a)–(d), the tetratic interaction is fixed at J₄ = 2, and the nematic interaction increases from J₂ = 0.2, 0.4, 0.8, to 1. On the sphere, in Fig. 5(e)–(h), the tetratic interaction is J₄ = 0.25, and the nematic interaction increases from J₂ = 0.05, 0.1, 0.2, to 0.5. For weak nematic interactions, the domain walls are only faint yellow, meaning that they have low energy, and hence they are quite long. For the smallest ratio J₂/J₄, in Fig. 5(e), we can see that the domain walls form a network connecting the point defects around the sphere. On both the disk and the sphere, the domain walls occasionally have kinks, which are presumably artifacts of the underlying mesh.

As the nematic interaction increases, the domain walls become brighter yellow, meaning higher energy, and hence they grow shorter. They bind the point defects into pairs, which can be regarded as string defects. Each string defect has a total topological charge of +1/2. The disk has two string defects, for a total topological charge of +1, as required by the boundary conditions. The sphere has four string defects, with a total topological charge of +2, as required by topology. For the highest ratio J₂/J₄, in Fig. 5(h), each string shrinks into a single nematic defect with an extended core.

Our results can be compared with recent simulations of colloidal crystals on spherical interfaces.²¹ That study uses a very different computational approach from our current paper; it is based on the packing of hard particles with different types of polyhedral shapes. In their results, the structure and defects depend on the particle shape. Cubical particles form a square lattice, which has tetratic orientational order as well as crystalline positional order. Around the sphere, this lattice has eight defects of topological charge +1/4 each, similar to our simulation of tetratic order in Fig. 2(e). By contrast, tetrahedral particles form a more complex woven motif, with edges of the tetrahedra oriented inward and outward on the sphere. In that case, the overall structure on the sphere has many extended defects, which are mostly connected with each other to form a network. It appears somewhat similar to our simulation of strong tetratic order with weak nematic order in Fig. 5(e). We speculate that their results might also be understood in terms of strong tetratic order with weak nematic order. The woven motif in their simulations certainly has four-fold symmetry, like a tetratic phase. This four-fold symmetry may be broken down to two-fold, like a nematic phase, by forcing the woven motif into a curved geometry. The combination of these two types of order may lead to a structure of point defects and domain walls, as we have found.

3.3 Polar and hexatic

For a further example of the same type of phenomenon, consider the combination of polar and hexatic order. One physical realization of this combination is in thin films of smectic liquid crystals.^8–10,12 In that case, the hexatic order corresponds to the orientations of bonds between neighboring molecules, which have six-fold symmetry. The polar order corresponds to the direction of molecular tilt, which has only one-fold symmetry. These two types of order are coupled, because the tilt generally has a preferred orientation with respect to the local bonds.

Fig. 6 shows simulations of the phase with polar and hexatic order, first on a disk with radial anchoring (parts (a) and (b)) and then on a sphere in the Mercator projection (parts (c) and (d)). At each mesh site i, the orientational order is represented by a symbol with six arms, decorated by a single arrow. The six arms show the hexatic order, and the arrow shows the polar order. For the disk, the arms are drawn at the angles θ_i, θ_i + π/3, …, θ_i + 5π/3, and the arrow at θ_i only. For the sphere, the arms are drawn along the vector ŝ_i, rotated by 0, π/3, …, 5π/3 about the local normal, and the arrow along ŝ_i only. The hexatic interaction favors alignment of neighboring symbols, regardless of the arrows, as if they were six-fold symmetric objects. The polar interaction favors alignment of neighboring arrows.

Fig. 6 Phase with combined polar and hexatic order. Simulations on disk: (a) J₆ = 2, J₁ = 0.5. (b) J₆ = 0.5, J₁ = 2. Simulations on sphere: (c) C = 10, J₆ = 1, J₁ = 0.8. (d) C = 5, J₆ = 0.25, J₁ = 1.

When the hexatic coefficient J₆ is strong and the polar coefficient J₁ is weak, as in Fig. 6(a) and (c), the simulations show point defects in the hexatic order, each with topological charge of +1/6, which are visible as red points of high energy density. The disk has six of these defects, for a total topological charge of +1, and the sphere has twelve, for a total of +2. The point defects are connected by domain walls, which are visible as yellow lines. Across each domain wall, the polar order rotates through an angle of π/3, so that it goes from one of the six directions of hexatic order to another. The domain walls are fairly low in energy, so they can be quite long, and they connect all of the point defects around the disk or sphere in a large network.

When the polar coefficient becomes stronger and the hexatic coefficient becomes weaker, as in Fig. 6(b) and (d), the domain walls become much higher in energy. As a result, the domain walls must become shorter, and hence the system forms “star defects.” Each star defect has a point defect at the center, connected to five domain walls as arms, with five more point defects at the ends of the arms. The point defect at the center is a +1/6 defect for the hexatic order, and also a +1 defect for the polar order. For that reason, it has a high energy density, indicated by the red color. The point defects at the ends of the arms are only +1/6 defects for the hexatic order. They have lower energy density, and hence only the yellow color. Each entire star defect has a topological charge of +1. There is one star defect on the disk, and two on the sphere.

The star defects in these simulations are very similar to star defects in tilted hexatic liquid crystals, which were observed experimentally and explained theoretically by Dierker, Pindak, and Meyer.¹² The main difference is that the experimental structures are more symmetrical, with straighter arms, while our simulated star defects are more disordered, with kinks in the arms. As in the previous case of nematic and tetratic order, the kinks are presumably artifacts of the underlying mesh.

By comparing the results of this section with the combinations of orientational order in the previous sections, we can see that star defects and string defects are both examples of the same general phenomenon. We will use the term “hybrid defects” to describe the class of topological structures that combine point defects and domain walls in different types of orientational order.

3.4 Tetratic and hexatic

In all of the examples so far, we have considered a combination of m-atic and n-atic order where n is a multiple of m. In those cases, the m-atic order is a more specific, lower-symmetry version of the n-atic order. In other words, the m-atic order selects m special directions among the n special directions associated with n-atic order. It reduces the symmetry by a factor of n/m, which is an integer.

As a final example, suppose that a system has a combination of tetratic (m = 4) and hexatic (n = 6) order. In this combination, n is not a multiple of m. The potential energy H_4,6 of eqn (5) and Fig. 3(d) might have either four or six minima, depending on the ratio J₄/J₆. Two of the minima are located at θ_i − θ_j = 0 and π. The other minima shift as functions of J₄/J₆.

In this case, what should we expect for the hybrid defects? To answer that question, we note that the greatest common divisor of m = 4 and n = 6 is 2, corresponding to nematic order. For that reason, nematic order is a more specific version of tetratic order, and nematic order is a more specific version of hexatic order. Hence, we speculate that the phase with combined tetratic and hexatic order should form hybrid defects with topological charge of +1/2. Each of these hybrid defects could combine two tetratic defects of charge +1/4 each, along with three hexatic defects of charge +1/6 each, connected by some combination of domain walls.

To test this speculation, we perform simulations of the model with combined tetratic and hexatic order. Fig. 7(a)–(c) shows the results on a disk with radial anchoring, and Fig. 7(d)–(f) on a sphere in the Mercator projection. At each mesh site i, the orientational order is represented by two arrows, which show the angles θ_i and θ_i + π on the disk, or equivalently the vector ŝ_i rotated by 0 and π about the local normal on the sphere. We do not draw additional decorations corresponding to tetratic or hexatic order, because they would be too cluttered on the figure.

Fig. 7 Phase with combined tetratic and hexatic order. Simulations on disk: (a) J₆ = 2, J₄ = 1. (b) J₆ = 2, J₄ = 1.5. (c) J₆ = 2, J₄ = 4. Simulations on sphere: (d) C = 50, J₆ = 1, J₄ = 0.5. (e) C = 50, J₆ = 1, J₄ = 0.75. (f) C = 50, J₆ = 1, J₄ = 2.5.

When the hexatic coefficient J₆ is large and the tetratic coefficient J₄ is small, as in Fig. 7(a) and (d), the phase is mainly hexatic, and tetratic order can be regarded as a small perturbation. Here, the red points of high energy density are hexatic defects of topological charge +1/6. There are six of these points on the disk, and twelve on the sphere. These points cannot be independent defects, because they are not compatible with the presence of some tetratic order. Hence, they must be joined by domain walls. Across each domain wall, the orientational order rotates by approximately π/3, as expected for a phase that is mainly hexatic. On the disk, in Fig. 7(a), the system forms two hybrid defects, each consisting of three +1/6 defects connected by domain walls, each with a total charge of +1/2. However, on the sphere, in Fig. 7(d), all of the defects are connected together in a network that spans the entire sphere.

When J₄ is large and J₆ is small, as in Fig. 7(c) and (f), the phase is mainly tetratic, and the hexatic order is small. The red points are now tetratic defects of topological charge +1/4. There are four of these points on the disk, and eight on the sphere. Once again, these points cannot be independent defects, because they are not compatible with the existence of hexatic order, and hence they must be joined by domain walls. Across each domain wall, the orientational order rotates by approximately π/2, as expected for a phase that is mainly tetratic. On the disk, in Fig. 7(c), the structure is organized in two hybrid defects, each involving two +1/4 defects connected by domain walls, each with a total charge of +1/2. By contrast, on the sphere, in Fig. 7(f), all of the defects are again connected together in a network over the whole sphere.

When J₄ and J₆ are of similar magnitude, as in Fig. 7(b) and (e), the results are more difficult to interpret. In these cases, the system has a large number of point defects, beyond the minimum number required by topology. These defects are connected by domain walls, which are longer than the minimum required for the connection. For both the disk and the sphere, the defects and domain walls form a network over the entire system. They do not break into hybrid defects with topological charge of +1/2.

Overall, the simulation results for tetratic and hexatic order show the same general phenomenon as the results for the other combined order parameters. As in the previous cases, the system forms point defects linked by domain walls, consistent with the topological requirements for each of the order parameters. The main difference in this new combination is that the topological structures often cover the entire geometry, instead of breaking into hybrid defects with the expected charge. It is possible that separate hybrid defects might still be the ground state, but that the simulation algorithm cannot consistently reach this ground state. We speculate that many arrangements of tetratic and hexatic order have similar energies, and the energy differences are too small to drive relaxation into a simple ground state.

4 Discussion

In this paper, we have simulated topological structures that occur when a system has two coupled orientational order parameters, m-atic and n-atic. The simulations show that this combination generally results in the formation of hybrid defects, such as the strings and stars found in earlier research.

In the simplest case where n is a multiple of m, the hybrid defects involve n/m point defects, each with topological charge +1/n. These point defects are connected by domain walls, across which the orientational order rotates through an angle of 2π/n. Each of the hybrid defects then has a total topological charge of +1/m.

If n is not a multiple of m, hybrid defects can form with a topological charge of +1/l, where l is the greatest common divisor of m and n. However, in this case, the energy landscape is more complex, and our simulations often show networks of defects extending over the entire geometry.

As we have discussed, combinations of orientational order parameters commonly occur in ferroelectric nematic liquid crystals, tilted hexatic liquid crystals, colloidal crystals on a spherical surface, and one can expect to find them in biological and active materials. Thus, our results should help to understand topological structures found in this wide range of systems.

Conflicts of interest

There are no conflicts to declare.

Data availability

Data for this article, including Python code for Monte Carlo simulations, VTK files with simulation output, and Mathematica notebooks to visualize the results, are available at https://doi.org/10.5281/zenodo.19226000.

Acknowledgements

We would like to thank H. Aharoni, S. C. Glotzer, G. N. Jones, P. W. A. Schönhöfer, and Q.-H. Wei for helpful discussions.

References

1. N. D. Mermin, Rev. Mod. Phys., 1979, 51, 591–648.
2. M. Kleman, Points, Lines and Walls: In Liquid Crystals, Magnetic Systems and Various Ordered Media, Wiley, 1983.
3. P. M. Chaikin and T. C. Lubensky, Principles of Condensed Matter Physics, Cambridge University Press, 1995.
4. M. Kleman and O. D. Lavrentovich, Soft Matter Physics: An Introduction, Springer, 2003.
5. G. P. Alexander, B. G.-G. Chen, E. A. Matsumoto and R. D. Kamien, Rev. Mod. Phys., 2012, 84, 497–514.
6. J. V. Selinger, Introduction to Topological Defects and Solitons: In Liquid Crystals, Magnets, and Related Materials, Springer, 2024.
7. D. R. Nelson and B. I. Halperin, Phys. Rev. B:Condens. Matter Mater. Phys., 1979, 19, 2457–2484.
8. D. R. Nelson and B. I. Halperin, Phys. Rev. B:Condens. Matter Mater. Phys., 1980, 21, 5312–5329.
9. J. V. Selinger and D. R. Nelson, Phys. Rev. Lett., 1988, 61, 416–419.
10. J. V. Selinger and D. R. Nelson, Phys. Rev. A:At., Mol., Opt. Phys., 1989, 39, 3135–3147.
11. D. H. Lee and G. Grinstein, Phys. Rev. Lett., 1985, 55, 541–544.
12. S. B. Dierker, R. Pindak and R. B. Meyer, Phys. Rev. Lett., 1986, 56, 1819–1822.
13. S. Yi, C. Zhou, Z. Hong, Z. Ma, M. Huang, S. Aya, R. Zhang and Q.-H. Wei, arXiv, 2025, preprint, arXiv:2511.22806 DOI:10.48550/arXiv.2511.22806.
14. H. Nishikawa, K. Shiroshita, H. Higuchi, Y. Okumura, Y. Haseba, S.-I. Yamamoto, K. Sago and H. Kikuchi, Adv. Mater., 2017, 29, 1702354.
15. X. Chen, E. Korblova, D. Dong, X. Wei, R. Shao, L. Radzihovsky, M. A. Glaser, J. E. Maclennan, D. Bedrov, D. M. Walba and N. A. Clark, Proc. Natl. Acad. Sci. U. S. A., 2020, 117, 14021–14031.
16. O. D. Lavrentovich, Proc. Natl. Acad. Sci. U. S. A., 2020, 117, 14629–14631.
17. N. Sebastián, M. Čopič and A. Mertelj, Phys. Rev. E, 2022, 106, 021001.
18. R. J. Mandle, Soft Matter, 2022, 18, 5014–5020.
19. B. Basnet, M. Rajabi, H. Wang, P. Kumari, K. Thapa, S. Paul, M. O. Lavrentovich and O. D. Lavrentovich, Nat. Commun., 2022, 13, 3932.
20. P. Kumari, B. Basnet, H. Wang and O. D. Lavrentovich, Nat. Commun., 2023, 14, 748.
21. G. N. Jones, P. W. A. Schönhöfer and S. C. Glotzer, Soft Matter, 2026, 22, 2558–2567.
22. P. K. Mishra, P. S. Mondal, P. Jena and S. Mishra, New J. Phys., 2025, 27, 074602.
23. A. Dinelli, L. Dumoulin and K. Kruse, arXiv, 2026, preprint, arXiv:2601.18307 DOI:10.48550/arXiv.2601.18307.
24. J.-M. Armengol-Collado, L. N. Carenza, J. Eckert, D. Krommydas and L. Giomi, Nat. Phys., 2023, 19, 1773–1779.
25. J. Eckert, B. Ladoux, R.-M. Mège, L. Giomi and T. Schmidt, Nat. Commun., 2023, 14, 5762.
26. J.-M. Armengol-Collado, L. N. Carenza and L. Giomi, eLife, 2024, 13, 1–24.
27. M. Chiang, A. Hopkins, B. Loewe, M. C. Marchetti and D. Marenduzzo, Proc. Natl. Acad. Sci. U. S. A., 2024, 121, 1–7.
28. L. Happel, G. Oberschelp, V. Grudtsyna, H. P. Jain, R. Sknepnek, A. Doostmohammadi and A. Voigt, eLife, 2025, 14, 1–35.
29. F. Vafa and A. Doostmohammadi, Europhys. Lett., 2025, 152, 57002.
30. T. C. Lubensky and J. Prost, J. Phys. II, 1992, 2, 371–382.
31. D. R. Nelson, Nano Lett., 2002, 2, 1125–1129.
32. H. Shin, M. J. Bowick and X. Xing, Phys. Rev. Lett., 2008, 101, 037802.
33. R. L. B. Selinger, A. Konya, A. Travesset and J. V. Selinger, J. Phys. Chem. B, 2011, 115, 13989–13993.
34. J. Geng and J. V. Selinger, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2009, 80, 011707.

---