# Dimensional crossover on multileg attractive-$U$ Hubbard ladders

Anastasia Potapova  Ian Pilé  HSE University, 101000 Moscow, Russia  Tian-
Cheng Yi  Beijing Computational Science Research Center, Beijing 100084, China
Rubem Mondaini  Beijing Computational Science Research Center, Beijing 100084,
China  Evgeni Burovski  HSE University, 101000 Moscow, Russia

(August 29, 2024)

###### Abstract

We study the ground state properties of a polarized two-component Fermi gas on
multileg attractive-$U$ Hubbard ladders. Using exact diagonalization and
density matrix renormalization group method simulations, we construct grand
canonical phase diagrams for ladder widths of up to $W=5$ and varying
perpendicular geometries, characterizing the quasi-one-dimensional regime of
the dimensional crossover. We unveil a multicritical point marking the onset
of partial polarization in those phase diagrams, a candidate regime of finite-
momentum pairing. We compare our findings with recent experimental and
theoretical studies of quasi-one-dimensional polarized Fermi gases.

## Introduction

Quantum liquids in low dimensions are predicted to support a variety of exotic
superfluid phases, the so-called Fulde-Ferrell-Larkin-Ovchinnikov state (FFLO)
[1, 2] being a prime example. Here the superconducting order parameter is
modulated in space due to Cooper pairing at finite center-of-mass momentum.
The FFLO state remains elusive despite being predicted more than fifty years
ago. There is a growing body of indirect evidence
[*[andreferencestherein]Agterberg2020_cuprates_review] for FFLO in experiments
in $\mathrm{CeCoIn_{5}}$ in rotating magnetic field [4], organic
superconductors [5, 6], pnictides [7] and cuprates [8]. Yet, direct
experimental observation of the FFLO state in condensed matter systems remains
challenging. One reason the FFLO phase is elusive is that the impurity
scattering can completely suppress the FFLO long-range ordering [9].

Ultracold gases provide a promising platform for the quest for FFLO-like
physics thanks to the ability to engineer clean systems and unprecedented
control over polarization, interactions, and geometry. Several geometries have
been explored: spherically symmetric three-dimensional (3D) traps [10, 11,
12], elongated cigar-shaped traps [13, 14] and arrays of tubes [15]. Due to an
applied trapping potential, the polarized gas cloud phase separates into the
unpolarized (equal densities, EQ) and partially polarized (PP) regions—and the
PP gas is a candidate for the FFLO state. The phase boundaries between regions
can be directly measured via an in situ imaging. By varying either the
confinement ratio or the intertube tunneling, experiments probe the 1D-to-3D
crossover of polarized Fermi gases. An immediate conclusion from the
experiments is that the relative locations of the PP and EQ states in an
external potential are inverted between the 1D and 3D limits: in 3D, the EQ
phase occupies the center of the trap, while in 1D the center is PP and EQ is
pushed to the wings.

The limiting cases of the 1D-to-3D crossover of a polarized two-component
Fermi gas are relatively well understood from the theory side. In the 3D case,
the uniform polarized state with FFLO long-range ordering is only stable in a
narrow range of interactions and polarizations [16]. In two dimensions (2D),
the FFLO state is believed to be stable in a wide range of parameters [17],
including the regime of light doping from half-filling with small spin
polarization [18]. Mean-field studies [19], and Monte Carlo simulations
indicate further possibilities for competing long-range orders. For a review,
see, e.g. [20] and references therein.

In 1D, the whole PP phase has a dominant FFLO-like algebraic ordering [21] and
is stable in a wide range of interactions. For a trapped 1D gas, the phase
separation scenario is consistent with the experiments in highly elongated
cigar-shaped traps [22]. Following the initial field-theoretical [21] and
Bethe ansatz [22] treatments, the 1D polarized Fermi gas was extensively
studied numerically, using DMRG [23, 24, 25, 26, 27] and QMC simulations [28,
29, 30].

The 1D-to-3D crossover for polarized Fermi gases has been studied in the mean-
field approximation (MFA) in Refs. [31, 32, 33]. While the MFA approach allows
tracing the main features of the crossover qualitatively, it fails to capture
several qualitative effects seen in experiments at low densities close to the
1D limits: the MFA misrepresents the topology of the grand canonical phase
diagram in the 1D limit [31, 33].

This paper addresses strong correlation effects beyond MFA in the quasi-1D
regime. To this end, we study the quasi-1D attractive Hubbard model on wide
ladders of up to five legs. We employ two numerically unbiased approaches:
exact diagonalization (ED) of the Hamiltonian on small clusters of varying
aspect ratio and density-matrix renormalization group (DMRG) simulations [34,
35] in ladder-geometries. Our numerical results for two and three legs are
consistent with previous simulations of Refs. [36] and [37]. By considering
wider ladders, we map out the evolution of the $T=0$ grand canonical phase
diagram for quasi-1D geometries from a strict 1D limit to higher dimensions.
Specifically, we find that (i) in presence of the external trapping potential,
the 3D scenario of the phase separation is a robust feature that sets in
immediately away from the strict 1D limit, and (ii) the stability region of
the FFLO-like phase shrinks upon increasing the ladder’s width, i.e., away
from the strict 1D limit.

We also note that considering wide ladders (with flat or cylinder geometry in
the perpendicular direction) is becoming feasible from the computational point
of view and has been recently used for various problems related to variants of
the repulsive Hubbard model [38, 39, 40].

## Model and Methods

We consider the attractive-$U$ Hubbard model of a two-component Fermi gas at
zero temperature on a $W\times L$ ladder defined by the Hamiltonian:

$\hat{H}=-t\sum_{\langle
ij\rangle,\sigma}\left(\hat{c}_{i,\sigma}^{\dagger}\hat{c}_{j,\sigma}^{\phantom{\dagger}}+{\rm
h.c.}\right)+U\sum_{i}\hat{n}_{i,\uparrow}\hat{n}_{i,\downarrow}\ .$ (1)

Here $\hat{c}_{i,\sigma}^{\dagger}$ ($\hat{c}_{i,\sigma}$) is a creation
(annihilation) operator of a fermion with pseudospin
$\sigma=\uparrow,\downarrow$ on site $i$;
$\hat{n}_{i,\sigma}=\hat{c}_{i,\sigma}^{\dagger}\hat{c}_{i,\sigma}^{\phantom{\dagger}}$
is the corresponding number operator. We set the hopping amplitude $t=1$ as
the energy scale, and $U<0$ is the local Hubbard attraction parameter between
two fermions with opposite spins. The summation over $i$ in Eq. (1) runs over
$L\times W$ sites, where $L$ is the length of an $W$-leg ladder.

In the canonical ensemble, the ground state energies,
$E_{0}(N_{\uparrow},N_{\downarrow})$, of Eq. (1) are characterized by the
cardinality of spin-up and spin-down fermions, $N_{\uparrow}$ and
$N_{\downarrow}$, or equivalently by the filling fractions,
$n_{\sigma}=N_{\sigma}/LW$. Considering the case where $n_{\uparrow}\geqslant
n_{\downarrow}$, the following ground states are of interest:

  1. 1.

Vacuum (V): $n_{\uparrow}=n_{\downarrow}=0$

  2. 2.

Equal densities (EQ): $n_{\uparrow}=n_{\downarrow}$. At $T=0$ this phase
supports the BCS regime.

  3. 3.

Partially polarized (PP) phase: $n_{\uparrow}>n_{\downarrow}>0$. This phase is
the FFLO candidate.

  4. 4.

Fully polarized (FP) phase: $n_{\uparrow}>0$ and $n_{\downarrow}=0$.

For a lattice model (1), we further distinguish whether the majority spin band
is filled ($n_{\uparrow}=1$—we call this the FP2 phase) or not (i.e.,
$n_{\uparrow}<1$—we call this the FP1 phase).

Changing the variables from the canonical to the grand canonical ensemble, we
introduce the effective magnetic field $h$ and the chemical potential $\mu$
via

$\mu={\left(\frac{\partial E_{0}}{\partial N}\right)}_{P}\;,\qquad
h={\left(\frac{\partial E_{0}}{\partial P}\right)}_{N}\;,$ (2)

where $N=N_{\uparrow}+N_{\downarrow}$ is the total particle number and
$P=N_{\uparrow}-N_{\downarrow}$ is the polarization.

The phase boundaries are found by approximating the derivatives (2) with
finite differences — this is equivalent to comparing the ground state energies
of different phases. For instance, the boundary between V and EQ is given by
comparing the energy of a two-particle state with the zero energy of an empty
band, i.e. $\mu=E_{0}(1,1)/2$. Likewise, the boundary between FP1 and V is
related to the bottom of the single-particle band, $\mu=-h+E_{0}(1,0)$.
Further details, including the finite-difference approximations to Eq. (2) for
computing the EQ-PP and FP-PP phase boundaries, are given in Appendix A

## Results and discussion

We resort to DMRG computations 111For that, we make use of both ALPS [45] and
iTensor [46] implementations. to extract the ground-state energies
$E_{0}(N_{\uparrow},N_{\downarrow})$ of the Hamiltonian (1) on ladders of
length $L$ and widths $W$ of up to $W=5$ with open boundary conditions in
_both_ dimensions. For that, we fix the length $L=40$, and check that using
larger values of $L$ does not significantly change the results. We also run
simulations of $W=3,4,5$ with the cylinder geometry (i.e., periodic boundary
conditions in the perpendicular direction). Appendix B highlights the
necessity of using DMRG techniques by showing that ED results are largely
impacted by finite-size effects in the small clusters amenable to
computations. For the interaction strength, we take $U=-7$. At these
parameters, we typically use up to 2000 states for the bond dimension, so that
the truncation error in the DMRG process is below $10^{-8}$ for $W\leqslant 3$
and below $10^{-6}$ for $W=4,5$.

For the strictly 1D model ($W=1$), our numerical results are consistent with
previous Bethe Ansatz results [22], and DMRG simulations [27]. For $W=2$ and
$W=3$ our results agree with previous DMRG simulations [37, 36]. Figure 1
shows the grand canonical phase diagram for the $W=4$ ladder—other $W>2$
diagrams are qualitatively similar. We also checked that imposing periodic
boundary conditions in perpendicular directions does not significantly alter
the results. Several features stand out in Fig. 1. First, the PP phase—the
candidate to support the FFLO-like physics—occupies a significant area on the
phase diagram. This fact was previously seen in simulations of $W\leqslant 3$
ladders, and our results for $W=4,5$ support this being a robust feature of
quasi-1D geometries.

Figure 1: The phase diagram for $W=4$, $U=-7$, and $L=40$. Phases considered
are vacuum (V) ($n_{\uparrow}=n_{\downarrow}=0$), equal densities
(EQ)($n_{\uparrow}=n_{\downarrow}$), partially polarized
(PP)($n_{\uparrow}>n_{\downarrow}>0$) and two fully polarized phases: FP1 with
$n_{\uparrow}>0,n_{\downarrow}=0$ and FP2 with $n_{\uparrow}=0$ ,
$n_{\downarrow}>0$. Here we only show the region $h>0$ and $\mu<U/2$ because
the phase diagram is symmetric for $h\to-h$ and $\mu\to U/2-\mu$. Filled
circles are DMRG results, and dashed lines connecting them are guides to the
eye. The marker colors are mapped by the value $N_{\uparrow}$ that identifies
the boundary line according to the corresponding colorbars – See Eqs. (4) and
(5).

Second, the EQ, PP and FP1 phases meet at a “multicritical” point $O$, whose
location is fixed by the few-body physics of the model,
$(\mu_{O},h_{O})=(E_{0}(1,1)/2,-E_{0}(1,1)/2+E_{0}(1,0))$. That the EQ-FP
phase boundary also terminates at the point $O$ can be traced to the fact that
the only bound states in the Hubbard model (1) are pairs of spin-up and spin-
down fermions [42]. We checked that the binding energies of three-body states
(“trimers”) are zero for all values of $W$, modulo finite-$L$ corrections. We
also note that one of the main artifacts of the mean-field approximation is
the prediction that the point $O$ splits into two tricritical points [31, 33].
A compilation of parts of the phase diagram close to the multicritical point
$O$ with a growing number of legs is shown in Fig. 2.

The PP-FP1 (“polaron”) line.— The phase boundary between FP1 and PP phases
corresponds to the so-called Fermi-polaron problem: a single spin-down
particle in a sea of $N_{\uparrow}$ spin-up fermions (the points along the
line correspond to the increasing $N_{\uparrow}$). While for $W=1$ the phase
boundary is smooth and monotonic, for $W>1$ the boundary has kinks, which can
be traced to the filling of non-interacting energy bands in the spectrum [36].
In general, for a $W$-leg ladder, there are $W$ (overlapping) branches of the
single-particle spectrum, and kinks on the “polaron” line in the $\mu$-$h$
plane correspond to onsets of partial filling of multiple branches. Whenever
multiple branches are partially filled, the “polaron” line is non-monotonic
and has additional chaotic oscillations. These are likely due to the
finite-$L$ level spacing of the motion along the ladder.

Another point to note is that for $W\geqslant 2$, the initial part of the
“polaron” line—which corresponds to the filling of the lowest single-particle
band— is close to linear. This initial linear behavior is followed by a cusp
and a sharp downturn for the range where the second branch starts filling.
This behavior was first observed on a two-leg ladder in Ref. [36] and our
simulations indicate that it persists for larger values of $W$ and varying the
boundary conditions in the perpendicular direction.

_The EQ-PP line.—_ We now turn our attention to the boundary between the EQ
and PP phases. The main difference between the strictly 1D case of $W=1$ and
$W>1$ is that in the vicinity of the “multicritical” point $O$, the slope of
the boundary $\partial h/\partial\mu<0$ for $W=1$ [22, 27] while for $W>1$ the
slope of the boundary has the opposite sign, $\partial h/\partial\mu>0$. This
has been first noted for $W=2$ in Ref. [36], and later confirmed for $W=2,3$
in Ref. [37]. Our results, illustrated in Figs. 1 and 2 for up to $W=5$
corroborate the conclusion that the sign change is a robust feature that
differentiates the strict 1D limit from quasi-1D geometries with $W\neq 1$.

In the presence of an external trapping potential, the shell structure (or
phase separation) of the atomic cloud can be directly read off the grand
canonical phase diagram in the local density approximation. For a spin-
independent trapping potential, $V(x)$, applied along the ladder direction,
the cloud structure corresponds to $\mu\to\mu-V(x)$, i.e., a vertical cut on
the $h-\mu$ plane. This way, the sign of the slope of the EQ-PP boundary
differentiates between the 1D behavior with a two-shell structure where the PP
phase occupies the center of the trap [14, 15], and a higher-dimensional
behavior where the shell structure is inverted: the center of the trap has
equal densities and the PP phase occupies an outer shell [10, 11, 15]. Our
results, illustrated in Fig. 1 and 2, show that the 3D-like shell structure is
a robust feature of the quasi-1D geometries. We stress that in our
simulations, the EQ-PP boundary is monotonic, in contrast to the mean-field
results [31, 33], where the MF approximation generates artifacts at low
density (equivalently, strong interactions) which lead to an apparently
reentrant behavior 222In fact, we do observe some reentrant behavior if we
allow the hopping amplitudes along the rungs of the ladder, $t_{\perp}$ to
differ from the hopping amplitudes along the ladder, $t_{\parallel}$. We
reserve a thorough discussion of these phenomena for a separate
investigation..

Figure 2: The vicinity of the “multicritical” point $O$ for $W=1$ (black
line), for $W=2$ (blue line), $W=3$ (green line), $W=4$ (purple line) and
$W=5$ (yellow line) and $L=40$. Here we shift the boundaries for each $W$
value so that for each $W$’s “multicritical” point $O$ is at the origin. We
also only show for each $W$ the EQ-PP and PP-FP1 lines. Filled symbols are
DMRG results, and dashed lines are to guide an eye. We also note that the
vicinity of the multicritical point corresponds to the continuum limit of the
lattice model (1) since the filling fractions are small, $N_{\sigma}\ll WL$.

_The PP phase stability region.—_ In Fig. 2, we compare the low-filling parts
of phase diagrams for ladders with $W=1$ to $W=5$, where we shift the values
of $\mu$ and $h$ for each $W$ so that their respective points $O$ coincide. In
Fig. 2 we only show the EQ-PP and PP-FP1 lines for clarity. It is immediately
seen from Fig. 2 that with increasing the width $W$, the EQ-FP boundary shifts
to the right while the “polaron” line shifts upwards.

We thus conclude that the stability region for the PP phase—the candidate
phase for the FFLO-like physics—shrinks on approach away from a strict 1D
limit _while in the quasi-1D regime_. The last comment is due to the fact that
the position of the kink and a downturn on the “polaron” line, visible in Fig.
2, shifts to smaller values of filling fraction, $N_{\uparrow}/WL$, for
increasing $W$. Since the position of the kink is related to the filling of
the second single-particle branch of the transversal motion, and the gap
between branches goes to zero in 2D or 3D, our approach does not allow us to
make definite statements about the behavior of the stability region of the
FFLO phase for $W\sim L$.

_The FLLO in the PP phase.—_ To probe the FFLO character of the PP phase, we
compute the superconducting correlation functions,
$\Gamma(\mathbf{x},\mathbf{x}_{0})=\langle\hat{\Delta}^{\dagger}(\mathbf{x}_{0})\hat{\Delta}(\mathbf{x}+\mathbf{x}_{0})\rangle$,
where
$\hat{\Delta}(\mathbf{x})=\hat{c}_{\mathbf{x}\uparrow}\hat{c}_{\mathbf{x}\downarrow}$
annihilates a pair of fermions at a lattice site with coordinates
$\mathbf{x}$. In the 1D limit, $W=1$, $\Gamma$ is expected [21] to decay
algebraically with $|\mathbf{x}|$ and oscillate with the typical FFLO momentum
$Q=\pi(n_{\uparrow}-n_{\downarrow})$.

Figure 3: Fourier-transformed pair-pair correlation functions
$\widetilde{\Gamma}(k_{x},k_{y})$, in momentum space on a $5\times 40$ ladder;
insets show the corresponding real-space decay of such correlations. The total
density $n=n_{\uparrow}+n_{\downarrow}=0.1$ is fixed, whereas the polarization
changes. In (a) $N_{\uparrow}=N_{\downarrow}=10$; (b) $N_{\uparrow}=15$,
$N_{\downarrow}=5$; and (c) $N_{\uparrow}=19$, $N_{\downarrow}=1$. The
vertical dashed lines depict the characteristic momenta from FFLO in one-
dimensional systems, $Q=\pi(N_{\uparrow}-N_{\downarrow})/L$. See text for
discussion. The different curves in the main panels show the different allowed
transverse momentum modes $k_{y}$, whereas in the insets the actual
transversal coordinate $y$.

To trace the behavior of the superconducting correlations in ladder
geometries, we consider a $W=5$ ladder and take $\mathbf{x}_{0}$ at the center
of the lattice, $\mathbf{x}_{0}=(L/2,0)$, where the legs of the ladder are
indexed from $-(W-1)/2$ to $(W-1)/2$. Figure 3 shows the Fourier transform,
$\widetilde{\Gamma}(k_{x},k_{y})$, of $\Gamma(\mathbf{x},\mathbf{x}_{0})$ for
$L=40$, and different polarizations $(N_{\uparrow},N_{\downarrow})=(10,10)$
[panel (a)], $(15,5)$ [panel (b)], and $(19,1)$ [panel (c)], with a fixed
total occupancy. The insets display the corresponding real-space correlations
along the ladder, resolved by the $y$-coordinate in the transverse direction.

For balanced occupancies, the real-space pair-pair correlations decay
monotonically with distance which translates to a
$\widetilde{\Gamma}(k_{x},k_{y})$ peaked at zero-momentum [Fig. 3(a)]. In
turn, once a finite density-imbalance is selected, oscillations on the
correlations in real space are clearly visible, resulting in peaks of the
Fourier-transformed correlations at finite momentum [Figs. 3(b) and (c)]. This
is direct evidence of FFLO-type physics occurring in wider ladders. The peak
location in $k_{x}$, however, which occurs at the same point irrespective of
the $k_{y}$ value, is not aligned with the expected typical FFLO momentum of
1D systems. In practice, a combination of modes sets in, a fact already
appreciated in $W=2$ ladders in Ref. [36].

Figure 4: (a) Fourier transformed pair-pair correlations
$\widetilde{\Gamma}(k_{x},k_{y})$ with zero-transverse momentum $(k_{y}=0)$,
and fixed particle-number $N=20$. Data at $k_{x}=0$ and $P=0$ is omitted for
clarity. (b) The corresponding real-space correlations along the central leg
($y=0$) of the $5\times 40$ ladder. (c) The compilation of the peak-locations
in panel (a) vs. the polarization $P=N_{\uparrow}-N_{\downarrow}$ accompanied
by a linear fitting and the typical momentum dependence in 1D FFLO problems.

Since this finite-momentum peak depends on the total polarization, we compile
in Fig. 4 the Fourier transformed $\widetilde{\Gamma}(k_{x},0)$ correlations
with growing polarization $P=0\ldots 18$, and fixed particle number $N=20$ on
the same $5\times 40$ ladder. A finite-momentum peak at $k_{x}^{\rm max}$
occurs only when $P>0$ [Fig. 4(a)], tied to the appearance of oscillations
with the distance of the real-space correlations [Fig. 4(b)]. This peak
position linearly follows the polarization $P$, as with the expected typical
FFLO momentum $Q=\pi(n_{\uparrow}-n_{\downarrow})$, albeit with a slightly
decreased slope [Fig. 4(c)].

Whether such oscillations and finite-momentum pairing remain in truly two-
dimensional systems, studied utilizing unbiased methods, is an open question.
Our results point out that still in the quasi-one-dimensional limit, they yet
occur. Lastly, it is important to emphasize that mean-field investigations
unveil the possible scenarios of normal-polarized phases [31, 33], that is,
without pairing formation. In studies that employed unbiased methods as ours
in $W=2$ ladders [36], such regimes have been ruled out but deserve future
investigation at wider ladders.

## Summary

We study the quasi-1D limit of the dimensional crossover of a polarized Fermi
gas. By simulating the attractive-$U$ Hubbard model on wide ladders
encompassing up to $W=5$ legs, we trace the evolution with $W$ of the grand
canonical phase diagram, which features BCS-like and FFLO-like states. Our
simulations complement earlier mean-field-based investigations of the
dimensional crossover, where the uncontrollable nature of the mean-field
approximation leads to artifacts at strong coupling for quasi-1D systems. We
show that a qualitative difference in a shell structure of a gas in an
external potential, observed in experiments with ultracold gases, is a robust
feature that differentiates a strict 1D limit and higher dimensions, including
quasi-1D systems. Our simulations indicate that the stability region of the
FFLO phase in the quasi-1D regime shrinks with increasing the ladder width
away from the strict 1D case.

###### Acknowledgements.

We heartily thank Th. Jolicoeur and G. Misguich for multiple fruitful
discussions. We acknowledge financial support by RFBR and NSFC, project number
21-57-53014. R.M. acknowledges support from the NSFC Grants No. U2230402,
12111530010, 12222401, and No. 11974039. Numerical simulations were carried
out using the HSE HPC resources [44], and the Tianhe-2JK at the Beijing
Computational Science Research Center.

##  Appendix A Extracting phase boundaries in ED and DMRG data

In the main text, we argue that an analysis of the energetics of different
filling sectors of the Hamiltonian can extract the different phases on the
$\mu-h$ plane. Here we present this analysis in an expanded level of detail.
In the canonical ensemble, the ground state energies,
$E_{0}(N_{\uparrow},N_{\downarrow})$, of Eq. (1) are characterized by the
numbers of spin-up and spin-down fermions, $N_{\uparrow}$, and
$N_{\downarrow}$. Considering the case where $N_{\uparrow}\geqslant
N_{\downarrow}$, the following ground states are of interest:

  1. 1.

Vacuum (V): $N_{\uparrow}=N_{\downarrow}=0$

  2. 2.

Equal Densities (EQ): $N_{\uparrow}=N_{\downarrow}$. At $T=0$ this phase
supports the BCS ordering.

  3. 3.

Partially polarized (PP) phase: $N_{\uparrow}>N_{\downarrow}>0$. This phase is
the FFLO candidate.

  4. 4.

Fully polarized (FP) phase: $N_{\uparrow}>0$ and $N_{\downarrow}=0$.

For a lattice model, we further distinguish whether the majority spin band is
filled ($N_{\uparrow}=L\times W$—we call this the FP2 phase) or not (i.e.,
$N_{\uparrow}<L\times W$—we call this the FP1 phase). Changing the variables
from the canonical to the grand canonical ensemble, we introduce the effective
magnetic field $h$ and the chemical potential $\mu$ via

$\mu={\left(\frac{\partial E_{0}}{\partial N}\right)}_{P}\;,\qquad
h={\left(\frac{\partial E_{0}}{\partial P}\right)}_{N}\;,$ (3)

where $N=N_{\uparrow}+N_{\downarrow}$ is the total particle number and
$P=N_{\uparrow}-N_{\downarrow}$ is the polarization.

The phase boundaries are found by approximating the derivatives (3) with
finite differences, leading to

  1. 1.

The boundary between EQ and PP:

$\mu,h=\frac{E_{0}(N_{\uparrow}+1,N_{\downarrow}\pm
1)-E_{0}(N_{\uparrow},N_{\downarrow})}{2}$ (4)

where on the right-hand side the plus sign gives $\mu$ and the minus sign
gives $h$. Note that the phase boundary corresponds to
$N_{\uparrow}=N_{\downarrow}$ in Eq. (4).

  2. 2.

Likewise, the PP-FP1 boundary between is given by

$\mu,h=\frac{E_{0}(N_{\uparrow}\pm 1,1)-E_{0}(N_{\uparrow},0)}{\pm 2}$ (5)

  3. 3.

The boundary between EQ and V

$\mu=\frac{E_{0}(1,1)}{2}$ (6)

  4. 4.

The boundary between V and FP1

$\mu=-h+E_{0}(1,0)$ (7)

  5. 5.

The boundary between FP1 and FP2

$\mu=-h-E_{0}(1,0)$ (8)

From the ED or DMRG data, we construct the grand canonical phase diagrams
using Eqs. (4)-(8) (where the phase boundaries are parameterized by
$N_{\uparrow}$). For the strictly 1D model ($W=1$), our results are consistent
with previous Bethe Ansatz results [22], and DMRG simulations [27]. For $W=2$
and $W=3$ our results agree with previous DMRG simulations [37, 36].

##  Appendix B ED results and comparison to global minimization

Figure 5: Similar to the main text, the ground phase diagrams of the model (1)
with $U=-7$ on 16-sites cluster with aspect ratios $l_{x}\times l_{y}$: (a)
$16\times 1$ and (b) $4\times 4$. Phases are labeled similarly to Fig. 1.
Filled markers give the loci of the EQ-PP transition, according to Eq. (4);
empty markers the PP-FP1 phase transition, Eq. (5). The latter display
substantially large oscillations stemming from finite-size effects. Dotted
lines give the EQ-V transition; dashed, the V-FP1 transition, and dash-dotted
ones, the FP1-FP2 transition for the corresponding system sizes. Overlaid to
it are colors that map the ground state according to the different phases when
inspecting all $(N_{\uparrow},N_{\downarrow})$ sectors of the Hamiltonian.

We benchmark our exploration in the main text by using the ED method in small
clusters. In particular, we contrast one- and two-dimensional geometries with
the same number of sites (16-site clusters), observing how the dimensionality
affects the $\mu-h$ phase diagram. While a characterization via finite
differences is sufficient to reliably extract the phase boundaries as shown in
the main text and further described in Appendix A, the small number of sites
that can be tackled within this technique makes this procedure yield only a
rough location of the limits between different phases.

This is exemplified in Fig. 5 for the two aforementioned cluster shapes. An
overall similar structure of the phase diagram is seen for the DMRG results in
$W=4$ ladders [Fig. 1], but the large oscillations on the phase boundaries
given by the markers here attest to the large finite-size effects in such
lattice sizes. Because all different sectors $(N_{\uparrow},N_{\downarrow})$
are readily available in this method, we can also draw the different
boundaries by direct comparison of the total ground-state energy:

$\displaystyle E_{0}(\mu,h)=\min_{\\{N_{\uparrow},N_{\downarrow}\\}}[$
$\displaystyle E_{0}(N_{\uparrow},N_{\downarrow})$ (9)
$\displaystyle-\mu(N_{\uparrow}+N_{\downarrow})-h(N_{\uparrow}-N_{\downarrow})]\
,$

where direct inspection of the $E_{0}$’s can be used to classify which phase a
given point $(\mu,h)$ the ground-state belongs to. By mapping the
corresponding five different phases (see main text) to colors, we can then
describe a refined phase diagram in Fig. 5, which follows in tandem with the
description of the finite difference estimation used in the main text. While
also possible in the DMRG method, such exploration does not give a much better
boundary between the phases because the finite-size effects are smaller in the
larger cluster sizes we can tackle. Therefore, the finite difference scheme
suffices in drawing the phase diagram there.

Finally, we emphasize that the finite differences procedure described in
Appendix A is numerically friendlier than the one given by Eq. (9) since only
a limited set of ground-state energy computations is sufficient to obtain the
phase boundaries. That is, other than the trivial cases, $E_{0}(1,1)$ and
$E_{0}(1,0)$ (giving EQ-V, V-FP1 and FP1-FP2 boundaries), the remaining
boundaries require $E_{0}(N_{\uparrow}+1,N_{\uparrow}\pm 1)$ and
$E_{0}(N_{\uparrow},N_{\uparrow})$ (EQ-PP); and $E_{0}(N_{\uparrow}\pm 1,1)$,
$E_{0}(N_{\uparrow},0)$ for the PP-FP1 transition. Thus with a limited number
of runs, we can satisfactorily build the grand-canonical phase diagram of the
model as shown in Fig. 1, for example.

## References

  * Fulde and Ferrell [1964] P. Fulde and R. A. Ferrell, Superconductivity in a strong spin-exchange field, Phys. Rev. 135, A550 (1964).
  * Larkin and Ovchinnikov [1964] A. I. Larkin and Y. N. Ovchinnikov, Nonuniform state of superconductors, Sov. Phys. JETP 20, 762 (1964).
  * Agterberg  _et al._ [2020] D. F. Agterberg, J. S. Davis, S. D. Edkins, E. Fradkin, D. J. Van Harlingen, S. A. Kivelson, P. A. Lee, L. Radzihovsky, J. M. Tranquada, and Y. Wang, The physics of pair-density waves: Cuprate superconductors and beyond, Annual Review of Condensed Matter Physics 11, 231 (2020).
  * Lin  _et al._ [2020] S.-Z. Lin, D. Y. Kim, E. D. Bauer, F. Ronning, J. D. Thompson, and R. Movshovich, Interplay of the spin density wave and a possible Fulde-Ferrell-Larkin-Ovchinnikov state in ${\mathrm{cecoin}}_{5}$ in rotating magnetic field, Phys. Rev. Lett. 124, 217001 (2020).
  * Wosnitza [2018] J. Wosnitza, FFLO states in layered organic superconductors, Annalen der Physik 530, 1700282 (2018).
  * Imajo  _et al._ [2021] S. Imajo, T. Kobayashi, A. Kawamoto, K. Kindo, and Y. Nakazawa, Thermodynamic evidence for the formation of a Fulde-Ferrell-Larkin-Ovchinnikov phase in the organic superconductor $\lambda-{(\mathrm{BETS})}_{2}{\mathrm{gacl}}_{4}$, Phys. Rev. B 103, L220501 (2021).
  * Cho  _et al._ [2017] C.-w. Cho, J. H. Yang, N. F. Q. Yuan, J. Shen, T. Wolf, and R. Lortz, Thermodynamic evidence for the Fulde-Ferrell-Larkin-Ovchinnikov state in the ${\mathrm{kfe}}_{2}{\mathrm{as}}_{2}$ superconductor, Phys. Rev. Lett. 119, 217002 (2017).
  * Du  _et al._ [2020] Z. Du, H. Li, S. H. Joo, E. P. Donoway, J. Lee, J. C. S. Davis, G. Gu, P. D. Johnson, and K. Fujita, Imaging the energy gap modulations of the cuprate pair-density-wave state, Nature 580, 65 (2020).
  * Song and Koshelev [2019] K. W. Song and A. E. Koshelev, Quantum FFLO state in clean layered superconductors, Phys. Rev. X 9, 021025 (2019).
  * Shin  _et al._ [2006] Y. Shin, M. W. Zwierlein, C. H. Schunck, A. Schirotzek, and W. Ketterle, Observation of phase separation in a strongly interacting imbalanced Fermi gas, Phys. Rev. Lett. 97, 030401 (2006).
  * Zwierlein  _et al._ [2006a] M. W. Zwierlein, A. Schirotzek, C. H. Schunck, and W. Ketterle, Fermionic superfluidity with imbalanced spin populations, Science 311, 492 (2006a).
  * Zwierlein  _et al._ [2006b] M. W. Zwierlein, C. H. Schunck, A. Schirotzek, and W. Ketterle, Direct observation of the superfluid phase transition in ultracold fermi gases, Nature 442, 54 (2006b).
  * Partridge  _et al._ [2006] G. B. Partridge, W. Li, R. I. Kamar, Y. an Liao, and R. G. Hulet, Pairing and phase separation in a polarized Fermi gas, Science 311, 503 (2006).
  * Liao  _et al._ [2010] Y.-a. Liao, A. S. C. Rittner, T. Paprotta, W. Li, G. B. Partridge, R. G. Hulet, S. K. Baur, and E. J. Mueller, Spin-imbalance in a one-dimensional Fermi gas, Nature 467, 567 (2010).
  * Revelle  _et al._ [2016] M. C. Revelle, J. A. Fry, B. A. Olsen, and R. G. Hulet, 1d to 3d crossover of a spin-imbalanced Fermi gas, Phys. Rev. Lett. 117, 235301 (2016).
  * Sheehy and Radzihovsky [2006] D. E. Sheehy and L. Radzihovsky, BEC-BCS crossover in “magnetized” Feshbach-resonantly paired superfluids, Phys. Rev. Lett. 96, 060401 (2006).
  * Gukelberger  _et al._ [2016] J. Gukelberger, S. Lienert, E. Kozik, L. Pollet, and M. Troyer, Fulde-Ferrell-Larkin-Ovchinnikov pairing as leading instability on the square lattice, Phys. Rev. B 94, 075157 (2016).
  * Vitali  _et al._ [2022] E. Vitali, P. Rosenberg, and S. Zhang, Exotic superfluid phases in spin-polarized Fermi gases in optical lattices, Phys. Rev. Lett. 128, 203201 (2022).
  * Datta  _et al._ [2019] A. Datta, K. Yang, and A. Ghosal, Fate of a strongly correlated $d$-wave superconductor in a Zeeman field: The Fulde-Ferrel-Larkin-Ovchinnikov perspective, Phys. Rev. B 100, 035114 (2019).
  * Kinnunen  _et al._ [2018] J. J. Kinnunen, J. E. Baarsma, J.-P. Martikainen, and P. Törmä, The Fulde–Ferrell–Larkin–Ovchinnikov state for ultracold fermions in lattice and harmonic potentials: a review, journal = Reports on Progress in Physics,,  81, 046401 (2018).
  * Yang [2001] K. Yang, Inhomogeneous superconducting state in quasi-one-dimensional systems, Phys. Rev. B 63, 140511 (2001).
  * Orso [2007] G. Orso, Attractive Fermi gases with unequal spin populations in highly elongated traps, Phys. Rev. Lett. 98, 070402 (2007).
  * Feiguin and Heidrich-Meisner [2007] A. E. Feiguin and F. Heidrich-Meisner, Pairing states of a polarized Fermi gas trapped in a one-dimensional optical lattice, Phys. Rev. B 76, 220508 (2007).
  * Rizzi  _et al._ [2008] M. Rizzi, M. Polini, M. A. Cazalilla, M. R. Bakhtiari, M. P. Tosi, and R. Fazio, Fulde-Ferrell-Larkin-Ovchinnikov pairing in one-dimensional optical lattices, Phys. Rev. B 77, 245105 (2008).
  * Tezuka and Ueda [2008] M. Tezuka and M. Ueda, Density-matrix renormalization group study of trapped imbalanced Fermi condensates, Phys. Rev. Lett. 100, 110403 (2008).
  * Heidrich-Meisner  _et al._ [2010a] F. Heidrich-Meisner, A. E. Feiguin, U. Schollwöck, and W. Zwerger, BCS-BEC crossover and the disappearance of Fulde-Ferrell-Larkin-Ovchinnikov correlations in a spin-imbalanced one-dimensional Fermi gas, Phys. Rev. A 81, 023629 (2010a).
  * Heidrich-Meisner  _et al._ [2010b] F. Heidrich-Meisner, G. Orso, and A. E. Feiguin, Phase separation of trapped spin-imbalanced Fermi gases in one-dimensional optical lattices, Phys. Rev. A 81, 053602 (2010b).
  * Batrouni  _et al._ [2008] G. G. Batrouni, M. H. Huntley, V. G. Rousseau, and R. T. Scalettar, Exact numerical study of pair formation with imbalanced fermion populations, Phys. Rev. Lett. 100, 116405 (2008).
  * Casula  _et al._ [2008] M. Casula, D. M. Ceperley, and E. J. Mueller, Quantum Monte Carlo study of one-dimensional trapped fermions with attractive contact interactions, Phys. Rev. A 78, 033607 (2008).
  * Wolak  _et al._ [2010] M. J. Wolak, V. G. Rousseau, C. Miniatura, B. Grémaud, R. T. Scalettar, and G. G. Batrouni, Finite-temperature quantum Monte Carlo study of the one-dimensional polarized Fermi gas, Phys. Rev. A 82, 013614 (2010).
  * Parish  _et al._ [2007] M. M. Parish, S. K. Baur, E. J. Mueller, and D. A. Huse, Quasi-one-dimensional polarized Fermi superfluids, Phys. Rev. Lett. 99, 250403 (2007).
  * Sun and Bolech [2013] K. Sun and C. J. Bolech, Pair tunneling, phase separation, and dimensional crossover in imbalanced fermionic superfluids in a coupled array of tubes, Phys. Rev. A 87, 053622 (2013).
  * Sundar  _et al._ [2020] B. Sundar, J. A. Fry, M. C. Revelle, R. G. Hulet, and K. R. A. Hazzard, Spin-imbalanced ultracold Fermi gases in a two-dimensional array of tubes, Phys. Rev. A 102, 033311 (2020).
  * White [1992] S. R. White, Density matrix formulation for quantum renormalization groups, Phys. Rev. Lett. 69, 2863 (1992).
  * Schollwöck [2005] U. Schollwöck, The density-matrix renormalization group, Rev. Mod. Phys. 77, 259 (2005).
  * Feiguin and Heidrich-Meisner [2009] A. E. Feiguin and F. Heidrich-Meisner, Pair correlations of a spin-imbalanced Fermi gas on two-leg ladders, Phys. Rev. Lett. 102, 076403 (2009).
  * Burovski  _et al._ [2019] E. Burovski, R. Ikhsanov, A. Kuznetsov, and M. Kagan, Phase diagrams of polarized ultra-cold gases on attractive-$U$ Hubbard ladders, Journal of Physics: Conference Series 1163, 012046 (2019).
  * Kantian  _et al._ [2019] A. Kantian, M. Dolfi, M. Troyer, and T. Giamarchi, Understanding repulsively mediated superconductivity of correlated electrons via massively parallel density matrix renormalization group, Phys. Rev. B 100, 075138 (2019).
  * Chung  _et al._ [2020] C.-M. Chung, M. Qin, S. Zhang, U. Schollwöck, and S. R. White (The Simons Collaboration on the Many-Electron Problem), Plaquette versus ordinary $d$-wave pairing in the ${t}^{{}^{\prime}}$-Hubbard model on a width-4 cylinder, Phys. Rev. B 102, 041106 (2020).
  * Huang  _et al._ [2022] K. S. Huang, Z. Han, S. A. Kivelson, and H. Yao, Pair-density-wave in the strong coupling limit of the Holstein-Hubbard model, npj Quantum Materials 7, 17 (2022).
  * Note [1] For that, we make use of both ALPS [45] and iTensor [46] implementations.
  * Orso  _et al._ [2010] G. Orso, E. Burovski, and T. Jolicoeur, Luttinger liquid of trimers in Fermi gases with unequal masses, Phys. Rev. Lett. 104, 065301 (2010).
  * Note [2] In fact, we do observe some reentrant behavior if we allow the hopping amplitudes along the rungs of the ladder, $t_{\perp}$ to differ from the hopping amplitudes along the ladder, $t_{\parallel}$. We reserve a thorough discussion of these phenomena for a separate investigation.
  * Kostenetskiy  _et al._ [2021] P. S. Kostenetskiy, R. A. Chulkevich, and V. I. Kozyrev, HPC resources of the higher school of economics, Journal of Physics: Conference Series 1740, 012050 (2021).
  * Bauer  _et al._ [2011] B. Bauer, L. D. Carr, H. G. Evertz, A. Feiguin, J. Freire, S. Fuchs, L. Gamper, J. Gukelberger, E. Gull, S. Guertler, A. Hehn, R. Igarashi, S. V. Isakov, D. Koop, P. N. Ma, P. Mates, H. Matsuo, O. Parcollet, G. Pawłowski, J. D. Picon, L. Pollet, E. Santos, V. W. Scarola, U. Schollwöck, C. Silva, B. Surer, S. Todo, S. Trebst, M. Troyer, M. L. Wall, P. Werner, and S. Wessel, The ALPS project release 2.0: open source software for strongly correlated systems, Journal of Statistical Mechanics: Theory and Experiment 2011, P05001 (2011).
  * Fishman  _et al._ [2020] M. Fishman, S. R. White, and E. M. Stoudenmire, The ITensor software library for tensor network calculations (2020).