Lifshitz transition in the phase diagram of two-leg $t$-$J$ ladder systems at low filling

We use a combination of numerical matrix product states (MPS) and analytical approaches to investigate the phase diagram of the two-leg $t$-$J$ ladder in the region of low to intermediate fillings. We choose the same coupling strength along the leg- and rung-directions, but study the effect of adding a nearest-neighbor repulsion $V$. We observe a rich phase diagram and analytically identify a Lifshitz-like band filling transition, which can be associated to a numerically observed crossover from s-wave to d-wave like superconducting quasi-long range order (QLRO). Due to the strong interactions, the Lifshitz transition is smeared into a crossover region which separates two distinct Luttinger theories with unequal physical meaning of the Luttinger parameter. Our numerically exact MPS results spotlight deviations from standard Luttinger theory in this crossover region and is consistent with Luttinger theory sufficiently far away from the Lifshitz transition. At very low fillings, studying the Friedel-like oscillations of the local density identifies a precursor region to a Wigner crystal at small values of the magnetic exchange interaction $J/t$. We discuss analytically how tuning parameters at these fillings modifies the phase diagram, and find good agreement with MPS results.


I. INTRODUCTION
A particularly interesting system of strongly correlated electrons is the t-J model, which was originally introduced as a simplification for the strong coupling limit of the Hubbard model, as already realized in the 1970s by Spa lek et al. [1,2] and later in the context of cuprate superconductivity [3,4].It possesses a rich phase diagram, and is believed to be a basic model for the study of hightemperature superconductivity [5] (see, e.g., [6]).In the context of the latter, it has been subject to numerous studies, but despite these efforts, its lack of integrability [7] allows the phase diagram to be well known only for one-dimensional systems [8].More recently, the advent of ultracold polar molecules [9][10][11][12][13][14][15][16] on optical lattices [17] inspired a generalization of the original t-J-model with fully tunable interactions [18,19].While its experimental realization is an ongoing challenge, progress has been made to demonstrate that spin exchange can be realized in these setups [20][21][22], so that future investigations of the t-J model in such experiments is envisaged.On a chain, the phase diagram of variants of this model are well studied numerically, e.g., using matrix product states (MPS) [8,23].The full tunability of interactions and the long-range nature of the dipolar interactions in the polar molecule setups show that interesting modifications can be achieved, e.g., an enhanced superconducting phase [18,24,25], or topological SC [26].However, the step towards exploring two-dimensional systems remains a major challenge.Recent progress has been reported using iPEPS [27][28][29][30][31][32][33][34], and also by treating multi-leg ladder systems (see, e.g., [35][36][37][38][39][40][41][42]).In this context, the simplest non-trivial extension of the chain system is to treat twoleg ladder systems (see, e.g., [43][44][45][46][47][48][49][50][51][52]), which lately have been proposed to describe organic crystals, e.g., doped crystals of terphenyl [53].The phase diagram has been investigated in some detail, but with an apparent focus on higher densities.Here, we aim at complementing these studies by considering two directions of interest: i) investigate in more detail the behavior at low densities; ii) investigate the effect of tuning interaction strengths as previously studied in chain systems.We pursue this goal by combining field theory, analytical considerations at very low fillings, and from a detailed numerical study using state-of-the art MPS methods.One important aspect is the 2-band nature of the two-leg ladder systems, leading to band-filling transitions [54].We show that such Lifshitz-type transitions are relevant for understanding the phase diagram at low fillings.Notably, here it leads to a more complicated field theory interpolating between two known low-energy field theories.In particular, in the two field theories the dependence of observables on the Luttinger parameter K c differs, so that in the crossover regime it is unknown how precisely the observables depend on K c .
It is necessary to numerically compare in detail the behavior of observables such as the algebraic decay of correlation functions in order to map out the phase diagram.This is done in the following for two variants of the t-J model, which allows us to study the effect of adding a Coulomb-type repulsion.This has been found to enhance superconducting phases [25,55], and here we can investigate the interplay with the Lifshitz-type transition.These considerations are complemented by an analytical treatment at very low densities, far below the Lifshitz transition, which allows us to estimate the extension of the phases when tuning parameters of the model, which can be useful in the context of future polar molecule experiments.
The paper is organized as follows.In Sec.II A, we introduce the model and briefly discuss the two phase diagrams which are presented in Fig. 1.This is followed arXiv:2211.13065v1[cond-mat.str-el]23 Nov 2022 up in Sec.II B by a presentation of analytical considerations regarding an appearing Lifshitz transition and the low energy field theory in means of bosonization.In Sec.II C, we introduce the observables that we investigate to explore the ground state phase diagram.In Sec.III, we present details to how we obtain the phase diagrams from the DMRG data.In addition, in Sec.IV we discus how an analytical treatment of two electrons can be used to estimate the size of the superconducting phase in the zero density limit.We conclude with a summary in Sec.V.The appendices contain further details of the discussions and results presented in the main part of the paper.

II. OUR SETUP: MODELS, METHODS, OBSERVABLES, AND PHASE DIAGRAM
A. Variants of the t-J-Model on the Two-Leg Ladder Geometry In the following model, Eq. ( 1), we consider the usual t-J model [1,4,6,56], but allow for a variable nearestneighbor Coulomb repulsion V , which has been studied before (see, e.g., [18,25,55]).On the two-leg ladder systems treated by us it reads with S i,l S i+1,l − V 4 n i,l n i+1,l , Here, l = 1, 2 is the leg index, i labels the rung, σ the spin direction, P s projects out double occupancies, c i,l,σ is a fermionic annihilation operator acting on site (i, l) and spin-direction σ, n i,l = σ c † i,l,σ c i,l,σ is the occupation operator, and S i,l = c † i,l,σ σ σσ c i,l,σ /2.We set the lattice constant equal to unity and work in units in which ≡ 1.
In the following, we treat antiferromagnetic spin exchange J > 0. For V = 1, we obtain the usual t-J model as obtained from 2nd order perturbation theory in the strong coupling limit of the Hubbard model.
In the following we will call this case the regular t-J ladder.Furthermore, we treat the system with V = 0, which is obtained by adding the corresponding nearestneighbor Coulomb-repulsion to the original t-J model.This case will be referred to as the V = 0 t-J ladder.One important effect of adding a Coulomb-repulsion is that it suppresses phase separation, and can lead to enhanced superconducting phases [55].Here, we revisit its effect in the low-density regime of the two-leg ladder system.
In Fig. 1 we present our main results and show the phase diagrams of both variants as obtained using MPS and the analytical approaches discussed further below.As can be seen, the phase diagrams are quite similar to each other: both show a sequence from a gapless Luttinger-liquid (LL) phase with central charge c = 2 and dominant spin correlations (SDW) at small values of J/t to a c = 1 LL with finite spin gap and dominant pairing correlation functions.At even larger values of J/t, phase separation sets in.These features and sequences of phases are very similar to the findings in t-J chains, in particular at very low densities.At the density n = 0.5, a charge gap opens for small J/t and a charge density wave insulator (CDWI) is formed.The V = 0 t-J ladder has an enhanced SC phase, which has been observed also for the corresponding chain system [25] (note that throughout the manuscript we denote by SC a phase with dominant pairing correlation functions, i.e., a phase with SC quasi-long-range order).In particular, the size of the spin-gap region before phase separation is substantially increased.However, in contrast to the chain, where the sequence of phases is not altered for V = 0, here at a density around n = 0.3 an additional c = 2 LL phase with dominant CDW-correlations is realized.This raises the question, if and how changing V in further t-J -systems (e.g.broader ladder systems or in 2D) can lead to new features in the phase diagram.Both systems show at low densities and small values of J/t a precursor-region to a Wigner crystal, which is similar to the findings in Hubbard chains reported in Ref. 57.In the spin-gap region, the s-wave SC at very low densities is identified to be a gas of electron pairs, which was previously reported also for the chain systems [8].
In the following, we explain in some detail how these phase diagrams were obtained and describe the field theoretical treatment around the Lifshitz transition line, which we estimate using a simple Hartree-Fock ansatz.

B. Analytical considerations: Lifshitz transition
In this section we summarize the analytical expectations for our system, Eq. ( 1).We consider a Hubbard-Heisenberg ladder instead of the constrained t-J ladder, i.e. we consider with and all other terms as in Eq. ( 1).The Hubbard-Heisenberg ladder is equivalent to the t-J ladder only in the limit U → ∞.Following a similar strategy to the one successfully applied in the literature [48] for high densities, we here study the Hubbard-Heisenberg ladder at small and intermediate U/t using perturbation theory and bosonization.Clearly, there is no guarantee that our results would reproduce the U → ∞ limit of the Hubbard-Heisenberg model, and of course perturbation theory and Gutzwiller projection do not commute.

Lifshitz transition at the Hartree-Fock-level
The non-interacting tight binding Hamiltonian on a ladder with nearest neighbor hopping can be solved exactly by the basis transformation which decouples the system into two independent chains and brings the Hamiltonian into the form where a = 0, π.Each chain can independently be diagonalized resulting in dispersion relations Thus, a Lifshitz transition occurs when the Fermi energy hits E F = −t, which corresponds to an electron density of n = n i = 1/2 per site (quarter filling).Hence, based on the non-interacting estimate, there is a Lifshitz transition at the upper edge of the phase diagram presented in Fig. 1.We now discuss interaction induced shifts of the position of this transition as displayed in Fig. 1.Using a Hartree-Fock calculation (details are relegated to Appendix A) we obtain an effective quadratic Hamiltonian (6) Under the assumption that only the lower band is occupied we find and the dispersion Note that the Hubbard U does not enter t * or t * ⊥ , but merely leads to an overall energy shift in the spectrum.This fact is used as we push the Hartree-Fock calculation beyond the limits of its validity when we send U → ∞ to enforce the equivalence between Hubbard-Heisenberg and t-J models.By assumption of filling only the lower band we further employ k F = πn for the Fermi energy so that the Lifshitz transition is defined by the condition This condition is displayed as a solid red curve in Fig. 1.
Effectively, the spin-interaction increases the splitting by t * ⊥ of the bands.

Lifshitz transition beyond Hartree-Fock
Lifshitz transitions in interacting one-dimensional systems are non-trivially affected by the interactions in the system [54,[58][59][60][61].We here briefly review this physics from different standpoints: First, when the chemical potential is far below the upper band, the effect of interband interactions only leads to virtual processes renormalizing the Luttinger liquid (LL) in the lower band, which can be treated using perturbation theory.There is a typical energy scale E p where this perturbation theory breaks down, i.e. the single band LL physics is inapplicable when the chemical potential is closer than E p to the upper band.As we explain in more detail below, this energy scale can also be estimated from the immediate vicinity of the Lifshitz transition.Second, when the chemical potential is far in the two-band regime, one may study the multi-band system using perturbative RG.Again, these RG equations break down inside a window of size E p above the bottom of the upper band.Third, we now discuss the physics inside this window, concentrating on a chemical potential µ π → 0 − right below the Lifshitz transition.Then, the free two-particle propagator in the upper subband takes the form In the presence of an intrasubband interaction V within the upper band (but for the moment neglecting interband interactions), particles repeatedly scatter off each other resulting in an exact inverse two-particle Green's function, Fig. 2 a), Clearly, the second term always dominates for on shell excitations ω = k 2 /4m as soon as |µ π | < E p ≡ mV 2 /8.Differently said, the particles at the Lifshitz transition are always strongly interacting and form a so-called "impenetrable electron gas", in which the two-electron wave function has nodes at equal particle positions even for opposite spin.As a consequence, when the lower band is coupled to this impenetrable gas, pair tunneling of singlet Cooper pairs is strongly suppressed and it was shown that also spin-spin interband interactions are marginally irrelevant.It was found in Ref. [54] that the lower subband remains a LL and ultimately screens the single particle excitations in the upper subband thereby suppressing their quasi particle weight.It is thus reasonable to think of a Lifshitz transition of polarons instead of electrons.
The interaction V is the most relevant perturbation at the critical point [54] and defines the crossover scale E p .We estimate this coupling for the model (2), keeping in mind that we want to ultimately push our theory to the U → ∞ limit.To leading order, V ∼ U and we dropped weaker intersite interactions.The interaction constant V is screened and for µ π → 0 − only the filled subband can contribute to static screening.Using the Dzyaloshinski-Larkin theorem [62,63], according to which the random-phase-approximation (RPA) is exact for a linearized spectrum in 1D, we obtain the following effective interaction, Fig. 2 b), and v F,0 = 2t * sin(πn).This interaction enters the ladder resummation, Eq. ( 10), and we exploit that for on-shell twoparticle excitations in the upper band the static screening approximation is justified, i.e.Π(ω = k 2 /4m, k) 1/πv F,0 .Physically, this follows from the much faster dynamics in the lower as compared to the upper subband.In summary we find V πv F,0 and E p = π 2 t * sin 2 (πn)/2.
We also comment on the regime of partially, but dilute filling of the upper subband.It is theoretically harder to describe, yet the same crossover scale E p is believed to limit the strongly interacting regime for positive µ π , as well [54].Approaching µ π → E p from above, we estimate the screening of the bare interaction by Eq. ( 11) and Π(ω, k) ≈ 1/πv F,0 + 1/πv F,π which is determined by the Fermi velocities of both filled bands.
Using E p = mV 2 /8 both above and below the Lifshitz transition, the window of impeneatrable electron gas is bounded by These conditions are displayed as dashed red lines in Fig. 1.We remark in passing that above the Lifshitz transition πn = k π F + k 0 F , we will use this relationship in the remainder of the paper.

Low-energy field theory
Interacting one-dimensional fermionic systems are field theoretically suitably captured by means of bosonization [63] in the limit when the important energy scales are small with respect to the Fermi energy counted from the edges of a given band (x i = αi, where we have introduced the lattice constant α for clarity), i.e., Here, Φ a,σ (x), Θ a,σ (x), are conjugate fields which are slow on the scale of the lattice constant and k a F is the Fermi momentum in band a = 0, π.It is convenient to introduce bosonic fields in the charge and the spin channels c and s, respectively, so that, in bosonic language, the kinetic part of the action takes the form Here, K ζ,a is the Luttinger parameter and v ζ,a the Luttinger velocity, respectively, of the corresponding degree of freedom.In the non-interacting limit K ζ,a = 1 and, of course, the bosonic fields in the a = π subband only exist above the Lifshitz transition.
In the presence of interactions, the situation is more subtle as some of the bosonic degrees of freedom gap out.It is custumary to summarize these interacting phases by the label CnSm where n and m denote the number of gapless bosonic modes in the charge and spin channel, respectively [64].A means to efficiently extract the total number of gapless modes numerically is to measure the central charge c = n + m.Numerically, this can be done by analyzing the spatial behavior of the von Neumann entanglement entropy, which for systems with open boundary conditions is given by [65][66][67] where x is the position of the bipartition and d a nonuniversal constant.S(x) is easily computed by MPS [23].
However, the open boundary conditions (OBC) lead to additional oscillations in S(x), which can be understood in terms of the oscillations of the local kinetic energy [68].These are of the form [69] B(t(l, i) − t), where B is a free constant that needs to be fitted, t(l, i) = (l,i),(h,j) c † (h,j) c (l,i) is the local kinetic energy and t is the mean value of t(l, i) in the bulk.The resulting central charge is indicated in the phase diagrams of Fig. 1; more details of our procedure can be found in Appendix B.
We now review the analytical expectations, and first consider µ π < −E p .In this limit, the a = π band may be disregarded and the effective low-energy theory is a spinful interacting single band model in the band of a = 0 orbitals.The charge mode remains gapless with nonuniversal K c,0 .The spin mode is gapless with K 0,s = 1 for small J/t (a C1S1 phase with c = 2) but gaps out for sufficiently large spin interactions leading to a C1S0 phase.Technically this is defined by the condition that the interaction constant mixing chiral and antichiral spin currents changes sign [63,70].
Next, consider the regime µ π > E p .In this limit, the bosonization of the two-leg Hubbard model is justified for moderate U and weak coupling.Renormalization group calculations [48,64] predict that out of the four bosonic modes, only the field describing the total charge, remains gapless.This corresponds to a C1S0 phase with action with K c,+ and v +,c the corresponding LL parameter and velocity.
Finally, in the regime −E p < µ π < E p of the impenetrable electron gas around the Lifshitz transition, a direct bosonization approach is challenging.Indeed, in the regime 0 < µ π < E p the interaction energies are large as compared to the chemical potential µ π of the upper subband, while in the regime −E p < µ 0 < π, the perturbative inclusion of virtual interband processes breaks down.Yet, one may hope that qualitative aspects of the physics above/below the Lifshitz transition persist in the impenetrable crossover regime and this work constitutes a numerical exploration thereof.Indeed, our numerical results indicate that the central charge continues to be either 1 or 2, so that also in this intermediate regime the phases can be characterized as CnSm phases.However, this regime realizes a crossover region of two different effective field theories, Eqs. ( 15) and ( 18), respectively.Usually, extracting the LL parameter from the numerics is possible and the results can be used to characterize the different phases.However, in such a crossover region of two different low-energy field theories, this approach is not as straight-forward in the following way: as discussed further below, the structure factor follows the expectation from bosonization (in particular a linear behavior at small momenta k), so that the usual procedure allows one to extract a numerical value for the LL parameter in the charge sector.This is the value plotted in the phase diagrams of Fig. 1.However, following Tab.I, the linear coefficient of the structure factor has to interpolate between K c,0 and K c,+ of the field theories Eqs. ( 15) and ( 18) with an unknown function.
We remark in passing that the commensurate filling at density n = 1/2 is special as it allows for a charge density wave state with gapped charge sector (charge density wave insulator, CDWI).

C. Observables and analytical expectations
In this section, we introduce the calculated observables we used to determine the phase diagrams, based on the approach of Refs.8 and 25.

Global observables
We define the spin gap as where E 0 is the groundstate energy of a finite system with L lattice sites and N particles in the corresponding spin sector S z total .We extrapolate to the thermodynamic limit (TL) by keeping the density n = N/L fixed and taking L → ∞.
In a similar fashion also a charge gap can be defined via where µ + and µ − are the energies needed to add or remove a particle, respectively.However, since adding or removing one particle would necessarily change S z total by one, this quantity can be influenced by the behavior of the spin gap.Therefore, we define the charge gap by adding and removing two particles to keep S z total =0, i.e., Again, we perform a finite-size extrapolation to obtain the charge gap in the TL.Furthermore, we test for the appearance of phase separation (PS) by computing the inverse compressibility, ) where e 0 = E 0 /N is the groundstate energy per particle computed for a system with L lattice sites and with S z total = 0. We choose ∆n = 0.05, which gives a sufficient approximation to the derivative in Eq. ( 21).

Correlation functions
In order to determine which kind of quasi-long-range order (QLRO) is predominant we compute the correlation functions in the charge, spin, and Cooper channels.
First, we consider the correlations in total and relative charge density with n± i = ni,1 ± ni,2 .The corresponding structure factor is computed via where R is the total number of rungs and i, j denote the ith and jth rung, respectively.Similarly, we consider correlations of the rung spins with S z,+ = S z i,1 + S z i,2 .Finally, we study singlet pairing correlations of rung and leg pair creation/annihilation operators ∆ r S (i) = ), and the leg pairing operator Since our ladder model is a quasi one-dimensional system, we assume that the correlation functions are predicted by bosonized field theory, i.e., that they are of the form [63] in our numerical approach, the parameters A , B , α , β , k , and ϕ can be used as free parameters to be fitted (they are then supplemented by charge spin Cooper Summary of expected coefficients in the power-law Eq. ( 28) in the low and high density limit, respectively.All power-laws in the Cooper channels (rr,rl, ll ) are expected to be equal.The low-density limit displays both a C1S1 phase and a C1S0 phase (results for the latter being quoted in brackets).subscripts c , s , rr , ll , rl in the corresponding charge, spin, and Cooper channels).In the field theory, the values can be expressed using the corresponding LL parameter, see Tab.I and the detailed discussion in the next Sec.II C 3. Note that in Eq. ( 28), for the sake of simplicity we ignore any kind of logarithmic correction and take only one harmonic into account.If one of the pairing correlation functions has the slowest decay, we identify SC quasi-long-range order and call the system superconducting.On the other hand, if either the density or the spin correlations are dominant, we will speak respectively of charge density wave (CDW) and spin density wave (SDW) quasi-long-range order.
In the SC phases, we distinguish between s-wave and d-wave pairing by means of the relative sign of leg-leg and leg-rung correlators.However, we emphasize that in the ladder geometry the s-and the d-wave channels correspond to the same irreducible representation of the symmetries of the system and consequently mix.

Analytical expectation for correlation functions
The exponents entering these correlation functions are related to the Luttinger parameter in charge and spin sector, respectively.In table I we summarize the analytical expectations for this relationship.In the following we comment on the results.
First, we consider correlations in the charge channel, Eq. (22).From the definition of fermionic creation and annihilation operators in the band of 0/π orbitals, Eq. ( 3), it is evident that n− i = c † i,0 c i,π + H.c. is gapped throughout the phase diagram.In contrast the correlations of the total charge are gapless and we use the field theory expression for the total density [63] Because of spin-interactions, the operator for 4k F oscillations is generally expected to be effectively independent of Φ a,s [63].The bosonized expression for the total density has multiple implications for the density correlator Eq. ( 22), where we use the form given by Eq. ( 28).
To begin with, we concentrate on the low density limit µ π < −E p below the lower dashed line of Fig. 1.The first term in Eq. ( 29) leads to α c = 2 and A c = −K c,0 /π 2 .As long as the spin sector is gapless, both k = 2k 0 F and k = 4k 0 F oscillations are possible and decay with β c = K c,0 + 1 and β c = K c,0 , respectively.In contrast, when the spin sector is gapped and Φ σ,0 ∈ 2πZ/ √ 8 (a singlet superconductor), the cosine in the density operator orders and 2k 0 F oscillations dominate the result.The situation is different when The density-density correlators do not display powerlaw correlated oscillations, yet, the correlator of (n + i ) 2 does [48].
Next, we review the spin channel [63].In full analogy to the above said, the z component of the total magnetization is In the single band limit µ π < −E p , we expect two phases.
In the C1S1 phase, α s = 2 with A s = −1/π 2 and predominant k = 2k 0 F oscillations with β s = K c,0 + 1.When the spin sector is gapped and Φ σ,0 ∈ 2πZ/ √ 8, the spin correlator quickly vanishes on length scales larger than the inverse spin gap.This gapped behavior is expected to qualitatively persist to the large density limit µ π > E p where we expect a single bosonic mode in the total charge sector.
Finally, we review the Cooper channel, see also Appendix C. The rung and leg pair annihilators take the bosonized form where (−1) a = 1 ((−1) a = −1) for a = 0, (a = π).Note that the operator content of these bosonized expressions is the same, but the relative sign changes below and above the Lifshitz transition: For low densities, cos(k 0 ) > 0, while for large densities when both bands are populated cos(k a ) ∝ −(−1) a .This accounts for the transition from s-wave to d-wave pairing.
In summary, we expect correlation functions with A = 0 for all Cooper correlators (rr, ll, rl ).We only considered k = 0, for which β = 1/K c,0 + 1 in the C1S1 phase and β = 1/K c,0 in the C1S0 phase.At large densities far above the Lifshitz transition, there is only a C1S0 phase with gapless Φ c,+ mode.Again, A = 0 = k but β = 1/(2K c,+ ).In contrast to the low density limit, where K c,0 > 1 implies the dominance of superconducting correlations, for large enough filling the condition becomes K c,+ > 1/2.

Determining the Luttinger parameter from the charge structure factor
In field theory, Luttinger parameters are defined by means of the prefactor of the action in Eqs. ( 15) and (18).While the Luttinger parameter of the spin modes is K s = 1 whenever a gapless spin mode is present in our phase diagram, we here briefly summarize how we determine the non-universal Luttinger parameter in the charge sector.
As we find that throughout the phase diagram there is only one gapless charge mode, we numerically obtain the Luttinger parameter in charge space by studying the structure factor of the density correlation functions, Eq. ( 22).Using the form Eq. ( 28) and α c = 2, for small values of the momenta k this leads to In the one band (two band) regime Hence, the Luttinger parameter can be obtained by fitting Eq. ( 32) for k → 0 to the numerically obtained structure factor.The slope is then equal to −πA c , and applying the corresponding relation gives the Luttinger parameter.Throughout the paper, we use the relation for A low to obtain the value of K c .This is the value displayed in the phase diagrams of Fig. 1.However, as discussed in Sec.II B 3, in the crossover region around the Lifshitz transition we expect the numerical value of the slope to be given by an unknown mixture of K c,0 and K c,+ , so that its interpretation is more involved than in the respective high-or low-density case, and the so-obtained numerical value cannot directly be used to characterize the phase diagram.We therefore suppress the index + or 0 in K c , since this aspect becomes only relevant outside the crossover region.

Friedel oscillations and precursor Wigner crystal
In Ref. 57 Friedel-like oscillations were used to identify a precursor region to a Wigner crystal in low-filled Hubbard chains.The 4k 0 F = 4πn term, which has a similar origin as the 4k 0 F in Eq. ( 29), is interaction induced and is prominent at smallest density and values of J/t in Fig. 1.Here, we apply the same analysis to our t-J-ladder systems and use the bosonization prediction for the density in the one-band regime [57] n It is also possible to estimate the value of the Luttinger parameter by fitting Eq. ( 33) at very low densities.For the C1S1 phase, the results are consistent with the ones obtained from the structure factor.However, in the C1S0 phase the fits are much more difficult to control, so that we will not further discuss this approach.Fitting Eq. ( 33) allows us to identify the precursor region of the Wigner crystal, for which the ratio F 2 /F 1 is finite, but vanishes outside this region.An exemplary fit of the local electron density in the Wigner crystal regime as well as F 2 /F 1 as a function of J/t are presented in Appendix D. As can be seen there, we do not obtain a sharp transition, but a crossover region in which the value F 2 /F 1 gradually decreases.

D. Details on the MPS calculations
We used the MPS-code contained in the SymMPSpackage [71] to calculate the ground state energies, local observables, and correlation functions.Both variants of the model were calculated with open boundary conditions on systems with up to 200 lattice sites.We used 22 sweeps to ramp up to the maximal bond dimensions χ max = 2000, and afterwards continue with further sweeps until convergence is reached using this value of χ max (a sweep is going once through the lattice).This is done by setting the SymMPS control parameter for the ground state energy (which roughly corresponds to the absolute error in the ground state energy [72]) to 10 −10 .Usually, this threshold is reached after 30 − 50 sweeps, in more difficult cases we went up to 200 sweeps.Also in these more difficult cases, a SymMPS-error of the energy ∼ 10 −9 is obtained.In addition, the parameters are set such that a discarded weight of 10 −12 is obtained.However, in particular at the higher values of the filling treated by us (n ∼ 0.4 − 0.5), the entanglement in the system is so high that the discarded weight in the final sweep is ∼ 10 −6 .As a reference, it would be interesting to compare our numerical results to approaches using the full SU(2) symmetry of the systems, or by reducing the entanglement, e.g., using the mode transform discussed in Refs.73 and 74.In the context of our paper, the results shown at the higher values of n will have a larger numerical error, and we will show only the results, for which we have the highest confidence.

A. Global Observables
Fig. 3 shows examples for our results for the inverse compressibility as well as for the spin and charge gap of the regular t-J ladder after extrapolating to the TL as a function of J/t at densities n = 0.15, n = 0.30, and n = 0.50.The value of J/t at which the inverse compressibility κ −1 vanishes marks the onset of phase separation.For n = 0.50 the system undergoes the transition to phase separation at J/t ≈ 2.3 while it happens at J/t ≈ 3.0 and J/t ≈ 3.4, respectively, for n = 0.30 and n = 0.15.Note that in the phase separation region the numerics become unstable due to the high degeneracy of the ground state, so that we do not display results inside the PS region.
The spin gap opens exponentially slowly, which makes it difficult to estimate the exact position of the critical point.We make a conservative estimate for the numerical accuracy of the values of the spin gap in the TL to be ∼ 10 −3 (see, e.g., 75 for an estimate of the accuracy that can be obtained when determining a gap in spin ladders using MPS).In Refs.18 and 25, a similar threshold was used to identify the line at which the spin gap opens.In this way, we estimate at n = 0.50 the spin gap to open at J/t ≈ 0.7 while for n = 0.15 it opens at J/t ≈ 2.4.In the case of the intermediate densities n = 0.30, the spin gap opens at a later point J/t ≈ 2.8.The inset in Fig. 3b illustrates our approach: it shows the spin gap for n = 0.15 close to the critical point determined by the spin gap.The estimated error is marked by a graycolored interval around 0. We find the value of the spin gap exceeds the error estimate for J/t ≥ 2.4 which is, thus, our estimate for the critical point.
The charge gap is zero for almost all investigated densities and all values of J/t, within an accuracy of ∼ 10 −3 .Only the density n = 1/2 forms an exception for which a commensurate filling allows for a charge gapped, CDW ordered insulating phase, which we denote by CDWI in order to distinguish from the CDW LL phase.In the regular t-J model, this phase extends at n = 0.5 from J/t = 0 to J/t ≈ 1.5 while for the V = 0 case the phase ends at J/t ≈ 2.4.
In order to extrapolate our data into the TL, we performed a finite-size scaling by fitting a second-order polynomial to the numerical results for different system sizes L as a function of 1/L.Examples of such fits are shown in the insets of Fig. 3a an Fig. 3c.Thereby, we encountered an even-odd effect.Every time L/2 is an odd number, we artificially remove one electron to ensure that there are as many spin up electrons as there are spin down electrons, and our system stays in the subspace of zero magnetization.Interestingly, the systems sizes for which this was necessary show a different scaling behavior than the other systems.However, in both cases, the numerical results scale to the same value in the TL within an accuracy of 10 −3 , consistent with the discussion above.
We compute the central charge according to the procedure outlined in Sec.II C 3 and Appendix B for finite systems up to L = 200 sites with a subsequent finite-size extrapolation.We obtain a value for the central charge of c ≈ 2 in the C1S1 phases and c ≈ 1 in the C1S0 regions of the phase diagrams of both systems treated.This transition is in good agreement with the position of the line at which the spin gap opens.These findings are in agreement with the results for the spin and charge gaps described before.A more detailed discussion of these FIG. 3. Comparison of the inverse compressibility (top), the spin gap (middle) and the two particle charge gap (bottom) for the regular t-J ladder in the TL as function of J/t for the densities n = 0.15 (blue), n = 0.3 (orange) and n = 0.5 (green).The inset in the central panel displays the spin gap at n = 0.15 close to its opening point.It also shows the estimated numerical accuracy of 10 −3 , indicated by the gray region.The spin gap is considered to be finite if its value exceeds 10 −3 .The inset in the lower and upper panel show, respectively, an example for the finite size scaling of the charge gap and inverse compresibility.Red and blue dots denote the values with J/t = 1 for systems where the number of lattice sites divided by two is even or odd, respectively.findings and of the accuracy estimated by us is found in Appendix B.

B. Luttinger Parameter
Fig. 4a shows an example for the structure factor of N (|i − j|) + for different values of J/t at density n = 0.20.As it can be seen, the structure factor goes linearly to 0 (with an accuracy ∼ 10 −11 ) for k → 0. This is in agreement with the expectations from bosonization and can be used to estimate the value of the Luttinger parameter according to Eq. ( 32).Furthermore, one can observe that the structure factor develops a kink at 2k F when the system approaches the spin gapped phase.An exception to this behavior forms the structure factor for very small densities and low J/t where we observe a discontinuity in the derivative of the structure factor at 4k F , indicating the existence of such oscillations in the correlation func-tion and therefore in the particle density.As Ref. [57] argued, strong interactions between the electrons cause 4k F oscillations, which indicate the appearance of a precursor phase towards a Wigner crystal.
In Fig. 4b we compare the structure factors of the N ± (|i − j|) correlation functions.The structure factor N − (|i − j|) does not vanish for k → 0 and has a vanishing derivative at k = 0. Therefore, in contrast to N + (|i − j|), the N − (|i − j|) correlations contain contributions of gapped modes and support the results obtained by the central charge that at least one charge mode must be gapped.
The Luttinger parameter of the ungapped mode can be obtained by a linear fit of the structure factor for values close to k = 0 and multiplying the slope with π (compare Eq. ( 32)).For our analysis, we waive a finite-size scaling and only use the structure factor of the largest system size we investigated (i.e.L = 200) since the lattice size seems not to play an important role, as further discussed in Appendix D.Even though the structure factors in Fig. 4 seem completely linear close to k = 0, they still contain slight oscillations.This might be an artifact caused by finite-size effects.Therefore, to average out these oscillations, the fit must be performed through more than just the closest two values to k = 0 as one would expect.It turns out that using the six closest values to k = 0 reproduce the results in Ref. [49] for n = 0.25 very well, so that we used this for all other cases as well.An exemplary fit is shown in Fig. 4a.
Usually a Luttinger parameter larger than one indicates attractive interactions in the low-energy field theory and one expects SC.However, in the crossover region around the Lifshitz transition indicated by the dashed lines in Fig. 1, this is not reproduced.We associate this discrepancy to the complications discussed in Sec.II B 3, so that the characterization of the phase diagram is only possible by carefully analyzing the behavior of correlation functions and the aforementioned global observables.

C. Correlation Functions and Superconductivity
An important aspect to characterize the different phases is to determine the dominant correlation function, i.e., the algebraically decaying correlation function with smallest exponent.Based on the bosonization results discussed in Sec.II C 2, we use the fit function where we neglected logarithmic corrections, and fit the numerical results for the correlation functions for the total density (c), the total spin (s), for rung-rung pairing (rr), leg-leg pairing (ll), and rung-leg pairing (rl).Note that in Eq. ( 28) we only took one harmonic into account, while in Eq. ( 34) we allow for two oscillating contributions with different k-values.However, it turned out that for all pairing and for the spin correlations, the main contributions in the oscillatory part come from 2k F oscillations.If we keep all three terms, the fit routine becomes more unstable even though it mostly finds only 2k F contributions.Therefore, these correlations are only fitted with one of the two oscillatory terms, as in Eq. ( 28).For the fit of the total density correlation function, we used, however, all three terms if also 4k F contributions are relevant, i.e., if we were in the precursor Wigner crystal region.
We introduce the notation for the smallest (i.e.dominant) exponent of all correlation functions via If the total density or the total spin correlation functions dominate, we consider the system to be in a CDW or SDW phase, respectively.If, on the other hand, one of the pairing correlations dominates, we call the system superconducting.
Note that for the leg-leg pairing correlation functions we only calculated intra-leg correlations.Furthermore, in order to average out possible contributions from the boundary, we apply a running average by choosing a window of 6 sites around the center of the system.The average is then taken for each possible distance |i − j| over the six different initial positions An exemplary fit of the total density correlation function is shown in Fig. 5.As can be seen, it deviates from a purely algebraic decay at large distances.This can hint to logarithmic corrections, boundary or finite size effects, or difficulties in the convergence of the MPS and shows FIG. 6. Smallest exponent x ξ of the total density (orange), the total spin (blue), the rung-rung pairing (green), and the leg-leg pairing (red) correlation functions as a function of J/t for the regular t-J model.(a) results for density n = 0.10; (b) results for n = 0.45.The graphs where the system enters phase separation.In addition, the total spin correlation function stops where the system develops a spin gap.The black line denotes the value of J/t for which the Luttinger parameter Kc = 1, determined by the structure factor of the total density correlation function.The dashed red line indicates where the spin gap opens and the spin correlations start to decay exponentially.
that it is difficult to obtain the value of the exponents with a high accuracy.Typically, in order to minimize the effect of these artifacts, we restrict our fit to consider only the first 10-30 data points.As further discussed below, we estimate the relative error in our values of the exponents to be between 10-20%.Fig. 6 shows our results for the regular t-J model for x ξ as a function of J/t for different correlation functions at densities n = 0.1 and n = 0.45.We do not show the exponents for the rung-leg correlation functions since the fits are numerically unstable and give different qualitative behavior depending on the number of oscillatory terms we keep.This is further discussed in Appendix D, where we find the Fourier spectrum of these correlation functions not to contain sharp 2k F oscillations, in contrast to the other pairing correlation functions, so that is much more difficult to obtain meaningful fits using expression (34).Note that for n = 0.45 we show our results for x ll only in the C1S1 phase.For larger values of J/t, we find the fits to become unstable.This can be due to the lower accuracy obtained in this region of the phase diagram, as discussed in Sec.II D. Fig. 6a demonstrates that below the crossover region of the Lifshitz transition the expectations from bosonization are fulfilled: for K c < 1, the spin correlations decay slowest, while for K c > 1 the pairing correlations dominate, irrespective of the presence of a spin-gap.This is the same behavior as observed for chains at low fillings [8,25].In addition, the different pairing channels are degenerate.
In contrast, at n = 0.45 the value of J/t ≈ 1.7 at which K c = 1 and the value of J/t ≈ 1.3, at which the pairing correlations become dominant, clearly disagree.The pairing correlations become dominant as soon as the spin gap opens.This further illustrates that in the vicinity of the Lifshitz transition the value of K c as obtained by us cannot be used to directly characterize the phase diagram.For small J/t the non-oscillating term in the density correlation function dominates (α c ) at both values of the filling shown.According to Tab.I, one expects x c = 2.While for n = 0.45 this is obtained in good accuracy, at n = 0.1 we find x c ∼ 1.85, clearly deviating from the expected value.This allows us to estimate a relative error of ∼ 10% for this quantity, and in a conservative estimate we assume the relative error for all exponents shown in Fig. 6 to be of the order 20% The phase boundaries shown in Fig. 1 are determined using this approach, and hence are affected by a small error bar, which is not displayed there for the sake of clarity.
We distinguish between s-wave-like and d-wave-like pairing based on the discussion in Ref. 76.We find that in the superconducting phases, the rung-rung correlation functions are strictly positive.Thus, the sign of the rungleg pairing reveals if a rung and a leg pair have a relative phase.If P rl s has the same sign as P rr s (i.e.positive), we call the phase (extended) s-wave SC, while for different sign (i.e.negative values of P rl s ) we call the phase d-wave SC.This classification is inspired by the definitions of s-and d-wave SC on a square lattice.As can be seen in Fig. 7, the relative sign between P rr s and P rl s depends on the filling n.At low filling (n = 0.  distances, illustrating the crossover from s-wave like to d-wave like SC when increasing the densities.

D. Gas of electron pairs
The formation of a spin gap at very low densities can be explained by electrons forming a gas of free electron pairs.We follow the description of such a gas given in Ref. [8].Fig. 8 shows a comparison of the ground state energy per particle for different densities and for the exact solution of a single two-particle bound state (cf.Sec.IV).For n = 0.05 both energies match well in the spin gapped region.This suggests that at n = 0.05 the groundstate consists of a set of bound electron pairs that are strongly diluted.Note, however, that the dashed lines show the transitions in the zero-density limit and the deviations already indicate weak interaction corrections.The noninteracting pair gas picture clearly breaks down at higher densities.For n = 0.15, one can argue that the qualitative behavior of the ground state energy as a function of J/t is still consistent with the low energy picture but with an already strong energy shift.

IV. ZERO-DENSITY LIMIT
In order to obtain a more detailed analytical picture of the lower edge of the phase diagram, we revisit an energetic argument first put forward by Emery et al. [77].It relies on the exact expression for the energy of a two- FIG. 8. Groundstate energy per particle as a of J/t for the densities n = 0.05, n = 0.15, and n = 0.25 in a ladder with 200 sites as a function of J/t.The black line is the exact groundstate energy per particle of two particles on an empty infinite lattice.The blue diagonal shows the groundstate energy per particle of the Heisenberg part of the regular t-J ladder for n = 1.The left and right dashed lines mark respectively the onset of the spin gapped phase and phase separation predicted by the low density limit.The blue and green array blow the blue graph denotes the SDW and spin gap phase for n = 0.05 respectively.
particle bound state [78][79][80].The argument goes like this: Consider the ground state of two electrons in an infinite t-J system.For small J/t the electrons behave like free particles up to a value of J (1) c /t, above which a bound two-particle state with energy E B becomes energetically favorable.In the (very) dilute limit and under the assumption that bound states of a higher number of particles play no role, this yields the boundary from SDW to the "gas of electron pairs" at the bottom of the phase diagram in Fig. 1.The transition point to phase separation, J c /t, is obtained from comparing the energy of two electrons in a hole-free "island" to the singlet bound state.Since in the TL these islands will still be very large, one can use the average electron energy of a Heisenberg model ground state.However, similar to the treatment in the polar molecule quantum simulators [18,25], we allow in the following also for an XXZ-anisotropy, and introduce the anisotropy parameter α = J z /J.In one spatial dimension the bound state energy of the conventional t-J model is given by E 1d B = −J − 4t 2 /J.Adapting the 2d calculation by Lin [78] to the ladder geometry (cf.Appendix E) yields an implicit equation of the form where E s is the energy of a singlet and We obtain J For the conventional t-J model this yields J (1) c = 2t.We note that the same formula holds true in extended 1d and 2d systems.The transition point to phase separation, J c , is obtained from the comparison of E B with the energy of two electrons in the ground state of the hole-free model.In this case the kinetic term does not contribute and the V -term only shifts the energy.We therefore calculated the ground state energy of the spinexchange term numerically with MPS.The dark blue line in Fig. 8 shows the resulting energy of two electrons in the phase separation, while the black line is the energy of the two-particle bound state (36).Hence, J (2) c is the point where both lines intersect.For values of J between J (1) c and J (2) c , it is most favorable for the electrons to form bound pairs giving rise to the gas of electron pairs in the phase diagram.The numerical results for density n = 0.05 are consistent with this prediction but show already some renormalization of the J c .
In Figure 9 we plot the width J (2) of the superconducting region in the phase diagram for varying values of V and the anisotropy α.At α = 1 (left plot) the lower edges of the phase diagrams in Fig. 1 are well reproduced.One can clearly see that a low or negative value of V enhances the superconducting parameter range.In addition, a value 0 < α < 1 can also broaden it further.
V. SUMMARY Using a combination of numerical MPS and analytical methods, we investigate the ground state phase diagram of two variants of the two-leg t-J ladder for densities below n = 0.5.The difference between both variants is the strength of the nearest-neighbor Coulomb repulsion V , which is known from previous work on chains [18,25,55] to enhance SC phases.
We numerically compute the spin and charge gaps, and the inverse compressibility, which we extrapolate to the TL.In addition, we compute for finite systems with up to L=200 sites the central charge, the Luttinger parameter, and correlation functions.We determined the Luttinger parameter by calculating the density correlation structure factor.We complement this by an analytical treatment of the band-filling Lifshitz-type transition using a simple Hartree-Fock ansatz.At J/t = 0, the Lifshitz transition, where the Fermi surface changes from two to four Fermi points, happens at n = 0.5.Using the Hartree-Fock approximation, we estimated the position of the Lifshitz transition as a function of J/t.This bandfilling transition is meaningful for the phase diagram of both models, it manifests itself as a crossover between two effective field theories in which the physical meaning of the bosonic charge mode and of the Luttinger parameter is of different nature.Around this Lifshitz transition, we estimate the size of the crossover region using RPA.While a field theory at the Lifshitz transition [54] had been developed before, important predictions such as the crossover function of the Luttinger parameter are still unknown.
For example, at low densities, the numerically obtained Luttinger parameter is in good agreement with the behavior of correlation functions and the phase boundary or crossover between phases, as known for t-J chain systems [8,25].In contrast, while in the crossover region it is still possible to numerically obtain values for a Luttinger parameter, its behavior and phase boundaries are incompatible with standard bosonization theory.We associate this to the aforementioned unknown crossover between the two different field theories at large and low densities, respectively, which makes it difficult to obtain meaningful numerical results for the characterization of the phase diagram.Therefore, we analyze in detail the decay of different correlation functions in order to characterize the phase diagram.
For small J/t, both variants of the t-J ladder possess dominant SDW correlations in a gapless C1S1 phase.Within this phase, at low fillings and for small J/t, we identify a precursor of a Wigner crystal, similar to the findings in Hubbard chains [57].For larger values of J/t both systems develop a spin gapped phase (C1S0).This observation is supported by the computed value of the central charge, which drops from two to one once the C1S0 phase is entered.For the lowest densities, the spin gapped phase can be understood as a gas of free electron pairs.The system with V = 0, however, has an additional CDW LL phase in the vicinity of the Lifshitz transition before opening the spin gap.Such a significant change of the phase diagram was not observed in chains [18,25].This opens the question if and how changing V in higher dimensional t-J -systems (e.g.broader ladder systems or in 2D) can lead to new features in the phase diagram.
The crossover region from s-wave-like SC, in which the sign of rung-rung and rung-leg pairing correlators is the same, to a d-wave-like SC phase in which these correlators have different sign, falls into the crossover regime of the Lifshitz transition.This raises the question if in a more accurate treatment of the Lifshitz transition a connection between both could be seen.This could be investigated, e.g. by computing the electronic spectral functions of the model as a function of the filling, which will allow to determine the precise point of the band filling transition.
One key observation of this paper is that the superconducting phase and the spin gapped phase seem to coincide in the n → 0 limit.Setting V = 0 leads to an enhanced SC phase, whose width can be obtained in this limit in an exact calculation of the binding energy of two electrons and a comparison between the exact ground state energy and the ground state energy of a Heisenberg chain.In addition, we exemplarily show for the lower edge of the phase diagram how an XXZ-anisotropy in the spin exchange can modify this phase.Both, tuning V and the XXZ-anisotropy can be relevant for experiments with ultracold polar molecules on optical lattices, which should allow full tunability of the interactions in the t-J model.
It will be interesting to study the Lifshitz transition and its impact on the phase diagram when increasing the number of legs in ladder systems.There might be multiple crossovers at which the system may qualitatively change its behavior.Also, it will be interesting to study such broad ladder systems for the Hubbard model or for magnetically frustrated systems like ladders with nextnearest-neighbour interactions.
the expected values c = 2 and c = 1, respectively.However, in particular for the largest system sizes, a deviation from the linear extrapolation is seen.This can have different reasons, one of them being the convergence being more difficult to control for the ladder system than for chains, in which this approach was used to obtain the central charge with high precision.Note also that the oscillations complicate the fitting procedure; it would be preferable to use periodic boundary conditions, in which these oscillatory terms do not appear.However, this case is more difficult to control with MPS, in particular also for the ladder systems treated here.It would be desirable to improve the MPS approach, e.g. using the mode transform of Refs.73 and 74, in order to reach a higher accuracy in the vicinity of the Lifshitz transition.Here, a small error margin remains in our treatment.Therefore, our results in the TL do not show a sharp drop in the value of c from 2 to 1 at the phase boundary, but more a smooth transition.Also note the exponentially small opening of the spin-gap, which indicates a large length scale in the vicinity of the critical points.This also influences the accuracy in the determination of the value of c, and much larger system sizes would be needed to obtain the central charge with higher precision, which again is difficult to achieve using MPS for the present ladder systems, so that we refrain from doing so at this point.Fig. 12 displays the so-obtained values of the central charge for both phase diagrams shown in Fig. 1.Note the difference between the region in which c = 1 and the line of opening of the spin-gap; while we cannot fully rule out additional effects playing a role, the aforementioned difficulties in reaching a higher precision in the estimation of c make it plausible that throughout the phase diagrams in the TL either a value of c = 2 (C1S1 phase) or c = 1 (C1S0 phase) is realized, and that at the phase boundaries this value shows a sharp drop in the TL, in accordance with the expectations from the field theory of Sec.II B 3. tinuum fields of right and left movers.We omit any pair density wave, which implies that the antisymmetric component of the leg-pairing operator is dropped.Using (−1) a = 1 ((−1) a = −1) for a = 0, (a = π) we find ∆ r S (i) −i √ 2 a (−1) a R a,σ (x i )σ y L a,σ (x i ) (C4) ∆ l S (i) −i √ 2 a cos(k a )R a,σ (x i )σ y L a,σ (i) .Fig. 13 displays the Fourier transform of the three pairing correlation functions treated by us.As can be seen, the rung-leg correlator shows a smooth behavior around 2k F , while the other two correlators show a kink in the Fourier transform, indicating oscillations in real space with wave number 2k F .This indicates that further effects influence the behavior of the rung-leg correlator in real space, which make it difficult to use the ansatz Eq. (34) for obtaining meaningful values of the exponents for this quantity.Fig. 14 shows the same quantities as Fig. 3, but for the case V = 0.As can be seen, the overall behavior is similar, but the values of the phase transition points differ between both models.for the V = 0 t-J ladder in the TL as function of J/t for the densities n = 0.15 (blue), n = 0.30 (orange) and n = 0.50 (green).
3. Charge structure factor for different system sizes Fig. 15 shows a comparison of the results for the symmetric part of the density structure factor for system sizes L = 80, . . ., 200.As can be seen, finite size effects are not playing an important role for the case displayed.Therefore, we focus on the behavior of the structure factor for the largest system sizes treated by us.

Precursor Wigner crystal
Figure 16 shows an exemplary fit of the local electron density according to Eq. ( 33) in the Wigner crystal regime at n = 0.10 and J/t = 1 where 4k 0 F oscillations play an important role.Fig. 17 shows the ratio F 1 /F 2 i.e. the dominance of the 4k 0 F oscillations for n = 0.05 as a function of J/t.Appendix E: Some details on the determination of the two-particle bound states at very low densities The general calculation idea was outlined by Lin [78].It makes use of a wave function ansatz and the explicit solution of the stationary Schrödinger equation (note the missing brackets in Eqs. ( 9) to (11) in [78]!)The energy of the two-particle bound state is determined by an implicit equation of the form where r = 8 for a chain and r = 16 for a square lattice.
The integral I 0 (E) is defined (in the s-wave case for total momentum Q = 0) as follows, Although in our case we cannot assume the 90 • rotation symmetry, we can make a connection to the original calculation by taking the parity eigenbasis in the y-direction (corresponding to momenta 0 and π) and halving of the respective matrix element t.This yields and consequently Eq. (37).Furthermore, we find that in the case of the ladder r = 12 in (E2).

FIG. 1 .
FIG. 1. Ground state phase diagrams of (top) the regular and (bottom) the V = 0 two-leg t-J ladder (1) obtained using MPS and the Hartree-Fock approach to the band-filling transition of Sec.II B 1. SDW stands for a 2-channel Luttinger liquid (LL; with central charge c = 2 as indicated) with dominant spin-density-wave correlations; CDW stands for a LL with dominant charge-density-wave correlations; SC stands for a c = 1 LL with dominant pairing correlation functions.s-wave and d-wave indicate for the corresponding type of SC as discussed in Sec.II C. The precursor region to a Wigner crystal is identified via 4kF -contributions to Friedel-like density oscillations.PS stands for phase separation, identified by a diverging compressibility.At very low fillings, the s-wave SC phase is identical to a gas of free electron pairs.The bold red line denotes the occurrence of a Lifshitz transition according to the Hartree-Fock ansatz (9).The dashed red lines estimate the crossover region around the Lifshitz transition connecting the two different field theories at high and low densities, respectively.In this region, the Luttinger parameter Kc(orange circles indicate the line at which Kc = 1 as obtained from the charge structure factor, see Secs.II C 3, II C 4) cannot uniquely be determined using standard approaches.All boundaries between the different LL phases are estimated by directly comparing the exponents of the different correlation functions.The blue line denotes the opening of the spin gap.The black horizontal line at n = 0.5 indicates the opening of a charge gap inducing a symmetry broken CDW insulator (CDWI) phase.The two points at zero density are taken from the calculation in Sec.IV.

FIG. 4 .
FIG.4.Examples for the structure factor of the total density correlation functions N+ at n = 0.2 for the regular t-J ladder for systems with 200 lattice sites and different values of J/t.(a) structure factor of N+ for J/t = 1 (SDW phase), J/t = 2 (close to the phase transition), and J/t = 3 (deep in the superconducting phase).The black dashed line denotes the wave-number 2kF for n = 0.2, while the red solid line is an example of a linear fit to the values corresponding to the six smallest k-values.(b) comparison between the structure factor for the N+ and N− correlation functions for J/t = 1.

FIG. 5 .
FIG. 5. Exemplary fit of the total density correlation function N+ at n = 0.15 and J/t = 2.8

45 FIG. 7 .
FIG. 7. Comparison of the different pairing correlation functions at parameters n = 0.10 and J/t = 2.5 (a), n = 0.35 and J/t = 2.0 (b), and n = 0.45 and J/t = 2.0 (c) in the superconducting phase of the regular t-J model.
FIG. 9. Width J (2) c −J (1) c of the superconducting region in the phase diagram in the zero-density limit.a) anisotropy α = 1 (standard case in this paper); b) width as a function of both α and V .A reduced value of α increases the superconducting region at negative V .

=
−6t, i.e. the kinetic energy of two free electrons.AsI 0 (E B → −6t) diverges, we need to solve −12t 2 /E s = −E B = 6t.This results in E s != −2t.In the general case with anisotropy, the energy of a singlet is given by

JFIG. 13 .
FIG.13.Real part of the Fourier transformation of the rungrung, rung-leg and leg-leg pairing correlation functions for J/t = 2.0 and n = 0.30.In contrast to the rung-leg, the rung-rung and leg-leg correlation functions show a clear kink at 2kF .

FIG. 14 .
FIG.14.Comparison of the inverse compressibility (top), the spin gap (middle) and the two particle charge gap (bottom) for the V = 0 t-J ladder in the TL as function of J/t for the densities n = 0.15 (blue), n = 0.30 (orange) and n = 0.50 (green).

FIG. 15 .
FIG.15.Comparison between the structure factor of the symmetric density correlation function for different system sizes for the regular two-leg t-J ladder.Note that for L < 200 the SymMPS energy convergence was set to 10 −8 .

FIG. 16 .
FIG.16.Expectation value of local electron density for n = 0.10 and J/t = 0.1 on one leg as function of the rung position.The black line denotes a fit applied to the data enclosed within the two dashed lines.

FIG. 17 .
FIG. 17. Ratio between the amplitudes, via fits of the local electron density, of the 2kF (F1) and 4kF (F2) term in Eq. (33) as function of J/t.The dashed line denotes a threshold where the 4kF term becomes negligible.