Emergent Potts Order in a Coupled Hexatic-Nematic XY model

Addressing the nature of an unexpected smectic-A' phase in liquid crystal 54COOBC films, we perform large scale Monte Carlo simulations of a coupled hexatic-nematic XY model. The resulting finite-temperature phase diagram reveals a small region with composite Potts $\mathbb{Z}_3$ order above the vortex binding transition; this phase is characterized by relative hexatic-nematic ordering though both variables are disordered. The system develops algebraic hexatic and nematic order only at a lower temperature. This multi-step melting scenario agrees well with the experimental observations of a sharp specific heat anomaly that emerges above the onset of hexatic positional order. We therefore propose that the smectic-A' phase is characterized by composite Potts order and bound-states of fractional vortices.


I. INTRODUCTION
The appearance of a sharp specific heat signal in the multiple-step melting sequence [1][2][3] of certain freestanding liquid crystal films ( Fig. 1 (a)) remains an outstanding mystery. Such transitions are typically associated with the unbinding of topological defects [4][5][6] that are not associated with acute thermodynamic signatures. Here we revisit this unsolved problem, bringing to it modern concepts and methods that have been developed in other areas. In a minimalist model [7][8][9][10], we explore whether vortex fractionalization can lead to an emergent three-state Potts phase above a Kosterlitz-Thouless binding temperature. Supported by analytical arguments, we perform this study using large-scale Monte Carlo simulations, finding a parameter region of the phase diagram consistent with experimental observations. Optical reflectivity, electron diffraction and specific heat measurements on free-standing twodimensional (2D) films of 54COOBC (n-pentyl-4n-pentanoyloxybiphenyl-4-carboxylate) [1][2][3] provide the experimental motivation for our work. Theoretically, a two-stage melting sequence was expected in these 2D films with an intermediate hexatic phase residing between the isotropic liquid and the 2D crystalline solid. By contrast, experimentalists observed three-step melting with two phases separating the solid and the liquid states. In particular they detected a hexatic phase and then at higher temperatures a mystery intermediate liquid phase with no long-range orientational order. The experimentalists referred to this unexpected "hidden order" (HO) in 2D films of 54COOBC as smectic-A (Sm-A ) phase. It is sandwiched between a (disordered) isotropic smectic-A (Sm-A) phase at higher temperature * porth@iastate.edu and a hexatic (Hex-B) phase with bond orientational quasi-long-range order (QLRO) at lower temperature (see Fig. 1 (a)). At even lower temperatures, the system develops positional QLRO in the hexagonal crystalline phase (Cry-B). The transition from disordered Sm-A into HO Sm-A phase is characterized by a pronounced specific heat anomaly that is in striking disagreement with the broad features predicted by the 2D melting theory of Kosterlitz, Thouless, Halperin, Nelson and Young (KTHNY) [4][5][6][11][12][13][14]. The scaling exponent associated with this specific heat was reported to be α 54COOBC = 0.30 ± 0.07 [1], suggesting that the HO Sm-A phase has 3-state Potts order at temperatures above the conventional vortex binding transition.
Theoretically this problem can be studied by a min- 2) (center), where T KT,hex and TKT,nem are the binding temperatures of q ϑ = 1/3 (blue tri-arrows denote the hexatic ϑ variable) and qϕ = 1 (red arrows denote the nematic ϕ variable) vortices, respectively. Here T3 refers to a three-state Potts transition that is well-established at temperatures below T KT,hex , but could persist above TKT,nem. The disordered phase is characterized in (a) by free hexatic vortices (blue dot), which correspond to fractional nematic vortices, where qϕ = 1/3 due to the presence of a Potts domain wall (in black) where ∆n = ±1. For J2 J6, these hexatic vortice then become bound in neutral pairs as one lowers the temperature below T KT,hex . The hexatic ordered phase is presented pictorially in (b), where the hexatic variable (blue) is still ordered, but the lack of a single Potts domain (i.e. presence of Potts domain walls where ∆n = ±1) prevents order in the ϕ (red). Only below T3 does a single domain cover the whole system. As one increases J2/J6, a new sequence of phase transition occurs. Firstly, there is the development of composite qϕ = 1 vortices, bound states of three q ϑ = 1/3 defects, at a confinement transition labeled as T conf . Below this, there is a network of local Potts domains (c) due to the absence of free q ϑ = 1/3 defects and their associated "dangling" Potts domain walls. Within these domains, neither ϕ nor ϑ is ordered, and a composite nematic vortex can sit there. A zoom of such vortex is presented in (d). Beyond the red halo corresponding to the vortex core, the vortex is perceived as one with qϕ = 1 and q θ = 1. Within this core, however, is a structure where three q ϑ = 1/3 vortices are bound together by the Potts domain walls joining them. Below T3, the system becomes a single domain of the local Potts order (e), and there are free composite nematic vortices. Motivated by experiment, we confirm the presence of a Potts phase at a temperature T3 > TKT,nem associated with relative ordering of the hexatic and nematic degrees of freedom and the vanishing of the Potts domain walls.
imalist coupled hexatic-nematic model on a 2D square lattice [7][8][9][10], where the hexatic degrees of freedom (θ, invariant underθ i →θ i + 2πn/6 for integer n) describes the orientational order of neighboring molecules' center of mass, and the nematic degrees of freedom (φ, invariant underφ i →φ i + πn for integer n) corresponds to the orientation of the rod shaped molecules of 54COOBC with respect to a fixed laboratory axis (see Fig. 1  (1.1) vortices are present with "charge" q associated with the phase winding 2πq around them; here q ϑ = ∆ϑ 2π = 1 3 and q ϕ = ∆ϕ 2π = 1. When the coupling λ is finite, the two vortex types are no longer independent, since now ϑ − ϕ ≡ 2π 3 n (mod 2π) where n is an integer (the mod 2π simply associatesñ i = −1, −2 to n i = 2, 1, respectively). There are therefore three inequivalent relative alignments of ϕ and ϑ (see Fig. 1 (c)). This suggests the presence of a well-defined 3-state Potts order parameter, which we write as (1.3) Sinceσ i = 2π 3 n i for finite λ , then a finite mσ = |Mσ| is the direct consequence of long-range order in the relative Potts variable.
Another drastic consequence of a finite λ , the vortex charges are now related by the relation where ∆n is the number of walls encountered that separate different Potts domains, whereσ →σ+ 2π 3 across the Potts domain wall. This expression has two important consequences: • q ϑ defects are bound to Potts domain walls (in black in Fig. 2 (a)) which follows from Eq. 1.4 as 1 3 = 0 + 1 3 . • Integer q ϕ = 1 vortices are composites formed of q ϑ = 1 3 defects bound by domain walls (Fig. 2 (e)) since, referring to Eq. 1.4, 3 × 1 3 = 1.
When there are no free q ϑ = 1 3 charges and thus no "dangling" Potts domain walls, a unique local Potts order parameter can be defined. Therefore, the binding of q ϑ = 1 3 vortices at the hexatic KT transition for J 2 J 6 results in a network of Potts walls ( Fig. 2 (b)) separating distinct local Potts domains. At some lower temperature one expects one such domain to dominate the system leading to long-range Potts order, and a finite value of the order parameter mσ from Eq. 1.3. Indeed previous computational work in this parameter regime confirms the presence of a Potts ordered phase below the hexatic Kosterlitz-Thouless transition [9,10] as indicated in the schematic phase diagram in Figure 2.
The situation on the nematic KT side of the phase diagram (J 2 J 6 ) is more subtle. From previous studies [16][17][18][19][20][21][22][23][24][25], we expect that there exists a parameter regime where bound states of three q ϑ vortices form above the nematic KT transition (Fig. 2 (d)). This is driven by confinement of q ϑ = 1/3 fractional vortices into a composite and extended q ϕ = 1 vortex. At this "confinement" temperature, local Potts order develops. If this scale was merely a crossover, then there could not be any lower temperature transition into Potts long-range order. By continuity with the J 2 J 6 side, the low-temperature ordered phase must have Potts LRO. Therefore, one has to embrace the results that the confinement is a true phase transition [22]. For lower temperatures than the confinement transition, the associated Potts domain walls disappear and there is long-range Potts order (Fig. 2 (e)), again reflected in a finite value of the order parameter mσ. This will certainly be the case when the composite vortices bind at T KT,nem ; but the coincidence of these two transition temperatures would surely be indicative of an underlying unknown symmetry. The 54COOBC measurements suggest Potts ordering at temperatures above that of the nematic vortex binding [1][2][3]. We therefore probe whether we can tune Eq. 1.2 to a parameter regime of its phase diagram where there is a Potts phase above the nematic Kosterlitz-Thouless transition (Fig. 2 (c), on the right side); this would emulate the experimental observation.
The proposed Potts phase involves the relative orientation of the disordered hexatic and the nematic phases. It is then composite in nature since it will necessarily involve higher-order correlations of the primary hexatic and nematic angles. A simple example of such composite order is nematicity in localized spin models such as the two-dimensional (2D) J 1 -J 2 Heisenberg square lattice model [26][27][28]. There, nematic order corresponds to a relative ordering of spins, which breaks fourfold lattice rotation symmetry even though the spins retain their full SU(2) rotation symmetry. A path towards the development of composite order is via partial melting of an underlying multi-component primary order parameter due to increasing fluctuations. Being intertwined with a primary order, composite orders often occur in proximity to other ordered phases, giving a natural explanation for the complexity of phase diagrams observed in correlated systems. Examples are frustrated magnets [26][27][28][29][30][31][32][33][34][35][36][37][38], unconventional superconductors [39][40][41][42][43][44][45][46][47][48], ultracold atoms [49] and liquid crystals [50].
Above the defect-binding transitions, vortices are the "elementary excitations" of the hexatic-nematic coupled system of Eq. 1.2 with charge associated with their phase winding. It is known that, as is the case here, integer charge vortices can lower their energy by splitting into multiple fractional charge vortices linked by domain wall strings [18], if it is energetically favored. Such fractionalization describes the phenomenon where the elementary excitations of a system carry fractional charges (or quantum numbers) that are different from ones of the microscopic degrees of freedom [51]. Fractional excitations are often confined into objects with integer charge by a strong string tension force, yet can unbind at a confinement-deconfinement transition [22,[52][53][54]. Confinement is well-known from elementary particle physics, where it describes the binding of quarks, which carry fractional electric charge, into integer charge baryons or mesons. It also occurs frequently in condensed matter systems, for example, in low-dimensional magnets [55][56][57], spin ice models [33,58], quantum Hall systems [59], coupled atomic-molecular superfluids [18], and generalized XY models with vector magnetic or nematic degrees of freedom [20,22,60]. Here we propose that the 54COOBC films provide another setting for this phenomenon, arguing that confinement of fractional nematic vortices, as shown in Fig 2 (e), leads to the composite Potts phase at temperatures above that of the nematic Kosterlitz-Thouless transition. Confinement of fractional nematic vortices drives the binding of three elementary q ϑ = 1/3 hexatic vortices. This bound state formation removes the dangling ends of Potts domain walls and thus enables the appearance of a well-defined Potts order parameter in the system.
In this paper, we use large-scale parallel-tempering classical Monte-Carlo simulations to investigate the finite temperature phase diagram of the model described in equation 1.2, which captures the relevant degrees of freedom in the experimental system and is consistent with its symmetries. We demonstrate that it contains a small region where Potts order exists even though the underlying hexatic and nematic angles remain disordered. We show that the transition to the Z 3 ordered phase lies in the 2D Potts universality class, which is characterized by a specific heat divergence with scaling exponent α = 1/3, in good agreement with the experimental scaling exponent. Hexatic and nematic QLRO only develop at slightly lower temperatures via a Kosterlitz-Thouless (KT) phase transition. We support our unbiased numerical findings by analytical arguments that describe how fractionalization of nematic vortices can lead to extended vortex cores and separate the Potts transition from the KT transition. The sequence of upper Potts and lower KT transitions serve as a natural explanation of the experimentally observed melting process from Hex-B to Sm-A' to Sm-A. This resolves the long-standing puzzle of why the two-step melting theory of KTHNY, which accounts well for the experimental observations at the Hex-B to Sm-A' (and the lower Cry-B to Hex-B) transitions, fails to explain the observed features at the Sm-A' (HO) to Sm-A transition.
The remainder of the paper is organized as follows: in Sec. II, we explore the thermodynamical phase diagram of the coupled hexatic-nematic model that describes 54COOBC films for zero and finite λ coupling, using analytical estimates. We highlight the interplay of vortex confinement, Potts domain-walls and fractionalization, which plays a key role. In Sec. III, we describe the Monte Carlo algorithm, and introduce the observables that we use to identify the different thermodynamic phases. In Sec. IV, we present results of our Monte Carlo simulations. This includes our main result: the finite temperature phase diagram of the coupled hexaticnematic model. We discuss the behavior of the system in different regions of the phase diagram. We close the paper in Sec. V with concluding remarks and outlook for future research directions.

II. HEXATIC-NEMATIC XY MODEL
For technical reasons, in the rest of the paper we study a model in the same universality class as that of Eqs. 1.1 and 1.2. This is done through a transformation of the hexatic and nematic degrees of freedom of the minimal coupled model of Eq. 1.1. Rescaling the angles θ i = 6θ i , φ i = 2φ i (equivalently, θ i = 3ϑ i , φ i = ϕ i ) so they both cover the range θ i , φ i ∈ [0, 2π), yields the dimensionless expression (2.1) Here, we have introduced J = 1 2 (J 2 + J 6 ), λ = λ /J and ∆ ≡ J 2 /J such that 0 ≤ ∆ ≤ 2 covers all exchange coupling ratios J 2 /J 6 . A value of ∆ = 1 corresponds to the isotropic limit J 2 = J 6 . This description of our minimal model in terms of two O(2) variables is most useful for the Monte-Carlo study we provide in the next sections, hence the change of variables.
In the remainder of this section, we provide an intuitive and semi-analytical description of the expected phase diagram of this model. Starting from the uncoupled model is subsection II A, we first show that the hexatic-nematic coupling term is a relevant perturbation that tends to induce relative Z 3 Potts order at temperatures larger than the KT transition temperature. Furthermore, we point out in subsection II B the important role of vortex excitations in the system, especially as to how they differ from those in the uncoupled model. In subsection II C, we show why the fractionalization of the nematic vortices leads to extended vortex cores, and that is a necessary requirement for the Potts transition to occur above the KT transition. Finally we review some numerical results of related models in subsection II D.

A. Uncoupled model and relevance of coupling term
In the absence of a coupling term, λ = 0, hexatic and nematic degrees of freedom undergo separate KT transitions at temperatures [5,[61][62][63][64]  by a sudden jump of hexatic (p = 6) or nematic (p = 2) spin stiffness K p = ρ p /T from zero to the universal value K p (T p,KT ) = 2/π. The background color in Fig. 3 shows that the specific heat c exhibits a broad hump above the transition at about T = 1.1T KT [61]. In the isotropic limit, ∆ = 1, the two KT transitions occur at the same temperature: T (λ=0) KT = 0.89J. In order to determine the effect of a finite hexaticnematic coupling term (λ > 0) on the phase diagram, we calculate its renormalization group (RG) scaling dimension D λ . This can be done straightforwardly at ∆ = 1 and T KT , where free vortex excitations are absent, to yield Here, we have used that K R,p = 2/π is the exact value of the renormalized spin stiffness at the KT transition temperature T (λ=0) p,KT [61]. Note that K p is independent of p only for the rescaled Hamiltonian in Eq. (2.1), where the minimal phase winding of both hexatic and nematic vortices is equal to 2π, but the final result (D λ = 3/4) is identical if one uses Eq. (1.2). Further details of this derivation are presented in appendix A. A positive scaling dimension D λ > 0 indicates that λ is a relevant perturbation at T (λ=0) KT and will drive the system away from the uncoupled KT fixed point. This suggests that the system develops long-range Potts order above the KT transition, T Z3 > T (λ=0) KT . Certainly as λ grows towards longer length scales, hexatic and nematic angles are forced to arrange into one of the three parallel configurations, shown in Fig. 1(d) for the equivalent ϑ and ϕ variables. The constraint due to λ, as expressed for Eq. (2.1), is then Here, n i = −1, 0, 1 is a Z 3 degree of freedom. Note that if the locking constraint is fulfilled, the model (2.1) becomes equivalent to a generalized XY model [60,65], in our case for q = 3 [66][67][68]. We will review the results of the fully locked model in subsection II D. The locking transition occurs as a crossover at a temperature T λ (λ) that depends on λ. Since D λ > 0, it follows that T λ > T KT for all nonzero values of λ. For small initial values of λ 1, T λ can be estimated from analyzing the RG equation of an XY model in a threefold potential, which yields T λ ≈ 8π 9 J [11]. For large initial values of λ 1, the locking occurs at a temperature T λ J. For T < T λ , the symmetry of the model is lowered

B. Vortex excitations in the coupled model
For T < T λ the relative angle between nematic and hexatic degrees of freedom is a Z 3 variable n i . Being discrete, n i can develop LR order at a finite transition temperature T 3 . The central question is whether T 3 lies above or below the KT transition of the "center-of-mass" O(2) variable. If T 3 > T KT the situation corresponds to the experimentally observed order of phase transitions in 54COOBC. This turns out to be a rather delicate issue [9,10] that requires a careful and unbiased large-scale computational effort, which is described below in Secs. III and IV. We note that several studies of related coupled XY models, obtained by taking the limit of λ → ∞ in Eq. (2.1), revealed intriguing behavior close to T KT at ∆ = 1 [67][68][69][70][71][72]. In this section, we give analytical arguments that reveal the subtleties which arise when addressing this question.
Vortex excitations in the system are extremely important, as their binding (unbinding) is related to the KT ordered (disordered) phases. Let us now discuss their role, and in particular the impact of nematic-hexatic phase locking at T < T λ as described by Eq. (2.4) on their formation. For simplicity, we focus on the isotropic point ∆ = 1, but similar arguments can be given for other values of ∆. It will be advantageous to refer to Hamiltonian 1.2, which lends itself to a clean interpretation of the vortex defects. We provide a lexicon for the vortex excitations for the Hamiltonian of Eq. 2.1, which is used extensively for technical reasons in the following sections. The same arguments can be made for the model of Eq. 1.1 in the original variables, for which the minimal winding around a nematic (hexatic) vortex is given by π (2π/6), hence their respective names. For simplicity, we use λ whenever we refer to the hexatic-nematic coupling term, irrespective of the particular model.
In the absence of the coupling term (λ = 0), hexatic and nematic systems undergo independent KT transitions (see Fig. 3). In the language of Hamiltonian 1.2, the minimal charge of a hexatic (nematic) vortex is q ϑ = 1/3 (q ϕ = 1) corresponding to a phase winding of 2πq ϑ and 2πq ϕ around a vortex, respectively. Thus, at the KT transition, point vortices of charge q ϕ = ±1 (q ϕ = ±1/3) unbind in the nematic (hexatic) system. Free vortices limit the correlation lengths ξ 6 (ξ 2 ) for ϑ (ϕ) to a finite value for all T > T KT . Vortices are bound into pairs of total charge zero in the critical phase for T < T KT , where the correlation lengths are infinite.
The situation is notably different at nonzero λ and temperatures T < T λ below the locking crossover. While the domain of the angles reads ϑ ∈ [0, 2π 3 ) and ϕ ∈ [0, 2π), the locking condition imposes a distinct constraint on the phase winding around nematic and hexatic vortices, such that where ∆n counts the changes in the Potts index n as one loops around a vortex core. Note that the modulo operation originates from the definition of the Potts variables n, such that n = −2 ≡ 1.
If the system is Potts disordered, one expects a sample to be swarmed by a network of domain walls of the σ = 2π 3 n variable, as domains get smaller and smaller for temperatures above a Potts ordering temperature T 3 . For such a system with global Potts disorder, ∆n = 0 generically. Hence, Eq. 2.5 leads to 1 3 ≡ 0+ 1 3 , i.e. domain walls where ∆n = 1 are necessarily attached to hexatic vortices. These domain walls are energetically costly due to the significant nematic gradient energy that arises across them. Proliferation of these vortices and their eventual unbinding leads to an hexatic KT transition above the T 3 transition. This is the mechanism at play at small J 2 /J 6 (or, alternatively, for small ∆ in Eq. 2.1), as it can be see in Fig. 2 (c), on the left side.
In the presence of Potts order, one has that globally, ∆n = 0. The solution to Eq. 2.5 is then 1 ≡ 1 or 3( 1 3 ) ≡ 1. Both cases correspond to nematic vortices, but there is a subtle difference between the two. In the first case, an hexatic vortex of charge three times its elementary charge is at the same site as a nematic vortex. This is a nematic point vortex.
The alternative is for the nematic vortex to correspond to a triad of hexatic vortices. Since each hexatic vortex is attached to a Potts domain wall, one simple way to solve such a triad is to link all three hexatic vortices via their domain wall, at it can be seen in Fig. 2 (e). In essence, the domain walls act as a binding force for the fractional vortices, i.e. the hexatic vortices, as they imply locally fractional nematic vortices. Such a split vortex can be much more extended than its point vortex counterpart. In subsection II C, we provide an analytical argument that extended nematic vortices are lower energy than their point vortex. Even though the domain walls can be quite energetic, vortices of charge triple that of their elementary charge are extremely costly due to the phase winding around the vortex core.
Note that if one uses the convention of Eq. (2.1), where both θ and φ are within the interval [0, 2π), the minimal phase winding around a vortex is given by 2π in both cases, corresponding to charge q θ = 1 and q φ = 1 elementary vortices. The relative locking between hexatic and nematic angles is described then by Eq. (2.4), leading to the constraint that Thus, a nematic vortex with minimal 2π phase winding must be accompanied by a hexatic vortex of charge q θ = 3 that exhibits a phase winding of 2πq θ = 6π. In both descriptions, a nematic vortex of minimal charge q φ necessarily pairs with a hexatic vortex whose charge is three times larger than its minimal charge. The hexatic vortex is then a situation where q θ = 1 and ∆n = 1, which means that there is a Potts domain wall due to the 2π/3 mismatch in the φ variable. Similarly, the composite nematic vortex has q φ = 1 and three q θ = 1 vortices, such that (1 + 1 + 1)/3 ≡ 1, with each hexatic vortex confined via the attractive domain wall. For all the following sections, we use the model presented at Eq. 2.1. Translation between the different vortex formalism (Eq. 1.2 vs 2.1) can be done through this subsection. In other words, in this model, the transition to the relative Potts phase corresponds to confinement transition of fractionalized q ϕ = 1/3 nematic vortices, which drives the formation of a bound state of the three attached elementary q ϑ = 1/3 hexatic vortices through the constraint of Eq. 2.5. In the absence of free elementary hexatic vortices, a Potts order parameter can be well defined even in the presence of nematic vortices. In contrast, no independent Z 3 degree of freedom exists below the nematic KT transition, as the potential λ acts as a uniaxial field for the hexatic degrees of freedom. Furthermore, elementary hexatic vortices are attached to "dangling" domain walls which, if they were not confined, would fully destroy the LR Potts order that is set in the system. The composite vortices of total elementary nematic charge are then the only ones allowed in the relative Z 3 ordered state. We cover the energetics of the extended nematic vortices in the following subsection. This transition is expected to lie in the 2D Potts universality class, which is characterized by the critical exponents α = 1/3, β = 1/9, γ = 13/9, ν = 5/6 [73]. Note in particular that the specific heat exponent α, such that c ∝ t −α , leads to a pronounced divergence of the specific heat distinct from the smooth behavior observed for KT transitions with a broad hump at about 1.1T KT [61]. We will exploit this notable difference to distinguish between Potts and KT transitions in our Monte-Carlo simulations in section IV.

C. Vortex fractionalization and the extended vortices
Extended nematic q ϕ vortices are the free "elementary excitations" in the system at temperatures below T λ . Lower charge vortices are held together by domain wall strings with finite tension, arising from the gradient energy cost of the domain wall due to the "fractional" hexatic vortex they contain. This observation immediately questions the possibility that a KT transition can take place below the Potts transition, because point-like vortices are expected to bind at a temperature above a possible Potts transition T 3 . Specifically, setting ∇ϑ = ∇ϕ once below T λ (alternatively, setting ∇θ = 9∇φ due to Eq. (2.4)) leads to T λ =0 KT,nem = 10 T λ=0 KT ≈ 10π 2 J. The estimated KT transition temperature is thus even larger than T λ ≈ 8π 9 J (which bounds T 3 in the small λ regime). This estimate, however, leaves out the possibility of vortex fractionalization and the emergence of extended vortices [16,17]. Indeed, it is well known that higher charge vortices can reduce their (gradient) energy by splitting into multiple lower charge objects. In addition, hexatic and nematic degrees of freedom are continuous, therefore the domain walls separating regions with different Potts variable n i can acquire a finite width ξ dw . The width ξ dw is determined by a balance between the cost of violating the locking constraint in Eq. (2.4) imposed by the λ term and the gain in gradient energy by distributing the angle mismatch over a finite length. The width of the domain walls approaches the minimal size of the lattice constant only as the renormalized λ J. For example, splitting a q θ = 3 vortex into three q θ = 1 vortices reduces the gradient energy by a factor of 3 2 − 3 = 6 [18] (similarly, splitting a q ϑ = 1 = 3( 1 3 ) vortex). Previous mean-field studies of continuous models with the same symmetry have shown that such an arrangement is energetically beneficial [20,[23][24][25]74].
Here, this splitting, however, implies fractionalization of the joint nematic vortex into three vortices of fractional charge q φ = 1/3. These are held together by domain wall strings where ∆n = ±1. The competition of gradient and domain wall energy results in an extended vortex core [18]. By comparing the energy of a point vortex (3 = 3(1)) to a split vortex (1 + 1 + 1 = 3(1), which looks like a point vortex at distances larger than the core), we arrive at the conclusion that it is energeti-cally favored to split the vortex. For a circular geometry as in Fig. 2 (e), this leads to an optimal radius of R 2.5a for the value of the coupling λ used in this work, or a split vortex with an area 6.25 times greater than its point vortex counterpart. Thorough derivation of this result is presented in Appendix B. Importantly, the unbinding of extended vortices is known to occur at a lower temperature than the unbinding of point-like vortices [19]. This can be understood from the fact that the initial value of the vortex fugacity increases by a factor of (R/a) 2 , where a is the microscopic lattic scale. For R a, the KT transition temperature scales as T KT ∼ 1/ ln[(R/a) 2 ] and is thus significantly reduced in the case of large vortex core sizes. The crucial open question is whether it is reduced to a value below the Potts transition temperature. To address this question in an unbiased way, we have performed large scale classical Monte-Carlo simulations, which will be discussed next.

D. Previous numerical studies
The particular model presented in Eq. 2.1 was first studied in two dimensions using Monte-Carlo techniques [9,10], yet with a focus on a different region of the phase diagram (∆ 1) and at substantially smaller system sizes than presented here. A three-dimensional version of the model was studied as well [75,76]. A related model, with O(2)×Z 2 symmetry instead of our O(2)×Z 3 symmetry, associated with atom-molecular mixing, was studied extensively in two dimensions [77,78]. The ∆ 1 regime was however not the focus of those studies.
The community has been more focused on the infinite coupled limit. Such generalized XY models can be written as where δφ ij = φ i − φ j , with p = 3 corresponding to the λ → ∞ limit of Eqs. 1.2 and 2.1, and p = 2 is obtained through the same limit for the atomic-molecular mixing model. These two limiting model were studied via renormalization group studies [21,22,60,65,79,80], Monte-Carlo techniques [67-70, 81, 82], matrix-product states [72], bosonization [71,83]. They all share a common phase diagram structure [84], with a well understood regime at ∆ 1 with a p-state discrete transition at T d (Ising or Potts) below a Kosterlitz-Thouless transition at T KT where fractional defects unbind. For ∆ 1, the available evidence supports a single transition temperature corresponding to both confinement/deconfinement transition of fractional vortices and an unbinding of integer charge vortices [21,22,71]. From our analysis of the composite vortex in section II C, we find that there is no gain at λ → ∞ for an extended vortex over a point one. This would likely collapse the two transitions we expect into one in that limit, such that they would be unseparable with current numerical accuracy. This is partly why our numerical study is focused on the coupled model in the intermediate regime, so as to see the postulated effect of the formation of composite vortices.
We note that in the case of the fully frustrated XY model (FFXY), which represents a periodic array of Josephson junctions with half a quantum flux through it, one rather has two O(2) order parameters coupled via a term ∼ cos [2(θ i − φ i )]. For this model, extensive analytical [80,[85][86][87][88] and Monte-Carlo simulations [62,[89][90][91][92][93][94][95] were able to distinguish a clear regime where the relative Z 2 symmetry breaking happens at a higher temperature than the global KT transition. It was also shown that the presence of the Ising line above the KT transition leads to the prospect of an emergent supersymmetry at the Ising-KT multicritical point [96]. This success was due to the power of modern implementation of Monte-Carlo algorithms, as well as the ever-increasing computing power available, further motivating our reexamination of this classic hexatic-nematic problem.

III. MONTE-CARLO SIMULATION: ALGORITHM AND OBSERVABLES
In this section, we describe the algorithm used for our large-scale parallel tempering Monte-Carlo (MC) simulations of the model in Eq. (2.1). We give details in subsection III A, and introduce the various observables that are measured to obtain the phase diagram of the model, presented in Fig. 4, in subsection III B. The detailed investigation of the full phase diagram is presented in Sec. IV. All data and code needed to generate the results in this paper are available online [97].

A. Technical details of Monte-Carlo algorithm
The employed MC algorithm combines a standard parallel tempering update [98,99] with a generalized Wolff cluster algorithm adapted to coupled XY models [77,100]. The simulations were performed on a square lattice with periodic boundary conditions of N = L × L sites, with L ranging from 8 to 420. For system sizes L ≤ 200 (L > 200), we simulate 40 (64) configurations in parallel at different, geometrically spaced temperatures between T min = 0.5J and T max = 1.6J to obtain a general overview of the phase diagram. When looking directly at the nature of a phase transition at T c , we always used a new temperature range such that (T min , T max ) = (0.95 T c , 1.05 T c ). The generalized Wolff algorithm takes into account that due to the coupling term, λ i cos(θ i − 3φ i ), one must flip the hexatic and nematic angles, θ i and φ i , in an anisotropic manner in order to explore both global O(2) and Z 3 symmetries. The Wolff clusters are therefore built as follows: one first starts by randomly selecting a flip direction η ∈ [0, 6π), a site i and the type of variable, θ i or φ i . Depending on which variable was chosen, we then apply the flip rule We perform the simulations for at least 4 × 10 5 MC steps, where each step consists of a parallel-tempering and a generalized Wolff move. To ensure thermalization, we discard the first half of obtained configurations and measure thermodynamic observables only during the second half of MC steps. We introduce the different observables that we measure next. Finally, error bars are obtained using the standard jackknife method [101].

B. Observables, phases and phase transitions
To distinguish the different phases of the model, we investigate several thermodynamic observables that can be grouped into three classes. First, we measure observables related to the energy fluctuations in the system, the specific heat per site c and the Binder cumulant of the energy B E : The second class of observables we study are magnetizations and their susceptibilities that characterize the different phases Here, we have introduced the Z 3 Potts variable as the following: such that the hexatic-nematic coupling Hamiltonian term reads λ cos[3σ i ]. We also measure the associated susceptibilities and Binder cumulants (a = θ, φ, σ): Below we perform a scaling analysis of m a , χ a and c in order to extract the critical exponents of the observed second order phase transition. The Binder cumulants of the magnetizations undergo step-like transitions from a value of 1/3 at high temperature to a value of 2/3 at low temperature at both KT and second-order phase transitions.
Finally, the third class of observables that we investigate are spin stiffnesses ρ α , which describe the free energy change of the system under a uniform phase twist (along a given lattice direction α). When hexatic and nematic variables are coupled to each other by nonzero λ, one cannot separate contributions arising from each individual variable and thus only extract the total stiff-ness of the system. In order to fulfill the potential term constraint (2.4), we apply a uniform phase twist of the form φ i+α,i , θ i+α,i → φ i+α,i + ψ, θ i+α,i + 3ψ , (3.6) where φ ji = φ j − φ i , θ ji = θ j − θ i , Here, ψ is uniform across the system andα corresponds to eitherx orŷ lattice direction. Note that the phase twist that is applied to the hexatic angle θ i across each bond is three times larger than the twist applied to the nematic angle. The spin stiffness for a twist along direction α is defined as ρ (α) = 1 N ∂ 2 F ∂ψ 2 and following a standard derivation [106], one arrives at the explicit expressions Here, we have defined the total stiffness ρ that is averaged over both lattice directions. The summation i, j α runs over all bonds along direction α. For λ = 0, cross correlations between hexatic and nematic variables are absent and the total stiffness decomposes into the sum ρ(λ = 0) = ρ 2 ∆ + 9ρ 6 (2 − ∆), where ρ 2 (ρ 6 ) is the stiffness of the nematic (hexatic) system for a twist angle of ψ. Note that the factor of 9 in the ρ 6 part is due to the fact that we apply a uniform twist angle of 3ψ to the hexatic angles. It is well known that the discontinuous jump of the spin stiffness at the KT transition can be directly associated with the charge q of the unbinding vortices (phase winding of 2πq around the vortex) [68,107,108]. In particular, one finds just below the transition, while the stiffness vanishes above the transition. The hexatic-nematic model of Eq. (2.1) supports two types of vortex excitations, associated with the hexatic and nematic angles. Because of the choice we made in Eq. (3.6) to apply the phase twist uniformly on the nematic φ variables, then the type of vortex unbinding in φ dictates the value of the stiffness jump. A KT transition dominated by the unbinding of q φ = ±1 vortices results in a normal jump of the total stiffness ρ by 2T KT,2 /π. We thus denote a KT transition characterized by this jump value as "nematic KT", T KT,2 ≡ T 2 . On the other hand, a KT transition driven by the unbinding of integer charge q θ = ±1 hexatic vortices, accompanied by their corresponding Potts domain walls, results in a jump of the total stiffness by 18T KT,6 /π. A KT transition characterized by this larger jump value is denoted as "hexatic KT", T KT,6 ≡ T 6 . To summarize, the location of stiffness jumps is used to determine the KT transition temperature T KT , and the height of the jump provides clear evidence on the type of vortex unbinding that occurs across the transition. In the following, we will use these observables to identify the different phases and phase transition universality classes presented in Fig. 4.

IV. COMPUTATIONAL RESULTS
Here we present results of extensive MC simulations of the model (2.1). Our main result is the phase diagram shown in Fig. 4 as a function of temperature T and coupling constant ratio ∆ = J 2 /J with J = 1 2 (J 2 + J 6 ). The phase diagram is obtained for fixed hexatic-nematic coupling strength λ = 2.1. Crucially, it shows a region around to ∆ = 1, where the Potts phase transition T 3 lies above the nematic KT transition T KT, nem , which is highlighted in Fig. 4(b). Although the two transitions occur close to each other and the ratio of transition temperatures reads T 3 /T KT (∆ = 1) = 1.007, a detailed numerical analysis presented below shows that T 3 > T KT with a high degree of statistical confidence. This numerically demonstrates that LR order in the emergent Potts Z 3 variable n i , as defined related to σ i in Eq. 3.4, exists even in a region with finite-range correlations of the hexatic and nematic degrees of freedom. This suggests that the HO phase that is experimentally observed in 54COOBC films is characterized by LR Potts order and the sharp specific heat divergence is associated with a 2D Potts phase transition. In the following subsections, we separately discuss the different regions in the phase diagram in the order of increasing values of ∆.
A. Upper KT and lower Potts transition at ∆ < 0.9 As shown in the phase diagram in Fig. 4(a), for ∆ < 0.9 and starting from the disordered high-temperature phase, the system first develops hexatic QLRO across a KT transition at temperature T 6 . This transition corresponds to an unbinding of hexatic vortices with charge q θ = ±1. Since T 6 < T λ , i.e. the temperature at which the λ coupling becomes relevant, each hexatic vortex is found at the end of a domain wall string, and can also be thought of as a fractional nematic vortex of q φ = 1/3. As T 6 > T 3 , Potts domain walls have proliferated already at lower temperatures. As shown in Fig. 5(a) the unbinding of fractional nematic vortices results in a jump of the total stiffness equal to ρ(T − 6 ) = 18 π T 6 . We determine T 6 by extracting T 6 (L) for different system sizes L using this criterion and then extrapolating to infinite system size. As shown in Fig. 5(b), the specific heat c exhibits only a broad maximum at about 1.1 T 6 , as expected from a KT transition. In contrast, c develops a pronounced peak at lower temperature T 3 , which grows with system size L. (a) Total spin stiffness ρ defined in Eq. (3.7) as a function of temperature T /J. Upturn at T6 signals KT transition into the hexatic phase. The transition temperature is determined most accurately using the universal jump criterion ρ(T − 6 ) = 18T6/π, which yields T6(∆ = 0.5) = 1.35±0.01. (b) Specific heat per site c as a function of temperature T /J for ∆ = 0.5. While c exhibits a broad maximum above the hexatic KT transition at 1.1 T6, it develops a sharp peak at the lower transition at T3, which increases with system size L. This is characteristic of Potts phase transition in 2D, where c ∼ t −α with exponent α = 1/3. We extract the transition temperature T3(∆ = 0.5) = 0.693 ± 0.005 from the location of the maximum of c for L = 60. Inset depicts the three Z3 domains with different relative hexatic-nematic ordering in the hexatic phase. Long-range Z3 order only develops at the lower T3 transition.
Such a power law singularity is expected, for example, at a 2D Potts phase transition, where c ∼ t −α with α = 1/3 and t = ( The inset of Fig. 5(b) schematically depicts the ordering that occurs at the lower transition: in the hexatic phase there exist different domains of the three-state Potts order parameter n i ∈ Z 3 (see Eq. (2.4)). Regions with different n i are separated by domain walls that also represent domain wall strings for the nematic angle, where φ winds by 2π/3. There exist fractional vortices with q φ = 1/3 at the end of these strings. In this regime, the cost of these domain walls is set by the nematic gradient energy, which is proportional to J 2 , which explains why T 3 → 0 as ∆ → 0. The transition at T 3 corresponds to the ordering of the Z 3 variable n i , which occurs via a 2D Potts phase transition. Since n i is a discrete degree of freedom, the system exhibits true LR Potts order below T 3 . This scenario of an upper hexatic KT transition, where θ variables develop QLRO, and a lower Potts transition, where n i (φ) develop LRO (QLRO, respectively) can also be observed in the Binder cumulants of the respective magnetizations, as shown in the left column of Fig. 10 in Appendix C. It is also clearly seen in the Binder cumulant of the energy B E , shown in Fig. 8. At the lower transition at T 3 , the Binder cumulant B E develops a sharp peak, which is a clear indication of a second-order phase transition. In contrast, close to the upper transition T 6 , B E only features a change in slope at about 1.1 T 6 , which is known to correspond to its behavior close to a KT transition. The Binder cumulant B E is, therefore, a convenient way to distinguish a second-order phase transition from a KT transition, as will be discussed more below. Finally, we note that such an order of phase transitions (upper KT, lower Potts) has previously been reported for the generalized XY model [9,68,77] that corresponds to the λ → ∞ limit of Eq. (2.1), as it is mentioned in section II D.

B. Upper Potts and lower KT transition at
0.9 < ∆ < 1.15 Let us now discuss the region of main interest in the phase diagram around ∆ = 1.0. As illustrated by the Binder cumulants B θ , B φ , B σ in the middle column of Fig. 10 in Appendix C, the ordering of hexatic, nematic and Potts degrees of freedom all occur at nearby temperatures in this region. As shown in Fig. 4, we find that the Potts (KT) transition temperature increases (decreases) monotonically with increasing ∆ until the three Binder cumulants cross at a value of ∆ = 0.9. For 0.9 < ∆ < 1.15, the Potts transition occurs above the KT transition. As the KT transition temperature exhibits a minimum at ∆ = 0.9, both transitions track each other with similar slope in that region. While the separation of the two transitions is small, (T 3 − T KT )/T 3 ≈ 1%, we can resolve them within error bars using extensive Monte Carlo simulations up to system sizes of L = 380. This is demonstrated in Fig. 6, which shows results for ∆ = 1.0. Figure 6(a) shows the spin stiffness ρ and the specific heat c that are used to extract KT and Potts transition temperatures, respectively. Specifically, we extract T 3 (L) from a fit of c(T, L) ∝ t −α with t = [T − T 3 (L)/T 3 (L)] and Potts exponent α = 1/3. As shown in detail in Fig. 6(b), we determine the KT transition temperature T 2 from the universal jump criterion ρ(T, L) = 2T 2 (L)/π. We note that the nature of the KT transition changes from hexatic to nematic at ∆ = 0.9. Figure. 6(c) displays the resulting system size dependent transition temperatures T 2 (L) and T 3 (L). We determine the transition temperatures in the thermodynamic limit by extrapolating to infinite system sizes using the ex- pected scaling forms T 2 (L) − T 2 (∞) = a/(ln L) 2 and T 3 (L) − T 3 (∞) = a L −1/ν with Potts correlation length exponent ν = 5/6. We find that T 3 (∞) > T 2 (∞) within a confidence of more than one standard deviation, specifically T 3 (∞) = 1.2022±0.0005 and T 2 (∞) = 1.194±0.002 at ∆ = 1.0. This is one of the main results of this work. In the thin region T 3 > T > T 2 the system exhibits LR order in the relative hexatic-nematic orientation n i even though both angles are still only short-range correlated with a finite correlation length. We believe that this scenario of inverted Potts and KT transitions can explain the experimental observations in thin films of 54COOB, as discussed above.
Once the system enters a region with LR Potts order, the only available asymptotically free vortex excitations are nematic q φ = 1 vortices, formed of three bound elementary q θ = 1 vortices, confined by the Potts domain wall's finite tension. As discussed in Sec. II C, the energy cost of these nematic vortices can be significantly lowered at finite values of λ by expanding the vortex core to host three hexatic vortices of q θ = 1. This decreases the transition temperature for the unbinding of such extended combined vortices by a factor of 1/ ln[R/a] 2 , where R > a is the vortex core size and a is lattice scale. As a result, the KT transition temperature for unbinding of nematic vortices may be pushed below T 3 , as we observe in our MC simulations. We emphasize that while we have focused on the value of λ = 2.1 in our work, we indeed find a region with T 3 > T 2 also for other values of λ = 0.5, 1.0, 4.0, demonstrating that our conclusions hold for an extended regime of hexatic-nematic couplings. We suggest that future work should explore the regime of small λ in more detail, since the splitting of T 3 and T 2 is potentially larger there, and one could observe them directly.
Still, we note that we cannot fully exclude a scenario with a single transition in the Potts universality class, but we consider this unlikely at finite λ based on our MC results. In other words, our numerical analysis shows that assuming that the stiffness makes a universal jump corresponding to an unbinding of nematic q φ = ±1 vortices, the associated KT transition would lie slightly below the Potts transition. In contrast, a single transition is expected in the λ → ∞ limit of our model (2.1), in analogy to the phase diagram obtained for the vector-nematic generalized XY model version of our model [21,22,69], where it is argued that vortex unbinding is suppressed by critical fluctuations of the discrete order parameter. Further work, both numerical and analytical, should be directed towards the λ → ∞ limit of Eq. (2.7), especially in the region of ∆ 1.0, to explore its phase diagram further and test this conjecture, especially in the region of ∆ 1.0.
To further study the nature of the phase transitions, we perform a finite size scaling analysis of specific heat c, Potts magnetization m σ , and susceptibility χ σ . As shown in Fig. 7, a proper rescaling of the axes leads to a data collapse of results for different system sizes 40 ≤ L ≤ 200 and temperatures. The collapse is most complete for the magnetization m σ , but also evident for FIG. 7. Finite size scaling analysis of (a) specific heat c, (b) Potts magnetization mσ, and Potts susceptibility χσ for system sizes between L = 40 and L = 380. We use the indicated scaling functions, which correspond to a 2D continuous Potts transition and include corrections to scaling. The resulting critical exponents and the Potts transition temperature T3 are collected in Table I. Data collapse is best for mσ, but also works fairly well for c and χσ.
c and χ σ . The critical exponents that we extract from the scaling analysis are consistent with a transition in the Potts universality class. The scaling analysis yields a transition temperature of T 3 = 1.20 ± 0.01 that agrees with the more precise value that is obtained from the scaling of the maximum of the specific heat with L. This confirms that the upper transition at T 3 lies in the Potts universality class, corresponding to the breaking of   Fig. 7. We use the scaling formÕ(t, , L]/(1 + a1L −ω 1 ), which includes corrections to scaling. The 2D Potts critical exponents are α = 1/3, β = 1/9 ≈ 0.11, γ = 13/9 ≈ 1.44, and ν = 5/6 ≈ 0.83, which is in fair agreement with our findings with the exception of ν. Note that when tracking only T3, and including small system sizes, we obtain a value of ν that lies much closer to its expected value (see Fig. 6 as well).
FIG. 8. Binder cumulant of the energy BE = E 4 E 2 2 − 1 as a function of temperature T /J for different values of ∆. The system size was set to L = 80. At ∆ = 0.5, the energy Binder develops a sharp peak at the Potts transition T3 (vertical dashed line) as obtained in Fig. 4(b). The apparent change in slope at around T = 1.5J coincides with the characteristic broad hump in the specific heat at around 1.1 T6 = 1.5 with T6 = 1.35J. At ∆ = 1.0, the sharp peak associated with the Potts transition has moved up in temperature and is located at T3 that we extract from the specific heat (see Fig. 6). The presence of a sharp peak, which persists up to the largest system sizes we consider, L = 300, is a clear indication of a Potts transition at ∆ = 1.0. Note that since BE ∝ L −2 we have rescaled the L = 300 curve by (300/80) 2 to lie on top of the L = 80 curve. In contrast, such a sharp peak is notably absent at ∆ = 1.5, where one only observes a change in slope around T = 1.6J associated with the nematic KT transition T2.
symmetry associated with the development of LR order in the relative orientation of hexatic and nematic angles. Finally, we briefly discuss the phase diagram at larger values of 1.15 ≤ ∆ ≤ 2, where only a single nematic KT transition occurs as a function of T (see Fig. 4). We identify the location and nature of the phase transition by the universal jump criterion of the stiffness, ρ(T 2 ) = 2T 2 /π, and find that the jump corresponds to an unbinding of nematic vortices of charge q φ = 1. Due to the locking constraint (2.4), these are equivalent to q θ = 3 hexatic vortices. Again, by fractionalizing the nematic vortices inside an extended vortex core so that there are three q θ = 1 vortices, the energy of such a composite vortex can be lowered, especially at smaller values of λ. As T 2 < T λ , QLRO of the nematic degrees of freedom immediately leads to QLRO of the hexatic variables, as they feel a uniaxial potential with a single minimum. Both O(2) and Z 3 symmetries are thus broken at the nematic KT transition for ∆ > 1. 15. This is confirmed clearly by the Binder cumulants B θ , B φ , B σ at ∆ = 1.5, which are shown in the right column of Fig. 10 in appendix C. Both B φ and B σ undergo a sharp transition at T 2 , while T 6 exhibits a smoother crossover towards the value of 2/3 consistent with QLRO as temperature is lowered. The Binder cumulant of the energy B E only shows a change in slope above T 2 and does not exhibit any additional feature at temperatures T < T 2 .
Finally, the specific heat exhibits two (broad) maxima, one at a temperature above T 2 , which corresponds to the usual broad KT maximum at 1.1 T 2 , where nematic vortices proliferate. In addition, there appears a maximum below T 2 , that roughly follows the hexatic phase transition in the uncoupled system at λ = 0 (see Fig. 3). This crossover is due to the production of q θ = 1 hexatic vortices that are bound to fractionalized and confined q φ = 1/3 vortices. These hexatic vortices are thus bound together by domain wall strings of finite tension and only release a small part of the entropy, leading to a less pronounced maximum of c.

V. CONCLUSION AND OUTLOOK
Motivated by an long-standing experimental mystery in liquid crystal films of 54COOBC, we have presented a detailed Monte Carlo simulation study of a generalized XY model with finite coupling strength λ between the two angular degrees of freedom. Our key finding is that there exists a parameter region in the phase diagram (for ∆ ≈ 1 in Eq. 2.1), where a sharp specific heat signature occurs at a temperature (T 3 ) above that where the spin stiffness jumps to its finite value; all the thermodynamic exponents are consistent with Z 3 symmetry. This phase has relative nematic-hexatic order, and can be characterized as one of free q φ = 1 nematic vortices formed of three bound q θ = 1 hexatic vortices. At lower temperatures, our numerical caling analysis demonstrates the presence of the Kosterlitz-Thouless transition (T KT,2 or T KT, nem ), where composite nematic-hexatic vortices become bound. For the coupling value (λ = 2.1) in our simulations we find the ratio of temperature scales to be T 3 /T KT,2 = 1.007. Both nematic and hexatic degrees of freedom only develop quasi-long range order at T KT,2 . Our simulations show no signs of first-order behavior anywhere in the phase diagram.
This scenario closely resembles the experimental observations in single-layer films of 54COOBC with a sharp specific heat signature between two smectic phases (Sm-A to Sm-A') at T Sm-A' = 66 • C and a lower KT transition into a hexatic phase (Sm-A' to Hex-B) at T Hex-B = 63 • C. The ratio of experimental transition temperatures is given by T Sm-A' /T Hex-B = 1.018, consistent with our findings. The transition into the Sm-A' phase is characterized by a sharp specific heat divergence with α ≈ 0.3, in agreement with the 2D Potts exponent α = 1/3 that we obtained in our investigation of this minimal model. While experimentally the hexatic correlation length makes a sudden jump into the Sm-A' phase, it remains finite and hexatic QLR order only develops at T Hex-B , again resembling our simulation results of a lower KT transition.
There are several additional features that emerge from our study. For example we reproduce previously reported results at smaller values of ∆, where a Potts transition [9] occurs below a hexatic KT transition, T 3 < T KT,6 (= T KT,hex ), that is characterized by the unbinding q θ = 1 vortices. Each hexatic defect is attached to a Potts domain wall (see Figure 2 (a)), so that the binding of these defects is accompanied by the formation of a network of local Potts domains. The prior establishment of a Potts transition at lower temperatures can be understood as the vanishing of these Potts walls and the formation of a single Potts domain.
Our analytic estimates suggest that confinement fractionalization play important roles in the emergence of the composite Potts phase. At a temperature T conf > T KT,2 , q φ = 1 vortices form composed of bound-states of three q θ = 1 defects (see Figure 2 (d)). Again the absence of "dangling" Potts walls means that a network of local Potts domains is formed. Then the question is whether the Potts ordering T 3 is above or at T KT,2 since there can be no Potts domain walls once the q φ defects bind. The extended vortex core sizes lower the KT transition to a value below T 3 , revealing a new phase. Though our numerical studies yield sharp thermodynamic features at a temperature T 3 , it is not possible to distinguish a crossover from a transition here since our system sizes are typically less than the average separation between free defects.
In order to maintain continuity and consistency with previous computational studies, we explored the hexaticnematic XY model (Eq. 2.1) with the same coupling value (λ = 2.1) as before [9]. Since the vortex core size increases with decreasing hexatic-nematic coupling strength λ, we expect that the two temperature scales T 3 and T KT,2 will be further separated at smaller values of λ, which should be explored in the future. Significantly larger system sizes should be accessible using techniques such as the worm algorithm [21,22,109] and matrix product states [72].
An alternative tuning parameter is hexatic ringexchange coupling [110][111][112][113][114], that is known to tune vortex core energies and could thus be used to suppress the lower KT transition, expanding the Potts phase. Tuning this parameter can be done numerically, and it has been shown experimentally that core energies and sizes are correlated with liquid crystal densities [114]. Adjusting these parameters may be done by varying local density, possibly by bending the freely-suspended liquid crystal films of 54COOBC. Because the proposed Potts phase involves relative hexatic-nematic ordering and hence higher-order correlations, it is challenging to probe experimentally since most probes measure twopoint correlation functions.
Taking our cue from intriguing results in a related model [21,22,71], it is an open question whether the two temperature scales T 3 and T KT,2 converge into one single transition at large values of λ. The two end points in ∆ of our composite Potts relative phase are worth further study. At ∆ = 0.9, the possible existence of a multicritical fixed point of O(2) × Z 3 symmetry suggests the exciting possibility of an emergent supersymmetry [96], in analogy to the FFXY [62]. Close to ∆ = 1.15, where the Potts transition disappears, there may be a deconfined critical point [21,22,71].
In conclusion we have presented the 54COOBC liquid crystal film problem as an experimental realization of the deconfinement of fractional vortices; here the "mystery" phase above the nematic KT transition is one of composite Potts ordering in the relative hexatic-nematic angles. We have reproduced the observed specific heat behavior as a function of temperature in our computational study of a minimalist model. The emergent Potts phase hosts extended nematic vortices that are analogous to baryons. Furthering the analogy, the three hexatic vortices that form such a bound state are akin to quarks. From this perspective, the experimentally observed Potts transition in the liquid crystal samples is a laboratory realization of quark confinement and baryon formation. Our work raises several further questions regarding the deconfinement of fractional vortices, in particular the need for a fully developed analytic theory for the interaction of domain walls with fractional vortices. We hope that our results will stimulate further interest in these problems.
We explore the relevance of the hexatic-nematic coupling ∝ λ cos [6(θ −φ)] in Eq. 1.1. The intermediate numerical factors are different for the models of Eqs. 1.2 and 2.1, but the end result is identical. We first turn to a long wavelength version of the model For K 2 = K 6 (i.e. ∆ = 1.0 in Eq. 2.1), one has that the hexatic T KT,6 and nematic T KT,2 are at identical temperatures, i.e. T KT,6 = T KT,2 . We can find whether the coupling constant λ 3 =λ/T is a relevant perturbation at the KT transition temperature by using the fact that at the transition T KT,p the exact value of the renormalized stiffnesses is known to be with bare K p = J p /T , from the model presented in Eq. 1.1. Following Refs. [61,115], we evaluate the correlation function associated with the coupling between ϑ and ϕ, assuming no correlations betweenθ andφ, i.e., [θ(x) −φ(0)] 2 = 0, valid for λ 3 = 0, via calculating where we used 1 2 [ψ(x) − ψ(0)] 2 ≈ 1 2πK ln x a0 for Gaussian field ψ with coupling K. The exponent is Via the Kadanoff construction [61,115], we find whether λ 3 is a relevant perturbation by determining the sign of the scaling dimension Using Eq. (A2) at T KT,6 = T KT,2 , we obtain This means that λ 3 is a relevant perturbation at T = T KT,6 = T KT,2 and the system will flow away from the KT transition at non-zero λ 3 , and hence the system will encounter a high-temperature T λ scale where the variables become locked with one another, satisfying the λ coupling. This completes the derivation of the result presented at equation 2.3, obtained for model 2.1 but otherwise identical.
Appendix B: Estimation of the size of composite vortices.
In this section, we provide a rough estimate of the size of a composite vortex of total charge q φ = 1 formed of three hexatic vortices (q θ = 1) bound through a domain wall of the σ variable, such that σ → σ+2π/3 across such wall. This procedure was done in the following references [20,74,116], with relevance to quark deconfinement, and in [18] for the Z 2 symmetric case. We largely reproduce the treatment of this last reference here for our Z 3 symmetric case.
We consider an effective long wavelength model analogous to equation 2.1, obtained by transforming cos (θ i − θ j ) ∼ dr(∇θ) 2 /2 into the continuum limit, with K 2 = ∆/T , K 6 = (2 − ∆)/T and h = λ/T : Note that one can relate K 2 and K 6 to those of the original model of Eq. 1.1, used in appendix A, such that K p = p 2 K p . We consider two types of point defects in the two XY variables It is clear that the two variables have their own independent behavior for h = 0. However, for any finite h, an extensive energy cost is incurred if the system does not lock the two phases on average with each other such that θ = 3 φ . This further leads to the winding numbers to be related to each other, as it is explained in Eq. 2.6. The first way to satisfy this is to simply have q φ = 1 and q θ = 3. However, as we will show later, this is very costly as the phase has to wind very tightly around the vortex,  2) variables. (c) Visualization of a hexatic vortex q θ = 1, with its adjacent domain wall of the relative variable σ = 2πn/3 (as defined in Eq. 3.4), of width ξ dw . This can also be viewed as a fractional q φ = 1 3 vortex. (d) A composite vortex of total charge q φ = 1 and q θ = 3, which is clearly seen outside the core. It is formed of three fractional ones bound together by their domain wall strings. These are arranged in this simple geometry where the walls are of length R and the distance between q θ = 1 vortices is √ 3R.
generating a large core energy for the vortex. The second way to satisfy this constraint is to have q θ = 1 and ∆n = 1 (or equivalently, q φ = 1/3) In the simplest picture, this fractional vortex in the φ generates a line defect originating at the vortex core (see Fig. 9), through which φ will rapidly wind by 2π/3. It is possible however that some part of the rapidly winding "leaks" into the θ variable. For h extremely large or K i very small, we expect that this wall defect will be very thin (ξ dw ξ i ), but as h is decreased, it should become wider and wider, leaving to φ more space to wind the remaining. Note that in both case, the energy of such a wall will scale linearly with the length of the wall, such that E dw ∝ L.
− K∇ 2 σ + 3h sin (3σ) = 0 We solve this equation for the situation where we have σ(x, y → ∞) → 0 and σ(x, y → −∞) → 2π/3 (this is done without loss of generality, as all domain walls have a jump of 2π/3), which leads to a domain wall along the x direction. The solution for this domain wall is then simply σ dw (y) = 4 3 arctan (e 2y/ξ dw ) with ξ dw = 2 3 2K 3h (B12) We get that α dw (y) = K ab Ka σ dw (y) and then, for our initial variables One has that the energy of a domain wall of length L is E dw = dw L with There is a caveat to this result. As K → 0 then the energy of the wall goes to 0 and so does its width. On the other hand, as h → ∞, the wall width goes to 0 but the energy increases dramatically. It is important in this context to bring back the lattice cutoff a ∼ 1 here. Hence, we have ξ dw = 2 3 2K 3h > 1. For h < 4K/9 the theory is well behaved. For larger h one has to use ξ dw = 1 and adjust Eq. B19. Note that putting back the initial parameters K 2 , K 6 and h in the domain wall energy we get K = 9(K 2 K 6 )/(K 2 + 9K 6 ) which leads to the following expression for the domain wall energy: dw = 4 6h K 2 K 6 K 2 + 9K 6 (B20) or, using ∆, With this new insight, we wish to compare the energy of a point defect with q φ = 1 to the assembly of three hexatic objects with each q θ = 1, linked by a domain wall, for the same total charge. This arrangement is shown in Fig. 9. We denote the first situation as a point vortex (pv) and the second as a split vortex (sv). The energy of a pure defect with no domain wall, i.e. one in which a q θ = 3 vortex is accompanied by a q φ = 1 vortex, is simply expressed as: E pv = E c,θ (3) + E c,φ (1) + 9πK 6 ln (L/ξ 6 ) + πK 2 ln (L/ξ 2 ) (B22) = π 2 2 (9K 6 + K 2 ) + 9πK 6 ln (R/ξ 6 ) + πK 2 ln (R/ξ 2 ) + 9πK 6 ln (L/R) + πK 2 ln (L/R) = E c pv (R) + 9πK 6 ln (L/R) + πK 2 ln (L/R) (B24) the second line uses the first approximation to the core energy of the vortices E c ∼ π 2 q 2 K/2, with q being the charge of the vortex. Note that in this formula, we have explicity put is an arbitrary radius R, and split the log parts. This will come in handy when we compare this new core energy E c pv (R) to the split core energy E c sv (R), as we consider R to be the radial length of a split combination of vortices. We will then be able to compare only the energetics inside the radius R, since for r > R, both combinations will look and act like a point vortex. This point vortex is thought to be the dominant kind as h K, since then the domain wall width becomes extremely small and its energy very large.
In the case of a split vortex, we need to include some extra energetic terms, i.e. the logarithmic repulsion between vortices themselves V q−q (r) = −2πKq 2 ln (r/ξ). There are three such terms here for each type of vortex, and they are separated by a distance of √ 3R for this simple geometry (see Fig. 9). We also include the attractive domain wall energy dw R for each of the domain walls. We then get E c sv (R) = 3(E c,θ (1) + E c,φ (1/3)) − 3 × 2πK 6 ln ( At this point, we stress that point configurations, the point vortex and the split vortex, are topologically equivalent at r > R. We can then compare the two configurations in order to determine the region in which a split vortex would be less energetic than the point vortex. We start by finding out the optimal size of a split vortex by finding R max such that ∂E c sv (R) ∂R | Rmax = 0. In the following, we used ξ 6 = ξ 2 = a. This leads to the following expression for the optimal radial size of a split vortex R max /a = 2π K 2 /9 + K 6 √ 3 dw (B27) One sees that driving the coupling h to infinity completely annihilates the split vortices structures (by dw → 0). There is however a regime with low h where split vortex will be quite extended. Note that this expression does not depend on the temperature, but only on the interaction parameters themselves.
In the case of strong coupling, i.e. h > 4K/9 (for K 2 = K 6 = 1, this is h > 4/10) then the domain wall width is ξ dw = a = 1 i.e. the lattice spacing. This is the case for the λ = 2.1 we have chosen for our simulations. In this case, we get σ dw (y) = 4 3 arctan (e 2y ) (B28) which then leads to dw = 16K/9 (B29) and then, for ∆ = 1.0, our computational value, we have R max /a 2.5 (B30) which completes the derivation of the vortex size estimate that is presented in the text, in section II C. Using this value for the radius of the composite vortex into the energy estimates from equations B24 and B26, we get that for such an extended vortex, proving that in the regime of our Z 3 relative order phase, composite vortices are favored with respect to point ones due to their extended nature, even with the cost of domain walls in the core of the vortex. sharply increase close to T = 1.2J, suggesting the presence of a Potts transition. The crossing B σ (L) yields a transition temperature T 3 1.202 that is consistent with our findings in Fig. 4. Note that the Binder cumulants do not allow us to easily address the question whether the nematic and Potts transitions are separated. To demonstrate this we rather rely on a refined scaling analysis, presented in Fig. 6 in the main text.