Nematicity-enhanced superconductivity in systems with a non-Fermi liquid behavior

We explore the interplay between nematicity~(spontaneous breaking of the sixfold rotational symmetry), superconductivity, and non-Fermi liquid behavior in partially flat-band models on the triangular lattice. A key result is that the nematicity (Pomeranchuk instability), which is driven by many-body effect and stronger in flat-band systems, enhances superconducting transition temperature in a systematic manner on the $T_{\rm c}$ dome. There, a $s_{x^2+y^2} - d_{x^2-y^2} - d_{xy}$-wave symmetry, in place of the conventional $d_{x^2-y^2}$-wave, governs the nematicity-enhanced pairing with a sharp rise in the $T_{\rm c}$ dome on the filling axis. When the sixfold symmetry is spontaneously broken, the pairing becomes more compact in real space than in the case when the symmetry is enforced. These are accompanied by a non-Fermi character of electrons in the partially flat bands with many-body interactions.


INTRODUCTION
Strongly correlated systems have become an epitome in the condensed-matter physics, as exemplified by the high-temperature superconductivity in the cuprate [1], and iron-based [2] families. These compounds exhibit rich phase diagrams as hallmarked by the emergence of unconventional superconductivity, and a plethora of symmetry-broken phases such as spin and charge nematicity and stripe orders.
Quest for finding novel high-temperature superconductors spurs interests in exploring many-body systems with short-range repulsions but with (nearly) flat subregions in the band dispersion arising from hopping beyond nearest neighbors or from lattice structures [3,4]. These systems with dispersionless band portions permit numerous scattering channels for the electrons and can give rise to various exotic quantum phases such as spin and charge density waves [5,6], Mott insulating [7], and bad-metallic phases [8], as well as the formation of spatially extended Cooper pairs [9,10].
Interaction and the flatness of the band structure can be intimately related to geometric and quantum frustration in producing strong correlation effects. The spin liquid behavior in hexagonal lattices, such as organic compounds [11,12]) and inorganic Herbertsmithites [13,14], are typical examples. In these exotic liquids, the classical picture is no longer valid, and their quantum phase transitions cannot be described within Landau's phase transition theory.
Aside from these many-body phenomena, the electron nematicity, i.e., spontaneous breaking of spatial rotational symmetry triggered by many-body interactions, is another manifestation of the correlation effects [15,16].
It is an intriguing direction to pursue the physical origins of these symmetry-broken phases [17][18][19], and to grasp the interplay between nematicity and other phases such as superconductivity [10,[20][21][22][23] and non-Fermi liquid [24]. Various studies report different roles of nematic fluctuations on the superconductivity including the competition between these two phases, e.g., in doped BaFe 2 As 2 [25], the assistance of nematicity to enhance the superconductivity transition temperature, e.g., in twisted bilayer graphene [26], and the negligible effect of nematicity on the superconducting phase, e.g., in FeSe [27]. One crucial aspect is figuring out which of these possibilities occur in systems that have a flat or partially flat band in their dispersions [5,6,28,29].
In this paper, we bring these features together to explore the interplay between the nematicity and superconductivity in partially flat band (PFB) models on the triangular lattice, effective model for Moire-produced basis states [30]. Here the lattice structure frustrates magnetic orders, thereby giving opportunities for nematic instabilities to arise. As a key finding, we shall demonstrate that nematicity can significantly enhance transition temperatures (T c ) in the superconducting phase, with a s x 2 +y 2 −d x 2 −y 2 −d xy -wave pairing symmetry. This occurs for an intermediate Hubbard repulsion and in a non-Fermi liquid regime.

MODEL
The Hubbard Hamiltonian on the isotropic triangular lattice reads where c † kσ (c kσ ) creates (annihilates) an electron with momentum k = (k x , k y ) and spin σ at site i, n iσ ≡ c † iσ c iσ . The repulsive Hubbard interaction is denoted as U (> 0), and µ is the chemical potential. The non-interacting band dispersion for the triangular lattice is given as where t is the nearest-neighbor hopping (taken as a unit of energy) and t is the second-neighbor hopping. Here, we consider (t, t ) = (1.0, 0.15), which possesses a nearly flat region along K − K − K , see dashed line in Fig. 1.
For the interaction, we set an intermediate U = 4.5t, with the inverse temperature set to be β ≡ 1/(k B T ) = 30/t except in Fig. 3(c).

NUMERICAL METHOD
To study paramagnetic phases with no spin imbalance, we employ the dynamical mean-field theory (DMFT) combined with the fluctuation exchange approximation (FLEX), known as the FLEX+DMFT [31]. This method comprises DMFT and FLEX double loops, solved self-consistently at each FLEX+DMFT iteration. In this work, we solve the DMFT impurity problem by the modified iterative perturbation theory [32,33]. The momentum-dependent FLEX self-energy is constructed from the bubble and ladder diagrams. Af-  ter removing the doubly-counted diagrams in the local FLEX self-energy, the FLEX+DMFT self-energy is updated. The momentum-dependent self-energy in the FLEX+DMFT incorporates vertex corrections generated from the DMFT iterations into the local part of the FLEX self-energy. Even though our FLEX+DMFT method does not deal with spatial vertex corrections, larger coordination number and frustrated magnetic fluctuations in the triangular lattice give rise to more local self-energies and less dominant spatial vertex corrections than in the square lattice [34]. As a result, the FLEX+DMFT is considered to be a reliable approach that incorporates local and nonlocal correlations.
When we start from the non-interacting tight-binding Hamiltonian, Eq. (1), that has the sixfold rotational (C 6 ) lattice symmetry, the solution of the many-body problem may exhibit a lower the symmetry. To study the phases with/without C 6 symmetry, we solve the FLEX+DMFT loops with/without imposing the C 6 constraints. To explore the Pomeranchuk instability with the broken C 6 symmetry, we take an initial self-energy as Σ in = 0.05[cos(k x ) − cos( √ 3k y /2) cos(k x /2)] which acts as a seed for distorting the Fermi surface for the FLEX+DMFT iterations. The FLEX+DMFT calculations are here performed on a 64 × 64 momentum grid and an energy mesh with 2048 points.

NEMATICITY AND NON-FERMI LIQUID BEHAVIOR
We start with presenting the momentum distribution function plotted in panels (a-c) in Fig. 2  the maxima of the momentum-dependent distribution function are below unity. The system exhibits a fillingdependent degrading of C 6 down to a twofold C 2 symmetry in n k . Namely, we have here an emergence of nematicity, or a Pomeranchuk instability. The breaking of C 6 is seen to occur even right at half-filling, while the electron-doped case shows a preserved C 6 . To quantify the broken C 6 symmetry, we introduce point-group resolved Pomeranchuk order parameters defined as ξ d x 2 −y 2 = k d x 2 −y 2 (k)n k and ξ dxy = k d xy (k)n k , with k = 1 [35]. The form factors, d x 2 −y 2 = cos(k x ) − cos( √ 3k y /2) cos(k x /2) and d xy = √ 3 sin( √ 3k y /2) sin(k x /2), describe the distortion of the Fermi surface in the point group C 6 , and ξ is a real number with values between zero (when C 6 is preserved) and unity. Figure 3(a) displays ξ against the band filling. We can see that, as the band filling is reduced, ξ starts to grow, and at a critical band filling n c1 = 1.02 [vertical blue line in Fig. 3(a)], ξ dxy undergoes a first-order phase transition [36,37]. At this filling, the onset of nematicity is accompanied by a Lifshitz transition, where the Fermi surface delineated by the ridges in |G k | 2 (not shown) is not only distorted but undergoes a topological change from closed to open structures.
We further notice that the filling dependence of the nematicity differs between ξ d x 2 −y 2 and ξ dxy in the PFB model; compare purple and magenta lines in Fig. 3(a). For 0.986 < n < n c1 , ξ d x 2 −y 2 is dominant, while ξ dxy takes over below n = 0.986, which we call the second characteristic band filling, n c2 [vertical dashed sky-blue line in Fig. 3(a)]. While ξ d x 2 −y 2 displays a first-order transition at n c1 , ξ dxy exhibits a crossover at n c2 . This suggests that thermodynamic parameters such as temperature at which ξ dxy and ξ d x 2 −y 2 experience the first-order transitions are different from each other.
To trace back the origin of the nematic phases, let us next present the momentum-dependent spin susceptibility χ s (k) for the PFB model in Fig.2(d-f). In the electron-doped regime where the Pomeranchuk instability is absent, χ s respects the sixfold rotational symmetry of the lattice, with peaks at k = ( √ 3π/2, 0) and its equivalent positions under C 6 . As band filling is decreased below the half-filling, the spin susceptibility develops spikes around k = ( √ 3π/3, 2π/3) and the equivalent places under a C 2 subgroup of the original C 6 rotational symmetry. The appearance of spikes at mid-momenta in the spin susceptibilities indicates the presence of long-range spin fluctuations in our systems [38].
In general, an electronic nematicity without breaking the translational symmetry can be driven by structural transitions, charge [39] or spin [40] fluctuations. Our Hamiltonian does not deal with the distortion of the lattice or phonons and thus precludes structural transitions. We have checked that the charge susceptibility is at least an order of magnitude smaller than the spin susceptibility. Thus the spin-mediated correlations should be responsible for the emergence of the Pomeranchuk instability [41]. Now let us turn to a non-Fermi liquid character of the present electronic systems, since the flat portions of the band may well exert peculiar effects. We can quantify this in terms of the impurity self-energy in DMFT by fitting the imaginary part of the self-energy on Matsubara axis to |ImΣ DMFT (iω n )| ∝ ω α n , and present the result for the exponent α in Fig. 3(b). In general, α = 1 at small ω n (c.f., α = 2 on real frequency axis as |ImΣ DMFT (ω)| ≈ max(ω 2 , T 2 ) at small T ) characterizes the Fermi liquid, while α < 0.5 will signify a non-Fermi liquid (bad metal) behavior [42][43][44][45]. Above the first order Pomeranchuk transition for n > n c1 , α's computed for systems with (dashed lines) and without (solid lines) the enforced C 6 constraint trivially coincide with each other. We can see that both systems display strong non-Fermi liquid behavior with α well below 1. If we turn to n < n c1 for which we have revealed the nematicity, Fig. 3(b) shows notable differences in α between the cases where C 6 is enforced or not. After a sharp drop at n c1 as the band filling is reduced, α gradually increases (decreases) in the presence (absence) of the imposed sixfold constraint. Eventually α starts to decrease with decreas-  ing n at n ≈ 0.85 in the PFB system. The persistent α < 0.5 for n < n c1 implies that the nematic phase resides in the non-Fermi liquid regime.

SUPERCONDUCTIVITY
Now let us come to our key interest in pairing instabilities, for which we solve the linearized Eliashberg , to find the largest eigenvalue λ for the spin-singlet, evenfrequency superconducting gap function ∆. Here, k ≡ (k, iω n ) with ω n the fermionic Matsubara frequency, and the effective interaction for singlets given as The pairing is identified when λ exceeds unity [46]. Figure 3(c) depicts λ in the presence (dashed green lines) or absence (solid blue) of imposed C 6 symmetry in PFB.
When the sixfold rotational symmetry is enforced, we get λ < 0.8 in PFB model, indicating that the singlet superconductivity does not arise for the temperature (k B T = t/30) considered here. We can still notice that λ displays a double-peak structure with a minimum at n min = 0.95. The dip is shown to occur at the band filling at which the d x 2 −y 2 gap function with two-nodal lines for n > n min changes into a more complicated multinodal-line gap functions for n < n min , see [47,48] for details. This behavior of the gap function reflects a crossover from the antiferromagnetic spin structure with a single nesting vector for n > n min , to a more complex spin structure for n < n min where single peaks in the spin susceptibility evolve into extended structures [see Fig.2(d-f)]. Thus the system for n < n min goes beyond the conventional nesting physics. Similar structure in λ and associated gap function have also been reported for PFB systems on the square lattice [9], again in the absence of nematicity.
In a dramatic contrast, if we allow the C 6 symmetry to be broken spontaneously, λ soars from those with C 6 restriction, as seen for n c < n < 1.15. This occurs concomitantly with the Pomeranchuk order parameters (ξ's), which grow precisely in this filling region. Just below n c , λ in the systems with broken C 6 (solid blue lines in Fig. 3 (c)) exhibits a rapid growth and exceeds unity. This is inherited in the superconducting transition temperatures (T c ), presented in Fig. 4(a).
T c , with the broken C 6 , exhibits a single-dome structure as a function of band filling. We can observe that the presence of a flat portion in PFB or a van Hove singularity for t = 0 have similar effects on the largest values of T c when ξ dxy ≥ ξ d x 2 −y 2 . One should note that, while a van Hove singularity at E F only occurs at a single point on the filling axis, a flat portion of the band can accommodate a range of band filling. This difference is reflected in the width of the T c dome at a given temperature; see Fig. 4(a) and SM [48]. The maximum of T c in the PFB is seen to take place close to n c2 at which ξ dxy exceeds ξ d x 2 −y 2 . Note that the superconducting transition temperature becomes almost doubled as we pass through n c1 , see Fig. 4(a), which should come from the interplay between nematicity, spin fluctuations, and superconductivity.
Let us now delve into the gap function in momentum space in Fig. 4(b-d). In the electron-doped regime, the PFB model exhibits a conventional d x 2 −y 2 paring symmetry [49]. This behavior of the gap function persists for n > n c . On the other hand, below n c1 where the C 6 symmetry is broken down to its C 2 subgroup, the dominant channel of instability is a mixture of s x 2 +y 2 , d x 2 −y 2 and d xy -wave symmetries, see Fig.S16 in SM [48].
To better understand the role of nematicity in superconducting phases, we can look at eff is the effective interaction with the imposed C 6 constraint. As shown in Fig. S19 in SM [48], V eff is much intensified when C 6 is lifted. Since χ c is much smaller than χ s , the effective pairing interaction reflects the momentum-dependence of the spin-susceptibility under the Pomeranchuk distortions. This effective interaction assists electrons to nonlocally form Cooper pairs [15,37,50]. The deformation in ∆V eff allows firstorder perturbation corrections in the distortion, which should be responsible for the drastic changes in λ below n c . This contrasts with the previous study on the interplay between nematicity and superconductivity, where the enhancement of λ originates from the second-order perturbation corrections and thus results in much smaller changes [41].

DISCUSSION AND SUMMARY
We have studied whether and how an emergent nematicity affects superconductivity in partially flat bands on the regular triangular lattice. We have shown with the FLEX+DMFT that nematicity dramatically affects pairing symmetry, and the T C is significantly enhanced by the lowered point-group symmetry in the electronic structure. This is shown to occur in a non-Fermi liquid regime, which is characterized by blurred Fermi surfaces, momentum-dependent fractional occupations of the band, and a fractional power-law in the self-energy.
In the presence of nematic order, the superconducting symmetry changes from an (extended) d x 2 −y 2 -wave to a s x 2 +y 2 −d x 2 −y 2 −d xy -wave, where unlike the conventional nesting-driven case, the pairing interaction is governed by an intricate spin susceptibility structure.
Future works should include the elaboration of the way in which the non-Fermi liquid property affects the superconductivity, and exploration of the interplay between Pomeranchuk instability and superconductivity in multiband/orbital systems with flat regions.
, is computed by analytically continuing the DMFT impurity Green's function on Matsubara frequency axis with the Padé approximation [51]. Similarly, the momentum-dependent FLEX spectral function, A(k, ω) = −1/π Im[G(k, ω)], is evaluated from the momentum-dependent Green's function with an analytical continuation.
Normalized double-occupancy, n ↑ n ↓ /( n ↑ n ↓ )), is the ratio between the number of doubly occupied sites n ↑ n ↓ to the uncorrelated value n ↑ n ↓ .
Momentum-dependent distribution function is given by Pomeranchuk order parameters, ξ d x 2 −y 2 and ξ dxy , which quantify the amount of symmetry breaking [36], are given on the triangular lattice as [35] ξ describe the component-resolved distortion of d-wave instabilities for the Fermi surface distortion in the point group C 6 . When the C 6 symmetry is preserved, both ξ d x 2 −y 2 and ξ dxy are equal to zero. Static spin and charge susceptibilities. The static spin susceptibility reads where τ denotes imaginary time and β is the inverse temperature, while the static charge susceptibility is given by The local spin susceptibility is then given by while the DMFT impurity spin susceptibility is calculated as where the polarization function is χ 0 (0) = − ωn G imp (iω n )G imp (−iω n ), with G imp being the DMFT impurity Green's function and ωn = 1. Exponent of the impurity self-energy. To quantify non-Fermi liquid behavior we can fit the imaginary parts of the DMFT self-energies at small Matsubara frequencies (ω n ) to For Fermi liquids, ImΣ(iω n ) ∝ iω n with the exponent α = 1 on the Matsubara axis (i.e., ImΣ(ω) ∝ ω 2 on the real frequency axis). If α < 0.5, the bad metallic behavior is indicated as a signature of non-Fermi liquid. Superconducting gap functions. For pairing instabilities we solve the linearized Eliashberg equation to find leading eigenvalues λ and superconducting gap functions ∆ as where k = (k, iω n ) with the fermionic Matsubara frequency ω n , G denotes the single-particle Green's function, and the effective interaction for the spin-singlet pairing, In this paper, we focus ouselves on even-frequency gap functions with even parity in the Matsubara frequency, which satisfy ∆(k, iω n ) = [∆(k, iω n ) + ∆(k, −iω n )] /2. Superconductivity is identified when the largest eigenvalue λ exceeds unity.
Appendix B: Ideal partially flat-band model In this section, we explore how the way in which the flat portion in the band dispersion is prepared affects the physics. So here we introduce the Hubbard model on a model [which we call the ideal partially flat-band (iPFB)] with a truncated dispersion given as Namely, the dispersion is made perfectly flat for energies below κ for a tight-binding model having a nearest-neighbor t and second-neighbor t with the dispersion, see Fig. S1(a). Thus the larger the κ, the more extended the flat region. We have here inverted the signs of t and t to put the van Hove singularity near the band bottom for the triangular lattice (occupied at small band fillings). The band width is that of Eq.(C2) minus κ. Here we take the band width to be W = 7.533t.

Results for iPFB
Let us present the numerical result for iPFB, here for U = 4.5t and, except in Figs. S4(d), S7(b), S12(a), the inverse temperature is β = 1/(k B T ) = 30/t. We take t as the unit of energy as in the main text.

a. Nematicity and non-Fermi liquid behavior
In search for footprints of non-Fermi liquid behavior, we start with the absolute value of Green's functions |G k | 2 (top rows) and the momentum-dependent occupation number (middle) in Figs. S2, S3 for the iPFB models with κ = 1.0 and 1.5. The iPFB model displays an absence of sharp peaks in |G k | 2 as for t = 0 (not shown), i.e., absence of welldefined Fermi surfaces, which hints that non-Fermi physics governs these systems. A further signature of non-Fermi liquid characters is evident in n k in which we observe max[ n k ] < 0.95 for all presented band fillings in the iPFB model.
In the hole-doped regime of the iPFB model with κ = 1.0, we observe a broken sixfold symmetry, see panels (a,b) and (e,f) in Fig. S2. For κ = 1.5 in Fig. S3, the reduction from C 6 to C 2 occurs even at half-filling, while the symmetry is resumed for smaller band fillings, namely n < 0.8. These results suggest that the onset of broken C 6 symmetry is shifted to higher band fillings as we increase the size of the flat portion. Thus we can conclude that the spontaneous nematicity occurs both in the PFB and iPFB, even though the detail of band dispersion is significantly different between the two models.
If we look at the normalized double occupancy against band filling in Fig. S4(a) to see how this quantity is correlated with the spontaneous breaking of C 6 . Above a particular filling, the normalized double occupancy in systems with (dashed curves) and without (solid curves) enforced C 6 symmetry are identical. Below this critical band filling, which we call n c , the Pomeranchuk instability occurs as a first-order phase transition. Figure S4  as n c = 0.91 for κ = 0, n c = 0.97 for κ = 1.0, and n c = 1.08 for κ = 1.5. For κ = 1.5, if we further decrease the filling to n < 0.8 the C 6 symmetry is restored as we have already noted in Fig.S3. We refer to n cl = 0.81 as the lower critical filling.
We can also note that we have a Lifshitz transition at the critical band filling n c , i.e., the closed ridges in |G k | 2 above n c are topologically changed to an open structure, c.f. panels(b,c) in Fig. S2 and (c,d) in Fig. S3. The second Lifshitz transition occurs at the lower n cl , below which the open ridges in |G k | 2 resume a closed structure, see Fig. S3(a).
If we actually examine in Fig. S4(b) the Pomeranchuk order parameters (ξs) against band filling for κ increased from 0 to 1.5 at β = 30. For κ = 0, ξ dxy (light green squares) and ξ d x 2 −y 2 (dark green circles) experience a first-order phase transition at a n c . Below this critical filling, both orders grow, with ξ dxy larger of the two. For κ = 1.0, we observe a sharp rise in ξ d x 2 −y 2 at a shifted n c accompanied by a sharp rise in ξ dxy at n c2 = 0.964 (vertical dotted orange line). When κ is increased to 1.5, the critical band filling n c for ξ d x 2 −y 2 further increases, while ξ dxy has a jump at n c2 = 1.027 (vertical dotted sky-blue line) accompanied by a drop in ξ d x 2 −y 2 . For n < n c2 , the dominant Pomeranchuk instability has d xy character. This continues until the lower critical filling n cl = 0.81 below which the iPFB system retrieves the C 6 symmetry and all instabilities are gone.
We now look at their temperature dependence of ξs at respective n c in Fig. S4(d). For κ = 0, the rapid growth of both ξ d x 2 −y 2 and ξ dxy simultaneously occur around T ≈ 0.05. It is evident that these nematic orders at n c survive at higher temperatures with ξ dxy remaining dominant. For κ = 1.0, we see a first-order phase transition in ξ d x 2 −y 2 at T ≈ 0.038 followed by a sharp increase in ξ dxy at T ≈ 0.055. At T = 0.066, the two Pomeranchuk instabilities become comparable. For κ = 1.5, we witness a first-order phase transition for ξ d x 2 −y 2 at T ≈ 0.038, dark blue line in Fig. S4(d). This order parameter remains the only instability at higher temperatures as at n c , with ξ dxy almost vanishing as the temperature is increased.
We turn to the spin susceptibility in Figs. S2, S3(bottom rows) for κ = 1.0 and 1.5, respectively. In the electrondoped regime where the Pomeranchuk distortion is absent, the spin susceptibility exhibits spikes corresponding to antiferromagnetic spin structure at k = ( √ 3π/2, 0) with it's equivalent sixfold k-points. Below n c , the spin susceptibility exhibits two streaks, reflecting the C 2 symmetry. These are centered around k = ( √ 3π/3, 2π/3), and its twofold equivalent places. At n = 0.9 for κ = 1.0, the streaks in the spin susceptibility are rotated by π/6 and located at k = (− √ 3π/3, 2π/3), and its twofold equivalent places. For κ = 1.5, we can confirm that the sixfold symmetry is resumed in the spin susceptibility as well below n = n cl . In agreement with the conclusion of the main text, we can thus infer that the nematicity is triggered by spin fluctuations.
We can further look at the impurity and local spin susceptibilities in Fig. S5(a,b), respectively. Right at n = n c , the χ imp s and χ loc s in systems without the enforced C 6 symmetry exhibit drastic deviations from the C 6 -imposed results. At n c2 with κ = 1.0 and 1.5 where the dominant instability changes from ξ d x 2 −y 2 to ξ dxy , the impurity and local spin susceptibilities exhibit small kinks.
In addition, at fillings below n c (solid vertical lines in Fig. S5), the difference between χ loc s in the presence (dashed lines) and absence (solid) of C 6 becomes more pronounced as the flat-band region is widened. This is as expected, since the many-body effects, reflected in χ s (k), are stronger in partially flat band systems. This property has also been noticed for partially flat-band systems studied with the determinant Monte-Carlo [8] and FLEX+DMFT [9]  calculations on the square lattice.
To further track the fingerprints of the nematic orders, we present in Fig. S6 the momentum-dependent spectral function, A(k, ω), in iPFB with κ = 1.0 (b) or κ = 1.5 (c) at n = 0.9. In the presence of the interaction, the flat region at the bottom of the bands seen in A(k) changes little, aside from a downward energy shift, except along the high symmetry lines where C 6 is broken as identified in Figs. S2, S3. More specifically, for κ = 1.0, the nematicity vector is along K → K, which changes to K → K for κ = 1.5, see labels in Fig. S1(b). For the momenta in these regions, the spectral function splits with a gap U/2.
To corroborate the non-Fermi liquid character of electrons in our models, we plot the exponent α of the impurity self-energies in systems with (dashed lines) and without (solid) the sixfold symmetry in Fig. S4(c) above. For the whole band-filling region studied, α exhibits values less than 0.5, indicating non-Fermi liquids. In the electron-doped regime for n > n c with unbroken C 6 , α grows as we approach n c as the band filling is decreased. Just below n c , α exhibits a sharp drop, followed by a gradual enhancement accompanying the Pomeranchuk instabilities. Close to n = n c2 , α in the presence of nematicity (solid lines) undergoes another drastic reduction.

b. Superconductivity
We now turn to superconductivity in the iPFB. We present the largest eigenvalue λ of the Eliashberg equation for the singlet pairing in Fig. S7(a). If we first look at the result when the C 6 symmetry is imposed (dashed lines), λ remains below 0.85, i.e., no superconducting instabilities. Similar to our finding in the main text for the PFB model, we can characterize the double-dome structure in λ as a signature of a change in the number of nodal lines from 2 in the right dome to ≥ 4 for the left dome. The band filling at which these two domes are connected is indicated by arrows with the same color as the corresponding dashed lines in Fig. S7(a). When the nematicity is allowed,  in a sharp contrast, λ becomes significantly higher, and specifically rapidly exceeds unity just below the respective n c . After this peak, λ decreases as the filling is decreased, until a second rise occurs at n = n c2 , followed by a second decrease. Corresponding superconducting transition temperatures is displayed in Fig. S7(b). The double-dome structure is inherited in T c , which endorses the interplay between various symmetries of Pomeranchuk distortions and the superconducting instabilities.
To identify the pairing symmetry, we present the singlet gap functions in momentum and real spaces for the regular band in Fig. S8, for iPFB model with κ = 1.0 in Fig. S9, and for κ = 1.5 in Fig. S10. Above n c , the gap functions have d x 2 −y 2 symmetry [see panels (d) in these figures], where the pairing in real-space extends to one lattice spacing [panels (h)]. As in Fig. 4 in the main text for PFB, the gap function changes from the d x 2 −y 2 -wave to a γs x 2 +y 2 − d x 2 −y 2 − d xy -wave for n cl < n < n c when Pomeranchuk distortions emerge, see panels (a,b) in Fig. S8 for the regular band, panel (a) in Fig. S9 and panels (a,b) in Fig. S10 for iPFB. The value of γ increases from zero to almost 0.3 when the filling is reduced from n c up to n c2 . γ acquires values larger than 0.3 below n c2 . In iPFB with κ = 1.0 at n = 0.9 the superconducting gap function exhibits a π/6 rotated s x 2 +y 2 − d x 2 −y 2 − d xy -wave pairing, i.e.,s x 2 +y 2 − d x 2 −y 2 + d xy -wave pairing, which reflects the momentum-dependency of its spin-susceptibility. Below n cl , the pair becomes spatially extended (panel (e) in Fig. S10), which accompanies the appearance of multiple nodal lines in the k-space (panel (a)). While this is expected for the systems with partially flat bands [9] because of a bunch of pair-scattering channels around the flat regions, suppression of these types of pairings in the presence of Pomeranchuk instabilities is interesting. This can be understood by recalling that a nematicity, which brings about some superstructures, degrades the flatness of the noninteracting band dispersion. Concomitantly, the spin susceptibility acquires spikes structures.
Aside from the discussed singlet superconductivity, we have also examined the triplet gap functions and associated λ. Our results show that the λ never exceeds 0.4 in the region studied. This is why we have not presented these data in the present work.  , exhibit a singular value as small energies. A large DOS is also observed at ω ≈ −1.68, slightly above the bottom of the band, for the regular band system. Note that while the bandwidths of these four systems are identical, the bottom of their band structure is not located at the same energy.  To further explore the temperature-dependence of the Pomeranchuk instabilities, we plot ξ d x 2 −y 2 and ξ dxy at n c against temperature in Fig. S12(a). While for T < 0.0425 both of the nematic orders are negligible, ξ d x 2 −y 2 undergoes a first-order transition at T = 0.0425, which indicates the broken C 6 symmetry at higher temperatures.

Impurity spin-susceptibility
Impurity spin-susceptibility for PFB model with t = 0.15 in Fig. S12(b) shows that the increasing spin fluctuations for n > n c sharply drops when sixfold rotational symmetry is broken for n < n c . Local spin-susceptibility in Fig. S12(c), on the other hand, displays a maximum at n = n c2 in PFB. We can note that we hardly detect any features at n c2 in the impurity spin-susceptibility which might be related to a crossover from ξ d x 2 −y 2 to ξ dxy , see Fig. S12(a). When the structure of ξ undergoes an abrupt change from ξ d x 2 −y 2 to ξ dxy , due to the first-order transition in ξ dxy , the impurity spin-susceptibility exhibits a kink at n c2 , see e.g. Fig. S4(b) and related discussion in Sec. B. 〈n〉=1.0 Figure S14. Difference, ∆V eff = V eff − V C 6 eff , in the effective pairing interaction between the absence and presence of the C6 symmetry for the PFB system with t = 0.15 for U = 4.5 and β = 30 for n = 0.9 (a) and 1.0(b). Black hexagons represent the Brillouin zone. Note that the scale of the color bar differs by orders of magnitude between the panels.

Superconductivity
Fig. S14 presents ∆V eff = V eff − V C6 eff , where V C6 eff is the effective interaction with the imposed C 6 constraint. We can see that ∆V eff displays not only the broken C 6 symmetry but also amplitudes significantly intensified than in the C 6 case. See also the discussion on ∆V eff in the main text.
Gap functions in the real-space with singlet, even-frequency pairing for the PFB model with t = 0.15 is plotted in Fig. S15. In the main text, we presented the γs − d xy -wave symmetry of the gap function below n c in the momentum space; see Fig. 4(b-d). Evidently, the electron pairing is short-range in this filling region. This observation is in contrast with the formation of extended Cooper pairs in systems where the C 6 constraint is imposed. ∆ dxy (k) = √ 3 sin( √ 3k y /2) sin(k x /2).
In the presence of C 6 symmetry, the above two d-waves, depicted in Fig. S16(a-b), are associated with two degenerate irreducible representations in the C 6 point group, so that the topological superconductivity with gap function ∆ d x 2 −y 2 (k) + ı∆ dxy is also possible. When C 6 symmetry is reduced down to C 2 , on the other hand, ∆ s x 2 +y 2 and ∆ d x 2 −y 2 are no longer irreducible representations, and instead their linear combination s x 2 +y 2 − d x 2 −y 2 , ∆ s−d , with a form factor, ∆ s−d (k) = 3 cos( √ 3k y /2) cos(k x /2), finds room to emerge, see Fig. S16(c). In the C 2 point group, s x 2 +y 2 −d x 2 −y 2 reads s y2 , and in this point group s y2 and d xy are not degenerate, so that this group does not permit topological superconductivity. Still, linear combinations of ∆ dxy and ∆ s−d may be realized, as we have indeed shown in Figs. 4(bottom rows), S8, S9, S10. Fig. S16(d-f) shows α(s x 2 +y 2 − d x 2 −y 2 ) + βd xy with (α, β) = (0.