#
Attraction from frustration in ladder systems

Ivan Morera  Departament de Física Quàntica i Astrofísica, Facultat de Física,
Universitat de Barcelona, E-08028 Barcelona, Spain  Institut de Ciències del
Cosmos, Universitat de Barcelona, ICCUB, Martí i Franquès 1, E-08028
Barcelona, Spain  Department of Physics, Harvard University, Cambridge,
Massachusetts 02138, USA<EMAIL_ADDRESS>Annabelle Bohrdt  Department of
Physics, Harvard University, Cambridge, Massachusetts 02138, USA  Department
of Physics and Institute for Advanced Study, Technical University of Munich,
85748 Garching, Germany  Munich Center for Quantum Science and Technology
(MCQST), Schellingstr. 4, D-80799 München, Germany  ITAMP, Harvard-Smithsonian
Center for Astrophysics, Cambridge, MA 02138, USA  Wen Wei Ho  Department of
Physics, Harvard University, Cambridge, Massachusetts 02138, USA  MIT-Harvard
Center for Ultracold Atoms, Cambridge, MA, USA  Department of Physics,
Stanford University, Stanford, CA 94305, USA  Eugene Demler  Department of
Physics, Harvard University, Cambridge, Massachusetts 02138, USA  Institute
for Theoretical Physics, ETH Zurich, 8093 Zurich, Switzerland

###### Abstract

We analyze the formation of multi-particle bound states in ladders with
frustrated kinetic energy in two component bosonic and two component fermionic
systems. We focus on the regime of light doping relative to insulating states
at half-filling, spin polarization close to 100 percent, and strong repulsive
interactions. A special feature of these systems is that the binding energy
scales with single particle tunneling $t$ rather than exchange interactions,
since effective attraction arises from alleviating kinetic frustration. For
two component Fermi systems on a zigzag ladder we find a bound state between a
hole and a flipped spin (magnon) with a binding energy that can be as large as
0.6$t$. We demonstrate that magnon-hole attraction leads to formation of
clusters comprised of several holes and magnons and expound on
antiferromagentic correlations for the transverse spin components inside the
clusters. We identify several many-body states that result from self-
organization of multi-particle bound states, including a Luttinger liquid of
hole - magnon pairs and a density wave state of two hole - three magnon
composites. We establish a symmetry between the spectra of Bose and Fermi
systems and use it to establish the existence of antibound states in two
component Bose mixtures with SU(2) symmetric repulsion on a zigzag ladder. We
also consider Bose and Fermi systems on a square ladder with flux and
demonstrate that both systems support bound states. We discuss experimental
signatures of multi-particle bound states in both equilibrium and dynamical
experiments. We point out intriguing connections between these systems and the
quark bag model in QCD.

##  I Introduction

###  I.1 Motivation

Quantum frustrated systems hold the prospect of multifarious exotic many-body
states, including spin liquids and unconventional superconductors (see Ref.
Lacroix (2010) and references therein). In this paper we consider two
component mixtures of fermionic or bosonic particles on lattices with
frustrated kinetic energy in the regime of light doping relative to insulating
states at half-filling, spin polarization close to 100%, and repulsive SU(2)
symmetric interactions. We demonstrate that such systems exhibit a common
tendency to form multi-particle bound states with effective attraction between
particles arising from kinetic frustration. In many cases we find that the
ground state of such systems is best understood as resulting from self-
organization of multi-particle clusters. The systems under consideration can
be realized using ultracold Bose and Fermi atoms in optical lattices with the
currently available experimental techniques Bakr  _et al._ (2009); Sherson
_et al._ (2010); Parsons  _et al._ (2015); Cheuk  _et al._ (2015); Omran  _et
al._ (2015); Haller  _et al._ (2015); Edge  _et al._ (2015); Greif  _et al._
(2016); Brown  _et al._ (2017); Weitenberg  _et al._ (2011); Yang  _et al._
(2021); Sørensen  _et al._ (2005); Dalibard  _et al._ (2011); Atala  _et al._
(2014); Fukuhara  _et al._ (2013); Struck  _et al._ (2012); Cooper and
Dalibard (2013); Aidelsburger  _et al._ (2013); Miyake  _et al._ (2013); Jotzu
_et al._ (2014); Celi  _et al._ (2014); Goldman  _et al._ (2016); An  _et al._
(2017); Yang  _et al._ (2021). They may also be relevant for understanding
twisted bilayer graphene close to ferromagentic insulating states Sharpe  _et
al._ (2019).

We consider two specific lattice geometries: zigzag ladders and square ladders
with flux (see Fig. 1), but we expect that our conclusions apply to broader
classes of systems. Our starting “vacuum” is an insulating state with one
fermion per site on a frustrated ladder (zigzag or regular ladder with a flux)
and full spin polarization. When one of the spins is flipped, we call it a
magnon. When a single fermion is removed altogether, we call it a hole. The
starting point of our analysis is the observation of attraction between holes
and magnons arising from frustration of the kinetic energy of holes (see
discussion in Secs. II and III, and a related discussion for triangular
lattices in Zhang  _et al._ (2018)). We demonstrate that this attraction leads
to the formation of multi-particle bound states, which we denote as nHmM,
where integer $n$ and $m$ indicate the numbers of holes and magnons in the
composite object (see Fig. 1 for illustration), respectively.

Figure 1: Schematic representation of multi-particle bound states of a single
hole and one or several flipped spins (magnons) on top of a fully polarized
insulating state. (a) corresponds to the non-bipartite zigzag ladder and (b)
to the square ladder with a transverse magnetic flux $\Phi$.

The central objective of our paper is understanding the interplay of few- and
many-body phenomena. Many-body systems are usually understood from the
perspective of two-body correlations, such as the pairing amplitude in
superconductors that arises due to formation of Cooper pairs. Traditional
mean-field approaches rely on the order parameters defined to describe two
particle correlations (superconductivity, magnetism, spin and charge density
wave states). On the other hand, powerful theoretical methods have been
developed for analyzing few-body states in vacuum, such as the celebrated
Skorniakov-Ter-Martirosian (STM) equation Braaten and Hammer (2006); Greene
_et al._ (2017). However, the analysis of multi-particle composites within a
many-body system remains an open problem in many areas of physics. For
example, the primary objective of quantum chromodynamics (QCD) is developing
accurate models of triplets of quarks binding into nucleons, which in turn
combine to form nuclei. Intriguingly, we find many analogies between our
system and the “bag model” of QCD Hasenfratz and Kuti (1978). While in QCD
gluons provide a “bag” that holds quarks together, in the ladder systems
discussed in this paper, magnons provide a “bag” which traps one or several
holes, and in turn holes are holding the magnons together.

###  I.2 Overview of results

Figure 2: Phase diagram of the multi-particle bound states in the fermionic
zigzag ladder for different $J/t$ and different ratios of the magnon to hole
densities. The vertical dashed line represents the point $J/t\sim 2$ at which
the hole-magnon bound states disappears. The $\infty$ line represents the case
of a single hole with an arbitrary number of magnons.

Before presenting detailed microscopic calculations we summarize the main
results of our work in a non-technical manner.

Fermions in a zigzag ladder. Kinetic frustration of holes in a zigzag ladder
gives rise to formation of several types of multi-particle magnon-hole bound
states, such as shown in Figs. 1 and 2. The nature of stable composite objects
can be tuned by varying the ratio of interaction to single fermion hopping and
the relative concentration of holes and magnons, see Sec. III. We find large
binding energies of composites, since attraction is driven by the kinetic
energy of holes rather than the superexchange energy setting the amplitude of
spin dynamics. For example, the binding energy of 1H1M (one hole and one
magnon) can be as large as $0.6t$, where $t$ is the fermion hopping. We also
point out that the largest binding energy of the 1H1M pair is achieved in the
limit of vanishing superexchange interactions, which corresponds to large on-
site repulsion of two particles. Fig. 2 shows several cuts through the phase
diagram for simple ratios of the magnon and hole densities. We show that these
bound states have a common tendency to exhibit antiferromagnetic correlations
around the hole positions. In the case of a single hole, most of the magnons
accumulate around the hole forming an antiferromagnetic region and the rest
are pushed out of this region and form a spiral winding with a net $S^{z}$
magnetization as shown in Fig. 3. When multiple composite objects are present
in the system their statistics and interaction between them determine the
nature of the ground state. In Fig. 4 we report the pair correlation function
of single holes and hole-magnon pairs. We observe that one closely follows the
other one indicating that each hole in the system is paired with a magnon thus
creating a fluid of hole-magnon pairs. In Sec. IV.1 we show that this fluid
corresponds to a Luttinger liquid of fermionic 1H1M bound states, and in Sec.
IV.2 we discuss the pair density wave state of 2H3M composites.

Fermions and bosons in a square ladder with a flux.
Quantum particles on a square ladder in the presence of static gauge field
also experience frustration of kinetic energy as shown in Sec. VII. Optical
lattices with synthetic gauge fields have been realized experimentally for
both fermions and bosons Atala  _et al._ (2014); Aidelsburger  _et al._
(2013); Miyake  _et al._ (2013); Struck  _et al._ (2012); Jotzu  _et al._
(2014); Jiménez-García  _et al._ (2012), including square flux ladders Atala
_et al._ (2014); Fukuhara  _et al._ (2013). In Fig. 18 we present results for
the binding energy of 1H1M for different interaction strengths and different
flux per plaquette. We find the largest binding energy of around $0.3t$ in the
regime of vanishingly small superexchange energy.

Mott state of Bose atoms, antibound states.

Figure 3: Schematic picture of the antiferromagentic spin bag surrounding a
single hole. Inside the spinbag z-components of the spins are suppressed,
while the xy components exhibit antiferromagnetic correlations between two
sites of different legs. For clarity of presentation, in the figure we
oriented transverse components to point along the y axis. Outside the spinbag,
we find spiral winding of the xy components of the spins and z-magnetization
approaching its asymptotic value. The two angles $\theta$ and $\varphi$
characterize the magnetic orders present in the system. They correspond to the
angles between two adjacent spins in different legs in the xy plane and with
respect to the z axis, respectively. We define two regions: Region (I) inside
the spinbag and region (II) outside of it. In order to perform numerical
calculations we map the ladder to a one-dimensional system in which the sites
are indexed following the zigzag pattern specified by the numbers.

An important feature of hole-magnon systems for all lattice geometries is a
symmetry between the slightly doped Mott states of bosons and fermionic
models. The low energy states of the hole-magnon Hamiltonian in the fermionic
system can be mapped into the high energy states of the bosonic model, see
Sec. V for a detailed examination. This correspondence works in the regime of
strong repulsion between particles, when we can exclude real “doublons” and
describe both systems using $t$-$J$ type Hamiltonians. The change of sign of
the kinetic energy of holes in the two cases arises from the fact that a hole
creation operator is an annihilation operator of the original particle. Thus
the hole hopping term has the opposite sign relative to the hopping of the
original particle for fermions and the same sign for bosons (for details see
discussion in Sec. II). Superexchange interactions in the two systems also
differ in sign. Bosonic systems have ferromagentic exchange interaction Duan
_et al._ (2003) (we assume SU(2) symmetric repulsion between particles) while
fermions have antiferromagentic superexchange interactions Anderson (1950);
Hofstetter  _et al._ (2002). The consequence of the ${\cal H}$ to $-{\cal H}$
mapping is that bound hole-magnon states of fermionic systems correspond to
antibound states in the case of bosons. These antibound state can be
understood as analogues of repulsively bound states observed in experiments in
Winkler  _et al._ (2006) with one important differences that they are non-
local, i.e. the bound “elementary” particles in our systems (i.e. magnons and
holes) do not occupy the same site. It is worth noting that the symmetry
between bound and antibound states applies not only to the basic 1H1M pairs
but also to multi-particle clusters.

Dynamical probes of binding.

Figure 4: (a) Schematic representation of a self-organization of hole-magnon
pairs in the fermionic zigzag ladder. Hole-magnon pairs form a Luttinger
liquid with a characteristic Fermi momentum $k_{F}^{(B)}=\pi n_{B}$ given by
the density of pairs $n_{B}$. (b) Pair correlation function of single holes
and hole-magnon pairs in the fermionic zigzag ladder for $J/t=0$ and
$N_{h}=N_{m}=10$.

An efficient method of exploring bound and antibound states of many-body
systems is to use coherent quantum dynamics. To be concrete we consider the
case of a 1H1M pair, however, this discussion can be generalized to arbitrary
nHmM composites. The idea is as follows: use external local potential and/or
local optical spin control to prepare an initial state in which a hole and a
magnon are localized on neighboring sites. Switch off the localizing potential
and let the system evolve under the many-body Hamiltonian. When there are no
bound (antibound) states, the two particles will spread out essentially
independently and the probability to find them next to each other should
eventually decrease to zero (in a finite system it will saturate to a finite
value set by inverse of the system size). On the other hand, when there is a
bound (or antibound) state in the spectrum, we will find a finite probability
that the magnon and the hole spread out together. The probability to find the
two particles together after a long evolution time is given by the overlap
between the bound state wavefunction and the initial state. In Fig. 5 we show
dynamics in the zigzag ladder starting from the initial state in which one
hole and one magnon have been prepared on neighboring sites. We observe that
the joint probability distribution function remains peaked where the two
particles are close to each other. In Sec. VI we present a detailed analysis
of the quantum dynamics of a hole-magnon pair in the fermionic zigzag ladder.
The hole-magnon correlation at different times could be measured in current
ultracold atom experiments with a quantum gas microscope Bakr  _et al._
(2009); Sherson  _et al._ (2010); Parsons  _et al._ (2015); Cheuk  _et al._
(2015); Omran  _et al._ (2015); Haller  _et al._ (2015); Edge  _et al._
(2015); Greif  _et al._ (2016); Brown  _et al._ (2017); Weitenberg  _et al._
(2011); Yang  _et al._ (2021). A similar experimental technique has been used
previously in the Mott insulating regime of two component bosons to
demonstrate the existence of two magnon bound states in the Heisenberg model
in one-dimension Fukuhara  _et al._ (2013).

Figure 5: Dynamics of a hole-magnon pair in the zigzag ladder. Hole-magnon
correlation function $\langle\hat{n}^{h}_{i}\hat{n}^{m}_{j}\rangle$ computed
at different times for $J/t=0$. The initial state corresponds to putting a
hole and a flipped spin in adjacent sites in the middle of the lattice. At
long times a diagonal part can be discerned showing the correlated expansion
of the hole and the magnon bound together.

##  II Microscopic model

###  II.1 Kinetic frustration

While the most familiar mechanism of bound state formation is attractive
interaction between particles, kinetic frustration provides an alternative
route. Kinetic frustration arises in lattice systems when different quantum
paths of a propagating particle interfere destructively, see Fig. 6 panel (a)
and (b). In this situation the particle is not able to obtain the full kinetic
energy corresponding to the number of nearest neighbor sites in a non-
frustrated case. This frustration can be relieved by the presence of another
distinguishable particle, see Fig. 6 panel (c). When the original frustrated
particle moves from the initial to the final point along the first trajectory,
the second particle remains fixed at the same position. For the other path the
“frustrated” particle moves from the same initial to the same final point, but
its motion causes the second particle to move to a different site. For the two
trajectories the second particle ends up in different positions, hence, these
paths do not interfere and the “frustrated” particle can gain more kinetic
energy. Moreover, since alleviation of the kinetic energy frustration is
greater for shorter paths, it is energetically favorable for the two particles
to stay close to each other in space.

The primary goal of this paper is to present several examples of kinetic
frustration giving rise to formation of multi-particle bound states in two
component Bose and Fermi mixtures.

Figure 6: A single hole injected into a fully polarized insulating state
exhibits kinetic frustration arising from destructive interference of
different paths. This occurs for non-bipartite geometries such as the zigzag
ladder (a) or for bipartite ones in the presence of magnetic flux (b). When a
flipped spin is added into the system (c), the two different paths cannot
interfere since they lead to different final states. These two final states
correspond to the permutation of the two spins.

###  II.2 $t$-$J$ models

The $t-J$ model is given by,

$\hat{H}_{t-J}=\hat{\mathcal{P}}\left[-t\sum_{\langle
i,j\rangle,\sigma}\left(\hat{c}^{\dagger}_{i\sigma}\hat{c}_{j\sigma}+\text{h.c.}\right)\right]\hat{\mathcal{P}}+J\sum_{\langle
i,j\rangle}\mathbf{S}_{i}\mathbf{S}_{j},$ (1)

where $\hat{c}^{\dagger}_{i\sigma}$ creates a fermion or a boson at site $i$
with spin $\sigma=\uparrow,\downarrow$ and $\hat{\mathcal{P}}$ project onto
the subspace of no double occupancies. The spin operators are defined by
$S_{i}^{\gamma}=\sum_{\alpha,\beta}\hat{c}^{\dagger}_{i\alpha}\sigma_{\alpha,\beta}^{\gamma}\hat{c}_{i\beta}$
where $\sigma^{\gamma}$ are the Pauli matrices with $\gamma=x,y,z$. The sum
$\sum_{\langle i,j\rangle}$ is taken over first neighbors.

The $t-J$ model contains two energy scales: the quantum tunneling (hopping) of
particles $t$ and the superexchange $J$. The hopping $t$ gives rise to the
kinetic energy of holes, and the superexchange term describes propagation and
interaction of magnons. The $t-J$ model (1) is an effective low energy model
for the Hubbard model at very strong Hubbard interaction $U$ with $J=\pm
4t^{2}/U$, where $+$ ($-$) corresponds to the fermionic (bosonic) case Duan
_et al._ (2003); Trotzky  _et al._ (2008); Brown  _et al._ (2015); Nichols
_et al._ (2019); Jepsen  _et al._ (2020).

In order to study kinetic frustration in the $t-J$ model (1) we have to
specify the underlying lattice. Kinetic frustration can be easily achieved for
fermions in lattice geometries where the number of legs enclosed in a minimal
closed loop is odd. Therefore ladders are promising candidates where this
condition can be satisfied. The minimal lattice fulfilling that condition is a
zigzag ladder, in which a unit cell corresponds to two stacked triangles, see
Fig. 1 (a). Another promising platform for realizing a one-dimensional model
with kinetic frustration is flux ladder systems, which have been recently
realized in experiments with ultracold atoms Atala  _et al._ (2014); Tai  _et
al._ (2017). In this system a static synthetic gauge field is used to provide
the analogue of the Aharonov-Bohm flux of the ordinary magnetic field in
condensed matter systems. A particularly attractive feature of such systems is
the possibility of reaching high flux density per unit cell. In order to
numerically simulate such systems we have performed exact diagonalization
simulations for the cases of one hole and one magnon and density matrix
renormalization group (DMRG) simulations for systems with more holes and
magnons.

Single particle dispersion for the zigzag ladder.
The dispersion relation for a single hole in a fully polarized zigzag ladder
is given by $\epsilon_{h}(k)/t=\pm\\{\cos(k)+\cos(2k)\\}$, where $+$ ($-$)
corresponds to the fermionic (bosonic) case. For fermions kinetic frustration
manifests in the ground state of the hole since its ground state energy is
given by $\epsilon_{h}(k_{0})=-2.25t$ and $k_{0}=\arccos(-1/4)$. This ground
state energy is higher than the kinetic energy we would expect from counting
the number of nearest neighbors $-4t$. For bosons kinetic frustration appears
not in the ground state but for the highest energy single particle states. In
this situation the highest energy obtained $2.25t$ is lower than the expected
value of $4t$. To understand this effect, one can change the sign of hopping
$t\rightarrow-t$, so that the highest energy states of the positive $t$ model
correspond to the lowest energy states of the negative $t$ model. Then one
observes that in the negative $t$ case there is a relative minus sign between
the two trajectories around a single triangle (blue and red paths in Fig. 6),
which implies a destructive interference and thus a reduction of the kinetic
energy gain due to hopping. This will play a crucial role when studying bound
states of holes and magnons for the fermionic and bosonic systems.

Single particle dispersion for the square ladder with a flux
For a square flux ladder we consider propagation of a single hole on top of a
fully polarized insulating state. In this situation the dispersion relation of
the hole reads
$\epsilon_{\pm}({k})/t=-2\cos(\phi/2)\cos(k)\pm\sqrt{1+4\sin(\phi/2)^{2}\sin(k)^{2}}$
where $\phi$ corresponds to the magnetic flux per plaquette (see e.g. Tai  _et
al._ (2017)). The minimum of the band ($-$) is at $k_{m}=0$
($k_{m}=\pm\sin(\phi/2)\sqrt{1-\cos(\phi/2)/(4\sin(\phi/2)^{4})}$) for
$\phi<\phi_{c}$ ($\phi>\phi_{c}$) where $\phi_{c}\approx 0.43\times\pi$,
respectively. In this case two different bands are present $(\pm)$ because the
system can be thought of as two coupled 1D systems. It is easy to check that
the kinetic energy is always higher than the non-frustrated situation
$\epsilon_{\pm}({k})>-3t$ for all momenta when the magnetic flux is turned on.
Moreover, the spectrum is symmetric with respect to a change of sign of $t$,
shifting the momentum by $\pi$, and exchanging the two bands ($\pm$). This
symmetry shows that in this case frustration does not arise from the sign of
the bare $t$ (recall that hole hopping has different signs for the bosonic and
fermionic systems) but originates from the magnetic flux $\phi$ and appears
for both lowest and highest energy single particle states. As a consequence,
the dispersion of a single hole is the same for bosonic and fermionic systems.

##  III Bound States in fermionic Zigzag ladder

The kinetic energy of a hole gets modified when one of the possible paths
contains a magnon. This induces an effective interaction between the hole and
the magnon, which can be either repulsive or attractive depending on the
lattice geometry and particle statistics. Surprisingly this has implications
beyond the formation of bound hole-magnon pairs. In this section we
demonstrate that as one varies $J/t$ and relative concentrations of holes and
magnons, different types of composite objects are formed. We will focus on the
case of a fermionic zigzag ladder and discuss other systems in subsequent
sections.

###  III.1 Overview of bound states in fermionic zigzag ladder

We present the phase diagram of multi-particle bound states in the fermionic
zigzag ladder in Fig. 2. We note that we perform calculations at fixed ratios
of the hole to magnon densities, and with system sizes that are much larger
than the number of holes and magnons. Thus the phase diagram shown in Fig. 2
should be understood as corresponding to the limit of vanishingly small
concentrations of holes and magnons, with tuning parameters being $J/t$ and
the ratio between the densities of holes and magnons. Going to finite density
of dopants may result in stabilization of other phases, which we plan to
address in future work.

To understand the results presented in Fig. 2 it is useful to consider cuts
through the phase diagram at fixed $N_{m}/N_{h}$ as we vary $J/t$. We find
that generally large $J/t$ favors the decomposition of hole-magnon bound
states into free holes and magnons. By contrast, for small $J/t$ we find
stabilization of large composite objects. For the case $N_{m}/N_{h}=1$ we find
that the bound hole-magnon state appears only when $J/t<2$. The origin of the
suppression of hole-magnon binding for larger values of $J$ is discussed in
the next subsection. We do not find formation of larger composites, which
suggests that hole-magnon pairs repel each other for equal and small densities
of holes and magnons.

For $N_{m}/N_{h}=3/2$ we observe that at $J/t\sim 0.5$ a trimer bound states
1H2M is formed. Between $J/t\sim 0.3$ and $J/t\sim 0.5$ we find that hole-
magnon pairs and trimers coexist, which suggests repulsive interaction between
trimers and pairs. For smaller values $J/t<0.3$ a pentamer 2H3M composed of
two holes and three magnons is formed. This indicates that the attraction
between holes and magnons induces an effective attraction between two holes.
This is the smallest bound state that we have found which includes a pair of
holes. For $N_{m}/N_{h}=2$ we find a similar scenario but the hole-magnon
trimer appears for a wider range of $J/t$ and an hexamer 2H4M, which includes
a pair of holes, appears for a narrow range of $J/t$. In Sec. V we discuss the
many-body phases associated with the pentamer and the hole-magnon pair.

Finally we study a single hole immersed in a fluid of magnons: this is the
$N_{m}/N_{h}=\infty$ line in Fig. 2. By lowering the superexchange interaction
$J/t$ we observe that the number of magnons bound to the hole increases.
Magnons bound to the hole form a finite size cloud around the hole, a
phenomenon that we will refer to as formation of a “spinbag”. Within the
spinbag region we observe strong antiferromagnetic correlations of the XY spin
components between the two legs of the ladder as shown in Fig. 3. This feature
is a counterpart of the ferromagnetic spinbag (Nagaoka polaron) found in
bipartite lattices Nagaev (1992); Alexandrov (2007); Auerbach (1994); White
and Affleck (2001). Furthermore, similar to the Nagaoka polaron in the limit
$J/t\rightarrow 0$ the bag size becomes of the order of the system size, which
in our case indicates that the full system will exhibit antiferromagnetic
correlations of the transverse spin components. Antiferromagnetic correlations
surrounding the hole alleviate its kinetic frustration by making all possible
paths distinguishable, thus lowering the energy. In this way, hole-magnon
binding provides an explanation to the origin of antiferromagnetism in non-
bipartite lattices Haerter and Shastry (2005); Sposetti  _et al._ (2014).

###  III.2 Analysis of a 1H1M pair

Figure 7: Binding energy $E_{B}/t$ of the hole-magnon pair (red dots) on the
fermionic zigzag ladder obtained in the limit of having an infinite ladder by
performing a finite size scaling (see App. A) and the hole-magnon pair
extension (blue squares) $\sqrt{\langle z^{2}\rangle}/a$ for a ladder of
$2\times 200$ sites.

In order to solve the hole-magnon problem of the fermionic $t-J$ model Eq. (1)
in a zigzag ladder we first map the ladder to a one-dimensional system in
which the sites are indexed as shown in Fig. 3. Then we pass from the original
particle operators to operators of holes and magnons. In this way we reduce
the problem to a one-dimensional two-body problem, see Sup  for a detailed
examination. Then the effective two-body problem can be solved and all the
hole-magnon bound state properties can be obtained. We set the energy of the
fully polarized state with one fermion on every site to be zero, i.e. we
measure energies of all states relative to the fully polarized band insulating
state. The hole-magnon binding energy is defined as $E_{B}\equiv
E_{1H1M}-E_{1H}-E_{1M}$ where $E_{1H1M}$, $E_{1H}$ and $E_{1M}$ are the
groundstate energies for the hole-magnon, the single hole, and the single
manon states respectively. We find the largest binding energy when $J/t=0$,
see Fig. 7. This demonstrates that kinetic frustration provides an effective
attraction between a hole and a magnon. At the same time the superexchange
interaction $J$ provides an effective repulsion between them. To see this
effect one can consider the difference of superexchange energies of having a
hole and a magnon at adjacent sites relative to the case of being at larger
distances, $E_{J,\text{adjacent}}-E_{J,\text{separated}}=J/2$. This shows that
in the fermionic system with $J>0$ having a hole and a magnon on adjacent
sites increases the exchange energy. At $J/t=2$ the repulsion exactly cancels
the attraction and the hole-magnon pair vanishes.

Figure 8: Hole-magnon correlation function on the fermionic zigzag ladder for
the case of $N_{h}=N_{m}=1$ and different values of the superexchange
interaction $J/t=0.1,2.2$ for (a) and (b), respectively.

The hole-magnon bound state is characterized by having the two particles close
to each other in real space. By construction they are hardcore particles
because they cannot occupy the same site. By defining a relative distance
between the hole and the flipped spin $z$ we compute the size of the bound
state relative to the size of the system $\langle z^{2}\rangle/\langle
L^{2}\rangle$, see App. A. In the regime of strongest binding, $J/t\rightarrow
0$, the size of the bound state is much smaller than the system size and in
fact, independent of it. With increasing $J$ and decreasing binding energy,
the bound state expands (see Fig. 7). At $J/t=2$ the binding energy goes to
zero and the size of the bound state reaches the system size. Beyond that it
becomes proportional to the system size, indicating a transition from a bound
to a non-bound state (see App. A for a detailed finite size scaling analysis
of the size of the hole-magnon composite object). The hole-magnon correlation
function $\langle n^{h}_{i}n^{m}_{j}\rangle$ can be used to see the real space
structure of this bound state. For small values of $J/t$ we see that the
positions of the magnon and the hole are staying close to each other, see Fig.
8. Moreover the correlation function has a peak on adjacent sites which also
indicates a small spatial extent of the bound state. With increasing $J/t$ the
distance between the hole and the magnon begins to increase and positions of
the hole and the magnon become less correlated. These results indicate that
analysis of the hole-magnon correlation function should provide direct
experimental evidence of the bound state formation.

###  III.3 Analysis of the 1H2M trimer

Figure 9: Hole-magnon, (a) and (b), correlation function and magnon-magnon
correlation function, (c) and (d), on the fermionic zigzag ladder for
$N_{h}=1$, $N_{m}=2$ and two values of the superexchange interaction
$J/t=0.2,0.6$.

In the previous subsection we discussed that superexchange interaction gives
rise to an effective repulsion between a hole and a magnon, which competes
with attraction arising from kinetic frustration. The relative strength of
these opposing interactions determines the transition between trionic and pair
binding on the $N_{m}/N_{h}=2$ line. For small values of $J/t$ the ground
state of the system has 1H2M trimers. As $J/t$ increases, the system undergoes
a transition in which trimers dissociate into 1H1M pairs and free magnons.
Numerically we identify this transition by computing
$E_{B}=E_{1H2M}-E_{1H1M}-E_{1M}$. When the binding energy $E_{B}$ is negative,
a trimer state is stable. When $E_{B}$ is positive, trions become unstable to
dissociating into hole-magnon pairs and free magnons. We observe that this
transition occurs around $J/t\sim 0.5$. Correlation functions can also be used
to detect the appearance of a trimer. Since these trimers contain two magnons,
these bound states can be revealed by observing two magnon bunching in the
correlation function $\langle n^{m}_{i}n^{m}_{j}\rangle$. On the other hand,
for large values of $J/t$, when the system decomposes in a pair and a free
magnon, the magnon-magnon correlation function shows that the two magnons are
separated by a large distance comparable to the system size, see Fig. 9.
Moreover the hole-magnon correlation function always shows that there is a
magnon close to a hole (see the enhanced probability close to the diagonal
$i_{h}=j_{m}$ line in Fig. 9 upper row). This can be understood from the
observation that for both trimers and 1H1M pairs, we should find a hole and a
magnon close to each other.

###  III.4 2H2M: Effective repulsive interactions

Figure 10: (a) Hole-magnon, (b) hole-hole and (c) magnon-magnon correlation
functions on the fermionic zigzag ladder for $N_{h}=2$, $N_{m}=2$ and
$J/t=0.1$.

When multiple hole-magnon pairs are present in the system effective
interactions between them appear. In particular, magnons can mediate effective
interactions between holes. We analyze whether these effective interactions
are enough to bind a pair of holes in zigzag ladders with nearly full spin
polarization. Surprisingly we do not find hole pairs for the cases of $2H1M$
and $2H2M$. In the case of two holes and two magnons we observe that the
system decomposes into two hole-magnon pairs. The binding energy
$E_{2H2M}-2E_{1H1M}$ approaches zero when we increase the system size for any
value of $J/t$. This indicates that there is an effective repulsive force
between hole-magnon pairs. By computing the hole-hole and magnon-magnon
correlation we observe that the two holes and the two magnons are separated in
space, see Fig. 10. On the other hand, the hole-magnon correlation still shows
that a hole and a magnon sit close together. This indicates that the system
forms two pairs which then stay away from each other.

###  III.5 2H3M: Hole pairing

Figure 11: Hole density (a) for the case of $N_{h}=2$, $N_{m}=3$ and different
values of the superexchange interaction on the fermionic zigzag ladder. Hole
density of the projected state $|\psi_{h}\rangle$ (b) which has a fixed hole
in the middle of the lattice, see main text. We removed the central point to
improve visibility.

In the zigzag ladder the smallest bound state containing a pair of holes is a
pentamer formed by two holes and three magnons. We checked that the binding
energy of the pentamer is negative with respect to all possible decompositions
into smaller objects when $J/t<0.3$. We observe however that the binding
energy with respect to the decomposition into a trimer and a pair is small.
Therefore, although the pentamer is the groundstate of the system, already
small temperatures will lead to the breaking of a pentamer into a 1H1M pair
and a 1H2M trimer. In the remainder of this subsection we will focus on the
properties of the ground state of the system. In order to detect the
appearance of hole pairs we measure the hole density, see Fig. 11. For large
values of $J/t$ the hole density exhibits two humps indicating that the holes
are separated by a large distance. With decreasing $J/t$ the two humps begin
to approach to each other and finally fuse into a single one. After the fusion
the two holes share the same region of space forming the pentamer together
with three magnons. In order to quantify how far apart the two holes are, we
fix the position of the first hole and compute the spatial distribution of the
probability density of the second hole. Technically this is done by first
projecting the groundstate wavefunction into the state which contains a hole
in the center of the lattice
$|\psi_{h}\rangle\equiv\hat{P}^{h}_{L}|\psi\rangle=\hat{c}_{L\uparrow}\hat{c}^{\dagger}_{L\uparrow}\hat{c}_{L\downarrow}\hat{c}^{\dagger}_{L\downarrow}|\psi\rangle$
and then computing the density of the second hole
$\langle\psi_{h}|\hat{P}^{h}_{i}|\psi_{h}\rangle$ White and Scalapino (1997).
We also normalize the state $\langle\psi_{h}|\psi_{h}\rangle=1$. For large
values of $J/t$ this probability has maxima at the positions $L/4$ and $3L/4$,
which is the expected result for two hardcore particles. When reducing $J/t$
the positions of the two maxima start to approach the center of the lattice.
This shows that the two holes start approaching each other, which indicates
effective attraction between them and formation of a bound state. As a
function of $J/t$ the distance between the two maxima exhibits a minimum value
of $2\times 10$ sites at $J/t\sim 0.05$, see Fig. 11. We conclude that the
pentamer exists but it is a loosely bound object with a relatively large size.

###  III.6 1HNM: Antiferromagnetic spinbag

Figure 12: Magnonic affinity energy $\mu_{mN}/t$ on the fermionic zigzag
ladder (a) as a function of the superexchange interaction for the cases with
$N_{h}=1$ and different numbers of magnons. Number of magnons forming the
antiferromagnetic spinbag $N_{m}^{*}$ (b) as a function of the superexchange
interaction $J/t$. Dashed line represents the analytical result for the
antiferromagnetic spinbag $N_{m}^{*}\propto 1.10(t/J)^{1/3}$, see App. B.
Dotted line is the fit to numerical data which gives $N_{m}^{*}\propto
1.63(t/J)^{1/3}$.

We now discuss the character of a magnetic polaron when a single hole makes a
bound state with several magnons. With this goal in mind we studied systems
that consist of a single hole and up to eight magnons. We define the binding
energy $E_{B}^{(N)}\equiv E_{1HNM}-E_{1H}-NE_{1M}$ and then compute the
difference $\mu_{mN}\equiv E_{B}^{(N)}-E_{B}^{(N-1)}$. This quantity plays the
role of a magnonic affinity energy. We present its dependence on the
superexchange interaction for different number of magnons in Fig. 12 panel
(a). The negative value of this quantity indicates that when the $N$-th magnon
is added to the system, it is energetically favorable to attach it to the
1H(N-1)M bound state as an extra bound particle. Addition of magnons will
continue to increase the size of the cluster bound to the hole as long as the
magnon affinity remains negative. We refer to such clusters of bound magnons
as a “spinbag” (see Fig. 3). When the magnon affinity becomes positive, the
excess magnons are pushed away from the bag. The point at which the magnonic
affinity energy changes sign indicates the optimal number of magnons forming a
spinbag $N_{m}^{*}$. For small values of $J/t$ the number of magnons inside
the spinbag increases which also results in the expansion of its size, see
Fig. 12 panel (b). A simple scaling analysis can be used to understand the
dependence of $N_{m}^{*}$ on the ratio $J/t$ (see Sup  for a detailed
analysis). For small values of $J/t$ the antiferromagnetic spinbag is frozen
and it acts as an effective potential trapping the hole inside it. Therefore
the kinetic energy of the hole decreases for larger spinbags $t/N_{m}^{2}$. On
the other hand, creating the antiferromagnetic spinbag costs an extra magnetic
energy $JN_{m}$. This competition between kinetic and magnetic energy gives
the optimal number of magnons present in the spinbag
$N_{m}^{*}\propto(t/J)^{1/3}$. Similar scalings have been found for the
magnetic polaron in a two-dimensional square lattice Shraiman and Siggia
(1988); White and Affleck (2001); Grusdt  _et al._ (2018a)

###  III.7 Phase diagram

As we discussed before the fermionic $t-J$ model in a zigzag ladder exhibits a
large zoo of bound states with different numbers of holes and magnons. We now
summarize how we obtained transitions between different types of bound states
as a function of $J/t$. Our first step was to compute binding energies defined
previously for different numbers of the holes and magnons. We use this
procedure to identify the lowest energy bound states for different values of
$J/t$ and different ratios of the hole to magnon densities. Then we double the
number of particles and compute the respective binding energy $E_{2N}-2E_{N}$.
A negative binding energy denotes an effective attraction and particles will
cluster and create a larger bound state. Instead a positive binding energy
indicates an effective repulsive interaction and two separated bound states
appear. This means that we have found the largest possible bound state for
this ratio of number of magnons and holes at this value of $J/t$. Iterating
this procedure for different values of $N_{m}/N_{h}$ and $J/t$ we obtain the
phase diagram presented in Fig. 2.

##  IV Many-body phases of fermionic zigzag ladders

Figure 13: Binding energies relative to the hole-magnon pair 1H1M
$E_{B}=E_{10H10M}-10E_{1H1M}$ (blue dots) and the pentamer 2H3M
$E_{B}=E_{6H9M}-3E_{2H3M}$ (orange squares) for the fermionic zigzag ladder as
a function of the lattice length $2L$.

In this section we review many-body states that emerge in systems with many
bound states. More formally we analyze zigzag ladders as we increase the
number of holes and magnons while keeping the ratio $N_{m}/N_{h}$ fixed. After
identifying “optimal” multi-particle bound states, we do not find further
clustering. This suggest that optimal composites of holes and magnons exhibit
repulsive interactions with each other. We begin by discussing the situation
of small densities, when the interaction between different bound states is
much smaller than the binding energy. In this regime the system can be
understood as a dilute fluid of weakly interacting bound states, and the total
energy of the system $E_{NB}$ should be equal to $N_{B}E_{1B}$, where $E_{1B}$
is the energy of a single bound state and $N_{B}$ is the number of bound
states (see Fig. 13). We present results of these calculations for the cases
of $10H10M$ and $6H9M$ where the energy of the system converges to
$10E_{1H1M}$ and $3E_{2H3M}$, respectively. In the first case the energy
converges slowly with the system size. In the second case the energy converges
to the asymptotic value already for small sizes. This indicates strong
repulsive interactions between pentamers which forces them to avoid spatial
overlap.

###  IV.1 Luttinger liquid of composite fermions

Figure 14: Correlation functions of holes (dots) and hole-magnon pairs (upper
triangles) on the fermionic zigzag ladder for $N_{h}=10$, $N_{m}=10$ and
$J/t=0$. We note that operators $c_{\sigma}$ denote the original fermionic
operators.

Our next step is to establish the best physical picture of the quantum fluids
of bound states. As our first example we analyze a system in which we increase
the densities of holes and magnons but fix the ratio $N_{m}/N_{h}=1$. 1H1M
pairs exhibit strong binding with the scale of the binding energy set by $t$
and effective repulsive interaction between pairs. This suggests that the
many-body state can be understood as a Luttinger liquid of 1H1M pairs. In
order to show the paired nature of the liquid we define the pair correlation
function of hole-magnon pairs $g_{2}=\frac{\langle
n_{i}^{c}n_{j}^{c}\rangle}{\langle n_{i}^{c}\rangle\langle n_{j}^{c}\rangle}$
for $i>j$ with $n_{i}^{c}=n_{i}^{h}n_{i+1}^{m}+n_{i}^{m}n_{i+1}^{h}$ being
$n_{i}^{h}$ and $n_{i}^{m}$ the hole and magnon density on site $i$,
respectively. This correlator shows that every hole in the system is
accompanied by a magnon on an adjacent site, see Fig. 4. This result strongly
supports the hypothesis of the liquid of pairs, but we need to provide
additional verification that the long range part of correlations is consistent
with the Luttinger liquid phase.

The Luttinger liquid of bare holes and the Luttinger liquid of hole-magnon
composites can be rigorously distinguished by identifying the operators that
decay most slowly in space. To this end we compute two types of correlation
functions: the correlation function of bare holes $\langle\hat{c}_{\uparrow
i}\hat{c}^{\dagger}_{\uparrow j}\rangle$ and the correlation function of pairs
$\langle\hat{c}_{\uparrow i}S^{+}_{i}\hat{c}^{\dagger}_{\uparrow
j}S^{-}_{j}\rangle$ (see Fig. 14, we also note that here operators
$c_{\sigma}$ denote the original fermionic operators). The hole correlation
function has an exponential decay indicating that we cannot characterize the
system as a Luttinger liquid of individual holes. On the other hand, the
correlation function for pairs exhibits a much slower decay, which for large
distances approaches a power-law decay with oscillations with a characteristic
length scale $1/(\pi n_{B})$ set by the pair density
$n_{B}=N_{h}/(2L)=N_{m}/(2L)$ as predicted by Luttinger liquid theory
Giamarchi (2004). This indicates that a paired Luttinger liquid appears when
the number of magnons is equal to the number of holes and $J/t<2$.

###  IV.2 Pair density wave

Figure 15: Hole density on the fermionic zigzag ladder (a) for the case of
$N_{h}=6$ and $N_{m}=9$ at $J/t=0.05$. Hole density of the projected
wavefunction $|\psi_{3h}\rangle$ (b) which has three holes fixed at the
positions specified by the dashed lines.

In the earlier discussion of energetics we pointed out evidence of strong
repulsive interactions between pentamers. The repulsive interactions between
pentamers provides the opportunity to explore a many-body phase of pairs of
holes. We will now discuss numerical evidence that pentamers form a crystal
like phase. We call this state a pair density wave phase because every
pentamer is a charge two object. The hole density presented in Fig. 15
provides a strong indication of the pair density wave. In a system with six
holes we find three strong peaks in the density with very low density between
the peaks. This suggests that 2H3M pentamers are strongly localized. Moreover
by fixing the position of three holes and computing the hole density for the
remaining three ones we observe that the latter are strongly localized close
to the positions of the three fixed holes, see Fig. 15 panel (b). In the
middle region between two bound pairs of holes the probability of finding
another hole is negligible. Furthermore these densities seem very similar to
the one found for a single pentamer in Fig. 11 panel (b). Before concluding
this discussion we point out that to rigorously define a crystal phase in 1d
in the thermodynamic limit requires spontaneous breaking of discrete
translational symmetry. Breaking of translational symmetry is only possible
for rational densities of pentamers. For irrational densities one can at most
find quasi-long range order. Presence or absence of the long range crystal
order is expected to depend strongly on the precise pentamer density. To
identify true long range order in numerics calls for analyzing large system
sizes. We postpone this investigation to subsequent publications.

##  V Bosons in zigzag ladder

Figure 16: Band structure of the hole-magnon problem on the zigzag ladder as a
function of the total quasimomentum $Q$ for $J/t=0.1$. The grey area denotes
the hole-magnon scattering continuum and the continuous line the bound state
band. The left (right) panel corresponds to the fermionic (bosonic) case
$t>0,J>0$ ($t<0,J<0$). The band structure is obtained by solving the effective
two-body problem of a single hole and a magnon in the zigzag ladder, see App.
A.

The fermionic $t-J$ model in a zigzag ladder presents the unique opportunity
to explore the attraction between holes and magnons and thus attractive bound
states coming from kinetic frustration. On the other hand, for the bosonic
case repulsive hole-magnon bound states (antibound states) driven by kinetic
frustration can be found. These states are also characterized by having the
hole and magnon close in real space but their energy is higher than the
scattering continuum, see Fig. 16. Thus these are not the groundstate of the
system. As we will show there is a direct transformation relating the hole-
magnon bound state in both models.

In order to obtain the relation between the fermionic and bosonic model it is
convenient to express the Hamiltonian (1) in terms of the hole and magnon
operators rather than the original particles. Since we are working close to
the fully polarized state we can introduce the transformations
$t\,\hat{c}^{\dagger}_{\uparrow i}\hat{c}_{\uparrow j}\rightarrow\pm
t\,\hat{h}^{\dagger}_{j}\hat{h}_{i}$ and $t\,\hat{c}^{\dagger}_{\downarrow
i}\hat{c}_{\downarrow j}\rightarrow\pm
t\,\hat{h}^{\dagger}_{j}\hat{h}_{i}S^{+}_{i}S^{-}_{j}$. The positive and
negative signs correspond to the bosonic and fermionic case, respectively.
These signs come from the commutation and anticommutation relations satisfied
by bosons and fermions respectively. Moreover the superexchange coupling also
changes sign when going from bosons to fermions $J=\pm 4t^{2}/U$ being
positive for fermions and negative for bosons Duan  _et al._ (2003); Anderson
(1950); Hofstetter  _et al._ (2002). This establishes a connection between the
fermionic and bosonic $t-J$ models
$\hat{H}_{t-J}^{F}\Longleftrightarrow-\hat{H}_{t-J}^{B}$. If a hole-magnon
bound state appears below the scattering continuum for the fermionic case it
also appears above the continuum for the bosonic case, see Fig. 16. This sets
a connection between attractive (fermionic) and repulsive (bosonic) bound
states for $t-J$ models.

##  VI Experimental probes with dynamics

Figure 17: Main panel: Probability of finding a hole and a flipped spin at a
relative distance less than three lattice sites $P(|z|/a<3)$ as a function of
time for a fermionic zigzag lattice of $2\times 50$ sites for two values of
$J/t$. Inset panel: Probability of finding the hole and magnon at a relative
distance $z=a$ as a function of $J/t$ obtained in the limit of having an
infinite ladder by performing a finite size scaling, see App. A.

So far we discussed hole-magnon and other composites at zero temperature. We
expect that at temperature smaller than the binding energies our conclusions
will remain accurate. However, when temperature exceeds the binding energy, we
expect composites to dissociate. Thus one of the practical challenges for
experimental observation of the bound states that we discussed in previous
sections is the requirement of cooling systems to temperatures below the
binding energies. This difficulty is exacerbated by the requirement of having
strongly polarized Fermi mixtures for which Pauli principle cuts off contact
interactions and suppresses thermalization. An alternative is to perform
dynamical experiments. These are based on starting a real time evolution of a
fully polarized state which has a hole and a magnon localized on adjacent
sites in the middle of the lattice. This state has an overlap with the hole-
magnon bound state (around $0.3$). In the inset of Fig. 17 we present the
probability of finding the hole and the magnon at nearest sites as a function
of the superexchange coupling $J/t$. A finite overlap of the initial state
with the bound state ensures that during the coherent evolution there is a
finite probability for the hole and the magnon to stay close to each other as
shown in Fig. 17. In order to avoid finite size effects we have performed
dynamical simulations up to the point where the hole or magnon touches the
ends of the system. These finite time and finite size simulations make the
probability to always saturate to a finite value even if no bound state is
present. To dynamically discern the appearance of a hole-magnon bound state we
require to have saturation values much larger than one over the system size.
For small values of $J/t$ we observe saturation values that can be ten times
larger than one over the system size. Hence such “quantum walk” like dynamics
Preiss  _et al._ (2015) starting from the adjacent hole-magnon configuration
can reveal the existence of the bound state. The dynamical observation of
hole-magnon bound states could be adressed in current ultracold atom
experiments with a quantum gas microscope Bakr  _et al._ (2009); Sherson  _et
al._ (2010); Parsons  _et al._ (2015); Cheuk  _et al._ (2015); Omran  _et al._
(2015); Haller  _et al._ (2015); Edge  _et al._ (2015); Greif  _et al._
(2016); Brown  _et al._ (2017); Weitenberg  _et al._ (2011); Yang  _et al._
(2021).

It is interesting to point out the similarity of the fermionic and bosonic
systems under this time evolution. Since the two Hamiltonians are related by a
minus sign $\hat{H}_{t-J}^{F}\Longleftrightarrow-\hat{H}_{t-J}^{B}$ but the
procedure satisfies dynamical symmetry analogous to the one discussed in Ref.
Schneider  _et al._ (2012) we expect the time evolution to be identical in
both cases. Therefore far-from-equilibrium experiments can be used either to
detect attractive (fermionc) or repulsive (bosonic) bound states.

##  VII Fermions and bosons in flux ladders

Figure 18: Binding energy of the hole-magnon for a square ladder as a function
of the transverse magnetic flux for different values of the superexchange
interaction. The fermionic case corresponds to (a) and the bosonic one to (b).

A single hole propagating in a ferromagnetic background in a non-bipartite
lattice always experiences kinetic frustration. However in a square lattice
the kinetic energy of a single hole can also become frustrated if different
paths contribute with different relative phases. Therefore we explore the
appearance of hole-magnon bound states in a square ladder with magnetic flux.

Hole-magnon bound states are found and they exhibit very similar properties to
the ones obtained for the zigzag ladder. Near the hole the spins tilt the XY
plane and antiferromagnetic correlations are found, see Fig. 1. In order to
quantify the parameter space where hole-magnon bound states can be observed we
compute the respective binding energy $E_{B}\equiv E_{1H1M}-E_{1H}-E_{1M}$,
see Fig. 18. The bound state appears for a wide range of values of the
magnetic flux but the range of $J/t$ is small compared with the zigzag ladder.

The single hole spectrum for the regular ladder with flux is symmetric under a
change of sign of the single particle hopping. Thus starting with a spin
polarized band insulator for fermions or a Mott insulator of bosons with
$n=1$, we expect to find identical single particle spectra of individual
holes. This proves that attractive bound states are present for both bosonic
and fermionic systems. This is the main difference between bipartite and non-
bipartite lattices. In Fig. 18. we present the binding energies for the
fermionic panel (a) ($J/t>0$) and the bosonic panel (b) ($J/t<0$) situation.
Since we still have the property that fermionic and bosonic models written in
terms of hole and magnon operators differ by the sign, we infer that in both
systems attractive and repulsive bound states are present. The wavefunction
that describes an attractive/repulsive bound state for fermions also describes
a repulsive/attractive bound state for bosons.

##  VIII Summary and Outlook

The main result of our work is a demonstration of the effective attractive
interaction arising from kinetic frustration of distinguishable particles.
These interactions lead to formation of several types of multi-particle
complexes, which can then result in several types of interesting many-body
phases. In this paper we focused on two component Bose and Fermi mixtures and
two types of frustrated lattices: zigzag ladders and square ladders with a
flux. We studied the regimes close to fully polarized insulating states, when
these systems can be naturally described in terms of doped holes and magnons,
i.e. flipped spins. We made concrete predictions for the phase diagrams of
these models including identifying the most stable multi-particle bound states
and understanding many-body states that emerge from the interaction of such
composites. We expect that qualitatively, the results presented in our paper
should be valid to general models with kinetic frustration.

We used analysis of binding energies and correlation functions to discern the
lowest energy multi-paritcle bound states for different values of the
interaction strength and several representative density ratios. Examination of
the correlation functions also allowed us to scrutinize the internal structure
of these multi-particle bound states. We show that the system has a tendency
to create an antiferromagnetic background around the position of the hole (see
discussion in App. B). This antiferromagnetic background acts as an effective
potential to the hole which confines it in this region. This object can be
understood as a multi-particle bound state of many magnons and a single hole
and we term it an antiferromagnetic spinbag. When adding a second hole to the
system an effective attraction between the two holes can be found. We show
that the minimal configuration exhibiting a pair of two holes is a pentamer
composed of two holes and three magnons.

We identified a generic dynamic symmetry between the fermionic and bosonic
Hubbard models close to spin polarized insulating states. We demonstrated that
expressing these models in terms of the hole and magnon operators gives rise
to effective Hamiltonians which differ only in the overall sign. This symmetry
maps the low energy states of the fermionic model to the high energy states of
the bosonic model. Therefore when a bound state is present in the fermionic
system we should find an antibound state for the bosonic system and vice
versa: a bound state for bosonic system implies an antibound state for the
fermionic system.

We discussed dynamical experiments which can be used to probe formation of
multiparticle bound states. In particular, we considered a protocol for a
fermionic system in a zigzag ladder in which a hole and a flipped spin are
initialized on neighboring sites. After the hole and the magnon are released,
we find that their coherent dynamics can be decomposed into two contributions:
in the first one the hole and the magnon expand essentially independently of
each other, and in the second one they move together. The latter part can be
understood by observing that the initial configuration has a large overlap
with the 1H1M bound state, which results in the hole-magnon pair expanding as
a whole. Because of the dynamic symmetry between the fermionic and bosonic
models, this procedure can be applied either for detecting a bound state for
fermions or an antibound state for bosons. Moreover this probe could be
generalized to systems with more holes or magnons in order to detect larger
multi-particle bound or antibound states.

Our work extends beyond analysis of individual multi-particle bound states. We
discuss several many-body phases that can be understood from the perspective
of self organization of the multi-particle bound states. We argue that higher
order correlation functions provide direct signatures of the emergent many-
body states of composite objects and review two concrete examples. Our first
example is a Luttinger liquid of 1 hole - 1 magnon pairs. We show that in this
case composite operators for the hole-magnon pairs exhibit much slower spatial
decay than operators of the original particles. The second example we discuss
is a crystal of pentamers, in which every pentamer is comprised of two holes
and three magnons. We refer to this phase as a pair density wave state since a
pentamer is a charge two object.

Our work on multi-particle bound state formation due to kinetic frustration
can be extended in several directions. Powerful experimental tools developed
for cold atomic ensembles in optical lattices allow to satisfy two criteria at
the same time: relevance to current technologies and richness of theoretically
expected phenomena. A promising direction is to analyze ladders in which
tunneling/interactions within the chains and between them are different.
Applying potential gradients also allows to realize mixed dimensional systems,
in which along certain directions there is exchange interaction but no single
particle tunneling Grusdt  _et al._ (2018b). In the case of bosonic systems
one can also consider Hubbard models in which interactions are not SU(2)
symmetric. This should translate into the anisotropy of exchange interactions
Jepsen  _et al._ (2020, 2021). We discussed in Sec. III that the Ising part of
the interaction contributes to the repulsive interaction between a magnon and
a hole. Thus by making the z-axis part of the exchange interaction stronger,
we expect to find suppression of magnon-hole binding, which should in turn
affect stability of all multi-particle composites. A more detailed analysis of
multi-particle formation and self-organization in systems with magnetic flux
is also an interesting future directions. With increasing magnetic flux we
expect to find stronger frustration and stronger binding. This should make
multi-particle systems, for which we find binding energies to be small, easier
to observe in experiments.

Our results indicate that hole pairing is strongly suppressed for ladders
relative to the 2D case. A possible explanation is that the hole pairs in two
dimensions form at finite angular momentum and in ladders they are suppressed
because they cannot fully expand in the perpendicular direction. From the
results of Ref. Zhang  _et al._ (2018) one expects that 2D systems exhibit
hole pairing arising from formation of 2 holes - 1 magnon trimers. An
interesting question then is to extend our work to a larger number of legs and
study the transition from the quasi-1D geometry to the full 2D one.

The antiferromagnetic spinbags that we find provide an intriguing analogy to
the quark bag model used to explain the asymptotic freedom of QCD. In our
system the antiferromagnetic background provides an effective static potential
to the holes. While holes are fast degrees of freedom, they become trapped
inside these antiferromagnetic regions. A possible future perspective would be
to exploit this analogy and study the backreaction of the bag to the hole
movement. We expect that this will allow to find excitations of composite
particles in our system that are analogous to excitations of hadrons
considered in the context of the bag model Hasenfratz and Kuti (1978). In
particular we anticipate vibrational modes, that can be visualized as a hole
“rattling” inside the AF bag, as well as surface modes, that correspond to
local displacements of the bag with respect to its equilibrium position around
the hole. The latter family of excitations should be particularly interesting
in the 2D setting. Finally, this analogy could be used to explore other
directions with the goal of providing new insight into interesting open
questions regarding the nature of confinement. One particularly intriguing
question is to use an experimental platform of frustrated quantum systems to
realize an analogue of the transition between the hadronic and paired phases,
expected in QCD at higher density. Many questions remain poorly understood
regarding this phase transition, including the possibility of some
intermediate exotic phases. In our frustrated system a similar situation could
be explored in the “pentamer phase” when the density of holes to magnons is
fixed at 2/3 but individual densities are being increased. As the density of
pentamers increases, they should start to overlap, and at some critical
density holes may become “liberated” from the bags. This will corresponds to
breaking of pentamers, however, we expect that residual interaction between
holes will still result in their pairing. Understanding this transition,
including exploring the possibility of other phases, is an interesting open
question.

###### Acknowledgements.

We acknowledge helpful discussions with P. Dolgirev, M. Greiner, F. Grusdt,
Eun-Ah Kim and P. Schauss. I. Morera thanks financial support from the MINECO
(Spain) Grants No. FIS2017- 87534-P and FIS2017-84114-C2-1-P and the European
Union Regional Development Fund within the ERDF Operational Program of
Catalunya (project QUASICAT/ QuantumCat). A. Bohrdt is supported by the NSF
through a grant for the Institute for Theoretical Atomic, Molecular, and
Optical Physics at Harvard University and the Smithsonian Astrophysical
Observatory. W. W. Ho is supported in part by the Stanford Institute of
Theoretical Physics (SITP). E. Demler thanks Harvard-MIT CUA, the ARO grant
number W911NF-20-1-0163, the NSF grant OAC-1934714, the NSF EAGER-QAC-QCH
award number 2037687. DMRG computations have been performed using TeNPy
Hauschild and Pollmann (2018).

## References

  * Lacroix (2010) C. Lacroix,  _Introduction to frustrated magnetism : materials, experiments, theory_ (Springer, Berlin London, 2010).
  * Bakr  _et al._ (2009) W. S. Bakr, J. I. Gillen, A. Peng, S. Fölling, and M. Greiner, Nature 462, 74 (2009).
  * Sherson  _et al._ (2010) J. F. Sherson, C. Weitenberg, M. Endres, M. Cheneau, I. Bloch, and S. Kuhr, Nature 467, 68 (2010).
  * Parsons  _et al._ (2015) M. F. Parsons, F. Huber, A. Mazurenko, C. S. Chiu, W. Setiawan, K. Wooley-Brown, S. Blatt, and M. Greiner, Phys. Rev. Lett. 114, 213002 (2015).
  * Cheuk  _et al._ (2015) L. W. Cheuk, M. A. Nichols, M. Okan, T. Gersdorf, V. V. Ramasesh, W. S. Bakr, T. Lompe, and M. W. Zwierlein, Phys. Rev. Lett. 114, 193001 (2015).
  * Omran  _et al._ (2015) A. Omran, M. Boll, T. A. Hilker, K. Kleinlein, G. Salomon, I. Bloch, and C. Gross, Phys. Rev. Lett. 115, 263001 (2015).
  * Haller  _et al._ (2015) E. Haller, J. Hudson, A. Kelly, D. A. Cotta, B. Peaudecerf, G. D. Bruce, and S. Kuhr, Nature Physics 11, 738 (2015).
  * Edge  _et al._ (2015) G. J. A. Edge, R. Anderson, D. Jervis, D. C. McKay, R. Day, S. Trotzky, and J. H. Thywissen, Phys. Rev. A 92, 063406 (2015).
  * Greif  _et al._ (2016) D. Greif, M. F. Parsons, A. Mazurenko, C. S. Chiu, S. Blatt, F. Huber, G. Ji, and M. Greiner, Science 351, 953 (2016).
  * Brown  _et al._ (2017) P. T. Brown, D. Mitra, E. Guardado-Sanchez, P. Schauß, S. S. Kondov, E. Khatami, T. Paiva, N. Trivedi, D. A. Huse, and W. S. Bakr, Science 357, 1385 (2017).
  * Weitenberg  _et al._ (2011) C. Weitenberg, M. Endres, J. F. Sherson, M. Cheneau, P. Schauß, T. Fukuhara, I. Bloch, and S. Kuhr, Nature 471, 319 (2011).
  * Yang  _et al._ (2021) J. Yang, L. Liu, J. Mongkolkiattichai, and P. Schauss, (2021), arXiv:2102.11862 .
  * Sørensen  _et al._ (2005) A. S. Sørensen, E. Demler, and M. D. Lukin, Phys. Rev. Lett. 94, 086803 (2005).
  * Dalibard  _et al._ (2011) J. Dalibard, F. Gerbier, G. Juzeliūnas, and P. Öhberg, Rev. Mod. Phys. 83, 1523 (2011).
  * Atala  _et al._ (2014) M. Atala, M. Aidelsburger, M. Lohse, J. T. Barreiro, B. Paredes, and I. Bloch, Nature Physics 10, 588 (2014).
  * Fukuhara  _et al._ (2013) T. Fukuhara, P. Schauß, M. Endres, S. Hild, M. Cheneau, I. Bloch, and C. Gross, Nature 502, 76 (2013).
  * Struck  _et al._ (2012) J. Struck, C. Ölschläger, M. Weinberg, P. Hauke, J. Simonet, A. Eckardt, M. Lewenstein, K. Sengstock, and P. Windpassinger, Phys. Rev. Lett. 108, 225304 (2012).
  * Cooper and Dalibard (2013) N. R. Cooper and J. Dalibard, Phys. Rev. Lett. 110, 185301 (2013).
  * Aidelsburger  _et al._ (2013) M. Aidelsburger, M. Atala, M. Lohse, J. T. Barreiro, B. Paredes, and I. Bloch, Phys. Rev. Lett. 111, 185301 (2013).
  * Miyake  _et al._ (2013) H. Miyake, G. A. Siviloglou, C. J. Kennedy, W. C. Burton, and W. Ketterle, Phys. Rev. Lett. 111, 185302 (2013).
  * Jotzu  _et al._ (2014) G. Jotzu, M. Messer, R. Desbuquois, M. Lebrat, T. Uehlinger, D. Greif, and T. Esslinger, Nature 515, 237 (2014).
  * Celi  _et al._ (2014) A. Celi, P. Massignan, J. Ruseckas, N. Goldman, I. B. Spielman, G. Juzeliūnas, and M. Lewenstein, Phys. Rev. Lett. 112, 043001 (2014).
  * Goldman  _et al._ (2016) N. Goldman, J. C. Budich, and P. Zoller, Nature Physics 12, 639 (2016).
  * An  _et al._ (2017) F. A. An, E. J. Meier, and B. Gadway, Science Advances 3 (2017), 10.1126/sciadv.1602685.
  * Sharpe  _et al._ (2019) A. L. Sharpe, E. J. Fox, A. W. Barnard, J. Finney, K. Watanabe, T. Taniguchi, M. A. Kastner, and D. Goldhaber-Gordon, Science 365, 605 (2019).
  * Zhang  _et al._ (2018) S.-S. Zhang, W. Zhu, and C. D. Batista, Phys. Rev. B 97, 140507 (2018).
  * Braaten and Hammer (2006) E. Braaten and H.-W. Hammer, Physics Reports 428, 259 (2006).
  * Greene  _et al._ (2017) C. H. Greene, P. Giannakeas, and J. Pérez-Ríos, Rev. Mod. Phys. 89, 035006 (2017).
  * Hasenfratz and Kuti (1978) P. Hasenfratz and J. Kuti, Physics Reports 40, 75 (1978).
  * Jiménez-García  _et al._ (2012) K. Jiménez-García, L. J. LeBlanc, R. A. Williams, M. C. Beeler, A. R. Perry, and I. B. Spielman, Phys. Rev. Lett. 108, 225303 (2012).
  * Duan  _et al._ (2003) L.-M. Duan, E. Demler, and M. D. Lukin, Phys. Rev. Lett. 91, 090402 (2003).
  * Anderson (1950) P. W. Anderson, Phys. Rev. 79, 350 (1950).
  * Hofstetter  _et al._ (2002) W. Hofstetter, J. I. Cirac, P. Zoller, E. Demler, and M. D. Lukin, Phys. Rev. Lett. 89, 220407 (2002).
  * Winkler  _et al._ (2006) K. Winkler, G. Thalhammer, F. Lang, R. Grimm, J. Hecker Denschlag, A. J. Daley, A. Kantian, H. P. Büchler, and P. Zoller, Nature 441, 853 (2006).
  * Trotzky  _et al._ (2008) S. Trotzky, P. Cheinet, S. Fölling, M. Feld, U. Schnorrberger, A. M. Rey, A. Polkovnikov, E. A. Demler, M. D. Lukin, and I. Bloch, Science 319, 295 (2008).
  * Brown  _et al._ (2015) R. C. Brown, R. Wyllie, S. B. Koller, E. A. Goldschmidt, M. Foss-Feig, and J. V. Porto, Science 348, 540 (2015).
  * Nichols  _et al._ (2019) M. A. Nichols, L. W. Cheuk, M. Okan, T. R. Hartke, E. Mendez, T. Senthil, E. Khatami, H. Zhang, and M. W. Zwierlein, Science 363, 383 (2019).
  * Jepsen  _et al._ (2020) P. N. Jepsen, J. Amato-Grill, I. Dimitrova, W. W. Ho, E. Demler, and W. Ketterle, Nature 588, 403 (2020).
  * Tai  _et al._ (2017) M. E. Tai, A. Lukin, M. Rispoli, R. Schittko, T. Menke, D. Borgnia, P. M. Preiss, F. Grusdt, A. M. Kaufman, and M. Greiner, Nature 546, 519 (2017).
  * Nagaev (1992) E. Nagaev, Journal of Magnetism and Magnetic Materials 110, 39 (1992).
  * Alexandrov (2007) A. S. Alexandrov,  _Polarons in advanced materials_ (Canopus Pub, Dordrecht, 2007).
  * Auerbach (1994) A. Auerbach,  _Interacting Electrons and Quantum Magnetism_ (Springer New York, New York, NY, 1994).
  * White and Affleck (2001) S. R. White and I. Affleck, Phys. Rev. B 64, 024411 (2001).
  * Haerter and Shastry (2005) J. O. Haerter and B. S. Shastry, Phys. Rev. Lett. 95, 087202 (2005).
  * Sposetti  _et al._ (2014) C. N. Sposetti, B. Bravo, A. E. Trumper, C. J. Gazza, and L. O. Manuel, Phys. Rev. Lett. 112, 187204 (2014).
  * (46) See Supplementary Material.
  * White and Scalapino (1997) S. R. White and D. J. Scalapino, Phys. Rev. B 55, 6504 (1997).
  * Shraiman and Siggia (1988) B. I. Shraiman and E. D. Siggia, Phys. Rev. Lett. 60, 740 (1988).
  * Grusdt  _et al._ (2018a) F. Grusdt, M. Kánasz-Nagy, A. Bohrdt, C. S. Chiu, G. Ji, M. Greiner, D. Greif, and E. Demler, Phys. Rev. X 8, 011046 (2018a).
  * Giamarchi (2004) T. Giamarchi,  _Quantum physics in one dimension_, International series of monographs on physics (Clarendon Press, Oxford, 2004).
  * Preiss  _et al._ (2015) P. M. Preiss, R. Ma, M. E. Tai, A. Lukin, M. Rispoli, P. Zupancic, Y. Lahini, R. Islam, and M. Greiner, Science 347, 1229 (2015).
  * Schneider  _et al._ (2012) U. Schneider, L. Hackermüller, J. P. Ronzheimer, S. Will, S. Braun, T. Best, I. Bloch, E. Demler, S. Mandt, D. Rasch, and A. Rosch, Nature Physics 8, 213 (2012).
  * Grusdt  _et al._ (2018b) F. Grusdt, Z. Zhu, T. Shi, and E. Demler, SciPost Phys. 5, 57 (2018b).
  * Jepsen  _et al._ (2021) P. N. Jepsen, W. W. Ho, J. Amato-Grill, I. Dimitrova, E. Demler, and W. Ketterle, (2021), arXiv:2103.07866 [cond-mat.quant-gas] .
  * Hauschild and Pollmann (2018) J. Hauschild and F. Pollmann, SciPost Phys. Lect. Notes , 5 (2018).
  * Ogata and Shiba (1990) M. Ogata and H. Shiba, Phys. Rev. B 41, 2326 (1990).
  * Ren and Anderson (1993) Y. Ren and P. W. Anderson, Phys. Rev. B 48, 16662 (1993).
  * Zaanen  _et al._ (2001) J. Zaanen, O. Y. Osman, H. V. Kruis, Z. Nussinov, and J. Tworzydlo, Philosophical Magazine B 81, 1485 (2001).
  * Kruis  _et al._ (2004) H. V. Kruis, I. P. McCulloch, Z. Nussinov, and J. Zaanen, Phys. Rev. B 70, 075109 (2004).
  * Hilker  _et al._ (2017) T. A. Hilker, G. Salomon, F. Grusdt, A. Omran, M. Boll, E. Demler, I. Bloch, and C. Gross, Science 357, 484 (2017).
  * Bohrdt  _et al._ (2018) A. Bohrdt, D. Greif, E. Demler, M. Knap, and F. Grusdt, Phys. Rev. B 97, 125117 (2018).

##  Appendix A Hole-magnon bound state

###  A.1 Two-body problem

By doping a fully polarized insulator in a frustrated geometry with vacancies
(i.e. holes) and flipped spins (i.e. magnons) effective interactions between
them appear. In this section we consider an insulator polarized along the $z$
axis in a fermionic zigzag ladder and study the hole-magnon problem. Holes and
magnons can rigorously be created by employing the transformations
$\hat{c}_{\uparrow i}=\hat{h}^{\dagger}_{i}$ and $\hat{c}_{\downarrow
i}=\hat{h}^{\dagger}_{i}S^{+}_{i}$. In this way the original $2L-N_{h}$
particles problem can be reduced to a problem of $N_{h}$ holes and $N_{m}$
magnons, being $2L$ the total number of sites in the ladder. Under these
transformations the $t-J$ Hamiltonian becomes,

$\displaystyle\hat{H}_{tJ}$ $\displaystyle=$ $\displaystyle\pm
t\sum_{i,\,e}\left(\hat{h}^{\dagger}_{i+e}\hat{h}_{i}+\text{h.c.}\right)+J\sum_{i,\,e}\mathbf{S}_{i}\mathbf{S}_{i+e}$
(2) $\displaystyle\pm$ $\displaystyle
t\sum_{i,\,e}\left(\hat{h}^{\dagger}_{i+e}\hat{h}_{i}S^{-}_{i}S^{+}_{i+e}+\text{h.c.}\right),$

where $e=a,2a$ being $a$ the lattice spacing and the $+$ and $-$ signs
correspond to the fermionic and bosonic case, respectively. They appear
because of the anticonmutation and commutation relations satisfied by the hole
operators in each case. Under these transformations the fully polarized
insulating state $|\text{FP}\rangle$ becomes the vacuum of holes and magnons.
Holes and magnons can be created by applying the respective creation
operators. In particular the hole-magnon state is defined as
$\hat{h}^{\dagger}_{i}S^{-}_{j}|\text{FP}\rangle\equiv|ij\rangle$. We can
solve this two-particle problem by separating the center $R/a=(i+j)/2$ and
relative $z/a=i-j$ motion using the set of states

$|\psi\rangle=\sum_{ij}e^{iQR}\psi_{Q}(z)|ij\rangle,$ (3)

where we introduce the total quasi-momentum $Q$ of the pair and the
wavefunction in relative position space $\psi_{Q}(z)$. Since the interaction
part of the Hamiltonian (2) does not couple states with different quasi-
momentum $Q$ we can obtain the equation of motion for the relative part,

$\displaystyle E_{Q}\psi_{Q}(z)=\sum_{e=\pm a,\pm 2a}\left(\pm
te^{iQe/2}+\frac{J}{2}e^{-iQe/2}\right)\psi_{Q}(z+e)$ (4) $\displaystyle+$
$\displaystyle\sum_{e=\pm a,\pm 2a}\delta(z-e)\left(\pm
t\psi_{Q}(-e)+\frac{J}{2}\psi_{Q}(e)\right),$

where the energy of the pair $E_{Q}$ depends parametrically on the total
quasi-momentum. We measure the hole-magnon energy $E_{Q}$ with respect to the
energy of the ferromagnetic background. The hole-magnon interaction contains
an exchange term of strength $t$ and a bare nearest-neighbor interaction of
strength $J/2$. For fermions (bosons) the nearest-neighbor interaction $J>0$
($J<0$) has a repulsive (attractive) effect on the pair. Moreover, since
parity is a good quantum number of our model the exchange term can be cast
into a diagonal form $\psi_{Q}(-e)=\pm\psi_{Q}(e)$ by considering a symmetric
or antisymmetric state with respect to the hole-magnon exchange. In this way,
we notice that the exchange term has an attractive (repulsive) effect for the
antisymmetric (symmetric) solution in the fermionic model. In the bosonic
model we observe the opposite. Finally, since by definition we cannot have a
hole and a magnon in the same position we impose a hard-core constrain in the
solution $\psi_{Q}(0)=0$ by including a very strong on-site repulsion.

We numerically solve Eq. (4) and we obtain a hole-magnon bound state for the
fermionic model and an antibound state in the bosonic one, see Fig. 16. The
bound and antibound states are antisymmetric with respect to the exchange of
the hole and magnon positions. Thus we do not observe binding when the
solution is symmetric. The appearance of an antisymmetric bound (antibound)
state can be explained by the presence of the interexchange interaction which
leads to an effective attraction (repulsion) between the hole and the magnon
in the fermionic (bosonic) model.

Figure 19: Main panel: Binding energy $E_{B}$ (continuous line) and effective
mass $m^{*}$ (dots) of the hole-magnon bound state in the fermionic zigzag
ladder as a function of $J/t$. Inset: Band structure of the hole-magnon
problem as a function of the total quasimomentum $Q$ for $J/t=0.1$. The grey
area denotes the hole-magnon continuum and the continuous line the bound state
band. Dotted lines correspond to the Born-Oppenheimer approximation, see main
text.

The minimum of the binding energy appears at $Q=0$ for any value of $J/t$.
Therefore we define the binding energy at this quasi-momentum
$E_{B}=E_{1H1M}-E_{1H}-E_{1M}$, see Fig. 19. The binding energy is larger for
small values of $J/t$ and it vanishes close to $J\approx 2t$ which can be
easily interpreted as the point where the bare repulsion cancels the
interexchange term. Since binding is stronger at $J/t=0$ we can neglect the
magnon kinetic energy allowing us to perform a Born-Oppenheimer approximation.
First, we solve Eq. (4) for a static magnon and we obtain the hole energy
$E^{h}$ and the wavefunction $\psi_{Q}^{h}(z)$. Then, we estimate the bound
state energy by adding the kinetic energy of the magnon which is given by
$E^{m}_{Q}=-J(\Gamma(e)\cos(Q)+\Gamma(2e)\cos(2Q))$, where we define the
overlap $\Gamma(e)=|\langle\psi^{h}(e)|\psi^{h}(0)\rangle|\approx|e|/3$ being
$e=\pm a,\pm 2a$. This gives us,

$E_{Q}=E^{h}-\frac{J}{3}\left(\cos(Qa)+2\cos(2Qa)\right).$ (5)

For small values of $J/t$ we observe a nice agreement between the Born-
Oppenheimer approximation and the full solution of Eq. (4), see Fig. 19.
Moreover we can compute the effective mass of the hole-magnon bound state at
the minimum of the band $(2m^{*})^{-1}=|\frac{\partial^{2}E_{Q}}{\partial
Q^{2}}|$. From the Born-Oppenheimer approximation we obtain
$m^{*}=1/(3Ja^{2})$ which diverges for $J/t\rightarrow 0$ indicating the
immobile nature of the bound state, see Fig. 19.

###  A.2 Finite size scaling

Figure 20: Extension of the hole-magnon pair $\sqrt{\langle z^{2}\rangle}/a$
(a) as a function of $J/t$ for different system sizes. Probability of finding
a hole and a flipped spin at a distance equal to the lattice spacing $a$ (b)
as a function of $J/t$ for different system sizes.

In order to ensure the bound state nature of the hole-magnon pair for a range
of $J/t$ we have performed a finite size scaling. We have obtained the binding
energy, the extension of the pair and the probability of finding the hole and
the magnon in nearest sites for large systems sizes and then we have
extrapolated to the infinite size limit. The binding energy in the infinite
size limit becomes negative for $J/t<2$ and the probability $a|\psi(z=a)|^{2}$
becomes zero for $J/t>2$ as shown in Fig. 7 and Fig. 20. This signals a
transition from bound to non-bound pair in the infinite size limit at $J/t=2$.
Moreover the size of the pair is independent of the system size for $J/t<2$ in
the large size limit and it is strongly dependent on it for $J/t>2$. This also
supports the appearance of a hole-magnon bound state at $J/t<2$ with a
characteristic size independent of the system size.

##  Appendix B Antiferromagnetic spinbag model

Figure 21: Magnetization and in-plane spin-spin correlations (a) computed over
the projected wavefunction $|\psi_{h}\rangle$ for which we fix the hole
position to be at the center of the lattice. Hole density (b) of the
groundstate obtained from DMRG (blue dots) and from the variational
calculation (dashed line). Both panels correspond to a computation of the
fermionic zigzag ladder with a single hole and 20 magnons for $J/t=0$.

We now consider the situation of binding many magnons to a single hole thus
creating a spin polaron or spinbag. An analytical treatment of this problem
using the Hamiltonian (2) is very hard since the hole hopping disturbs the
spin background. Therefore we change to a representation where the hole lives
in the links of a new lattice with $2L-1$ sites i.e. the squeezed space Ogata
and Shiba (1990); Ren and Anderson (1993); Zaanen  _et al._ (2001); Kruis  _et
al._ (2004); Hilker  _et al._ (2017); Bohrdt  _et al._ (2018). In this way the
hole hopping between the two legs of the ladder simply becomes
$\left(t\sum_{i}\hat{h}^{\dagger}_{i+1}\hat{h}_{i}+\text{h.c.}\right)$ since
it does not disturb the spin order in squeezed space. On the other hand, for
the hole hopping between two sites of the same leg we have to take into
account that two spins are exchanged in squeezed space. The operator doing
such process is
$\left(\hat{h}^{\dagger}_{i+2}\hat{h}_{i}+\hat{h}^{\dagger}_{i}\hat{h}_{i+2}\right)\left(\frac{1}{2}+2\,\mathbf{S}_{i}\mathbf{S}_{i+1}\right)$.
Concerning the spin interaction we have to take into account that two spins
will not interact if a hole was occupying one of the two sites in the original
lattice, giving the term
$J\mathbf{S}_{i}\mathbf{S}_{i+1}(1-\hat{h}^{\dagger}_{i}\hat{h}_{i})$ and
$J\mathbf{S}_{i}\mathbf{S}_{i+2}(1-\hat{h}^{\dagger}_{i}\hat{h}_{i})(1-\hat{h}^{\dagger}_{i+1}\hat{h}_{i+1})$.
Moreover if a hole is present in the original lattice between two spins these
two will be nearest neighbor in squeezed space and they will interact through
the the superexchange coupling
$J\mathbf{S}_{i}\mathbf{S}_{i+1}\hat{h}^{\dagger}_{i}\hat{h}_{i}$. Therefore
the Hamiltonian in squeezed space becomes,

$\displaystyle\hat{H}_{tJ}$ $\displaystyle=$ $\displaystyle
t\sum_{i}\left(\hat{h}^{\dagger}_{i+1}\hat{h}_{i}+\text{h.c.}\right)+J\sum_{i}\mathbf{S}_{i}\mathbf{S}_{i+1}$
$\displaystyle+$ $\displaystyle
t\sum_{i}\left(\hat{h}^{\dagger}_{i+2}\hat{h}_{i}+\hat{h}^{\dagger}_{i}\hat{h}_{i+2}\right)\left(\frac{1}{2}+2\,\mathbf{S}_{i}\mathbf{S}_{i+1}\right)$
$\displaystyle+$ $\displaystyle
J\sum_{i}\mathbf{S}_{i}\mathbf{S}_{i+2}(1-\hat{h}^{\dagger}_{i}\hat{h}_{i})(1-\hat{h}^{\dagger}_{i+1}\hat{h}_{i+1}).$

From this Hamiltonian we notice that the hole hopping in the same leg
($\hat{h}^{\dagger}_{i+2}\hat{h}_{i}$) provides a strong coupling of strength
$t$ between the hole and the magnons. Notice that this coupling induces strong
antiferromagnetic correlations between the two legs even though the original
Ising interaction is not present for $J/t=0$. These antiferomagnetic
correlations between the two legs will appear when a hole is present and the
direct hopping dominates over the superexchange interaction $t/J\gg 1$. The
original superexchange interaction presents magnetic frustration and at the
semiclassical level it favours the formation of a spiral order of angle
$\theta=\arccos(-1/4)$. Since we work at a fixed total mangetization (i.e. a
fixed number of magnons) we expect that the system can be understood as
consisting of two regions: Region (I) is adjacent to the hole and has spins in
the XY plane with antiferromagnetic correlations between the legs. Region (II)
extends outside of region (I), has nearly full $S^{z}$ polarization and spiral
winding of the small XY components of the spins, see Fig. 3. We name this
object antiferromagnetic spinbag or antiferromagnetic spin polaron Nagaev
(1992) since it is analogous to the ferromagnetic one Auerbach (1994) but its
origin resides in the frustration of the system.

To obtain qualitative expressions for the size and energy of the
antiferromagnetic spinbag we consider a semiclassical approximation of the
spin part of the wavefunction
$\mathbf{S}_{i}\mathbf{S}_{i+1}=\tilde{S}^{2}\left(\sin(\varphi)^{2}\cos(\theta)+\cos(\varphi)^{2}\right)$
with $\tilde{S}^{2}=S(S+1)=3/4$ The angles $\theta$ and $\varphi$ determine
the spin orientation in the XY plane and the angle respect to the $z$ axis,
respectively (see top panel of Fig. 21). Therefore we can relate the total
number of magnons $N_{m}$ and the size of the antiferromagnetic spinbag $2R$
using the angle $\varphi_{II}$,

$\cos(\varphi_{II})\left(L-R\right)=L-N_{m}.$ (7)

By considering the limit $J/t\ll 1$ we assume a frozen magnetic order and
solve the equation of motion for the hole with a fixed spin background along
the lines of the Born-Oppenheimer approximation. The hole equation of motion
in each region is given by,

$\displaystyle(I):\,E_{h}h_{i}$ $\displaystyle=$ $\displaystyle
t\left(h_{i+1}+h_{i-1}\right)-\frac{3J}{4}\left(h_{i}+h_{i+1}\right)$ (8)
$\displaystyle-$ $\displaystyle t\left(h_{i-2}+h_{i+2}\right),$
$\displaystyle(II):\,E_{h}h_{i}$ $\displaystyle=$ $\displaystyle
t\left(h_{i+1}+h_{i-1}\right)$ $\displaystyle+$
$\displaystyle\frac{3J}{32}\left(15\cos(\varphi_{II})^{2}-7\right)\left(h_{i}+h_{i+1}\right)$
$\displaystyle+$
$\displaystyle\frac{t}{8}\left(15\cos(\varphi_{II})^{2}+1\right)\left(h_{i-2}+h_{i+2}\right).$

These equations show that the hole energy is reduced when it is localized
inside the antiferromagnetic region (I). Thus the antiferromagnetic background
acts as an effective square well potential to the hole. In order to obtain an
analytical expression for the hole energy $E_{h}$ we solve the hole equation
of motion by taking the continuum limit and proposing a gaussian wavefunction
localized in the antiferromagnetic region
$h_{i}=\exp\\{-x_{i}^{2}/(R^{2})\\}/\sqrt{R\sqrt{\pi/2}}$, where we employ the
size of the bag $2R$. Moreover Eq. (7) establishes the relation between the
size of the bag $R$ and the angle $\varphi_{II}$ when the total number of
magnons is specified. Taking the limit $L\ll R\sim N_{m}\ll 1$ we obtain a
hole energy that scales with the number of magnons as,

$E_{h}\propto 4.28\frac{t}{N_{m}^{2}}-0.04\frac{J}{N_{m}}.$ (10)

We now consider the magnetic energy cost of creating an antiferromagnetic
spinbag. We compare the magnetic energies with the same number of magnons when
there is an antiferromagnetic spinbag and when there is not. In the first case
the antiferromagnetic region has zero magnetic energy and the spiral region
has an energy $E_{M1}$ which is determined by the angle $\varphi_{II}$ (see
Eq. (7)) and the size of the spinbag $R$. In the second case the full ladder
presents a spiral order with energy $E_{M2}$ and an angle satisfying
$\cos(\varphi)L=L-N_{m}$. The difference between these two magnetic energies
in the limit $L\gg R\gg 1$ is,

$E_{M1}-E_{M2}\propto\frac{51}{8}JN_{m}.$ (11)

The magnetic energy increases when creating an antiferromagnetic spinbag
because the angle $\varphi$ of the spiral order is reduced. The total energy
of the antiferromagnetic spinbag is given by,

$E_{ASB}=4.28\frac{t}{N_{m}^{2}}-0.04\frac{J}{N_{m}}+\frac{51}{8}JN_{m},$ (12)

where we discard the constant terms with respect to the number of magnons.
From Eq. (12) we observe that increasing the size of the antiferromagnetic
spinbag reduces the hole kinetic energy but increases the magnetic one.
Therefore there is an optimal number of magnons forming the antiferromagnetic
spinbag $N_{m}^{*}$. Considering the situation of a small antiferromagnetic
region compared with the spiral region $L\gg R\sim N_{m}\gg 1$ we can estimate
the number of magnons,

$N_{m}^{*}\propto 1.10\,(t/J)^{1/3},$ (13)

considering the condition $(J/t)\ll 1$.

The semiclassical spinbag theory predicts that under doping the fermionic
zigzag ladder with a single hole and magnons, a certain number of magnons,
given by Eq. (13), accumulate around the hole with an antiferromagnetic
ordering and the remaining ones are pushed away from this region forming a
spiral order with a net magnetization pointing in the $z$ axis, see panel (a)
of Fig. 21. In order to benchmark our theory we have performed DMRG
simulations of the fermionic zigzag ladder with a single hole and different
number of magnons. We observe that for small values of $J/t$, the hole density
accumulates in a certain region of the lattice, see panel (b) of Fig. 21. This
can be attributed to the large effective mass of the antiferromagnetic spinbag
and the fact that we are working with open boundary conditions thus explicitly
breaking the translational invariance of the system. By computing the
projected wavefunction $\hat{P}^{h}_{L}|\psi\rangle=|\psi_{h}\rangle$ which
has a hole fixed at the center of the lattice we obtain the spin ordering
surrounding the hole White and Scalapino (1997). Around the hole the
magnetization suddenly decays indicating that the spins tilt to the XY plane.
Moreover the in-plane spin-spin correlations show an antiferromagnetic order
inside this region. Outside this region we recover a perfect ferromagnetic
spin order. This confirms the formation of an antiferromagnetic spinbag for
small values of $J/t$. Notice that the simulation is performed at $J/t=0$.
Therefore all the magnons are contained in the antiferromagnetic spinbag and
the spiral order simply becomes a ferromagnetic order since $\varphi=0$.