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

Steffen Bollmann  Institut für Theoretische Physik, Georg-August-Universität
Göttingen, Friedrich-Hund-Platz 1, D-37077 Göttingen, Germany  Max Planck
Institute for Solid State Research, 70569 Stuttgart, Germany  Alexander
Osterkorn  Institut für Theoretische Physik, Georg-August-Universität
Göttingen, Friedrich-Hund-Platz 1, D-37077 Göttingen, Germany  Elio J. König
Max Planck Institute for Solid State Research, 70569 Stuttgart, Germany
Salvatore R. Manmana  Institut für Theoretische Physik, Georg-August-
Universität Göttingen, Friedrich-Hund-Platz 1, D-37077 Göttingen, Germany

###### Abstract

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łek 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 high-temperature 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 two-leg 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 [55, 25], 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.1, we introduce the model and
briefly discuss the two phase diagrams which are presented in Fig. 1. This is
followed up in Sec. II.2 by a presentation of analytical considerations
regarding an appearing Lifshitz transition and the low energy field theory in
means of bosonization. In Sec. II.3, 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.

Figure 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.2.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.3. The precursor region to a
Wigner crystal is identified via $4k_{F}$-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
$K_{c}$(orange circles indicate the line at which $K_{c}=1$ as obtained from
the charge structure factor, see Secs. II.3.3, II.3.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.

##  II Our Setup: Models, Methods, Observables, and phase diagram

###  II.1 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,
56, 6], but allow for a variable nearest-neighbor Coulomb repulsion $V$, which
has been studied before (see, e.g., [55, 18, 25]). On the two-leg ladder
systems treated by us it reads

$\mathcal{H}^{tJV}_{\text{ladder}}=\mathcal{H}^{t}+\mathcal{H}^{JV}_{\text{leg}}+\mathcal{H}^{JV}_{\text{rung}}$
(1)

with

$\displaystyle\mathcal{H}^{t}=$
$\displaystyle-t\sum_{i,l,\sigma}P_{s}\left[c^{\dagger}_{i,l,\sigma}c^{\phantom{\dagger}}_{i+1,l,\sigma}+\text{H.c.}\right]P_{s}$
$\displaystyle-t\sum_{i,\sigma}P_{s}\left[c^{\dagger}_{i,1,\sigma}c^{\phantom{\dagger}}_{i,2,\sigma}+\text{H.c.}\right]P_{s},$
$\displaystyle\mathcal{H}^{JV}_{\text{leg}}=$ $\displaystyle
J\sum_{i,l}\left[\vec{S}_{i,l}\vec{S}_{i+1,l}-\frac{V}{4}n_{i,l}n_{i+1,l}\right],$
$\displaystyle\mathcal{H}^{JV}_{\text{rung}}=$ $\displaystyle
J\sum_{i}\left[\vec{S}_{i,1}\cdot\vec{S}_{i,2}-\frac{V}{4}n_{i,1}n_{i,2}\right].$

Here, $l=1,2$ is the leg index, $i$ labels the rung, $\sigma$ the spin
direction, $P_{s}$ projects out double occupancies, $c_{i,l,\sigma}$ is a
fermionic annihilation operator acting on site $(i,l)$ and spin-direction
$\sigma$,
$n_{i,l}=\sum_{\sigma}c^{\dagger}_{i,l,\sigma}c^{\phantom{\dagger}}_{i,l,\sigma}$
is the occupation operator, and
$\vec{S}_{i,l}=c^{\dagger}_{i,l,\sigma}\vec{\sigma}_{\sigma\sigma^{\prime}}^{\phantom{\dagger}}c^{\phantom{\dagger}}_{i,l,\sigma^{\prime}}/2$.
We set the lattice constant equal to unity and work in units in which
$\hbar\equiv 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 nearest-neighbor 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.

###  II.2 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

$\mathcal{H}^{UJV}_{\text{ladder}}=\mathcal{H}^{t}_{0}+\mathcal{H}^{U}+\mathcal{H}^{JV}_{\text{leg}}+\mathcal{H}^{JV}_{\text{rung}}$
(2a)

with

$\displaystyle\mathcal{H}^{t}_{0}=$
$\displaystyle-t\sum_{i,l,\sigma}\left[c^{\dagger}_{i,l,\sigma}c^{\phantom{\dagger}}_{i+1,l,\sigma}+\text{H.c.}\right]$
(2b)
$\displaystyle-t\sum_{i,\sigma}\left[c^{\dagger}_{i,1,\sigma}c^{\phantom{\dagger}}_{i,2,\sigma}+\text{H.c.}\right],$
(2c) $\displaystyle\mathcal{H}^{U}$
$\displaystyle=U/2\sum_{i,l}n_{i,l}(n_{i,l}-1),$ (2d)

and all other terms as in Eq. (1). The Hubbard-Heisenberg ladder is equivalent
to the t-J ladder only in the limit $U\rightarrow\infty$. 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\rightarrow\infty$
limit of the Hubbard-Heisenberg model, and of course perturbation theory and
Gutzwiller projection do not commute.

####  II.2.1 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

$c_{i,0/\pi,\sigma}=\frac{1}{\sqrt{2}}\left(c_{i,1,\sigma}\pm
c_{i,2,\sigma}\right)$ (3)

which decouples the system into two independent chains and brings the
Hamiltonian into the form

$\mathcal{H}^{t}=-t\sum_{i,a,\sigma}\left[c^{\dagger}_{i,a,\sigma}c^{\phantom{\dagger}}_{i+1,a,\sigma}+\text{H.c.}\right]-t\sum_{i}(n_{i,0}-n_{i,\pi}),$
(4)

where $a=0,\pi$. Each chain can independently be diagonalized resulting in
dispersion relations

$\varepsilon_{0/\pi}(k)=-2t\cos(k)\mp t.$ (5)

Thus, a Lifshitz transition occurs when the Fermi energy hits $E_{F}=-t$,
which corresponds to an electron density of $n=\langle n_{i}\rangle=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

$\mathcal{H}^{t}=-t_{\|}^{*}\sum_{i,a,\sigma}\left[c^{\dagger}_{i,a,\sigma}c^{\phantom{\dagger}}_{i+1,a,\sigma}+\text{H.c.}\right]-t_{\perp}^{*}\sum_{i}(n_{i,0}-n_{i,\pi}).$
(6)

Under the assumption that only the lower band is occupied we find

$\displaystyle t_{\|}^{*}$ $\displaystyle=t+J\frac{\sin(\pi
n)}{8\pi}\left(3-{V}\right)$ (7a) $\displaystyle t_{\perp}^{*}$
$\displaystyle=t+J\frac{n}{8}\left(3-V\right)$ (7b)

and the dispersion

$\varepsilon_{0/\pi}(k)=-2t_{\|}^{*}\cos(k)\mp t_{\perp}^{*}.$ (8)

Note that the Hubbard $U$ does not enter $t_{\|}^{*}$ or $t_{\perp}^{*}$, 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\rightarrow\infty$ to enforce the equivalence between Hubbard-
Heisenberg and $t$-$J$ models. By assumption of filling only the lower band we
further employ $k_{F}=\pi n$ for the Fermi energy so that the Lifshitz
transition is defined by the condition

$-2t_{\|}^{*}\cos(\pi n)-t_{\perp}^{*}=-2t_{\|}^{*}+t_{\perp}^{*}.$ (9)

This condition is displayed as a solid red curve in Fig. 1. Effectively, the
spin-interaction increases the splitting by $t_{\perp}^{*}$ of the bands.

####  II.2.2 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
$\mu_{\pi}\rightarrow 0^{-}$ right below the Lifshitz transition. Then, the
free two-particle propagator in the upper subband takes the form
$\mathcal{D}_{\pi}^{(0)}(k,\omega)=-i\sqrt{m/4\epsilon^{+}}$, with
$\epsilon^{+}=\omega-k^{2}/4m+2\mu_{0}+i0$ (in our case, $2m=1/t_{\|}^{*}$).
In the presence of an intrasubband interaction $\mathcal{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),

Figure 2: Two-particle Green’s function in the upper subband. a) Diagrammatic
representation of Eq. (10). b) Diagrammatic representation of Eq. (11)

$[D_{\pi}(k,\omega)]^{-1}=[D_{\pi}^{(0)}(k,\omega)]^{-1}-\mathcal{V}.$ (10)

Clearly, the second term always dominates for on shell excitations
$\omega=k^{2}/4m$ as soon as $|\mu_{\pi}|<E_{p}\equiv m\mathcal{V}^{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 $\mathcal{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\rightarrow\infty$ limit. To leading order,
$\mathcal{V}\sim U$ and we dropped weaker intersite interactions. The
interaction constant $\mathcal{V}$ is screened and for $\mu_{\pi}\rightarrow
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),

$\mathcal{V}_{\rm eff}(\omega,k)^{-1}=\mathcal{V}^{-1}+\Pi(\omega,k),$ (11)

with $\Pi(\omega,k)=v_{F,0}k^{2}/\pi[(\omega+i0)^{2}-v_{F,0}^{2}k^{2}]$ and
$v_{F,0}=2t^{*}_{\|}\sin(\pi n)$. This interaction enters the ladder
resummation, Eq. (10), and we exploit that for on-shell two-particle
excitations in the upper band the static screening approximation is justified,
i.e. $\Pi(\omega=k^{2}/4m,k)\simeq 1/\pi v_{F,0}$. Physically, this follows
from the much faster dynamics in the lower as compared to the upper subband.
In summary we find $\mathcal{V}\simeq\pi v_{F,0}$ and
$E_{p}=\pi^{2}t_{\|}^{*}\sin^{2}(\pi 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
$\mu_{\pi}$, as well[54]. Approaching $\mu_{\pi}\rightarrow E_{p}$ from above,
we estimate the screening of the bare interaction by Eq. (11) and
$\Pi(\omega,k)\approx 1/\pi v_{F,0}+1/\pi v_{F,\pi}$ which is determined by
the Fermi velocities of both filled bands.

Using $E_{p}=m\mathcal{V}^{2}/8$ both above and below the Lifshitz transition,
the window of impeneatrable electron gas is bounded by

$-2t_{\|}^{*}\cos(\pi n)-t_{\perp}^{*}=-2t_{\|}^{*}+t_{\perp}^{*}\pm E_{p}.$
(12)

These conditions are displayed as dashed red lines in Fig. 1. We remark in
passing that above the Lifshitz transition $\pi n=k_{F}^{\pi}+k_{F}^{0}$, we
will use this relationship in the remainder of the paper.

####  II.2.3 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}=\alpha i$, where we have introduced the lattice constant
$\alpha$ for clarity), i.e.,

$\displaystyle c_{i,a,\sigma}\sim\sum_{\pm}e^{\pm
ik_{F}^{a}x_{i}-i(\pm\Phi_{a,\sigma}(x_{i})-\Theta_{a,\sigma}(x_{i}))}.$ (13)

Here, $\Phi_{a,\sigma}(x),\Theta_{a,\sigma}(x)$, are conjugate fields which
are slow on the scale of the lattice constant and $k_{F}^{a}$ is the Fermi
momentum in band $a=0,\pi$. It is convenient to introduce bosonic fields in
the charge and the spin channels $c$ and $s$, respectively,

$\Phi_{a,c/s}=\frac{\Phi_{a,\uparrow}\pm\Phi_{a,\downarrow}}{\sqrt{2}}$ (14)

so that, in bosonic language, the kinetic part of the action takes the form

$S_{\rm kin}=\sum_{\begin{subarray}{c}\zeta=c,s\\\
a=0,\pi\end{subarray}}\frac{1}{2\pi K_{\zeta,a}}\int d\tau
dx\;\left[\frac{\dot{\Phi}^{2}_{a,\zeta}}{v_{\zeta,a}}+v_{\zeta,a}{\Phi^{\prime}_{a,\zeta}}^{2}\right].$
(15)

Here, $K_{\zeta,a}$ is the Luttinger parameter and $v_{\zeta,a}$ the Luttinger
velocity, respectively, of the corresponding degree of freedom. In the non-
interacting limit $K_{\zeta,a}=1$ and, of course, the bosonic fields in the
$a=\pi$ 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 C$n$S$m$ 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]

$\mathcal{S}(x)=\frac{c}{6}\ln\left[\frac{L}{\pi}\sin\left(\frac{\pi
x}{L}\right)\right]+d,$ (16)

where $x$ is the position of the bipartition and $d$ a nonuniversal constant.
$\mathcal{S}(x)$ is easily computed by MPS [23]. However, the open boundary
conditions (OBC) lead to additional oscillations in $\mathcal{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)-\bar{t})$, where $B$ is a free
constant that needs to be fitted,
$t(l,i)=\sum_{\braket{(l,i),(h,j)}}\braket{c^{\dagger}_{(h,j)}c^{\phantom{\dagger}}_{(l,i)}}$
is the local kinetic energy and $\bar{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
$\mu_{\pi}<-E_{p}$. In this limit, the $a=\pi$ 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 $\mu_{\pi}>E_{p}$. In this limit, the bosonization
of the two-leg Hubbard model is justified for moderate $U$ and weak coupling.
Renormalization group calculations[64, 48] predict that out of the four
bosonic modes, only the field describing the total charge,

$\Phi_{+,c}=\frac{\Phi_{0,c}+\Phi_{\pi,c}}{\sqrt{2}},$ (17)

remains gapless. This corresponds to a C1S0 phase with action

$S_{\rm C1S0}=\frac{1}{2\pi K_{c,+}}\int d\tau
dx\;\left[\frac{\dot{\Phi}^{2}_{+,c}}{v_{+,c}}+v_{+,c}{\Phi^{\prime}_{+,c}}^{2}\right],$
(18)

with $K_{c,+}$ and $v_{+,c}$ the corresponding LL parameter and velocity.

Finally, in the regime $-E_{p}<\mu_{\pi}<E_{p}$ of the impenetrable electron
gas around the Lifshitz transition, a direct bosonization approach is
challenging. Indeed, in the regime $0<\mu_{\pi}<E_{p}$ the interaction
energies are large as compared to the chemical potential $\mu_{\pi}$ of the
upper subband, while in the regime $-E_{p}<\mu_{0}<\pi$, 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 C$n$S$m$ 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. 1, 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.

How does this affect observables? Schematically, the electronic operators are
$c_{i,\sigma}\sim\sum_{\pm}e^{i[\mp\Phi_{0,c}+\Theta_{0,c}]/\sqrt{2}}$ below
the Lifsitz transition, but
$c_{i,\sigma}\sim\sum_{\pm}e^{i[\mp\Phi_{+,c}+\Theta_{+,c}]/2}$ above the
transition. Thus, exponents stemming from
$\langle\Phi_{c}(x)\Phi_{c}(0)\rangle$ correlators, which are $\sim 1/K_{c,0}$
below the transition become $\sim 1/2K_{c,+}$ above the transition, while
those exponents stemming from $\langle\Theta_{c}(x)\Theta_{c}(0)\rangle$
correlators change from $\sim K_{c,0}$ to $K_{c,+}/2$. The latter correlators
however do not show up in the observables we study, see below.

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).

###  II.3 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, 25.

####  II.3.1 Global observables

We define the spin gap as

$\Delta_{S}=E_{0}(N,S^{z}_{\text{total}}=1)-E_{0}(N,S^{z}_{\text{total}}=0)\,,$
(19)

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}_{\text{total}}$. We extrapolate to the thermodynamic limit (TL) by
keeping the density $n=N/L$ fixed and taking $L\to\infty$.

In a similar fashion also a charge gap can be defined via

$\Delta_{C}=\mu_{+}-\mu_{-},$

where $\mu_{+}$ and $\mu_{-}$ are the energies needed to add or remove a
particle, respectively. However, since adding or removing one particle would
necessarily change $S^{z}_{\text{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}_{\text{total}}$=0,
i.e.,

$\displaystyle\Delta_{C}=$ $\displaystyle
E_{0}(N+2,S^{z}_{\text{total}}=0)+E_{0}(N-2,S^{z}_{\text{total}}=0)$
$\displaystyle-E_{0}(N,S^{z}_{\text{total}}=0).$ (20)

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,

$\kappa^{-1}=n^{2}\frac{\partial^{2}e_{0}}{\partial n^{2}}\approx
n^{2}\frac{e_{0}(n+\Delta n)+e_{0}(n-\Delta n)-2e_{0}(n)}{\Delta n^{2}}\,,$
(21)

where $e_{0}=E_{0}/N$ is the groundstate energy per particle computed for a
system with $L$ lattice sites and with $S^{z}_{\text{total}}=0$. We choose
$\Delta n=0.05$, which gives a sufficient approximation to the derivative in
Eq. (21).

####  II.3.2 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

$\displaystyle N(i,j)_{\pm}$
$\displaystyle=\braket{\hat{n}_{i}^{\pm}\hat{n}_{j}^{\pm}}-\braket{\hat{n}_{i}^{\pm}}\braket{\hat{n}_{j}^{\pm}},$
(22)

with $\hat{n}_{i}^{\pm}=\hat{n}_{i,1}\pm\hat{n}_{i,2}$. The corresponding
structure factor is computed via

$\mathcal{N}_{\pm}(k)=\frac{1}{R}\sum_{i,j=1}^{R}\text{e}^{ik(i-j)}N_{\pm}(|i-j|),$
(23)

where $R$ is the total number of rungs and $i,j$ denote the $i$th and $j$th
rung, respectively.

Similarly, we consider correlations of the rung spins

$S(i,j)=\braket{S^{z,+}_{i}S^{z,+}_{j}}-\braket{S^{z,+}_{i}}\braket{S^{z,+}_{j}}$
(24)

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
$\Delta_{S}^{r}(i)=\frac{1}{\sqrt{2}}\left(c_{i,1,\downarrow}c_{i,2,\uparrow}-c_{i,1,\uparrow}c_{i,2,\downarrow}\right)$,
and the leg pairing operator
$\Delta_{S}^{l}(i,l)=\frac{1}{\sqrt{2}}\left(c_{i,l,\downarrow}c_{i+1,l,\uparrow}-c_{i,l,\uparrow}c_{i+1,l,\downarrow}\right)$,

$\displaystyle P^{rr}_{S}(i,j)$
$\displaystyle=\braket{\left(\Delta^{r}_{S}(i)\right)^{\dagger}\Delta^{r}_{S}(j)},$
(25) $\displaystyle P^{ll}_{S}(i,j)$
$\displaystyle=\braket{\left(\Delta^{l}_{S}(i,l)\right)^{\dagger}\Delta^{l}_{S}(j,l)},$
(26) $\displaystyle P^{rl}_{S}(i,j)$
$\displaystyle=\braket{\left(\Delta^{r}_{S}(i)\right)^{\dagger}\Delta^{l}_{S}(j,l)}.$
(27)

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]

$\displaystyle
C(|i-j|)=\frac{A}{|i-j|^{\alpha}}+B\frac{\cos(k|i-j|+\varphi)}{|i-j|^{\beta}};$
(28)

in our numerical approach, the parameters
$A\,,B\,,\alpha\,,\beta\,,k\,,\text{and}\,\varphi$ can be used as free
parameters to be fitted (they are then supplemented by 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. 1 and the detailed discussion in the next
Sec. II.3.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.

####  II.3.3 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 1 we
summarize the analytical expectations for this relationship. In the following
we comment on the results.

| charge | spin | Cooper
---|---|---|---
$A_{\rm low}$ | $-K_{c,0}/\pi^{2}$ | $-1/\pi^{2}$ | $0$
$\alpha_{\rm low}$ | $2$ | $2$ | N/A
$\beta_{\rm low}$ |  $K_{c,0}+1$ ($K_{c,0}$) |  $K_{c,0}+1$ (gap) |  $1/K_{c,0}+1$ ($1/K_{c,0}$)
$A_{\rm high}$ | $-2K_{c,+}/\pi^{2}$ | $0$ | $0$
$\alpha_{\rm high}$ | $2$ | N/A | N/A
$\beta_{\rm high}$ | gap | gap | $1/2K_{c,+}$
Table 1: 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).

First, we consider correlations in the charge channel, Eq. (22). From the
definition of fermionic creation and annihilation operators in the band of
$0/\pi$ orbitals, Eq. (3), it is evident that
$\hat{n}_{i}^{-}=c^{\dagger}_{i,0}c^{\phantom{\dagger}}_{i,\pi}+\text{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]

$\displaystyle\hat{n}_{i}^{+}$ $\displaystyle\simeq
2n-\sum_{a=0,\pi}\Big{[}\frac{\sqrt{2}}{\pi}\Phi_{a,c}^{\prime}$
$\displaystyle+2n[e^{i(2k_{F}^{a}x_{i}-\sqrt{2}\Phi_{a,c})}\cos(\sqrt{2}\Phi_{a,s})+\text{H.c.}]$
$\displaystyle+2n[e^{i(4k_{F}^{a}x_{i}-\sqrt{8}\Phi_{a,c})}+\text{H.c.}]\Big{]}.$
(29)

Because of spin-interactions, the operator for $4k_{F}$ oscillations is
generally expected to be effectively independent of $\Phi_{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 $\mu_{\pi}<-E_{p}$ below the
lower dashed line of Fig. 1. The first term in Eq. (II.3.3) leads to
$\alpha_{c}=2$ and $A_{c}=-K_{c,0}/\pi^{2}$. As long as the spin sector is
gapless, both $k=2k_{F}^{0}$ and $k=4k_{F}^{0}$ oscillations are possible and
decay with $\beta_{c}=K_{c,0}+1$ and $\beta_{c}=K_{c,0}$, respectively. In
contrast, when the spin sector is gapped and $\Phi_{\sigma,0}\in
2\pi\mathbb{Z}/\sqrt{8}$ (a singlet superconductor), the cosine in the density
operator orders and $2k_{F}^{0}$ oscillations dominate the result.

The situation is different when $\mu_{\pi}>E_{p}$ – here only $\Phi_{c,+}$ is
gapless. While $\alpha_{c}=2$ persists, $A_{c}=-2K_{c,+}/\pi^{2}$. The
density-density correlators do not display power-law correlated oscillations,
yet, the correlator of $(\hat{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

$\displaystyle S_{i}^{z,+}$
$\displaystyle\simeq-\sum_{a=0,\pi}\Big{[}\frac{\sqrt{2}}{\pi}\Phi_{a,s}^{\prime}$
$\displaystyle-2n[e^{i(2k_{F}^{a}x_{i}-\sqrt{2}\Phi_{a,c})}\sin(\sqrt{2}\Phi_{a,s})+\text{H.c.}].$
(30)

In the single band limit $\mu_{\pi}<-E_{p}$, we expect two phases. In the C1S1
phase, $\alpha_{s}=2$ with $A_{s}=-1/\pi^{2}$ and predominant $k=2k_{F}^{0}$
oscillations with $\beta_{s}=K_{c,0}+1$. When the spin sector is gapped and
$\Phi_{\sigma,0}\in 2\pi\mathbb{Z}/\sqrt{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
$\mu_{\pi}>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

$\displaystyle\Delta_{S}^{r}(i)$
$\displaystyle\simeq\frac{-i}{\sqrt{2}}\sum_{a}(-1)^{a}e^{-i\sqrt{2}\Theta_{a,c}}\cos(\sqrt{2}\Phi_{a,s}),$
(31a) $\displaystyle\Delta_{S}^{l}(i)$
$\displaystyle\simeq\frac{-i}{\sqrt{2}}\sum_{a}\cos(k_{a})e^{-i\sqrt{2}\Theta_{a,c}}\cos(\sqrt{2}\Phi_{a,s}),$
(31b)

where $(-1)^{a}=1$ ($(-1)^{a}=-1$) for $a=0$, ($a=\pi$). 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})\propto-(-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
$\beta=1/K_{c,0}+1$ in the C1S1 phase and $\beta=1/K_{c,0}$ in the C1S0 phase.
At large densities far above the Lifshitz transition, there is only a C1S0
phase with gapless $\Phi_{c,+}$ mode. Again, $A=0=k$ but $\beta=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$.

####  II.3.4 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 $\alpha_{c}=2$, for small values of the momenta
$k$ this leads to

$\mathcal{N}(k)=-\pi A_{c}|k|.$ (32)

In the one band (two band) regime $A_{c}=-K_{c,0}/\pi^{2}$
($A_{c}=-2K_{c,+}/\pi^{2}$). Hence, the Luttinger parameter can be obtained by
fitting Eq. (32) for $k\to 0$ to the numerically obtained structure factor.
The slope is then equal to $-\pi A_{c}$, and applying the corresponding
relation gives the Luttinger parameter. Throughout the paper, we use the
relation for $A_{\rm 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.2.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.

####  II.3.5 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_{F}^{0}=4\pi n$
term, which has a similar origin as the $4k_{F}^{0}$ in Eq. (II.3.3), 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]

$\displaystyle n(x_{i})$ $\displaystyle=n-F_{1}\frac{\sin(2\pi
nx_{i})}{[\sin(\frac{\pi
x_{i}}{L+1})]^{\frac{K_{c,0}+1}{2}}}-F_{2}\frac{\sin(4\pi
nx_{i}-\phi)}{[\sin(\frac{\pi x_{i}}{L+1})]^{2K_{c,0}}}.$ (33)

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.

###  II.4 Details on the MPS calculations

We used the MPS-code contained in the SymMPS-package [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 $\chi_{\rm max}=2000$, and afterwards continue with further sweeps
until convergence is reached using this value of $\chi_{\rm 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 $\sim 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\sim 0.4-0.5$), the
entanglement in the system is so high that the discarded weight in the final
sweep is $\sim 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, 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.

##  III Numerical Results

###  III.1 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 $\kappa^{-1}$ vanishes marks the onset of phase separation.
For $n=0.50$ the system undergoes the transition to phase separation at
$J/t\approx 2.3$ while it happens at $J/t\approx 3.0$ and $J/t\approx 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 $\sim
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, 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\approx 0.7$ while
for $n=0.15$ it opens at $J/t\approx 2.4$. In the case of the intermediate
densities $n=0.30$, the spin gap opens at a later point $J/t\approx 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 gray-colored interval around 0. We find the value of the spin gap
exceeds the error estimate for $J/t\geq 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 $\sim 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\approx 1.5$ while for the $V=0$ case the phase ends at
$J/t\approx 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.

Figure 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.

We compute the central charge according to the procedure outlined in Sec.
II.3.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\approx 2$ in the C1S1 phases and $c\approx 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 findings and of the accuracy estimated by us
is found in Appendix B.

###  III.2 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 $\sim 10^{-11}$) for
$k\to 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 function 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_{\pm}(|i-j|)$
correlation functions. The structure factor $\mathcal{N}_{-}(|i-j|)$ does not
vanish for $k\to 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 $\pi$ (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.2.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.

Figure 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 $2k_{F}$ 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$.

###  III.3 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.3.2, we use the fit function

$D(|i-j|)=\frac{A}{|i-j|^{\alpha}}+B\frac{\cos(k_{1}|i-j|-\phi)}{|i-j|^{\beta}}\\\
+C\frac{\cos(k_{2}|i-j|-\psi)}{|i-j|^{\gamma}},$ (34)

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

$x_{\xi}=\text{min}(\alpha_{\xi},\beta_{\xi}\,,\gamma_{\xi}),\quad\xi\in\\{c,s,rr,ll,rl\\}\,.$
(35)

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
$i\in[L/2-3,L/2+2]$.

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 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%.

Figure 5: Exemplary fit of the total density correlation function $N_{+}$ at
$n=0.15$ and $J/t=2.8$ Figure 6: Smallest exponent $x_{\xi}$ 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 stop 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 $K_{c}=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.

Fig. 6 shows our results for the regular $t$-$J$ model for $x_{\xi}$ 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.4.

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\approx 1.7$ at which $K_{c}=1$ and
the value of $J/t\approx 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 ($\alpha_{c}$) at both
values of the filling shown. According to Tab. 1, one expects $x_{c}=2$. While
for $n=0.45$ this is obtained in good accuracy, at $n=0.1$ we find $x_{c}\sim
1.85$, clearly deviating from the expected value. This allows us to estimate a
relative error of $\sim 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 $\lesssim 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 rung-
leg pairing reveals if a rung and a leg pair have a relative phase. If
$P_{s}^{rl}$ has the same sign as $P_{s}^{rr}$ (i.e. positive), we call the
phase (extended) s-wave SC, while for different sign (i.e. negative values of
$P_{s}^{rl}$) 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_{s}^{rr}$ and $P_{s}^{rl}$ depends on the
filling $n$. At low filling ($n=0.1$), both have positive sign. At
intermediate filling ($n=0.35$), $P_{s}^{rr}$ remains strictly positive, but
$P_{s}^{rl}$ oscillates around zero. At higher filling ($n=0.45$), we see that
for large distances $|i-j|$ $P_{S}^{rr}$ and $P_{S}^{rl}$ have a different
sign for all distances, illustrating the crossover from s-wave like to d-wave
like SC when increasing the densities.

Figure 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.

###  III.4 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 non-interacting 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.

Figure 8: Groundstate energy per particle as a function 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.

##  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-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_{c}^{(1)}/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}^{(2)}/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,

$\hat{H}_{XXZ}=\frac{J}{2}\sum_{j}\left(S_{j}^{+}S_{j+1}^{-}+\text{H.c.}\right)+J_{z}\sum_{j}S_{j}^{z}S_{j+1}^{z}\,,$

and introduce the anisotropy parameter $\alpha=J_{z}/J$.

In one spatial dimension the bound state energy of the conventional $t$-$J$
model is given by $E_{B}^{\text{1d}}=-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

$\displaystyle\begin{split}-\frac{12t^{2}}{E_{s}}=\frac{1}{I_{0}(E_{B})}-E_{B},\end{split}$
(36)

where $E_{s}$ is the energy of a singlet and

$\displaystyle\begin{split}I_{0}(E)\quad&\overset{L\rightarrow\infty}{=}\quad-\frac{1}{2}\Big{(}\frac{1}{\sqrt{(E-2t)^{2}-(4t)^{2}}}\\\
&\qquad\qquad+\frac{1}{\sqrt{(E+2t)^{2}-(4t)^{2}}}\Big{)}.\end{split}$ (37)

We obtain $J_{c}^{(1)}$ upon equating $E_{B}\overset{!}{=}-6t$, i.e. the
kinetic energy of two free electrons. As $I_{0}(E_{B}\rightarrow-6t)$
diverges, we need to solve $-12t^{2}/E_{s}=-E_{B}=6t$. This results in
$E_{s}\overset{!}{=}-2t$. In the general case with anisotropy, the energy of a
singlet is given by
$E_{s}=\big{(}-\frac{1}{2}-\frac{\alpha}{4}-\frac{V}{4}\big{)}J$ yielding

$J_{c}^{(1)}=\frac{8t}{2+\alpha+V}\,.$ (38)

For the conventional $t$-$J$ model this yields $J_{c}^{(1)}=2t$. We note that
the same formula holds true in extended 1d and 2d systems. The transition
point to phase separation, $J_{c}^{(2)}$, 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
spin-exchange 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_{c}^{(2)}$
is the point where both lines intersect. For values of $J$ between
$J_{c}^{(1)}$ and $J_{c}^{(2)}$, 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}$.

Figure 9: Width $J_{c}^{(2)}-J_{c}^{(1)}$ of the superconducting region in the
phase diagram in the zero-density limit. a) anisotropy $\alpha=1$ (standard
case in this paper); b) width as a function of both $\alpha$ and $V$. A
reduced value of $\alpha$ increases the superconducting region at negative
$V$.

In Figure 9 we plot the width $J_{c}^{(2)}-J_{c}^{(1)}$ of the superconducting
region in the phase diagram for varying values of $V$ and the anisotropy
$\alpha$. At $\alpha=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<\alpha<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 [55, 18, 25] 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 band-filling 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\rightarrow 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
next-nearest-neighbour interactions.

###### Acknowledgements.

We acknowledge insightful discussions with R.M. Noack, Ö. Legeza, D.
Chakraborty, A. Schnyder, J. Mitscherling, P. M. Bonetti, D. Vilardi, R.
Scholle and V. Mishra. S.R.M. acknowledges financial support by ESI Vienna,
where part of this work was accomplished, and by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) - 217133147/SFB 1073,
project B03. A.O. acknowledges financial support by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) - 217133147/SFB 1073,
project B07.

## References

  * Spałek and Oleś [1977] J. Spałek and A. M. Oleś, Physica B+C 86-88, 375 (1977).
  * Chao  _et al._ [1977] K. A. Chao, J. Spałek, and A. M. Oleś, Journal of Physics C: Solid State Physics 10, L271 (1977).
  * Anderson [1987] P. W. Anderson, Science 235, 1196 (1987).
  * Zhang and Rice [1988] F. C. Zhang and T. M. Rice, Phys. Rev. B 37, 3759 (1988).
  * Bednorz and Müller [1986] J. G. Bednorz and K. A. Müller, Zeitschrift für Physik B Condensed Matter 64, 189 (1986).
  * Dagotto [1994] E. Dagotto, Rev. Mod. Phys. 66, 763 (1994).
  * Caux and Mossel [2011] J.-S. Caux and J. Mossel, Journal of Statistical Mechanics: Theory and Experiment 2011, P02023 (2011).
  * Moreno  _et al._ [2011] A. Moreno, A. Muramatsu, and S. R. Manmana, Phys. Rev. B 83, 205113 (2011).
  * Ni  _et al._ [2008] K.-K. Ni, S. Ospelkaus, M. H. G. de Miranda, A. Pe’er, B. Neyenhuis, J. J. Zirbel, S. Kotochigova, P. S. Julienne, D. S. Jin, and J. Ye, Science 322, 231 (2008).
  * Ospelkaus  _et al._ [2010] S. Ospelkaus, K.-K. Ni, D. Wang, M. H. G. de Miranda, B. Neyenhuis, G. Quéméner, P. S. Julienne, J. L. Bohn, D. S. Jin, and J. Ye, Science 327, 853 (2010).
  * Ni  _et al._ [2010] K. K. Ni, S. Ospelkaus, D. Wang, G. Quemener, B. Neyenhuis, M. H. G. de Miranda, J. L. Bohn, J. Ye, and D. S. Jin, Nature 464, 1324 (2010).
  * Carr  _et al._ [2009] L. D. Carr, D. DeMille, R. V. Krems, and J. Ye, New Journal of Physics 11, 055049 (2009).
  * Gadway and Yan [2016] B. Gadway and B. Yan, Journal of Physics B: Atomic, Molecular and Optical Physics 49, 152002 (2016).
  * Carr and Ye [2009] L. D. Carr and J. Ye, New Journal of Physics 11, 055009 (2009).
  * Wall  _et al._ [2014] M. L. Wall, K. R. A. Hazzard, and A. M. Rey, Quantum magnetism with ultracold molecules, in  _From Atomic to Mesoscale_ (World Scientific, 2014) Chap. Chapter 1, pp. 3–37.
  * Bohn  _et al._ [2017] J. L. Bohn, A. M. Rey, and J. Ye, Science 357, 1002 (2017).
  * Bloch  _et al._ [2008] I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80, 885 (2008).
  * Gorshkov  _et al._ [2011a] A. V. Gorshkov, S. R. Manmana, G. Chen, J. Ye, E. Demler, M. D. Lukin, and A. M. Rey, Phys. Rev. Lett. 107, 115301 (2011a).
  * Gorshkov  _et al._ [2011b] A. V. Gorshkov, S. R. Manmana, G. Chen, E. Demler, M. D. Lukin, and A. M. Rey, Phys. Rev. A 84, 033619 (2011b).
  * Yan  _et al._ [2013] B. Yan, S. A. Moses, B. Gadway, J. P. Covey, K. R. A. Hazzard, A. M. Rey, D. S. Jin, and J. Ye, Nature 501, 521 (2013).
  * Hazzard  _et al._ [2014] K. R. A. Hazzard, B. Gadway, M. Foss-Feig, B. Yan, S. A. Moses, J. P. Covey, N. Y. Yao, M. D. Lukin, J. Ye, D. S. Jin, and A. M. Rey, Phys. Rev. Lett. 113, 195302 (2014).
  * Li  _et al._ [2022] J.-R. Li, K. Matsuda, C. Miller, A. N. Carroll, W. G. Tobias, J. S. Higgins, and J. Ye, arXiv:2208.02216 (2022).
  * Schollwöck [2011] U. Schollwöck, Annals of Physics 326, 96 (2011).
  * Cheng  _et al._ [2015] C. Cheng, B.-B. Mao, F.-Z. Chen, and H.-G. Luo, EPL 110, 37002 (2015).
  * Manmana  _et al._ [2017] S. R. Manmana, M. Möller, R. Gezzi, and K. R. A. Hazzard, Phys. Rev. A 96, 043618 (2017).
  * Fazzini  _et al._ [2019] S. Fazzini, L. Barbiero, and A. Montorsi, Phys. Rev. Lett. 122, 106402 (2019).
  * Jordan  _et al._ [2008] J. Jordan, R. Orús, G. Vidal, F. Verstraete, and J. I. Cirac, Phys. Rev. Lett. 101, 250602 (2008).
  * Corboz  _et al._ [2010] P. Corboz, R. Orús, B. Bauer, and G. Vidal, Phys. Rev. B 81, 165104 (2010).
  * Verstraete  _et al._ [2008] F. Verstraete, V. Murg, and J. I. Cirac, Advances in Physics 57, 143 (2008).
  * Cirac and Verstraete [2009] J. I. Cirac and F. Verstraete, Journal of Physics A: Mathematical and Theoretical 42, 504004 (2009).
  * Corboz  _et al._ [2011] P. Corboz, S. R. White, G. Vidal, and M. Troyer, Phys. Rev. B 84, 041108 (2011).
  * Poilblanc  _et al._ [2014] D. Poilblanc, P. Corboz, N. Schuch, and J. I. Cirac, Phys. Rev. B 89, 241106 (2014).
  * Li  _et al._ [2021] J.-W. Li, B. Bruognolo, A. Weichselbaum, and J. von Delft, Phys. Rev. B 103, 075127 (2021).
  * Hubig  _et al._ [2020] C. Hubig, A. Bohrdt, M. Knap, F. Grusdt, and J. I. Cirac, SciPost Phys. 8, 021 (2020).
  * Dodaro  _et al._ [2017] J. F. Dodaro, H.-C. Jiang, and S. A. Kivelson, Phys. Rev. B 95, 155116 (2017).
  * Cheng  _et al._ [2018] C. Cheng, R. Mondaini, and M. Rigol, Phys. Rev. B 98, 121112 (2018).
  * Jiang  _et al._ [2018] H.-C. Jiang, Z.-Y. Weng, and S. A. Kivelson, Phys. Rev. B 98, 140505 (2018).
  * Bohrdt  _et al._ [2020a] A. Bohrdt, E. Demler, F. Pollmann, M. Knap, and F. Grusdt, Phys. Rev. B 102, 035139 (2020a).
  * Bohrdt  _et al._ [2020b] A. Bohrdt, F. Grusdt, and M. Knap, New Journal of Physics 22, 123023 (2020b).
  * Jiang  _et al._ [2021] S. Jiang, D. J. Scalapino, and S. R. White, Proceedings of the National Academy of Sciences 118, e2109978118 (2021).
  * Gong  _et al._ [2021] S. Gong, W. Zhu, and D. N. Sheng, Phys. Rev. Lett. 127, 097003 (2021).
  * Jiang  _et al._ [2022] S. Jiang, D. J. Scalapino, and S. R. White, arXiv:2206.07812 10.48550/ARXIV.2206.07812 (2022).
  * Tsunetsugu  _et al._ [1994] H. Tsunetsugu, M. Troyer, and T. M. Rice, Phys. Rev. B 49, 16078 (1994).
  * Dagotto and Rice [1996] E. Dagotto and T. Rice, Science 271, 618 (1996).
  * Troyer  _et al._ [1996] M. Troyer, H. Tsunetsugu, and T. M. Rice, Phys. Rev. B 53, 251 (1996).
  * Hayward and Poilblanc [1996] C. A. Hayward and D. Poilblanc, Phys. Rev. B 53, 11721 (1996).
  * Scalapino  _et al._ [2001] D. J. Scalapino, S. R. White, and I. Affleck, Phys. Rev. B 64, 100506 (2001).
  * White  _et al._ [2002] S. R. White, I. Affleck, and D. J. Scalapino, Phys. Rev. B 65, 165122 (2002).
  * Feiguin  _et al._ [2008] A. E. Feiguin, S. R. White, D. J. Scalapino, and I. Affleck, Phys. Rev. Lett. 101, 217001 (2008).
  * White  _et al._ [2015] S. R. White, D. J. Scalapino, and S. A. Kivelson, Phys. Rev. Lett. 115, 056401 (2015).
  * Musser and Senthil [2022] S. Musser and T. Senthil, arXiv preprint arXiv:2207.07128 (2022).
  * Lu  _et al._ [2022] X. Lu, D.-W. Qu, Y. Qi, W. Li, and S.-S. Gong, arXiv:2211.03778 10.48550/ARXIV.2211.03778 (2022).
  * Sous  _et al._ [2021] J. Sous, N. A. Gadjieva, C. Nuckolls, D. R. Reichman, and A. J. Millis, Nano Letters 21, 9573 (2021), pMID: 34761676\.
  * Meng  _et al._ [2011] T. Meng, M. Dixit, M. Garst, and J. S. Meyer, Phys. Rev. B 83, 125323 (2011).
  * Troyer  _et al._ [1993] M. Troyer, H. Tsunetsugu, T. M. Rice, J. Riera, and E. Dagotto, Phys. Rev. B 48, 4002 (1993).
  * Auerbach [2012] A. Auerbach,  _Interacting Electrons and Quantum Magnetism_, Graduate Texts in Contemporary Physics (Springer New York, 2012).
  * Söffing  _et al._ [2009] S. A. Söffing, M. Bortz, I. Schneider, A. Struck, M. Fleischhauer, and S. Eggert, Phys. Rev. B 79, 195114 (2009).
  * Takahashi [1971] M. Takahashi, Progress of Theoretical Physics 46, 1388 (1971).
  * Ogata and Shiba [1990] M. Ogata and H. Shiba, Phys. Rev. B 41, 2326 (1990).
  * Göhmann  _et al._ [1998] F. Göhmann, A. G. Izergin, V. E. Korepin, and A. G. Pronko, International Journal of Modern Physics B 12, 2409 (1998).
  * Essler  _et al._ [2005] F. H. L. Essler, H. Frahm, F. Göhmann, A. Klümper, and V. E. Korepin, in  _The One-Dimensional Hubbard Model_ (Cambridge University Press, Cambridge, 2005) pp. xi–xvi.
  * Dzyaloshinskii and Larkin [1973] I. E. Dzyaloshinskii and A. I. Larkin, Zh. Eksp. Teor. Fiz. 65, 411 (1973).
  * Giamarchi [2003] T. Giamarchi, in  _Quantum Physics in One Dimension_ (Oxford University Press, 2003) pp. 70–99.
  * Lin  _et al._ [1997] H.-H. Lin, L. Balents, and M. P. A. Fisher, Phys. Rev. B 56, 6569 (1997).
  * Holzhey  _et al._ [1994] C. Holzhey, F. Larsen, and F. Wilczek, Nuclear physics b 424, 443 (1994).
  * Korepin [2004] V. E. Korepin, Physical review letters 92, 096402 (2004).
  * Calabrese and Cardy [2004] P. Calabrese and J. Cardy, J. Stat. Mech.: Theor. and Exp. , P06002 (2004).
  * Laflorencie  _et al._ [2006] N. Laflorencie, E. S. Sørensen, M.-S. Chang, and I. Affleck, Phys. Rev. Lett. 96, 100603 (2006).
  * Roux  _et al._ [2008] G. Roux, S. Capponi, P. Lecheminant, and P. Azaria, The European Physical Journal B 68, 293 (2008).
  * Gogolin  _et al._ [2004] A. O. Gogolin, A. A. Nersesyan, and A. M. Tsvelik,  _Bosonization and strongly correlated systems_ (Cambridge University Press, Cambridge, England, 2004).
  * [71] S. Paeckel and T. Köhler, Scipal sym-mps, https://www.symmps.eu/.
  * [72] S. Paeckel, Oral communication.
  * Krumnow  _et al._ [2016] C. Krumnow, L. Veis, O. Legeza, and J. Eisert, Phys. Rev. Lett. 117, 210402 (2016).
  * Krumnow  _et al._ [2021] C. Krumnow, L. Veis, J. Eisert, and O. Legeza, Phys. Rev. B 104, 075137 (2021).
  * Michaud  _et al._ [2010] F. Michaud, T. Coletta, S. R. Manmana, J.-D. Picon, and F. Mila, Phys. Rev. B 81, 014407 (2010).
  * Noack  _et al._ [1994] R. M. Noack, S. R. White, and D. J. Scalapino, Phys. Rev. Lett. 73, 882 (1994).
  * Emery  _et al._ [1990] V. J. Emery, S. A. Kivelson, and H. Q. Lin, Physical Review Letters 64, 475 (1990).
  * Lin [1991] H. Q. Lin, Physical Review B 44, 4674 (1991).
  * Ogata  _et al._ [1991] M. Ogata, M. U. Luchini, S. Sorella, and F. F. Assaad, Physical Review Letters 66, 2388 (1991).
  * Hellberg and Manousakis [1995] C. S. Hellberg and E. Manousakis, Physical Review B 52, 4639 (1995).

##  Appendix A Analytical Phase diagram

In this appendix we present details on the interaction induced suppression of
the Lifshitz transition, as summarized in Sec. II.2.1 of the main text.

In all of the following, we will use the expansion in 0 and $\pi$ orbitals
accross the rungs, Eq. (3), to find

$\contraction{}{c^{\dagger}}{{}_{i,l,}}{c}c^{\dagger}_{i,l,\sigma}c_{i,l^{\prime},\sigma^{\prime}}=(-1)^{l+l^{\prime}}n\delta_{\sigma\sigma^{\prime}}/2,$
(39)

under the assumption that we are below the Lifshitz transition and thus
$\langle\sum_{\sigma}c^{\dagger}_{i,0,\sigma}c_{i,0,\sigma}\rangle=2n$ with
$n$ the density per site.

We start by studying the Hubbard U term, Eq. (2d). The calculation of the
Hartree-Fock self-energy leads to (normal ordering is taken before
contractions)

$\displaystyle\mathcal{H}^{\rm eff}_{U}$
$\displaystyle=U\sum_{i,l,\sigma,\sigma^{\prime}}:\contraction{c^{\dagger}_{i,l,\sigma}}{c}{{}_{i,l,\sigma}}{c}c^{\dagger}_{i,l,\sigma}c_{i,l,\sigma}c^{\dagger}_{i,l,\sigma^{\prime}}c_{i,l,\sigma^{\prime}}:$
$\displaystyle+U\sum_{i,l,\sigma,\sigma^{\prime}}:\contraction{}{c^{\dagger}}{{}_{i,l,}}{c}c^{\dagger}_{i,l,\sigma}c_{i,l,\sigma}c^{\dagger}_{i,l,\sigma^{\prime}}c_{i,l,\sigma^{\prime}}:$
$\displaystyle=\frac{Un}{2}\sum_{i,l,\sigma}c^{\dagger}_{i,l,\sigma}c_{i,l,\sigma}.$
(40)

We employ Einstein summation convention for spin indices. Note that we neglect
possible mean-field decoupling in the spin and pairing channels. Thus, the
Hubbard $U$ term only leads to an overall shift in the chemical potential for
both bands. This effect is trivial and thus disregarded in the main text.

Next, we study the term stemming from intrarung interactions

$\displaystyle\mathcal{H}^{\rm eff}_{\rm rung}$
$\displaystyle=\frac{J}{4}[\vec{\sigma}_{\sigma_{1},\sigma_{1}^{\prime}}\cdot\vec{\sigma}_{\sigma_{2},\sigma_{2}^{\prime}}-V\delta_{\sigma_{1},\sigma_{1}^{\prime}}\delta_{\sigma_{2},\sigma_{2}^{\prime}}]$
$\displaystyle\times\sum_{i}\Big{[}\left(:\contraction{c^{\dagger}_{i,1,\sigma_{1}}}{c}{{}_{i,1,\sigma_{1}^{\prime}}}{c}c^{\dagger}_{i,1,\sigma_{1}}c_{i,1,\sigma_{1}^{\prime}}c^{\dagger}_{i,2,\sigma_{2}}c_{i,2,\sigma_{2}^{\prime}}:+1\leftrightarrow
2\right)$
$\displaystyle+\left(:\contraction{}{c^{\dagger}}{{}_{i,1,\sigma_{1}}}{c}c^{\dagger}_{i,1,\sigma_{1}}c_{i,1,\sigma_{1}^{\prime}}c^{\dagger}_{i,2,\sigma_{2}}c_{i,2,\sigma_{2}^{\prime}}:+1\leftrightarrow
2\right)\Big{]}$
$\displaystyle=-\frac{JVn}{4}\sum_{i,l}c^{\dagger}_{i,l,\sigma}c_{i,l,\sigma}$
$\displaystyle+\frac{Jn}{8}\left[3-V\right]\sum_{i}\left(c^{\dagger}_{i,1,\sigma}c_{i,2,\sigma}+\text{H.c.}\right).$
(41)

Apart from yet another contribution to the overall chemical potential, this
result is the derivation of Eq. (7b).

Finally, we evaluate nearest neighbor intrachain interactions,

$\displaystyle\mathcal{H}^{\rm eff}_{\rm intra}$
$\displaystyle=\frac{J}{4}[\vec{\sigma}_{\sigma_{1},\sigma_{1}^{\prime}}\cdot\vec{\sigma}_{\sigma_{2},\sigma_{2}^{\prime}}-{V}\delta_{\sigma_{1},\sigma_{1}^{\prime}}\delta_{\sigma_{2},\sigma_{2}^{\prime}}]\sum_{i,l}$
$\displaystyle\Big{[}\left(:\contraction{c^{\dagger}_{i,l,\sigma_{1}}}{c}{{}_{i,l,\sigma_{1}^{\prime}}}{c}c^{\dagger}_{i,l,\sigma_{1}}c_{i,l,\sigma_{1}^{\prime}}c^{\dagger}_{i+1,l,\sigma_{2}}c_{i+1,l,\sigma_{2}^{\prime}}:+i\leftrightarrow(i+1)\right)$
$\displaystyle+\left(:\contraction{}{c^{\dagger}}{{}_{i,l,\sigma_{1}}}{c}c^{\dagger}_{i,l,\sigma_{1}}c_{i,l,\sigma_{1}^{\prime}}c^{\dagger}_{i,l,\sigma_{2}}c_{i,l,\sigma_{2}^{\prime}}:+i\leftrightarrow(i+1)\right)\Big{]}$
$\displaystyle=-\frac{JVn}{4}\sum_{i,l}c^{\dagger}_{i,l,\sigma}c_{i,l,\sigma}$
$\displaystyle-\frac{J}{4}\mathcal{C}\left[3-{V}\right]\sum_{i,l}\left(c^{\dagger}_{i,l,\sigma}c_{i+1,l,\sigma}+\text{H.c.}\right),$
(42)

where

$\displaystyle\mathcal{C}$ $\displaystyle=\frac{1}{4}\sum_{l,\sigma}\langle
c^{\dagger}_{i,l,\sigma}c_{i+1,l,\sigma}\rangle=\frac{1}{2}\int_{-k_{F,\pi}}^{k_{F,\pi}}\frac{dk}{2\pi}e^{ik}$
$\displaystyle=\frac{\sin(k_{F,\pi})}{2\pi}=\frac{\sin(\pi n)}{2\pi}.$ (43)

This is the origin of Eq. (7a) of the main text.

Figure 10: Exemplary fit of the entanglement entropy in the regular $t$-$J$
ladder at $n=0.10$ and $J/t=2.0$. The Plot only shows the entanglement entropy
for bisections between the rungs.

##  Appendix B Entanglement entropy and central charge

Since we treat systems with OBC, for which the MPS is more efficient, we need
to take into account oscillatory terms in the entanglement entropy. A typical
example is shown in Fig. 10, together with the fitting as explained in Sec.
II.2.3. In the shown case, a good quality of the fit is obtained for the
largest system size treated by us ($L=200$), and the numerical value of the
central charge is $c\approx 1.9$, which is close to the expected value $c=2$.

Figure 11: Examplary finite size scaling for the central charge in the regular
$t$-$J$ ladder. The right panel presents a finite size scaling in the C1S1
phase while the left panel shows a finite size scaling in the C1S0 phase.

In Fig. 11 we present a finite size extrapolation of the obtained values for
$c$ in the C1S1 and the C1S0 phases, respectively. We see that in both cases
finite-size effects are rather small, and that an extrapolation leads closer
to 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, 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.2.3.

Figure 12: Colormap plot of the central charge in the regular (top) and $V=0$
(bottom) two-leg $t$-$J$ ladder. The red line denotes where the spin gap
becomes finite while the black line marks the onset of phase separation. In
the phase separated region, the numerical data was set to zero. The filling
$n=0.05$ is left out since the fits were difficult to control.

##  Appendix C Fourier Transformation of Pairing Correlation Functions

In this appendix we summarize details about the pairing correlation functions
quoted in Sec. II.3.3.

The rung pair annihilator is

$\displaystyle\Delta_{S}^{r}(i)$
$\displaystyle=\frac{-i}{\sqrt{2}^{3}}\left(c_{i,1}^{T},c_{i,2}^{T}\right)\sigma_{y}\left(\begin{array}[]{cc}0&1\\\
1&0\end{array}\right)\left(\begin{array}[]{c}c_{i,1}\\\
c_{i,2}\end{array}\right)$ (48)
$\displaystyle=\frac{-i}{\sqrt{2}^{3}}\left(c_{i,\pi}^{T},c_{i,0}^{T}\right)\sigma_{y}\left(\begin{array}[]{cc}-1&0\\\
0&1\end{array}\right)\left(\begin{array}[]{c}c_{i,\pi}\\\
c_{i,0}\end{array}\right).$ (53)

In contrast, the symmetric leg-pairing operator is

$\displaystyle\sum_{l}\frac{1}{2}\Delta_{S}^{l}(i)$
$\displaystyle=\frac{-i}{\sqrt{2}^{3}}\left(c_{i,1}^{T},c_{i,2}^{T}\right)\sigma_{y}\left(\begin{array}[]{cc}1&0\\\
0&1\end{array}\right)\left(\begin{array}[]{c}c_{i,1}\\\
c_{i,2}\end{array}\right)$ (58)
$\displaystyle=\frac{-i}{\sqrt{2}^{3}}\left(c_{i,\pi}^{T},c_{i,0}^{T}\right)\sigma_{y}\left(\begin{array}[]{cc}1&0\\\
0&1\end{array}\right)\left(\begin{array}[]{c}c_{i,\pi}\\\
c_{i,0}\end{array}\right),$ (63)

while its antisymmetric counter part is

$\displaystyle\sum_{l}\frac{(-1)^{l}}{2}\Delta_{S}^{l}(i)$
$\displaystyle=\frac{i}{\sqrt{2}^{3}}\left(c_{i,1}^{T},c_{i,2}^{T}\right)\sigma_{y}\left(\begin{array}[]{cc}1&0\\\
0&-1\end{array}\right)\left(\begin{array}[]{c}c_{i+1,1}\\\
c_{i+1,2}\end{array}\right)$ (68)
$\displaystyle=\frac{i}{\sqrt{2}^{3}}\left(c_{i,\pi}^{T},c_{i,0}^{T}\right)\sigma_{y}\left(\begin{array}[]{cc}0&1\\\
1&0\end{array}\right)\left(\begin{array}[]{c}c_{i+1,\pi}\\\
c_{i+1,0}\end{array}\right).$ (73)

We next use the continuum limit $c_{i,a,\sigma}\simeq
R_{a,\sigma}(x)e^{ik_{a}x_{i}}+L_{a,\sigma}(x)e^{-ik_{a}x_{i}}$, where $R,L$
are continuum 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=\pi$) we find

$\displaystyle\Delta_{S}^{r}(i)$
$\displaystyle\simeq\frac{-i}{\sqrt{2}}\sum_{a}(-1)^{a}R_{a,\sigma}(x_{i})\sigma_{y}L_{a,\sigma}(x_{i})$
(74) $\displaystyle\Delta_{S}^{l}(i)$
$\displaystyle\simeq\frac{-i}{\sqrt{2}}\sum_{a}\cos(k_{a})R_{a,\sigma}(x_{i})\sigma_{y}L_{a,\sigma}(i)\,.$
(75)

##  Appendix D Further Numerical Results

###  D.1 Fourier transform of the pairing correlation functions

Figure 13: Real part of the Fourier transformation of the rung-rung, 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 $2k_{F}$.

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.

###  D.2 Compressibility, spin gap, and charge gap for $V=0$

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.

Figure 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).

###  D.3 Charge structure factor for different system sizes

Figure 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. 15 shows a comparison of the results for the symmetric part of the
density structure factor for system sizes $L=80,\,\ldots,\,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.

Figure 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.

###  D.4 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_{F}^{0}$ oscillations for $n=0.05$
as a function of $J/t$.

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

##  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

$\ket{\Psi}=\sum_{i_{1},i_{2}}\Phi(i_{1},i_{2})c_{i_{1}\uparrow}^{\dagger}c_{i_{2}\downarrow}^{\dagger}\ket{0}$
(76)

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

$\frac{rt^{2}}{J}=\frac{1}{I_{0}(E)}-E,$ (77)

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 $\vec{Q}=0$) as
follows,

$I_{0}(E)=\frac{1}{2\pi}\int_{-\pi}^{\pi}\frac{\text{d}\vec{q}}{E-2\epsilon_{\vec{q}}}\,.$
(78)

Although in our case we cannot assume the $90^{\circ}$ 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 $\pi$) and
halving of the respective matrix element $t$. This yields

$\displaystyle\begin{split}I_{0}(E)=\frac{1}{2\pi}\int_{-\pi}^{\pi}\text{d}q_{x}\;\Big{(}&\frac{1}{E+2t+4t\cos(q_{x})}\\\
+&\frac{1}{E-2t+4t\cos(q_{x})}\Big{)}\end{split}$ (79)

and consequently Eq. (37). Furthermore, we find that in the case of the ladder
$r=12$ in (77).


* (12) H. Goldman and E. Fradkin, _Loop models, modular invariance, and three-dimensional bosonization_ , _Physical Review B_ 97 (2018) 195112.
  * (13) R. Dijkgraaf and E. Witten, _Topological gauge theories and group cohomology_ , _Communications in Mathematical Physics_ 129 (1990) 393.
  * (14) E. Witten, _Fermion path integrals and topological phases_ , _Reviews of Modern Physics_ 88 (2016) 035001.
  * (15) D. Gaiotto and A. Kapustin, _Spin tqfts and fermionic phases of matter_ , _International Journal of Modern Physics A_ 31 (2016) 1645044\.
  * (16) L. Bhardwaj, D. Gaiotto and A. Kapustin, _State sum constructions of spin-tfts and string net constructions of fermionic phases of matter_ , _Journal of High Energy Physics_ 2017 (2017) 1.
  * (17) J. Fröhlich and P. Marchetti, _Soliton quantization in lattice field theories_ , _Communications in Mathematical Physics_ 112 (1987) 343.
  * (18) A. Cappelli, E. Randellini and J. Sisti, _Three-dimensional topological insulators and bosonization_ , _Journal of High Energy Physics_ 2017 (2017) 1.
  * (19) E. Fradkin and S. Kivelson, _Modular invariance, self-duality and the phase transition between quantum hall plateaus_ , _Nuclear Physics B_ 474 (1996) 543.
  * (20) E. Witten, _Three lectures on topological phases of matter_ , _La Rivista del Nuovo Cimento_ 39 (2016) 313.
  * (21) A.P.O. Chan, T. Kvorning, S. Ryu and E. Fradkin, _Effective hydrodynamic field theory and condensation picture of topological insulators_ , _Phys. Rev. B_ 93 (2016) 155122.
  * (22) H.-G. Zirnstein and B. Rosenow, _Cancellation of quantum anomalies and bosonization of three-dimensional time-reversal symmetric topological insulators_ , _Physical Review B_ 88 (2013) 085105.
  * (23) S.D. Geraedts and O.I. Motrunich, _Line of continuous phase transitions in a three-dimensional u (1) loop model with 1/r 2 current-current interactions_ , _Physical Review B_ 85 (2012) 144303.
  * (24) A. Polyakov, _Quark confinement and topology of gauge theories_ , _Nuclear Physics B_ 120 (1977) 429.
  * (25) W.-H. Hsiao and D.T. Son, _Duality and universal transport in mixed-dimension electrodynamics_ , _Physical Review B_ 96 (2017) 075127.
  * (26) W.-H. Hsiao and D.T. Son, _Self-dual $\nu$= 1 bosonic quantum hall state in mixed-dimensional qed_, _Physical Review B_ 100 (2019) 235150\.
  * (27) D. Dudal, A.J. Mizher and P. Pais, _Exact quantum scale invariance of three-dimensional reduced qed theories_ , _Physical Review D_ 99 (2019) 045017.
  * (28) A. Cappelli and G.R. Zemba, _Modular invariant partition functions in the quantum hall effect_ , _Nuclear Physics B_ 490 (1997) 595.
  * (29) E. Witten, _Non-abelian bosonization in two dimensions_ , _Communications in Mathematical Physics_ 92 (1984) 455.
  * (30) R. Floreanini and R. Jackiw, _Self-dual fields as charge-density solitons_ , _Physical Review Letters_ 59 (1987) 1873.
  * (31) A. Cappelli, G.V. Dunne, C.A. Trugenberger and G.R. Zemba, _Conformal symmetry and universal properties of quantum Hall states_ , _Nucl. Phys. B_ 398 (1993) 531 [hep-th/9211071].
  * (32) P. Ginsparg, _Applied conformal field theory_ , _arXiv preprint hep-th/9108028_ (1988) .
  * (33) A. Vishwanath and T. Senthil, _Physics of three-dimensional bosonic topological insulators: Surface-deconfined criticality and quantized magnetoelectric effect_ , _Physical Review X_ 3 (2013) 011016\.
  * (34) D. Gaiotto, A. Kapustin, N. Seiberg and B. Willett, _Generalized global symmetries_ , _Journal of High Energy Physics_ 2015 (2015) 1\.
  * (35) R. Thorngren, _Anomalies and bosonization_ , _Communications in Mathematical Physics_ 378 (2020) 1775.
  * (36) J. McGreevy, _Generalized symmetries in condensed matter_ , _Annual Review of Condensed Matter Physics_ 14 (2023) 57.
  * (37) P.R. Gomes, _An introduction to higher-form symmetries_ , _arXiv preprint arXiv:2303.01817_ (2023) .
  * (38) L. Bhardwaj, L.E. Bottini, L. Fraser-Taliente, L. Gladden, D.S. Gould, A. Platschorre et al., _Lectures on generalized symmetries_ , _Physics Reports_ 1051 (2024) 1.
  * (39) R. Luo, Q.-R. Wang and Y.-N. Wang, _Lecture notes on generalized symmetries and applications_ , _arXiv preprint arXiv:2307.09215_ (2023) .
  * (40) A. Kapustin and N. Seiberg, _Coupling a qft to a tqft and duality_ , _Journal of High Energy Physics_ 2014 (2014) 1.
  * (41) A. Karch, D. Tong and C. Turner, _A web of 2d dualities: $\mathbb{Z}_{2}$ gauge fields and arf invariants_, _SciPost Physics_ 7 (2019) 007.
  * (42) A. Kapustin and R. Thorngren, _Fermionic spt phases in higher dimensions and bosonization_ , _Journal of High Energy Physics_ 2017 (2017) 1.
  * (43) X. Chen, A. Tiwari and S. Ryu, _Bulk-boundary correspondence in (3+ 1)-dimensional topological phases_ , _Physical Review B_ 94 (2016) 045113.
  * (44) P. Ye and Z.-C. Gu, _Vortex-line condensation in three dimensions: a physical mechanism for bosonic topological insulators_ , _Physical Review X_ 5 (2015) 021029.
  * (45) B. Binegar, _Relativistic field theories in three dimensions_ , _Journal of Mathematical Physics_ 23 (1982) 1511.
  * (46) R. Jackiw and V.P. Nair, _Relativistic wave equation for anyons_ , _Physical Review D_ 43 (1991) 1933.
  * (47) T. Hansson, V. Oganesyan and S.L. Sondhi, _Superconductors are topologically ordered_ , _Annals of Physics_ 313 (2004) 497\.
  * (48) M. Metlitski and A. Vishwanath, _Particle-vortex duality of 2d dirac fermion from electric-magnetic duality of 3d topological insulators_ , in _APS March Meeting Abstracts_ , vol. 2016, pp. K51–006, 2016.
  * (49) A.G. Abanov and A. Cappelli, _Hydrodynamics, anomaly inflow and bosonic effective field theory_ , _arXiv e-prints_ (2024) arXiv:2403.12360 [2403.12360].