Spin-Triplet Pairing Induced by Near-Neighbor Attraction in the Cuprate Chain

In quantum materials, the electronic interaction and the electron-phonon coupling are, in general, two essential ingredients, the combined impact of which may drive exotic phases. Recently, an anomalously strong electron-electron attraction, mediated by phonons, has been unveiled in one-dimensional copper-oxide chain Ba$_{2-x}$Sr$_x$CuO$_{3+\delta}$. Yet, it is unclear how this strong near-neighbor attraction $V$ influences the superconductivity pairing in the compound. Here we perform accurate many-body calculations to study the extended Hubbard model with on-site Coulomb repulsion $U>0$ and attraction $V<0$ that well describes the cuprate chain and likely other similar transition-metal materials with both strong correlations and lattice effects. We find a rich quantum phase diagram containing an intriguing Tomonaga-Luttinger liquid phase -- besides the spin density wave and various phase separation phases -- that can host dominant spin-triplet pairing correlations and divergent superconductive susceptibility. Upon doping, the spin-triplet superconducting regime can be further broadened in the parameter space and extends to larger $U$, offering a feasible mechanism to realize $p$-wave superconductivity in realistic cuprate chains.


Introduction
Strongly correlated materials, where the electronic structure cannot be approximated by the reductive band theory, have become a research frontier. In particular, two types of unconventional superconductivity have attracted considerable attention. One of them is the high-T c superconductivity discovered in cuprates [1]. Although this class of materials has been investigated for nearly 40 years, the pairing mechanism remains an enigma [2,3]. The other type of unconventional superconductivity is the topological triplet-pairing superconductivity [4][5][6], where electron fractionalizes into Majorana excitations [7,8] and is the foundation for topological quantum computing [9,10]. Therefore, pursuing such exotic superconductivity in realistic compounds constitutes a stimulating research topic.
The single-band Hubbard model, as the prototypical model carrying the strong correlation effects, has been widely em-ployed in the studies of many-body electron systems [11][12][13][14][15][16] as variants of this model are relevant to the two-dimensional (2D) cuprate superconductors. Besides, quasi-1D cuprate chains also constitute important class of strongly correlated materials that host intriguing correlated electron states and effects, e.g., the Tomonaga-Luttinger liquid (TLL) with spincharge separation [17][18][19]. On the other hand, most theoretical studies of the ground-state and dynamical properties [20][21][22][23] also lie in 1D as rigorous many-body simulations are more accessible using analytics [24], exact diagonalization, density matrix renormalization group (DMRG) [25] and quantum Monte Carlo [26][27][28]. Since both the on-site interaction U and near-neighbor (NN) interaction V correspond to the electronic repulsion at different distances, previous numerical studies focused on the cases with repulsive U, V > 0 as supposed relevant to real materials [29][30][31][32][33].
Most recently, a paradigm shift occurs as an in situ ARPES experiment on the 1D cuprate chain Ba 2−x Sr x CuO 3+δ (BSCO) has revealed an anomalously strong attraction V < 0 between NN electrons [34]. In contrast to the intrinsic electron-electron Coulomb repulsion, this attractive interaction is likely to be mediated by the strong electron-phonon coupling in transition metal oxides [35]. Such an effective attraction largely missed previously may serve as a key ingredient in both understanding the high-T c superconductivity and enabling exotic quantum phases in correlated materials [36][37][38][39][40][41][42]. Therefore, an interesting question naturally arises: Does such an effective attraction V help establish superconductivity pairing between the strongly correlated electrons?
To address this question, and also motivated by the recent experimental realization of such attractive-V extended Hubbard model (EHM, see Fig. 1a), we employ large-scale DMRG simulations and systematically explore its phase diagram. We especially focus on the possible realization of spin-triplet superconductivity while identifying all phases. At both half and quarter fillings, we have numerically determined the ground-state phase diagrams of the EHM, from which we identify a robust gapless TLL phase with a prominent spintriplet superconducting pairing (TS) with algebraic singularity. In two dimensions, the triplet superconducting (SC) state is topologically non-trivial where the fractional excitation can emerge on the boundary [8,[43][44][45][46]. However, quantum fluctuations are usually too strong in 1D such that interacting electrons in a Hubbard-type chain usually behave as a TLL, con-tradicting the mean-field and small-cluster predictions. Therefore in this paper, we refer this emergent TLL phase with divergent superconducting susceptibility to as a gapless TS phase.
Our main findings are summarized in Fig. 1. At half filling (see Fig. 1b), the TS phase survives only up to a finite U c /t 2.3 and is absent when U > U c . At quarter filling (see Fig. 1c), this TS phase extends to larger U s comparable to those in cuprates [34]. Between this TS phase and the regular PS phases with singly (PS 1 ) and doubly (PS 2 ) occupied clusters, we further identify an exotic PS x phase where the clustered electrons form the TLL and even TS states. With the model parameters determined from fitting dynamical data of BSCO, our study reveals a close proximity of this doped cuprate chain to the p-wave superconductivity, and provide theoretical guide for realizing such gapless TS phase in 1D cuprate chains.

Results
EHM with NN attraction. The BSCO chain can be described by the EHM with on-site U > 0 and NN attraction V < 0, whose Hamiltonian reads where c † iσ (c iσ ) is the electron creation (annihilation) operator, σ =↑, ↓ labels the electron spin, and n i = n i↑ + n i↓ is the particle number operator at site i. Throughout the study, we set hopping amplitude t = 1 as the energy unit, and focus on the ground state phase diagrams at both half and quarter fillings. In this work, we employ DMRG method with non-Abelian symmetry implemented [47,48] (see Methods and Supplementary Note 1).
To characterize various quantum phases, we compute the spin, charge, and pairing correlation functions. The spin-spin correlation is defined as F (r) = S i · S j , with S i(j) the spin operator at site i(j) and r ≡ j − i. The charge density correlation is defined as D(r) = n i n j − n i n j , where n i(j) is the particle number operator at site i(j). To characterize the superconducting pairing correlation, we consider both the spin-singlet (s-wave) pairing ,↓ for s = 1, 0, −1, respectively. Note that the EHM in Eq. (1) is SU(2) invariant so the above three components are degenerate in the spin-triplet channel, and we thus take the averaged Φ T (r) = 1 3 s Φ T,s (r) from our SU(2) DMRG calculations and compare it with Φ S .
Analytical results from the TLL theory. The TLL theory puts rigorous constraints [49][50][51][52] on our numerical results, which we always compare with and make use of in the analysis of our numerical data. In TLL, two-point correlation functions including the spin, charge and pairing correlations all de-  Fig. 1. Extended Hubbard model and phase diagrams. a illustrates the BSCO compound and corresponding extended t-U -V Hubbard model with NN hopping t, on-site repulsive U , and NN attractive V terms. b and c show the quantum phase diagrams of the EHM at half and quarter fillings, respectively. The solid black line represents the asymptotic phase boundary V = − U 2 − 8 ln 2 3U in the strong coupling limit, and the dashed line for V = −U/2 in the weak coupling limit [36,41]. The blue circle in c represents the parameters U = 8 and V = −1 of the doped 1D cuprate chain BSCO [34]. cay in power law ∼ r −α , with exponents α determined by two basic Luttinger parameters K σ and K ρ , respectively related to the spin and charge degrees of freedom (see more details in the Supplementary Note 2). To accurately evaluate these intrinsic parameters, one can calculate the momentum-dependent spin structure factor S m (k) and charge structure factor S c (k), and then extract K σ and K ρ .
For the current EHM in Eq. (1) with SU(2) spin symmetry, K σ = 1 for the spin density wave (SDW), TLL, and TS phases with gapless spin excitations, while K σ = 0 in the spin gapped phase PS 2 . Therefore, K ρ uniquely determines the power-law exponents α of various correlations: for charge and spin correlations there exist a uniform mode with exponent α 0 = 2 and a 2k F mode with α 2k F = 1 + K ρ ; for the pairing correlations Φ S and Φ T , they both have uniform modes with the same exponent α SC = 1 + 1/K ρ , which dominates over the spin and charge correlations when K ρ > 1. Consequently, the low-T behaviors of the staggered magnetic, charge, and pairing susceptibilities are also controlled by K ρ , i.e., χ SDW ∼ T Kρ−1 , χ CDW ∼ T Kρ−1 , and χ SC ∼ T 1/Kρ−1 . For K ρ > 1 or < 1, these susceptibilities exhibit apparently distinct behaviors as T → 0. Thus, the Luttinger parameter constitutes an essential quantity characterizing the underlying phases of a 1D system. In practice, we extract the Luttinger parameter K ρ via a second-order polynomial fitting of S c (k) in the small k regime [31,33,50,53] (see Supplementary Note 3 for details). To minimize the boundary effect, we evaluate the correlation functions using sites away from both ends.
Quantum phase diagram at half filling. We summarize our main findings at half filling in the phase diagram of Fig. 1b, where the SDW, phase separation PS 2 with doubly occupied sites clustered, and most remarkably, a TLL phase with promi-  nent superconductive pairing is uncovered. To show the distinction of these phases, we present simulations along two typical paths in Fig. 2, namely, the U = 1.6 and U = 4 vertical cuts in the phase diagram.
The Luttinger parameter K ρ clearly separates the U = 1.6 systems into three regimes. As the interaction strength increases to |V | > |V c | 1 (but smaller than the phase separation transition strength |V s |, which will be discussed later), in Fig. 2a1 there exists an intermediate regime with K ρ > 1. We also compute the central charge c by fitting the entanglement entropy (see more details in Supplementary Note 4), and from Fig. 2b1 c is found to change from c 1 to about 2 for |V c | < |V | < |V s |, confirming that the intermediate phase has both gapless spin and charge modes. On the other hand, also as shown in Fig. 2, for the U = 4 case K ρ remains small for all values of V and does not exceed 1 (see Fig. 2a2) and the central charge remains c = 1 (Fig. 2b2), showing the absence of such intermediate phase.
With further increase of the attractive interaction for either U = 1.6 or 4, the system eventually exhibits phase separation for |V | > |V s |. The critical strength V s dependent on U is shown in Fig. 1 (see the detailed estimation of V s in Supplementary Note 1). Specifically for the two selected cuts, we found V s −1.55 for U = 1.6 (see Fig. 2a1-c1) and V s −2.42 for U = 4 (see Fig. 2a2-c2). In such a PS state, the clustered part consists of doubly-occupied sites and no singularity can be observed in various correlations. Therefore, we denote it as PS 2 to distinguish from other PS phases discussed later.
Among these three phases in the U = 1.6 case (and for other interactions U < U c 2.3, c.f., Fig. 1b), we are particularly interested in the intermediate one due to the signature of triplet pairing. As evidenced by the charge correlation results in Fig. 3a, the charge gap is closed by the attractive V term, and the Luttinger parameter K ρ can be fitted to be greater than 1 (see the inset of Fig. 3a, and more details in Supplementary Note 3). According to the TLL theory, the superconductive paring decays r −αSC with the exponent α SC = 1 + 1/K ρsmaller than the algebraic exponent (1 + K ρ ) of both the charge and spin correlations when K ρ > 1and thus constitutes the dominant correlation in the charge-2e channel, with an algebraically diverging pairing susceptibility χ SC (T ) for low temperature T .
In the weak attraction regime |V | < |V c |, K ρ vanishes in the thermodynamic limit [54] and K σ = 1 due to the spin SU(2) symmetry. In Fig. 3b, an quasi-long range spin order with an algebraic exponent of α SDW = 1 appears, which has logarithmically diverging spin structure factor of S m (k = π) (see Fig. 2c1,c2 and the insets). This is well consistent with the SDW scenario with a finite charge gap and quasi-long range spin order (see Supplementary Note 5). On the other hand, for the intermediate phase in Fig. 2c1 S m (π) ceases to increase vs. L, as the 2k F mode spin correlation decays faster than ∼ r −2 shown in Fig. 3b, which reveals a non-diverging magnetic susceptibility and thus rather distinct magnetic properties from that of the SDW phase.
Gapless triplet superconducting phase. As shown in Fig. 3c,d, it can be observed that both the singlet-(Φ S ) and triplet-pairing (Φ T ) exhibit power-law decay behaviors, and the latter with p-wave pairing symmetry clearly dominates over the former with the s-wave pairing symmetry. This is clearly demonstrated in Fig. 4a, where the strengths of the two correlations Φ T (r) and Φ S (r) are compared at a fixed distance r = 20. Though two pairing correlations are comparable in the SDW regime, Φ T (r) clearly surpasses Φ S (r) once entering the intermediate-V phase: the latter turns to decreasing, while Φ T (r) keeps increasing and becomes over one order of magnitude greater than Φ S (r).
Such a dominance of the triplet pairing in the TS phase holds for different distances r other than the fixed distance r = 20 in Fig. 4a. This dominance is reflected in the spatial distribution of both pairing correlations in Fig. 3c,d. There we find Φ T firstly decay exponentially in the SDW phase (the blue dots), then exhibits power-law behaviors for |V c | < |V | < |V s | (the red dots), and decays again exponentially for |V | > |V s | (the grey dots). We notice there is virtually no uniform but only 2k F mode in Φ S , as reflected in the smooth curves Φ S (r) × (−1) r−1 in Fig. 3d. For the gapless TS phase where we are most interested in, the dominance of Φ T is reflected by the comparison of Figs. 3c and d: Φ T (r) decays slower   than r −2 , while Φ S (r) decays faster than r −2 (Fig. 3d). More quantitatively, the ratio between these two pairing correlations |Φ T (r)/Φ S (r)| scales in power law r Kρ−1 , since the leading scaling in Φ T and Φ S is 1/r 1+1/Kρ and 1/r Kρ+1/Kρ , respectively (see Supplementary Note 2). We present such a powerlaw scaling extracted from our DMRG simulations in Fig 4b. Therefore, in the intermediate regime the pairing correlation Φ T dominates over Φ S not only in magnitude but actually in long-distance scaling, making it a rather unique gapless TS phase.
When compared to the phase diagram obtained in Ref. 36, our DMRG results in Fig. 1b show some agreement on the existence of three phases, yet there are still noticeable differences. Particularly, our DMRG calculations identify the upper boundary of the TS phase in agreement with V = −U/2 obtained from the perturbation theory in the small U regime while it deviates from this line in the strong coupling regime. Consequently, in contrary to Ref. 36 where the TS phase was shown extending to infinite U , our results in Fig. 1b suggest it can only survive up to U c 2.3, located in a much narrower regime. On the other hand, when compared to more recent studies [55,56] where the phase diagrams are only schematic, here we pinpoint the numerically accurate phase boundaries with large-scale DMRG calculations and reveal the predominant triplet quasi-long range TS pairing relevant to the realistic cuprate chain BSCO, decades after such a TS instability was proposed [36,40].
Finite doping. Besides half filling, we have also explored the phases in the doped EHM systems. We first focus on the quarter filling, where the triplet pairing instability is approximately maximized, as will be discussed later. The extracted phase diagram is presented in Fig. 1c. Here, we select a cut along U = 4 and explain the properties of each phase in Fig. 5. Similar to half filling, the Luttinger parameter K ρ > 1 characterizes the intrinsic nature of the correlations and separates the U = 4 systems into four phases (see Fig. 5a). Particularly, for V < V c −0.8, we identified a TS regime following the same principle as half filling, manifested as enhanced triplet and singlet pairing correlations. Between these two correlations, we evaluated their ratio |Φ T /Φ S | and found its envelop increasing monotonically as |V | enhances and exceeding 1 for V < V c (see Fig. 5b), despite some oscillations with distance r. Note the two pairing correlations now show the same scaling at long distance. Importantly, the TS phase at quarter filling is significantly wider than that at half filling, particularly in the large U regime.
Besides the TS regime, there are three different inhomogeneous PS phases, i.e., PS 1 , PS x , and PS 2 in Fig. 1c, in the doped system. The real-space charge distributions n(i) are shown in Fig. 5c, from which we see that in the PS phases the electrons cluster with filling n = 1, 2 or x ∈ (1/2, 1]. To track the evolution among these PS phases when V changes, we pick the center of the system as a representative, which always lies in the filled domain in a PS state due to the open boundary, and extract n(i = L/2) for different U and V strengths in Fig. 5d. This filling density starts with n(L/2) = 0.5 (i.e., the TLL and TS phases) and deviates from the uniform quarter filling when |V | is stronger than certain transition value. As n(L/2) = x is not a fixed integer value but varies between 0.5 and 1, we denote this regime as PS x . For small U , like U = 2, the system jumps from PS x to PS 2 at a second transition point. In contrast, this transition is preceded by a third PS phase for large U > U c 2.3 (the same as that of half filling). Taking U = 4 as an example, PS x firstly transits into an n(L/2) = 1 phase (denoted as PS 1 ), and then jumps into PS 2 as |V | further increases. For the doped cases with filling factors other than 1/4, the quantum phase diagram is qualitatively similar to that of Fig. 1c. The phase boundaries of PS 1 and PS 2 actually remain intact for other doping since they reflect the local energy relation between singly and doubly occupied states. The quantum many-body states in the clustered part of the three PS phases -PS 1 , PS 2 , and PS x -only depends on the interaction parameters U and V .
The existence of the PS x phase was missed in early studies on the same model [36,40], and the distinct feature of PS x is the clustered electrons that constitute a TLL liquid with fractional filling. With x continuous tuned by V , the clustered part of PS x can also become close to half filling in terms of density, i.e., x = 1. Nevertheless, it is distinct from that in the PS 1 phase, as the clustered electrons in the latter form a charge gapped SDW instead of a gapless TLL. Even more interestingly, we can also identify a K ρ > 1 regime and significant TS pairing correlations in the clustered part of PS x , showing the existence of gapless TS cluster in (at least part of) the PS x phase (see more details in Supplementary Note 6).
TS pairing in the 1D cuprate BSCO. Although the phase boundaries, i.e., the critical strengths of V , are U -dependent, they can be determined analytically at quarter filling in the U → ∞ limit [57,58]. In this limit, the Luttinger parameter  identified as U 8 and V −1 [34]. Despite anomalously strong, the effective attraction V is still slightly below this threshold.
To search for TS in larger parameter space, we further explore the full doping dependence. To approximate the realistic materials, we fix U = 8 and three different values of V , and evaluate the Luttinger parameter K ρ for a wide range of doping. As shown in Fig. 6, the uniform TS phase characterized as K ρ > 1 can be realized only if |V | > 1.2 (and |V | 1.7 before PS x sets in), in order to exhibit prominent superconducting instability below 40% doping, the maximal accessible doping at current experimental conditions. Therefore, the doped BSCO resides on the boundary to a TS phase (as also indicated in Fig. 1c), and a slight reduction of on-site U or enhancement of near-neighbor attraction V may drive it into the TS phase -both can be achieved by manipulating the electron-phonon coupling either inside the crystal or via a substrate [35].

Discussion
Our simulation is based on the recently extracted attractive extended Hubbard model for 1D cuprate BSCO from experiments [34]. Although this newly demonstrated model and its parameters have been theoretically reproduced from the electron-phonon coupling [35], its impact on emergent phases, especially unconventional superconductivity phases, remains unknown. In this work, we employ DMRG -the method of choice for 1D correlated systems -to investigate the EHM with both on-site repulsive and near-neighbor attractive interactions. At both half and quarter fillings, we identify a prominent gapless TS phase with the p-wave pairing induced by the attractive interactions. Different from the long-range order (hidden) assumption in the context of mean-field theory, the p-wave superconducting order identified in this correlated 1D chain is quasi-long-ranged: the triplet pairing correlation Φ T decays as a power-law at long distance and presents as the dominant charge-2e excitations in the gapless TLL, and specially, at half filling it dominates over the singlet pairing Φ S also in large distance scaling. Such dominance results in divergent triplet superconductive susceptibility at low temperature. This phenomenon can be detected by the spectral depletion in ARPES or the Drude peak in optical conductivity, both of which are accessible for in situ synthesized quasi-1D materials.
As the experimentally extracted model parameters for cuprates [34] are close to, though not within, the TS phase identified in our simulations, our finding may motivate further investigation and manipulation of cuprates towards a p-wave topological superconductor. Couplings between the cuprate chains may open a charge gap and introduce edge modes that can be very useful in future quantum technologies. Due to the chemical and structural similarity between 1D and 2D cuprates, our results of the TS phase in the attractive EHM here shed light on and call for further many-body studies of the superconductivity in the EHM of higher dimensions [59][60][61][62].
Lastly, our conclusion on the cuprate chain can be extended to other related electronic materials. Considering the widely existing electron repulsion and electron-phonon coupling, this model with a repulsive U and an attractive V may also be applicable, as a low-energy approximation, for other transition-metal oxides. There different cuprate compounds and other materials may exhibit different microscopic parameters (U and V ) due to their distinct chemical environments, and the rich quantum phases revealed in the EHM model studies here may find their interesting materialization.

Methods
Density matrix renormalization group. We perform DMRG calculations with the charge U(1) and spin SU(2) symmetries implemented through the tensor library QSpace [47,48], and compute system sizes up to L = 512 to obtain the spin, charge, and superconductive correlations, etc, with high precision. In the calculations, we retain up to m * = 2048 multiplets, equivalent to m ≈ 4000 U(1) states, which render small truncation errors 10 −7 . We use the open boundary conditions as in conventional DMRG calculations. Due to the existence of attraction V , particularly near the PS phase one needs to introduce pinning fields at both ends and perform sufficient numbers of sweeps (even over 100 times) to fully converge the results, e.g., the charge distribution along the chain (see Supplementary Note 1).

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.

CODE AVAILABILITY
All numerical codes in this paper are available upon request to the authors. for their technical support and generous allocation of CPU time.

Supplementary Note 1. DMRG Techniques in Simulating the Extended Hubbard Model
Relieving the boundary effects by pinning fields. At half filling and especially near the PS 2 phase boundary, there is strong boundary effects in the DMRG calculations, which significantly lower the simulation efficiency. To relieve this problem and accelerate the calculations, we apply a pinning field term −V (n i − 1) 2 to both boundaries i = 1, L, where n i = n i↑ + n i↓ is the local particle number operator at site i. Here we set the strength of the pinning term proportional to V , so that the boundary effect can be well reduced for various cases with different attraction strengths.
Such pinning terms do help the DMRG to relieve the strong boundary effects and converge to the uniform charge distribution in the bulk. For example, in Supplementary Figure 1 we show the results of half-filled EHM at U = 3 and V = −1.8 (i.e., in SDW phase) with and without pinning. There, the entanglement entropy is defined as where ρ(l) is the reduced density matrix of a subsystem with length l. In Supplementary Figure 1(a), we see although the S E (l) distribution is strongly affected by boundary effects without pinning, such modulation can be "healed" by the pinning terms. An alternative way to control the boundary effects is to push the calculations to longer system sizes. In Supplementary Figure 1(b), we compute a very large system with L = 512, where S E (l) shows the expected dome shape deep in the bulk, despite that the boundary effect still penetrates into the bulk with over 50 sites.
According to the conformal field theory [63,64], for 1+1 dimensional critical systems, the entanglement S E should follow a linear scaling with the conformal distance π sin π(2l+1) 2(L+1) with L the total system size. In the SDW phase with only one gapless spin mode, we have central charge c = 1 that can be obtained by fitting the S E data. However, a naive analysis fails for the L = 128 case due to strong finite-size effects as shown by the blue symbol and lines in Supplementary Figure 1(d). When the pinning term is turned on, the S E (l) curve exhibits a usual dome shape and is consistent with central charge of c = 1 even for a moderate system size L = 128, as shown in Supplementary Figure 1(a) and (d). Moreover, the pinning field makes the charge density distribution n(i) more homogeneous (thus closer to half-filling throughout the chain), as shown in Supplementary Figure 1(c). By pushing the calculations to larger sizes, like L = 512 in Supplementary Figure 1(e) we plot S E vs.l and find it falls into a linear relation in the bulk (largel regime), also confirming the central charge c = 1 in the SDW phase.
In conclusion, the pinning fields on the boundaries of the attractive-V EHM can significantly reduce the boundary effects, making the results much more accurate and well-behaved for analysis. Fortunately, at other fillings the problem is less severe, and it is also quite hard to find a proper pinning strength for every U , V and doping. Therefore in the quarter-filling calculations we did not employ the trick.
Wave function initialization and DMRG sweeps. Since there are multiple phase separation phases and first-order phase transitions in the phase diagram Fig. 1 of the main text, it is good to start with a proper initial wave function for later DMRG sweeps, so as to avoid being trapped in local minima. Through careful numerical tests, we found the so-called PS 2 initialization constitutes a very good option in our practical calculations, where the particle number is doubly occupied (n = 2) in the middle and n = 0 towards the ends (c.f., the n sw = 0 lines in Supplementary Figure 2).
In Fig. 5c of the main text we have shown the charge density distribution n(i) at quarter filling, in various phases of the EHM. Specifically, in Supplementary Figure 2 we show how the initial PS 2 state evolves vs. DMRG sweeps n sw , which finally converges to the distribution that corresponds to the true ground state. For example, in Supplementary Figure 2(a) we find, for V = −1 in the TS phase, the initially clustered part -the n = 2 plateau -gradually smears out and expands to the whole system, recovering finally a virtually uniform distribution. For V = −1.7 in the PS x phase, the n = 2 plateau firstly shrinks into a n = 1 one and then to the final n = x plateau in the central of the chain, as shown in Supplementary Figure 2  Determination of the PS 2 Boundary. Here we show the details of determination of condensation transition boundary in the half-filling phase diagram in Fig. 1b of the main text. As discussed above, the PS 2 initialization works very well in practice. Therefore, with this initial state we perform the DMRG calculations of different V values, and obtain an estimation of the transition point V s for fixed U . However, when V is close to the condensation transition point V s , the convergence becomes quite slow. To determine V s more efficiently, we compare the energy e g (per site) versus V computed with two different initial states, i.e., the PS 2 and a uniform initial state. For example, in Supplementary Figure 3(a) we show the energy e g versus V , obtained with the two different initial states. We find in the TS phase the two e g values obtained with different initial states are in excellent agreement, while in the PS 2 phase they are different: the PS 2 initial state leads to a significantly lower energy. Therefore we extrapolate the energy curve from two sides of the transition through a linear fitting, and find the crossing point [the black dot in Supplementary Figure 3 can be remained in the analysis. In Supplementary Figure 3(b) we plot the estimated V s from different system sizes, and compare the two schemes, i.e., including or excluding the pinning term contributions in e g . We see the results of the two schemes approach each other as L is increased. By a second-order polynomial fitting and extrapolation to L = ∞, the two schemes give a consistent result. However, the V s estimated by subtracting the pinning term are found to converge faster, and in practice the L = 128 result is already precise enough with a relative error only about 10 −3 . Therefore, in the practical calculations we stick to such a fast scheme and use the L = 128 data to obtain an accurate estimation of V s .

Supplementary Note 2. Spin, Charge, and Pairing Correlations in the Tomonaga-Luttinger Liquid Theory
The Tomonaga-Luttinger liquid (TLL) theory can be used to describe a large family of one-dimensional (1D) quantum critical states [49][50][51][52]. Here we briefly recapitulate certain results from the TLL theory, relevant to the current work, on the spin, charge, and pairing correlations used in the analysis of our DMRG results in the main text. For the results listed below, we focus on the cases with spin SU(2) symmetry and have the Luttinger parameter K σ = 1 for the spin gapless states.
In general, there exist multiple modes for a given two-point correlation functions. For example, consider the charge densitydensity correlation D(r) = n i n j − n i n j with n i(j) the particle number operator at site i(j) and r ≡ j − i, and the spin correlation F (r) = S i · S j with S i(j) the spin operator at site i(j). Up to the first two dominant modes, both correlations have a similar form as [50,58,65] where A 1 and B 1 are model-dependent parameters. Generally, when K ρ < 1 the 2k F oscillation term dominates in both the charge and spin correlations, while for K ρ > 1 the uniform r −2 term becomes the leading one. These different modes are reflected in the singularities of the structure factor S(k). The uniform mode results in S(k) K ν |k|/π for k → 0, wherẽ K ν=ρ = K ρ for charge andK ν=σ = 3 4 K σ for spin; while for the 2k F mode, it leads to S(k) ∼ c 1 + c 2 |k ∓ 2k F | Kρ for k → ±2k F .
Next we consider the two superconductive pairing correlations. The singlet pairing correlation function is defined as , whose two leading algebraic modes take the scaling form [65,66] Φ S (r) = C 0 r 1+1/Kρ + C 1 where the logarithmic correction is neglected. Therefore, for K ρ < 1 the 2k F term dominates, while for K ρ > 1 the uniform term takes over and the overall pairing correlation decays algebraically with an exponent 1 + 1/K ρ < 2. On the other hand, the (averaged) triplet pairing correlation Φ T (r) = 1 3 s=±1,0 Φ T,s (r) scales as with the logarithmic corrections also omitted. Different from Φ S , here the 2k F term is much weaker and the uniform term is always playing a dominant role in long distance r. This is indeed what we have observed in our DMRG calculations of Φ T , as can be seen in Fig. 3 of the main text as well as Supplementary Figures 8 and 9 below. In addition, we note the coefficients A 1 , B 1 , C 0 , C 1 , etc. are model-dependent, and can take very different values or even be absent. For example, in the main text we have found Φ S only has 2k F mode in the half-filled TS phase, i.e. C 0 0, based on our accurate DMRG calculations.
Besides the above charge-2e correlations, the single-particle Green's function G(r) = σ c † iσ c jσ also shows a power-law behavior as r −1−α G with α G = 2 ν γ ν [67], where As mentioned above, for the EHM with SU(2) spin symmetry, we always have K σ = 1. Therefore, G(r) has the slowest decaying power −1 when K ρ = 1 (at the SDW-TS transition, TLL-TS crossover, etc). By taking a Fourier transformation of G(r), we obtain the occupation number distribution n(k) in momentum space. In the TLL phase, although in general n(k) is continuous at ±k F (except for the K ρ = 1 case), there is nevertheless singularity in the form of n(k) ∼ c 3 + c 4 |k ± k F | α G sign(k ± k F ) right at the Fermi vector k F .

Supplementary Note 3. Determination of the Luttinger Parameter K ρ
As discussed in Supplementary Note 2, the Luttinger parameter K ρ can be extracted from the charge structure factor, i.e., S c (k) K ρ |k|/π for k → 0. Here S c (k) is the Fourier transformation (FT) of the charge correlation function where D(i, j) = n i n j − n i n j . In the standard FT, the momentum k takes discrete value 0, 2π L , 4π L , ... , 2π(L−1) L . To collect more momentum data points for the purpose of fitting, we extend the value of k to the continuum in [0, 2π] (or equivalently [−π, π]). That is, we still compute S c (k) using Eq. (S6) but now the k can take the continuum of values in the Brillouin zone, which constitutes a natural and smooth interpolation method.
For the charge-gapless phase like TLL, the asymptotic behavior of S c (k) approaching k = 0 is linear with k, while in the charge-gapped phase like the SDW we have K ρ = 0 and the small-k scaling is instead quadratic. This offers us an efficient way to extract K ρ from the charge structure factor S c (k). In practice, to reduce the boundary effects, when we compute S c (k) the left-and right-most L E edge points are skipped and only the bulk correlation data are used. By doing so, we observe that the the determined K ρ results are not very sensitive to L E and thus numerically accurate, as shown Supplementary Figure 4(d).
At half filling, as the SDW and TS phases are respectively charge gapped and gapless, we employ the second-order polynomial fitting S c (k) = Ak 2 + K ρk + C, wherek = k/π is the (normalized) momentum. The intercept C is introduced as the discarded edge sites can break the particle number conservation, making S c (0) = 0. The S c (k) results for U = 1.6 at half filling are shown in Supplementary Figure 4(a). As the charge gap in the SDW phase is small in the weak coupling regime of U = 1.6, the change of S c (k) from quadratic to linear behaviors can be seen clearly if we zoom in into the small k regime in Supplementary  Figure 4(b). Since the quadratic region in SDW phase is narrow, we need to carefully choose the momentum range [k min , k max ] for our fittings. We set the lower bound of the fitting range to k min = 2π/L (to avoid finite-size gap effects), and different upper bounds k max ranging from 0.03π to 0.06π have also been chosen in practical fitting. We take average of the fitting results to provide a reliable estimate of K ρ with error bars. As shown in Supplementary Figure 4(c), indeed we observed that in the SDW phase the quadratic coefficient A is predominant while in the TS phase (V < V c −1) the linear coefficient K ρ dominates.
At quarter filling the charge sector is always gapless in the uniform TLL and TS phases, and thus we can employ a linear fitting S c (k) = K ρk + B throughout. Besides, we are also interested in the clustered regime in phase separation phase PS x and extract the Luttinger parameter K ρ of the clustered electrons. To obtain that, we select the bulk correlations measured on the central regime, with the criteria that the local particle density being very close to and at lease 99% of the density in the For the 1+1 dimensional quantum critical states described by the conformal field theory [63,64], the central charge c constitutes an important characteristic of the universality class of quantum criticality, which can be extracted by fitting the entanglement data by the scaling form wherel = ln 4(L + 1) π sin π(2l + 1) 2(L + 1) is the conformal distance.  problem, we use the OBC scaling form of S E [64] S E (l) = c 6 ln 4(L + 1) π sin π(2l + 1) 2(L + 1) To further reduce the boundary effects, we use the bulk L/2 sites for fitting. In Supplementary Figure 5 Here we provide more supportive data in the determination of quantum phase diagram of the EHM at half filling.
The U = 1.6 case. As discussed in Supplementary Note 2, the single-particle Green's function G(r) decays slowest as r −1 when K ρ = 1. Indeed, in Supplementary Figure 6 decays in power law and departs from the r −1 scaling as K ρ > 1.
Regarding the spin correlations, although there is always 2k F = π singularity in the spin correlations in both the SDW and TS phases, there are still clear distinctions. In Supplementary Figure 6(b) we see there is divergent behavior in the spin structure factor S m (k) at k = π for V V c . This is because the long-range spin correlation in the SDW phase scales as (−1) r /r, and S m (π) thus diverges as ln(L) in the thermodynamic limit. However, within the TS phase the exponent of the 2k F = π correlation is 1 + K ρ > 2, meaning the spin correlation decays faster than 1/r and corresponds to a non-divergent S m (π). In Supplementary Figure 6(b) we indeed observe these features, and in particular find there only a cusp in S m (π) in TS phase with V V c .
The U = 4 case. In the large U = 4 case, there exists only SDW and PS 2 but no intermediate TS phase. The transition point to PS 2 can be accurately determined by the method described in Supplementary Note 1 above (c.f., Supplementary Figure 3), and it is estimated as V s −2.42. In Supplementary Figure 7 we plot various correlation functions and structure factors throughout V s < V ≤ 0. Except the quasi-long-range AFM correlation with scaling (−1) r /r in Supplementary Figure 7 Figure 7(g,h)], all decay exponentially. In TLL, although there is no true Fermi point existing, the momentum-space distribution n(k) nevertheless exhibits singularity at ±k F [50][51][52]. In Supplementary Figure 7  Below we provide supportive data for determination of the ground-state phase diagram of EHM at quarter filling.
Various correlations and structure factors at quarter filling. In Supplementary Figure 8, we consider the case of U = 4 and show various correlation and structure factor results. Similar to the half-filling case in Supplementary Figure 7, here we again have K σ = 1 for the spin gapless cases. At the crossover point V c −0.8 we have K ρ = 1 and the G(r) decays algebraically with a power of −1 [c.f. Eq. (S5)]. Moreover, for V < V c it still decays in power-law but the exponent becomes smaller than −1 as K ρ > 1. In Supplementary Figure 8(b) the momentum-space particle number distribution n(k) shows singularity at ±π/4 (i.e., ±k F ) for different interactions V within the TLL and TS phases.
For the spin and charge correlations in Supplementary Figure 8(c-f), we find there exists weak 2k F singularity in the TLL phase, which becomes very weak (and even negligible) in the TS phases [c.f., Supplementary Figure 8(c,e)]. The 2k F singularity can be clearly seen in the spin structure factor results shown in Supplementary Figure 8(d), where the structural peak at ±π/2 rather prominent in the TLL phase becomes smeared out in the TS regime. This is due to the fact that the dominant correlations in both the spin and charge channels in the TS phase belong to the uniform k = 0 mode.
Lastly, in Supplementary Figure 8(g,h) we show the singlet and triplet pairing correlations, both of which exhibit algebraic behaviors in the TLL and TS phases. It can be noticed that Φ T gets significantly enhanced as |V | increases in the TS phase, where the power-law exponent exceeds −2 (meaning decaying more slowly), while the intensity of Φ S barely changes as |V | enhances.
Pairing correlations in the clustered region of PS x . In the main text, we have mentioned the clustered electrons in PS x constitute a gapless TLL state, which even has prominent TS pairing in a parameter regime not far from the uniform TS regime (c.f., the regime with K ρ > 1 in Fig. 5a of the main text). To be more specific, in Supplementary Figure 9 Supplementary Figure 9(a,b)]. Interestingly, we can further add some electrons to the system and make the whole system away from quarter filling, which, however, does not alter the electronic states in the clustered plateaus [see Supplementary Figure 9(b)]. When computing the pairing correlations in these wider plateaus, the boundary effects are further reduced, and the results in Supplementary Figure 9(c,d,e) are in practice computed in this manner. The results confirm that the properties of the clusters depend only on U ,V but not on the total filling factor.
In Supplementary Figure 9(c,d), we find both pairing correlations decay slower than r −2 in the electron clustered regime, and the triplet pairing Φ T again dominates. In Supplementary Figure 9(e) we show the ratio |Φ T /Φ S | directly, and find at long distances Φ T is significantly stronger than Φ S , confirming the TS nature of the gapless clustered electrons.