# Magnonic Superconductivity

Khachatur G. Nazaryan  Liang Fu  Department of Physics, Massachusetts
Institute of Technology, Cambridge, MA 02139

###### Abstract

We uncover a new superconducting state with partial spin polarization induced
by a magnetic field. This state, which we call ”magnonic superconductor”,
lacks a conventional pairing order parameter, but is characterized instead by
a composite order parameter that represents the binding of electron pairs and
magnons. We rigorously demonstrate the existence of magnonic superconductivity
with high transition temperature in a triangular lattice Hubbard model with
repulsive interaction. We further show that magnonic Cooper pairs can attract
to form higher-charge bound states, which can give rise to charge-$6e$
superconductivity.

Unraveling the mechanism of high-temperature superconductivity in doped
cuprates and other correlated insulators is a long-standing challenge in
condensed matter physics [1, 2]. An appealing idea is the resonating valence
bond (RVB) theory [3, 4], originally proposed for parent Mott insulators that
are quantum spin liquids formed from highly entangled spin-singlet pairs. More
recently, unconventional superconductivity has been discovered in a variety of
two-dimensional (2D) materials, which emerges from doping insulating states.
These include ZrNCl [5], WTe2 [6], and crystalline graphene multilayers [7,
8]. While undoped cuprates are antiferromagnetic Mott insulators, the parent
states in these 2D materials are of a different nature. Undoped ZrNCl, WTe2
and crystalline graphite are all band insulators with an even number of
electrons per unit cell. Upon doping, the insulating state in ZrNCl and WTe2
evolves directly into the superconducting state. Unconventional
superconductivity is evidenced by doping-induced BEC/BCS crossover in ZrNCl
and strong violation of Pauli limit in WTe2.

The emergence of unconventional superconductivity from doping band insulators
in 2D materials invites new ideas on the pairing mechanism [9, 10, 11, 12, 13,
14]. An interesting possibility is that two doped charges form a bound state
in the insulating background. Counter-intuitively, recent works have shown
that such two-particle bound state—often with non-$s$-wave pairing
symmetry—can be formed due to interband fluctuations in band insulators with
repulsive electron-electron interaction [15, 16, 17, 18, 19]. This mechanism
of electron attraction from repulsion gives rise to unconventional
superconductivity upon doping.

In this work, we demonstrate a novel kind of pairing and superconductivity
that emerges from doping spin-polarized insulators under a magnetic field. In
our theory superconductivity is formed from the condensation of charge-$2e$
multiparticle bound states made of two doped charges and one or more magnons
(spin-flips), which carry spin $S\geq 2$ distinct from Cooper pairs. We call
such bound state a ”magnonic Cooper pair”, and the corresponding
superconductor a “magnonic superconductor”. Importantly, a magnonic
superconductor does not have a standard pairing order parameter $\langle
c^{\dagger}_{i\alpha}c^{\dagger}_{j\beta}\rangle=0$, but is characterized by a
composite order parameter, for example, $\langle
c^{\dagger}_{i\uparrow}c^{\dagger}_{j\uparrow}c^{\dagger}_{k\uparrow}c_{l\downarrow}\rangle\neq
0$.

As a concrete example, we demonstrate the emergence of composite pairing and
magnonic superconductivity in a triangular lattice Hubbard model over a wide
range of magnetic field. Starting with the fully spin-polarized insulator at
half filling, we study this system doped with two holes as a function of the
magnetic field and find numerically exact solution for the ground state. Below
a critical field, the system transitions into a non-fully polarized state with
one or more magnons present, the regime of our interest.

Over a wide range of magnetic field below the saturation field, we find that
the lowest energy state can be a two-holes-one-magnon (2H1M) bound state,
i.e., with the involvement of one magnon, two doped holes attract. This
attraction occurs when the tight-binding energy dispersion has three
degenerate maxima and one minimum. Remarkably, this magnonic Cooper pair
features a large binding energy comparable to the hopping amplitude and a
small size of a few lattice constants. The condensation of these spin-$2$
charge-$2e$ bound states gives rise to magnonic superconductivity.

We note that related ideas on composite pair superconductivity have been
considered in Kondo lattice systems [20], where spin-$0$ Cooper pair may
couple to a spin-$1$ magnon to form a spin-$1$ charge-$2e$ order parameter.
This is different from our case of spin-$2$ composite pairing. Most closely
related to our work is Ref.[21], which found similar composite pairing in a
different model, albeit only in a very narrow range of magnetic field. Our
work further identifies a general mechanism for 2H1M bound state formation due
to holes having three valley flavors. Finally, by analyzing the interactions
between these 2H1M charge-$2e$ bosons, we demonstrate the formation of three-
boson bound states in a certain parameter regime, which gives rise to
charge-$6e$ superconductivity.

Triangular lattice Hubbard model under magnetic field. For concreteness, we
study the Hubbard model on the triangular lattice with nearest and next
nearest neighbor hoppings $t_{1}$ and $t_{2}$, subject to a magnetic field $h$
that couples to electron spin:

$\displaystyle H=$ $\displaystyle-
t_{1}\sum_{\sigma;\langle\textbf{i},\textbf{j}\rangle}\left(c_{\textbf{i}\sigma}^{\dagger}c_{\textbf{j}\sigma}+h.c.\right)-t_{2}\sum_{\sigma;\langle\langle\textbf{i},\textbf{j}\rangle\rangle}\left(c_{\textbf{i}\sigma}^{\dagger}c_{\textbf{j}\sigma}+h.c.\right)$
$\displaystyle+U\sum_{\textbf{i}}n_{\textbf{i}\uparrow}n_{\textbf{i}\downarrow}-\frac{h}{2}\sum_{\textbf{i}}\left(n_{\textbf{i}\uparrow}-n_{\textbf{i}\downarrow}\right).$
(1)

As we will show later, the essential physics of magnonic Cooper pairing
remains robust against additional perturbations such as longer-range hopping
or interaction.

Figure 1: Phase diagram of the infinite-$U$ Hubbard model on the triangular
lattice at half-filling with two doped holes, as a function of the magnetic
field $h$ and tight-binding parameters $t_{2}/t_{1}$. Various multi-particle
bound states of holes and magnons are found below the saturation field.

Exactly at half-filling (with one electron per site), the ground state is a
fully spin-polarized insulator above a critical field denoted as $h_{J}$. In
the strong-coupling regime (large $U$), $h_{J}$ is determined by the spin
exchange interaction $J_{a}=4t^{2}_{a}/U$ ($a=1,2$) and therefore is small.
When one or more holes are introduced, however, a much larger saturation field
(denoted as $h_{1}$) comparable to the hopping energy scale is required to
fully polarize all the spins in the system. This is because the kinetic motion
of holes is affected by the presence of magnons. This kinetic mechanism also
leads to various bound states of holes and magnons, as we show below.

We first review the case of one doped hole. Starting with the seminal work of
Zhang, Zhu and Batista [21], the problem of one doped hole in two-dimensional
triangular lattice Hubbard model has been thoroughly studied in recent years
[22, 23]. Below a critical field, the ground state with one hole is a non-
fully polarized state that contains a spin polaron, a one-hole-one-magnon
bound state with spin $S=3/2$ and a large binding energy on the order of the
hopping amplitude.

For the triangular lattice Hubbard model with only nearest-neighbor hopping
($t_{2}=0$), the spin polaron is shown to be the lowest-energy-per-charge
excitation over a wide range of field. Consequently, the ground state at
finite doping density is a spin polaron liquid that exhibits a magnetization
plateau [22, 24], which has been experimentally observed in semiconductor
moiré superlattices [25].

For the purpose of studying pairing and superconductivity, this work is
concerned with the system of two doped holes. Importantly, we consider the
Hubbard model with both nearest and second nearest neighbor hoppings. As we
will show below, $t_{2}\neq 0$ plays a crucial role in allowing for the
binding of two holes and magnons.

Figure 2:  (a), (b) Many-body spectrum of the infinite-$U$ Hubbard model in
the presence of two doped holes and one spin flip as a function of COM
momentum along $M$-$\Gamma$-$K$ direction for $t_{2}/t_{1}=0.15$ and $0.11$,
respectively. The lowest energy branch (blue color) corresponds to the two-
holes-one-magnon bound state, which is separated from the high-energy
continuum (orange color) by a large gap. Also shown on the right are the
ground state wavefunction at fixed positions of the spin-flip (down arrow) and
the hole (white circle), and the hole-hole correlation function $c^{hh}$.

We use the exact diagonalization method to calculate the energy spectrum of
the system with two holes and $m$ magnons, denoted as $2HmM$. First, we
diagonalize the many-body Hamiltonian $H$ at zero field $h=0$ to find the
lowest energy state with a fixed number of magnons, whose energy is denoted as
$\epsilon^{2HmM}$. By utilizing the center-of-mass momentum as a good quantum
number, we obtain results for systems up to $N=36^{2}=1296$ sites ($m=1$),
$N=6\times 12=72$ sites ($m=2$), and $N=6\times 6=36$ sites ($m=3$). In the
presence of magnetic field, the energies of states with two holes and $m$
magnons are given by

$\displaystyle E^{2HmM}=\epsilon^{2HmM}+(m+1)h,$ (2)

where the second term comes from the Zeeman energy cost associated with holes
and spin-flips.

Figure 3:  (a) The inverse effective mass in units of $t_{1}a^{2}/\hbar^{2}$
(red) and average distance between holes (blue) of the magnonic Cooper pair
(2H1M). The average distance was calculated on lattices $12\times 12$,
$18\times 12$ and $18\times 18$, showing no sensitive dependence on the size,
indicating a well-localized state. (b) Energy gap of 4H2M state vs 2 $\times$
2H1M. For $t_{2}/t_{1}$ within $(0.09,0.15)$ the magnonic Cooper pairs attract
and form a composite particles, while for larger $t_{2}/t_{1}$ they repel.

For a given magnetic field $h$, by comparing these energies of different
states, we determine the number of magnons $m$ which minimizes the energy. The
full phase diagram of the ground state with two holes thus obtained is
presented in the Fig. 1, as a function of $t_{2}/t_{1}$ and the magnetic field
$h$. As the magnetic field is reduced, distinct hole-magnon bound states with
different numbers of magnons are found. For a wide range of field below the
saturation field denoted as $h_{2}$, the ground state of two holes is either
two unbound spin polarons (1H1M) or a 2H1M bound state, depending on the tight
binding model parameters. Below a lower critical field, bound states with two
or more magnons appear. Below we will analyze various regions in detail.

Our findings reveal that for small enough $t_{2}/t_{1}$, below the saturation
field the fully polarized state (FP) transitions into a state with two
separate charge-$e$ spin polarons, consistent with the previous works [22,
24]. In contrast, a small second-neighbor hopping $t_{2}\geq 0.09t_{1}$ is
enough to drive the system into a new phase: over a wide range of field, we
find a bound state of two holes and one magnon ($2H1M$). This establishes the
formation of a magnonic Cooper pair, a composite boson of charge $2e$ and spin
$S=2$.

We further calculate the many-body energy spectra of the system with two holes
and one magnon as a function of the many-body momentum, shown in Fig.2. For
two representative values of $t_{2}/t_{1}=0.11$ and $0.15$, the lowest energy
branch is well separated from the the continuum by a gap. This branch
corresponds to the energy dispersion of the $2H1M$ bound state as a function
of its center-of-mass (COM) momentum. Note that when $t_{2}/t_{1}$ surpasses
$0.15$, the ground state has COM momentum at $\Gamma$, while within the region
$0.09<t_{2}/t_{1}<1/8$, the COM momentum shifts to $M$ point. The binding
energy is remarkably large and on the order of the hopping energy. This
implies that the binding of two holes is robust against perturbations to the
model Hamiltonian, e.g., longer-range hoppings or electron-electron
repulsions.

The internal structure of the $2H1M$ bound state can be visualized by
examining its wavefunction, which depends on the positions of two holes and
one magnon. By fixing the magnon at the origin and one of the holes at a
neighboring site, the dependence of the wavefunction on the coordinate of the
second hole is depicted in Fig.2. Additionally, the hole-hole correlation
function $c^{hh}(\textbf{r})=\langle
n_{h}(\textbf{0})n_{h}(\textbf{r})\rangle$ clearly shows that the two holes
are tightly bound. Remarkably, the bound state size as defined by the average
distance between the two constituent holes is only around $3$ lattice
constants over a wide range of tight binding parameters, see Fig. 3.

Figure 4:  (a), (b) Spin up and spin down electron dispersions (separated by
Zeeman energy) for $t_{2}/t_{1}=0,0.3$, respectively. The electron band
minimum is at $\Gamma$ point. Meanwhile, the band maximum is double-degenerate
at $K,K^{\prime}$ points for $t_{2}/t_{1}<1/8$; it becomes triple-degenerate
for $t_{2}/t_{1}>1/8$ shifting to $M$ points. In (a) two holes are efficiently
placed at band maxima allowing for polaron formation. In (d), system can host
three holes at band maxima, leading to the formation of a magnonic Cooper
pair.

Insights into the origin of 2H1M bound state can be gained from considering
the kinetic energy dispersion, illustrated in the Fig. 4. We observe that
while the band minimum occurs at the $\Gamma$, the location of the band
maximum highly depends on the parameter $t_{2}/t_{1}$. For weak $t_{2}$, the
maximum is doubly degenerate, occurring at $K$ and $K^{\prime}$. However, when
$t_{2}/t_{1}$ exceeds 1/8, the maximum shifts to $M$ points leading to a
three-fold valley degeneracy. As a result, the system can host two holes in
any of the 3 band maxima. In the absence of any magnons (i.e. the fully
polarized phase), the energy of two-hole state in the Hubbard model is simply
$2\epsilon_{M}$, with $\epsilon_{M}$ being the energy of a hole at $M$ point.
In contrast, in the presence of a magnon the translation symmetry of the two
holes in the COM frame is broken, and the state of two holes can be in a
coherent superposition of three valleys of the form:
$|M_{1},M_{2}\rangle+|M_{2},M_{3}\rangle+|M_{3},M_{1}\rangle$. And such
coherent superposition lowers the energy of the system. This picture also
applies to the small $U$ limit. Here 2H1M state is a bound state of 3 spin
$\downarrow$ holes and 1 spin $\uparrow$ electron. We have numerically
verified that the transition to such 4 particle bound state occurs precisely
at $t_{2}/t_{1}=1/8$, and the results are presented in the supplementary
materials.

As the magnetic field is further reduced, the system continues to reveal its
richness. Notably, when $t_{2}/t_{1}$ surpasses $0.13$, the phase
characterized by two separate polarons (2 x 1H1M) transitions into a state
where the polarons combine into a composite bound state (2H2M). In this phase,
the charge carriers are bosons with a charge of 2$e$ and spin $S=3$. By
further decreasing magnetic field, the system transitions into states with an
increasing number of magnons. Figure 1 shows that the 2H2M phase has a
significantly narrower width compared to the 2H1M phase. As the field further
decreases, we anticipate additional transitions to states with a larger number
of magnons, eventually leading to magnon condensation.

Having established the phase diagram of two holes, we now discuss the effects
of finite doping density. The existence of magnonic Cooper pairs over a wide
range of magnetic fields naturally suggests the possibility of
superconductivity. At low density our system maps to an interacting Bose gas
where the boson corresponds to 2HmM bound states. If the interaction between
these bosons is repulsive, the ground state of the 2D Bose gas will be a
superfluid.

We now determine the sign of interaction of the magnonic Cooper pairs. We
perform ED study for a system with 4 holes and 2 magnons on a $6\times 6$
lattice and compare the energy with that of two separate 2H1M bound states.
The Fig. 3b) reveals that the interaction is attractive for $t_{2}/t_{1}$
within $(0.09,0.15)$ and then shifts to repulsive for larger $t_{2}/t_{1}$.

Therefore, for $t_{2}/t_{1}>0.15$, the ground state at low density is a
superconductor formed by Bose condensation of 2H1M magnonic Cooper pairs. It
is well known that the BKT superfluid transition temperature for a two
dimensional Bose gas with weak repulsive interactions can be estimated as [26,
27],
$k_{B}T_{\mathrm{BKT}}=\frac{n^{\text{b}}}{m}\frac{2\pi\hbar^{2}}{\ln\left(380\hbar^{2}/mU_{0}\right)}$.
Here $n^{\text{b}}=n/2$ is the boson density and $n$ the hole doping density,
$m$ the effective mass of bosons, $k_{B}$ is the Boltzmann constant and
$U_{0}$ is the contact repulsion strength between the bosons in the continuum.
For our Hubbard model $U_{0}$ can be extracted from the microscopic
parameters, and we expect it to be of order of hopping amplitude up to a
numerical prefactor. But because the dependence on $U_{0}$ is through a $\ln$,
the exact value should not affect the order of magnitude estimate of
$T_{BKT}$. The important dependence is on the effective mass $m$ of the 2H1M
state, which we determine from the many-body spectrum. The results on the Fig.
3a) show that for $t_{2}>0.15t_{1}$ the inverse mass $1/m$ is of order
$0.1t_{1}a^{2}/\hbar^{2}$. Therefore, the superconducting transition
temperature $k_{B}T_{BKT}$ is of the order of $0.1t_{1}\cdot\nu/\sqrt{3}$
where $\nu=na^{2}\sqrt{3}/2$ is the number of holes per site. We compare
$T_{BKT}$ with the Fermi energy, $E_{F}$, of the hole system in the fully
polarized state, which for low doping can be found as $E_{F}=2\pi/3\cdot
t_{1}\nu\sqrt{(9t_{2}/t_{1}-1)(1-t_{2}/t_{1})}$. Here the factor $1/3$ comes
from the three-fold degeneracy of the band minima ($M$ points). At a ratio of
$t_{2}/t_{1}=0.2$, $E_{F}$ is approximately $1.6\nu t_{1}$. Consequently, the
superconducting transition temperature $T_{BKT}$ is approximately $5\%$ of
$E_{F}/k_{B}$.

We now compare our findings to other mechanisms of strong coupling
superconductivity from repulsive interaction. Refs. [16, 17, 18] introduced a
mechanism of superconductivity mediated by virtual excitons in doped
insulators. In the narrow band limit, a controlled analysis shows that the
pairing energy scale is $t^{2}/\Delta$, where $\Delta$ represents the
insulating gap. In contrast, the pairing energy and $T_{BKT}$ in our case are
on the order of the hopping amplitude $t$, which is parametrically larger.

Similar to the spirit of our work, Ref. [21] found two-hole-one-magnon bound
state in a different triangular lattice Hubbard model with nearest and third-
neighbor hoppings. However, in their model this bound state only exists in a
narrow range of magnetic field. Our work reveals the microscopic origin of the
2H1M bound state formation due to the threefold valley degeneracy of hole
dispersion. Moreover, thanks to its large binding energy and stability range,
the 2H1M bound state in our model is robust against the introduction of
perturbations such as longer range hopping and longer range interaction.

For $t_{2}/t_{1}<0.15$, our ED calculation indicates that the two 2H1M bosons
attract, yielding three degenerate ground states with COM momentum at $M$
points, separated by a finite energy gap to higher energy states. We attribute
this to the fact that the ground state is a bound state of two 2H1M bosons in
different valleys, and therefore carries a total momentum $M_{1}+M_{2}=M_{3}$
etc. The presence of attraction between charge-$2e$ 2H1M bosons leads us to
conjecture that three magnonic Cooper pairs with total momentum at
$M_{1}+M_{2}+M_{3}=\Gamma$ will form a stable composite bound state. This will
give rise to a magnonic superconductivity through a Bose condensation of
charge $6e$ and spin $S=6$ bosons.

To explicitly illustrate this interesting scenario, we study an effective
model of dilute boson gas with three flavors corresponding to the three $M$
valleys. The Hamiltonian for this model takes the following form dictated by
symmetry:

$\displaystyle\mathcal{H}=\frac{1}{2}\sum_{\alpha,ij}\int
d\textbf{r}\;m_{\alpha,ij}^{-1}\psi_{\alpha}^{\dagger}\bm{\partial}_{i}\bm{\partial}_{j}\psi_{\alpha}$
$\displaystyle+\frac{1}{2}\int
d\bm{r}\;\left(\sum_{\alpha}U(\psi_{\alpha}^{\dagger})^{2}(\psi_{\alpha})^{2}-\sum_{\alpha\neq\beta}Vn_{\alpha}n_{\beta}\right)$
$\displaystyle+\frac{1}{2}\int
d\bm{r}\sum_{\alpha\neq\beta}J(\psi_{\beta}^{\dagger})^{2}(\psi_{\alpha})^{2}$
(3)

Here $\alpha,\beta=1,2,3$ denote the different $M$ valleys; the first term
describes the kinetic energy with anisotropic mass $m_{\alpha,ij}$ around the
$M_{\alpha}$ point. The second and third terms represent the intravalley and
intervalley density-density interactions. The third term is allowed by
symmetry and results in the scattering of two particles from one valley to
another valley. In the low density limit, only contact interactions are
relevant and their strengths are denoted as $U,V$ and $J$ (which are
determined by our microscopic Hubbard model parameters). Importantly, the
intervalley interaction $V$ is attractive and responsible for the binding
energy of two 2H1M bosons from different valleys. We further analyze the
effective model of bosons; details can be found in the Supplementary Material.
To determine what kind of boson bound states is energetically favorable, we
compare the binding energy per particle for two-boson and three-boson states.
Assuming that intravalley repulsion $U$ is repulsive and large, we find that
the three-boson bound states are more favorable than the two-boson ones and
even infinitesimal attraction $V<0$ is enough to form such bound states.

Therefore, we conclude that over a wide range of magnetic field, our Hubbard
model at low density is a dilute system of magnonic Cooper pairs, charge-$2e$
bosons consisting of two holes and one magnon. The boson energy dispersion may
have a unique minimum at $\Gamma$ or three degenerate minima at $M$. In the
former case the system is a charge-$2e$ magnonic superconductor, while in the
latter case the system exhibits charge-$6e$ superconductivity due to boson
triple condensation. Both magnonic superconductors are characterized by
composite order parameters involving both electrons and magnons. Last but not
the least, the binding between hole pairs and magnons into magnonic Cooper
pairs implies that the magnetization as a function of the magnetic field
exhibits a plateau where the total spin $S$ depends only on the doping
$\delta$:

$\displaystyle S=N(1-2\delta)/2$ (4)

with $N$ is the number of unit cells. This value differs from the fully spin-
polarized state which has a total spin $N(1-\delta)/2$ as well as the spin-
polaron metal, where the total spin is $N(1-3\delta)/2$ [24].

###### Acknowledgements.

We thank Aidan Reddy and Margarita Davydova for helpful discussions. This work
was supported by a Simons Investigator Award from the Simons Foundation.

## References

  * Lee  _et al._ [2006] P. A. Lee, N. Nagaosa, and X.-G. Wen, Doping a mott insulator: Physics of high-temperature superconductivity, Rev. Mod. Phys. 78, 17 (2006).
  * Norman [2011] M. R. Norman, The challenge of unconventional superconductivity, Science 332, 196 (2011).
  * Anderson [1973] P. W. Anderson, Materials Research Bulletin 8, 153 (1973).
  * Fazekas and Anderson [1974] P. Fazekas and P. W. Anderson, On the ground state properties of the anisotropic triangular antiferromagnet, The Philosophical Magazine: A Journal of Theoretical Experimental and Applied Physics 30, 423 (1974).
  * Nakagawa  _et al._ [2021] Y. Nakagawa, Y. Kasahara, T. Nomoto, R. Arita, T. Nojima, and Y. Iwasa, Gate-controlled bcs-bec crossover in a two-dimensional superconductor, Science 372, 190 (2021).
  * Fatemi  _et al._ [2018] V. Fatemi, S. Wu, Y. Cao, L. Bretheau, Q. D. Gibson, K. Watanabe, T. Taniguchi, R. J. Cava, and P. Jarillo-Herrero, Electrically tunable low-density superconductivity in a monolayer topological insulator, Science 362, 926 (2018).
  * Zhou  _et al._ [2021] H. Zhou, T. Xie, T. Taniguchi, K. Watanabe, and A. F. Young, Superconductivity in rhombohedral trilayer graphene, Nature 598, 434 (2021).
  * Zhou  _et al._ [2022] H. Zhou, L. Holleis, Y. Saito, L. Cohen, W. Huynh, C. L. Patterson, F. Yang, T. Taniguchi, K. Watanabe, and A. F. Young, Isospin magnetism and spin-polarized superconductivity in bernal bilayer graphene, Science 375, 774 (2022).
  * Dong  _et al._ [2023] Z. Dong, L. Levitov, and A. V. Chubukov, Superconductivity near spin and valley orders in graphene multilayers, Phys. Rev. B 108, 134503 (2023).
  * Cea  _et al._ [2022] T. Cea, P. A. Pantaleón, V. o. T. Phong, and F. Guinea, Superconductivity from repulsive interactions in rhombohedral trilayer graphene: A kohn-luttinger-like mechanism, Phys. Rev. B 105, 075432 (2022).
  * Chou  _et al._ [2022] Y.-Z. Chou, F. Wu, J. D. Sau, and S. Das Sarma, Acoustic-phonon-mediated superconductivity in bernal bilayer graphene, Phys. Rev. B 105, L100503 (2022).
  * Ghazaryan  _et al._ [2021] A. Ghazaryan, T. Holder, M. Serbyn, and E. Berg, Unconventional superconductivity in systems with annular fermi surfaces: Application to rhombohedral trilayer graphene, Phys. Rev. Lett. 127, 247001 (2021).
  * Chatterjee  _et al._ [2022] S. Chatterjee, T. Wang, E. Berg, and M. P. Zaletel, Inter-valley coherent order and isospin fluctuation mediated superconductivity in rhombohedral trilayer graphene, Nature Communications 13, 6013 (2022).
  * Sous  _et al._ [2023] J. Sous, C. Zhang, M. Berciu, D. R. Reichman, B. V. Svistunov, N. V. Prokof’ev, and A. J. Millis, Bipolaronic superconductivity out of a coulomb gas, Phys. Rev. B 108, L220502 (2023).
  * Slagle and Fu [2020] K. Slagle and L. Fu, Charge transfer excitations, pair density waves, and superconductivity in moiré materials, Phys. Rev. B 102, 235423 (2020).
  * Crépel and Fu [2021] V. Crépel and L. Fu, New mechanism and exact theory of superconductivity from strong repulsive interaction, Science Advances 7, eabh2233 (2021).
  * Crépel and Fu [2022] V. Crépel and L. Fu, Spin-triplet superconductivity from excitonic effect in doped insulators, Proceedings of the National Academy of Sciences 119, e2117735119 (2022).
  * Crépel  _et al._ [2022] V. Crépel, T. Cea, L. Fu, and F. Guinea, Unconventional superconductivity due to interband polarization, Phys. Rev. B 105, 094506 (2022).
  * He  _et al._ [2023] Y. He, K. Yang, J. B. Profe, E. J. Bergholtz, and D. M. Kennes, Superconductivity of repulsive spinless fermions with sublattice potentials, Phys. Rev. Res. 5, L012009 (2023).
  * Abrahams  _et al._ [1995] E. Abrahams, A. Balatsky, D. J. Scalapino, and J. R. Schrieffer, Properties of odd-gap superconductors, Phys. Rev. B 52, 1271 (1995).
  * Zhang  _et al._ [2018] S.-S. Zhang, W. Zhu, and C. D. Batista, Pairing from strong repulsion in triangular lattice hubbard model, Phys. Rev. B 97, 140507 (2018).
  * Davydova  _et al._ [2023] M. Davydova, Y. Zhang, and L. Fu, Itinerant spin polaron and metallic ferromagnetism in semiconductor moiré superlattices, Phys. Rev. B 107, 224420 (2023).
  * Morera  _et al._ [2023] I. Morera, M. Kanász-Nagy, T. Smolenski, L. Ciorciaro, A. m. c. Imamoğlu, and E. Demler, High-temperature kinetic magnetism in triangular lattices, Phys. Rev. Res. 5, L022048 (2023).
  * Zhang and Fu [2023] Y. Zhang and L. Fu, Pseudogap metal and magnetization plateau from doping moiré Mott insulator, SciPost Phys. Core 6, 038 (2023).
  * Tao  _et al._ [2024] Z. Tao, W. Zhao, B. Shen, T. Li, P. Knüppel, K. Watanabe, T. Taniguchi, J. Shan, and K. F. Mak, Observation of spin polarons in a frustrated moiré hubbard system, Nature Physics (2024).
  * Prokof’ev  _et al._ [2001] N. Prokof’ev, O. Ruebenacker, and B. Svistunov, Critical point of a weakly interacting two-dimensional bose gas, Phys. Rev. Lett. 87, 270402 (2001).
  * Roy  _et al._ [2021] A. Roy, M. Ota, A. Recati, and F. Dalfovo, Finite-temperature spin dynamics of a two-dimensional bose-bose atomic mixture, Phys. Rev. Res. 3, 013161 (2021).

Supporting Material for “High-temperature magnonic superconductivity under
magnetic field”

by Kh. G. Nazaryan, L. Fu

##  A Free particle spectrum

In the main text, we discussed the main mechanism that gives an origin to the
magnonic Cooper pair formation. In this section, we detail the single particle
band structure.

Following our main discourse, we start our from a sufficiently strong magnetic
field, which imposes a ferromagnetic ordering into the system. In the case of
half-filling, the ferromagnetic ground state is expressed as

$\displaystyle\left|f\right\rangle=\prod_{\textbf{m}}c_{\textbf{m}\uparrow}^{\dagger}|0\rangle.$
(A.5)

Here, $|0\rangle$ symbolizes the vacuum state. We have oriented the magnetic
field, and consequently, all the spins up.

The system is then doped with a hole or an extra electron. The extra electron
will have a spin down and be located at a site already occupied by a spin down
electron.

The dispersion of a single hole/electron doped on top of the half-filled
system can be easily obtained as

$\displaystyle\epsilon(\bm{\kappa})_{h/e}=h/2+U\delta_{e}\pm\mu\mp
2\left(t_{1}\epsilon_{1}(\textbf{k})+t_{2}\epsilon_{2}(\textbf{k})\right),$
(A.6)
$\displaystyle\epsilon_{1}(\textbf{k})=\cos(k_{x}a)+\cos\left(k_{x}a/2+k_{y}a\sqrt{3}/2\right)+\cos\left(k_{x}a/2-k_{y}a\sqrt{3}/2\right),$
(A.7)
$\displaystyle\epsilon_{2}(\textbf{k})=\cos(k_{y}a\sqrt{3})+\cos\left(3k_{x}a/2-k_{y}a\sqrt{3}/2\right)+\cos\left(3k_{x}a/2+k_{y}a\sqrt{3}/2\right),$
(A.8)

where $\delta_{e}$ indicates that the $U$ term emerges only for the doped
electron (due to double occupancy); $\mu$ is the chemical potential. The
electron dispersion is illustrated in the Fig. S1. We infer from this that
while the band minimum is always located at $\Gamma$ point, the band maximum
is strongly dependent on the parameter $t_{2}/t_{1}$. For weak $t_{2}$, the
maximum is double-degenerate and is at $K,K^{\prime}$, but above a threshold
value of $\tau=1/8$, it becomes triple degenerate and moves to $M$ points.

Figure S1: (a) Schematic represeantation of the triangular-lattice Hubbard
model with the nearest-neighbor, $t_{1}$, and next-nearest-neighbor, $t_{2}$,
hopping parameters. (b), (c) Are contour plots, illustrating the free electron
dispersion within the first BZ for hopping ratios $t_{2}/t_{1}=0,0.3$,
respectively. The electron band minimum is at $\Gamma$ point. Meanwhile, the
band maximum is double-degenerate at $K,K^{\prime}$ points for
$t_{2}/t_{1}<1/8$; it becomes triple-degenerate for $t_{2}/t_{1}>1/8$ shifting
to $M$ points.

##  B Many-body problem

We are interested in the interplay of two holes doped into the system. In the
context of supercritical magnetic fields, the ground state is expressed as
follows,

$\displaystyle|\Phi^{2h}\rangle=\sum_{\textbf{k},\textbf{l}}\Phi_{\textbf{kl}}|\Phi_{\textbf{kl}}\rangle,\quad|\Phi_{\textbf{kl}}\rangle=c_{\textbf{k}\uparrow}c_{\textbf{l}\uparrow}|f\rangle.$
(B.9)

It is an arbitrary linear superposition of states each containing two holes
located at cites k and l; and the coefficients $\Phi_{\textbf{kl}}$ embody the
two-particle wave function.

Upon lowering the magnetic field, in alignment with the principal discourse,
the system transitions to a new ground state. In this state, one of the spin
up electrons flips its spin and hops to another site, creating a hole and a
doublon. The wave function corresponding to this excitation is given by,

$\displaystyle|\Psi^{3h1e}\rangle=\sum_{\textbf{i},\textbf{j},\textbf{k},\textbf{l}}\Psi_{\textbf{ijkl}}|\Psi_{\textbf{ijkl}}\rangle,\quad|\Psi_{\textbf{ijkl}}\rangle=c_{\textbf{i}\downarrow}^{\dagger}c_{\textbf{j}\uparrow}c_{\textbf{k}\uparrow}c_{\textbf{l}\uparrow}|f\rangle$
(B.10)

###  B.1  Shrodinger’s Equation

Here we detail the analysis of the 4 particle (1 doublon and 3 holes) state
discussed in the main text.

We act at this state with our Hamiltonian and derive,

$\displaystyle-
t_{1}\sum_{\textbf{i},\textbf{j},\textbf{k},\textbf{l}}\sum_{\bm{\delta}\in\bm{\delta}^{(1)}}\Psi_{\textbf{ijkl}}\bigg{[}|\Psi_{\textbf{i}-\bm{\delta}\textbf{j}\textbf{k}\textbf{l}}\rangle-|\Psi_{\textbf{i}\textbf{j}+\bm{\delta}\textbf{k}\textbf{l}}\rangle-|\Psi_{\textbf{i}\textbf{j}\textbf{k}+\bm{\delta}\textbf{l}}\rangle-|\Psi_{\textbf{i}\textbf{j}\textbf{k}\textbf{l}+\bm{\delta}}\rangle+\left[\bm{\delta}\to-\bm{\delta}\right]\bigg{]}$
$\displaystyle-
t_{2}\sum_{\textbf{i},\textbf{j},\textbf{k},\textbf{l}}\sum_{\bm{\delta}\in\bm{\delta}^{(2)}}\Psi_{\textbf{ijkl}}\bigg{[}|\Psi_{\textbf{i}-\bm{\delta}\textbf{j}\textbf{k}\textbf{l}}\rangle-|\Psi_{\textbf{i}\textbf{j}+\bm{\delta}\textbf{k}\textbf{l}}\rangle-|\Psi_{\textbf{i}\textbf{j}\textbf{k}+\bm{\delta}\textbf{l}}\rangle-|\Psi_{\textbf{i}\textbf{j}\textbf{k}\textbf{l}+\bm{\delta}}\rangle+\left[\bm{\delta}\to-\bm{\delta}\right]\bigg{]}$
$\displaystyle+U\sum_{\textbf{i},\textbf{j},\textbf{k},\textbf{l}}\left(1-\delta_{\textbf{i}\textbf{j}}-\delta_{\textbf{i}\textbf{k}}-\delta_{\textbf{i}\textbf{l}}\right)\Psi_{\textbf{ijkl}}|\Psi_{\textbf{i}\textbf{j}\textbf{k}\textbf{l}}\rangle=\left(\epsilon^{3h1e}(t_{1},t_{2})+2\mu+2h\right)\sum_{\textbf{i},\textbf{j},\textbf{k},\textbf{l}}\Psi_{\textbf{ijkl}}|\Psi_{\textbf{i}\textbf{j}\textbf{k}\textbf{l}}\rangle$
(B.1)

Here $\bm{\delta}^{(1)}$ are the three basis vectors on triangular lattice
(vectors to nearest neighbors), and $\bm{\delta}^{(2)}$ are the vectors to
next-nearest neighbors.

We note that the Pauli exclusion principle leads the antisymmetry of the wave
function under the swap of two holes,

$\displaystyle|\Psi_{\textbf{ijkl}}\rangle=-|\Psi_{\textbf{ikjl}}\rangle=-|\Psi_{\textbf{ijlk}}\rangle$
(B.2)

We introduce a linear coordinate that encompasses the position of each
particle on the two-dimensional lattice,

$\displaystyle\alpha_{m}=L(x_{m}-1)+(y_{m}-1)+1,\quad\alpha_{m}\in[1,L^{2}],\quad
m=i,j,k,l.$ (B.3)

The initial 2D coordinates of the given particles are then identically
recovered from these linear coordinates. Then we can generate basis states for
which the linear coordinate of the hole j is smaller than that of hole k,
$\alpha_{j}<\alpha_{k}$, and similarly, $\alpha_{k}<\alpha_{l}$,

$\displaystyle\alpha_{j}<\alpha_{k}<\alpha_{l}$ (B.4)

By Pauli exclusion principle, these states span the entire Hilbert space.

We then proceed by projecting the Eq. (B.1) onto the bra- state to obtain a
set of equations on the coefficients $\Psi_{\textbf{ijkl}}$:

$\displaystyle-
t_{1}\sum_{\bm{\delta}\in\bm{\delta}^{(1)}}\bigg{[}\left(\Psi_{\textbf{i}+\bm{\delta}\textbf{j}\textbf{k}\textbf{l}}-\Psi_{\textbf{i}\textbf{j}-\bm{\delta}\textbf{k}\textbf{l}}-\Psi_{\textbf{i}\textbf{j}\textbf{k}-\bm{\delta}\textbf{l}}-\Psi_{\textbf{i}\textbf{j}\textbf{k}\textbf{l}-\bm{\delta}}\right)+\left[\bm{\delta}\to-\bm{\delta}\right]\bigg{]}$
$\displaystyle-
t_{2}\sum_{\bm{\delta}\in\bm{\delta}^{(2)}}\bigg{[}\left(\Psi_{\textbf{i}+\bm{\delta}\textbf{j}\textbf{k}\textbf{l}}-\Psi_{\textbf{i}\textbf{j}-\bm{\delta}\textbf{k}\textbf{l}}-\Psi_{\textbf{i}\textbf{j}\textbf{k}-\bm{\delta}\textbf{l}}-\Psi_{\textbf{i}\textbf{j}\textbf{k}\textbf{l}-\bm{\delta}}\right)+\left[\bm{\delta}\to-\bm{\delta}\right]\bigg{]}$
$\displaystyle+U\sum_{\textbf{i},\textbf{j},\textbf{k},\textbf{l}}\left(1-\delta_{\textbf{i}\textbf{j}}-\delta_{\textbf{i}\textbf{k}}-\delta_{\textbf{i}\textbf{l}}\right)\Psi_{\textbf{ijkl}}=\left(\epsilon^{3h1e}(t_{1},t_{2})+2\mu+2h\right)\Psi_{\textbf{ijkl}}$
(B.5)

In the limit of infinite $U$, the double occupancy is forbidden. And thus, the
doublon must be located at the same site with one of the holes (spin flip
channel). In other words, the condition for one of the Kronecker deltas must
always be satisfied,

$\displaystyle\alpha_{i}=\alpha_{j},\quad\text{or}\quad\alpha_{i}=\alpha_{k},\quad\text{or}\quad\alpha_{i}=\alpha_{l}.$
(B.6)

This significantly reduces the Hilbert space dimension, allowing to simulate
larger system sizes.

We also note that a similar approach will yield equations for any number of
holes and doublons. The resulting equations can be easily recovered in the
general case by noticing that the doublons hop with a coefficient $-t_{1}$ and
the holes with $+t_{1}$.

###  B.2 Fixed Center of Mass Momentum

In the discussed problem the center of mass (COM) momentum serves as a good
quantum number. To exploit this, we construct the basis of Fock states in the
site occupation representation. Then construct a COM translation operator
$T_{\textbf{R}}$ which shifts all the particles in the system by R, and thus
preserves the relative coordinates (mod s). For every Fock state we apply
$T_{\textbf{R}}$ to see what states it is related to and sort all Fock states
into cycles. The eigenstates of the translation operator are then constructed
from each of these cycles as

$\displaystyle\Psi^{n}_{\textbf{P}}(\textbf{r}_{ij},\textbf{r}_{ik},\textbf{r}_{il})=\frac{1}{N_{n}}\sum_{\textbf{R}}e^{-i\textbf{P}\textbf{R}}T_{\textbf{R}}\Psi(0,\textbf{r}_{ij},\textbf{r}_{ik},\textbf{r}_{il}),$
(B.7)

where $N_{n}$ is the length of the cycle $n$. Each cycle has $N_{n}$ allowed
values for P which are determined from the periodicities under applying the
translation operator $T_{\textbf{R}}$.

We then decompose the wave function into a superposition of these states and
apply the Hamiltonian to construct the corresponding matrix and determine the
spectrum. This approach is applicable in general case but often is
computationally intensive.

In the case when we have a single doublon, we can employ a more optimized
version of this approach by fixing the COM coordinate with the doublon and
introducing the relative ones as,

$\displaystyle\textbf{R}=\textbf{i},\quad\textbf{r}_{ij}=\textbf{j}-\textbf{i},\quad\textbf{r}_{ik}=\textbf{k}-\textbf{i},\quad\textbf{r}_{il}=\textbf{l}-\textbf{i}.$
(B.8)

We carry out a Fourier transform over R,

$\displaystyle\Psi_{\textbf{P}}(\textbf{r}_{ij},\textbf{r}_{ik},\textbf{r}_{il})=\frac{1}{L}\sum_{\textbf{R}}e^{-i\textbf{P}\textbf{R}}\Psi(\textbf{R},\textbf{r}_{ij},\textbf{r}_{ik},\textbf{r}_{il}),$
(B.9)

to obtain equations for a fixed center of mass momentum P. They can be written
as,

$\displaystyle-
t_{1}\sum_{\bm{\delta}\in\bm{\delta}^{(1)}}\bigg{[}\left(e^{i\textbf{P}\bm{\delta}}\Psi_{\textbf{P}}\left(\textbf{r}_{ij}-\bm{\delta},\textbf{r}_{ik}-\bm{\delta},\textbf{r}_{il}-\bm{\delta}\right)-\Psi_{\textbf{P}}\left(\textbf{r}_{ij}+\bm{\delta},\textbf{r}_{ik},\textbf{r}_{il}\right)-\Psi_{\textbf{P}}\left(\textbf{r}_{ij},\textbf{r}_{ik}+\bm{\delta},\textbf{r}_{il}\right)-\Psi_{\textbf{P}}\left(\textbf{r}_{ij},\textbf{r}_{ik},\textbf{r}_{il}+\bm{\delta}\right)\right)+\text{h.c.}\bigg{]}$
$\displaystyle-
t_{2}\sum_{\bm{\delta}\in\bm{\delta}^{(2)}}\bigg{[}\left(e^{i\textbf{P}\bm{\delta}}\Psi_{\textbf{P}}\left(\textbf{r}_{ij}-\bm{\delta},\textbf{r}_{ik}-\bm{\delta},\textbf{r}_{il}-\bm{\delta}\right)-\Psi_{\textbf{P}}\left(\textbf{r}_{ij}+\bm{\delta},\textbf{r}_{ik},\textbf{r}_{il}\right)-\Psi_{\textbf{P}}\left(\textbf{r}_{ij},\textbf{r}_{ik}+\bm{\delta},\textbf{r}_{il}\right)-\Psi_{\textbf{P}}\left(\textbf{r}_{ij},\textbf{r}_{ik},\textbf{r}_{il}+\bm{\delta}\right)\right)+\text{h.c.}\bigg{]}$
$\displaystyle\quad+U\left(1-\delta_{\textbf{r}_{ij}0}-\delta_{\textbf{r}_{ik}0}-\delta_{\textbf{r}_{il}0}\right)\Psi_{\textbf{P}}\left(\textbf{r}_{ij},\textbf{r}_{ik},\textbf{r}_{il}\right)=\left(\epsilon^{3h1e}(t_{1},t_{2})+2\mu+2h\right)\Psi_{\textbf{P}}\left(\textbf{r}_{ij},\textbf{r}_{ik},\textbf{r}_{il}\right)$
(B.10)

We can again introduce the linear variables
$\alpha_{ij},\alpha_{ik},\alpha_{il}$. The condition to cover all the Hilbert
space is

$\displaystyle\alpha_{ij}<\alpha_{ik}<\alpha_{il}$ (B.11)

We note that because of this condition the second and third Kronecker deltas
in the above equation will always be 0.

The fixed CM momentum becomes even more useful when dealing with the infinite
$U$ limit. In this case, the constriction becomes even more strict, limiting
the Hilbert space to $\alpha_{ij}=1$.

###  B.3 Exact Diagonalization (ED)

To employ ED we construct the hopping matrices described above. Though their
linear size increases rapidly with the system size, the matrices remain very
sparse. This allows us to employ algorithms specifically designed for storing
(compressed sparse row CSR) and manipulating large sparse matrices. Namely,
the first several eigenvalues and eigenvectors are found by employing Lanczos
algorithm.

We implemented our code both in Wolfram Mathematica and Python to assure an
exact agreement of the numerical results. Additionally, we developed two
different methods for ED in real space, a general code which deals with any
given number of holes and doublons, as well as specific algorithm which works
more efficiently for the 3 hole 1 doublon problem, allowing to scale the
system size. They both agreed exactly on every simulated scenario.

The both algorithms were implemented for finite and infinite $U$ cases. We
made sure that the results in the large $U$ limit converged to the infinite
$U$ results seamlessly. For infinite $U$ limit, we were able to scale the
system size up to $N=18\times 18$.

Separately, we implemented the ED code in the basis of fixed CM momentum as
well. Thanks to the reduced Hilbert space dimension, we were able to reach a
rather large system size of $N=36\times 36$ for the 2H1S problem. We checked
that the results for smaller system sizes agree exactly with the direct ED
calculations.

##  C Robustness of results

In this section we examine the 2H1M phase, discussed in the main text, and
show that it is a robust property of the our model, resilient within infinite
lattice systems. To conclusively establish this we explore the energy gap
observed in the Fig. 2a),b) of the main text. We analyze the width of the gap
as a function of the system size and extend of our findings from finite
lattices to infinite lattice configurations. The extrapolated data
corresponding to finite systems is graphically represented in Fig. S2a). For
small enough $t_{2}/t_{1}$, when the band minimum is at $K$ point, the
corresponding gap shrinks as $\Delta_{K}\propto L^{-2}$, eventually vanishing
for an infinite system. Conversely, the gap width demonstrates a different
pattern for larger $t_{2}$, when the band minima shift to $M$ or $\Gamma$
points. While here the gap again follows a $L^{-2}$ scaling with a reasonable
accuracy, it saturates to a finite value in the infinite size limit,

$\displaystyle\Delta_{M}\approx 0.146t_{1},\quad t_{2}=0.09t_{1},$ (C.1)
$\displaystyle\Delta_{\Gamma}\approx 0.375t_{1},\quad t_{2}=0.15t_{1}.$ (C.2)

This delivers compelling evidence that this energy gain is a phenomenon
inherent to infinitely large systems as well, thus affirming the robustness of
the predicted phase as an intrinsic feature of the Hubbard model.

Additionally, we illustrate the dependency of the average distance between the
holes on the parameter $t_{2}/t_{1}$ in the Fig. S2b), and compare the results
for systems of increasing linear size $L$. we observe that, initially, at
small $t_{2}$, the holes are far apart (of order of system size), indicating
that one of the holes bounds to the spin flip, forming a polaron, while the
other hole is repelled from that pair. Here the average size scales with the
system size.

When $t_{2}/t_{1}$ exceeds $0.09$, the localization length of the ground state
abruptly drops, reaching $|\Delta r|/a\sim 3$. This occurs because the GS
transitions to the $P=M$ sector, which was already localized even for smaller
$t_{2}$. Importantly, here the pair size shows very little sensitivity to the
system size, proving the state to be well-localized.

Figure S2:  (a) Explores the energy gaps at fixed CM momentum as a function of
inverse linear system size $1/L$, when the corresponding state is the GS. The
extrapolation reveal that, for $t_{2}=0$ the gap vanishes for infinite
systems, while for $t_{2}/t_{1}=0.09,0.15$ saturates to finite values. (b) The
average distance between the holes as a function of $t_{2}/t_{1}$ for
increasing system size. The state becomes well-localized for
$t_{2}/t_{1}\gtrsim 0.09$, and the GS momentum shifts from $K$ to $M$, then t
$\Gamma$. (c) Effective masses of the 2H1M particle $m_{x},m_{y}$. The masses
are anisotropic when the GS is at $M$ points.

##  D Finite on-site interaction $U$

Our main text focused on the limit of infinite on-site interaction energy
$U/t_{1}$. As discussed above, the Hilbert space dimension is significantly
reduced in that case, allowing for more reliable calculations for larger
system sizes. In this section, we present the results for finite $U/t_{1}$
case. From our main discussion we know that the $K$ and $M$ points in momentum
space play a pivotal role in our system. Thus when fixing the finite system
size $L$, we need to choose it in a way that both $K$ and $M$ points are
reachable, i.e. $L$ should be divisible by $6$. And since $12\times 12$ is
very computationally intensive, the largest system size that we can have is
$N=6\times 6$.

We illustrate that the remarkable property of 2H1M bound state manifests in
the case of finite interaction $U$ as well. In the Fig. S3 we illustrate the
energy gap between the 2H1M state and the (1H1M + 1H) state (polaron and
hole):

$\displaystyle\Delta=\epsilon^{3h1e}-(\epsilon^{pol}+\epsilon^{h}).$ (D.3)

We observe that for any interaction $U$, the localized magnonic Cooper pair
state becomes more favorable for large enough second-neighbor hopping
parameter. Interestingly, while for large enough $U$ the localization takes
place near $t_{2}/t_{1}\approx 0.09$, in agreement with the infinite $U$
results, by lowering the $U$, the threshold value is gradually pushed towards
exactly $t_{2}=t_{1}/8$. This result is in a perfect agreement with the
intuitive picture of the pairing mechanism involving the holes shifting to the
$M$ points, described in the main text. Additionally, we see that the gap has
a local maximum for intermediate interactions of $U/t_{1}\sim 15$.

Figure S3: Energy gap between the 2H1M and the polaron and hole (1H1M + 1H)
state, Eq. (D.3) as a function of $U/t_{1}$ and $t_{2}/t_{1}$.

##  E Interactions between 2H1M particles

As discussed in the main text, we determined the sign of interactions of the
magnonic Cooper pairs by performing an ED study for a system of 4 holes and 2
magnons on a $6\times 6$ lattice. We showed that the net interaction is
attractive for these particles.

Here we present the analysis for the effective model described in the main
text for studying the interaction of effective bosons. To illustrate of the
localization of the bosons into a composite particles in the low doping
regime, we consider a triangular lattice model with nearest neighbor hopping.
In the long-wavelength regime the internal structure of the lattice is
rendered unimportant, but the choice of the triangular lattice facilitates the
modeling of the anisotropic spectrum of the 2H1M bosons. We obtained that
these bosons have anisotropic masses, typically with $m_{y}/m_{x}\sim 2$. We
orient eigen-axis for the type 1 boson dispersion, $\alpha=1$, along the $x$
and $y$ axis with masses $m_{x}$ and $m_{y}$. For the boson 2 and 3 we can
then rotate these axes by $\pi/3$. In other words, for the boson of type 1, as
an example, we model hopping with parameter $t$ in the $x$ direction, and
$t^{\prime}$ in the two other directions. In the long wavelength limit, the
spectrum is expanded as

$\displaystyle\epsilon_{b,1}=\frac{p_{x}^{2}}{2m_{x}}+\frac{p_{y}^{2}}{2m_{y}},\quad
m_{x}=\frac{\hbar^{2}}{a^{2}(2t+t^{\prime})},\quad
m_{y}=\frac{\hbar^{2}}{3a^{2}t^{\prime}}$ (E.4)

To obtain $m_{y}/m_{x}\sim 2$ for 2H1M boson, we choose $t^{\prime}/t=1/2$.
The anisotropic spectra for the other two type bosons are modeled in a similar
manner.

In order to show the localization, it is enough to perform ED analysis with at
most a single boson of each type. Thus, in the Hamiltonian from the Eq.(3) of
the main text, the $U,J$ terms are left out as they involve interactions of
two bosons of the same type. We model the attractive $V$ term on the lattice
as on-site interaction with magnitude $V_{0}$.

###  E.1 Two boson localization

We start from the two boson localization on the example of $\alpha=1,2$ type
bosons. Similarly to the analysis in the Sec. B.1, we can write the Shrödinger
equation for the two particle wave function $\Phi_{\textbf{ij}}$ in the form

$\displaystyle-t\left(\Phi_{\left(\textbf{i}+\bm{\delta}_{1}\right)\textbf{j}}+\Phi_{\left(\textbf{i}-\bm{\delta}_{1}\right)\textbf{j}}\right)-t^{\prime}\sum_{m=2,3}\left(\Phi_{\left(\textbf{i}+\bm{\delta}_{m}\right)\textbf{j}}+\Phi_{\left(\textbf{i}-\bm{\delta}_{m}\right)\textbf{j}}\right)$
(E.5)
$\displaystyle\quad-t\left(\Phi_{\textbf{i}\left(\textbf{j}+\bm{\delta}_{2}\right)}+\Phi_{\textbf{i}\left(\textbf{j}-\bm{\delta}_{2}\right)}\right)-t^{\prime}\sum_{m=1,3}\left(\Phi_{\textbf{i}\left(\textbf{j}+\bm{\delta}_{m}\right)}+\Phi_{\textbf{i}\left(\textbf{j}-\bm{\delta}_{m}\right)}\right)-V_{0}\Phi_{\textbf{ij}}\delta_{\textbf{ij}}=E^{(2)}\Phi_{\textbf{ij}}$
(E.6)

The two-particle problem can be conveniently studied by changing the initial
variables to the center of mass R and relative motion l,

$\displaystyle\Phi_{\textbf{ij}}=\Phi_{\textbf{Rl}},\quad\textbf{R}=\textbf{i},\quad\textbf{l}=\textbf{j}-\textbf{i}$
(E.7)

This will rewrite the equation above in the form

$\displaystyle-t\left(\Phi_{\left(\textbf{R}+\bm{\delta}_{1}\right)\left(\textbf{l}-\bm{\delta}_{1}\right)}+\Phi_{\left(\textbf{R}-\bm{\delta}_{1}\right)\left(\textbf{l}+\bm{\delta}_{1}\right)}\right)-t^{\prime}\sum_{m=2,3}\left(\Phi_{\left(\textbf{R}+\bm{\delta}_{m}\right)\left(\textbf{l}-\bm{\delta}_{m}\right)}+\Phi_{\left(\textbf{R}-\bm{\delta}_{m}\right)\left(\textbf{l}+\bm{\delta}_{m}\right)}\right)$
$\displaystyle\quad-t\left(\Phi_{\textbf{R}\left(\textbf{l}+\bm{\delta}_{2}\right)}+\Phi_{\textbf{R}\left(\textbf{l}-\bm{\delta}_{2}\right)}\right)-t^{\prime}\sum_{m=1,3}\left(\Phi_{\textbf{R}\left(\textbf{l}+\bm{\delta}_{m}\right)}+\Phi_{\textbf{R}\left(\textbf{l}-\bm{\delta}_{m}\right)}\right)-V_{0}\Phi_{\textbf{Rl}}\delta_{\textbf{l}}=E^{(2)}\Phi_{\textbf{Rl}}$
(E.8)

We can then carry out a Fourier transform over R,
$\Phi_{\textbf{P}}(\textbf{l})=\frac{1}{\sqrt{N}}\sum_{\textbf{R}}e^{-i\textbf{P}\textbf{R}}\Phi_{\textbf{R}\textbf{l}}$,
to obtain equations for a fixed center of mass momentum P. The resulting
equations read as follows:

$\displaystyle-\left(te^{i\textbf{P}\bm{\delta}_{1}}+t^{\prime}\right)\Phi_{\textbf{P}}\left(\textbf{l}-\bm{\delta}_{1}\right)-\left(te^{-i\textbf{P}\bm{\delta}_{1}}+t^{\prime}\right)\Phi_{\textbf{P}}\left(\textbf{l}+\bm{\delta}_{1}\right)$
$\displaystyle\quad-\left(t^{\prime}e^{i\textbf{P}\bm{\delta}_{2}}+t\right)\Phi_{\textbf{P}}\left(\textbf{l}-\bm{\delta}_{2}\right)-\left(t^{\prime}e^{-i\textbf{P}\bm{\delta}_{2}}+t\right)\Phi_{\textbf{P}}\left(\textbf{l}+\bm{\delta}_{2}\right)$
$\displaystyle\quad-t^{\prime}\left(e^{i\textbf{P}\bm{\delta}_{3}}+1\right)\Phi_{\textbf{P}}\left(\textbf{l}-\bm{\delta}_{3}\right)-t^{\prime}\left(e^{-i\textbf{P}\bm{\delta}_{3}}+1\right)\Phi_{\textbf{P}}\left(\textbf{l}+\bm{\delta}_{3}\right)-V_{0}\Phi_{\textbf{P}}\left(\textbf{l}\right)\delta_{\textbf{l}}=E\Phi_{\textbf{P}}({\textbf{l}})$
(E.9)

This equation has two different families of solution, depending on whether
$\Phi_{\textbf{P}}(\textbf{l}=0)$ is zero or finite. In the first case, there
is no overlap between the wave functions of the two bosons, while in the
second case the bosons are bound. Here we will analyze them and quantify the
transition from one to another depending on the parameters of the model.

(1) For the first case, $\Phi_{\textbf{P}}(\textbf{l}=0)=0$, we can proceed
easily by carrying out another Fourier transform over the relative coordinate
l, as
$\tilde{\Phi}_{\textbf{P}}(\textbf{k})=\frac{1}{\sqrt{N}}\sum_{\textbf{l}}e^{-i\textbf{k}\textbf{l}}\Phi_{\textbf{P}}(\textbf{l})$.
Then the equations are algebraic and can be solved to yield the dispersion for
this case:

$\displaystyle E^{(2)}_{\text{unb}}(\textbf{P},\textbf{k})$
$\displaystyle=\left(te^{i\textbf{P}\bm{\delta}_{1}}+t^{\prime}\right)e^{-i\textbf{k}\bm{\delta}_{1}}+\left(te^{-i\textbf{P}\bm{\delta}_{1}}+t^{\prime}\right)e^{i\textbf{k}\bm{\delta}_{1}}$
$\displaystyle+\left(t^{\prime}e^{i\textbf{P}\bm{\delta}_{2}}+t\right)e^{-i\textbf{k}\bm{\delta}_{2}}+\left(t^{\prime}e^{-i\textbf{P}\bm{\delta}_{2}}+t\right)e^{i\textbf{k}\bm{\delta}_{2}}$
$\displaystyle+t^{\prime}\left(e^{i\textbf{P}\bm{\delta}_{3}}+1\right)e^{-i\textbf{k}\bm{\delta}_{3}}+t^{\prime}\left(e^{-i\textbf{P}\bm{\delta}_{3}}+1\right)e^{i\textbf{k}\bm{\delta}_{3}}$
(E.10)

In this case the bosons move absolutely independently. It is easy to show that
the minimum energy is achieved for COM momentum $\textbf{P}=0$, and the long-
wavelength expansion for the dispersion reads as

$\displaystyle
E^{(2)}_{\text{unb}}(\textbf{k})\approx-4(t+2t^{\prime})+\frac{\hbar^{2}q_{x}^{2}}{2M_{1}}+\frac{\hbar^{2}q_{y}^{2}}{2M_{2}},\quad
M_{1}=\frac{\hbar^{2}}{3a^{2}(t+t^{\prime})};\quad
M_{2}=\frac{\hbar^{2}}{a^{2}(t+5t^{\prime})},$ (E.11)
$\displaystyle\begin{pmatrix}q_{x}\\\
q_{y}\end{pmatrix}=\begin{pmatrix}\cos\phi&-\sin\phi\\\
\sin\phi&\cos\phi\end{pmatrix}\begin{pmatrix}k_{x}\\\
k_{y}\end{pmatrix},\quad\phi=-\pi/6.$

With the minimal energy being $E^{(2)}_{\text{unb,min}}=-4(t+2t^{\prime})$. We
note that this is simply twice the GS energy of a single boson
$E^{(1)}_{\text{min}}=-2(t+2t^{\prime})$.

(2) In the case when $\Phi_{\textbf{P}}(\textbf{l}=0)\neq 0$, we again carry
out a Fourier transform. Here we have an additional term arising from
$\delta_{\textbf{l}0}$. It can be Fourier transformed as,
$\frac{1}{\sqrt{N}}\sum_{\textbf{l}}\delta_{\textbf{l}0}f(\textbf{l})e^{-i\textbf{k}\textbf{l}}=\frac{1}{\sqrt{N}}f(\textbf{l}=0)$.
The resulting equation is then written as follows,

$\displaystyle
E^{(2)}_{\text{b}}\tilde{\Phi}_{\textbf{P}}(\textbf{k})=E^{(2)}_{\text{unb}}(\textbf{P},\textbf{k})\tilde{\Phi}_{\textbf{P}}(\textbf{k})-\frac{V_{0}}{\sqrt{N}}\Phi_{\textbf{P}}(\textbf{l}=0)$
(A3.12)

We then express,

$\displaystyle\tilde{\Phi}_{\textbf{P}}(\textbf{k})=-\frac{V_{0}}{\sqrt{N}}\frac{\Phi_{\textbf{P}}(\textbf{l}=0)}{E^{(2)}_{\text{b}}-E^{(2)}_{\text{unb}}(\textbf{P},\textbf{k})}$
(E.13)

Carrying out an inverse Fourier transform we obtain,

$\displaystyle\Phi_{\textbf{P}}(\textbf{l})=\frac{1}{\sqrt{N}}\sum_{\textbf{k}}e^{i\textbf{k}\textbf{l}}\tilde{\Phi}_{\textbf{P}}(\textbf{k})=-\frac{V_{0}}{N}\Phi_{\textbf{P}}(\textbf{l}=0)\sum_{\textbf{k}}\frac{e^{i\textbf{k}\textbf{l}}}{E^{(2)}_{\text{b}}-E^{(2)}_{\text{unb}}(\textbf{P},\textbf{k})}$
(E.14)

Now we fix $\textbf{l}=0$ to obtain a self-consistency equation on the energy
$E^{(2)}_{\text{b}}$,

$\displaystyle
1=-\frac{V_{0}}{N}\sum_{\textbf{k}}\frac{1}{E^{(2)}_{\text{b}}-E^{(2)}_{\text{unb}}(\textbf{P},\textbf{k})}$
(E.15)

Depending on the parameters of the system, either of the two above solutions
can have lower energy. Let us find the critical value of the on-site
interaction term $V^{\text{c}}_{0}$, after which the bound solution is more
favorable than the unbound one.

The lowest energy from the unbound solution is, as we found,
$E^{(2)}_{\text{unb, min}}=-4(t+2t^{\prime})$. At $V_{0}=V^{\text{c}}_{0}$, it
should coincide with the energy of bound solution. So, we plug it into the Eq.
(E.15), evaluate the sum and express $V^{\text{c}}_{0}$.

The summation can be replaced by integration as
$\sum_{\textbf{k}}\to\frac{A}{(2\pi)^{2}}\int_{BZ}d^{2}\textbf{k}$, where $A$
is the sample surface area. Additionally, since the second term in the
denominator Eq. has a minimum at $\textbf{k}=0$, we expand the denominator up
to the second order to arrive to the following integral,

$\displaystyle
1\approx\frac{V^{\text{c}}_{0}}{N}\frac{A}{(2\pi)^{2}}\int_{BZ}d^{2}\textbf{q}\frac{1}{\frac{\hbar^{2}q_{x}^{2}}{2M_{1}}+\frac{\hbar^{2}q_{y}^{2}}{2M_{2}}}$
(E.16)

We notice that the last integral diverge logarithmically at small $q$. To
avoid the divergence, we introduce a cutoff at the inverse lattice size $q\sim
1/(La)$. Within the logarithmic precision, the integral is evaluated as,

$\displaystyle
1\approx\frac{V^{\text{c}}_{0}}{N}\frac{A}{(2\pi)^{2}}\int_{1/Na}^{1/a}\frac{dq}{q}\int_{0}^{2\pi}\frac{d\varphi}{\frac{\hbar^{2}\cos^{2}\varphi}{2M_{1}}+\frac{\hbar^{2}\sin^{2}\varphi}{2M_{2}}}=\frac{V^{\text{c}}_{0}}{N}\frac{A}{(2\pi)^{2}}\frac{2\pi\ln
L}{3ta^{2}}J(t^{\prime}/t)$ (E.17)

Here we introduced a dimensionless function

$\displaystyle
J(t^{\prime}/t)=\int_{0}^{2\pi}\frac{d\varphi}{2\pi}\frac{6}{3(1+t^{\prime}/t)\cos^{2}\varphi+(1+5t^{\prime}/t)\sin^{2}\varphi}=\frac{2\sqrt{3}}{\sqrt{5(t^{\prime}/t)^{2}+6t^{\prime}/t+1}},\quad
J(1)=1$ (E.18)

The factor $A/N$ is the surface area of a unit cell of a triangular lattice,
$s_{0}=a^{2}\sqrt{3}/2$. This finally gives,

$\displaystyle\frac{V^{\text{c}}_{0}}{t}\approx\frac{2\pi}{\ln
L}\sqrt{5(t^{\prime}/t)^{2}+6t^{\prime}/t+1}\propto\frac{1}{\ln L}.$ (E.19)

Importantly, this result shows that the critical value of the on-site term,
$V^{\text{c}}_{0}$, tends to zero with the increasing lattice size. This
signals that for an infinite system, the two bosons will form a bound state
even at infinitesimal interaction strength.

We illustrate the two boson localization in the Fig.S4 by solving the two-body
problem numerically and looking into the average relative distance between the
holes as a function of interaction strength $V_{0}$.

Figure S4: Mean relative distance $\langle r\rangle/a$ between two bosons with
$t^{\prime}/t=1/2$ for $L=30,12$ lattice sizes versus interaction strength
$V_{0}/t$. For weak $V_{0}$ the relative distance is of order of the system
size. With increasing $V_{0}$, the state tends to become localized.

###  E.2 Three boson localization

A key property that our models allows to study is the three-boson bound state.
We perform an ED analysis with three bosons of different types. We start by
looking into the average distance between the bosons in the GS. The
calculations reveal that the relative distances between the bosons are equal
and they decrease with increasing interaction $V_{0}$. The results are
depicted in Fig. S5a) for two different system sizes, indicating the emergence
of a three-boson bound state.

In order to determine which type of localization will take place in the case
of finite doping, the two-boson or three-boson binding, we look into the
binding energies per bosons in both cases. Let us first introduce the binding
energies as

$\displaystyle\Delta^{(2)}(V_{0})=E^{(2)}(V_{0})-2E^{(1)},\quad\Delta^{(3)}(V_{0})=E^{(3)}(V_{0})-3E^{(1)},$
(E.20)

for two and three particle states, respectively, and
$E^{(1)}=-2(t+2t^{\prime})$ is the single boson GS energy. And the binding
energies per boson are then naturally defined as

$\displaystyle\Delta^{(2)}_{b}=\Delta^{(2)}/2,\quad\Delta^{(3)}_{b}=\Delta^{(3)}/3.$
(E.21)

We compare these values in the Fig. S5b). The results indicate that the
binding per boson in the three-boson bound state is always stronger than in
the two-boson bound state. This establishes, that the bosonic system in case
of finite doping will condense into bosonic trios. Therefore, the 2H1M bosons
will bound into 6H3M states with charge 6e and spin $S=6$.

Figure S5: (a) Mean relative distance $\langle r\rangle/a$ between three
bosons with $t^{\prime}/t=1/2$ for $L=18,12$ lattice linear sizes versus
interaction strength $V_{0}/t$. For weak $V_{0}$ the relative distance is of
order of the system size. With increasing $V_{0}$, the state tends to become
localized. (b) The comparison of binding energies per boson for three-boson
bound state $|\Delta^{(3)}_{b}|/t$ and two-boson bound state
$|\Delta^{(2)}_{b}|/t$. The three-boson bound state is more favorable

##  F Hardcore boson model for 2H1M particles

In order to demonstate that the three-boson localization is a robust property,
qualitatively independent on the model, here we present the analysis of
another toy model leading to three 2H1M boson localization into a composite
particles. We consider hard core bosons (forbidden double occupancy) with
nearest neighbor attraction $V_{0}$. To obtain a lattice model with band
minima at the $M$ points, we take into consideration a triangular lattice with
both first and second nearest neighbor hoppings, similar to our main model.
This model is written as

$\displaystyle
H=t_{1}^{\prime}\sum_{\langle\textbf{i},\textbf{j}\rangle}\left(b_{\textbf{i}}^{\dagger}b_{\textbf{j}}+h.c.\right)+t_{2}^{\prime}\sum_{\langle\langle\textbf{i},\textbf{j}\rangle\rangle}\left(b_{\textbf{i}}^{\dagger}b_{\textbf{j}}+h.c.\right)-V_{0}\sum_{\langle\textbf{i},\textbf{j}\rangle}n_{\textbf{i}}n_{\textbf{j}},\quad
t_{2}^{\prime}\geq 0.125t_{1}^{\prime}$ (F.22)

In this model, the effective masses of the particles near the band minima are
anisotropic. We can determine them by expanding the dispersion around the
$M_{1}$ point

$\displaystyle
m_{x}(M_{1})=\frac{\hbar^{2}}{a^{2}(9t^{\prime}_{2}-t^{\prime}_{1})};\quad
m_{y}(M_{1})=\frac{\hbar^{2}}{3a^{2}(t^{\prime}_{1}-t^{\prime}_{2})};$ (F.23)

This allows to tune the hopping parameters to obtain the actual effective
masses of the 2H1S bosons.

We then dope three bosons into the system and explore their interactions on a
finite lattice of linear size $L$. In the absence of $V_{0}$, the bosons are
delocalized. The Fig. S6 shows the localization of the bosons as a function of
interaction strength $V_{0}$. We compare the average distance between the
bosons for different system sizes of $L=12,18$. For $V_{0}\geq 0.3$, the both
results agree with $20\%$ accuracy, providing a compelling evidence of three-
particle localization. Finally, we look into the binding energy $\Delta$ of
the bound state by comparing the energy of the composite boson to thrice that
of a single boson at $M$ point. The inset to the Fig. S6 illustrates that
while $\Delta$ is slightly positive for vanishing $V_{0}$, it already becomes
negative for $V_{0}/t_{1}^{\prime}\sim 0.1$, indicating attraction.

Figure S6: Localization of hard-core bosons for
$t^{\prime}_{2}=0.15t^{\prime}_{1}$. The figure illustrates the average
distance between the bosons in the ground state as a function of $V_{0}$ for
$L=12,18$. The inset shows the binding energy of the bound state.