Stable Flat Bands, Topology, and Superconductivity of Magic Honeycomb Network

Jongjun M. Lee, Chenhua Geng, Jae Whan Park, Masaki Oshikawa, Sung-Sik Lee, 5 Han Woong Yeom, 1 and Gil Young Cho ∗ Department of Physics, Pohang University of Science and Technology (POSTECH), Pohang 37673, Republic of Korea Institute for Solid State Physics, The University of Tokyo, Kashiwa, Chiba 277-8581, Japan Center for Artificial Low Dimensional Electronic Systems, Institute for Basic Science (IBS), Pohang 37673, Korea Department of Physics & Astronomy, McMaster University, 1280 Main St. W., Hamilton ON L85 4M1, Canada Perimeter Institute for Theoretical Physics, 31 Caroline ST. N., Waterloo ON N2L 2Y5, Canada (Dated: August 6, 2019)

1. Introduction: Topology of electron wavefunctions and strong correlations are the two main resources for modern condensed matter physics to develop and grow. [1][2][3] The strong electronic correlation is the hidden root of the central enigma in the field, namely the microscopic coexistence of competing orders. In such a material, spin/charge orders and superconductivity are coherently intertwined 3 to form a self-organized superstructure, which is against the conventional wisdom. On the other hand, the discovery of the non-trivial topology in electronic states has drastically changed our view on quantum states and left huge impacts on various fields of physics. 4,5 The two fields have dramatically evolved when a new theoretical/experimental platform, which can show the desired phenomena, has been discovered. Recently, the novel electronic states emerging from moire superstructures have received huge attention as a potential platform for realizing surprising varieties of topology and correlation-driven phenomena. [6][7][8][9][10][11][12][13][14] In this manuscript, we uncover that yet another superstructure electronics, namely a honeycomb network of one-dimensional metals, is also an ideal system to explore both the physics of the diverse topology and many-body phenomena and that it may be a hidden driver of the emerging superconductivity in charge-density wave materials [15][16][17][18][19] such as 1T-TaS 2 , which we focus here. However, the theory will be generically applicable to the broader class of materials supporting the network superstructure such as 3 He and 4 He mixtures absorbed to graphite, 20 several chemical compounds, 21-25 twin boudary networks in MoSe 2 grown on MoS 2 , 26 and artificially-built network. 27 Experimentally, 1T-TaS 2 is known to have a rich phase diagram of charge-density wave (CDW) orders 28 and su-perconductivity (SC), [15][16][17] which are accesible by tuning temperature, pressure, and doping. Depending on the degree of the charge ordering, there are three distinct regimes: commensurate CDW (C-CDW), nearlycommensurate CDW (NC-CDW), and incommensurate CDW (IC-CDW) . The C-CDW phase has a long-ranged charge ordering and appears at the lowest temperature with the ambient pressure. This state is the correlationdriven Mott insulator, 29 which has been suspected to be the long-sought spin liquid. 30,31 When the C-CDW ordering is slightly suppressed by pressure or doping, then the domain walls appear inbetween the locally chargeordered domains, i.e., it enters into the NC-CDW state. If the pressure or doping increase further, the SC emerges from this NC-CDW state. Note that a similar phase diagram is obtained in another CDW material, namely 1T-TiSe 2 . 18 This suggests that the NC-CDW [15][16][17][18][19] is somehow essential for realizing SC in these CDW materials though how this actually happens has been unclear.
We will point out that this NC-CDW state and its associated honeycomb network of metallic domain walls 32-34 host a series of robust flat bands, which can greatly help for the SC to emerge. A recent STM study 33 combined with ab-initio DFT calculation provided an unprecedented detail of the electronic structure of the network, where the metallic nature of the domain walls and tri-junctions in the network is clearly exposed. Motivated by the experiments, 32-34 we consider the emergent electronic states out of this network superstructure, which has been missing so far, and thoroughly examine the stability of the flat bands. We note that the emergence of entirely new electronic states from the superstructure is strongly reminiscent of the Moire physics 10 of arXiv:1907.00012v2 [cond-mat.str-el] 5 Aug 2019 twisted bilayer graphene. Basing on the robustness of the flat bands, we consider the correlation-driven many-body phenomena. The effect of the correlation is particularly drastic when the Fermi level is near the flat bands. In this case, we show that the network can strongly enhance the superconducting T c . Hence we clarify the long-sought cooperative mechanism between the NC-CDW state and the emerging SC in 1T-TaS 2 . On top of this, we show that the honeycomb network is an outstanding system to realize a long list of topological band structures including two-component Dirac/quadratic band touching, 35 threecomponent spin-1 Dirac fermion, 36 Chern insulator and zero-dimensional corner state 37-40 which is closely related to that of higher-order topology.
2. Model and Flat Bands: We first present a simple tight-binding model which exposes the essential physics. It consists of the modes inside the CDW gap, which are entirely from the domain wall regions. See Fig 1 (A).
where the second sum over {r, r } ∈ J runs over the sites around the junctions J and t J ≤ t 0 . Here σ represents the spin and the sum over it is assumed. Here we have assumed the spin rotational symmetry for simplicity. The model also has time-reversal symmetry T , reflection R x , and six-fold rotation C 6 , which are the symmetries observed in the STM experiment. 33 Diagonalizing Eq.(1), we find a cascade of flat bands, Dirac and quadratic band crossings. See Fig 1 (B). On top of this, a three-component spin-1 Dirac fermion 36 can appear when t J /t 0 → 0. The number of the flat bands and topological band crossings is proportional to the number of the sites between the junctions. In 1T-TaS 2 , 33 we expect that the energy difference between the neighboring flat bands is roughly 30 meV and there are at least 4-5 flat bands (or more) inside the Mott gap. The topological band crossings are protected by crystalline and Tsymmetries. On the other hand, the flat dispersion cannot be generally protected because it requires infinitely many parameters to be tuned. Hence, they are generically fragile, e.g., against the second neighbor hoppings. 41 Despite of the fragile nature, flat bands are an ideal stage to realize correlation-dominated physics, such as ferromagnetism, and thus have been studied vigorously. 21 A well known general mechanism which gives rise to flat bands is an imbalance between sublattice sites in bipartite lattices. This was also applied 42 to hyperhoneycomb systems which have some similarities with the systems we study in this paper. However, such flat bands can be generically removed by inclusion of short-ranged further neighbor hoppings, which exists in real materials.
Remarkably, our flat bands from the network defy this standard phenomenology and are stable against D 6symmetric local perturbations. For example, addition of the third neighbor hoppings t 3 (or even the fourth neighbor hoppings) near the nodes does not disperse the flat bands. Note that t J in Eq.(1) is already the second neighbor hoppings. We can even include the insulating electronic states from the domain area (described by H dom ). See Fig.1

(A). The full Hamiltonian is now
where c r,σ (d r,σ ) electrons belong to the network (to the domains). Here H dom is described by the band insulator which has the two energy-splitted states per site and the different sites are connected by the small hopping t D [ Fig. 1 (c)]. From the perturbative reasoning, we expect that this model includes all the possible symmetric local perturbations. Notably, the flat bands and overall band shapes remain almost intact inside the insulating gap [ Fig.1 (C)]. This result implies that in the NC-CDW 1T-TaS 2 , even if the bulk bands from the domain region rise up, flat bands inside the gap are almost intact. Finally, we checked that the flat bands survive 41 against the Rashba spin-orbit coupilng. Such stability is unique to the honeycomb network and is absent for other network such as triangular and square networks. 41 The above finding is consistent with our previous study. 33 There, we have considered a single-particle scattering problem of the same network where an electron propagates ballistically along 1d wires and scatters at the junctions. This description reflects faithfully the essential features of the band structure of Eq.(1), i.e., the repeated appearance of the flat bands and topological band crossings. Within this description, independent of the specific values of the parameters, we always find the flat bands (the summary is in SI. 41 ). This unusual stability of the flat bands requests some explanation. For this, we look carefully into the structure of the wavefunctions inside the flat bands. We find that those wavefunctions vanish at the junctions. 41 This means that when the low-energy degrees of freedom are entirely from the network and only the nearest neighbor hoppings are included, i.e., t J = 0 in Eq.(1), the wavefunction is a standing wave ψ(l), e.g., ψ(l) ∝ sin( πl L ) [ Fig  2 (A)]. From such standing waves, we can construct a set of localized states 43 which consist of the flat bands: we consider a linear combination around the honeycomb plaquette with a sign oscillation, i.e., Ψ ∼ a (−1) a ψ(l) for a ∈ {1, 2, · · · 6} labeling the six links around the plaquette. See the pictorial representation in Fig. 2 (A). It is straightforward to see that this state cannot disperse into the neighboring plaquettes because of the destructive interference. See Fig 2 (A) and its caption. Such destructive interference persists as far as the hopping distance is shorter than the length of the wire. We can also understand that the cascade of the flat bands must appear because there are many standing wave solutions per wire. This illustrates the importance of the symmetry and the locality on the stability of the flat bands. By passing, we note that the energy of the flat bands is related to the energy of a single wire subject under the open boundary condition. 44 From the above, we expect the flat bands to be removed either by long-ranged direct hopping processes across the CDW domains or by breaking symmetries. Indeed, we can show that the flat bands are lifted when the symmetries are broken [Fig 2 (C) and SI 41 ], and/or when such long-ranged hoppings [e.g., t f term in Fig 2  (A)] are included [Fig 2 (B)]. This means that, for the NC-CDW state of 1T-TaS 2 , we need a sizable hopping between the two sites apart by ∼ O(80)Å 33 to remove the flat bands, which is not realistic. This explains why the dispersion of the "flat bands" are very small but still finite when the domain sites are included(the "dispersion" is induced by the exponentially suppressed but nonzero hoppings through insulating domain sites.). We also note that when the time-reversal symmetry is broken, the band touchings are gapped out and it results in dispersive Chern bands [Fig 2 (C)]. Here, the resulting band width can be still small and thus we find a long-desired relatively-flat Chern band, which is a natural platform for fractional Chern physics. 45 Similarly, we can straightforwardly generalize this to the quantum spin Hall bands.
We comment on the effect 46-49 of the interlayer coupling in 1T-TaS 2 . The interlayer interaction has been suggested to be important in metallization and SC in 1T-TaS 2 . However, we remark that there are plenty of experimental data and theory 31,33 suggesting that the main physics is essentially two-dimensional. For example, the resistivity along c-axis is much larger, e.g., by 500 times, 50 than the intralayer resistivity 50,51 and anisotropic two-dimensional charge transfer is observed for the NC-CDW state. 52 Further, the SC T c is almost insensitive to the pressure 15 and not much affected under the dimensional reduction. 16 Basing on these, we focus on the two-dimensional physics here.
3. Superconducting States: Having established the stability of the flat bands, we now discuss the many-body physics when the Fermi level is near one of the flat bands. In general, such system is unstable toward various particle-hole (e.g., ferromagnetism 53 ) and particleparticle channels. However, guided by experiments, we mainly focus on SC in this paper.
First, we perform the simplest BCS mean-field theory with the phenomenological interaction Projecting to the BCS channel, we obtain: 41 where g l is the interaction strength along the pairing channel of the angular momentum l, which we compute numerically.ˆ l;p is the corresponding pairing order parameter. Note that, within the mean-field decomposition of Eq.(2), only the spin-singlet sector appears. Below we consider only the s-wave and (d + id)-wave pairing channels, i.e., g 0 and g 2 = g −2 of Eq.(3). We find that typically {g 0 , g 2 } are positive when U is large and positive. However, as the negative V is added, g 0 and g 2 decrease toward negative and signal the instability toward the superconducting states. As the negative V increases, generically g 2 first becomes negative and thus it opens up a window for the (d ± id)-channel pairing. If U and V are both negative, then typically g 0 is more negative than g 2 and the s-wave pairing is favored. The phase diagram including the ferromagnetism (which appears when the interactions are dominantly repulsive) is in SI. 41 The magic of the flat bands appears when the gap equation is solved. The standing waves receive ±1 around the plaquette Ψa(l) ∼ (−1) a ψ(l). Here the red colored waves are with (−1) sign and the blue colored waves are with (+1) sign. When the state attempts to hop to the next plaquette, there are two contributions: one from the blue one t × (−ψ(l)) and the other from the red one t × (+ψ(l)). When summed up, they exactly cancel. (B) Typical band structure when the hopping across the domain is added. For this process, there is no destructive interference and so the flat bands become dispersive. (C) Band structure with the broken T -symmetry. In this case, all the band touchings are removed and the flat bands become dispersive. Furthermore, all the bands carry the non-zero Chern numbers.
where F l (k) is the form factor for each pairing channels, 41 e.g., F 2 (k) → (k x + ik y ) 2 for |k| 1, and k = µ. Hence, we find that the SC gap is linearly enhanced, i.e., | l | ∼ |g l | (if g l < 0), instead of the standard exponential suppression ∼ exp(−1/|g l |ν l ). Because there is no other scale, the mean-field energy of the SC is linear in g l and so is the T c . 41 Hence, the honeycomb network strongly enhances the SC T c .
In the s-wave SC, the BdG fermion is fully gapped as usual. On the other hand, the (d ± id)-wave SC is not fully gapped and exhibits the doubled quadratic band touching. 41 Since this state is gapless, it does not support any topologically-protected edge mode. Because the quadratic band touching is marginally unstable 35 against the short-ranged repulsive interactions toward the chiral or nematic states, there will be a successive transition at the temperatures below the SC T c .
On the other hand, it also instructive to consider the strong-coupling limit of the network, 54 where the interaction is bigger than the hopping integrals. For this, we start from the decoupled strongly-correlated wires, each of which may be described by a Tomonaga-Luttinger liquid (TLL) where the Luttinger parameters {K c , K s } capture the correlation non-perturbatively. 2 This is the correlated version of the scattering problem in our previous work, 33 which faithfully represented the band structures.
The two-dimensional SC can be preempted by the spin gap 2,55 in each wire. Once the spin gap forms, the low-energy physics of each wire is described by a single-component TLL of θ a,c describing the fluctuating SC pairs. 2 Since the SC pair is bosonic, the junction of three TLLs at each vertex of the honeycomb network corresponds to the bosonic Y-junction, 56 rather than the fermionic one 57,58 where the fermion statistics plays an important role. When each wire has sufficient attraction, i.e., K c ≤ 1, 2 then the interwire coupling between the SC fluctuaions, namely the Josephson coupling J, becomes relevant. 2 This will lock the spatial pattern of the phases of the SC fluctuations between different wires. Only interested in the pattern of the phases, we simply note that the problem is symmetrically equivalent to the XY model on the Kagome lattice, while leaving the full quasi-1d interchain mean-field treatment 55 to the future study Depending on the sign of J, either the two-dimensional s-wave or (d ± id)-wave SCs can emerge from these fluctuations. When J > 0, then the so-called √ 3 × √ 3 order appears, 59 which translates as the (d ± id)-SC. If J < 0, the conventional s-wave SC emerges. When the repulsive-U dominates the junction region and the region becomes Mott insulating, 60 J > 0 can appear. From this strong-coupling limit, we can learn how the 2k F -density wave of 1d wires competing with SC is suppressed. For generic filling of each wire, the momentum k F will not be commensurate with the wire length L, i.e., φ L = k F L is not a rational number. This will frustrate the phases of the density waves, and thus their true two-dimensional order is strongly suppressed to develop. This gives a natural favor on the SC over insulators.
When the filling per wire is commensurate, then the network can develope a charge density wave order and become an insulator. This can lead to the emergent corner states at the tri-junctions and corners of the sample boundary, which are akin to those of higher-order topological insulators. 37-39 Motivated by references 61,62 where the nature of the insulating domain walls of C-CDW 1T-TaS 2 were discussed, we consider the half-filled wires which can develop a period-2 charge density modulation. When the tri-junction frustrates the charge order and traps solitons, we can show that there is an associated crystal symmetry-protected bound state at the junction (In SI, 41 we provide a few concrete realizations.). By passing, we note that gapped domain wall junctions in the C-CDW state 1T-TaS 2 61 and similar junctions in twin boudary networks of MoSe 2 grown on MoS 2 26 are found experimentally although the presence of higherorder topology in those networks is not clear at this moment.
4. Conclusions: In this Letter, we have considered the electronic structure of a conducting honeycomb network. We have constructed a lattice model and uncovered the emergence of the cascade of the flat bands. Combining the scattering theory, 33 tight-binding calculation and also the analytical reasoning, we have established and explained the emergence of many flat bands stable against various local perturbations. Compared to the previous studies, the discovery and analysis of many stable flat bands are unprecedented. Not only the flat bands, we also demonstrated that the honeycomb network is an ideal place to find diverse topological band structures: the regular Dirac and spin-1 Dirac fermions, quadratic band touchings, Chern insulators, topological band insulators, and also zero-dimensional corner states.
When the chemical potential is tuned near one of the flat bands, the many-body effects are enhanced. Motivated by CDW materials with coexisting SC, we have concentrated on the SC states and have discovered the anomalous enhancement of SC T c . In both the weakcoupling mean-field and strong-coupling Luttinger liquid limits, we found that the nodal (d ± id)-SC can emerge when the repulsive-U plays significant role, e.g., suppressing s-wave SC to emerge. When the on-site U is attractive, the system enters into the s-wave state. Given that the domain wall states of 1T-TaS 2 presumably experience small correlation effect and that the junction regions are quite metallic from the STM data, 33 our result suggests that the SC state of 1T-TaS 2 is likely a regular s-wave SC (with some p-wave component due to the spin-orbit coupling). We remark that the strong enhancement of the SC by the NC-CDW network discovered in this paper provides a new example of the so-called "intertwinement" of the orders. This is clearly beyond the previous bosonic analysis, e.g., Landau-Ginzburg theory 19 or nonlinear σ models, which are blind to these emergent electronics from the network superstructure. Also our work here clarifies the precise role of the domain wall network of the CDW materials, for example 1T-TaS 2 , in realizing SC and explains how the two seemingly-competing orders coexist and self-organize. It is also intriguing to note the similarity to the charge-density wave network and moire twisted-bilayer-grphane physics. In both cases, understanding low-energy electronic states, which have very different characteristics from microscopic states, are critical to decipher the origin of the low-temperature phases.
The signature of the network and its role in SC can be detected in various experiments. First, inside the normal state, photo-emission experiment can in principle access the emergent band spectrum. Within the currently-available experimental data, 63-65 we find that it is unclear either to confirm or to rule out the existence of the cascade of the flat bands in 1T-TaS 2 . At least one nearly flat band is observed right below the chemical potential. This calls for further investigation. Second, magneto-transport and oscillations in resistivity of the normal state and SC state can provide the information about the primary conducting and SC channels. They have been applied in the small twist-angle bilayer graphene 66,67 and textured superconducting states of 1T-TiSe 2 18,19 to identify the network geometry. Also the flat bands show a few characteristic behaviors in thermodynamic quantities, which we summarize in SI. 41 -Note Added : After the completion of this work, we became aware of independent work on decorated star lattices 68 where flat bands and higher-order topology are also discussed. Here we summarize the theory part of our previous work 1 . Specifically, we will introduce the scattering description of the honeycomb network, which reproduces the key structure of the tight-binding model in the main text.
The STM experiment 1 essentially found that the domains of the nearly-commensurate charge-density wave form a regular honeycomb lattice, and thus the domain walls are the links of this regular honeycomb lattice. Furthermore, the domain walls trap finite local density of states near the Fermi level. We note that the domain walls are generically expected to trap some in-gap modes due to the topological solitonic modes (though these modes may appear away from the Fermi level.). Motivated from these findings, we consider a regular array of one-dimensional metals living on the links of a honeycomb lattice. Similar network models of one-dimensional metals have been studied in the context of quantum Hall plateau transition, known as "Chalker-Coddington model" 2 , and also in the twisted bilayer graphene at a small twisting angle 3 .

Model
To capture the physics of the network, we introduce the two wavefunctions on the links of the honeycomb network: ψ a and ψā. Here ψ a represent the chiral mode propagating from an A-sublattice (of the network) to a B-sublattice (of the network) and ψā for the mode propagating from a B-sublattice to an A-sublattice. See Fig. 1.
Hence, we can associate ψ a to the node at an A-sublattice and ψā to the node at a B-sublattice, i.e., ψ a is an out-going mode along the link a = x, y, z from an A-sublattice, and ψā is an out-going mode along the link a = x, y, z from a B-sublattice (See Fig. 1).
We further assume that these modes propagate ballistically within each link and scatter only at the nodes of the honeycomb lattice. We further assume that there are six-fold rotation C 6 , mirror R x , and T symmetries, and the scattering between the modes respects the symmetries. (As in the main text, we assume SU(2) spin rotational symmetry and suppress the spin indices here.) With all of these in hand, we can write down the scattering problem at an A-sublattice.
Here, the left-hand side ψ a (R), a = x, y, z represents the out-going modes from the A-sublattice, which is related by a scattering matrixT A to the in-coming modes ψā(R), a = x, y, z appearing on the righ-hand side (See Fig. 1). The additional phase factor ∼ exp(−i E v F L) is the phase accumulated by the incoming modes while it propagates coherently from the neighboring B-sublattices to the A-sublattice at R. Here v F is the Fermi velocity within the one-dimensional metal, which is expected to be similar to that of the bulk electron, and L is the length of the link. The scattering matrixT A is fixed by the combination of the crystal symmetries and T -symmetry. With the unitarity of the scattering matrix, we find . Similary we have the following scattering problem at the B-sublattice.  whereT B =T A by the crystal symmetries. Now we can perform the Fourier transformation and solve these scattering problems. On performing the Fourier transformation, we find whereV q = diag [exp(iq ·ê x ), exp(iq ·ê y ), exp(iq ·ê z )]. Hence, the energy spectrum can be obtained by diagonalizinĝ T q , which is unitary. In terms of the eigenvalues e i j (q) , j = 1, 2, · · · 6 ofT q , we have Here n ∈ Z and thus the minibands are repeating in the energy in period of 2π v F L . Mathematically this ambiguity in n originates from the ambiguity of j (q) by 2π appearing in the eigenvalues e i j (q) , j = 1, 2, · · · 6. Physically this repetition can be traced back to the excitation energy of the microscopic one-dimensional modes with the same momentum q, i.e., for a given q, there are different one-dimensional modes with energy 2π v F L n, n ∈ Z. Thus we expect that this repetition will fill up within the band width of the original parent 1d band. Indeed, this is the band structure that we find from the tight-binding problem, where the certain unit structure repeats in energy.

Spectrum of Network
Next we analyze the band structure out of this scattering description. As apparent from the Fig 2, the spectrum features (i) Dirac cones at K and K , (ii) flat bands, and (iii) quadratic band touchings at Γ-point, which are the features of the tight-binding band structure. Now the crucial question is if these features are stable against the symmetric deformation of the scattering matrices. We see that the parameters that we can tune are {t A = t B = t, v F , χ}, which determine the scattering amplitudes at the nodes and the phase accumulated by the modes while they travel along the links. We find that the overall band structures remain the same. In particular, the flat bands always survive. See Fig. 2. momentum space, where ψ q (r) is the Bloch function of the flat bands at the momentum q and the site r in a unit cell. From Fig.3, we can clearly see that the Bloch function has zeros at the tri-junctions upto the numerical error. This confirms that the flat band states can be generated by the standing waves living in each wires. In fact, we can do better: from the Bloch state, we can even read off the sign structures of the non-dispersing states in the Fig 2 (A) of the main text. The Bloch states have a staggered ±1 signs around the nodes, which is precisely the same structure as the wavefunction of the main text.

Proof of Existence of A Single Flat Band
For a few fine-tuned cases, we can in fact prove the existence of a single flat band. The case is that, there are odd number of sites between the tri-junctions and the Hamiltonian has only the nearest neighbor hoppings. For example, for the case with a single site inbetween the tri-junctions Fig. 4, we can factorize the Hamiltonian into the following: Since the Hamiltonian as a matrix has the rank less than its dimension, it must have an eigenvalue 0 for all the q. This is the flat band because the corresponding eigenvalue is a constant zero over all the momentum q. In fact, such factorization can be easily generalized into the cases with the odd number of sites between the tri-junctions. Suppose that (2m + 1)-sites exist between the tri-junctions, where m is an integer. Hence, there will be 6m + 5 sites per unit cell. We can let the adjacent lattice points be included in different sets A and B since the lattice is bipartite. Assuming that the lattice point at the junction is included in the set A without loss of generality, we find that the set A will have 3m + 3 sites and the set B will have 3m + 2 sites. Now, we index the lattice points in the set A with the integers from 1 to 3m + 3 and the lattice points in the set B with the integers from 3m + 4 to 6m + 5. Then, since the matrix representation of the tight-binding Hamiltonian is zero for all (i, j) except when i-th and j-th lattice points are adjacent. This makes the Hamiltonian block-off-diagonal, where the blocks are (3m + 3) × (3m + 3) and (3m + 2) × (3m + 2) sized zero matrices. Hence, the rank of the matrix is lower than the dimension of the matrix by 1, which signals the emergence of the flat band. Fig.5 is an example. This honeycomb network has three sites between tri-junctions, and the sites are indexed following the above method. The Hamiltonian, which is 11 × 11, can be written as the two rectangular off-block-diagonal matrices.
However, we remark that such factorization is absent (at least we couldn't find it after some trials) for the generic cases with further neighbor hoppings. Hence, this approach cannot be used to prove the existence of many stable flat bands. generalizes. From the previous work, we also have observed that C 3 × R x cannot lift the flat bands. 1

Appendix F: Details of BCS Calculations
Here we present the details of the BCS mean-field calculation. This involves the projection of the bare lattice-scale interactions to the BCS channels and the calculation of the mean-field gap equation and energy.

Coupling Constants
The microscopic pairing interaction that we introduced phenomenologically in the main text can be generally expressed as where U and V are real-valued constants,n R;a is the number operator at position R and a-th site and a, b and "end" indicate the adjacent sites in real space. We are going to project this interaction terms to the BCS pairing term. First, we consider the spin degree of freedom then the first term can be expanded.
We only select n ↑ n ↓ pairs in this expansion since we are interested in the conventional BCS channel. By using the Fourier transformation, the first pairing term becomes where the form factor is defined as the following.
F U (p, q; l, k) = a u k;a u * p;a u * q;a u l;a (F4) Restricting the momentum summation into the pairing channels only, we obtain the BCS channel.
Similarly, we calculated the second and the third terms also. Thus, Next, we further perform the expansion of g(p, k) in terms of the angular momentum sectors, i.e., the s-wave and the d-wave channels.

Mean-Field Solution and Energy
We define the pairing order parameter.
The form factor F l (k) and the g l is defined as for the s-wave case, and   TABLE II. Table of {f l U , f l V , l = s, d} for the system with foive sites between the junctions. Note that the band index means the "n-th" lowest band.
negative to be attractive.
Then, the BCS Hamiltonian is given as the following, ignoring the constant term for a moment.
We find the unitary matrix Q which diagonalize the Hamiltonian.
Then we define states, which are number-conserving, using the eigenvectors.
The normalized BCS ground state is given as the following.
By definition of the pairing order parameter, we calculate it with the BCS ground state.
Assuming that the Fermi energy is exactly at the flat band, we have ξ k = 0.
Finally, we can get the BCS mean-field energy.
As the result of the numerical calculation, the BCS mean-field energies for each possible superconductor type are where A is area of the Brillouin zone. With the BCS mean-field energy, we draw the phase diagrams for the tightbinding honeycomb network models. Only the nearest-neighbor hoppings were considered. For the completeness, we also have included the ferromagnetism in Fig.11. The phase diagrams are drawn at each flat band of the model. We may turn on the next-nearest-neighbor hoppings, but the phase diagrams were slightly changed but they do not induce much difference. Note that, when U is small, then the window for the d-wave SC is also small. However, when U is sizable and positive, then the window for the d-wave is also large.

BdG Fermion Spectrum
We draw the BdG fermion spectrum.
which has the flat band with the quadratic band touching at k = 0. With the BCS pairing interaction, the fermion spectrum becomes the following.
hence there are four bands at the low-energy limit. The BdG fermion spectra of the s-wave and the d-wave are plotted in Fig.12 near the Γ-point. Starting from the QBT with a flat band, the band gap are opened by the pairing interaction symmetrically. Note that they are gapped for the s-wave and gapless with doubled quadratic band touchings for the d-wave at the Γ-point. This state potentially has an interesting quantum critical behavior, which we leave for the future study.

Appendix G: Corner States at Tri-Junctions
In this supplemental information, we provide a few concrete realizations of the zero-dimensional states localized at the tri-junctions (or nodes) of the honeycomb network. We concentrate on the half-filled per wire case, in which the leading insulating instability is the period-2 charge-density waves.
The zero-dimensional states is the soliton of the charge-denstiy waves, and they carry the fractional quantum numbers, e.g., the electric charge. Here we assume SU(2) symmetry and so consider effectively spinless fermions. Now, imagine that the tri-junction induces the frustrations for the charge-density wave order parameters between the neighboring links of the honeycomb plaquette. Due to the commensurability of the filling, the only allowed frustration is the π-phase shift between the neighboring links and is trapped at the junction. Once such mode is trapped at the junction, it cannot move around when the crystalline symmetries are protected. [Here we do not allow the Hilbert space to be changed when we consider the symmetric deformations of the model.] Then, we find that the 2d domains are insulating and its first-order boundaries are also insulating, but only its second-order boundary, i.e., the corner, has "in-gap" modes, which can be protected by the crystal symmetries. This is very parallel with the corner modes in the higher-order topology, (or more precisely, obstructed atomic insulator with the non-trivial nested Wilson loop topology). [4][5][6][7] Note that, in the reference 8 , the junctions of gapped wires are considered in an entirely different context, and the emergent corner states are discovered. Our finding is consistent with theirs.
It is well-known that, in 1d, a soliton of the period-2 charge density waves is the same as the boundary mode of the Su-Schieffer-Heeger model. In fact, our construction here is essentially the charge-density wave verions of the higher-order topology. In this paper, we will not attempt to present the full theory and classification. Instead, we will present only the minimal contents, and the precise connection and classification of the "higher-order topological" domain wall states will be reported elsewhere.

Tight-Binding Models
Two examples of tight binding models exhibiting the localized corner state is shown in Fig.13. Because of the period-2 modulation, we modulate t and g in each link. Infinitesimal on-site energies, ±m ,were given at each junctions to break a symmetry weakly (this is commonly done in the investigation of the corner charge in higher-order topology and polarization chain). [4][5][6][7] They are explicitly written in Eq.G1 and Eq.G2. Then, we calculated the localized charges at each junction at half-filling. The charges at A(Red) and B(Blue) site showed opposite signs but their amplitude was nearly 0.5 when g is smaller enough than t (i.e., small correlation length limit). The better localization of electric charge is expected when more sites are assumed on the wires as the Su-Schrieffer-Higger model does.  The flat bands are split into two flat bands at (E 0 + h) for spin-down and (E 0 − h) for spin-up so the density of states is given as the following.
The number of occupied electrons for each flat band is where σ =↑, ↓ and f ( , µ) is the Fermi-Dirac distribution. Magnetization is proportional to the difference of number of the spin-up and spin-down electrons. The spin susceptibility near the flat band shows a similar tendency with the specific heat as indicated in Fig.14. The spin susceptibility is zero at the zero-temperature. As the fermi energy approach to the flat band energy, the spin susceptibility diverges at the zero-temperature.