Classification of spin-$1/2$ fermionic quantum spin liquids on the trillium lattice
Ming-Hao Li Sounak Biswas S. A. Parameswaran





	Rudolf Peierls Centre for Theoretical Physics, Parks Road, Oxford, OX1 3PU, UK


	Rudolf Peierls Centre for Theoretical Physics, Parks Road, Oxford, OX1 3PU, UK
	Institut für Theoretische Physik und Astrophysik, Universität Würzburg, 97074 Würzburg, Germany

	Rudolf Peierls Centre for Theoretical Physics, Parks Road, Oxford, OX1 3PU, UK





§ ABSTRACT 

We study fermionic quantum spin liquids (QSLs) on the three-dimensonal trillium lattice of corner-sharing triangles. We are motivated by recent experimental and theoretical investigations that have explored various classical and quantum spin liquid states on similar networks of triangular motifs with strong geometric frustration. Using the framework of Projective Symmetry Groups (PSG), we obtain a classification of all symmetric 𝖹_2 and 𝖴(1) QSLs on the trillium lattice.  
We find 2 𝖹_2 spin-liquids, and a single 𝖴(1) spin-liquid which is proximate to one of the 𝖹_2 states. The small number of solutions reflects the constraints imposed by the two non-symmorphic symmetries in the space group of trillium.  Using self-consistency conditions of 
the mean-field equations, we obtain the spinon band-structure and spin structure factors corresponding to these states. All three of our spin liquids are gapless at their saddle points: the 𝖹_2 QSLs are both nodal, while the  𝖴(1) case hosting a spinon Fermi surface. One of our 𝖹_2 spin liquids hosts a stable gapless nodal star, that is protected by projective symmetries against additions of further neighbour terms in the mean field ansatz.
 We comment on directions for further work.




Classification of spin-1/2 fermionic quantum spin liquids on the trillium lattice 
    S.A. Parameswaran
    September 9, 2024
==================================================================================






§ INTRODUCTION
 

Spin liquids are magnetic systems that fail to order at the temperatures expected on the basis of their exchange energy scale, while  exhibiting cooperative behaviour that  distinguish them from  high-temperature paramagnet <cit.>. A natural ingredient
leading to such lack of ordering  is geometric frustration, wherein the lattice structure eliminates  simple ground state configurations
which minimise the exchange interaction energy between magnetic
moments <cit.>,  and which would typically lead to symmetry-breaking in the thermodynamic limit. This is  manifest, for instance, in a
system of three classical spins with pairwise antiferromagnetic Heisenberg interactions: the lowest-energy state of the triangle cannot be described in terms of minimal-energy configurations of each of the individual bonds. Frustrated lattices can be assembled by tiling space with such elementary
units — typically triangles or tetrahedra — in order to form
edge-sharing or corner-sharing structures: common examples are the triangular and
kagome lattices in two dimensions (2D), and the pyrochlore and hyperkagome lattices in 3D.
Classical ground states of antiferromagnets on such lattices are macroscopically
degenerate <cit.>. These degeneracies can often be
understood in the exactly solvable large-N limit: frustration is signaled by  a
macroscopically degenerate manifold of continuously connected ordering
wavevectors <cit.>. Such “classical spin liquids” often
order at very low temperatures T (much lower than the scale set by exchange
couplings), in accord with the third law of thermodynamics, that forbids the finite T→ 0 entropy associated with an extensive ground-state degeneracy. Typically, thermal or quantum fluctuations select an ordering
wave vector out of this manifold, in a phenomenon termed “order by
disorder" <cit.>. 
However, in some cases the system is sufficiently frustrated that the quantum mechanical ground state selected  as T→ 0 continues to exhibit no broken symmetries, and instead is a quantum spin characterized by an emergent deconfined gauge structure. The resulting quantum spin liquid (QSL) is often strikingly characterized by the appearance of fractionalized excitations, whereas its gauge structure is more subtly encoded in certain long-range entanglement properties of the ground-state wavefunction. 



A powerful framework to understand  QSL
ground states of quantum spin systems is provided by the projective symmetry group <cit.>. This framework, which  makes the emergent gauge structure especially transparent, builds on the so-called “parton
construction" <cit.>, and proceeds by representing  each spin in terms of auxiliary fermionic `spinons',
S⃗=1/2f^†_i σ⃗_ij f_j. The physical
Hilbert space of quantum spins is recovered via the 
projection i.e. by imposing the constraint that there is exactly one fermion per site. The resulting  Hamiltonian of these
auxiliary (or Abrikosov) fermions is generically quartic and can then be studied within a mean field decoupling wherein 
operators corresponding to fermion hopping, fermion-pair creation, and
fermion-pair annihilation are self-consistently determined, leading to a quadratic mean-field “ansatz". By construction, ground states of such
ansatzes correspond to symmetric, disordered wavefunctions, i.e.,
candidate QSL states. 

This parton (or “projective”) construction suggests low-energy effective descriptions for these phases in terms
of spinons coupled to gauge fields. Of course, the resulting strongly-coupled problem can be challenging to treat in a controlled fashion, particularly in cases where the spinon degrees of freedom are gapless. Nevertheless, a key feature of the parton approach is that it provides a systematic framework to enumerate and classify candidate variational  wavefunctions in terms of their topological structure, in much the same way that the Landau-Ginzburg formalism provides a useful starting point to investigate broken symmetries as captured by conventional  mean-field wavefunctions. Such a classification is facilitated by the projection of the mean-field Hamiltonian from the large Hilbert space
of auxiliary fermions back to the physical spin Hilbert space — essential in order to obtain a sensible spin wave function —
which confers a 
“gauge structure" to the fermion Hilbert space. Specifically, mean field
ansatzes which correspond to the same physical wave function after projection
are related by a gauge transformation. Consequently, the mean field fermion
ansatzes are
only  required to be invariant under physical symmetries up to an associated gauge transformation. In other words, the mean-field ansatz is invariant under a  projective symmetry
group (PSG) which is usually larger
than the physical symmetry group of the QSL wave function. 
However, there exists a group of pure gauge transformations —  typically 𝖹_2, 𝖴(1) or
𝖲𝖴(2)) —  termed the invariant gauge group (IGG),  which leaves the mean field
anstaz invariant. The IGG and PSG together  characterize the low-energy, long-wavelength fluctuations around the mean field ansatz: these involve
fermions coupled to gauge fields, with the gauge group specified precisely by
the IGG, and  fermions in a mean-field dispersion  classified by representations of the PSG.

In other words, different PSGs capture distinct
“quantum orders" of QSL phases with a specified IGG, in much the same way that the physical symmetry groups
characterise broken symmetries. 

Notably, there can be distinct PSGs corresponding to the same physical symmetry manifest in the spin wavefunction, underscoring the fact that these “quantum orders” can be richer than their classical counterparts.


Experimental evidence for QSLs and the resulting need to characterize their emergent low-energy properties has driven a systematic program of applying the projective construction to a variety of frustrated lattices <cit.>. The resulting  mean-field
ansatzes provide starting points for more refined calculations where the
projected fermion wavefunctions can be calculated variationally using Monte
Carlo approaches <cit.>. (Alternative parton constructions that split
the spins into bosons <cit.> offer a complementary set of insights into the phenomenology of possible QSLs and  their possible proximate phases.)




In this work, we continue this program by classifying  symmetric spin liquid states on
the trillium lattice <cit.>, a three-dimensional network of corner sharing
triangles displayed in Fig. <ref>. A natural theoretical motivation of this problem is that the motif of corner-sharing triangles is expected to seed significant magnetic frustration,  like the better-known kagome and hyperkagome lattices. At a more experimentally-grounded level, trillium is  the magnetic lattice of MnSi, or that of the Ce moments in
CeIrSi, which has been considered before in the context of frustrated
magnetism <cit.>. Recent characterisations of the quantum
spin liquid material
K2Ni2(SO4)3 <cit.> show that the magnetic
Ni^2+ ions, with S=1, lie on two interconnected trillium lattices, having
exactly the same set of symmetries as single trillium lattice— implying that
these structures share the classification of symmetric spin liquid states in
terms of projective symmetry groups. Another compound KSrFe2(PO4)3 with structures similar to K2Ni2(SO4)3, with S=5/2, has been shown to exhibit no long range order down to T=0.19K <cit.>. Our interest in trillium is also seeded
by its remarkable similarity to the hyperhyperkagome (HHK) lattice which
describes the network of coupled Cu2+ ions in
PbCuTe2O6, which was  shown to host a QSL
ground state in a series of recent
experiments <cit.>, leading to theoretical
work on spin liquid states on the underlying HHK
structure <cit.>. Both trillium and HHK are
three-dimensional networks of corner-sharing triangles with a cubic Bravais
lattice where each site belongs to three corner-sharing triangles. Classical
frustrated magnets on these lattices share similar
phenomenology <cit.>: a large regime with classical spin
liquid behaviour,  eventually yielding to co-planar ordering at very low
temperatures. For both lattices, large-N approaches yield “partial
ordering" <cit.>, characterized by a
macroscopic but sub-extensive set of ordering wave vectors, whose manifold forms a line (HHK) or surface (trillium) in three
dimensional reciprocal space. This is distinct from the large-N signatures of a classical spin liquid, where this manifold would be extensive <cit.>, as obtained, e.g. for the pyrochlore, kagome and hyperkagome lattices. [Note that a recent study of classical Ising models with  three-spin interactions on the trillium and HHK showed that both  host very similar  classical fractal spin liquids, with “fractonic”  glassy behaviour arising out of kinetic constraints <cit.>; however this is unlikely to be directly relevant to the QSL problem studied here.]
  
  

 The  HHK lattice has the same space group and hence the same classification of PSGs as the three-dimensional hyperkagome lattice <cit.>, in which each site is shared by two (rather than three)
  corner sharing triangles. The latter has been the subject of  numerous theoretical
  investigations of its ordered and spin-liquid
  states <cit.> motivated by its relevance to the candidate QSL material
Na4Ir3O8 <cit.>. The corresponding 𝖯4_1 32 space group  has both 3-fold rotations and a 4-fold
  non-symmorphic screw rotation, with the latter  known to
  cause a drastic reduction of total number of QSL states <cit.>. In contrast, 
  the 𝖯2_1 3 space group of trillium  has a three-fold rotation, along with
  two non-symmorphic screw rotations <cit.>. In the light
  of the preceding discussion, it is natural to ask what QSL phases are
  consistent with symmetries of the trillium lattice. To this end, in this
  paper we undertake a classification of PSGs for spin-1/2 QSLs on this lattice. Although  experiments <cit.> indicate that on K2Ni2SO4 is best understood as an effective spin-1 system, understanding the simpler spin-1/2 case is an important first step towards a more comprehensive study of the higher-spin problem. Accordingly, we
  hope that the present work will guide the
  interpretation of results of future experiments, and add to our theoretical
  understanding of QSL phases in three dimensions.

  The rest of this paper is organised as follows. In Sec. <ref> we introduce the crystal structure of the trillium lattice and the symmetry 
generators of its space group. In Sec. <ref> we present the symmetry group relations of
trillium, and outline the classification of its PSGs using them. We also present the gauge transformations accompanying
physical symmetries for all of the PSGs. The details are relegated to Appendices  <ref> and  <ref>.






§ BACKGROUND: PROJECTIVE SYMMETRY GROUP FORMALISM
 
We briefly review the projective symmetry group classification of parton mean-field theories, as applied to spin models with Heisenberg exchange interactions. Readers familiar with the parton approach can jump ahead to the next section, but may wish to quickly skim this section to orient themselves with our notation and conventions. We begin with the Heisenberg model on a given spatial lattice,

    H=∑_{i, j}J_ijS⃗_i·S⃗_j.


In order to implement the PSG, we first enlarge the Hilbert space by decomposing spins into Abrikosov fermions as follows:

    S⃗_i=∑_α,β1/2f^†_iασ⃗_αβf_iβ.

The above equation maps the spin Hilbert space to the subspace of the Abrikosov fermion Hilbert space in which the fermion occupation number on each site is 1. This means that, on the operator level, we strictly have ∑_αf^†_iαf_iα=Id. Indeed, by using the identity, we can verify that [S^m,S^n]=iϵ_lmnS^l. In fact, a second constraint, is also introduced as a consequence of the first: ∑_α, β f_iαf_iβϵ_αβ = 0.  (One can verify by considering ∑_α, β f_iαf_iβϵ_αβ∑_γf^†_iγf_iγψ, where ∑_γf^†_iγf_iγψ = ψ by virtue of single-occupancy.)


In terms of the Abrikosov fermions, the Heisenberg Hamiltonian reads (up to some constants)

    H   =∑_{i, j}∑_αβμνJ_ij1/4(f^†_iασ⃗_αβf_iβ)·(f^†_jμσ⃗_μνf_jν)
       =∑_{i, j}∑_αβ-1/2J_ij(f^†_iαf_jαf^†_jβf_iβ+1/2f^†_iαf_iαf^†_jβf_jβ),

 We now study H within a mean-field approximation, by introducing parameters for expectation values of operators

    η_ijϵ_αβ=-2⟨ f_iαf_jβ⟩,   χ_ijδ_αβ=2⟨ f^†_iαf_jβ⟩;

where η_ij=η_ji and χ_ij=χ^†_ji.

As is usual, we expand operators in  Eq. <ref> in terms of fluctuations about their expectation values and ignore terms which are quadratic in fluctuations,  leading to 

    H_MFT =    - ∑_⟨ i, j⟩3/8J_ij(χ_ji f^†_iμf_jμ 
       + η_ij f^†_iμf^†_jμ + h.c.
       - |χ_ij|^2 - |η_ij|^2)  
       +∑_i (μ^3_i(f^†_i↑f_i↑ - f_i↓f^†_i↓) 
       + 1/2(μ^1_i+i μ^2_i)f_iμf_iνϵ_μν+ h.c.).
        ,

where we have introduced the Lagrange multipliers μ^1,2,3_i to impose the  
one-fermion-per-site   constraint at a mean-field level. These, as well as the parameters χ_ij and η_ij, are determined self consistently.

To discuss the 𝖲𝖴(2) gauge structure of the mean-field Hamiltonian, it is convenient to introduce a spinor representation

    ψ≡[ ψ_1; ψ_2 ]≡[   f_↑; f^†_↓ ].



In terms of these spinors, the  H_MFT can be be compactly rewritten as:

    H_MFT   = ∑_⟨ i, j⟩3/8 J_ij[1/2Tr(U^†_ijU_ij) - (ψ^†_i U_ijψ_j + h.c.)] 
                   + ∑_i,lμ^l_iψ^†_iτ^l ψ_i ,

where the matrix U_ij captures both mean-field parameters via 

    U_ij≡[ χ^†_ij   η_ij; η^†_ij  -χ_ij ] .

The constraint implementing projection into the spin Hilbert space at the mean-field level now has the form:

    ⟨ψ^†_i τ^l ψ_i ⟩ = 0,      l=1,2,3.

The {U_ij} and {μ^m_i } together constitute variational parameters that specify the mean-field “ansatz”for the  Hamiltonian and the  corresponding ground state wavefunction. Variationally optimizing the parameters to obtain the lowest energy ground state is equivalent to determining the parameters self-consistently.

It is crucial to realize that the ground state of the mean field spinon Hamiltonian is not a valid spin wave function, since the on-site constraints are only enforced on average. The final  wavefunction in terms of the physical spin degrees of freedom is constructed from the mean-field spinon state by Gutzwiller projection, i.e. Ψ_spin = P_GΨ_MFT.

The spinor representation makes the 𝖲𝖴(2) gauge redundancy  of the
mean-field Hamiltonian manifest. The Hamiltonian is invariant, trivially, under
the site-dependent gauge transformation ψ_i↦ W_iψ_i, 
U_ij↦ W_i U_ij W_j^† and μ^m_i ↦ W_i μ^m_i W_i^†, with W_i∈𝖲𝖴(2) since this leaves physical spin operator invariant [cf.  <ref>]. Therefore, the mean-field anstaze
parametrised by U_ij,μ^m_i and  W_i U_ij
W_j^†, W_i μ^m_i W_i^† share the same physical spin wavefunctions, i.e. after projection into the
 spin Hilbert space. This has significant consequences for what we
require of symmetric mean-field ansatzes. Consider the action of a symmetry g:
U_ij↦ U_g(i)g(j). For a symmetric ansatz we no longer
require U_g(i)g(j)=U_ij; instead, we only require  that   there exist transformations G_g(n) ∈𝖲𝖴(2) for all sites n, such that G_g(g(i)) U_g(i)g(j)G^†_g(g(j)) = U_ij. The gauge redundancy then implies that physical properties of the state represented
by the ansatz has not changed. The physical transformations g together with the gauge transformation, (G_g(i),g), which leaves
the ansatz invariant, constitute the projective symmetry group (PSG).
The PSG characterises the symmetries of the ansatz, and serves to classify and
characterise different mean field spin liquid states. The PSG also determines the low-energy  description of fluctuations about the
mean field states.  From the preceding discussion on the gauge structure it is
clear that not all fluctuations of the mean-field parameters { U_ij} are
physical: the unphysical fluctuations between gauge inequivalent states must be
described by gauge fields in the effective theory. The effective theories,
then, are likely to be fermions coupled to gauge fields. The gauge structure of
the low energy theory is in general not   given by the high energy
gauge group 𝖲𝖴(2), but is instead determined by the
“invariant gauge group” (IGG) <cit.>. The IGG is a subgroup of the PSG
comprised of pure gauge transformations which leave the ansatz invariant,
i.e., 𝒢 = { W_i|W_i U_ijW^†_j=u_ij, W∈𝖲𝖴(2)}.  Given the central importance of the IGG, one usually labels
QSLs by the IGG, leading to the terminology of “𝖹_2,
𝖴(1), or 𝖲𝖴(2)" QSLs. The PSGs, therefore, play a role
for mean-field QSL phases akin to that of ordinary symmetry groups for broken-symmetry phases, distinguishing quantum
disordered states with the same physical symmetries but different emergent properties.This is a good place to flag one final complication: namely, that that certain PSGs do not correspond to non-zero mean field ansatzes. Therefore, simply tabulating the list of PSGs does not conclude the classification of PSGs; it is essential to investigate the physical constraints that each imposes on the mean field ansatzes. Despite this complication, it is nevertheless useful to organize the investigation of symmetric spin liquid ground states on a
given spatial lattice in terms of an enumeration of all PSGs, for a given set of physical
symmetries (typically, the full lattice space group as well as time reversal symmetry) and the IGG. This
allows the construction of the corresponding mean-field ansatzes and  spin liquid wavefunctions.  In the balance of this paper, we implement such a program
 for the trillium lattice.



§ PSGS OF THE TRILLIUM LATTICE
 









 §.§ The trillium lattice



We begin by describing the trillium lattice and its spatial symmetries. These, along with time reversal, will constitute the physical symmetries that our QSL ground states (after projection to the correct Hilbert space) must respect, and are hence central to the PSG construction.

The trillium lattice, shown in Fig. <ref>,  has a simple cubic Bravais lattice with four sub-lattices: α, β, γ and  δ. The positions of the sublattice sites relative to the unit cell center are given by:

    r⃗^0_α   = (κ , κ , κ ), r⃗^0_β = (1/2+κ , 1/2-κ , 1-κ), 
    r⃗^0_γ   = (1-κ , 1/2+κ , 1/2-κ), r⃗^0_δ = (1/2-κ , 1-κ , 1/2+κ),

where κ is a free parameter. As mentioned before, the nearest neighbor bonds on the lattice form a network of corner sharing triangles, with each site participating in three triangles, which are the elementary frustrated motifs.

We  denote the  position of a unit cell i by the vector u⃗_⃗i⃗=(x,y,z), where x,y,z are integers.  
A generic lattice site i  is referred to by specifying its unit cell position and sublattice as i≡ (x,y,z;s); such a site lies at position u⃗_i+r⃗^0_s. Since the mean-field parameters {U_ij} specifying the ansatz are associated with the links, it is  convenient to uniquely label all links for the purpose of further discussion. We do so by exploiting lattice translation invariance: there are 12 links per unit cell, all of which are translationally inequivalent. We introduce the labels ζ=(1,2… 12) for these links, and specify each of these links in Tab. <ref>.

Trillium has  space group 𝖯2_1 3, with the symmetry generators {T_x, T_y, T_z, g_a, g_b, g_c }. Here, T_is are the three translational generators. g_c is a threefold rotation about the (1,1,1) axis passing through the origin of an unit cell. g_a and g_b are the generators of the 2-fold non-symmorphic screw rotations. g_a involves a π rotation about an axis in the (0,0,1) direction passing through the point (1/2,0,0), followed by a translation by 1/2 of the unit-cell distance along the rotation axis. g_b involves a π rotation about an axis in the direction (0,1,0) 
passing through the point (0,0,1/2) followed by a translation of 1/2 of the untit-cell distance along the rotation axis.  It has been noted previously <cit.> that non-symmorphic symmetries generally lead to strong constraints on possible PSGs, and a consequent reduction of their number.

The generators g_a,g_b and g_c act on a lattice site i≡(x,y,z,s) via

    g_a:     (x,y,z;α)↦ (-x, -y-1, z;δ), 
        (x,y,z;β)↦ (-x-1, -y-1, z+1;γ),
        (x,y,z;γ)↦ (-x-1, -y-1, z;β),
        (x,y,z;δ)↦ (-x, -y-1, z+1;α),
    
    g_b:     (x,y,z;α)↦ (-x-1, y, -z;γ),
        (x,y,z;β)↦ (-x-1, y, -z-1;δ),
        (x,y,z;γ)↦ (-x-1, y+1, -z;α),
        (x,y,z;δ)↦ (-x-1, y+1, -z-1;β),
    
    g_c:     (x,y,z;α)↦ (z, x, y;α),
        (x,y,z;β)↦ (z, x, y;γ),
        (x,y,z;γ)↦ (z, x, y;δ),
        (x,y,z;δ)↦ (z, x, y;β).




 §.§ PSG classification on the trillium lattice


The PSG involves the group of the transformations (G_g(n),g) that leaves
the mean-field ansatz invariant. Here g is a physical symmetry
transformation, and G_g(n) ∈𝖲𝖴(2) is the associated
site-dependent gauge transformation, with n denoting the physical site.  (G_g(n),g) acts on a mean-field
parameter  U_ij as 

    (G_g(n),g): U_ij↦ G_g(g(i)) U_g(i)g(j)G^†_g(g(j)).

It follows from consecutive action on the ansatz that the product of two PSG elements is given by the group compatibility condition

    (G_g_1(n),g_1) ∘ (G_g_2(n),g_2) = (G_g_1(n)G_g_2(g^-1_1 n),g_1 g_2),

and the group inverse by

    (G_g_1(n), g_1 )^-1  = (G^†_g_1(g_1(n)),g^-1_1 ).

Since the
elements of the IGG 𝒢 are pure gauge transformations which leave
the ansatz invariant, it is clear that whenever (G_g(i),g) is an element of
the PSG, (W G_g(i),g), for all W ∈𝒢, is also an element of the
PSG. If one considers a gauge-equivalent ansatz, W_i U_ij W^†_j, the
PSG element (G_g(i),g) changes to (W_i G_g(i) W^†_g(i)). PSGs related
  by such gauge transformations are equivalent; they are associated with
  gauge-equivalent ansatzes and represent the same QSL phase. Our task is to
  find all such  equivalence classes; in other words, to single out one 
  representative from each class by fixing the gauge freedom.

  It is convenient to carry out this task purely “algebraically”,
  i.e., by making no reference to the ansatz. To do this, we note that given a
  physical symmetry group and the IGG 𝒢, the PSG can be viewed as a
  group equipped with a projection 𝒫 to the physical symmetry
  group, such that 𝒫: (G_g(i),g) ↦ g. From the discussion in
  the previous paragraph, 𝒫: (W G_g(i),g) ↦ g for W ∈𝒢 . As a corollary, 𝒫 projects pure gauge
  transformations in the IGG back to the identity element, 𝒫: (W ,e)
  ↦ e for W ∈𝒢. 

The projection map between the PSG and the physical symmetry group implies that
the gauge transformation G_g(i) associated with the symmetry transformation
g is constrained by the relations between symmetry group elements g.  
These constraints on
G_g can be used to enumerate all gauge-inequivalent choices of G_g for
all symmetry transformations g, and hence enumerate all PSGs.

To see this, one begins with the relations between the symmetry generators {T_x,
T_y, T_z, g_a, g_b, g_c, 𝒯} which completely specify the group.
Each such relation will lead to an equation constraining the associated PSG
elements. The minimal set of such relations that specify the group is called
the “presentation” of the group. Using the GAP computer algebra package <cit.>, we obtain the
finite presentation of trillium space group of : 

    g^3_c=e,
        T^-1_zg^2_a=e,
        T^-1_yg^2_b=e,
        T^-1_xT^-1_yT_xT_y=e,
        T^-1_yT^-1_zT_yT_z=e,
        T^-1_zT^-1_xT_zT_x=e,
        g^-1_aT_xg_aT_x=e,
        g^-1_aT_yg_aT_y=e,
        g^-1_aT^-1_zg_aT_z=e,
        g^-1_bT_xg_bT_x=e,
        g^-1_bT^-1_yg_bT_y=e,
        g^-1_bT_zg_bT_z=e,
        g^-1_cT^-1_yg_cT_x=e,
        g^-1_cT^-1_zg_cT_y=e,
        g^-1_cT^-1_xg_cT_z=e,
        g^-1_ag^-1_cg^-1_bT^-1_xT_yg_ag_c=e,
        g^-1_bg^-1_aT_xT^-1_yT_zg_bg_a=e,
        g^-1_cg^-1_bT^-1_xT_yg_ag_bg_cg_b=e,


where e denotes the identity of the symmetry group.We focus further on QSLs on the trillium lattice which respect time-reversal symmetry (TRS).  
The TRS operator 𝒯 acts on the mean-field ansatz by complex conjugating  the mean-field parameters
U_ij and μ_i. It is convenient to include a global gauge transformation i τ_2 in the definition of
G_𝒯, such that we have 

    (G_𝒯,𝒯):  U_ij↦   G_𝒯(i) i τ_2 U^*_ij (-i τ_2) G^†_𝒯(j)
     
     =    -G_𝒯(i) U_ij G^†_𝒯(j).

Including  TRS introduces the following additional relations, which express the fact that 𝒯 commutes with generators in the space group:


    𝒯^2=e, 
       𝒯^-1T^-1_x𝒯T_x=e,
       𝒯^-1T^-1_y𝒯T_y=e,
       𝒯^-1T^-1_z𝒯T_z=e,
       𝒯^-1g^-1_a𝒯g_a=e,
       𝒯^-1g^-1_b𝒯g_b=e,
       𝒯^-1g^-1_c𝒯g_c=e.


Chiral spin liquids, which break TRS and some lattice symmetries separately
while preserving their combinations,  have also been considered in the
literature <cit.>. For chiral PSGs, one considers the symmetry
group generated by g 𝒯^ϵ_g instead of the usual symmetry group generated by {g} <cit.>. 
 ϵ_g={0,1} specifies whether the lattice symmetry g is preserved on its own (ϵ_g=0), or 
preserved only up to TRS (ϵ_g=1). 
The trillium SG relations given by Eqs. <ref>-<ref> 
impose the constraint ϵ_g=0 for all generators. This can be easily seen from the fact that
for each generator g, there exists one SG relation which has only an odd number of appearances of that generator, which forces ϵ_g=0. One could still consider spin liquids which respect all lattice symmetries but not TRS. In all our PSG calculations, we first derive the PSG classification without TRS, and then impose TRS at the end.  While this immediately gives us  the PSGs without TRS, we forego a consideration of mean-field ansatzes corresponding to such PSGs, restricting ourselves to the study of physical fully symmetric spin liquids.
Ground states for classical spins on the trillium lattice <cit.>
are also known to be non-chiral (which, for classical spin configurations, is equivalent to co-planarity).


The projective
relation between the symmetry group elements and the corresponding PSG elements
allow us to translate the above symmetry relations
(Eqs <ref>- <ref>) into constraint equations for the  PSG elements.
Consider a general symmetry group relation among a set of elements, ∏_νg_ν= e. The product of the corresponding 
PSG elements are given by (G̃, ∏_νg_ν=e), where 
  G̃ can be constructed from the matrices G_g_ν(i) using
  Eq. <ref>.
Under the projection 𝒫 to the symmetry group elements (G̃,e)↦ e); this
  immediately implies a constraint equation expressing that G̃ must be a member of the IGG, G̃∈𝒢.

The unknowns in these
equations are of two kinds: first, the site-dependent 
gauge transformation matrices { G_x , G_y, G_z, G_a, G_b, G_c, G_𝒯} 
accompanying each symmetry transformation in { T_x, T_y,
T_z, g_a, g_b, g_z, 𝒯}; 
and second, an
element of the IGG W∈𝒢 corresponding to each symmetry group relation in 
Eqs <ref>-<ref>. Solving
these equations, along with choice of gauge described earlier, leads to
the different inequivalent PSGs. 

Explicit   procedures for solving these equations in a fixed gauge are detailed for specific lattices in
Ref. <cit.>, as well as several later works that classifying PSGs in
different spatial lattices <cit.>. We have
undertaken this procedure to enumerate and classify all symmetric spin liquids
with the IGG set to both 𝖹_2 and 𝖴(1). The calculations
are tedious, and  hence we have relegated their details to the
Appendices <ref> and  <ref> for conciseness. Each inequivalent PSG is uniquely specified by  the 
expressions for site-dependent
gauge transformations { G_x , G_y, G_z, G_a, G_b, G_c, G_𝒯}
which accompany the symmetry transformation. We now summarize our results by specifying these gauge transformations for all the PSGs that we identify.

When the IGG is fixed to 𝖹_2, we find 4 inequivalent PSGs.  Once the global
gauge freedoms are fixed, the gauge transformation matrices associated with
lattice translations are uniform, with no position or sublattice
dependence for all PSGs, i.e., G_x=G_y=G_z=1. The PSGs can be
uniquely indexed by constraints on the gauge transformation matrices obtained
from the PSG equations. First, the transformation corresponding to
time-reversal G_𝒯 takes the values τ_0 or i τ_z, though
the PSGs corresponding to G_𝒯=τ_0 do not lead to any non-zero
mean-field ansatzes. Gauge transformation matrices associated with other
symmetry generators are also unit-cell independent, although they retain a
sublattice dependence. Second, the gauge transformation associated with the
rotation g_c acting on sites of sublattice α,
G_c(α)=𝒜,  takes the 2 values exp(i k(2π/3)  τ_z)
for k={0,1}. All other gauge transformations can be specified in
terms of these three, as detailed in Tab. <ref>. The 2
possible values of G_c(α) and
the 2 possible values of G_𝒯 lead to 4 inequivalent PSGs, out
of which only 2 (corresponding to G_𝒯=i τ_z) lead to
non-zero mean field ansatzes. 


Next, we fix the IGG to 𝖴(1). As in the case
of 𝖹_2, the gauge transformations corresponding to the three translations
are uniform, G_x=G_y=G_z=1. The other gauge transformations, however, acquire 
both a unit-cell and a sublattice dependence. The PSGs can again be indexed by the parameters
specifying certain gauge transformations. The gauge transformation associated with the
rotation g_c acting on sites of sublattice α,
G_c(α)=𝒜,  takes the 2 values exp(i k(2π/3)  τ_z)
for k={0,1}. On including time reversal,
we find two possibilities for the associated gauge-transformation G_𝒯. First, G_𝒯 can be iτ_z uniformly, and this case does not lead to any physical spin-liquid ansatz with non-zero mean-field parameters, and so we do not consider these PSGs further. Second, G_𝒯 can acquire a space-dependent form depending on A, which leads to physical spin liquids. Following the second possibility, therefore, we have 
1 𝖴(1) spin liquid, as we will later show that only the k=0 case leads to nearest neighbor mean field spin liquid states. We specify the PSGs by expressing all gauge transformations in terms of the parameters A in Tab. <ref>. In the next section, we will  construct mean-field ansatzes for spin liquids corresponding to these PSGs and proceed to investigate them.





§ MEAN-FIELD SPIN LIQUID PHASES
 

    



In this section, we construct the mean field QSL solutions to the Heisenberg Hamiltonian on the trillium lattice. For this purpose, we will consider only  nearest-neighbor interactions and set J=8/3  henceforth.  The form of the ansatzes are constrained by the PSG. 
Concretely, a given PSG   (G_g,g),g∈𝖯2_13×𝖹^𝒯_2 requires that,

    ∀ g:      G_g(g(i)) U_g(i)g(j)G^†_g(g(j)) = U_ij,
        G_g(g(i)) μ_g(i)G^†_g(g(i)) = μ_i.


The  ansatz for each PSG is derived by systematically imposing Eq. <ref> using the 
gauge transformations detailed in Tab. <ref> and Tab. <ref>; the results of this procdure, detailed in Appendix <ref>, are  tabulated in Table <ref>. The labeling scheme in the table is such that, Given a generic label 𝖹_2x, we can read off the phases in our PSG solutions: 𝒜 = exp(ix2π/3τ_z). There is only 1 𝖴(1) QSL, which is labeled as such. The primed parameters, i.e. quantities like U^x', are defined as follows:

    [ U^x'; U^y' ]   = 1/2[    -1  √(3); -√(3)    -1 ][ U^x; U^y ],
    [ U^x”; U^y” ]   = 1/2[    -1 -√(3);  √(3)    -1 ][ U^x; U^y ].





Concretely, this means we need to find {U_ij, μ_i} such that the following self-consistency equations and on-site constraints are satisfied:

    χ_ij   =⟨ f^†_i↑f_j↑⟩+ ⟨ f^†_i↓f_j↓⟩, 
    η_ij   =⟨ f_i↓f_j↑⟩- ⟨ f_i↑f_j↓⟩, 
    
        1     = ⟨ f^†_i↑f_i↑⟩+ ⟨ f^†_i↓f_i↓⟩, 
    
        0     = ⟨ f_i↑f_i↓⟩.

We can impose symmetry conditions on the mean-field solutions to reduce the number of parameters, and the PSG determines these symmetry conditions by requiring that Eq. <ref> is respected. 
We discuss how the symmetry conditions constrain the mean-field ansatzes in detail in
Appendix <ref> and Appendix <ref>, and the results are tabulated in Tab. <ref>.
Once we obtain the form of the mean-field Hamiltonians, we assemble the
Hamiltonian using our variational parameters, and solve the non-linear equation set Eq. <ref> using the NLSolve package <cit.> available in
julia <cit.>. We set up our system with periodic boundary conditions (PBC),
with L=99 in the three directions. 

The numerical values of mean-field parameters and the corresponding energies obtained from the self-consistent solutions are 
summarized in Table <ref>. Note that the  mean-field energies of 𝖹_20 and 𝖴(1) state are the same, while 𝖹_21 has a higher energy. A comment on the close energies of the 𝖹_20 and 𝖴(1) state is in order. The ansatzes corresponding to these states are uniform, i.e, U_ij=U^x τ_x +U^y τ_y (μ^i= μ^x τ_x + μ^y τ_y) for all links (sites) of the 𝖹_20 state while U_ij=U^z τ_z (μ_i = μ^z τ_z) uniformly on all links (sites) of the 𝖴(1) state. The saddle-point  solutions for the mean-field parameters of the 𝖹_20 state have the property |μ^x/μ^y| = |U^x/U^y|. Thus we can use a global 𝖲𝖴(2) transformation to change this ansatz — only at the saddle point — to the 𝖴(1) ansatz displayed in Table  <ref>. Similar phenomenon has been noted in previous works such as <cit.>.
Therefore, while the nearest-neighbor ansatzes for both QSLs are the same at the saddle point, the general ansatzes displayed in Table  <ref> refer to different QSL states with different IGGs. 


    



 §.§ Relations between the QSLs 

We can also infer connections between the 3 distinct QSLs. First, we note that the 𝖴(1) QSL is the parent state of the 𝖹_20 QSL.  This can be seen by first performing a gauge transformation on 𝖴(1):

    W(α)= W(β)= W(γ)= W(δ)= e^-iπ/4τ_y,

so that for this QSL we now have

    U_ζ = U^z τ_x,   ζ∈{1,…,12 },
       μ_s  = μ^z τ_x,    s∈{α,…,δ}.

From this, we see that introducing the perturbations Δ U_ζ∼τ_y and Δμ_s ∼τ_y breaks the 𝖴(1) symmetry down in a manner that results in the 𝖹_20 state. The 𝖹_21 QSL state cannot be obtained by perturbing around the 𝖴(1) QSL state.



 §.§ Spinon spectra and the nodal star

Our mean-field ansatz allow us to determine the structure of excitations on the mean-field ground state. We compute the eigenvalues of the mean-field Hamiltonian for different wave-vectors to obtain the spinon dispersion spectra in the Brillouin zone.

Fig. <ref> shows the spinon band structures for the different mean-field QSL states. As explained earlier, the nearest-neighbour ansatzes for the 𝖹_20 and 𝖴(1) states are related by a global 𝖲𝖴(2) at the saddle point, leading to identical band structures.
To illustrate the structure of possible gapless modes, we plot the set of gapless points in the Brillouin zone (BZ) for 𝖴(1) and 𝖹_21 states in Fig. <ref>. The collection of gapless points at the saddle point for the 𝖹_20 QSL is identical to that of the 𝖴(1) QSL.

We note that all the mean field states we obtained are gapless at saddle point. The 𝖴(1) possesses a spinon fermi surface at the center of the BZ, with another sheet of spinon fermi surface at the corners.
The 𝖹_2 mean-field states are also gapless. The 𝖹_20 state has similar spectrum to those of the 𝖴(1).  We note that the gaplessness of the 𝖹_20 state seems to not be protected by symmetries and one might generically expect a gap to open up when our nearest-neighbor ansatzes are extended to include further neighbor terms.


The 𝖹_21 state hosts a spectrum with a “nodal star" of gapless points, with dispersion-less bands running from the center of the BZ to its 8 corners. This can be seen from Fig. <ref>.  This nodal star is not a specific property of the short-range ansatz we use to display the bands in Fig. <ref>; rather it is robust to the addition of arbitrary links in the ansatz. In App. <ref> we prove that the gapless nodal star is protected by projective symmetries of the 𝖹_21 phase. 
Such gapless nodal stars have received significant attention in the pyrochlore lattice <cit.>, where two gapless bands along the nodal star were recently proven to be protected by the projective symmetries <cit.>. Such lines were also observed in FCC structures in Ref. <cit.>, where the whole mean field Hamiltonian vanishes along the nodal star. Gapless nodal loops were observed in diamond lattice  <cit.> where strong evidence of symmetry-protection was provided by showing that the  gapless nodal loops  persist despite longer range bond amplitudes being  included in the ansatz.

Our proof of the protected nodal star is algebraic and close in spirit to that of Ref. <cit.>. We  look at the symmetries of the mean-field Hamiltonian directly in momentum-space

    H_ MFT=∑_k⃗ψ^†(k⃗) H_ MFT (k⃗) ψ (k⃗),

When the spinors (Eq. <ref>) are arranged as ψ(k⃗)=(ψ^α_1(k⃗),ψ^β_1(k⃗),ψ^γ_1(k⃗),ψ^δ_1(k⃗),ψ^α_2(k⃗),ψ^β_2(k⃗),ψ^γ_2(k⃗),ψ^δ_2(k⃗)), time-reversal already implies that H_ MFT(k⃗) takes the form 
  
    H_ MFT(k⃗)= 
      [       0_4×4 h_4 × 4(k⃗); h^†_4×4(k⃗)       0_4×4 ].

This block off-diagonal hermitian structure implies that the eigenvalues come in symmetric pairs of ± E(k⃗) everywhere in the BZ. Next, we work out the most general form of h_4 × 4(k⃗) allowed by the projective representations of
the symmetries (G_a,g_a), (G_b,g_b) and (G_c,g_c). 
Restricting the general form of h_4×4(k⃗) to the “nodal star” wavevectors  k⃗=(± k, ± k, ± k), we show, using elementary linear algebraic techniques,  that it has a maximum rank of 3. This implies that H_ MFT has a maximum rank of 6 along the nodal star, proving the existence of two gapless bands. The complete proof  involves explicit expressions for the most general h_4 × 4(k⃗) allowed by projective symmetries, and is fleshed out in Appendix <ref>.By computing the equilibrium state energy of the mean field spinon models, we estimate the temperature dependence of the specific heat for the nodal star spin liquid state. Specifically we see that for 𝖹_20, C_v∼ T^1.22, where as for 𝖹_21, C_v∼ T^0.73. The numerical results are given in Fig. <ref>. The above analysis is not performed for the 𝖴(1) spin liquid, since we expect the gauge field excitations at low energies to modify the results from the calculations of the non-interacting model.

    







 §.§ Spin structure factors


    


In Fig. <ref>, we plot the static structure factor of 𝖹_20, 𝖹_21 and 𝖴(1) mean field states in the k_y-k_z plane. The definition of the static structure factor is:

    𝔖^s_i,s_j(q⃗) ≡1/N∑_R⃗ e^-iq⃗· (R⃗ + d⃗^0_ij)⟨S⃗_(0;s_i)·S⃗_(R⃗;s_j)⟩,

where s_i and s_j are the sub-lattice indices of site i and j, R⃗ is the distance between the two unit cells, and d⃗^0_ij is the distance between the two sub-lattice sites within in the unit cell. And we compute the sum of all these components:

    𝔖(q⃗) ≡∑_s_i, s_j𝔖^s_i,s_j(q⃗),

then plot the normalized results.The structure factor plots indicate that we have obtained remarkable quantum spin liquid states, especially the 𝖴(1) QSLs, which are visibly featureless, implying a sharp departure from the ordered states. We also note that the 𝖹_21 state is the most featureful among the three. 



§ CONCLUSIONS AND OUTLOOK
 
In this work, we have computed the PSGs for the trillium lattice both with and without time reversal symmetry. In the former case we implement the full  construction of the nearest neighbor mean field (fermionic) parton Hamiltonian of the corresponding quantum spin liquid states.  We find two distinct such QSLs with a  𝖹_2 gauge group, and a single example of a QSL with a  𝖴(1) gauge group. We also obtained the corresponding mean-field spinon  band structures and static structure factors, providing some basic thermodynamic and spectral information on these states.  Our main results  are reported in Tab. <ref> and Tab. <ref>.

As noted in the introduction, one of our principal motivations is the recent report of QSL-type behaviour in 
K2Ni2SO4; our work represents a stepping stone towards a parton mean-field analysis of this system, which hosts a double trillium lattice with spin-1 moments. Accordingly, a natural next step in this program is to modify the PSG analysis to account for these differences, and 
perform variational Monte Carlo studies of the Gutzwiller projected mean field QSL wavefunctions to compare with the available experimental data. These tasks are currently underway, and we hope to report on them in the near future. 


SAP acknowledges useful conversations with Yasir Iqbal, and the hospitality of the Max Planck Institute for Complex Systems (MPI-PKS) Dresden, where part of this work was completed. SB acknowledges useful discussions with Atanu Maity.
This work was supported by the European Research Council under the European Union Horizon 2020 Research and Innovation Programme, Starting Grant [Agreement No. 804213-TMCS],  a Buckee Scholarship at Merton College (ML), the Deutsche Forschungsgemeinschaft via Grant No. AS120/16-1,  Project No. 493886309 (SB) and the Gutzwiller Fellowship of the MPI-PKS (SAP)














§ IGG=𝖹_2
 
One can translate the SG relations to the PSG relations:


    G_c(g_c^3(i))G_c(g_c^2(i))G_c(g_c(i))=η_c,
        G^†_z(g_a^2(i))G_a(g_a^2(i))G_a(g_a(i))=η_a,
        G^†_y(g_b^2(i))G_b(g_b^2(i))G_b(g_b(i))=η_b,
        G^†_x(T^-1_yT_xT_y(i))G^†_y(T_xT_y(i))G_x(T_xT_y(i))G_y(T_y(i))=η_xy,
        G^†_y(T^-1_zT_yT_z(i))G^†_z(T_yT_z(i))G_y(T_yT_z(i))G_z(T_z(i))=η_yz,
        G^†_z(T^-1_xT_zT_x(i))G^†_x(T_zT_x(i))G_z(T_zT_x(i))G_x(T_x(i))=η_zx,
        G^†_a(T_xg_aT_x(i))G_x(T_xg_aT_x(i))G_a(g_aT_x(i))G_x(T_x(i))=η_ax,
        G^†_a(T_yg_aT_y(i))G_y(T_yg_aT_y(i))G_a(g_aT_y(i))G_y(T_y(i))=η_ay,
        G^†_a(T^-1_zg_aT_z(i))G^†_z(g_aT_z(i))G_a(g_aT_z(i))G_z(T_z(i))=η_az,
        G^†_b(T_xg_bT_x(i))G_x(T_xg_bT_x(i))G_b(g_bT_x(i))G_x(T_x(i))=η_bx,
        G^†_b(T^-1_yg_bT_y(i))G^†_y(g_bT_y(i))G_b(g_bT_y(i))G_y(T_y(i))=η_by,
        G^†_b(T_zg_bT_z(i))G_z(T_zg_bT_z(i))G_b(g_bT_z(i))G_z(T_z(i))=η_bz,
        G^†_c(T^-1_yg_cT_x(i))G^†_y(g_cT_x(i))G_c(g_cT_x(i))G_x(T_x(i))=η_cyx,
        G^†_c(T^-1_zg_cT_y(i))G^†_z(g_cT_y(i))G_c(g_cT_y(i))G_y(T_y(i))=η_czy,
        G^†_c(T^-1_xg_cT_z(i))G^†_x(g_cT_z(i))G_c(g_cT_z(i))G_z(T_z(i))=η_cxz,
        G^†_a(g^-1_cg^-1_bT^-1_xT_yg_ag_c(i))G_c^†(g^-1_bT^-1_xT_yg_ag_c(i))
       × G_b^†(T^-1_xT_yg_ag_c(i))G^†_x(T_yg_ag_c(i))
       × G_y(T_yg_ag_c(i))G_a(g_ag_c(i))G_c(g_c(i))=η_acb,
        G^†_b(g^-1_aT_xT^-1_yT_zg_bg_a(i))G^†_a(T_xT^-1_yT_zg_bg_a(i)) 
       × G_x(T_xT^-1_yT_zg_bg_a(i))G^†_y(T_zg_bg_a(i))
       × G_z(T_zg_bg_a(i))G_b(g_bg_a(i))G_a(g_a(i))=η_ab,
        G^†_c(g^-1_bT^-1_xT_yg_ag_bg_cg_b(i))G^†_b(T^-1_xT_yg_ag_bg_cg_b(i))
       × G^†_x(T_yg_ag_bg_cg_b(i))G_y(T_yg_ag_bg_cg_b(i))G_a(g_ag_bg_cg_b(i)) 
       × G_b(g_bg_cg_b(i))G_c(g_cg_b(i))G_b(g_b(i))=η_cba.


The Gs in the above relations are 𝖲𝖴(2) matrices, and are associated with the 𝖲𝖴(2) gauge symmetry, which transforms the G in the following way:

    G_g(i)↦ W(g^-1(i))G_g(i)W^†(i),    W∈𝖲𝖴(2);

This equation can be understood as follows: under gauge transformation, we have:

    U_ij↦Ũ_ij≡ W(i)U_ijW^†(j),

and the requirement for the gauge transformed PSG is:

    G̃_g(g(i))Ũ_g(i)g(j)G̃^†_g(g(j)) = Ũ_ij.

From the above relations we derived the gauge transformation of Gs.Aside from the gauge symmetry, we note that we can replace a generic element G with 𝔤G, where 𝔤∈IGG={τ_0, -τ_0 }. Wisely making use of this fact is going to help us reduce the number of phases on the right hand side of the PSG equations.By performing:

    G_x↦η_cyx G_x,    G_z↦η_czy G_z, 
       G_a↦η_acbη_cbaG_a,    G_b↦η_acbη_cyx G_b, 
       G_c↦η_c G_c,

we eliminate the phases on the right hand side (RHS) of Eq. <ref>, Eq. <ref>, Eq. <ref>, Eq. <ref> and Eq. <ref>.



 §.§ Solving for the translational Elements

Let us start by considering the following equations Eq.<ref>, Eq.<ref> and Eq.<ref> that arise because of the commutation of translational generators. Canonically, this gives us the following expressions of G_x, G_y, G_z:

    G_x(x,y,z;s)   =τ_0, G_y(x,y,z;s)=η^x_xyτ_0, 
     G_z(x,y,z;s)   =η^x_zxη^y_yzτ_0.





 §.§ Solving for G_c

Using the IGG 𝖹_2 gauge symmetry, we had eliminated the phases on the RHS of Eq.<ref>, Eq.<ref>. To solve for G_c, one then plug the canonical expressions of the translational PSG elements into Eq.<ref>, Eq.<ref> and Eq.<ref>. One arrives at the following expressions:


    G^†_c(T^-1_y(i))η^-x_xyG_c(i)   =τ_0,
    
        G^†_c(T^-1_z(i))η^-x_zxη^-y_yzG_c(i)G_y(g^-1_c(i))   =τ_0,
    
        G^†_c(T^-1_x(i))G_c(i)G_z(g^-1_c(i))   =η_cxz.


Further simplifying the expressions, one arrives at:


    G_c(x,y,z)   =η^x_xyG_c(x,y-1,z),
    
        G_c(x,y,z)   =η^x_zxη^y_yzη^-y_xyG_c(x,y,z-1),
    
        G_c(x,y,z)   =η^-y_zxη^-z_yzη_cxzG_c(x-1,y,z).


The above expressions are valid for all sub-lattice indices, and we have suppressed the s indices. One then assumes that the following form is valid for G_c: G_c≡ f_c(x,y,z;s)𝔐_c(s). Because of the mentioned reason, we have f_c(x,y,z;s)=f_c(x,y,z). Then the separation of variables allows one to arrive at:


    f_c(x,y,z)   =η^x_xyf_c(x,y-1,z),
    
        f_c(x,y,z)   =η^x_zxη^y_yzη^-y_xyf_c(x,y,z-1),
    
        f_c(x,y,z)   =η^-y_zxη^-z_yzη_cxzf_c(x-1,y,z).


For f_c to be a path-independent function, there are certain constraints that the phases have to satisfy. For example, one considers two paths to arrive at f_c(x+1, y+1, z): 1.) f_c(x, y, z)↦ f_c(x+1, y, z)↦ f_c(x+1, y+1, z); 2.) f_c(x, y, z)↦ f_c(x, y+1, z)↦ f_c(x+1, y+1, z). One then compares the phases resulting from the two paths, and enforces them to be identical. Such a process produces the relevant constraints on the phases. We check the path independence on the xy, yz and zx planes respectively, and arrive at the following constraint:

    η_xy=η^-1_zx;η_yz=η_xy; η_zx=η^-1_yz.

It follows then η_xy=η_yz=η_zx=η_1. The previous equations on f_c become:


    f_c(x,y,z)   =η^x_1f_c(x,y-1,z),
    
        f_c(x,y,z)   =η^x_1f_c(x,y,z-1),
    
        f_c(x,y,z)   =η^-(y+z)_1η_cxzf_c(x-1,y,z).


Therefore, at this point we claim that G_c=η^xy+xz_1η^x_cxz𝔐_c(s).Let us take a look at Eq.<ref>. One can eliminate the phase on the RHS by making use of the IGG gauge symmetry. Plugging the above expression into Eq.<ref>, we arrive at:


    𝔐^3_c(α)   =τ_0; 
    𝔐_c(δ)𝔐_c(γ)𝔐_c(β)   =τ_0; 
    η_cxz    =1.


It is useful to make a summary before we close this subsection:

1.) η_c=η_cyx=η_czy=η_cxz=1;

2.) η_xy=η_yz=η_zx=η_1;
 3.)G_c=η^xy+xz_1𝔐_c(s), for which the following relations are satisfied:

    𝔐^3_c(α)   =τ_0; 
    𝔐_c(δ)𝔐_c(γ)𝔐_c(β)   =τ_0;




 §.§ Solving for G_a

To solve for G_a, one plugs the simplified expressions of the translational PSG elements into Eq.<ref>, Eq.<ref> and Eq.<ref>. One arrives at the following expressions:


    G^†_a(T_x(i))G_a(i)   =η_ax,
    
        G^†_a(T_y(i))G_y(T_y(i))G_a(i)G_y(g^-1_a(i))   =η_ay,
    
        G^†_a(T^-1_z(i))G^†_z(i)G_a(i)G_z(g^-1_a(i))   =η_az.


One makes the usual ansatz G_a(i)≡ f_a(x,y,z;s)𝔐_a(s), only this time one does not have f_a(x,y,z;s)=f_a(x,y,z), for the evaluation of G_y/z(g^-1_a(i)) is not s-independent.

We have, for s=α/δ, the following conditions for f_a:


    η^-1_axf_a(x,y,z;α/δ)   =f_a(x+1,y,z;α/δ),
    η^-1_ayf_a(x,y,z;α/δ)   =f_a(x,y+1,z;α/δ),
    η_azη_1f_a(x,y,z;α/δ)   =f_a(x,y,z+1;α/δ);


and for s=β/γ, the following conditions for f_a:


    η^-1_axf_a(x,y,z;β/γ)   =f_a(x+1,y,z;β/γ),
    η^-1_ayη^-1_1f_a(x,y,z;β/γ)   =f_a(x,y+1,z;β/γ),
    η_azf_a(x,y,z;β/γ)   =f_a(x,y,z+1;β/γ).


Note that this time we do not have to check the path independence of f_a, as the phases appearing in the above equations are constants. We then arrive at the following expressions:

    f_a(x,y,z;α/δ)=η^-x_axη^-y_ayη^z_azη^z_1,
       f_a(x,y,z;β/γ)=η^-x_axη^-y_ayη^-y_1η^z_az,

from which we write:

    G_a(x,y,z;α/δ)=η^-x_axη^-y_ayη^z_azη^z_1𝔐_a(α/δ),
       G_a(x,y,z;β/γ)=η^-x_axη^-y_ayη^-y_1η^z_az𝔐_a(β/γ).

In plugging these expressions into Eq.<ref>, we first consider i≡ (x,y,z;α). The condition we arrive at is:

    G^†_z(-x, -y-1, z;δ)G_a(-x, -y-1, z;δ)G_a(x,y,z;α)=η_a,

which further simplifies to:

    η^x+y+1_1η_ay𝔐_a(δ)𝔐_a(α)=η_a.

The above equation dictates that η_1=1. Consequently, one has:

    𝔐_a(δ)𝔐_a(α)=η_aη^-1_ay.

We then consider other sublattice sites, and they give us:


    G_a(-x-1, -y-1,z+1;γ)G_a(x,y,z;β)   =η_a,
    
    G_a(-x-1, -y-1,z;β)G_a(x,y,z;γ)   =η_a,
    
    G_a(-x, -y-1,z+1;α)G_a(x,y,z;δ)   =η_a.


Plugging the explicit forms into the above equations, and we arrive at:


    𝔐_a(γ)𝔐_a(β)   =η_aη^-1_axη^-1_ayη^-1_az,
    𝔐_a(β)𝔐_a(γ)   =η_aη^-1_axη^-1_ay,
    𝔐_a(α)𝔐_a(δ)   =η_aη^-1_ayη^-1_az.


It is useful to make a summary again before we close this subsection:

1.) η_xy=η_yz=η_zx=η_1=1, and as a result, G_x=G_y=G_z=τ_0;
 2.) We have G_a(x,y,z;s)=η^-x_axη^-y_ayη^z_az𝔐_a(s); for which the following relations are satisfied:


    𝔐_a(δ)𝔐_a(α)   =η_aη^-1_ay, 
    𝔐_a(γ)𝔐_a(β)   =η_aη^-1_axη^-1_ayη^-1_az,
    𝔐_a(β)𝔐_a(γ)   =η_aη^-1_axη^-1_ay,
    𝔐_a(α)𝔐_a(δ)   =η_aη^-1_ayη^-1_az.





 §.§ Solving for G_b

Now we attack Eq.<ref>, Eq.<ref> and Eq.<ref>. Since G_x, G_y and G_z are trivial now, the equations are reduced to the following form:


    G^†_b(T_x(i))G_b(i)   =η_bx,
    
        G^†_b(T^-1_y(i))G_b(i)   =η_by,
    
        G^†_b(T_z(i))G_b(i)   =η_bz.


We make the ansatz G_b≡ f_b(x,y,z)𝔐_b(s), where the separation of variables is possible because the above conditions are s-independent. One can quickly arrive at the condition that G_b=η^-x_bxη^y_byη^-z_bz𝔐_b(s).Plugging the expression into Eq.<ref>, one ends up with:

    G_b(g_b(i))G_b(i)=η_b.

Considering the individual sublattice sites respectively, one arrives at:


    𝔐_b(γ)𝔐_b(α)   =η_bη^-1_bx,
    𝔐_b(δ)𝔐_b(β)   =η_bη^-1_bxη^-1_bz,
    𝔐_b(α)𝔐_b(γ)   =η_bη^-1_bxη^-1_by,
    𝔐_b(β)𝔐_b(δ)   =η_bη^-1_bxη^-1_byη^-1_bz.


A quick summary:

We have G_b(x,y,z;s)=η^-x_bxη^y_byη^-z_bz𝔐_b(s); for which the following relations are satisfied:


    𝔐_b(γ)𝔐_b(α)   =η_bη^-1_bx,
    𝔐_b(δ)𝔐_b(β)   =η_bη^-1_bxη^-1_bz,
    𝔐_b(α)𝔐_b(γ)   =η_bη^-1_bxη^-1_by,
    𝔐_b(β)𝔐_b(δ)   =η_bη^-1_bxη^-1_byη^-1_bz.




 §.§ Solving Eq.<ref>, Eq.<ref> and Eq.<ref>

To attack the remaining three equations, note that we can use the IGG gauge symmetry of G_a and G_b to eliminate η_acb and η_cba. The equations are reduced to:


    G^†_a(g^-1_cg^-1_b(i))G^†_c(g^-1_b(i))G^†_b(i) 
       × G_a(T^-1_yT_x(i))G_c(g^-1_aT^-1_yT_x(i))=τ_0, 
       G^†_b(g^-1_a(i))G^†_a(i)G_b(T^-1_zT_yT^-1_x(i)) 
       × G_a(g^-1_bT^-1_zT_yT^-1_x(i))=η_abτ_0, 
       G^†_c(g^-1_bT^-1_xT_y(i))G^†_b(T^-1_xT_y(i))G_a(i)
       × G_b(g^-1_a(i))G_c(g^-1_bg^-1_a(i))G_b(g^-1_cg^-1_bg^-1_a(i))=τ_0.


We consider first i=(x,y,z;α). Plugging this into Eq.<ref> gives us the following:

    G^†_a(y-1,-z,-x-1;β)𝔐_c(γ)G^†_b(x,y,z;α)
       × G_a(x+1,y-1,z;α)𝔐_c(δ)=τ_0.

Recalling G_a=η^-x_axη^-y_ayη^z_az𝔐_a(s) and G_b=η^-x_bxη^y_byη^-z_bz𝔐_b(s), the LHS of the above expression is evaluated as:

    η^y-1_axη^-z_ayη^x+1_az𝔐^†_a(β)𝔐^†_c(γ)
       ×η^x_bxη^-y_byη^z_bz𝔐^†_b(α)η^-x-1_axη^-y+1_ayη^z_az𝔐_a(α)𝔐_c(δ)
    
    =   (η_azη_bxη^-1_ax)^x(η_axη^-1_ayη^-1_by)^y(η_azη_bzη^-1_ay)^zη_azη_ay
       ×𝔐^†_a(β)𝔐^†_c(γ)𝔐^†_b(α)𝔐_a(α)𝔐_c(δ).

Since the RHS of the previous expression is unit cell independent, we would have the following equations:

    η_azη_bx   =η_ax, 
    η_byη_ay   =η_ax, 
    η_azη_bz   =η_ay.

Naming η_2≡η_ax, η_3≡η_ay and η_4≡η_az, we would have η_bx=η_2η^-1_4, η_by=η_2η^-1_3 and η_bz=η_3η^-1_4. Also, we now have G_a=η^-x_2η^-y_3η^z_4𝔐_a(s) and G_b=η^-x+y_2η^-y-z_3η^x+z_4𝔐_b(s). The previous constraint becomes:

    𝔐^†_a(β)𝔐^†_c(γ)𝔐^†_b(α)𝔐_a(α)𝔐_c(δ)=η^-1_3η^-1_4.

What happens now for Eq.<ref>? Plugging i=(x,y,z;α) into the equation gives us:

    𝔐^†_c(γ)G^†_b(x-1,y+1,z;α) 
       × G_a(x,y,z;α)G_b(-x,-y-1,z-1;δ) 
       ×𝔐_c(β)G_b(-y-1,-z,x-1;δ)=τ_0.

Evaluating the LHS of the above equation gives us:

    η^x-y_2η^y+z+1_3η^-x-z+1_4η^-x_2η^-y_3η^-z_4η^x-y-1_2
       ×η^y-z_3η^-x+z-1_4η^y+1-z_2η^z-x+1_3η^-y+x_4
    ×   𝔐^†_c(γ)𝔐^†_b(α)𝔐_a(α)𝔐_b(δ)𝔐_c(β)𝔐_b(δ)
    
    =   (η^-1_2η^-1_3η_4)^x(η^-1_3η^-1_4η_2)^y(η_4η^-1_2η_3)^z
       ×𝔐^†_c(γ)𝔐^†_b(α)𝔐_a(α)𝔐_b(δ)𝔐_c(β)𝔐_b(δ).

Again, the unit cell independence gives us an extra condition η_2=η_3η_4. The original equation becomes:

    𝔐^†_c(γ)𝔐^†_b(α)𝔐_a(α)𝔐_b(δ)𝔐_c(β)𝔐_b(δ)=τ_0.

No further conditions on the phases can be derived from the three equations. We are in the position to write G_a=η^-x-y_3η^-x+z_4𝔐_a(s) and G_b=η^-x-z_3η^y+z_4𝔐_b(s). Iterating scenarios with different s for i≡ (x,y,z;s), we arrive at the following constraints:


    𝔐^†_a(β)𝔐^†_c(γ)𝔐^†_b(α)𝔐_a(α)𝔐_c(δ)=η^-1_3η^-1_4, 
       𝔐^†_a(γ)𝔐^†_c(δ)𝔐^†_b(β)𝔐_a(β)𝔐_c(γ)=η_4, 
       𝔐^†_a(α)𝔐^†_c(α)𝔐^†_b(γ)𝔐_a(γ)𝔐_c(β)=τ_0, 
       𝔐^†_a(δ)𝔐^†_c(β)𝔐^†_b(δ)𝔐_a(δ)𝔐_c(α)=η_3, 
       𝔐^†_b(δ)𝔐^†_a(α)𝔐_b(α)𝔐_a(γ)=η_abη_3η_4, 
       𝔐^†_b(γ)𝔐^†_a(β)𝔐_b(β)𝔐_a(δ)=η_abη_3η^-1_4, 
       𝔐^†_b(β)𝔐^†_a(γ)𝔐_b(γ)𝔐_a(α)=η_abη_3η^-1_4, 
       𝔐^†_b(α)𝔐^†_a(δ)𝔐_b(δ)𝔐_a(β)=η_abη_3η^-1_4, 
       𝔐^†_c(γ)𝔐^†_b(α)𝔐_a(α)𝔐_b(δ)𝔐_c(β)𝔐_b(δ)=τ_0,
       𝔐^†_c(δ)𝔐^†_b(β)𝔐_a(β)𝔐_b(γ)𝔐_c(α)𝔐_b(α)=η^-1_3,
       𝔐^†_c(α)𝔐^†_b(γ)𝔐_a(γ)𝔐_b(β)𝔐_c(δ)𝔐_b(γ)=η_3η_4,
       𝔐^†_c(β)𝔐^†_b(δ)𝔐_a(δ)𝔐_b(α)𝔐_c(γ)𝔐_b(β)=η_4.





 §.§ Solving for the 𝔐s

We start by performing some 𝖲𝖴(2) gauge transformations. Let us consider a type of gauge transformation W(i)≡ W(s). Had we started at a generic gauge, we could always make the following gauge transformation:

    W(α) = 𝔐_a(α), W(β)=𝔐_c(β), 
       W(γ)=𝔐^†_c(δ), W(δ) = τ_0.

so that:


    g.t.:      𝔐_c(δ)↦ W(β)𝔐_c(β)W^†(β)=τ_0, 
       𝔐_c(γ)↦ W(δ)𝔐_c(δ)W^†(δ)=τ_0, 
       𝔐_a(δ)↦ W(α)𝔐_a(α)W^†(α)=τ_0.


Now we make use of Eq.<ref>, and arrive at 𝔐_c(β)=𝔐_c(γ)=𝔐_c(δ)=τ_0. It should be noted that we are silent on 𝔐_c(α). Indeed, it is not possible to use gauge symmetry alone to trivialise 𝔐_c(α). However, as we shall see, other equations will bring enough restrictions on the form of 𝔐_c(α).Before we take a step further, let us note that taking traces over Eq.<ref> and Eq.<ref> dictates that η_4=1. After the simplification, we have:


    𝔐^3_c(α)   =τ_0, 
    𝔐_a(δ)𝔐_a(α)   =η_aη^-1_3,
    𝔐_a(γ)𝔐_a(β)   =η_a,
    𝔐_a(β)𝔐_a(γ)   =η_a,
    𝔐_a(α)𝔐_a(δ)   =η_aη^-1_3, 
    𝔐_b(γ)𝔐_b(α)   =η_bη^-1_3,
    𝔐_b(δ)𝔐_b(β)   =η_b,
    𝔐_b(α)𝔐_b(γ)   =η_bη^-1_3,
    𝔐_b(β)𝔐_b(δ)   =η_b,
    𝔐^†_a(β)𝔐^†_b(α)𝔐_a(α)   =η^-1_3, 
    𝔐^†_a(γ)𝔐^†_b(β)𝔐_a(β)   =τ_0, 
    𝔐^†_a(α)𝔐^†_c(α)𝔐^†_b(γ)𝔐_a(γ)   =τ_0, 
    𝔐^†_a(δ)𝔐^†_b(δ)𝔐_a(δ)𝔐_c(α)   =η_3, 
    𝔐^†_b(δ)𝔐^†_a(α)𝔐_b(α)𝔐_a(γ)   =η_abη_3, 
    𝔐^†_b(γ)𝔐^†_a(β)𝔐_b(β)𝔐_a(δ)   =η_abη_3, 
    𝔐^†_b(β)𝔐^†_a(γ)𝔐_b(γ)𝔐_a(α)   =η_abη_3, 
    𝔐^†_b(α)𝔐^†_a(δ)𝔐_b(δ)𝔐_a(β)   =η_abη_3, 
    𝔐^†_b(α)𝔐_a(α)𝔐_b(δ)𝔐_b(δ)   =τ_0,
    𝔐^†_b(β)𝔐_a(β)𝔐_b(γ)𝔐_c(α)𝔐_b(α)   =η^-1_3,
    𝔐^†_c(α)𝔐^†_b(γ)𝔐_a(γ)𝔐_b(β)𝔐_b(γ)   =η_3,
    𝔐^†_b(δ)𝔐_a(δ)𝔐_b(α)𝔐_b(β)   =τ_0.


Let us first look at Eq.<ref>, which can be rewritten as:

    𝔐_b(α)𝔐^†_b(β)𝔐_a(β)𝔐_b(γ)𝔐_c(α)=η^-1_3.

The above expression, when combined with Eq.<ref> and Eq.<ref>, gives us:

    𝔐_b(α)𝔐_b(γ)=η^-1_a,

that, when combined with Eq.<ref>, gives us:

    η^-1_a=η_bη^-1_3.

Now we look at Eq.<ref>, which can be rewritten as:

    𝔐_a(δ)𝔐^†_b(γ)𝔐^†_a(β)𝔐_b(β)=η_abη_3.

The above expression, when combined with Eq.<ref> and Eq.<ref>, gives us:

    𝔐_a(δ)𝔐_a(α)=η_a,

that, when combined with Eq.<ref>, gives us:

    η_3=1.

Looking at Eq.<ref> and Eq.<ref>, with the help of Eq.<ref>, we reach:

    𝔐_a(β)=η_abη_a𝔐_b(δ).

Similarly, Eq.<ref> and Eq.<ref>, with the help of Eq.<ref>, give us:

    𝔐_a(γ)=η_abη_a𝔐_b(β).

Consider, now, Eq.<ref> and Eq.<ref>. The coupled equations can be manoeuvred to give us:

    𝔐_a(α)𝔐_b(β)𝔐_b(γ)=τ_0.

The above equation, when paired with Eq.<ref>, gives us (note that η_a=η_b now):

    𝔐_b(β)𝔐_a(β)=η^-1_a,

which can be coupled with Eq.<ref> and Eq.<ref> to give us:

    η_ab=η_a.

At this point, there is only one phase left in the problem: η_a=η_b=η_ab. Let us look at Eq.<ref>, which gives us:

    𝔐^3_b(β)=η^-1_a,

which then implies that:

    𝔐^3_b(δ)=τ_0.

We note that, at this point, there are four independent 𝖲𝖴(2) matrices, denoted as follows:

    𝒜≡𝔐_c(α), ℬ≡𝔐_b(δ), 𝒞≡𝔐_a(α), 𝒟≡𝔐_b(α).

More completely, the 𝔐_as and 𝔐_bs are represented as follows:

    𝔐_a(α)   =𝒞,
    𝔐_a(β)   =ℬ,
    𝔐_a(γ)   =η_a ℬ^†,
    𝔐_a(δ)   = η_a 𝒞^†;

and 

    𝔐_b(α)   =𝒟,
    𝔐_b(β)   =η_a ℬ^†,
    𝔐_b(γ)   =η_a 𝒟^†,
    𝔐_b(δ)   =ℬ.

There are, in fact, only five independent constraints for these matrices:


    𝒜^3=τ_0,
       ℬ^3=τ_0,
       𝒜=𝒞ℬ𝒞^†,
       𝒜=𝒟ℬ𝒟^†,
       𝒜=𝒟ℬ^† 𝒞^†.


However, note that 𝒞=τ_0 from Eq. <ref>, we immediately have:

    𝒜 = ℬ,   𝒟=𝒜^†.

In summary:

1.) G_x=G_y=G_z=τ_0;

2.) G_a/b/c(x,y,z;s)=𝔐_a/b/c(s), where the 𝔐s have the following forms:

    𝔐_a(α)=τ_0, 𝔐_a(β)=𝒜, 𝔐_a(γ)=η_a 𝒜^†, 𝔐_a(δ)=η_a τ_0;
       𝔐_b(α)=𝒜^†, 𝔐_b(β)=η_a 𝒜^†, 𝔐_b(γ)=η_a 𝒜, 𝔐_b(δ)=𝒜;
       𝔐_c(α)=𝒜, 𝔐_c(β)=τ_0, 𝔐_c(γ)=τ_0, 𝔐_c(δ)=τ_0.

Since we have η_a=± 1, and 𝒜=τ_0, e^i2π/3τ_z, e^i4π/3τ_z, corresponding to apparently 6 solutions. However, let us note that a further gauge transformation W(x,y,z;s)≡η_a^(x+y+z) leads us to:

1.) G_x=G_y=G_z= η_aτ_0;

2.) G_a/b/c(x,y,z;s)=𝔐_a/b/c(s), where the 𝔐s have the following forms:

    𝔐_a(α)=τ_0, 𝔐_a(β)=𝒜, 𝔐_a(γ)=𝒜^†, 𝔐_a(δ)=τ_0;
       𝔐_b(α)=𝒜^†, 𝔐_b(β)=𝒜^†, 𝔐_b(γ)=𝒜, 𝔐_b(δ)=𝒜;
       𝔐_c(α)=𝒜, 𝔐_c(β)=τ_0, 𝔐_c(γ)=τ_0, 𝔐_c(δ)=τ_0.

Due to the fact that η_a now becomes global signs which are elements of the IGG, we conclude that η_a=± 1 PSGs are equivalent. We also remark that a gauge transformation W(x,y,z;s)≡ iτ_x maps the PSG solutions in which 𝒜=e^i2π/3τ_z to that in which 𝒜=e^i4π/3τ_z.In conclusion, we have:

1.) G_x=G_y=G_z= τ_0;

2.) G_a/b/c(x,y,z;s)=𝔐_a/b/c(s), where the 𝔐s have the following forms:

    𝔐_a(α)=τ_0, 𝔐_a(β)=𝒜, 𝔐_a(γ)=𝒜^†, 𝔐_a(δ)=τ_0;
       𝔐_b(α)=𝒜^†, 𝔐_b(β)=𝒜^†, 𝔐_b(γ)=𝒜, 𝔐_b(δ)=𝒜;
       𝔐_c(α)=𝒜, 𝔐_c(β)=τ_0, 𝔐_c(γ)=τ_0, 𝔐_c(δ)=τ_0,

where 𝒜=τ_0, e^i2π/3τ_z, corresponding to 2 solutions.



 §.§ Adding Time-Reversal Symmetry

Having arrived at the PSG solutions given the space group for the trillium lattice, we are at a position to add time-reversal symmetry (TRS) to the story. 
The extra relations are translated into the corresponding PSG equations:


    G_𝒯(i)G_𝒯(i)=η_𝒯, 
        G^†_𝒯(T^-1_x(i))G^†_x(i)G_𝒯(i)G_x(i)=η_x𝒯, 
        G^†_𝒯(T^-1_y(i))G^†_y(i)G_𝒯(i)G_y(i)=η_y𝒯, 
        G^†_𝒯(T^-1_z(i))G^†_z(i)G_𝒯(i)G_z(i)=η_z𝒯, 
        G^†_𝒯(g^-1_a(i))G^†_a(i)G_𝒯(i)G_a(i)=η_a𝒯, 
        G^†_𝒯(g^-1_b(i))G^†_b(i)G_𝒯(i)G_b(i)=η_b𝒯, 
        G^†_𝒯(g^-1_c(i))G^†_c(i)G_𝒯(i)G_c(i)=η_c𝒯.


We make the ansatz that G_𝒯≡ f(x,y,z; s)𝔐_𝒯(s), and Eq.<ref>, Eq.<ref> and Eq.<ref> immediately tell us that:

    G_𝒯(i)=η^x_x𝒯η^y_y𝒯η^z_z𝒯𝔐_𝒯(s).

The above form, when plugged into Eq.<ref>, gives us the following constraint:

    𝔐^2_𝒯(s)=η_𝒯.

We now discuss about the consequences of Eq.<ref>, Eq.<ref> and Eq.<ref>.


  §.§.§ Solving Eq.<ref>

First, we consider i≡(x,y,z;β), since in this case G_c(i)=τ_0. It is straightforward to reach the following constraint on the phases:

    η_x𝒯=η_y𝒯=η_z𝒯≡η_5.

Also the following constraints for 𝔐_𝒯(s) arise if we iterate the sub-lattice indices:


    𝔐^†_𝒯(δ)𝔐_𝒯(β)=η_c𝒯, 
       𝔐^†_𝒯(β)𝔐_𝒯(γ)=η_c𝒯, 
       𝔐^†_𝒯(γ)𝔐_𝒯(δ)=η_c𝒯, 
       𝔐^†_𝒯(α)𝔐^†_c(α)𝔐_𝒯(α)𝔐_c(α)=η_c𝒯.


From Eq.<ref> and Eq.<ref>, we can reach 𝔐^†_𝒯(δ)𝔐_𝒯(γ)=τ_0. This statement, when coupled with Eq.<ref> gives us η_c𝒯=1.

A quick summary:

1.) η_x𝒯=η_y𝒯=η_z𝒯≡η_5, η_c𝒯=1;

2.) We have G_𝒯(x,y,z;s)=η^x+y+z_5𝔐_𝒯(s); for which the following relations are satisfied:

    𝔐^†_𝒯(α)𝔐^†_c(α)𝔐_𝒯(α)𝔐_c(α)=τ_0,
       𝔐_𝒯(β)=𝔐_𝒯(γ)=𝔐_𝒯(δ).



  §.§.§ Solving Eq.<ref> and Eq.<ref>

Iterating the sub-lattice indices, we arrive at the following constraints:


    𝔐^†_𝒯(δ)𝔐^†_a(α)𝔐_𝒯(α)𝔐_a(α)   =η_a𝒯,
    𝔐^†_𝒯(γ)𝔐^†_a(β)𝔐_𝒯(β)𝔐_a(β)   =η_a𝒯, 
    𝔐^†_𝒯(β)𝔐^†_a(γ)𝔐_𝒯(γ)𝔐_a(γ)   =η_a𝒯η_5,
    𝔐^†_𝒯(α)𝔐^†_a(δ)𝔐_𝒯(δ)𝔐_a(δ)   =η_a𝒯η_5,
    𝔐^†_𝒯(γ)𝔐^†_b(α)𝔐_𝒯(α)𝔐_b(α)   =η_b𝒯,
    𝔐^†_𝒯(δ)𝔐^†_b(β)𝔐_𝒯(β)𝔐_b(β)   =η_b𝒯η_5,
    𝔐^†_𝒯(α)𝔐^†_b(γ)𝔐_𝒯(γ)𝔐_b(γ)   =η_b𝒯η_5,
    𝔐^†_𝒯(β)𝔐^†_b(δ)𝔐_𝒯(δ)𝔐_b(δ)   =η_b𝒯.




  §.§.§ Collection of Constraints

We further specify that 𝔐_𝒯(α)≡ℰ and 𝔐_𝒯(β)=𝔐_𝒯(γ)=𝔐_𝒯(δ)≡ℱ. We note that Eq.<ref> and Eq.<ref> immediately imply that η_5=1 since 𝔐_a(γ)=η_a𝔐^†_a(β). Furthermore, comparing Eq.<ref> and Eq.<ref> gives us η_a𝒯=η_b𝒯≡η_6, since 𝔐_a(β)=𝔐_b(δ).In the end, we reach the following five independent constraints:


    ℰ^2=η_𝒯,
       ℱ^2=η_𝒯,
       ℰ^† 𝒜^† ℰ𝒜=τ_0,
       ℱ^† ℬ^† ℱℬ=η_6,
       𝒞^† ℰ^† 𝒞ℱ=η_6.


Since 𝒞=τ_0, and 𝒜=ℬ, we can determine that η_6=1 and ℰ=ℱ. Also, when η_𝒯=1, we have ℰ=ℱ=τ_0; when η_𝒯=-1, we have ℰ=ℱ=iτ_z. Since without TRS, we had 2 solutions, now we have 4 solutions, as collected in Tab. <ref>.






§ IGG=𝖴(1)
 
In this section, our target is to find the PSG solutions with 𝖴(1) IGG. The PSG relations listed in Section <ref> still hold, only with the signs on the RHS being replaced as η_𝔤≡exp[iϕ_𝔤]. To proceed, we mention a fact which is a blessing for us. In <cit.>, Wen proved that for PSG solutions with 𝖴(1) IGG, the Gs always have the following canonical forms:

    G_g(i)≡ (iτ_x)^n_ge^iθ_g(i)τ_z,    n_g=0  or 1,

where θ_g∈ [0, 2π). Another thing we would like to mention before moving on is that, in this section θ always stands for a function which depends on position i, whereas ϕ always stands for a constant phase.Let us first look at Eq.<ref>, which can be rewritten as:

    G_a(g_a(i))G_a(i)=G_z(g_a(i))e^iϕ_aτ_z.

The above equation already dictates that n_z=0. Why? This is straightforward to see if n_a=0. Now supposing n_a=1, we would have:

    LHS_<ref>=e^i(-θ_a(g_a(i))+θ_a(i))τ_z,

which also implies that n_z=0 on the RHS_<ref>. Similarly, due to Eq.<ref> and Eq.<ref>, we can conclude that n_y=n_c=0.There is a valuable lesson from the above operation. Given a PSG equation G_𝔤G^†_𝔥…=e^iϕ_𝔦τ_z, we demand that (n_𝔤-n_𝔥… )=0  mod 2. Using the above lesson, we see from Eq.<ref> that n_x=0. And Eq.<ref> tells us that n_b=0, whereas Eq.<ref> implies that n_a=0. This is remarkable, for n_x=n_y=n_z=n_a=n_b=n_c=0! What was for us originally a set of coupled 𝖲𝖴(2) matrix equations now reduces to a set of coupled 𝖴(1) matrix equations, which are equations of compact 𝖴(1) phases.


 §.§ Solving for the Translational Elements and the Simplification

Before we take a step further, let us rewrite the remaining PSG equations in terms of θs:


    θ_c(g^2_c(i))+θ_c(g_c(i))+θ_c(i)=ϕ_c, 
       -θ_z(g_a(i))+θ_a(g_a(i))+θ_a(i)=ϕ_a, 
       -θ_y(g_b(i))+θ_b(g_b(i))+θ_b(i)=ϕ_b, 
       -θ_a(T_x(i))+θ_x(T_x(i))+θ_a(i)+θ_x(g^-1_a(i))=ϕ_ax, 
       -θ_a(T_y(i))+θ_y(T_y(i))+θ_a(i)+θ_y(g^-1_a(i))=ϕ_ay, 
       -θ_a(T^-1_z(i))-θ_z(i)+θ_a(i)+θ_z(g^-1_a(i))=ϕ_az, 
       -θ_b(T_x(i))+θ_x(T_x(i))+θ_b(i)+θ_x(g^-1_b(i))=ϕ_bx, 
       -θ_b(T^-1_y(i))-θ_y(i)+θ_b(i)+θ_y(g^-1_b(i))=ϕ_by, 
       -θ_b(T_z(i))+θ_z(T_z(i))+θ_b(i)+θ_z(g^-1_b(i))=ϕ_bz, 
       -θ_c(T^-1_y(i))-θ_y(i)+θ_c(i)+θ_x(g^-1_c(i))=ϕ_cyx, 
       -θ_c(T^-1_z(i))-θ_z(i)+θ_c(i)+θ_y(g^-1_c(i))=ϕ_czy, 
       -θ_c(T^-1_x(i))-θ_x(i)+θ_c(i)+θ_z(g^-1_c(i))=ϕ_cxz, 
       -θ_a(g^-1_cg^-1_bT^-1_xT_yg_ag_c(i))-θ_c(g^-1_bT^-1_xT_yg_ag_c(i)) 
       -θ_b(T^-1_xT_yg_ag_c(i))-θ_x(T_yg_ag_c(i))
       +θ_y(T_yg_ag_c(i))+θ_a(g_ag_c(i)) +θ_c(g_c(i))=ϕ_acb,
       -θ_b(g^-1_aT_xT^-1_yT_zg_bg_a(i))-θ_a(T_xT^-1_yT_zg_bg_a(i))
       +θ_x(T_xT^-1_yT_zg_bg_a(i))-θ_y(T_zg_bg_a(i))
       + θ_z(T_zg_bg_a(i))+θ_b(g_bg_a(i))+θ_a(g_a(i))=ϕ_ab,
        -θ_c(g^-1_bT^-1_xT_yg_ag_bg_cg_b(i))-θ_b(T^-1_xT_yg_ag_bg_cg_b(i)) 
       -θ_x(T_yg_ag_bg_cg_b(i))+ θ_y(T_yg_ag_bg_cg_b(i))
       +θ_a(g_ag_bg_cg_b(i))+θ_b(g_bg_cg_b(i))
       +θ_c(g_cg_b(i))+θ_b(g_b(i))=ϕ_cba.


These θs and ϕs in the above equations are compact 𝖴(1) phase factors, and an equation θ = ϕ means θ = ϕ mod 2π. The θs are associated with the 𝖲𝖴(2) gauge symmetry like before, specifically the 𝖴(1) subgroup of 𝖲𝖴(2) transforms the θ in the following way:

    θ_U(i)↦θ_U(i)-φ(i)+φ(U^-1(i)),   φ∈ [0,2π);

note that since we do not want to spoil the choice of τ_z, we consider only the 𝖴(1) subgroup of 𝖲𝖴(2).Similar to the 𝖹_2 case, we eliminate the phases on the right hand side of Eq. <ref>, Eq. <ref>, Eq. <ref>, Eq. <ref> and Eq. <ref>.



 §.§ Solving for the translational Elements

Let us start by considering the equations that arise because of the commutation of translational generators. Canonically, this gives us the following expressions of G_x, G_y, G_z after a gauge fixing:

    G_x(x,y,z;s)=τ_0, G_y(x,y,z;s)=e^ixϕ_xyτ_z, 
       G_z(x,y,z;s)=e^i(xϕ_zx+yϕ_yz)τ_z.

In other words, we have the following representation:

    θ_x(x,y,z;s)=0, θ_y(x,y,z;s)=xϕ_xy, 
       θ_z(x,y,z;s)=xϕ_zx+yϕ_yz.




 §.§ Solving for θ_c

Using the IGG gauge symmetry, one can eliminate the phases on the RHS of Eq.<ref>, Eq.<ref> and Eq.<ref>, as each of θ_x, θ_y and θ_z appears only once in these equations. To solve for θ_c, one then plug the canonical expressions of the translational PSG elements into Eq.<ref>, Eq.<ref> and Eq.<ref>. We make an ansatz analogous to the one we made in the 𝖹_2 case: θ(i)≡ f(x,y,z;s)+𝔪(s). We realise that the equations under attention are valid for all sub-lattice indices, therefore we have f_c(x,y,z;s)≡ f_c(x,y,z), and:


    f_c(x,y,z)   =f_c(x,y-1,z)+xϕ_xy, 
    
    f_c(x,y,z)   =f_c(x,y,z-1) 
       +xϕ_zx+yϕ_yz-yϕ_xy, 
    
    f_c(x,y,z)   =f_c(x-1,y,z) 
       -yϕ_zx-zϕ_yz + ϕ_cxz.


Checking the path-independency of f_c, we arrive at the following constraint:

    ϕ_xy=ϕ_yz=-ϕ_zx≡ϕ_1.

Eventually we arrive at the conclusion that f_c(x,y,z)=(xy-xz)ϕ_1.Let us look at Eq.<ref>. We had eliminated the phase on the RHS by making use of the IGG gauge symmetry. Plugging the above expression into Eq.<ref>, we arrive at:


    3𝔪_c(α)   =0,
    𝔪_c(β)+𝔪_c(γ)+𝔪_c(δ)   =0, 
    ϕ_cxz    = 0.


Before moving on, we make a brief summary:

    θ_x(i)=0, θ_y(i)=xϕ_1, θ_z(i)=(y-x)ϕ_1, 
       θ_c(i)=(xy-xz)ϕ_1+𝔪_c(s).



 §.§ Solving for θ_a

To solve for θ_a, one plugs the simplified expressions of the translational PSG elements into Eq.<ref>, Eq.<ref> and Eq.<ref>. One arrives at the following expressions:


    θ_a(T_x(i))   =θ_a(i)-ϕ_ax,
    θ_a(T_y(i))   =θ_a(i)+θ_y(T_y(i)) 
       +θ_y(g^-1_a(i))-ϕ_ay,
    θ_a(i)   =θ_a(T^-1_z(i))+ϕ_az
       +θ_z(i)-θ_z(g^-1_a(i)).


One makes the usual ansatz θ_a≡ f_a(x,y,z;s)+𝔪_a(s), only this time one does not have f_a(x,y,z;s)=f_a(x,y,z), for the evaluation of θ_y/z(g^-1_a(i)) is not s-independent.

We have, for s=α/δ, the following conditions for f_a:


    f_a(x+1,y,z;α/δ)   =f_a(x,y,z;α/δ)-ϕ_ax,
    
    f_a(x,y+1,z;α/δ)   =f_a(x,y,z;α/δ)-ϕ_ay,
    
    f_a(x,y,z;α/δ)   =f_a(x,y,z-1;α/δ) 
       +ϕ_az+(2y-2x+1)ϕ_1;


and for s=β/γ, the following conditions for f_a:


    f_a(x+1,y,z;β/γ)   =f_a(x,y,z;β/γ)-ϕ_ax,
    
    f_a(x,y+1,z;β/γ)   =f_a(x,y,z;β/γ)-ϕ_ay-ϕ_1,
    
    f_a(x,y,z;β/γ)   =f_a(x,y,z-1;β/γ)+ϕ_az
       +(2y-2x)ϕ_1.


Checking the path-independency of f_a, we arrive at the following constraint:

    2ϕ_1=0⇒ϕ_1=0  or π.

After the path-independency is guaranteed, we arrive at the following expressions:


    f_a(x,y,z;α/δ)   =-xϕ_ax-yϕ_ay+z(ϕ_az+ϕ_1),
    
    f_a(x,y,z;β/γ)   =-xϕ_ax-y(ϕ_ay+ϕ_1)+zϕ_az.


Plugging the forms of θ_a≡ f_a(x,y,z;s)+𝔪_a(s) into Eq.<ref>, further constraints can be derived. Specifically, we iterate the sub-lattice index. Let us consider i≡ (x,y,z;α), we have:

    -θ_z(-x,-y-1,z;δ) + θ_a(-x,-y-1,z;δ) + θ_a(x,y,z;α) = ϕ_a,

which gives us:

    ϕ_a     = -(-y-1+x)ϕ_1 + (xϕ_ax + (y+1)ϕ_ay
       + z(ϕ_az+ϕ_1))+(-xϕ_ax - yϕ_ay
       + z(ϕ_az+ϕ_1)) + 𝔪_a(α) + 𝔪_a(δ) 
        = -(-y-1+x)ϕ_1 + ϕ_ay + 2zϕ_az
        + 𝔪_a(α) + 𝔪_a(δ).

The implication from the above equation is that:

    ϕ_1 = 0,    2ϕ_az=0;   𝔪_a(α) + 𝔪_a(δ) = ϕ_a - ϕ_ay.

For i≡ (x,y,z;β), we have:

    θ_a(-x-1,-y-1,z+1;γ) + θ_a(x,y,z;β) = ϕ_a,

which gives us:

    ϕ_a     = ((x+1)ϕ_ax + (y+1)ϕ_ay + (z+1)ϕ_az) 
       +(-xϕ_ax - yϕ_ay + zϕ_az) 
       + 𝔪_a(β) + 𝔪_a(γ) 
        = ϕ_ax + ϕ_ay + ϕ_az + 𝔪_a(β) + 𝔪_a(γ).

For i≡ (x,y,z;γ), we have:

    θ_a(-x-1,-y-1,z;β) + θ_a(x,y,z;γ) = ϕ_a,

which gives us:

    ϕ_a     = ((x+1)ϕ_ax + (y+1)ϕ_ay + zϕ_az) 
       +(-xϕ_ax - yϕ_ay + zϕ_az) 
       + 𝔪_a(β) + 𝔪_a(γ) 
        = ϕ_ax + ϕ_ay+ 𝔪_a(β) + 𝔪_a(γ).

Combining with the equation for i≡ (x,y,z;β), we have:

    ϕ_az=0;   𝔪_a(β) + 𝔪_a(γ) = ϕ_a - ϕ_ay - ϕ_ax.

Lastly, the case for i≡ (x,y,z;δ) does not give us new relations.

Before moving on, we make a brief summary:

    θ_x/y/z(i)=0, θ_a(i)= -xϕ_ax - yϕ_ay +𝔪_a(s).




 §.§ Solving for θ_b

To solve for θ_b, one plugs the simplified expressions of the translational PSG elements into Eq.<ref>, Eq.<ref> and Eq.<ref>. One arrives at the following expressions:

    θ_b(T_x(i))   =θ_b(i)-ϕ_bx,
    θ_b(i)   =θ_b(T^-1_y(i))+ϕ_by, 
    θ_b(T_z(i))   =θ_b(i)-ϕ_bz.

One makes the usual ansatz θ_b≡ f_b(x,y,z;s)+𝔪_b(s), we arrive at

    f_b(x,y,z;s)=-xϕ_bx + y ϕ_by - zϕ_bz.

Plugging the forms of θ_b≡ f_b(x,y,z;s)+𝔪_b(s) into Eq.<ref>, further constraints can be derived. Specifically, we iterate the sub-lattice index. Let us consider i≡ (x,y,z;α), we have:

    θ_b(-x-1,y,-z;γ) + θ_b(x,y,z;α) = ϕ_b,

which gives us:

    ϕ_b     =  ((x+1)ϕ_bx +yϕ_by + zϕ_bz) 
       + (-xϕ_bx + y ϕ_by - zϕ_bz) 
       + 𝔪_b(α) + 𝔪_b(γ) 
        = ϕ_bx + 2yϕ_by + 𝔪_b(α) + 𝔪_b(γ).

The implication from the above equation is that:

    2ϕ_by=0;   𝔪_b(α) + 𝔪_b(γ) = ϕ_b - ϕ_bx.

For i≡ (x,y,z;β), we have:

    θ_b(-x-1,y,-z-1;δ) + θ_b(x,y,z;β) = ϕ_b,

which gives us:

    ϕ_b     =((x+1)ϕ_bx +yϕ_by + (z+1)ϕ_bz) 
       + (-xϕ_bx + y ϕ_by - zϕ_bz) 
       + 𝔪_b(β) + 𝔪_b(δ) 
        = ϕ_bx + ϕ_bz + 𝔪_b(β) + 𝔪_b(δ).

For i≡ (x,y,z;γ), we have:

    θ_b(-x-1,y+1,-z;α) + θ_b(x,y,z;γ) = ϕ_b,

which gives us:

    ϕ_b     =((x+1)ϕ_bx +(y+1)ϕ_by + zϕ_bz) 
       + (-xϕ_bx + y ϕ_by - zϕ_bz) 
       + 𝔪_b(α) + 𝔪_b(γ) 
        = ϕ_bx + ϕ_by + 𝔪_b(α) + 𝔪_b(γ).

Combining with the equation for i≡ (x,y,z;α), we have:

    ϕ_by=0;   𝔪_b(α) + 𝔪_b(γ) = ϕ_b - ϕ_bx.

Lastly, the case for i≡ (x,y,z;δ) gives us the same relation from the case for i≡ (x,y,z;β), which is:

    𝔪_b(β) + 𝔪_b(δ) = ϕ_b - ϕ_bx - ϕ_bz.

Before moving on, we make a brief summary:

    θ_b(i)= -xϕ_bx- zϕ_bz +𝔪_b(s).




 §.§ Solving Eq. <ref>

The equation Eq. <ref> is then reduced to:

    -θ_a(g^-1_c g^-1_b T^-1_x T_y(i)) - θ_c(g^-1_b T^-1_x T_y(i)) 
       - θ_b(T^-1_x T_y(i)) + θ_a(i) + θ_c(g^-1_a(i)) = 0.

For i≡ (x,y,z;α), the above equation is:

    0   =-θ_a(y,-z,-x;β) - θ_c(-x,y,-z;γ) 
       - θ_b(x-1,y+1,z;α) 
       + θ_a(x,y,z;α) +θ_c(-x,-y-1, z-1;δ)
        = x(ϕ_bx - ϕ_ax)+ y(ϕ_ax - ϕ_ay)+ z(ϕ_bz - ϕ_ay)-ϕ_bx
       -𝔪_a(β)-𝔪_c(γ)-𝔪_b(α)+𝔪_a(α)+𝔪_c(δ).

We can conclude from the above equation that ϕ_2≡ϕ_ax = ϕ_ay = ϕ_bx =  ϕ_bz. Thus we have f_a(x,y,z) = -(x+y)ϕ_2 and f_b(x,y,z) = -(x+z)ϕ_2.

For i≡ (x,y,z;β), the above equation is:

    0   =-θ_a(y,-z-1,-x;γ) - θ_c(-x,y,-z-1;δ) 
       - θ_b(x-1,y+1,z;β) + θ_a(x,y,z;β) 
       +θ_c(-x-1,-y-1, z;γ)
        = -2ϕ_2 -𝔪_a(γ)-𝔪_c(δ)-𝔪_b(β)
       +𝔪_a(β)+𝔪_c(γ).

For i≡ (x,y,z;γ), the above equation is:

    0   =-θ_a(y+1,-z,-x;α) - θ_c(-x,y+1,-z;α) 
       - θ_b(x-1,y+1,z;γ) + θ_a(x,y,z;γ) 
       +θ_c(-x-1,-y-1, z-1;β)
        = -𝔪_a(α)-𝔪_c(α)-𝔪_b(γ)+𝔪_a(γ)+𝔪_c(β).

For i≡ (x,y,z;δ), the above equation is:

    0   =-θ_a(y+1,-z-1,-x;δ) - θ_c(-x,y+1,-z-1;β) 
       - θ_b(x-1,y+1,z;δ) + θ_a(x,y,z;δ) 
       +θ_c(-x,-y-1, z;α)
        = -ϕ_2-𝔪_a(δ)-𝔪_c(β)-𝔪_b(δ)
       +𝔪_a(δ)+𝔪_c(α).




 §.§ Solving Eq. <ref>

The equation is reduced to:

    -θ_b(g^-1_aT_xT^-1_yT_z(i))-θ_a(T_xT^-1_yT_z(i))
       +θ_b(i)+θ_a(g^-1_b(i))=ϕ_ab.

For i≡ (x,y,z;α), the above equation is:

    ϕ_ab   = -θ_b(-x-1,-y,z;δ) - θ_a(x+1,y-1,z+1;α)
       +θ_b(x,y,z;α)+θ_a(-x-1,y-1,-z;γ)
       =ϕ_2 -𝔪_b(δ)-𝔪_a(α)+𝔪_b(α)+𝔪_a(γ).

For i≡ (x,y,z;β), the above equation is:

    ϕ_ab   = -θ_b(-x-2,-y,z+1;γ) - θ_a(x+1,y-1,z+1;β)
       +θ_b(x,y,z;β)+θ_a(-x-1,y-1,-z-1;δ)
       =ϕ_2 -𝔪_b(γ)-𝔪_a(β)+𝔪_b(β)+𝔪_a(δ).

For i≡ (x,y,z;γ), the above equation is:

    ϕ_ab   = -θ_b(-x-2,-y,z;β) - θ_a(x+1,y-1,z+1;γ)
       +θ_b(x,y,z;γ)+θ_a(-x-1,y,-z;α)
       =-ϕ_2 -𝔪_b(β)-𝔪_a(γ)+𝔪_b(γ)+𝔪_a(α).

For i≡ (x,y,z;δ), the above equation is:

    ϕ_ab   = -θ_b(-x-1,-y,z+1;α) - θ_a(x+1,y-1,z+1;δ)
       +θ_b(x,y,z;δ)+θ_a(-x-1,y,-z-1;β)
       =ϕ_2 -𝔪_b(α)-𝔪_a(δ)+𝔪_b(δ)+𝔪_a(β).




 §.§ Solving Eq. <ref>

The equation Eq.<ref> is then reduced to:

    -θ_c(g^-1_bT^-1_xT_y(i))-θ_b(T^-1_xT_y(i))+θ_a(i) 
        +θ_b(g^-1_a(i))+θ_c(g^-1_bg^-1_a(i)) 
        +θ_b(g^-1_cg^-1_bg^-1_a(i)) =0.

For i≡ (x,y,z;α), the above equation is:

    0   =-θ_c(-x,y,-z;γ)-θ_b(x-1,y+1,z;α)
       +θ_a(x,y,z;α)+θ_b(-x,-y-1,z-1;δ)
       +θ_c(x-1,-y-1,-z;β)+θ_b(-y-1,-z,x-1;δ)
       = 2ϕ_2 -𝔪_c(γ)-𝔪_b(α)+𝔪_a(α)
       +𝔪_b(δ)+𝔪_c(β)+𝔪_b(δ).

For i≡ (x,y,z;β), the above equation is:

    0   =-θ_c(-x,y,-z-1;δ)-θ_b(x-1,y+1,z;β)
       +θ_a(x,y,z;β)+θ_b(-x-1,-y-1,z;γ)
       +θ_c(x,-y-1,-z;α)+θ_b(-y-1,-z,x;α)
       = ϕ_2 -𝔪_c(δ)-𝔪_b(β)+𝔪_a(β)
       +𝔪_b(γ)+𝔪_c(α)+𝔪_b(α).

For i≡ (x,y,z;γ), the above equation is:

    0   =-θ_c(-x,y+1,-z;α)-θ_b(x-1,y+1,z;γ)
       +θ_a(x,y,z;γ)+θ_b(-x-1,-y-1,z-1;β
       )+θ_c(x,-y-2,-z;δ)+θ_b(-y-2,-z,x;γ)
       = 3ϕ_2 -𝔪_c(α)-𝔪_b(γ)+𝔪_a(γ)
       +𝔪_b(β)+𝔪_c(δ)+𝔪_b(γ).

For i≡ (x,y,z;δ), the above equation is:

    0   =-θ_c(-x,y+1,-z-1;β)-θ_b(x-1,y+1,z;δ)
       θ_a(x,y,z;δ)+θ_b(-x,-y-1,z;α)
       +θ_c(x-1,-y-2,-z;γ)+θ_b(-y-2,-z,x-1;β)
       = 2ϕ_2 -𝔪_c(β)-𝔪_b(δ)+𝔪_a(δ)
       +𝔪_b(α)+𝔪_c(γ)+𝔪_b(β).



 §.§ Collected equations for 𝔪s

In this subsection, we summarize the coupled equations to solve for 𝔪s. Before doing so, we note that we can use the 𝖲𝖴(2) gauge symmetry to fix certain 𝔪s. Recall that the action of the gauge transformation is:

    g.t.:θ_U(i)↦ w(i) -  θ_U(i) + w(U^-1(i)).

We start off in a generic gauge where all 𝔪s are non-trivial. We first perform the gauge transformation w(β) = 𝔪_c(β) and w(γ) = -𝔪_c(δ). The consequence is that 𝔪_c(β) = 𝔪_c(δ) = 0. And because 𝔪_c(β)+ 𝔪_c(γ) + 𝔪_c(δ)=0 from one of our constraints, we have 𝔪_c(γ)=0. We then perform the gauge transformation w(α) = 𝔪_a(α), such that 𝔪_a(α)=0. And because 𝔪_a(α) + 𝔪_a(δ) = ϕ_a - ϕ_2 from one of our constraints, we have 𝔪_a(δ) = ϕ_a - ϕ_2.The remaining equations after the reductions are:


    3𝔪_c(α)=0, 
       𝔪_a(β) + 𝔪_a(γ) = ϕ_a - 2ϕ_2, 
       𝔪_b(α) + 𝔪_b(γ) = ϕ_b - ϕ_2, 
       𝔪_b(β) + 𝔪_b(δ) = ϕ_b - 2ϕ_2, 
       -𝔪_a(β)-𝔪_b(α) = ϕ_2 , 
       -𝔪_a(γ)-𝔪_b(β)+𝔪_a(β) = 2ϕ_2, 
       -𝔪_c(α)-𝔪_b(γ)+𝔪_a(γ) =0, 
       -𝔪_b(δ)+𝔪_c(α)  = ϕ_2, 
       -𝔪_b(δ)+𝔪_b(α)+𝔪_a(γ)  = ϕ_ab-ϕ_2, 
       -𝔪_b(γ)-𝔪_a(β)+𝔪_b(β) = ϕ_ab-ϕ_a, 
       -𝔪_b(β)-𝔪_a(γ)+𝔪_b(γ)  = ϕ_ab+ϕ_2, 
       -𝔪_b(α)+𝔪_b(δ)+𝔪_a(β)  = ϕ_ab+ϕ_a-2ϕ_2, 
       -𝔪_b(α)+𝔪_b(δ)+𝔪_b(δ)  = - 2ϕ_2, 
       -𝔪_b(β)+𝔪_a(β)+𝔪_b(γ)
         +𝔪_c(α)+𝔪_b(α) = - ϕ_2, 
       -𝔪_c(α)-𝔪_b(γ)+𝔪_a(γ)+𝔪_b(β)+𝔪_b(γ) = - 3ϕ_2, 
       -𝔪_b(δ)+𝔪_b(α)+𝔪_b(β) = -ϕ_a- ϕ_2.


We now set A≡𝔪_c(α), B≡𝔪_a(β), C≡𝔪_b(α) and D≡𝔪_b(β). Eq. <ref> tells us that -B-C=ϕ_2. Also, Eq. <ref> tells us that D = 2B - ϕ_a. Thus all the 𝔪s can be represented using A and B, as deduced from Eq. <ref> to Eq. <ref>:

    𝔪_a(α)    = 0,
    𝔪_a(β)    = B,
    𝔪_a(γ)    = ϕ_a-2ϕ_2 - B,
    𝔪_a(δ)    = ϕ_a - ϕ_2,
    𝔪_b(α)    = -ϕ_2-B,
    𝔪_b(β)    = 2B-ϕ_a,
    𝔪_b(γ)    = ϕ_b+B,
    𝔪_b(δ)    = ϕ_b -2ϕ_2+ϕ_a-2B,
    𝔪_c(α)    = A,
    𝔪_c(β)    = 0,
    𝔪_c(γ)    = 0,
    𝔪_c(δ)    = 0.

Eq. <ref> tells us that:

    A+2B = ϕ_a - ϕ_b -2ϕ_2,

whereas Eq. <ref> tells us that:

    A+2B = ϕ_a + ϕ_b -ϕ_2.

From which we can see that ϕ_2=-2ϕ_b. Now if we look at Eq. <ref>, we have ϕ_ab = -ϕ_b. In fact, Eq. <ref> to Eq. <ref> do not tell us more than this. Eq. <ref> tells us that 3B = 4ϕ_b+2ϕ_a, Eq. <ref> tells us that A-B = -ϕ_a-ϕ_b, Eq. <ref> tells us that A-B = -2ϕ_b and finally  Eq. <ref> tells us that 3B = ϕ_a + ϕ_b-2ϕ_2. The above relations allow us to assert that:

    ϕ_3≡ -ϕ_ab=ϕ_a=ϕ_b,   ϕ_2 = -2ϕ_3,    B = A+2ϕ_3,

and:

    𝔪_a(α)    = 0,
    𝔪_a(β)    = A + 2ϕ_3,
    𝔪_a(γ)    =3ϕ_3-A,
    𝔪_a(δ)    = 3ϕ_3,
    𝔪_b(α)    = -A,
    𝔪_b(β)    = -A+3ϕ_3,
    𝔪_b(γ)    = 3ϕ_3+A,
    𝔪_b(δ)    = 2ϕ_3+A,
    𝔪_c(α)    = A,
    𝔪_c(β)    = 0,
    𝔪_c(γ)    = 0,
    𝔪_c(δ)    = 0.

We also have:

    f_a(i) = 2(x+y)ϕ_3,    f_b(i) = 2(x+z)ϕ_3.

Similar to the 𝖹_2 case, we now consider a further gauge transformation:

    w(x,y,z;α)   =-2ϕ_3+ϕ_3(x+y+z),
    
        w(x,y,z;β/γ/δ)   =ϕ_3(x+y+z),

we see that:

1.) G_x=G_y=G_z= e^-iϕ_3τ_z;

2.) G_a/b/c(x,y,z;s)=e^i𝔪_a/b/c(s)τ_z, where the 𝔪s have the following forms:

    𝔪_a(α)=0, 𝔪_a(β)=A, 𝔪_a(γ)=-A, 𝔪_a(δ)=0;
       𝔪_b(α)=-A, 𝔪_b(β)=-A, 𝔪_b(γ)=A, 𝔪_b(δ)=A;
       𝔪_c(α)=A, 𝔪_c(β)=0, 𝔪_c(γ)=0, 𝔪_c(δ)=0.

Due to the fact that ϕ_3 now becomes global signs which are elements of the IGG, we conclude that ϕ_3 is redundant. We also remark that a gauge transformation W(x,y,z;s)≡ iτ_x maps the PSG solutions in which 𝒜=e^i2π/3τ_z to that in which 𝒜=e^i4π/3τ_z.In conclusion, we have:

1.) G_x=G_y=G_z= τ_0;

2.) G_a/b/c(x,y,z;s)=e^i𝔪_a/b/c(s)τ_z, where the 𝔪s have the following forms:

    𝔪_a(α)=0, 𝔪_a(β)=A, 𝔪_a(γ)=-A, 𝔪_a(δ)=0;
       𝔪_b(α)=-A, 𝔪_b(β)=-A, 𝔪_b(γ)=A, 𝔪_b(δ)=A;
       𝔪_c(α)=A, 𝔪_c(β)=0, 𝔪_c(γ)=0, 𝔪_c(δ)=0,

where 𝒜=τ_0, e^i2π/3τ_z.



 §.§ Adding Time-Reversal Symmetry

We firstly write the algebraic relations:


    G_𝒯(i)G_𝒯(i)=e^iϕ_𝒯τ_z, 
        G^†_𝒯(T^-1_x(i))G^†_x(i)G_𝒯(i)G_x(i)=e^iϕ_x𝒯τ_z, 
        G^†_𝒯(T^-1_y(i))G^†_y(i)G_𝒯(i)G_y(i)=e^iϕ_y𝒯τ_z, 
        G^†_𝒯(T^-1_z(i))G^†_z(i)G_𝒯(i)G_z(i)=e^iϕ_z𝒯τ_z, 
        G^†_𝒯(g^-1_a(i))G^†_a(i)G_𝒯(i)G_a(i)=e^iϕ_a𝒯τ_z, 
        G^†_𝒯(g^-1_b(i))G^†_b(i)G_𝒯(i)G_b(i)=e^iϕ_b𝒯τ_z, 
        G^†_𝒯(g^-1_c(i))G^†_c(i)G_𝒯(i)G_c(i)=e^iϕ_c𝒯τ_z.


As usual, the canonical form of G_𝒯(i) = (iτ_x)^n_𝒯e^iθ_𝒯(i)τ_z. When n_𝒯=0, we can show that G_𝒯 = iτ_z uniformly much like the case for 𝖹_2. Since the derivation is very similar to the 𝖹_2 case, it is not included here. This group of PSG solutions does not produce mean field 𝖴(1) spin liquids if we consider the constraint imposed by TRS. For the rest of the appendix, let us focus on the case when n_𝒯=1.Let us first look at Eq. <ref>, we straightforwardly conclude that ϕ_𝒯=π. We denote θ_𝒯≡ f_𝒯(x,y,z;s) + 𝔪_𝒯(s). Then Eq. <ref> to Eq. <ref> tell us that:

    f_𝒯(x,y,z;s)=xϕ_x𝒯 + yϕ_y𝒯 + zϕ_z𝒯.

We would like to plug the above results into Eq. <ref>. We arrive at the folllowing constraint:

    -θ_𝒯(g^-1_c(i)) + 2θ_c(i) + θ_𝒯(i) = ϕ_c𝒯.

For the case with i=(x,y,z;α), we have:

    ϕ_c𝒯   =-θ_𝒯(y,z,x;α) + 2A +θ_𝒯(x,y,z;α)
       =x(ϕ_x𝒯 - ϕ_y𝒯) + y(ϕ_y𝒯-ϕ_z𝒯) 
       + z(ϕ_z𝒯 - ϕ_x𝒯) + 2A,

we arrive at ϕ_4≡ϕ_x𝒯 = ϕ_y𝒯 = ϕ_z𝒯, and ϕ_c𝒯=-A, where we used 3A=0.

For the case with i=(x,y,z;β), we have:

    ϕ_c𝒯   =-θ_𝒯(y,z,x;δ) +θ_𝒯(x,y,z;β)
       = -𝔪_𝒯(δ) + 𝔪_𝒯(β).

For the case with i=(x,y,z;γ), we have:

    ϕ_c𝒯   =-θ_𝒯(y,z,x;β) +θ_𝒯(x,y,z;γ)
       = -𝔪_𝒯(β) + 𝔪_𝒯(γ).

For the case with i=(x,y,z;δ), we have:

    ϕ_c𝒯   =-θ_𝒯(y,z,x;γ) +θ_𝒯(x,y,z;δ)
       = -𝔪_𝒯(γ) + 𝔪_𝒯(δ).

Thus if we denote 𝔪_𝒯(β)≡ E, we have 𝔪_𝒯(γ)= E-A and 𝔪_𝒯(δ)= E+A.We now look at Eq. <ref>. Similar to the case before:

    -θ_𝒯(g^-1_a(i)) + 2θ_a(i) + θ_𝒯(i) = ϕ_a𝒯.

For the case with i=(x,y,z;α), we have:

    ϕ_a𝒯   =-θ_𝒯(-x,-y-1,z-1;δ) 
       +2𝔪_a(α)+θ_𝒯(x,y,z;α)
       = - (-x-y-1+z-1)ϕ_4 -𝔪_𝒯(δ) 
        + (x+y+z)ϕ_4 + 𝔪_𝒯(α) 
        = x(2ϕ_4) + y(2ϕ_4) 
       +2ϕ_4 + 𝔪_𝒯(α) - 𝔪_𝒯(δ),

from which we have:

    2ϕ_4 = 0,   𝔪_𝒯(α) - 𝔪_𝒯(δ) = ϕ_a𝒯.

For the case with i=(x,y,z;β), we have:

    ϕ_a𝒯   =-θ_𝒯(-x-1,-y-1,z;γ)  
       +2𝔪_a(β)+θ_𝒯(x,y,z;β)
        = 𝔪_𝒯(β) - 𝔪_𝒯(γ) + 2A
        = 0.

For the case with i=(x,y,z;γ), we have:

    ϕ_a𝒯   =-θ_𝒯(-x-1,-y-1,z-1;β) 
       +2𝔪_a(γ)+θ_𝒯(x,y,z;γ)
        = 3ϕ_4 + 𝔪_𝒯(γ) - 𝔪_𝒯(β) -2A  
        = 3ϕ_4.

Note that this relation combined with the one before, gives us that ϕ_4=0.

For the case with i=(x,y,z;δ), we have:

    ϕ_a𝒯   =-θ_𝒯(-x,-y-1,z;α) 
       +2𝔪_a(δ)+θ_𝒯(x,y,z;δ)
        = ϕ_4 + 𝔪_𝒯(δ) - 𝔪_𝒯(α).

Combining the above constraints, we arrive at a set of relations summarised here:

    ϕ_a𝒯 = ϕ_4 = 0
         𝔪_𝒯(α) = 𝔪_𝒯(δ),  𝔪_𝒯(β) = 𝔪_𝒯(γ)+A.

Let us now look at Eq. <ref>:

    -θ_𝒯(g^-1_b(i)) + 2θ_b(i) + θ_𝒯(i) = ϕ_b𝒯.

For the case with i=(x,y,z;α), we have:

    ϕ_b𝒯   =-θ_𝒯(-x-1,y-1,-z;γ)  
       +2𝔪_b(α)+θ_𝒯(x,y,z;α)
        = - 2A + 𝔪_𝒯(α) - 𝔪_𝒯(γ).

For the case with i=(x,y,z;β), we have:

    ϕ_b𝒯   =-θ_𝒯(-x-1,y-1,-z-1;δ) 
       +2𝔪_b(β)+θ_𝒯(x,y,z;β)
        = - 2A  + 𝔪_𝒯(β) - 𝔪_𝒯(δ).

For the case with i=(x,y,z;γ), we have:

    ϕ_b𝒯   =-θ_𝒯(-x-1,y,-z;α)  
       +2𝔪_b(γ)+θ_𝒯(x,y,z;γ)
        = 2A  + 𝔪_𝒯(γ) - 𝔪_𝒯(α).

For the case with i=(x,y,z;δ), we have:

    ϕ_b𝒯   =-θ_𝒯(-x-1,y,-z-1;β) 
       +2𝔪_b(δ)+θ_𝒯(x,y,z;δ)
        = 2A  + 𝔪_𝒯(δ) - 𝔪_𝒯(β).

Combining the above constraints, we arrive at ϕ_b𝒯 = 0 and no new relations.We can then summarise:

    𝔪_𝒯(α)    = E+A ,
    𝔪_𝒯(β)    = E,
    𝔪_𝒯(γ)    = E-A,
    𝔪_𝒯(δ)    = E+A .

In the above, E is a free 𝖴(1) phase. However, note that we did not make use of the IGG gauge degrees of freedom associated with TRS. Recalling that G_𝒯∼ G_𝒯W_𝒯, where W_𝒯∈𝖴(1). We choose W_𝒯≡exp (-iE), thus eliminating the free phase in our solutions. We collect the 𝖴(1) PSG solutions into Tab. <ref>. Thus we have:

    G_𝒯(r⃗,α) = iτ_x e^iAτ_z,
       G_𝒯(r⃗,β) = iτ_x,
       G_𝒯(r⃗,γ) = iτ_x e^-iAτ_z,
       G_𝒯(r⃗,δ) = iτ_x e^iAτ_z.








§ MEAN-FIELD ANSATZES FOR THE PSG SOLUTIONS


Our PSG classification  obtains a set of gauge-inequivalent transformations G_g for all g ∈𝖯2_13 ×𝖹_2. In this appendix, we 
derive the constraints imposed on the mean-field parameters U_ij and μ_i by requiring that an element of the PSG leaves the ansatz invariant. We repeat
this condition for convenience:

    ∀ g:     G_g(g(i)) U_g(i)g(j)G^†_g(g(j)) = U_ij,
        G_g(g(i)) μ_g(i)G^†_g(g(i)) = μ_i.
        <ref>



 §.§ 𝖹_2



Here we specify the ansatzes for the PSGs corresponding to the IGG being 𝖹_2.
As derived in Appendix <ref> and displayed in Table <ref> in the main text, the four 𝖹_2 PSGs can be
indexed by , 𝒜=exp(i 2 π n/3) for i=0,1, and ℰ=τ_0 or ℰ=i τ_z, in terms 
of which all gauge transformations are listed in Tab. <ref>.

We note that if G_𝒯= τ_0  then the invariance of the ansatz under time-reversal requires U_ij=-U_ij and μ_i = -μ_i, 
leading to no non-zero mean-field ansatzes.

For G_𝒯=i τ_z, the invariance under TRS requires the following
form for all links and sites:

    U_ij= U^x_ijτ_x + U^y_ijτ_y 
    μ_i = μ^x_iτ_x + μ^y_iτ_y.


Finally, since G_x=G_y=G_z=τ_0 in our solutions, we must require
U_ij and μ_ij to be translationally invariant. We then encode
the dependence of the parameters U_ij on the link (ij) by 
determining U for each of the unique links in Tab. <ref>, and determining the functions U^x_i, U^y_i, where i∈{1,2 ⋯ 12} specifies the link in Table <ref>. The on-site parameters
are described as μ_α, μ_β, μ_γ and μ_δ, where the subscripts denote the sublattice dependence.

Imposing the invariance of the ansatz under the action of (G_c, g_c), we get the following relations between 4 groups of links that are closed under the 
application of g_c

    U_1    = U_5             = U_9 
      U_2    = U_6             = U_10
      U_3    = 𝒜 U_7 = 𝒜^2 U_11
    
                U_4    = 𝒜 U_8 = 𝒜^2 U_12


The relations between different groups of links are obtained by the invariance under (G_b,g_b) and (G_a,g_a). The action of (G_a,g_a) gives us

    U_1    = U_2      
      U_3    = U_4        
      U_5    = 𝒜^2 U_7  
      U_6    = 𝒜^2 U_8  
      U_9    = 𝒜 U_12
    
     U_10   =  𝒜 U_11

Similarly, the invariance of all links under (G_b,g_b) give us the conditions 

    U_1    =  U_3      
      U_2    =  U_4        
      U_5    = 𝒜^2 U_8  
      U_6    = 𝒜^2 U_7  
      U_9    = U_10
    
               U_11   = U_12


Combining the conditions in Eqs. <ref>,  <ref> and  <ref> we find that the U_ij for all links 
can be specified in terms of only two parameters U^x and U^y:

    U_1 = U^x τ_x + U^y τ_y;   
       U_2 =U_1;      U_3 = U_1 ;      U_4 =U_1;
       U_5 = U_1;      U_6 = U_1;      U_7 = 𝒜^2 U_1;
       U_8 = 𝒜^2 U_1;      U_9 = U_1;       U_10= U_1; 
       U_11 = 𝒜^2 U_1 ;           U_12 =𝒜^2 U_1


Similarly, demanding the invariance of μ_i under (G_c,g_c) gives us μ_γ=μ_δ=μ_β,
and μ_α=𝒜^2μ_α. Under (G_a, g_a), we have μ_α= μ_δ and
μ_β=𝒜^2 μ_γ. This already implies that when 𝒜≠ 1, 
μ =0 on all sites. When 𝒜=1, site-independent on-site terms of the form μ^x τ_x + μ^y τ_y are allowed 
in the ansatz. 



 §.§ PSG-protected gapless nodal star in 𝖹_21 spin liquid

In this Appendix, we prove that the mean-field Hamiltonian H_ MFT(k⃗) (Eqs. <ref> and Eqs. <ref>) for the 𝖹_21 QSL has 
two zero-energy eigenvalues for k⃗=(± k, ± k, ± k). To this end, we first work out the most general PSG-allowed H_ MFT(k⃗)  for the 
𝖹_21 QSL. The rest of the discussion assumes the translation invariance of the ansatzes, which is true for all our QSLs.
First, we use the basis (ψ^α_1, ψ^α_2,⋯ψ^δ_1 ψ^δ_2 ) to express 
the H_ MFT(k⃗) in terms of 2×2 blocks as 

    H_ MFT(k⃗)= 
      [  h_α,α(k⃗)  h_α,β(k⃗)  h_α,γ(k⃗)  h_α,γ(k⃗); h_β ,α(k⃗) h_β ,β(k⃗) h_β ,γ(k⃗) h_β ,γ(k⃗);  h_γ,α(k⃗)  h_γ,β(k⃗)  h_γ,γ(k⃗)  h_γ,γ(k⃗);  h_δ,α(k⃗)  h_δ,β(k⃗)  h_δ,γ(k⃗)  h_δ,γ(k⃗) ]

As just demonstrated in the previous section, when 𝒜≠ 1 and G_𝒯=iτ_z, we have 
h_α,α=h_β,β=h_γ,γ=h_δ,δ=0 for all k⃗. 
The 
block matrices have the  form  U_x τ^x +U_y τ^y (Eq. <ref>) in real space, leading to

    h_α,β(r⃗,r⃗')= 
      [               0   U_α,β(r⃗-r⃗'); U^*_α,β(r⃗-r⃗')               0 ]

for a complex amplitude U(r⃗).
The fourier-transformed equivalent is given by 

    h_α,β(k⃗)    =1/N∑_r⃗,r⃗' h_α,β(r⃗,r⃗')e^ik⃗.(r⃗-r⃗')
       =  [            0    U_α,β(k⃗); U^*_α,β(-k⃗)            0 ]

Foreseeing repeated appearances of the off-diagonal form in Eq. <ref>, we introduce the shorthand M[u(k)], defined by 

    M[u(k⃗)]
    =  [        0    u(k⃗); u^*(-k⃗)        0 ]

Also note that from Eq. <ref> we know that h_α,β(r⃗,r⃗') = h_β,α(r⃗',r⃗). So we have from Eq. <ref>

    h_α,β(k⃗) = h_α,β(-k⃗)




The action of symmetries on the block matrices in k-space can be worked from their real-space equivalents, given by
Eq. <ref>. To show this explicitly for a general symmetry transformation g, we assume the unit-cell independence of gauge transformations 
which has been shown for all our QSLs. To reduce cumbersome expressions, we   introduce the shorthand α̅ and g_α(r⃗)  to denote the sublattice index and unit-cell position of the operation  
  g(r⃗;α).
  We have, from Eq. <ref>,

    (G_g,g): ( h_α,β(k⃗) =1/N∑_r⃗,r⃗' h_α,β(r⃗,r⃗')e^ik⃗.(r⃗-r⃗')) 
       ↦1/N∑_r⃗,r⃗'G_g(α̅) h_α̅β̅(g_α(r⃗),g_β(r⃗'))G^†_g(β̅) e^ik⃗.(r⃗-r⃗')
       = 1/N∑_r⃗,r⃗'G_g(α̅) h_α̅β̅(r⃗,r⃗')G^†_g(β̅)e^ik⃗.(g^-1_α̅r⃗-g^-1_β̅(r⃗'))
       = 1/N∑_r⃗,r⃗'G_g(α̅) h_α̅β̅(r⃗,r⃗')G^†_g(β̅) e^ik⃗'.(r⃗-r⃗')+ϕ_g(α,β)
       = G_g(α̅)h_α̅,β̅(k⃗')G^†_g(β̅)e^iϕ_g(α,β)
       (G_g,g): h_α,β(k⃗↦ G_g(α̅)h_α̅,β̅(k⃗')G^†_g(β̅)e^iϕ_g(α,β)

From the third line to the fourth, we have used the fact that one can always write k⃗.(g^-1_α̅r⃗-g^-1_β̅(r⃗')) as
k⃗'.(r⃗-r⃗')+ϕ_g(α,β) for some k⃗' and a constant ϕ_g(α,β) independent of r⃗-r⃗'— this is always true for symmetry operations which are linearly represented on the lattice sites. 

Now, let us consider the symmetry transformation g=g_b . g_a acting on h_α,β. Using Eq. <ref> followed by Eq. <ref> , we find 

    h_α,β(k_x,k_y,k_z)     = h_β,α(k_x,-k_y,-k_z) e^i k_x 
        = h_α,β(-k_x,k_y,k_z) e^i k_x

Eq. <ref> can only be satisfied if the real space amplitudes u(r⃗,r⃗') in Eq. <ref> satisfy

    U_α,β(r⃗)= u(y,z)(δ_x,1+ δ_x,0),

where δ is the Kronecker delta not to be confused with the sublattice index, and u(y,z) is any complex function of the coordinates y and z. This form implies that,

    h_α,β = M[u(k_y,k_z)ζ(k_x)], 
    ζ(k_x)= (1+exp(i k_x))

The function u(k_x,k_y) is the fourier transform of u(y,z) defined in Eq. <ref>. All other block matrices in Eq. <ref> can be expressed in 
terms of u(k_x,k_y) by applying symmetry transformations to h_α,β.
Applying g_c to  h_α,β using Eq. <ref> gives us

    h_α,γ=𝒜^2 M[u(k_z,k_x)ζ(k_y)],
    
    h_α,δ=𝒜 M[u(k_x,k_y)ζ(k_z)], 
    where 𝒜 =exp(i (2π/3)τ^z).

We note that 

    𝒜^n M[u]=M[ω^n u], 
    where ω = 2 π/3.

This gives us

    h_α,γ= M[ω^2 u(k_z,k_x)ζ(k_y)],
    
    h_α,δ= M[ω u(k_x,k_y)ζ(k_z)], 
    where ω =exp(i 2π/3).

Applying g_b to h_γ,δ gives us


    h_γ,δ =M[ω u(k_y,-k_z) ζ(-k_x)exp(i k_x)]

Finally, applying g_c and g^2_c to h_γ,δ gives us 

    h_β,γ =M[ω u(k_x,-k_y) ζ(-k_z)exp(i k_y)],
    
    h_δ,β =M[ω u(k_z,-k_x) ζ(-k_y)exp(i k_x)].


Eqs. <ref>, <ref>, <ref>, <ref>, along with Eq. <ref> specify the most general form 
of all block matrices appearing in H_ MFT(k⃗)  in Eq. <ref> that is allowed by projective symmetries of the 𝖹_21 state.

Now we express H_ MFT(k⃗) in the basis 
(ψ^α_1, ⋯ψ^δ_1, ψ^α_2, ⋯ψ^δ_2) to  have  

    H_ MFT(k⃗)   =
    [       0   h(k⃗); h^†(k⃗)       0;         ].

The matrix the most general  h(k⃗) allowed by the PSG is given by


    h(k⃗)    = 
    (
    [                             0             ζ(k_x) u(k_y,k_z)        ω ^2 ζ(k_y) u(k_z,k_x)         ωζ(k_z)
       u(k_x,k_y);          ζ(-k_x) u(-k_y,-k_z)                             0 ω  e^i k_yζ(-k_z) u(k_x,-k_y) ω  e^-i k_xζ(k_y) u(-k_z,k_x);     ω ^2 ζ(-k_y) u(-k_z,-k_x) ω  e^-i k_yζ(k_z) u(-k_x,k_y)                             0 ω  e^i k_zζ(-k_x) u(k_y,-k_z);         ωζ(-k_z) u(-k_x,-k_y) ω  e^i k_xζ(-k_y) u(k_z,-k_x) ω  e^-i k_zζ(k_x) u(-k_y,k_z)                             0;                               ])



Now, we proceed to show that h(k⃗) has a maximum rank of 3 on the points (± k, ± k, ± k).
First consider k⃗=(k,k,k). With the further shorthands

    u = u(k,k),       v = u(-k,-k),       ρ = u(k,-k)
       λ = u(-k,k),       ζ = ζ(k),        r = e^i k


we have

    h(k,k,k)=
    (
    [           0        ζ  u   ζ  u ω ^2      ζ  u ω;      ζ ^* v           0         ζρω     ζ ^* λω; ζ ^* v ω ^2     ζ ^* λω           0         ζρω;    ζ ^* v ω         ζρω     ζ ^* λω           0;             ])

h(k,k,k) can now be brought to the row-echelon form by the sequence of elementary row-transformations given by

    R_2 ↔ R_1,   R_3 →ω R_3 -R_1,  R_4 →ω^2 R_4  -R_1 
       R_3 → uωζ R_3   - λζ^* R_2 ,  R_4 → R_4 u - R_2 ρ
       
     R_4    = -u ω(ζρ  (ω +1)-ζ ^* λω)R_4  
       +uζω ^2 (ζ ^* λ +ζρ) R_3   .

in conjunction with the using the identities ω^3=1, ζ r* =ζ^* and r r^*=1.The result of the row transformations is 

    (
    [                   ζ ^* v                        0              ζ ^* ρ  r ω                ζλ  r^* ω;                        0                     ζ  u                ζ  u ω ^2                   ζ  u ω;                        0                        0 ζ  u (ω +1) (ζ ^* λ +ζρ)        ζ  u (ζ ^* λ +ζρ);                        0                        0                        0                        0;                          ])

The appearance of the 0 in the last diagonal element establishes that the maximum rank of h(k,k,k) can be 3.

The analysis need not be repeated for the H_MFT at other gapless points 
like (-k,k,k), (k,-k,k) etc.  All of h(± k, ± k, ± k) are related to 
h(k,k,k) by elementary “rank-preserving" row and column transformations. h(-k,k,k) 
can be obtained from h(k,k,k) by the transformations 

    R_1 ↔ R_2,  R_3 ↔ R_4,
     C_1 ↔ C_2,  C_3 ↔ C_4, 
    
    R_4 → r^* C_4, C_4 → r C_4,

followed by two  re-identifications: u ↔ v, which have been considered independent complex numbers in the proof; and ω↔ω^2 which survive the important properties ω^3=1 and 1+ω+ω^2=0. h(k,-k,k) can be obtained in turn from h(-k,k,k) by the transformations 

    R_2 → R_3, R_3 → R_4, R_4 → R_2, 
    
    C_2 → C_3, C_3 → C_4, C_4 → C_2, 
    
    R_1 →ω^2 R_1, C_1 →ω^2  C_1.


h(k,k,-k) can be obtained from h(-k,k,k) by the transformations 

    R_2 → R_4, R_4 → R_3, R_3 → R_2, 
    
    C_2 → C_4, C_4 → C_3, C_3 → C_2, 
    
    R_1 →ω R_1, C_1 →ω  C_1.



This completes the proof that h(± k, ± k, ± k) has a maximum rank of 3,
and consequently, H_MFT(± k, ± k, ± k) has two gapless bands which are protected by projective symmetries against the addition of arbitrarily long-ranged terms
in the mean-field ansatz.

















 §.§ 𝖴(1)

The 𝖴(1) spin liquid mean field ansatz has the following form:

    U_ij = iU^0_ijτ_0 + U^z_ijτ_z,

dictated by the fact that the ansatz is invariant under the 𝖴(1) IGG gauge transformation. We would like to investigate how the PSG solutions we obtained constrain the nearest neighbor mean field ansatz by subjecting them to the following test:

    ∀ g∈𝖯2_13×𝖹^𝒯_2: Ĝ_g ĝ(U_ij) = U_ij.

Among the PSGs we have, we first study the class in which A=0. After enumerating all the conditions imposed by the PSG, we arrive at the following nearest neighbor mean field ansatz in the class where A = 0.:


    U_i = λτ_z,   where  i ∈{1,…,12}.

Next we would like to argue that, when A≠ 0, there would be no nearest neighbor mean field ansatz. We write U_1/3≡ i𝔘_1/3exp [iφ_1/3τ_z]. The TRS conditions on these two bonds give:

    2φ_1 + A = π,   φ_3 = π/2.

However, the condition Ĝ_b ĝ_b(U_1)=U_1 gives us:

    φ_3 + 2A + 4ϕ_3 = φ_1.

We time the above equation by 2, and combined with the last two relations, we would immediately arrive at 5A = 0. Note that we had 3A=0. Thus the ansatz does not vanish only when A=0.We conclude that we obtain only one nearest neighbor 𝖴(1) mean field ansatz:


    U_ζ   = λτ_z,   ζ∈{1,…,12 }
       a_s     = ωτ_z,    i∈{α,…,δ}.





