Chiral spin density wave, spin-charge-Chern liquid and d+id superconductivity in 1/4-doped correlated electronic systems on the honeycomb lattice

Recently two interesting candidate quantum phases --- the chiral spin density wave state featuring anomalous quantum Hall effect and the d+id superconductor --- were proposed for the Hubbard model on the honeycomb lattice at 1/4 doping. Using a combination of exact diagonalization, density matrix renormalization group, the variational Monte Carlo method and quantum field theories, we study the quantum phase diagrams of both the Hubbard model and t-J model on the honeycomb lattice at 1/4-doping. The main advantage of our approach is the use of symmetry quantum numbers of ground state wavefunctions on finite size systems (up to 32 sites) to sharply distinguish different quantum phases. Our results show that for $1\lesssim U/t<40$ in the Hubbard model and for $0.1<J/t<0.80(2)$ in the t-J model, the quantum ground state is either a chiral spin density wave state or a spin-charge-Chern liquid, but not a d+id superconductor. However, in the t-J model, upon increasing $J$ the system goes through a first-order phase transition at $J/t=0.80(2)$ into the d+id superconductor. Here the spin-charge-Chern liquid state is a new type of topologically ordered quantum phase with Abelian anyons and fractionalized excitations. Experimental signatures of these quantum phases, such as tunneling conductance, are calculated. These results are discussed in the context of 1/4-doped graphene systems and other correlated electronic materials on the honeycomb lattice.

The reliable determination of quantum phase diagrams of correlated electronic systems has been one of the central issues in quantum condensed matter physics. In the past decades, different analytic and numeric methods have been developed to attack this problem, including renormalization group methods 1-4 , quantum Monte Carlo methods (for a review see Ref. 5), variational Monte Carlo methods 5,6 , the density matrix renormalization group (DMRG) method 7,8 and the recently developed tensor-network methods [9][10][11][12][13][14] . Although each method has its advantages and disadvantages, this growing list of theoretical techniques has enabled careful investigations, and sometimes reliable determinations of quantum phase diagrams of correlated systems. In particular, in presence of strong correlations, a reliable understanding of the quantum phase diagrams of realistic model Hamiltonians usually strongly relies on unbiased numerical techniques. For instance, the DMRG method has successfully determined quantum phase diagrams of various quantum spin systems, and exotic quantum spin liquid phases were revealed [15][16][17][18][19] .
However, in the presence of doping, due to the larger dimension of Hilbert space and stronger quantum entanglement, a reliable determination of quantum phases remains challenging. The challenge is partially due to the fact that competing quantum phases cannot be sharply distinguished on finite-size systems in an obvious fashion, while most cutting-edge numerical simulations can only be performed on finite-size systems.
In this work, we show that a combination of different quantum many-body techniques allows, to a certain level, a sharp determination of quantum phases of correlated electronic systems at some commensurate dopings 20 . Particularly, we demonstrate our approach in the 1/4-doped correlated systems on the honeycomb lattice. Our approach is based on our ability of analytically writing down symmetric quantum wavefunctions of different candidate quantum phases on finite size systems, studying their characteristic quantum numbers and other properties, and comparing with results from unbiased numerical simulations such as exact diagonalization and DMRG. When different quantum phases can be analytically shown to have different lattice quantum numbers, this approach has the power to sharply distinguish them even on small lattices.
Recently, interesting quantum phases were proposed for the 1/4-doped Hubbard model on the honeycomb lattice. Considering the nearest neighbor single-band tightbinding model on the honeycomb lattice, both the 1/4 electron-doped and hole-doped systems feature a Fermi surface of hexagonal shape (see Fig.1 even in the presence of weak interactions. There are two important features of the hexagonal Fermi surface: The opposite sides of the Fermi surface are nested by three wavevectors Q 1,2,3 , and three Van Hove singularities are located at the mid-points of the Brillouin Zone boundary M 1,2,3 . Previous studies have revealed two interesting candidate quantum phases: The chiral spin density wave state 21 and the d+id superconductor 22,23 , both of which can be understood starting from the two features of the hexagonal Fermi surface.
It is well-known that nested Fermi surfaces can cause magnetism.
Based on Hartree-Fock meanfield calculations 21 and functional renormalization group calculations 24,25 , it has been shown that the three nested wavevectors together could give rise to a rather exotic type of magnetic ordering at intermediate coupling strengths: The tetrahedral magnetic order which quadruples the unit cell (see Fig.2a). Due to the non-coplanar magnetic ordering pattern, electrons pick up Berry's phase when hopping around the lattice, similarly to the effect of a non-uniform magnetic field. Consequently the electronic band structure is found to carry a non-zero Chern number. This magnetically ordered phase, termed the chiral spin density wave (c-SDW) state, is a topological phase featuring gapless electronic edge states and anomalous quantum Hall effect σ xy = e 2 /h. 21,24 In addition, analytical renormalization group calculations, focusing on scattering involving electronic states at the Van Hove singularities, show a d+id superconductor as the ground state of the system, which in principle could be a high temperature phenomenon 23 (see Fig.2b for the pairing order parameter in real space). The same d+id superconductor (SC) has also been proposed for the Hubbard model and t-J model on the honeycomb lattice over a large range of doping levels based on renormalization group studies 22-26 , variational Monte Carlo approaches 27 , and tensor-network numerical simulations 28 . The d+id SC phase turns out to be a topological superconductor hosting a spin-quantum-Hall effect. 29 Apart from these two phases, in this work we propose yet another candidate quantum phase, denoted as spincharge-Chern liquid (SCCL), which could also be realized in correlated electronic systems on the honeycomb lattice at 1/4 doping. SCCL can be viewed as the resulting phase after the long-range magnetic order in the c-SDW phase is quantum melted. In the past, quantum melting of long-range magnetic order was discussed in the context of undoped quantum spin systems, and the resulting exotic phases, quantum spin liquids, have attracted considerable interest (see, e.g., 15,16,[30][31][32][33][34][35]. It is known that strong quantum fluctuations are necessary to stabilize such liquid phases. In fact, most candidate quantum spin liquid materials are spin-1/2 systems, where quantum fluctuations are strong. Intuitively, quantum fluctuations of spin degrees of freedom are likely to be even stronger in doped spin-1/2 systems, which can be justified by slave-fermion mean-field arguments (see Sec.II). This suggests that liquid phases such as SCCL may have a better chance to be stabilized in doped correlated electronic systems.
Unlike c-SDW, the SCCL phase respects spin-rotation and lattice translation and rotation symmetries, while breaking the time-reversal symmetry; nevertheless, both the charge and spin excitations are gapped in the bulk. This violation of Luttinger's theorem is due to the fact that SCCL is a fractionalized phase with topological order. For example, in the bulk SCCL features charge-1/2, spin-neutral anyon excitations with θ = π/4 exchange statistics. On the boundary, SCCL hosts chiral gapless edge states of charge-1, spin-neutral fermions. We show that although the electromagnetic response in the bulk of the SCCL is described by an anomalous quantum Hall response similar to the c-SDW phase: j x = σ xy E y , where σ xy = e 2 /h, the SCCL and c-SDW have very different signatures in transport experiments, which can be used to identify them in candidate materials. One important result of the current work is that the conductance through a weakly coupled tunneling junction with a metallic lead (namely, G e 2 /h) in a SCCL phase should vanish as G(T ) ∝ T 4 at low temperatures, while in the c-SDW phase this should obey G(T ) ∝ constant.
Experimentally, single-band correlated electronic models on the honeycomb lattice are relevant for many materials. For instance, doped graphene may be described by the Hubbard model in the intermediate correlated regime U/t = 2 ∼ 3. 36 More candidate materials, including certain transition metal oxide heterstructures, will be discussed later in this paper. Although experimental realization of 1/4-doping on these materials has not yet been reported, with the fast developing material science techniques on thin film synthesis, this doping level may be achievable within foreseeable future. This motivates us to carefully investigate the phase diagrams of the corre-lated electronic systems on the honeycomb lattice at 1/4 doping, especially over the intermediate to strong correlation strengths. The previous studies are either based on mean-field theories 21 which is biased, or renormalization group techniques [22][23][24][25] which presumably are under control only for the weak coupling regime.
We study both the 1/4-doped Hubbard model and t-J model on the honeycomb lattice: Here P G is the usual Gutzwiller projection operator removing the double occupancies in the t-J model. Due to the particle-hole symmetry of these nearest neighbor (t only) models, our study applies for both the electrondoped and hole-doped systems. We analytically constructed quantum wavefunctions, performed exact diagonalization on the 8-site sample, DMRG simulations on the 24-site and 32-site samples, and variational Monte Carlo simulations. These allow us to construct reliably, at least to a certain extent, the quantum phase diagrams from intermediate to strong coupling regimes. Our main results are summarized in Fig.3. The limitation of our calculations are twofold. First, we cannot address the phase diagram reliably for the weak coupling regime in the Hubbard model: U/t 1, because correlation lengths of the competing phases can be much larger than the investigated system sizes. Second, we cannot sharply distinguish the c-SDW phase from the SCCL phase, because they are distinguished only by long-range physics which requires careful finite size scaling and larger system sizes.
Despite these limitations, we find that c-SDW or SCCL is stablized in the majority of the physically realistic parameter regime: 1 U/t < 40 in the Hubbard model and 0.1 < J/t < 0.80(2) in the t-J model. Between the c-SDW and SCCL phases, the measurements of correlation functions suggest that SCCL is more likely to be realized in the small U/t and large J/t regimes within these parameter windows. The d+id SC phase is found in the t-J model at 0.80(2) < J/t. The sharp distinction between the c-SDW/SCCL phase and the d+id SC phase becomes possible because they have different lattice quantum numbers on the 32-site sample.
The remaining parameter regimes will be briefly discussed but will not be the focus of the present paper, since it is unclear whether these regimes are relevant for correlated materials. For instance, in the t-J model with J/t < 0.1 and in the Hubbard model with U/t 50, we find some inconclusive evidence for a different homogeneous phase. This parameter regime is adjacent to the infinite-U Hubbard problem, 37,38 and will be left for future study. In addition, when J/t 3 in the t-J model, FIG. 3. The phase diagrams of the correlated electronic systems on the honeycomb lattice at 1/4 doping. (a) In the Hubbard model, the ground state is found to be either in a c-SDW phase or in a SCCL phase over the majority of the parameter range 1 U/t < 40. (b) In the t-J model, the c-SDW phase or the SCCL phase is identified in the regime 0.1 < J/t < 0.80(2), seperated from the d+id superconductor phase at larger J/t by a first order phase transition. evidence of charge inhomogeneity is observed, which is likely due to phase separation and consistent with physical intuition.
There is a useful by-product of our investigation. It was proposed that, on finite size lattices, the rotational symmetry eigenvalues in the ground state manifold of a topologically ordered phase can be determined by the modular transformation matrices. [39][40][41][42][43] However, the SCCL phase here serves as a counterexample of this claim, because the rotational eigenvalues are system-size dependent (see Section IV).
The paper is organized as follows. In Section II we analytically construct symmetric quantum wavefunctions of the three phases: c-SDW, SCCL and d+id SC. We show that they have characteristic signatures in quantum numbers on finite size lattices. In Section III we present the results from a combination of different numerical simulations, which convincingly justify the phase diagrams in Fig.3. Because SCCL is a new topologically ordered quantum phase, its low energy effective theory, fundamental properties, and experimental signatures are studied in Section IV. Finally, we discuss our methodology and results, in particular in the context of a few candidate materials in Section V.

II. SYMMETRIC WAVEFUNCTIONS OF COMPETING QUANTUM PHASES
To identify fingerprints of these candidate quantum phases in unbiased numerical simulations, we explicitly write down the symmetric quantum wavefunctions of the competing c-SDW, SCCL and d+id SC phases on finite-size lattices. Note that there is no sense of spontaneous symmetry breaking on finite-size lattices and these quantum wavefunctions are symmetric (i.e, forming irre-ducible representations) of the full symmetry group involving both the SU (2) spin-rotation group and the lattice space group. This is also why the c-SDW and SCCL cannot be sharply distinguished on finite-size lattices because they share the same quantum numbers.
We construct the c-SDW/SCCL quantum wavefunctions using the slave-fermion approach, 44-48 and construct the d+id SC wavefunctions by the slave-boson approach. 49,50 Note that these wavefunctions are constructed in the Hilbert space of the t-J model; however, their quantum numbers on finite size lattices are unchanged in the Hubbard model. This is because the small J/t regime in the t-J model and large U/t regime of the Hubbard model are smoothly connected, and if two quantum states have different quantum numbers on finite size lattices, they cannot represent the same quantum phase. In addition, although these quantum wavefunctions are not the exact ground states of simple model Hamiltonians, they have the same universal properties of the quantum phases that they belong to, including symmetry quantum numbers.

A. d+id SC
Here we briefly describe these symmetric quantum wavefunctions. The details can be found in Appendix B). A d+id SC can be constructed using the slaveboson approach, 49,50 in which the electrons are split into fermionic spinons and bosonic holons: This parton construction enlarges the Hilbert space and has a U (1) gauge redundancy. This gauge redundancy is broken by boson condensation | b i | = √ x at zero temperature (x is the doping fraction), which is required to accomodate the doped charge at the mean-field level. A d+id SC can be represented if f iα fermions form a d+id SC band structure, and the bosons are condensed at the Γ-point. The associated physical wavefunction is obtained after projecting out the unphysical states of the t-J model, i.e., it is a simple Gutzwiller projected d+id SC wavefunction: where P N is the projector into a fixed fermion number sector, enforcing that the total number of fermions equals 3/4 of the total number of sites. |Ψ M F d+id (χ, ∆) is the ground state of the d+id SC mean-field Hamiltonian: Here χ is the real hopping, while singlet pairing ∆ ij has the real space pattern shown in Fig.2b. Namely, a) b) Red: −e iφ . Here the real number φ = φ b for bosons and φ = φ f for fermions. φ b and φ f can be viewed as two variational parameters. The above pattern is for one of the two degenerate ground states while the other one is its time-reversal image, which can be obtained by sending these amplitudes to their complex conjugates: Sites numbered 1 to 8 label the quadrupled unit-cell used in Appendix C.
∆ ij = ∆, ∆e 2πi/3 , ∆e −2πi/3 depending on the orientations of the bonds (∆ is chosen to be real). For simplicity we include the nearest neighbor amplitudes only. The µ f is tuned to satisfy f † iα f iα = 3/4 at the mean-field level, so it is not a variational parameter. Note that after the sign of χ is fixed, only the ratio ∆/χ is a variational parameter of the constructed wavefunction |Ψ d+id (χ, ∆) . This single-parameter variational wavefunction will be used in the variational Monte Carlo study in Section III B, which reproduces ∼ 97 − 99% of the ground state energy in the d + id phase shown in Fig.3b.

B. c-SDW/SCCL
To construct c-SDW/SCCL wavefunctions, we use the slave-fermion approach, [44][45][46][47][48] in which electrons are split into bosonic spinons and fermionic spinless holons: which also enlarges the Hilbert space and has a U (1) gauge redundancy. At the mean-field level, the spin dynamics is described by a bosonic superconductor, and the charge dynamics is described by a spinless fermion band structure: where B ij and A ij are boson singlet hopping and pairing on bond ij, and χ ij is the spinless fermion hopping. Nonzero B ij and A ij , which are required to describe c-SDW/SCCL, break the U (1) gauge redundancy down to Z 2 . The mean-field boson(fermion) wavefunc- is the ground state of the corresponding Hamiltonian in Eq.6, which can be mathematically represented as a permanent(determinant). The associated physical wavefunction |Ψ c−SDW/SCCL is obtained by gluing two parts together and going back to the Hilbert space of the t-J model.
More precisely, note that any physical state in the t-J model can be expanded in the spin-occupation basis {|s 1 , s 2 , s 3 ...s N }, where N is the number of sites, and s i =↑, ↓, 0 depending on whether the site-i is spin-up, spin-down or empty: where a certain ordering of sites is required in the last product to take care of the fermion sign. The physical wavefunction |Ψ c−SDW/SCCL is defined as: i.e., |Ψ c−SDW/SCCL is a product of a permanent(the second line) and a determinant(the third line). It turns out that the real space pattern of A ij , B ij , χ ij as shown in Fig.4 is describing the c-SDW/SCCL phases (see Appendix B). For simplicity we plot these amplitudes only on the nearest neighbor (NN) and next nearest neighbor (NNN) bonds. This complicated pattern ensures that the wavefunction is symmetric under lattice space group while capturing the tetrahedral spin correlation.
One can see that the unit-cell of the amplitudes doubles the original unit-cell of the honeycomb lattice, which indicates that the mean-field states |Ψ M F b (|Ψ M F f ) break translational symmetry. However, the physical state |Ψ c−SDW/SCCL is fully translationally symmetric, as shown in Appendix B. Similar states having doubled unit-cell of the mean-field amplitudes are often called πflux states in the context of quantum spin liquids.
In addition, this doubling of unit-cell is physically important. This is why the spinless fermion filling f † i f i = 1/4, required by the 1/4-doping, actually corresponds to a fully filled lowest f -fermion band, which is separated from higher bands by an energy gap generated by the imaginary part of the NNN hopping e iφ f . Similarly to the Haldane model of spinless fermions, 51 which preserves the original unit-cell of the honeycomb lattice, the lowest energy band of f -fermion here is found to carry nonzero Chern-number C = 1. Because f -fermion describes the charge dynamics, the electromagnetic response of c-SDW/SCCL features an anomalous quantum Hall response, σ xy = e 2 /h. Now we describe the difference between the c-SDW phase and the SCCL phase in the above slave-fermion formulation. At the mean-field level, µ b is chosen so that b † iα b iα = 3/4 to be consistent with the doping level. On a finite size lattice, this is always achieved by tuning µ b so that the boson band minima are close enough to, but not touching, zero. Note that when the bosonic band minima touch zero, boson condensation occurs and longrange tetrahedral magnetic order is established (see Appendix B). This is the c-SDW phase in the slave-fermion formulation. However, because boson condensation never occurs on finite size lattice due to the presence of boson pairing, the difference between the two phases appears only in the thermodynamic limit (L → ∞). In this limit, if the boson band minima separate from zero by a finite gap, the resulting phase is a SCCL; however, if the gap closes the resulting phase is a c-SDW.
The SCCL phase is thus a fully gapped phase in the bulk, which will be studied in detail in Section IV. Nevertheless it is helpful to mention some of its basic properties here. Because the bosons do not condense, there is a remaining Z 2 gauge dynamics which dictates the existence of a topological order. However, the topological order in SCCL is fundamentally different from a usual Z 2 topological order such as the one in Kitaev's toric code. 52 In a usual Z 2 topological order, there are three types of nontrivial quasiparticles: bosonic Z 2 -gauge-charge e, bosonic vison (π-gauge-flux) m, and the fermionic bound state em. But in SCCL, the three nontrivial quasiparticles are: spin-1/2 and charge-neutral bosonic Z 2 -gaugecharge b α (which can be identified with the spinons), spin-neutral and charge-1/2 vison v with statistical angle θ = ±π/4, and their bound states: spin-1/2-charge-1/2 anyon b α v with statistical angle ±5π/4. Here the two signs of the statistical angles correspond to the two degenerate ground states which are time-reversal images of each other. The charge-1/2 vison v is simply due to the fact that fermion-f fills a Chern band. Thus the vison, a π-gauge flux, will be bound with 1/2-charge.
The charge-1, spin-neutral fermionic holon-f differs from a spinon only by an electron. Therefore mathematically it is not a new type of quasiparticle. However, there are gapless chiral edge states formed by f on the boundary, which is clear at the mean-field level. This means that although the spin-gap is opened everywhere in the SCCL phase, the charge gap is closed on the boundary. Because of the spin gap, single electron tunneling into the edge states is forbidden at low energy. However, singletpairs of electrons can still tunnel into the edge, which is the origin of the T 4 power-law tunneling conductance at low energy.
Finally, the slave-fermion formulation of long-range magnetic order allows us to argue that the spin liquid phases, such as SCCL, may be easier to be stablized in the doped systems compared with the undoped spin-1/2 systems. In the past, a great number of spin-1/2 models were investigated in a search for quantum spin liquids. Only few of these models can host spin liquid phases. [53][54][55][56] In the slave-fermion formulation (which in the undoped case, is the same as the Schwinger-boson formulation), this can be understood as follows.
For a given value of mean-field parameters A ij /B ij , b † iα b iα increases as µ b increases and the boson quasiparticle gap decreases. In most cases, boson condensation is required to accommodate the boson density b † iα b iα = 1 in the undoped systems. For example, in the Q 1 = Q 2 state on the Kagome lattice, only for a rather small parameter window of A ij /B ij , a spin liquid state is stabilized. 44,57 Interestingly, this small window appears to be energetically favored in a variational Monte Carlo study of the J 1 -J 2 Heisenberg model, 56 which could explain the quantum spin liquid phase discovered in DMRG simulation. 15 However, in the doped case, b † iα b iα = 1−x where x is the doping level, suggesting a larger parameter range in which the liquid phase is stabilized. This is also consistent with physical intuition. In the slave-fermion mean-field description, in terms of spin dynamics, doping only means replacing S = 1/2 by S = 1/2(1 − x). Therefore doping effectively reduces the spin and increases the effects of quantum fluctuations.

C. Quantum numbers
After the symmetric wavefunctions are constructed on finite-size lattices, their symmetry quantum numbers can be analytically computed. In Table I we summarize the quantum numbers of the three competing phases on symmetric samples (see Appendix C for details). All wavefunctions are SU (2) spin singlets. We find that on 2N × 2N × 2 lattices 58 , the ground state wavefunctions of all the three competing phases always form two-fold irreducible representations(irreps) of symmetry group.
In particular, the two degenerate states in the angular momentum basis (rotational symmetry eigenbasis) exactly form time-reversal images of each other. This is a rather special case of time-reversal symmetry breaking phenomena. Although all three competing phases break time-reversal symmetry in the thermodynamic limit, without the analysis of lattice symmetries, naively one may expect that the time-reversal-related two-fold ground state sector is nondegenerate on finite size lattices due to tunneling. Here the quantum tunnel- ing between the two ground states is forbidden by the lattice rotational symmetry.
One may wonder that in the thermodynamic limit, apart from two-fold degeneracy induced by time-reversal symmetry breaking, there should also be a topological order induced degeneracy in the SCCL phase. In fact we will show in Section IV A that there will be four-fold degeneracy induced by topological order and totally we have eight-fold degeneracy. The ground states of SCCL shown in Table I correspond to a particular one of the four topologically degenerate sectors. The center of mass momentum of the other three sectors is at the three M 1,2,3 points. We believe that on the 2N × 2N × 2 finite lattices the energies of these three sectors are higher than the one shown in Table I, because only in the sector with the center of mass momentum at Γ the minima of spin-1/2 boson dispersion coincide with available momenta in Brillouin Zone; other three sectors are obtained by insertion of π-fluxes, which moves the momenta away from the position of boson minima and should lead to higher energy (see Appendix B).
From Table I one learns that the 32-site sample is the smallest system allowing a sharp distinction 59 between the c-SDW/SCCL phase and the d+id SC phase. However, it is likely that exact diagonalization on the 32site sample is beyond the currently available computing power. This motivated us to perform the 32-site DMRG calculations in Section III C.

III. NUMERICAL SIMULATIONS
Numerical calculations were performed on three samples shown in Fig. 5, each defined with periodic boundary conditions. In all numerics, t = 1 is fixed while J or U are varied.

A. Exact diagonalization on the 8-site sample
We first describe the results for the t-J model on the 8-site sample ( Fig. 5) with 6 fermions. Fig. 6 shows the ground state (GS) energy throughout the physically interesting parameter regime 0.1 ≤ J/t ≤ 2. 60 We find it is twofold degenerate and a spin singlet S = 0. Evaluating directly the matrix elements ψ i |O sym |ψ j of symmetry operators in the GS doublet i, j ∈ {1, 2}, we find that the translation, rotation, inversion, mirror and time inversion properties of the GS doublet match the ones shown in Table Ib to all available digits. The inset shows a more detailed scan revealing a level crossing to a singly degenerate GS below J/t = 0.089(1) which forms a trivial irrep of the symmetry group, but the relevance and nature of this very-low-J state will be studied in future work.
Turning to the Hubbard model, the GS energy on this sample is shown in Fig. 7, in comparison to results for the 32-site sample obtained using the DMRG method. Using ED we find a doubly degenerate ground state in the regime 0 < U/t < 61.31, matching the irrep shown in Table Ib. It is well known that the small-J regime in the t-J model and the large-U regime in the Hubbard model are related by perturbative analysis. Indeed we find that at U/t > 61. 31(1) the ground state forms a onedimensional trivial irrep of the symmetry group, which is consistent with the related level crossing in the t-J model at J/t = 0.089(1).
Finally, we use ED on this small sample as a benchmark for DMRG calculations which perfectly matched the ED energies.

B. Variational Monte Carlo calculations of the d+id superconductor phase in the t-J model
The Hilbert space on which |Ψ d+id (χ, ∆) (Eq. (3)) is defined is too large for direct computation. Therefore we use the Variational Monte Carlo (VMC) technique, within which the expectation values of observables in this state are calculated using: 6,61 where |φ is the considered many-body state, while |R are states in the appropriate Hilbert space which are probabilistically sampled using the first fraction in Eq. (9) as the distribution in a Metropolis algorithm. Concretely, the states |R in the t-J model Hilbert space are given by the spin-occupation basis: where N is the number of sites, s i =↑, ↓, 0 depending on whether the site-i is spin-up, spin-down or empty, c iα annihilates electron of spin α =↑, ↓ at site i, and |0 is the vacuum. There are exactly N F non-empty sites, enforcing the fixed fermion number, and obviously there is no double occupancy. We choose to order the c † iα operators according to site label i, thereby fixing the fermion signs in the |R(s 1 , s 2 , ...s N ) tJ basis. Similarly, in the Hubbard model we have s i =↑↓, ↑, ↓, 0, and where again there are in total exactly N F operators c † iα , and in the obtained |R(s 1 , s 2 , ...s N ) Hubbard we order them according to site label i, keeping the c † i↑ before the c † i↓ for each doubly occupied site i. We focus on the total S z equal to zero sector (in both models), by additionally choosing an equal number of spin-up and spin-down electrons. Note that the DMRG calculation conserves this spin quantum number of a state, so we can work in an S z sector. As discussed in detail in the following, we also measured quantities after projecting the wavefunction to a certain symmetry sector using a projector P , and note here that both the action of the operatorÔ and P are dealt with by acting directly on the R| in Eq. (9).
The optimal value of the single variational parameter, the pairing ∆/χ ∈ R, which minimizes the variational energy, is shown in Fig. 8. For smaller J/t the pairing is too small and harder to determine precisely.
The main signature of the d+id phase is the complex phase of pairing, Fig. 2b. We therefore calculate the pairpair correlation function: the singlet pairing. The pattern from Fig. 2b should be revealed in the long-range physics, so the most interest lies in pairs of nearest-neighbor bonds ij and kl which are as far from each other as possible. Table II reveals that the pattern indeed occurs and becomes weaker with decreasing J/t. The spin-spin correlation function is very short-ranged as expected, so we do not present it in detail. 62

C. DMRG simulations on the 32-site sample
We have used DMRG to obtain the ground state (GS) of the periodic 32-site sample (Fig. 5). Our calculation is based on the open-source DMRG software package ITensor, 8,63 where the periodic two-dimensional shape of our samples was implemented simply by introducing long-range hopping (of same size t) in the native DMRG one-dimensional representation of the system. The limit on dimension of MPS matrices was between 10.000 and 11.000. We find truncation errors around (2 ∼ 7) · 10 −4 , depending on model and parameter regimes. Although such error values seem too large in view of general DMRG performance, in this work we found it appropriate to apply a different physical criterion for convergence, namely, that the expectation values of symmetry transformations allow a clear assignment of quantum numbers to the ground state; additionally, when appropriate, in measurements we projected the GS to a sector having some quantum numbers fixed, to effectively get closer to the true GS. This approach will be described in detail below. Appendix F presents further details on our DMRG setup and convergence.
Focusing first on the t-J model, we find a very precise quantization of the inversion operator expectation value in the GS, as shown in Table III. For this sample there is a sharp transition at J/t = 0.80(2) at which the low-J ground states (blue phase in Fig. 3), having inversion −1, switch to high-J ground states (red in Fig. 3), which are in the +1 representation of inversion. Due to change of symmetry quantum numbers, we expect this to be a first order phase transition in thermodynamic limit. Given that GS is in a representation having inversion +1 (−1), and since there is no reason for additional degeneracy except due to time-reversal, the 60-degree rotation operator (C 6 , with C 3 6 =Inversion) should be represented by one of numbers A crucial subtlety here is that the DMRG calculation automatically provides a real-valued wavefunction for our real Hamiltonians. This DMRG wavefunction will be denoted as |ψ in the following discussion. If |ψ gives the converged true ground state, it must be an equal superposition of two conjugate partners in a two-dimensional irrep when C 6 is represented by a complex number. Simple calculation shows that generally the C 6 expectation value for a converged real ground state wavefunction must be one of {−1, −1/2, 1/2, 1}, corresponding to the four possible irreps of the symmetry group respectively: the C 6 -odd one-dimensional irrep, the two-dimensional irrep as shown in Table Ib, the two-dimensional irrep as the c-SDW/SCCL shown in Table Ia, and the trivial onedimensional irrep. Note that the DMRG we applied here can be viewed as a variational wavefunction technique in real space, in which lattice symmetry is not implemented at all.
Therefore we use the C 6 expectation value as a physical criterion for successful convergence of the DMRG wavefunction. Namely, if Ψ|C 6 |Ψ , with |Ψ defined shortly, is found to be one of the four values: {−1, −1/2, 1/2, 1}, the DMRG has successfully converged. On the 24-site sample (see Sec.III D), we find that using |Ψ = |ψ the C 6 expectation value is well converged in the parameter regimes of interest. However on the 32-site sample, in order to improve convergence, we project |ψ to the sector with center of mass momentum equal to Γ; namely, we use |Ψ = P Γ |ψ as the wavefunction in MC measurement of the C 6 expectation value, Eq. (9), where P Γ is the projection operator into the Γ-sector. (We also check that |ψ has a big portion in the Γ-sector for all parameter values so that this projection is not creating unphysical artifacts.) Fig. 9a demonstrates the result that in the low-J regime (0.1 < J/t < 0.8), the rotation expectation value is indeed consistent with 1/2 on the 32-site sample. Therefore the GS irrep in this regime is the same as the c-SDW/SCCL phase as shown in Table Ia. However, for the lowest values, 0 < J/t < 0.1, the C 6 does not converge to either of {−1, −1/2, 1/2, 1}, and this also happens for the 24-site sample for 0 < J/t < 0.07; on the other hand, the 8-site exact diagonalization shows a singlet ground state for 0 < J/t < 0.089(1). All this evidence suggests the existence of a different quantum phase in this lowest J regime. Given that such lowest J regime is not the most interesting for correlated materials, we leave it for future work, and focus exclusively on values J/t > 0.1.
In the high-J regime (J/t > 0.8) unfortunately the C 6 is close to zero and far from any of the {−1, −1/2, 1/2, 1}, which indicates that the 32-site DMRG GS for J/t > 0.8 has not converged well enough; it cannot give reliable information about correlations. Nevertheless the inversion quantum number for J/t > 0.8 is found to be accurately +1, consistent with the d+id SC and sharply distinguished from the 0.1 < J/t < 0.8 value −1 (see TableIII). In the following discussion and in the next Section, the J/t > 0.8 phase is actually confirmed to be the d+id SC using complementary variational Monte Carlo results as well as DMRG on the smaller 24-sample which has no such issues with convergence.
Energetics of the DMRG GS of t-J model are shown in Fig. 10. The energy of the single-parameter variational wavefunction discussed in Section II A is quantitatively compared to the DMRG energy, showing that the d+id candidate wavefunction captures more than 97 ∼ 99% of DMRG GS energy throughout the high-J phase. In addition, the energy of the d+id variational state deviates significantly in the low-J phase.
To further identify the nature of the DMRG GS, we consider spin-spin and pair-pair correlation functions.
The expectation values are obtained using the Monte Carlo (MC) technique, Eqs. (9), (10), (11), using between 300 and 1000 MC measurements with 40 MC steps between each measurement and with a 500 MC step thermalization. Further, the measurements are averaged across 64 independent MC runs, and the measurement errors in this paper represent the error of the mean. To correctly calculate observables we need to choose a particular rotation sector from the |Ψ , since the DMRG mixes rotation sectors by selecting a real wavefunction as discussed above. According to Table I this projection to a rotation eigenstate means breaking the timereversal symmetry, which should naturally happen in the thermodynamic limit. For all measurements on this 32site sample, in the phase with inversion −1 we choose the exp(−iπ/3) sector of C 6 . (Note that for the +1 inversion phase (J/t > 0.8) in the t-J model the C 6 is not converged, so we do not use it.) More precisely, the correlation functions we study next are obtained as ψ|P Γ P C6Ô P Γ P C6 |ψ in Eq. (9), with P C6 projecting into the desired rotation eigenspace.
We calculate the spin-spin correlation function by setting the observableÔ = S z (i)S z (j) with some sites i, j. Although the short-range physics dictates that nearest neighbor i, j correlations grow with J/t (this is indeed observed), we are interested in long-range physics and therefore choose the farthest pair i, j, Fig. 12a, finding that this correlation grows with going deeper into the low-J phase. Fig. 11 demonstrates the spin-spin correlation pattern for i, j bonds of all lengths, revealing a pattern consistent with tetrahedral spin correlations (Fig. 2a) in the low-J phase. The overlap of spin vectors in the tetrahedron predicts a ratio of −1/3 in the correlation when the spins at sites i, j are parallel compared to when they are not. Our measurement of this ratio for the farthest possible site pairs i, j is consistent with the prediction, Fig. 13a. We also calculate the spin chirality S i ·( S j × S k ) for the smallest triangle in the honeycomb lattice, see Fig. 14.
The magnitude of chirality of around 0.01 is consistent with magnitude of nearest neighbor spin-spin correlation of ∼ 0.05. 64 In the low-J phase, where spin indicates the c-SDW/SCCL state, the pair-pair correlation function is extremely short-ranged and beyond nearest bond pairs hard to distinguish from zero within our numerical precision (see Appendix F).
Let us now turn to the Hubbard model on the 32site sample, having ground state energy presented in Figure 7. Table III demonstrates our result that for a very wide range of parameters 1 < U/t < 40 the expectation value of inversion operator is very accurately quantized to −1. Figure 9b shows that for all U/t 6 we find a satisfying agreement of 60-degree rotation expectation value C 6 with +1/2. The same figure shows the influence of projection to Γ momentum, i.e., using |Ψ = P Γ |ψ , which significantly improves this agreement. It is not surprising that convergence worsens for low U/t, due to the existence of many low energy states, but we believe it is limited by our maximal available m. For instance, at U/t = 4 the C 6 expectation with projection to Γ momentum improves from 0.10(2) at m = 8.000 to 0. 16(2) at m = 10.500 (see Appendix F). In fact, using degenerate perturbation theory on the 32site sample around free electron state t = 1, U = 0 (Appendix G), we find the same quantum numbers as for the c-SDW/SCCL state. We therefore expect that the c-SDW/SCCL ground state quantum number persists through the whole range 0 < U/t < 40 on this sample.
The spin-spin correlation function (again, for this sample we take P Γ P C6 |ψ in the exp(−iπ/3) sector of C 6 ) throughout the entire well-converged and physically interesting regime U/t 4 is qualitatively the same as in the c-SDW/SCCL phase of t-J model, Fig. 11a (see also Appendix F). Quantitatively, Figs. 12c, 13b show how the long-range tetrahedral spin pattern describes this phase very well, and strengthens with growing U/t. This is consistent with the mapping between low-J and large-U models, confirming the c-SDW/SCCL nature of the phase in both models.

D. DMRG simulations on the 24-site sample
The fully symmetric 24-site sample, Fig. 5b, is large enough to provide some longer-range physics information, but small enough to allow excellent DMRG con- vergence and precise measurements (see general discussion of our DMRG convergence criteria in previous Section). It may even be suitable for exact diagonalization numerical simulations using currently available computing power. We therefore investigated the quantum numbers of the three competing states, c-SDW, SCCL and d+id superconductor, on this sample, and found that unfortunately all these phases share the same quantum numbers as in Table Ib. Therefore, a smooth crossover takes place in the t-J model. To support the claim that the high-J phase observed on the 32-site sample is the d+id SC, in this section we will consider the t-J model on the 24-site sample and show that it clearly exhibits a change in its correlation properties from the characteristic c-SDW/SCCL to the d+id SC behavior as J/t is increased within the 0.1 < J/t < 2 parameter region. (We will not discuss the Hubbard model on this sample.) Table III shows the very precise quantization of inversion to +1 in the DMRG GS in the entire region 0.1 < J/t < 2 (as explained in previous section, we do not further discuss the 0 < J/t < 0.1). The GS is almost entirely in the Γ momentum sector, so we use |Ψ = |ψ and find that C 6 is very close to −1/2 (Fig. 9a) in the entire considered parameter region. This corresponds to quantum numbers in Table Ib. The energetics in Fig. 10 shows 99% agreement with variational d+id wavefunction at larger J/t, which significantly worsens as we go to lower J/t, indicating the crossover to c-SDW/SCCL state.
Due to smaller sample size and the fact that Γ momentum projection is unnecessary, we could use 10.000 MC measurements in correlation functions, significantly reducing the statistical error. The correlation measurements are all done in the exp(−i2π/3) rotation sector, corresponding to the +1 value of inversion.
Figs. 11c,d contrast the spin-spin correlation at J/t = 0.2 and 2, respectively. The former is clearly consistent with tetrahedral spin correlations. On the other hand, the J/t = 2 case exemplifies a completely different, and much shorter ranged, spin correlation pattern. Fig. 12c quantifies the weakening of the tetrahedral pattern, which rapidly drops to zero with J/t growing towards 1, indicating the existence of the crossover.
Complementary information is found in the pair-pair correlation function (see Eq. (12)), presented in Table II. At largest value, J/t = 2, the correlation pattern matches the ideal pattern of Fig. 2b with percent precision. By the time we reach the lowest value J/t = 0.1, the overall correlation amplitude drops fivefold, the different pairs' correlation varies in amplitude significantly, and their relative phase of 2π/3 drops to 0.07·2π. The table shows that these results match the evolution of pair-pair correlation in the variational d+id wavefunction, up to an overall amplitude difference in the correlation function. Altogether, the existence of crossover between c-SDW/SCCL and d+id SC in the t-J model on this sample is clearly confirmed.

IV. THE SPIN-CHARGE-CHERN LIQUID
A. Low energy effective theory: Parton construction and the K-matrix formulation In two spatial dimensions, a description of Abelian topological order can be given by Abelian U (1) N Chern-Simons theory. [65][66][67] The low energy effective Lagrangian relevant for us has the following generic form where µ, ν, λ = 0, 1, 2 in 2+1D and summation over repeated indices is implied. K is a symmetric N × N matrix with integer entries. A quasiparticle in this theory is described by an N component integer vector l, whose components determine the N U (1) gauge charges of the excitation. The particle couples to internal gauge field a µ as −a I µ l I j µ . Here, j µ is the 3-current for a single quasiparticle.
The quasiparticle statistics can be easily read out by integrating out a I µ . The self(exchange) statistics of a quasiparticle l is given by its statistics angle while the mutual(braiding) statistics of a quasiparticle l and l is characterized by Quasiparicles generally have anyonic statistics and are thus nonlocal. However, there is a special type of quasiparticlel = Kl, where l ∈ Z N .l is mutual boson to all other quasiparticles, so it can be viewed as a local excitation, in the topologically trivial sector. Examples include electron excitations of fractional quantum Hall systems and spin-1 magnons in Z 2 spin liquids. Two quasiparticles whose difference is in the trivial topological sector should be considered as being in the same topological sector. Further, the ground state degeneracy (GSD) on a torus is 68,69 which is equal to the number of topological sectors (quasiparticle types). In the following we will construct the effective field theory for SCCL state. In the slave-fermion approach (5), the electron is separated into a bosonic spinon and a fermionic holon. The fermionic holons fill a C = 1 Chern band, which can be described by a Chern Simons term where a 2π flux (vortex) of gauge field a f µ is a holon particle. On the other hand, a pair of bosonic spinons can be described as a 2π flux of an internal gauge field a p µ . (In the liquid phase, there is a superfluid of spinon pairs, not of spinons.) Finally, the holon and spinon are glued together to form the electron by a U (1) gauge field a c µ . This U (1) gauge field acts as a constraint in the Lagrangian where the factor 2 accounts for pair of spinons having twice the internal gauge charge of a single spinon. Now, we define a I µ = (a f µ , a p µ , a c µ ), leading to and we find  = (0, 1, 1). Notice that the holon f = (1, 0, 0) and spinon b differ by an electron, so they belong to the same topological sector.
SCCL is however not fully described by its topological properties. Symmetry interplays with topological order, leading to symmetry fractionalization (see, e.g., Refs. 49,[70][71][72][73]. Within the K-matrix formulation, it is possible to assign quantum numbers of onsite symmetries, e.g., charge and spin, to quasiparticles. 74 Namely, we define the charge vector t c = (1, 0, 0) and S z vector t Sz = (1/2, −1, 0), so that a I µ couples to external test gauge fields as where A c µ is the gauge field that couples to electric charge, while A Sz µ couples to S z . Quasiparticle l carries electric charge t t c K −1 0 l and carries S z = t t Sz K −1 0 l. We can now identify spinons b ↑/↓ as (0, 0, ±1), while holon f remains just (1, 0, 0). It is straightforward to see that holon indeed carries electric charge 1 and S z = 0, while b ↑ (b ↓ ) carries no charge and S z = 1 2 (− 1 2 ). Electron e ↑(↓) is simply the bound state of f and b ↑(↓) , and it is in the topologically trivial sector. The vison, expressed by (0, ±1, 0), carries charge ± 1 2 , and since it has statistical angle π 4 , the vison can be viewed as 'half holon'. Bound state of spinon and vison carries both charge ± 1 2 and spin ± 1 2 , with statistical angle 5π 4 . There exists another state, described byK 0 ≡ −K 0 , which is related to the above state by time reversal. In this state, vison excitation has statistical angle − π 4 while bound state of spinon and vison has statistical angle − 5π 4 .

B. Modular Transformations and Rotation Quantum Numbers
S and T matrices obtained from modular transformations of ground states on torus are believed to encode quasiparticle braiding and exchange statistics. 68 Additionally, as pointed out by Refs.39 and 42, it seems that if system has C 6 rotation symmetry the ground state quantum numbers of C 6 equal the eigenvalues of ST.
The relation between modular S,T matrices and the rotational symmetry of a topologically ordered phase may be understood as follows. First note that S,T matrices are in principle measurable quantities in practical model Hamiltonians. In particular, given a topologically ordered phase in 2+1D with its topologically degenerate ground sector on torus T 2 , one can firstly find a minimally entangled state (MES) basis 39 . For instance, for the S-matrix element between two MES |Ξ i and |Ξ j : S ij , one can perform the following thought numerical measurement. Because the topological properties do not depend on local geometry, we can assume that these ground states live on a square with periodic boundary conditions. Then one can consider the state rotated by 90 • around the square center: R 90 • |Ξ i . Because R 90 • |Ξ i and |Ξ j belong to the same topological phase, in the absence of symmetry there should exist a Hamiltonian path H(τ ) (τ ∈ [0, 1]) such that |Ξ j (|Ξ j ) are the ground state of H(0)(H(1)), and the ground state sectors of H(τ ) are adiabatically connected. One can then define a projection operatorP τ into the ground state sector of H(τ ) for any given τ .
The many-body quantum amplitude related to the adiabatic time-evolution process of the S-transformation can be computed as s ij ≡ Ξ j |P (N −1)/N · ... ·P 2/N · P 1/N R 90 • |Ξ i as N → ∞. This computation is a realization of the topological quantum field theory timeevolution. In particular, if the system has a 90 • rotational symmetry, the Hamiltonian path H(τ ) can be con- 15. Symmetry quantum numbers for C6 rotation, calculated analytically (Appendix C) for the fourfold topologically degenerate ground state sector of the SCCL phase on different lattice sizes. C6 eigenvalues are plotted on the unit-circle in complex plane. The sets of eigenvalues differ by overall phase between different system sizes (N is integer, and X × Y labels number of unit-cells along a1, a2), although all systems are topologically a torus.
veniently chosen to be a constant: H(τ ) = H(0). In this case, s ij can be simply computed as the R 90 • transformation matrix in the MES basis: We expect that this quantum amplitude s ij is related to the S-matrix elements S ij at most by an overall ambiguity U (1) phase e iθ , which is due to the nonuniversal local physics in the time-evolution, and a phase e iφi−iφj which is due to the gauge choice of |Ξ i ,|Ξ j . Even with these ambiguities, based on the above argument, it is clear that in a 90 • rotational symmetric system, the R 90 • eigenvalues in the topologically degenerate ground state sector can be determined by the eigenvalues of the Smatrix, up to an overall U(1) phase factor. Similar consideration for a 60 • rotational symmetric system leads to the conclusion that the C 6 eigenvalues in the topologically degenerate ground state sector can be determined by the eigenvalues of the matrix product ST, up to an overall U(1) phase factor. In addition, it has been proposed that this U(1) phase factor is simply unity 39 which is consistent with numerical simulations on several model Hamiltonians 17,41,43 .
However, we find in the SCCL phase on the honeycomb lattice, the C 6 eigenvalues in the ground state sector and the eigenvalues of ST differ by an overall U(1) phase factor that is system-size dependent. In particular, one can obtain S and T matrices from our K-matrix. According to Ref 42 , using Eq. (20) and choosing four quasiparticle vectors as (0, 0, 0), (0, 0, 1), (0, 1, 0), (0, 1, 1), one obtains where ξ,η are U (1) phase factors. Although in Ref. 42 these phase factors are fully determined using modular transformations on fractional quantum Hall liquid analytic wavefunctions on torus, their values are not important for the following discussion.

C. Gapless edge states and experimental signatures
We will first derive the edge theory of SCCL using the effective field theory from previous subsection. We consider two cases of symmetry on the edge: 1) Charge conservation and spin-rotations around S z (group U (1) c × U (1) z ); and to capture more of the spin-rotation symmetry 2) Charge conservation, S z rotations, and πrotation around a perpendicular axis. The second case is detailed in in Appendix D, but in both cases we find a gapless chiral holon edge mode, which differs from the gapless chiral electron mode of the c-SDW state. We therefore propose several experimental signatures for distinguishing c-SDW and SCCL states in the last two subsections.
The effective action describing edge excitations of Abelian Chern-Simons theory can be derived from gauge invariance of Lagrangian Eq. (13) expanded by higher order (Maxwell) terms, on a manifold with boundary. 75 The edge physics is captured by N chiral boson fields {φ I φ I + 2π|1 ≤ I ≤ N }: Here, V I,J is positive definite constant matrix, which depends on system details. The number of right movers n + and left movers n − are given by the signature of K. The commutation relations between these chiral boson fields are fixed by the first term, and describe the following Kac-Moody algebra 75 : There is a one-to-one correspondence between quasiparticles in the bulk and chiral boson fields living on the edge. Operator V l = exp(i I l I φ I ) creates quasiparticle l on the edge. Generic action for scattering takes the form of Higgs terms: wherel are local bosonic excitations, which can be expressed as Kl for some integer vector l. However, in the presence of symmetry, chiral boson fields may transform nontrivially under symmetry operations, and some Higgs terms may be forbidden in the symmetry-preserving edge. 76

Edge modes with U (1)c × U (1)z symmetry
Now let us turn to edge theory for SCCL. As the set of independent local excitationsl we choose the columns of K 0 . The 1st column of K 0 matrix is an electron, carrying charge 1 and spin 1/2, while the 2nd column is boson pair with spin 1. They both transform nontrivivally under S z rotations. The 3rd column is bound state of two vison and a charge -1 holon, which is a trivial boson carrying trivial quantum number of S z and charge. Thus, the only Higgs term allowed by this symmetric boundary is This term gaps out two counter-propagating edge modes, and leaves the gapless chiral boson mode φ 1 on the edge. Thus, the edge theory of SCCL can be modeled as 1d chiral fermion liquid of spinless holons. Even with added π-rotation around an axis perpendicular to S z , as shown in Appendix D, φ 1 remains the only gapless edge mode. We therefore found that the edge with charge conservation and any of above spin-rotation symmetries has a chiral fermion liquid of spinless holons. The c-SDW edge on the other hand has a chiral fermion liquid of electrons. We therefore next propose tunneling experiments to distinguish these two phases.

Point junction
In this part, we will discuss the experimental signature of transport through a weak tunnel junction connecting a metallic/singlet SC lead to c-SDW/SCCL (Fig. 16). Our results of tunneling conductance are listed in a table in Fig. 16. Below we will find the same exponents for the voltage dependence of the conductance.
These scaling forms, and therefore the experimental signatures, should hold in the regime of weak tunneling, G e 2 /h. More formally, the weak-tunneling condition corresponds to the assumption T, V T K , where T K is a characteristic energy scale of the junction depending on details of the point contact.
The total Hamiltonian can be modeled as a sum of three pieces where H 0 is Hamiltonian for c-SDW/SCCL, H lead is Hamiltonian for SC/Metal lead, and H tunn describes tunneling through point contact. For the most general case, we can write where O 0 is electron or singlet pair annihilation operator on c-SDW/SCCL side, while O lead is the corresponding operator in the lead. Before calculating tunneling conductance, it is instructive to consider a simple renormalization group (RG) transformation, which tells us how the tunneling amplitude t varies with the energy (or temperature) scale. 77 Assume O 0 ∼ τ −δ0 and O lead ∼ τ −δ lead , where τ is imaginary time. Consider an RG step which integrates out Matsubara frequencies between Λ/b and Λ, where Λ is a high frequency cut-off. Then the RG equation for t is given to leading order by where δ = δ 0 + δ lead . At nonzero temperature, the RG flows are cut off by T (T V ), leading to t ef f ∼ tT δ−1 . One expects tunneling conductance to vary as t 2 ef f , which gives the result We now present the case of metal/SCCL junction in detail, referring the reader to Appendix E for the other cases listed in table of Fig. 16. (Note that the scaling for c-SDW/SC junction follows directly from Ref. 78.) Due to the spin gap on the boundary of SCCL, single electron tunneling will be exponentially suppressed at low temperatures. So, the leading contribution is from singlet pair tunneling, and in Eq. (28) we have: where the product of holon operators f † in O 0 represents annihilation of a local singlet pair of electrons on SCCL edge due to the presence of bosonic spinon pairing (see Eqs. (5), (6)), while coherence length ξ appears due to the Pauli principle. So, δ lead = 2δ F L = 1, where we used that the scaling dimension δ F L = 1/2 for Fermi liquid system in any dimension. 1 The operator f (x = ξ)f (x = 0) has the same scaling dimension as operator f (x = 0)∂ x f (x = 0), giving δ 0 = 1 + 2δ F L = 2, where the holon operator on the edge scales with δ F L since it forms a chiral fermion liquid analogous to the one on the edge of integer quantum Hall systems. 79 This leads to the announced G(T ) ∼ T 4 for this junction.
One generally expects that the voltage-dependent conductance G(V ) scales in the same way as G(T ). We checked that this is true using a perturbative calculation (i.e., the Fermi golden rule) of the nonlinear currentvoltage (I − V ) characteristic in the regime T V . The calculation details for all junctions are presented in Appendix E.

Line junction
A more common setup in experiments is the line junction, 80 which can be viewed as a large number of weakly coupled point junctions, as sketched in Fig. 17. Here, and throughout our discussion on the line junction, weakly coupled means where N is the total number of point junctions, while T (n) K is the characteristic energy scale determined by details of the n-th junction 80 . Physically, this weakcoupling condition in the line junction means that the regime of weak tunneling, G e 2 /h, is available, at least for low enough voltage (see further below). The special case of c-SDW/metallic lead junction is left for the end of this subsection, since it is much simpler to analyze and does not require such assumptions.
The number α = 0 which appears below is simply the value of exponent in Fig. 16 for the considered combination of quantum state and lead. The special case of c-SDW/metallic lead junction has exponent α = 0, and is discussed at the end.
First, let us consider the T V regime. We will find that for small voltages, the scaling of conductance can distinguish the c-SDW and SCCL in the same way as table in Fig. 16.
The expression for current-voltage characteristic we obtain (see Appendix E) is where the voltage difference V ≡ V R −V 0 (Fig. 17), α = 0 is the exponent in the point junction scaling G ∼ V α (table in Fig. 16), and the effective T K is the single parameter describing the line junction and incorporating all the T (n) K as well as their fluctuations: (Note that the definition of T K also depends on scaling exponent α.) The above α = 0 conductance result holds for all values of V, T K at T V , as long as the assumptions used to derive the expression hold, namely, each individual point contact is weakly coupled. This just means V T

(n)
K for all n. However, the effective T K can be much smaller than all T (n) K in a long line junction (large N ). Therefore, let us examine the tunneling conductance G in two regimes: T V T K and T, T K V . For the first regime we get: manifesting the same scaling form as that in point contact junction. (Note that still α = 0.) On the other hand, in the regime T, T K V we get The derivation for V T regime is similar, and we reach the same final conclusions as for previous case. The current-voltage characteristic in this regime is: see Appendix E. In this regime a characteristic energy scale T K , analogous but different from T K , describes a point junction having α = 0. For a given combination of quantum phase and lead forming the junction, we expect the ratio T K /T K to be a universal number of order 1.
With this in mind, we again consider two regimes: V T T K and V, T K T . We get while in the other regime: recalling that α = 0 is the same exponent found for the point junction, Fig. 16. Details are in Appendix E.
Concerning the c-SDW/metallic lead junction, which has scaling exponent α = 0 in the point junction, a simple calculation reveals that the line junction conductance is constant, G = e 2 h · c, with 0 < c ≤ 1 a non-universal constant describing the line junction (see Appendix E). c = 1 corresponds to the regime in which the chiral electron edge modes are equilibriated with the lead.
We conclude that the line junction tunneling conductance G(V )(G(T )) can distinguish between c-SDW and SCCL in the regime where G e 2 /h and T V (V T ), which corresponds to weakly coupled line junctions. In this case the edge modes are not thermally equilibriated with the lead. For example, in this regime, the zero bias tunneling conductance for the c-SDW/metallic lead line junction is temperature-independent while for the SCCL/metallic lead line junction it is expected to scale as T 4 . On the other hand, although G = e 2 /h is one experimental signature of the quantum anomalous Hall effect in the c-SDW phase, we find that even in the SCCL phase the universal value G = e 2 /h can be realized, e.g., in the regime of Eqs. (37), (40), where it represents the thermal equilibriation of chiral holon edge modes with the lead. Therefore the G = e 2 /h is not a unique property of the c-SDW phase.

V. DISCUSSION AND CONCLUSIONS
This paper studies the phase diagrams of correlated electronic models on the honeycomb lattice at 1/4 doping, using a combination of analytical construction of quantum wavefunctions and various numerical simulations. Interestingly, all phases appearing in our main results, the phase diagrams in Fig.3, are interaction-driven topological phases. In particular, we find that either the c-SDW state or the SCCL phase occupies the majority of the realistic parameter regimes for correlated materials. In the present study, due to the limitation of sample sizes, we cannot sharply distinguish these two phases in the phase diagrams. Distinguishing them in numerical simulations requires careful finite size scaling, which we leave as a subject of future investigation. However, we study the sharp signatures of c-SDW and SCCL phases in transport experiments, which can be used to identify and distinguish these phases in candidate materials.
The method applied here, namely using lattice quantum numbers to sharply distinguish competing quantum phases, is not limited to the models studied in this paper. In particular, in time-reversal symmetry breaking phases, the ground states often form non-trivial multidimensional irreps of the lattice symmetry groups. When this happens, the analytical understanding of the nontrivial irreps can be used to identify/distinguish candidate quantum phases in numerical simulations.
From a general point of view, what are the possible candidate phases in correlated electronic systems at generic fillings? First, charge inhomogeneity is always a possibility. For instance, stripe-like charge modulations have been observed in numerical simulations of the t-J model on the square lattice. 81 Assuming charge being homogeneous, incommensurately filled systems and commensurately filled systems are quite different at the conceptual level. Generally speaking, in order to accommodate an incommensurate filling, the system could either develop superconductivity, or the doped charges could form a Fermi surface. 82 In any case the system is expected to be a charge conductor in the bulk. However, at commensurate fillings, the system has a third option: The charges could condense into many-body states without causing a charge inhomogeneity or superconductivity, and a bulk energy gap of charge-excitations can be generated. We term this third scenario as the charge-insulator scenario.
In conventional quantum phases in which Luttinger's theorem 83 is valid, the charge-insulator scenario must be accompanied with translational symmetry breaking such as long-range magnetic ordering. The c-SDW phase belongs to this situation. However in exotic quantum phases in which fractionalization occurs, translational symmetry does not need to be broken. For instance, the SCCL phase is a translationally invariant charge insulator. Other examples include the recently studied Fractional Chern insulators, [84][85][86][87][88][89][90][91] which are symmetric manybody states that exist in models with commensurately filled nearly-flat bands in the presence of strong interactions.
One goal of this paper is to investigate the competition between the superconductivity and the charge-insulator phases in commensurately doped correlated systems. Exactly at the 1/4 doping, we find that the charge-insulator phase(c-SDW/SCCL) occupies the majority of realistic regimes of the models that we investigated. Meanwhile, although the d+id superconductor phase is found only at J/t > 0.8 in the t-J model, as a variational state, it captures > 90% of the ground state energy even in the regime 0.1 < J/t < 0.8(see Fig. 10). Therefore the d+id state serves as a nearby competing phase.
As doping deviates away from 1/4 slightly, the extra electric charges need to be absorbed by excitations in a charge-insulator. In the SCCL phase, these charge excitations form a finite density of anyons: v, b α v or fermionic chargeon-f ; while in the c-SDW phase, these excitations would be a finite density of electronic quasi-particles. However, the d+id superconductor state, as a charge superfluid, can absorb extra electric charges without causing excitations. In the regimes in which c-SDW/SCCL phase is realized at 1/4 doping, we expect that the ground state is likely to be the d+id state as the doping is tuned away from 1/4 by a finite amount.
One may wonder that due to the Mermin-Wagner theorem, the long-range magnetic order cannot be observed at finite temperatures in the c-SDW phase. In addition, our discussion of the low temperature tunneling conductance in the c-SDW phase did not consider this thermal fluctuation effect. However, in an ideal SU (2) symmetric system the correlation length of the magnetic order diverges exponentially at low temperatures. Thus even a tiny spin-orbit coupling strength would pin the magnetic order at low temperatures in realistic materials, which justifies our treatment.
Recently there has been a lot of interest in the understanding of interplays between global symmetry and topological order, which have been named as "symmetry enriched" phenomena (see, e.g., Refs.49,[70][71][72][73][92][93][94][95]. In the present study, the SCCL phase serves as a new example of a symmetry enriched topological phase which could be realized in materials. In the SCCL phase, the symmetry enriched phenomena include the charge-1/2 spin-neutral anyons with statistical angle π/4 and the gapless chargeon chiral edge states. And the latter one has direct experimental signature as G ∼ T 4 in tunneling conductance experiments. Our results are relevant for many correlated materials on the honeycomb lattice. Doped graphene, in which the long-range Coulomb interaction is screened, is an intermediately correlated material that may be modeled by the Hubbard model on a honeycomb lattice with U/t = 2 ∼ 3. 36 InCu 2/3 V 1/3 O 3 is a strongly correlated spin-1/2 antiferromagnet on the honeycomb lattice. 96 However, doping these systems up to 1/4 remains experimentally challenging but may be achievable in a foreseeable future due to the progress of experimental techniques on thin films. 97,98 In addition, recently a new route for realization of honeycomb lattice thin films was proposed, based on growth of (111) bilayers of perovskites. [99][100][101] For instance, after trigonal lattice distortion is included, a 1g -active compounds may be realizations of single-band correlated systems on the honeycomb lattice. 102 In addition, cold-atom optical lattices can be used to realize the Hubbard model on the honeycomb lattice. 103,104 We thank Fa Wang, Yuan-Ming Lu and Satoshi Okamoto for helpful discussions.
The DMRG calculations were performed using the ITensor library, http://itensor.org/. This study is supported by the Alfred P. Sloan foundation and National Science Foundation under Grant No. DMR-1151440. We thank Boston College Research Service for providing the computing facilities where the numerical simulations were performed. The symmetry group (SG) of our honeycomb lattice model is generated by the following symmetry operations (shown in Fig.1):(1) Translations T 1,2 by Bravais lattice vectors a 1,2 ; (2) The π 3 -rotation C 6 around theẑ axis through the honeycomb plaquette center; (3) Mirror reflection with respect to thex −ẑ plane combined with the time-reversal operation, labeled asσ. Note thatσ is an antiunitary symmetry since it includes time-reversal operation. It acts on the Hamiltonian through a combination of a unitary symmetry operation and complex conjugation C.

(A1)
The multiplication rules of the above SG are completely determined by the following algebraic relations: σC 6σ C 6 = e, where e represents the identity element of SG.

Appendix B: Parton construction of symmetric quantum wavefunctions
In this section, we use the slave-fermion method to construct the projective wavefunction of c-SDW/SCCL, and a slave-boson method to describe d+id SC.

c-SDW/SCCL states
In this Section we will consider all mean-field Ansätze allowed by the projective symmetry group construction, and pick out one that correctly describes the c-SDW/SCCL states.

a. Projective symmetry group analysis
The projective symmetry group (PSG) 44,49,50 classifies different mean field Ansätze, and we will briefly review it and apply it here. Although projective wavefunctions are invariant under the symmetry group action (listed in Appendix A), the mean field wavefunction before projection can still explicitly break symmetry. In fact, due to the U (1) gauge field that glues spinon and holon together, a mean-field wavefunction only needs to be invariant under a combined symmetry and gauge transformation. Also, there is a many-to-one correspondence between mean-field states and physical electron states: Any two parton mean-field states related to each other by a U (1) gauge transformation e iφ(r) correspond to the same electron state.

b. Wavefunction for c-SDW/SCCL
There are further constraints on a mean field Ansatz of the c-SDW/SCCL. First, we want quarter doped holons to fill a Chern band, which will lead to the anomalous quantum Hall response. This requires at least doubling of the unit cell. So we only consider the π-flux states having p 1 = 1. In this case, we double the unit-cell in x direction. When two Ansätze are a time reversal pair, we only need to consider one of them.
It turns out that p 1 = 1, p 6 = 0, p 7 = 0 gives the mean-field Ansatz for c-SDW/SCCL. We first construct the mean-field Hamiltonian with NN and NNN hopping/pairing. This particular PSG solution partially fixes the phases of mean-field parameters. The pattern is shown in Fig. (4). After solving Bogoliubov equations for bosons (spinons), we find that the boson band minima lie at ±(π/2, π) of the reduced Brillouin zone. Now, we are able to construct the wavefunction from the mean-field Hamiltonian. Let us consider the most general form of the Hamiltonian in momentum space. Spinon Hamiltonian has BCS form where β k = (b 1k↑ , . . . , b nk↑ , b † 1−k↓ , . . . , b † n−k↓ ) t is Nambu spinor in momentum space, and n is the number of sublattices. A(k) and B(k) are n × n matrices, the Fourier transforms of pairing and hopping, respectively.
We can use M (k) ∈ SU (n, n) for diagonalizing D(k) to get the spectrum of spinons. Expressing it is not hard to derive the BCS-type wavefunction for bosonic spinons as For the c-SDW/SCCL ansatz in Fig.4, there are four sites in one unit cell so n = 4 in this case. The boson condensation occurs (i.e., long-range magnetic order is established) when the boson band minima at ±(π/2, π) touch zero. When this happens, the zero energy modes satisfying D(±(π/2, π))Ψ(±(π/2, π)) = 0 determine the magnetic ordering pattern. They are found to be (in one of the two degenerate ground states): The general boson condensate takes the form: β (π/2,π) = c 1 Ψ[(π/2, π)] and β −(π/2,π) = c 2 Ψ[(π/2, π)], where c 1 , c 2 are two complex numbers. Here among the four real parameters in c 1 , c 2 , one of them, |c 1 | 2 + |c 2 | 2 , controls the magnitude of the magnetic order parameter. A different choice of the other three real parameters can be shown to generate a global SU (2) spin rotation in the spin space. The real space magnetic order pattern is nothing but the tetradedral pattern with the chirality shown in Fig.2. The other degenerate state can be obtained by time-reversal transformation. Now let us look at the fermionic holon part. Hamiltonian of holons is free fermion hopping model, where ψ(k) = (f 1k , . . . , f nk ) t , and n is band index. Using W (k) ∈ SU (n) to diagonalize h(k), we get where d ik = W ij (k) † f jk . Fermions fill bands from the lowest to the i-th, depending on doping. In the case of c-SDW/SCCL phases, the doped holon fills the lowest band. On the mean-field level, it is straightforward to show that the holon real hoppings on the nearest neighbor and second neighbor give a band structure with Dirac points located at ±(π/2, π) in the lowest two bands. The imaginary hoppings (see Fig.4) on the second neighbor open energy gaps at the two Dirac points and the resulting lowest band carries Chern number one. The wavefunction of c-SDW/SCCL is obtained from projection to physical Hilbert space as shown in Eq. (8). The wavefunction of c-SDW/SCCL is obtained from projection to phyical Hilbert space as shown in Eq. (8).
Finally, there is an important subtlety in the PSG construction related to finite samples. Although we explicitly construct a mean-field Ansatz which is invariant under a combination of symmetry and local gauge transformations, it is possible that we can not achieve this consistently on some finite lattice samples with PBC, i.e., having no open boundary. In c-SDW/SCCL case only 4N × 4N lattice sample supports the PSG pattern. However, when π-flux is included in both directions, the resulting Ansatz is symmetric (up to a gauge transformation) in (4N + 2) × (4N + 2) lattice samples. Wavefunctions obtained by π-flux insertion are related to topolog-ically degenerate ground states in thermodynamic limit. We discuss this further in Appendix C.

d+id SC state
Construction of d+id SC state is much simpler. Mean field Ansatz is given in Sec.II. Consider the bosonic holon part first. For the uniform hopping model, bosons will condense at Γ point, and only contribute a constant number after projection. For the fermionic spinon part, the mean field wavefunction is of BCS type: and v(k) are 2 × 2 matrices and ( u(k) v(k) ) are eigenvectors corresponding to positive eigenvalues of H M F d+id,f (k) in Eq. (4). Note that due to vanishing of pairing at the Γ point, |k = 0 = c † k=0,↑ c † k=0,↓ |0 k=0 is not a BCS type wavefunction, and only contributes a constant number (similarly to the bosonic part), so we can omit it in the following analysis.

Appendix C: Understanding quantum numbers
In this section, we use projective wavefunctions to analytically understand quantum numbers of c-SDW/SCCL and d+id SC on different lattice samples. The results are not limited to projective wavefunctions but hold throughout the quantum phase.

c-SDW/SCCL state
We will consider four wavefunctions formed from the considered Ansatz by flux insertion, as they represent the topologically degenerate ground state manifold (the flux is inserted through the handles of the torus formed by the periodic system). To understand quantum numbers for various lattice sizes, it is convenient to use momentum space. The mean-field Ansatz of c-SDW/SCCL already has a doubled unit-cell in x direction, and to make the Brillouin zone more symmetric we double the unit-cell in the other direction too. This enlarged unit-cell contains 8 sites and in this entire Section we will call it the "quadrupled UC" to avoid any confusion (see Fig. 4b).
Thus Brillouin zone becomes a hexagon, and it is simpler to consider C 6 rotation in momentum space.
It turns out that all further calculations are greatly simplified if we immediately insert a π-flux through both directions of every quadrupled UC in the c-SDW/SCCL Ansatz. Then we consider two types of samples analogous to Fig. 5a: The 4N × 4N × 2 = 2N × 2N × 8, to which the 32-site sample belongs; and the (4N + 2) × (4N + 2) × 2 = (2N + 1) × (2N + 1) × 8, to which the 8-site belongs. (Note that the latter family experiences the above π-flux insertion as an insertion through the entire system, and the Ansatz is changed to a topologically degenerate one; for the former family the flux insertion is a simple redefinition of gauge.) All PSG transformations can be performed consistently on all above samples in this redefined Ansatz. We label the state described by the redefined Ansatz as [0,0]. The other three topologically degenerate states are obtained by adding π-flux through entire system in different directions, and the states are labeled as [0,π], [π,0] and [π,π]. In the following, we will analyze the quantum numbers of these four states.
a. [0,0] state Because the quadrupled unit-cell is doubled comparing to unit-cell of mean-field Ansatz, we get double degeneracy for every band. Boson band minimum is moved to Γ point due to the insertion of π-flux through every quadrupled unit-cell. The special property of this [0,0] state is that the mean-field Ansatz is indeed invariant up to a gauge transformation defined by PSG on all 2N × 2N × 2 lattice sizes. Further, gauge transformation G U associated with symmetry operation U turns out to be independent of unit-cell, but only depends on sublattice index. For the other three states, we find that it is impossible to write a consistent mean field Ansatz invariant under all PSG operations (especially the C 6 rotation). In other words, the other three states break (rotation) symmetry explicitly.
In momentum space, PSG transformation is defined as . . , f nk ) t , while n = 8 is number of bands(sublattices). U • k is symmetry transformation for k points while S U (k) is an n × n unitary matrix which represents action of symmetry on sublattice. G U is the associated gauge transformation, with (G U ) ij = δ ij exp(iφ U (i)). Note that in general the gauge transformation of fermions has more freedom, and we can choose a different G U than for bosons. Here, for simplicity, we assume fermions have the same PSG as bosons. The mean-field Hamiltonian is invariant under PSG.
First we analyze the contribution to quantum numbers from fermionic (holon) part. For symmetry U and associated gauge transformation G U , the invariance of holon Hamiltonian can be expressed as Setting α(k) as an eigenvector of H f (k) with eigenvalue λ, we can define α(U It is easy to see that α(U • k) is indeed an eigenvector of H f (U • k) with eigenvalue λ. In this way, we can generate where we assume U m U • k = k, and m U can vary for different k. Note that However, since there is a two-fold degeneracy, it is always possible to choose appropriate α(k) such that Now we apply symmetry on this set of states α(U i−1 • k), i = 0, . . . , m U − 1. By definition, Using the definition of α(U i • k), it is straightforward to derive So under symmetry operation, this set of eigenstates will pick up a θ α(k) phase plus a gauge transformation. It is clear that θ α(k) is directly related to Berry phase of symmetry operation, which is independent of our choice of basis. (To be more precise, this phase is invariant under U (1) phase choice of α(U i • k)). From the above transformation law, it is not hard to get the contribution to quantum numbers from holons. Examples will be presented below. Let us now do a similar analysis on bosonic (spinon) part. For BCS-type Hamiltonian, the invariance of Hamiltonian under PSG transformation can be expressed as Assuming H b (k)β(k) = λβ(k), and using a similar method to above, we can generate β(U • k), . . . , β(U m U −1 • k) as eigenvectors of H b (U • k), . . . , H b (U m U −1 • k) with eigenvalue λ. By appropriately choosing these vectors, it is possible to make In the following, we will show that the additional U (1) phase exp[iθ β(U i •k) ] is unimportant for the BCS-type wavefunction. We only need to focus on G U in the BCStype wavefunction.
Applying symmetry U on M (k) defined in Eq.(B11), we get where Θ 1 (k) and Θ 2 (k) are n × n diagonal matrices, and their elements are additional U (1) phases for different eigenvectors discussed above. Particularly, According to Eq.(B12), Cooper pair creation operator is only picking up a gauge transformation defined by PSG. We can view this as b † ikα → G U (i) * b † ikα under symmetry transformation U , where G U (i) is the i-th diagonal element of G U . Since BCS-type wavefunction is formed by condensation of Cooper pairs, when acted on by symmetry, the only contribution comes from gauge transformation G U . It is worth mentioning that this result also applies to fermionic singlet superconductor, which appears in the case of fermionic spinon in d+id SC.
In the following, we will apply the above results to symmetry group defined in Appendix A. First, let us consider the quantum number of T 1 . Written in momentum space, its gauge transformation can be expressed as a diagonal matrix depending only on sublattice index, while Assuming α(k) is eigenstate of H f (k), then = e iθα(k) G † T1 α(k) (after choosing a convenient α(k)). It is easy to show that θ α(−k) = −θ α(k) . Thus the phase apart from G † T1 will always cancel. The holon wavefunction will transform as where i = 1, 2 for case of one quarter doping. We can view this as f † j → G * T1 (j)f † j . Now we turn to the spinon wavefunction. According to previous analysis, spinon b † ikα picks up phase G * T1 (i) under T 1 . For the total projective wavefunction, we have a constraint on Hilbert space: There is only one spinon or holon per site. Due to this constraint, the total phase obtained from T 1 is simply the product of G * T1 (i) for all lattice sites. So the T 1 quantum number of c-SDW/SCCL is 1.
For translation T 2 , we do a similar procedure as for T 1 , and find that the quantum number of T 2 also equals 1. So, we can conclude that the center of mass of [0,0] state is at Γ point for 2N × 2N × 2 lattice size, i.e., for both sample families introduced in this Section.
Let us turn to C 6 symmetry. It is straightforward to get the sublattice transformation matrix: while for bosonic spinon, where β = π/6. Note that although G b,C6 is also a consistent gauge transformation for fermion, we choose G f,C6 different from G b,C6 for simplicity. We have three classes of k points in Brillouin zone according to their transformation rule under C 6 : 1) Γ point, which transforms back to itself under C 6 , so m C6 = 1; 2) Three M points, which transform back to themselves under C 3 6 (inversion), so m C6 = 3; 3) Other k points, which are invariant only under C 6 6 , so m C6 = 6. Using the method developed above, we calculate the additional U (1) phase under C 6 for the 1st and 2nd holon band (in quarter doped case, holons always fill these 2 bands). The result is listed below: We checked this numerically for various mean-field parameter values.
It is easy to see that only Γ point contributes additional phase, which equals e iπ/3 . Under C 6 symmetry, holon wavefuncion transforms as For spinon part, the transformation law is We can view this as if every spinon picks up factor β = e iπ/6 (plus fermion gauge transformation) under C 6 . It is straightforward to calculate that for 4N × 4N × 2 lattice size C 6 quantum number equals e iπ/3 , while for (4N + 2) × (4N + 2) × 2 lattice size C 6 quantum number is e −2iπ/3 . For the state related by time reversal, quantum numbers are obtained by conjugation.

b. Other three states
Using the method developed above, we calculated translation quantum numbers of the three other states. It turns out that the center of mass of these three states are three M points ((0,π), (π,0) and (π,π)). While calculation details will not be presented in this paper, there is a simple physics picture. Consider adding π-flux in x direction to the [0,0] state, and then translating in the same direction. This corresponds to every fermion hopping one lattice spacing in x direction, and they will see this additional π-flux. Thus, compared to original state, the translation quantum number in x direction is multiplied by −1, so the center of mass will change from Γ point to M point [π, 0].
For rotation, we note that C 6 is not a symmetry for these three states. But the three states are symmetric under inversion symmetry C 3 6 . Applying the above method, we find that these three states have opposite inversion quantum number to [0,0] state, which is consistent with our field theory analysis in Section IV B.

d+id SC state
Understanding the quantum numbers of d+id SC is much simpler. Firstly, holons always condense at Γ point, and contribute an overall constant, thus can be neglected. Secondly, two spinons that occupy Γ point will also have no contribution, as discussed in Appendix B. For other spinons, which have a BCS-type wavefunction, the analysis of quantum numbers is similar to bosonic spinon part above: Under lattice symmetry, only gauge transformations contribute to quantum numbers.
For translation T 1 and T 2 , associated gauge transformations G T1 and G T2 are trivial. So, the center of mass is Γ for any lattice size.
Under C 6 rotation, mean-field wavefunction changes as We can view this as if every fermion picked up e −iπ/3 after C 6 (except for fermions at Γ point). So C 6 quantum number for 2N × 2N × 2 lattice size is independent of lattice size. Next consider the inversion C 3 6 quantum number. For d+id SC, it is always 1. For c-SDW/SCCL, on 4N ×4N × 2 lattice size, inversion quantum number equals -1, while on (4N + 2) × (4N + 2) × 2, inversion quantum number is 1. This provides a sharp signature to distinguish c-SDW/SCCL state and d+id SC in finite samples. To consider spin rotation symmetry in x and y directions, one must enlarge the K matrix by adding degrees of freedom that are in a topologically trivial phase 0 1 1 0 . Then we get Any such transformation by a matrix X in GL(5, Z), the group of unimodular N × N matrices, can be seen as a relabeling of topological degrees of freedom since l → l = X t l, and the physics remains unchanged. However, the particular choice of X, inspired by Ref. 74, allows an easier identification of physical properties. The charge vector t c = (1, 0, 0, 0, 0) and S z vector t Sz = (1/2, −1, 0, 0, 0) are direct extensions of the original ones. We identify holon f as (1, 0, 0, 0, 0), spinon b ↑ as (0, 0, 0, 1, 0) and spinon b ↓ as (0, 0, 0, 0, 1). Electron e ↑(↓) is simply the bound state of f and b ↑(↓) , and it is in the topologically trivial sector. Vison can be viewed as 'half holon', and is expressed as (0, 1, 0, 0, 0) with charge 1 2 and (0, 0, 1, 0, 0) with charge − 1 2 . It is easy to check that the statistical angles and quantum numbers of these quasiparticles are correct.
The general consideration of symmetry in this Kmatrix formulation has been considered in Ref. 74. Under symmetry g ∈ G s , chiral boson field φ I will transform as Notice that the above symmetry transformations {W g , δφ g |g ∈ G s } must be compatible with group structure of symmetry group G s . More precisely, the local bosonic degree of freedom {M I ≡ e il I J K I,J φ J } must form a linear representation of symmetry group G s while nonlocal quasiparticles can transform projectively.
It is not yet known how to incorporate the full SU (2) spin rotation symmetry in the K-matrix formulation. However, we can choose a subgroup of SU (2), which is generated by g 1 , the π rotation around S x direction, and rotations around S z direction, U θ , θ ∈ [0, 4π). They satisfy the following algebra: In fact, we can view this group as a projective representation of a SO(2) z Z 2 subgroup of SO(3).
Following Ref. 74, we find a consistent solution for {W g , δφ g |g ∈ G s } that describes SCCL, namely Explicitly, the quasiparticles transform as Notice that φ 2 and φ 1 − φ 2 only differ by a trivial boson, so they are actually the same quasiparticle with the same quantum numbers. Further, each spinon (φ 4,5 ) acquires −1 Berry phase after 2π-spin-rotation, while holon φ 1 and vison φ 2,3 transform trivially, as expected.
There are several Higgs terms allowed on a symmetric boundary by this transformation law. However, one should consider the largest subset such that all terms can condense simultaneously, meaning that the arguments of these terms commute. Furthermore, the condensed fields must not break the symmetry. Thus we arrive at the following Higgs terms: Define θ 1 = φ 1 + 2φ 3 , θ 2 = φ 1 − 2φ 2 and their conjugate variables are ϕ 1 = φ 4 + φ 5 , ϕ 2 = φ 4 − φ 5 . It is easy to show that {θ 1 , ϕ 1 }, {θ 2 , ϕ 2 } form two decoupled Luttinger liquids, so they can be gapped by the Higgs term L Higgs = C 1 cos θ 1 + C 2 cos θ 2 . The only remaining gapless degree of freedom is φ 1 . So the edge of SCCL is chiral fermion liquid of holons.
Appendix E: Tunneling conductance calculation for different junctions

Point junctions
For completeness we repeat the metal lead/SCCL case from the main text here.

Tunneling Hamiltonian is
For Fermi liquid systems, the scaling dimension δ F L = 1/2 in any dimension. 1 So, δ = 2δ F L = 1. Using Eq.(30), we get that G(T ) is constant in this case.

• Metal and SCCL
Due to the spin gap on the boundary of SCCL, single electron will decay exponentially when tunneling to edge of SCCL. So, the major contribution is from singlet pair tunneling.

• SC and c-SDW
The tunneling Hamiltonian is where Cooper pair operatorĉ is a complex number inside SC. Therefore δ = 2 in this case, and we get G(T ) ∼ T 2 .
• SC and SCCL Since singlet Cooper pairs are not influenced by spin gap, the result should be the same as for SC and c-SDW, namely G(T ) ∼ T 2 . We next present the perturbative Fermi golden rule calculations for various point junctions.
The tunneling current is where So, ρ 0 (V ) and ρ M (V ) are constant numbers. We get I(V ) ∼ V 0 ρ 0 (V )ρ M (V ) dV ∼ V , and tunneling conductance G(V ) = dI/dV is constant.
For SC lead, the main contribution is from tunneling of singlet Cooper pairs. Therefore, G(V ) scales in the same way for c-SDW and SCCL, and we get G(V ) ∼ V 2 . Comparing with the above results, the perturbative calculation is indeed consistent with simple RG analysis.

Line junction
Voltage on the metal/SC side is a constant number, labeled by V R (Fig. 17). Electron scattered from lead will lose its phase and always keep at the same voltage. On the c-SDW/SCCL side, voltage is maintained between scattering events and is accumulated, as shown in Fig. 17.
We first completely derive the case of junctions for which the point contact scaling exponent α = 0, and deal with the α = 0 case (c-SDW/metallic lead) at the end.
The voltage at n-th point junction is labeled by V n , while the tunneling current is I n . Due to anomalous quantum Hall response of electron/holon, we get V n − V n−1 = I n e 2 /h .
According to the result for point contact having α = 0, in the regime T V , where α is the scaling exponent for tunneling conductance obtained in point junction case. Eq.(E9) can also be viewed as definition of T (n) K . Define x n = V R − V n , to get which we can transform into a differential equation: Integrating the above equation from the initial x 0 to the final x N yields in which we defined the effective T K from the individual T After integration, one obtains where we define V ≡ x 0 = V R − V 0 . The total current flowing from metal/SC to c-SDW/SCCL is obtained from the voltage difference V N − V 0 , and is given by so the tunneling conductance is The result expressed holds for all values of V, T K at T → 0, as long as the assumptions used to derive the expression holds, namely, each individual point contact junction is weakly coupled, V T (n) K for all n. Note that the effective T K for a long line junction (large N ) can be very small compared to all T (n) K . Now, let us consider the small voltage regime, namely, V much smaller than temperature T . However, we still require the weak coupling condition for single point junctions, namely, T T Following similar steps as above, we get By solving this equation, it is straightforward to get the tunneling conductance as a function of T : where we define (E19) Finally we consider the c-SDW/metallic lead line junction, i.e., the case of α = 0. The derivation procedure is the same as for the above case, and starts from the point junction result: K , which physically corresponds to low enough temperatures and voltages. Using Eqs. (E8), (E20), and the same procedure as above, we get: where we defined Here we discuss the precise DMRG setup, convergence to true ground state with limiting MPS matrix size, and also present some measurements for parameter values not shown explicitly in the main text.
To represent the two-dimensional periodic samples in the DMRG in a way that eases convergence, we labeled the sites 1 . . . N such that the longest necessary hopping range is minimized. For present samples, which have aspect ratio of 1, it was sufficient to sequentially order site labels 1 . . . N from, say, left to right within each row and then from one row to the next. With larger twodimensional samples in lattices with higher coordination, it is advantageous to avoid labeling rows sequentially, but instead, starting from one row, sequence the one below it, then the one above it, and so on in an alternating fashion. We have checked for some parameter values that the labeling which minimizes the longest range hopping indeed allows faster sweeps and better convergence in the same amount of time. FIG. 19. Spin-spin correlation function Sz(i)Sz(j) in the DMRG ground state projected into sector with center of mass momentum Γ and C6 rotation eigenvalue exp(−iπ/3). Blue is positive, red negative, and disk radius is proportional to amplitude. Site i is fixed at green circle, while every bond i, j is averaged over translations and rotations to reduce statistical error by increasing the number of sampled observable values in MC. "Max" labels absolute amplitude of largest shown disk. All four parameters are chosen within the c-SDW/SCCL phase: (a),(b) Hubbard model; (c),(d) t-J model.
The convergence of DMRG energy is however limited in practice by the maximal size of MPS matrices, m, which does not exceed 11.000 in our calculations. In Fig. 18 we present a typical convergence of DMRG energy as a function of 1/m, with a linear fit extrapolation towards infinite m. This is not the common way of considering DMRG convergence, but it is informative given our m value limitations.
As discussed in Section III C, in this paper we quantify the DMRG convergence to the true ground state by using the expectation value of C 6 symmetry operation (60 • rotation), which should be one of {−1, 1/2} ({1, −1/2}) when the inversion is −1 (+1). (The inversion is always numerically very precisely quantized.) As shown in Fig. 9, the C 6 measurement on the 32-site sample indicates the convergence failure in the J/t > 0.8 phase; in the Hubbard model, the convergence progressively worsens with lowering U/t below the value 5. It is not surprising that convergence worsens for low U/t, but we believe it is mainly due to our m limitation. For instance, at U/t = 4 the C 6 expectation with projection to Γ momentum improves from 0.10(2) at m = 8.000 to 0. 16(2) at m = 10.500.
Next, we present additional details about correlation functions on the 32-site sample.
In Fig. 19 we show the spin-spin correlation function in the c-SDW/SCCL phase for several values of parameters, as addition to Fig. 11. The values are chosen to demonstrate how the longer-range spin correlations match the tetrahedral pattern even better as U/t grows and as J/t decreases. On the other hand, the magnitude of shortrange spin correlations grows with both U/t and J/t as expected. Let us here emphasize again that we use the total S z equal to zero sector in both models throughout this paper. The DMRG calculation conserves this quantum spin number of a state, as well as the total number of fermions.
In Section III C we claimed that the pair-pair correlation function on 32-site sample in the c-SDW/SCCL phase is very short ranged. Here we provide a numerical example to compare to 24-site sample results in Table II. On 32-site sample we consider the J/t = 0.78 DMRG GS projected to exp(−iπ/3) eigenspace of C 6 , and pairs of nearest-neighbor bonds separated exactly as in figure of Table II. Every correlation value is obtained using 64 MC runs of 10.000 measurements, and averaged over translations of the bond pair to additionally reduce statistical error. (The usual 500 measurements give a statistical error that overwhelms the value of correlations.) Out of the three bond pairs, the maximal correlation magnitude is 0.00059(2), to be compared with 0.00163(2) and 0.00447(2), the values for J/t = 0.78 and J/t = 2.0, respectively, for 24-site sample from Table II. The complex phases of the three bond-pair correlations in d+id state are 0, 1/3, −1/3 in units of 2π, but in the considered 32site measurement we find 0.5(1), 0.40(5), −0.45 (4).
Finally, we explained in Section III C that on the 32site sample the DMRG GS did not converge well in the large-J phase J/t > 0.8, so the correlation measurements are not trustworthy, but we note for completeness that in that regime the obtained DMRG GS with projection to Γ center of mass momentum and either exp(−iπ/3) or exp(−i2π/3) eigenvalue of C 6 , completely loses resemblance to tetrahedral spin pattern without developing a d+id pair-pair correlation pattern.