Quantum criticality in many-body parafermion chains

We construct local generalizations of 3-state Potts models with exotic critical points. We analytically show that these are described by non-diagonal modular invariant partition functions of products of $Z_3$ parafermion or $u(1)_6$ conformal field theories (CFTs). These correspond either to non-trivial permutation invariants or block diagonal invariants, that one can understand in terms of anyon condensation. In terms of lattice parafermion operators, the constructed models correspond to parafermion chains with many-body terms. Our construction is based on how the partition function of a CFT depends on symmetry sectors and boundary conditions. This enables to write the partition function corresponding to one modular invariant as a linear combination of another over different sectors and boundary conditions, which translates to a general recipe how to write down a microscopic model, tuned to criticality. We show that the scheme can also be extended to construct critical generalizations of $k$-state clock type models.


Introduction
At a quantum critical point two distinct phases of matter coexist. A remarkable feature of 1D systems is that such special points in the phase diagram are in general described by a field theory with conformal symmetry -a conformal field theory (CFT) [1,2]. In other words, the system exhibits a universal behavior regardless of the underlying microscopic model, i.e. what are the local degrees of freedom and how they interact. This universal description at the critical point enables to determine what phases of matter appear in the vicinity of the critical point when the system is perturbed away from it. However, the relation of the universality class to criticality of a microscopic model is a one-way problem. If a given microscopic model exhibits a critical point with diverging correlation length, given a catalog of possible CFTs it is comparatively easy to make ansatzes to verify which CFT describes it. The converse is not true though. While symmetries present in the CFTs constrain the possible microscopic models, there is in general no recipe to write down a model that exhibits a critical point described by a given CFT. The lack of framework to write down a microscopic model for a given CFT is an outstanding problem from several perspectives. From an academic perspective, it is of interest to understand the minimal microscopic conditions that can give rise to a given universal behavior. This enables to experimentally search for such behavior in existing or synthesized materials. While CFT is a standard theoretical tool, experiments probing CFT predictions beyond measuring critical exponents are still few [3,4]. Due to the coexistence of phases of matter at the critical points, critical models with universal behavior that are perturbed away from the critical point can also serve as starting points for models of gapped, possibly topological phases of matter. Recently, the interplay of symmetry-protected topological order and quantum criticality has been explored in spin-1/2 chains [5,6,7], with a unified picture emerging how the protecting symmetries and the possible fractionalized edge states determine the universality classes of transitions [8,9].
While the universal description of a critical point is a property of 1D systems (the nature of conformal critical points in 2D is an open question and a subject of cutting-edge numerical studies [10,11]), it has been proposed that critical 1D systems can serve as building blocks of 2D topologically ordered states of matter [12,13,14,15]. Critical 1D systems can be arranged into 2D array and when coupled together in a designed fashion, the system can become a gapped topologically ordered state whose nature depends only on the couplings and on the CFT describing the critical 1D systems. Thus new microscopic models with exotic critical points enable also the construction of new models for exotic 2D topological states of matter. In particular, since the advent of topological quantum computation, there is much interest to construct states that harbour non-Abelian anyons that could be employed for topologically protected quantum information processing [16].
Motivated by these open questions, in this work we advance the program initiated in our earlier works, where we constructed exactly solvable spin-1/2 chains for all criticalities in the so(N ) 1 universality class [17,18,6]. Instead of spin chains, we focus here on 3-state Potts models, where the local degrees of freedom are not spins, but clock variables, and construct generalizations that exhibit critical points whose critical behavior has not been previously discussed in the literature. 3-state Potts models themselves have recently attracted attention due to their relation to parafermion modes [19], that have been proposed to be realized, following the same principles that lead to the recent experimental discovery of Majorana modes [20,21], by inducing superconductivity via the proximity effect on fractional quantum Hall edge states, such as the ν = 1/3 Laughlin state [22,23]. A uniform array of such parafermion modes is unitarily equivalent to a 3-state Potts model via the parafermion version of the Jordan-Wigner transformation, the Fradkin-Kadanoff transformation [24]. The recent focus on parafermions arises not from themselves though, but from their collective behavior. Were they to hybridize in a 2D array, they could realize a state that hosts the coveted Fibonacci anyons that are universal for quantum computation [14]. This prospect, while speculative, has motivated research into 1D collective states of parafermion modes. Different parafermion phases have been understood from the symmetry protection perspective [25], the microscopics of Z 3 parafermion CFT describing the criticality of the 3-state Potts model has been analyzed [26] and the phase diagram in the presence of longer range parafermion tunneling has been studied [27].
While our primary motivation is the construction of new Potts-like models with exotic criticalities, our models turn out to be unitarily equivalent to parafermion chains with manybody interactions between the parafermion modes. Thus as a by-product of our construction we address the nature of the critical behavior of parafermion chains in the presence of such many-body terms that may also arise when parafermion modes hybridize. To our understanding, the effect of such many-body terms has not been considered previously in literature. We find that when many-body terms are comparable to local chemical potential-like terms, the chains are critical and described by the non-diagonal modular invariant partition functions of a CFT that is product of some number of Z 3 parafermion or u(1) 6 CFTs.
Usually, when talking about a critical point being described by a given CFT, one refers to the diagonal invariant. Our work explicitly shows that this is too restrictive. We demonstrate that a particular class of non-diagonal modular invariant partition functions -the permutation invariants -that are usually overlooked when talking about physical systems, also have local microscopic models that realize them. Before proceeding to the actual constructions and the physics underlying them, we summarize our main results -the microscopic generalizations of the 3-state Potts models and the exotic critical points they exhibit.

Summary of results: Generalized critical 3-state Potts models
The starting point for our construction is the critical nearest neighbour 3-state Potts model that is described by the Hamiltonian (see for instance Ref. [28]) The local clock operators X i and Z j commute on different sites, while on the same site they obey For future reference, we also define Y = ZX. Like a spin-1/2 chain can be written in terms of fermionic operators via a Jordan-Wigner transformation, 3-state Potts models can be written in terms of parafermion operators via a Fradkin-Kadanoff transformation [24,14]. Introducing the lattice parafermion operators that satisfy α 3 j = 1, α † j = α 2 j and α i α j = ω sign(i−j) α j α j , the Hamiltonian (1) takes the form While this model is quadratic in terms of the parafermion operators, it is strongly interacting and not solvable with a Fourier transformation. Due to the formal similarity to free, quadratic fermions, we refer to it describing 'quadratic' parafermions. 1 When higher order terms terms appear, such as four or six parafermion terms in the models we construct, we refer to the parafermions exhibiting many-body interactions. The nature of the criticality of the 3-state Potts model depends on the overall sign. As shown in Table 1, H − is in the universality class of Z 3 parafermion CFT with central charge c = 4/5 [29,30], while the criticality of H + is described by the u(1) 6 CFT with c = 1 [31,32,33]. The table also summarizes the new critical models we analytically construct and the CFTs that occur in them. Since they all involve terms such as that involve four or more parafermion operators, all the new critical models correspond to parafermions with many-body interactions.
To explain how these models were constructed, the paper is structured as follows. In Section 2.1 we introduce the key concept of modular invariant partition functions of CFTs that characterize the possible distinct critical behaviors. To illustrate these concepts and pave the  (2) is critical and described by a modular invariant CFT that has central charge c and the listed number of primary fields. In the CFT column Z 3 , u(1) 6 and su(2) 3 = c(Z 3 × u 1 (6)) denote the diagonal invariants corresponding to these CFTs, while π(G) and c(G) denote the non-trivial permutation and condensation invariants of the CFT G, respectively, and G ×n denotes a product of n CFTs of type G. Apart from the model with the su(2) 3 critical point that has a two-site unit cell, all the models are translationally invariant.
way for constructing the generalized 3-state Potts models, in Section 2.2 we revisit the relation between two critical transverse field Ising chains and the XY chain from the perspective of modular invariant partition functions. This example, as well as its generalizations, has been studied in our previous works [17,18] from the anyon condensation perspective, whose connection to the current approach of modular invariants is explained in Section 2.3. In Section 3 we review the Z 3 parafermion and u(1) 6 criticalities that appear in the 3-state Potts model and then in Section 4 we provide the detailed derivations of the new models listed in Table 1. Section 5 discusses the possible generalization of our construction to kstate Potts models and we conclude with Section 6. A detailed discussion on the fusion rule symmetries that are useful in constructing non-diagonal permutation invariants can be found in Appendix A.

The Method: Modular invariant partition functions
Our construction is based on the properties of the partition function of the CFT describing a critical point. The partition function of every two-dimensional CFT on the torus must be invariant under symmetry operations known as modular transformations [2,28]. In the current setting of one-dimensional chains, this means that the spectrum of a chain with periodic boundary conditions is described (in the thermodynamic limit) by a modular invariant partition function. It is possible that a given CFT admits several distinct modular invariant partition functions, in which case each describes critical behavior of distinct type (see [34,35,36,37] for some early references). Often in literature, due to these appearing most often in physically relevant systems, one implicitly refers to the diagonal partition function when referring to a particular critical point being described by a given CFT. However, in general this is too restrictive. There can also be other non-diagonal modular invariant partition functions, to which we refer to as permutation or condensation invariants for reasons to be explained below. We show that also these can in fact be realized in microscopic critical models. Before we do so, we should mention that many important properties, do not depend on the precise modular invariant, or even the boundary conditions. Perhaps the most important such property is the specific heat, which is determined by the central charge of the CFT. Nevertheless, many properties do depend on the details of (at least the low-lying part of) the spectrum, such as finite temperature and dynamical properties.
To show how different modular invariants can be realized, we start from a known critical microscopic model described by the diagonal invariant of a CFT and find all the non-diagonal modular invariant partition functions by employing additional symmetries of the CFT. Then we express these partition functions as linear combinations of the diagonal partition function when it is restricted to different symmetry sectors and have different boundary conditions. This linear combination is subsequently translated into a Hamiltonian term that, when added to the initial critical microscopic system described by the diagonal invariant, changes the nature of the criticality to the non-diagonal invariant. While the resulting system is in general non-local, in all cases we are able to find canonical duality transformations that give a local and translational invariant representative of the new critical model.
To demonstrate this method in detail, in this section we first review the modular invariance of CFT partition functions, and comment on the relation with extended algebras and the orbifold construction. Then, by revisiting our earlier work [17,18], we illustrate the method by showing how the critical XY model can be derived from two decoupled critical transverse field Ising chains. We also review the connection to anyon condensation that provides a simple criterium to predict when a given CFT admits particular type of non-diagonal invariant known as a condensation invariant.

Modular invariant partition functions
The partition function of a system described by the Hamiltonian H is given by Z = Tr(e −βH ). When the system is critical, the Hamiltonian can be expanded in the Virasoro algebra generators L i satisfying the Virasoro algebra with central charge c [2]. A special role is played by the operator L 0 . Its eigenvalues take the form ǫ = h φ i + n, where n a non-negative integer, that gives the energies of each state in the spectrum. All the energies are shifted away from integer values by h φ i , the scaling dimensions of a each chiral primary field φ i . The distinct primary fields φ i , whose number is finite in all the cases we consider, are specific to the CFT describing the critical point. Akin to usual symmetry sectors, every state that descends from φ i belongs to the same 'conformal tower' that can be viewed as a sector of the CFT. The partition functions of these sectors are the chiral characters associated with different primary fields φ i . Formally they are defined as the polynomials where the trace is over all the eigenstates of L 0 belonging to the conformal tower associated with the primary field φ i . Here q = e 2πiτ is a formal variable of the modular parameter τ whose powers are the energies ǫ of the states and the non-negative integers a n encode their degeneracy.
The primary field operators are in general chiral, which means that in 1+1D they can either propagate left or right. This means that a single chiral CFT can describe chiral states, such as edge states of 2+1D topologically ordered states, but not genuine 1D states that arise from local Hamiltonians (chiral operators do not have local representations). Instead, 1+1D critical systems are described by combining the two chiral halves of the CFT. The full partition function describing such systems can be written as Z = (qq) −c/24 Tr q L 0qL0 , where the bar denotes Hermitian conjugation. In terms of the left and right chiral characters χ φ andχ φ , the partition function of a 1+1D critical system takes the general form where φ 1 and φ 2 are summed over all primary fields of the CFT and M φ 1 ,φ 2 are non-negative integers.
Possible reparametrizations of the torus set stringent constraints on the allowed partition functions of a given CFT described by the matrices M with elements M φ 1 ,φ 2 . These modular transformations formally transform the modular parameter τ as Modular invariance of the partition function is then equivalent to demanding that where S and T are n × n modular matrices specific to a given CFT with n primary fields. Formally, they are obtained by studying the behavior of the chiral characters (5) under the modular transformations (7), but for most CFTs they can be found in the literature. The matrix T is always diagonal with entries T φ 1 ,φ 2 = e 2πi(h φ 1 −c/24) δ φ 1 ,φ 2 , while S takes a form that does not have a simple expression in terms of only the scaling dimensions h φ i . To find all matrices M subject to the constraints (8) amounts to classifying all modular invariant partition functions and hence all distinct critical behaviors that can occur in 1+1D systems. This, however, is a challenging task. A complete classification has been achieved only in a limited set of cases, such as the minimal models [37], but a general classification for arbitrary CFTs is still lacking. For a CFT with a small number of primary fields, one can perform a numerical brute-force search, which is the approach we take here. The different possible modular invariant partition functions fall into different classes. First, clearly the identity matrix M φ 1 ,φ 2 = δ φ 1 ,φ 2 commutes with both S and T . This is known as the diagonal modular invariant, which is the most common partition function associated with a given CFT. In the literature, a 1+1D critical systems being described by a CFT usually refers to this invariant. The second case corresponds to the matrix M being block-diagonal. Using the original characters χ φ i , one constructs new ones, χ ′ j = χ φ j 1 + χ φ j 2 + · · · , which are combined 'in a diagonal way' to form a new partition function Z ′ = j χ ′ jχ ′ j . Typically, not all the original characters appear in the block diagonal partition function, and it can happen that certain characters appear several times. For the purpose of this paper, we call these block-diagonal invariants condensation modular invariants due to their connection to anyon condensation explained in Section 2.3. One can view the new chiral characters χ ′ j as corresponding to a primary field of a CFT with the same central charge as the original one, but with a different field content with generically less primary fields. Thus criticality described by a condensation invariant of some CFT is always equivalent to the diagonal invariant of some other CFT. We note that there is direct relation between condensation invariants and CFTs with an extended chiral algebra. If the new vacuum character χ ′ 0 = χ 0 + j χ j , the primaries related to the χ j (which in the current paper, are chiral products of Virasoro characters) have integer scaling dimensions, that is, they are currents. Indeed, in the CFT corresponding to the condensation invariant, the chiral algebra has been extended with these currents. For more information about CFTs with extended algebras, we refer to [38].
Given the diagonal invariant, it is sometimes possible to construct a permutation modular invariant. It is of the form M φ 1 ,φ 2 = δ φ 1 ,π(φ 2 ) , where π(φ) denotes a permutation of the primary fields, which leaves the fusion rules of the primary fields invariant (we review fusion rule symmetries in Appendix A). In other words, all chiral characters of both chiral halves appear, but M is no longer the identity matrix, but a permutation matrix (thus, in this case, the chiral algebra is not extended). Depending on the symmetries of the CFT, it may or may not give rise to the same partition function as the diagonal one when (6) is written out as a polynomial in q andq. If not, then while the primary field content in each chiral half is the same as for the diagonal invariant, the local physical observables and the energy spectrum are different due to them being constructed from the left and right moving components of different primary fields. While criticality corresponding to permutation invariants has been little studied in the context of physical models, in this work we show that they can indeed also arise in local systems. We note that to construct permutation invariants, one does not have to start from the diagonal invariant. One can also obtain different invariants by performing a permutation of the fields on a block-diagonal (or condensation) invariant.

Example: Modular invariants for products of Ising CFTs
To illustrate these concepts, we give an explicit example of a case where distinct modular invariants appear -a product theory of two Ising CFTs. This example also illustrates how to deal with product CFTs and how the partition functions depend on twisted boundary conditions, which are the key methods for constructing new critical Potts-like models later. We refer to [39] for more information on the relation between the Ising ×2 and u(1) 4 CFTs that we will now discuss.
The Ising CFT has three primary fields, denoted by 1, σ and ψ, with scaling dimensions h 1 = 0, h σ = 1/16 and h ψ = 1/2. The S and T matrices are represented by For the Ising CFT, there is only a a single modular invariant partition function given by the diagonal invariant M = 1 3×3 . The product of two Ising CFTs, that we call Ising ×2 CFT, has c = 1 and nine primary fields that we will denote by pair of labels (1, 1); (1, σ); (1, ψ); . . . spanning all the possible combinations of the primary fields. The scaling dimensions for these product primary fields are additive in the constituent fields, i.e. h (φ 1 ,φ 2 ) = h φ 1 +h φ 1 . Similarly, their chiral characters are simply products χ (φ 1 ,φ 2 ) = χ φ 1 χ φ 2 and the representations of the modular matrices are given as the tensor products S Ising ×2 = S ⊗ S and T Ising ×2 = T ⊗ T . The matrix M is thus now a 9 × 9 matrix, but it is still straightforward to find all the modular invariant partition functions satisfying (8) by a brute force calculation. It turns out that there are three solutions given by The Z Ising ×2 is the diagonal invariant partition function of the Ising ×2 CFT, while Z π is a permutation invariant. In the case of Ising ×2 it is not independent, but gives rise to a partition function that is identical to the diagonal invariant. This is due to the 'layer' symmetry of Ising ×2 that leaves the theory invariant under relabeling of the primary fields On the other hand, the block diagonality and the absence of some primary fields in Z one finds that it should correspond to a diagonal invariant of a theory with three non-trivial primary fields with scaling dimensions hψ = 1 2 and h λ = hλ = 1 8 . This CFT, which is a c = 1 CFT of a compactified boson, is called u(1) 4 . In fact, this was to be expected since it is well known that the Ising ×2 CFT is obtained from the u(1) 4 CFT by means of an orbifold construction [40]. Going in the other direction, one can obtain the u(1) 4 CFT form the Ising ×2 theory by 'adding the bosonic field (ψ, ψ) to the chiral algebra'. We describe this relation between these CFTs from the perspective of anyon condensation below.
We have thus found that starting from the chiral Ising ×2 CFT, one can construct both the diagonal invariant and condensation invariant partition functions. Since it is well known that the critical transverse field Ising (TFI) chain is described by Ising CFT [28], the criticality corresponding to the diagonal invariant of Ising ×2 is clearly realized in a system of two decoupled TFI chains. The corresponding microscopic Hamiltonian is given, for instance, as the critical TFI chain with nest-nearest exchange interactions where σ α i are the usual spin-1/2 Pauli matrices. To write down a microscopic model for the condensation invariant, we begin by studying how the diagonal partition function of a single Ising CFT depends on the symmetry sectors and boundary conditions. The symmetry sectors are inherited from the TFI chain that has Z 2 spin flip symmetry and thus two symmetry sectors that we label by Q = 0, 1. A TFI chain can also have either periodic or anti-periodic boundary conditions, which we denote bỹ Q = 0, 1, respectively. Denoting the partition function in symmetry sector Q with boundary conditionsQ by ZQ Q , it is well known that [41] where there holds ZQ Q = Z Q Q as required by the duality between the symmetry sectors and boundary conditions when solving the transverse field Ising model with a Jordan-Wigner transformation [42]. Clearly, summing over the symmetry sectors with periodic boundary conditions gives the diagonal invariant partition function of the Ising CFT, Likewise, the Hamiltonian (14) for two decoupled TFI chains has Z 2 ×Z 2 symmetry and the boundary conditions can be independently chosen for both chains. Employing the property χ (φ 1 ,φ 2 ) = χ φ 1 χ φ 2 for product CFTs, one can then verify that the diagonal invariant of the Ising ×2 CFT is obtained by summing over the four symmetry sectors with periodic boundary conditions Instead of decoupled TFI chains, consider a coupled system where the boundary condition of one the Ising CFTs is given by the symmetry sector of the other, and vice versa. Employing the explicit forms of the partition functions (15), one can verify that one obtains now precisely the condensation invariant partition function (12) when summing over all symmetry sectors To implement this coupling in the decoupled Hamiltonian (14), we introduce the coupling boundary Hamiltonian where are the independent symmetry operators on even and odd sites, respectively. The Hamiltonian H Ising ×2 + H B correlates then the boundary conditions and symmetry sectors of the two chains precisely as done in the linear combination (17) for the condensation invariant partition function. By construction the resulting system must thus be critical and described by the u(1) 4 CFT. While this Hamiltonian is manifestly non-local, we have shown in earlier work [17,18] that by using non-local duality transformations it can be mapped precisely to the critical XY chain that is known to be described by the diagonal invariant of the u(1) 4 CFT. We mention that the term H B effectively causes the boundary conditions to be 'twisted'. In contrast to the usual, well-studied, twisted boundary conditions, in our case, the 'twist' in the boundary conditions is symmetry sector dependent.
The example we just described, allows us to comment on the relation between the construction and orbifolding. It is well known that the Ising ×2 CFT is the Z 2 orbifold of the u(1) 4 CFT. So, by constructing the condensation modular invariant, one is effectively doing an orbifold construction in reverse. On the level of microscopic hamiltonians, the 'condensation direction' is the more straightforward one, because one starts by a simple doubling. Realizing that H XY can be written as H Ising ×2 + H B after a non-local canonical transformation is in general much harder.
We have demonstrated how finding out all possible modular invariant partition functions enables to determine whether a given critical system enables to construct new microscopic models for the non-diagonal partition functions. Before applying this strategy to construct microscopic models for criticality in Potts-type models, we review briefly the connection to anyon condensation that provides an easy way to predict when such constructions might be possible without knowing all the modular invariant partition functions (which in general is a hard task for CFTs with increasing number of primary fields).

Anyon condensation perspective
There is an intriguing connection between the modular invariant partition functions and anyon condensation [43,44]. Chiral CFTs describe also the 1D edge states of 2D topologically ordered states of matter, with the possible anyonic quasi-particle excitations being in one-toone correspondence with the primary fields of the CFT. Via this bulk-edge correspondence the anyonic statistics of the quasi-particles are given by the scaling dimensions that are interpreted as their topological spins. If for some quasi-particle a the spin h a is an integer, then it may condense, which results in general in a transition to another topologically ordered state of matter. As this state has a different quasi-particle content in the bulk, it must also have a different CFT describing the edge and hence the framework of anyon condensation can in general be expected to also predict relations between CFTs themselves without referring to a particular realization of 2D topological order.
Indeed, in our previous works [17,18] we have shown that there is a precise counterpart. This insight enabled us to derive microscopic spin chains for all criticalities in the so(N ) 1 universality class. The example considered in Section 2.2 is the simplest case of this hierarchy (u(1) 4 ≃ so(2) 1 ) that also has a realization in generalized cluster models [6]. In particular, we showed that by correlating boundary conditions with symmetry sectors, those primary fields that are predicted to be confined in a condensation transition can be removed from the spectrum and that this changes the critical description precisely as predicted by the framework of anyon condensation [43,44]. In our example above, we employed this same correlating of symmetry sectors and boundary conditions to construct the modular invariant partition function (12). In the condensation language the block diagonality of Z u(1) 4 Ising ×2 resulting in the diagonal partition function of u(1) 4 CFT via the identifications (13) can be understood as follows. The boson (ψ, ψ) with h (ψ,ψ) = 1 condensed and thus became indistinguishable from the vacuum (1, 1). This results in the confinement of the primary fields (σ, 1), (σ, ψ), (1, σ) and (ψ, σ) that do not appear in the block-diagonal partition function. The consistency of the resulting theory requires that (ψ, 1) and (1, ψ) are treated as the same new fermionic field ψ and that the field (σ, σ) must split into two distinct fields λ andλ.
The advantage of the condensation picture is that even if not all modular invariant partition functions are known, the existence of a primary field with an integer scaling dimension that is also a simple current implies that a condensation is possible and hence a block-diagonal modular invariant partition function exists 2 . As it is easy to figure out which fields should be confined [43,44], this gives a principle how to construct condensation invariant partition functions by correlating symmetry sectors and boundary conditions and subsequently verify that modular invariance is indeed satisfied. We turn now to apply these ideas to construct generalized 3-state Potts models whose criticalities are described by non-diagonal modular invariants of either the permutation or the condensation type.

Building block: Critical behavior of the 3-state Potts model
In the example of Section 2.2, we started from decoupled TFI chains to construct the XY model corresponding to the u(1) 4 condensation invariant. To construct generalized critical 3-state Potts models, we follow the same strategy, but start from some number of decoupled 3-state Potts models. In this section we review the Z 3 parafermion and u(1) 6 criticality that appear in the 3-state Potts model (1) for different overall signs. We also summarize how the partition function depends on symmetry sectors and boundary conditions, which is needed for our construction.
For negative overall sign, the 3-state Potts model is described by the diagonal invariant of the Z 3 parafermion CFT [30], while for positive sign it is described by the diagonal invariant of the u(1) 6 CFT [33]. For both cases these are the only independent modular invariants. To describe the field content of these CFTs on equal footing, it is convenient to view the Z k parafermion CFT as the coset su(2) k /u(1) 2k , which holds for any positive integer k [29]. This enables to describe the primary field content of Z k theories in terms of the simpler field content of su(2) k and u(1) 2k CFTs. For generality, we provide these properties for generic k, although in this work we mainly focus on the case k = 3.
The su(2) k CFT has central charge c = 3k k+2 and contains k + 1 primary fields [2]. We denote these by ϕ l for l = 0, 1, . . . , k. The scaling dimension of each primary field is given by h l = l(l+2) 4(k+2) and they obey the fusion rules For example, one of the new models we construct has a su(2) 3 critical point, as shown in Table  1. This c = 9/5 CFT has four primary fields with scaling dimensions {0, 3/20, 2/5, 3/4}. On the other hand, the u(1) 2k CFT is the c = 1 compactified chiral boson theory with 2k primary fields [2]. The primary fields are labeled ϕ m , where m is an integer defined modulo 2k. Here we choose the convention that m lies in the range −k < m ≤ k, which means that the scaling dimensions are given by h m = m 2 2k and the fusion rules can be written as For example, for the u(1) 6 CFT describing the 3-state Potts model for a positive overall sign, there are six primary fields with scaling dimensions {0, 1/6, 1/6, 2/3, 2/3, 3/2}. The Z k parafermion CFT has the central charge c = 2(k−1) k+2 and it is equivalent to the coset su(2) k /u(1) 2k . Its primary fields are then labeled by the labels of the su(2) k and u(1) 2k theories, namely ϕ l m . In addition to the constraints on l and m we already introduced, due to the coset the labels also have to satisfy the additional constraint l + m = 0 mod 2 and the identification ϕ l m ≡ ϕ k−l m+k . Using this, one can choose labels φ l m such that |m| ≤ l, in which case the scaling dimensions of the k(k + 1)/2 primary fields are given by h l m = h l − h m . The fusion rules follow directly from the fusion rules of the su(2) k and u(1) 2k theories and are given by For the Z 3 parafermion CFT describing the 3-state Potts model for a negative overall sign, there are six primary fields with scaling dimensions {0, 2/3, 2/3, 2/5, 1/15, 1/15}. In the literature these are often also labeled as For completeness, we note that the modular T matrix is determined by the scaling dimensions as discussed in Section 2.1, and the S matrix is given by where ω = e 2πi/3 and φ = (1 + √ 5)/2. The ordering of the fields is (ϕ 0 0 , ϕ 0 2 , ϕ 0 −2 , ϕ 2 0 , ϕ 2 2 , ϕ 2 −2 ). While the fusion rules are not of paramount importance to our construction, we have have included them to highlight a possible ambiguity in using the labels l and m to label the primary fields ϕ l m . If a permutation of the labels leaves the fusion rules invariant, i.e. there is a fusion rule symmetry, then the partition function corresponding to the permuted labeling may or may not be identical to the diagonal one. This can be checked by writing out the partition functions (6) as polynomials in q andq and see whether two different labelings give rise to the same partition function or not. Thus the fusion rule symmetries provide easy ansatzes to try and see whether a given CFT admits non-diagonal partition functions. For instance, for the Z 3 parafermion CFT with primary fields ϕ l m , where l = 0, 2 and m = −2, 0, 2, the only permutation of the fields that leaves the fusion rules invariant is ϕ l m → ϕ l −m . However, because the characters associated with ϕ l m and ϕ l −m are identical, this permutation gives again the diagonal partition function. The fusion rule symmetries for Z 3 and u(1) 6 CFTs, and for their product theories relevant to the present work, are discussed in more detail in Appendix A.

Partition functions for different symmetry sectors and boundary conditions
The dependence of the partition functions on the symmetry sectors and twisted boundary conditions is known for the Z k parafermion and the u(1) 2k CFTs [41,47]. A k-state Potts model always has a Z k symmetry described by the operators P = L i=1 Z i , where Z k = 1 and the phase factor ω appearing in the commutation relations of the 3-state clock variables is replaced by ω = e 2πi k . This means that there are always k symmetry sectors, that we label by Q = 0, 1, . . . , k − 1, corresponding to the eigenvalues ω Q of P. Similarly, there are k possible twisted boundary conditions ψ(x + L) = ωQψ(x) that we label byQ = 0, 1, . . . , k − 1.
Using this notation, the partition function for symmetry sector Q with twisted boundary conditionsQ is given compactly by [47] where the factor 1 2 accounts for the double counting due to the identification ϕ k−l m+k ≡ ϕ l m . We note that the effect of twisting the boundary conditions only involves changing the labels m of the parafermionic fields ϕ l m . In the case of the u(1) 2k CFT, the partition functions for fixed symmetry sector and boundary condition take a completely analogous form. One simply drops the label l from the field to obtain ZQ Q = −k<m≤k χ ϕmχϕ m−2Q , where the sum is constraint by m −Q mod k = Q.
Analogous to the example presented in Section (2.2), to construct critical generalizations of the Potts models we start from several copies of either the 3-state Potts models at either Z 3 or u(1) 6 critical point, or mixtures of them. The criticality of such systems is then described by products of the respective CFTs. The primary field content of product theories is obtained simply by taking all possible products of the primary fields of each CFT in the product. The scaling dimensions of these products are given as the sum of the scaling dimensions of the constituent fields and the fusion rules are just direct products of their fusion rules. Since the symmetries sectors Q n and boundary conditionsQ n of each CFT in the product are independent, the partition function of a product theory of N CFTs in a given symmetry sector and for given boundary conditions is given simply as the product ZQ 1 ,...,Q N Q 1 ,...,Q N = ZQ 1 Q 1 · · · ZQ N Q N . For instance, a product of two Z 3 parafermion CFTs has 36 primary fields (a 1 , a 2 ) with scaling dimensions h (a 1 ,a 2 ) = h a 1 + h a 2 . There are nine symmetry sectors labeled by Q 1 , Q 2 = 0, 1, 2 with nine independent boundary conditionsQ 1 ,Q 2 = 0, 1, 2. For each choice of them the partition function is given by ZQ 1 Q 1 ZQ 2 Q 2 .

Generalized 3-Potts models for non-diagonal modular invariants
This section contains the detailed derivation of the critical generalizations of 3-state Potts models listed in Table 1. We consider products of up to three Z 3 parafermion and u(1) 6 CFTs that can be realized by the 3-state Potts model for different overall signs. To obtain the non-diagonal modular invariant partition functions, for products of two CFTs we determine all the distinct invariants by numerically solving the equations (8). For products of three CFTs, a brute force solution becomes challenging due to the size of the modular matrices. Instead, guided by the anyon condensation perspective (see Section 2.3) and the fusion rule symmetries (see Appendix A), we consider all possible boundary terms and determine which of these give rise to non-diagonal modular invariants. For each case, we write the invariant as a linear combination over symmetry sectors and boundary conditions (see Section 3.1) and, similar to the example discussed in Section 2.3, construct a Hamiltonian term that implements the change in the nature of criticality in the microscopic setting. Finally, we introduce canonical duality transformations to find unitarily equivalent local models for each non-diagonal modular invariant.

The permutation invariants of Z ×2
3 and u(1) ×2 6 CFTs We begin by considering the doubled Z ×2 3 parafermion CFT that is realized by the nextnearest neighbour Hamiltonian that describes two critical decoupled 3-state Potts chains. By brute force calculation we find that there are altogether 16 different modular invariants M . As expected from Z ×2

3
CFT containing no primary fields with integer scaling dimensions (so there are no bosons to condense from the anyon condensation perspective), all the invariants are permutation invariants. However, because of the fusion rule symmetries, there are only two distinct modular invariant partition functions (see Appendix A for details). All invariants related by the 8 permutations where s 1 = ±1 and s 2 = ±1, give rise to the same a partition function as the diagonal invariant of the Z ×2 3 CFT. On the other hand, the 8 permutations give rise to a permutation invariant corresponding to the partition function where again l + m = 0 mod 2.
Employing the knowledge how the partition functions depend on the symmetry sectors and boundary conditions (23), it is straightforward to verify that this permutation invariant can be expressed as the linear combination In other words, the boundary condition of the first chain is given by the symmetry sector of the second chain, and vice versa, but in a twisted manner. The Hamiltonian term implementing this coupling between the two 3-state Potts chains is given by where P n = L/2−1 j=0 Z 2j+n are the independent symmetry operators for each chain. While the Hamiltonian H Z ×2 3 + H B is manifestly non-local and breaks translation invariance, it can be brought into a translationally invariant form with the duality transformations where j = 1, 2, . . . , L/2. It is straightforward to verify that the dual operatorsX i andZ i satisfy again the clock algebra (2). In terms of them one obtains the Hamiltonian which is critical by construction and described by the non-trivial permutation invariant (27) of the doubled Z 3 parafermion CFT.
There is a clear difference in the energy spectrum of the criticality described by the diagonal or the permutation invariants Z ×2 3 . Recalling that the lowest energy states of each conformal  tower labelled by the scaling dimensions (h l , h r ) have energy h l + h r , Table 2 shows the degeneracies of the lowest lying states. For the diagonal invariant there are four conformal towers with scaling dimensions (h l , h r ) = ( 1 15 , 1 15 ) and four with (h l , h r ) = ( 16 15 , 16 15 ). These are absent for the permutation invariant and replaced by eight towers with (h l , h r ) = ( 1 15 , 16 15 ) and (h l , h r ) = ( 16 15 , 1 15 ), that correspond to non-diagonal combination of the chiral halves. This means that if one compares the spectra of these models, the lowest lying excitations that are present in H Z ×2 3 are absent in H π Z ×2

3
. We confirmed this behaviour by explicitly diagonalizing the Hamiltonian (31) for system sizes up to L = 12 sites, showing that the rescaled finite-size spectrum indeed clearly differs from the diagonal invariant and precisely corresponds to the non-trivial permutation invariant (27).
The permutation invariant of u(1) ×2 6 CFT The u(1) ×2 6 CFT is realized by two decoupled critical Potts chains described by the Hamiltonian −H Z ×2 3 , i.e. just by changing the overall sign of the Hamiltonian (24). Again, there are no fields with integer scaling dimensions, so there are no condensation invariants, but after accounting for the fusion rule symmetries we find again a single independent permutation invariant partition function Z π u(1) ×2

6
. While writing it explicitly out is not particularly illustrative, its form can be seen from Table 2, which shows how the energy spectrum of a critical system corresponding to it differs from one described by the diagonal invariant. As above, some of the lowest lying states of the diagonal invariant are missing for the permutation invariant, which exhibits also conformal towers corresponding to non-diagonal combinations of the chiral halves.
To construct a microscopic model that realizes the permutation invariant, we find that the same coupling (28) between the two chains is required. Since changing the overall sign in (24) does not affect the duality transformations, the construction in the previous subsection applies directly also here. Thus the permutation invariant partition function Z π u(1)

The su(2) 3 condensation invariant of the Z 3 × u(1) 6 CFT
The remaining case that involves only a product of two CFTs is the Z 3 × u(1) 6 criticality that is realized by two decoupled 3-state Potts models with opposite signs. This case is slightly more tricky than the previous two cases. When comparing the spectrum of a microscopic Hamiltonian to the predictions of a conformal field theory, one has to take a non-universal velocity into account, which is accomplished by an overall rescaling of the spectrum. If one couples two identical microscopic chains, the velocities are identical and one does not have to worry about them. However, in general the velocities of the two chains can be different. The Bethe ansatz solution of the 3-state Potts model indicates that the velocity of the Z 3 model is twice the velocity of the u(1) 6 model [32].
We can correct for this difference in velocity, without changing the universal description by Z 3 × u(1) 6 CFT, by hand and include a relative factor of two in the Hamiltonian, so that both criticalities have the same velocity. Our starting point thus is the Hamiltonian Unlike the two preceding cases, this CFT contains a primary field with an integer scaling dimension suggesting that a condensation invariant exists. Indeed, by numerically evaluating the modular invariants, we find only the condensation invariant Identifying the blocks as new primary fields with scaling dimensions derived from the constituent fields, we find that this condensation invariant corresponds to the diagonal invariant of the su(2) 3 CFT described in Section 3. This was to be expected, since the Z 3 parafermion CFT is equivalent to the coset su(2) 3 /u(1) 6 . Multiplying in an additional u(1) 6 factor essentially reverses the coset construction.
To obtain a local microscopic Hamiltonian whose critical point is described by the su(2) 3 CFT, we find again that one has to use the same coupling (28) and hence the same boundary term (29) as in the case of the two coupled Z 3 parafermion CFTs. Thus the same canonical transformations can be employed and we find the local Hamiltonian (34) which is translationally invariant with respect to a unit cell of two sites. Having constructed microscopic representatives for all different modular invariant partition functions that can be constructed from two critical 3-state Potts models, we now turn to consider non-diagonal modular invariants of three decoupled 3-state Potts models described by either tripled Z ×3 3 parafermion CFT or u(1) ×3 6 CFTs. Any mixed product of three CFTs yields only invariants that are products of those obtained above for doubled theories with an additional Z 3 or u(1) 6 theory. For instance, the only non-diagonal invariants of Z ×2 3 × u(1) 6 correspond to the condensation invariant su(2) 3 × Z 3 or the permutation invariant π(Z ×2 3 ) × u(1) 6 .

The condensation invariants of Z
We begin with the Z ×3 3 CFT, which is realized by the Hamiltonian The presence of third nearest neighbour interactions only means this Hamiltonian corresponds to three decoupled critical 3-state Potts chains. The Z ×3 3 CFT contains primary fields with scaling dimension h = 2 and hence one expects to find a condensation invariant. Indeed, we find, in addition to a permutation invariant to be considered in the next section, a condensation invariant that gives rise to the partition function This partition function corresponds to the diagonal invariant of a c = 12/5 CFT that has 24 primary fields with scaling dimensions given in Table 3. To our knowledge this CFT has not been considered before in literature and we call it here c(Z ×3 3 ). Using the properties (23), we find that Z c Z ×3 3 can be written as the linear combination Thus the boundary condition of one of the chains should depend on the product of the sectors of the other two chains, while the boundary conditions of the other two chains should depend on the sector of this chain only. This can be achieved by adding to (35) the boundary coupling where P n = L/3−1 j=0 Z 3j+n , for n = 1, 2, 3. To write down a local model for the condensation invariant, we employ again canonical duality transformations that now take the form (shown here for a chain of nine sites, the generalization to larger chains is straightforward) In terms of the dual operators, one finds the translationally invariant chain The condensation invariant of u(1) ×3 6 CFT Switching the sign of the Hamiltonian (35), i.e. considering −H Z ×3

The permutation invariants of Z ×3
3 and u(1) ×3 6 CFTs Apart from the condensation invariant, the Z ×3 3 theory also allows for a non-trivial permutation invariant, which is similar to the permutation invariant of the Z ×2 3 theory. To be explicit, in this case it is given by which can be expressed as the linear combination This coupling can be implemented in the three decoupled 3-state Potts models by introducing the boundary term where as above P n = L/3−1 j=0 Z 3j+n . To rewrite this Hamiltonian in an translationally invariant form, one can use the following canonical transformation (given here again for nine sites for clarity) In terms of the dual variables one obtains the local and translationally invariant Hamiltonian The permutation invariant of u(1) ×3 6 CFT The u(1) ×3 6 CFT admits also a permutation invariant corresponding to the partition function χ ϕm 1 χ ϕm 2 χ ϕm 3χ ϕ 3m 1 +2m 2χ ϕ 3m 2 +2m 3χ ϕ 3m 3 +2m 1 The analysis goes again through identically and thus the microscopic model realizing this invariant is given by H π(u(1) ×3 . This applies also to the generic case of π(u(1) ×n 6 ), which means that a permutation invariant can be realized in perturbed cluster clock models H π(u(1) ×n 6 ) = −H π(Z ×n 3 ) . CFTs for generic N . As N increases, the number of possible nondiagonal modular invariants grows quickly. We therefore limit ourselves to two cases, that are straightforward generalizations of the condensation and permutations invariants, as described in Sections 4.3 and 4.4.

Condensation and permutation invariants for Z ×N
The starting point is the Hamiltonian that describes N decoupled 3-state Potts models. As evaluating the modular invariants (8) explicitly become intractable for generic N , we instead consider generalizations of the boundary terms (38) and (43) as ansatzes for non-diagonal modular invariants. These are given by and couples the chains such that the boundary condition of each chain depends on the symmetry sectors of all the other chains, in analogy to the permutation boundary term (43). The can be transformed into a translationally invariant and local form by means of a generalization of the canonical transformations given in Sec. 4.3.
The resulting Hamiltonian is given by which is a generalization of the Hamiltonian (39) for the condensation invariant c(Z ×3 3 ). We verified that this model is again critical and described by a non-diagonal modular invariant for all N . Interestingly, we find that the type of invariant depends on N . If N is a multiple of three, the invariant is equal to a condensation invariant of the CFT consisting of the direct product of N Z 3 parafermion CFTs. This condensation invariant is a straightforward generalization of the condensation invariant (36). The field that condenses is (ϕ 0 2 , ϕ 0 2 , . . . , ϕ 0 2 ), which has scaling dimension h = 2N/3 and which is thus a boson for N mod 3 = 0. If N is not a multiple of three, the invariant is a permutation invariant obtained via the permutation (71) discussed in Appendix A.
By changing the overall sign of the Hamiltonian (49), we obtain a model that realizes a non-diagonal modular invariant of the u(1) ×N 6 CFT. Again, for N a multiple of three, the invariant is a straightforward generalization of the condensation invariant (40) for three 3-state Potts chains. The field that condenses is (ϕ 2 , ϕ 2 , . . . , ϕ 2 ), which has the scaling dimension h = N/3. In the other cases, the invariant is equal to the permutation invariant obtained via the permutation (66) discussed in Appendix A.
Like above, the Hamiltonian H Z ×N can be transformed into a translationally invariant and local form by means of a generalization of the canonical transformation given in Sec. 4.4. The resulting Hamiltonian is given by which is now a generalization of the Hamiltonian (44) for the permutation invariant π(Z ×3 3 ). We have again verified that this Hamiltonian is critical and described by a non-diagonal modular invariant. However, unlike (49), the corresponding invariant is now always a permutation invariant for all values of N . The permutation describing this invariant is given by (70), as discussed in Appendix A. Finally, changing the sign of (50) gives a critical model described by permutation invariants of u(1) ×N

6
CFT for all values of N . The corresponding permutation is given by (65) in Appendix A.
There are intriguing parallels of the termsZ iỸi+1 · · ·Ỹ i+n−2Zi+n−1 to those appearing in generalized spin-1/2 cluster models [6]. To be precise, were theZ i andX i operators replaced by Pauli matrices, then (50) would realize a hierarchy of critical models described by all so(N ) 1 CFTs (that correspond to condensation invariants of Ising ×N CFTs [18]). Whether this similarity is accidental, or whether there is some deeper connection, remains an open question. This similarity suggests also the existence of a clock cluster phase. Were the perturbationsX iX † i+1 +h.c. small compared to the cluster termsZ iỸi+1 · · ·Ỹ i+n−2Zi+n−1 +h.c.
(named as such due them all commuting with each other), one expects a gapped, possibly symmetry-protected topological phase with some number of parafermion edge states [25]. To our understanding such models have not been explored in the literature and it would be interesting to investigate whether they exhibit interesting entanglement properties akin to spin-1/2 cluster states [51].

Beyond 3-state Potts models: Non-diagonal invariants of Z ×2 k CFTs
A natural way to generalize our construction is to go from 3-state Potts models to clock models with Z k symmetry [29]. Such models are obtained by redefining the clock operators Z and X to describe k-state systems. Instead of (2), we demand that they now satisfy Z k = X k = 1, ZX = ωXZ, Z k−1 = Z † = Z −1 and X k−1 = X † = X −1 , where ω = e 2πi/k . The phase diagrams of k-state Potts models were already considered in [29]. It is known that when the relative magnitudes of the different powers of the clock operators are tuned suitably, they exhibit critical points described described by Z k parafermion CFTs, that we discussed in Section 3. For every k, the critical microscopic Hamiltonians take the form [29] which is Hermitian due to the sum over p.
While the k-state models exhibit also other critical points, we focus on this case to show that our construction is not limited to Z 3 parafermions only. Even with the restriction to the Hamiltonian (51), there are many ways in which we can use our construction to create further microscopic chains with particular CFTs. One option is to consider k decoupled chains and condense the boson that is present in the spectrum (see Section 3 for the field content of Z k parafermion CFTs). However, the number of primary fields grows very quickly upon increasing k. Hence, we focus on two decoupled k-state clock chains, described by the diagonal invariant of Z ×2 k CFT, and look for a generalization of the permutation invariant discussed in Section 4.1. As in the preceding Section 4.5, solving for all the modular invariants for a generic k becomes intractable. Thus we again make an ansatz for a suitable boundary term by generalizing the 3-state construction and verify that the resulting model is indeed critical and described by a non-diagonal modular invariant.
We start with the Hamiltonian, which has the Z ×2 k CFT describing its critical point, Motivated by (29), we introduce the boundary term where P n = We verified that this Hamiltonian is indeed critical and described by a non-diagonal modular invariant. Interestingly, we observe that for k odd, the new invariant is a permutation invariant of the Z ×2 k CFT, while for k even, the new invariant is a condensation invariant instead. This behaviour is consistent with the behaviour observed for k = 2 (condensation in Ising ×2 [17]) and k = 3 (Section 4.1). To conclude this section, we give these non-diagonal invariants of Z ×2 k CFT explicitly.

Permutation invariants for odd k
When k is odd, the modular invariant realized by the Hamiltonian Eq. (54) is a permutation invariant, which is a direct generalization of the k = 3 case discussed in Section 3. The chiral fields of the Z ×2 k CFT can be labeled as (ϕ l 1 m 1 , ϕ l 2 m 2 ), where l 1 , l 2 = 0, 2, . . . , k − 1 and m 1 , m 2 = −k + 1, −k + 3, . . . , k − 3, k − 1. The permutation that describes the invariant is given by which gives rise to the permutation modular invariant This modular invariant partition function clearly has a distinct spectrum from the diagonal modular invariant.

Condensation invariants for even k
In the case of k even, the chiral fields of the Z ×2 k CFT can be labeled as (ϕ l 1 m 1 , ϕ l 2 m 2 ), where now l 1 , l 2 = 0, 1, . . . , k/2. In addition, m 1 , m 2 must have the same parity as l 1 , l 2 , and lie in the range 0, 1, . . . , 2k − 1 for l i < k/2 and 0, 1, . . . , k − 1 for l i = k/2. The modular invariant realized by the microscopic model (54) corresponds to a condensation invariant that is obtained by condensing the boson (ϕ 0 k , ϕ 0 k ) with scaling dimension h = k/2, followed by a permutation of the fields. The closed form of the modular invariant partition function depends on the parity of k/2.
Explicitly, for k = 0 mod 4, this condensation invariant takes the following form (we remind the reader that the m label is defined modulo 2k) For both cases the corresponding CFT has (k/2) 2 ((k/2) 2 + k/2 + 2) primary fields and central charge c = 4(k−1) (k+2) . To obtain the modular invariant realized by the Hamiltonian (54), denoted by Z π•c Z ×2 k , one has to perform a further permutation of the anti-chiral fields on these condensation invariants. The needed permutation is the same as the permutation (55) used to construct the permutation invariants for k odd. Applying it to the partition functions above is equivalent to replacing all the anti-chiral charactersχ (ϕ l 1 ) . This does not change the central charge or the number of primary fields and the resulting modular invariants Z π•c Z ×2 k completely specify criticality of the Hamiltonian (54) for all k.

Conclusion
We have constructed several new microscopic clock-type models that by construction exhibit exotic critical points described by CFTs not considered previously in the literature. Some of these correspond to non-trivial permutation invariant partition functions, which shows explicitly that such modular invariants allowed by theory also have realizations in local and translational invariant models. While not limited to them, we focused on generalizations of 3-state Potts models. When these are mapped into parafermion chains with the Fradkin-Kadanoff transformation, our models correspond to parafermion chains with many-body terms. Thus our work addresses also the open question of phase transitions that such many-body terms can drive if they were to appear.
First and foremost, our main result is to demonstrate that exotic CFTs can appear in local and translational invariant models and to construct representative models for each, as summarized in Table 1. The physical realization of any Potts model is challenging, but proposals exist suggesting that domain walls in Abelian fractional quantum Hall edges could realize the parafermion modes [22,23], whose hybridization can result in parafermion chains that are unitarily equivalent to Potts models [14]. If the hybridization were also to give rise to many-body terms, then our critical models could also be realized in this setting. While speculative, if they were, following the same principles to construct a 2D topologically ordered phase from coupled critical 3-state Potts models with Z 3 parafermion criticality, our models could serve as building blocks for even more exotic states of topological 2D matter [14,15].
As we have shown by considering also critical generalizations of k-state Potts models, our method is versatile and lends itself for construction of local models for ever more exotic CFTs. While this is a rather academic exercise with the models increasing in complexity, ideally one would like to have local microscopic realizations for all modular invariants of all CFTs. Admittedly, this is a formidable task as the classification of all modular invariant partition functions is an open question. Still, progress towards this goal can be made with methods of presented here. To go beyond 3-state Potts models, we considered a few examples of models that can be constructed by coupling together two k-state clock models. A particularly interesting case to consider would be to coupling together two such models with opposite sign, generalizing the case of Sec. 4.2. Using the results in [49], i.e. that the central charge for the antiferromagnetic Z k parafermion chains is c = 1, one would expect that the resulting critical point is described by the su(2) k CFT 3 . These in turn can be used to describe minimal models via the coset construction [48], thus possibly enabling a different construction of microscopic realizations for all CFTs in this important class.
While our focus has been on the criticality of the models constructed, it would also be interesting to understand how the criticality relates to the symmetry protected topological phases in the vicinity of the critical points [25,8]. The models with commuting cluster like terms (50) should also be investigated. In spin-1/2 models such terms give rise to entangled ground states that are universal resources for quantum computation [51]. It would be interesting to investigate what, if any, new interesting phenomena appear in the clock model counterparts of cluster states.
Finally, we would like to point out that a generalization of the 3-state Potts model, that is slightly different from the one considered in Sec. 4.1, Eq. (31) was studied in [52]. That model can be obtained simply by exchanging the Z i Z i+1 term in Eq. (31) by Z i Z † i+1 . This cosmetically small change leads to rather different behaviour of the model, for which it turned out hard to fully characterize the critical point. It would be interesting to investigate the relations between these models, if any.

A Fusion rule symmetries of CFTs
In this appendix we discuss the possible permutations of the primary fields in the Z 3 parafermion and u(1) 6 CFTs that leave their fusion rules invariant. Such fusion rule symmetries enable to predict possible non-diagonal modular invariants that may exist in product theories. To this end we treat Z k parafermion CFTs as the cosets Z k ≃ su(2) k /u(1) 2k .
On the other hand, u(1) 2k contains 2k primary fields that we denote by ϕ m , where m is an integer defined modulo 2k. We adopt a convention that m lies in the range −k < m ≤ k. These primary fields obey the fusion rules ϕ m 1 × ϕ m 2 = ϕ m 1 +m 2 .
Treating the Z k parafermion CFT as the coset su(2) k /u(1) 2k means that its primary fields are labeled by the labels of the su(2) k and u(1) 2k theories, namely ϕ l m . The coset means that the labels are no longer independent though, but subject to the constraint l + m = 0 mod 2 as required. and to the field identification ϕ l m ≡ ϕ k−l m+k . With these constraints, the fusion rules follow directly from the fusion rules of the su(2) k and u(1) 2k theories, ϕ l 1 m 1 × ϕ l 2 m 2 = ϕ |l 1 −l 2 | m 1 +m 2 + ϕ |l 1 −l 2 |+1 m 1 +m 2 + · · · + ϕ max(l 1 +l 2 ,2k−l 1 −l 2 ) m 1 +m 2 .
For direct product theories the fusion rules are simply direct products of the fusion rules of the constituent fields. For instance, a product of two Z 3 parafermion CFTs has 36 primary fields (ϕ l 1 m 1 , ϕ l 2 m 2 ). They satisfy the fusion rules where the right hand side is given by all possible product fields compatible with the Z 3 fusion rules.

fusion rules
Let us consider first the symmetries of u(1) 6 fusion rules and product theories of them. For a single u(1) 6 CFT, the only possible permutation invariant is ϕ m → ϕ −m due to the cyclic nature of the fusion rules. This permutation does not lead to a distinct partition function. For a product of two such CFTs, there are more alternatives. All the permutations that either swap the two theories or apply the single theory permutations, leave the fusion rules invariant, but again do not change the partition function. On the other hand, there is also a further allowed permutation given by (ϕ m 1 , ϕ m 2 ) → (ϕ m 1 +2(m 1 +m 2 ) , ϕ m 2 +2(m 1 +m 2 ) ), which results in the permutation invariant partition function discussed in Section 4.1.
The fusion rules of any product theories of three or more u(1) 6 CFTs are invariant when the copies are permuted in an arbitrary way, or under inversion symmetries (63). These permutations do not lead to a change in the partition function. When the number of copies increases, the number of non-trivial permutation invariants grows. We do not give an exhaustive list, but mention two types of non-trivial permutation invariants, relevant for the models considered in this paper.