\section*{}}

\begin{document}
\normalem

\title{Pairing of holes by confining strings in antiferromagnets}

\author{F. Grusdt}
\email[Corresponding author email: ]{fabian.grusdt@lmu.de}
\affiliation{Department of Physics and Arnold Sommerfeld Center for Theoretical Physics (ASC), Ludwig-Maximilians-Universit\"at M\"unchen, Theresienstr. 37, M\"unchen D-80333, Germany}
\affiliation{Munich Center for Quantum Science and Technology (MCQST), Schellingstr. 4, D-80799 M\"unchen, Germany}

\author{E. Demler}
\affiliation{Institut f\"{u}r Theoretische Physik, ETH Zurich, 8093 Zurich, Switzerland}

\author{A. Bohrdt}
\address{ITAMP, Harvard-Smithsonian Center for Astrophysics, Cambridge, MA 02138, USA}
\affiliation{Department of Physics, Harvard University, Cambridge, Massachusetts 02138, USA}

\pacs{}

\date{\today}

\begin{abstract}
In strongly correlated quantum materials, the behavior of charge carriers is dominated by strong electron-electron interactions. These can lead to insulating states with spin order, and upon doping to competing ordered states including unconventional superconductivity. The underlying pairing mechanism remains poorly understood however, even in strongly simplified theoretical models. Recent advances in quantum simulation allow to study pairing in paradigmatic settings, e.g. in the $t-J$ and $t-J_z$ Hamiltonians. Even there, the most basic properties of paired states of only two dopants, such as their dispersion relation and excitation spectra, remain poorly studied in many cases. Here we provide new analytical insights into a possible string-based pairing mechanism of mobile holes in an antiferromagnet. We analyze an effective model of partons connected by a confining string and calculate the spectral properties of bound states. Our model is equally relevant for understanding Hubbard-Mott excitons consisting of a bound doublon-hole pair or confined states of dynamical matter in lattice gauge theories, which motivates our study of different parton statistics. Although an accurate semi-analytic estimation of binding energies is challenging, our theory provides a detailed understanding of the internal structure of pairs. For example, in a range of settings we predict heavy states of immobile pairs with flat-band dispersions -- including for the lowest-energy $d$-wave pair of fermions. Our findings shed new light on the long-standing question about the origin of pairing and competing orders in high-temperature superconductors.
\end{abstract}

\maketitle

\section{Introduction}
The history of physics has repeatedly taught us that nature tends to realize richer structures than one might first suggest. The most important example, by far, is constituted by the theory of atoms which has evolved from Thomson's featureless plum pudding model to our current picture of precisely quantized energy shells and a nucleus with structures down to the level of individual quarks. The only way to reveal such structures is to perform increasingly more precise measurements, or, with the benefit of hindsight and a microscopic Hamiltonian at hand, perform increasingly more precise numerical simulations. Here we address the question how much, and which, structure paired charge carriers in correlated quantum matter may have.

To date, important aspects of strongly correlated electrons remain poorly understood. Part of the reason is that the nature and microscopic structure of the emergent charge carriers is not fully understood. Among the most famous puzzles is the origin of pairing in high-temperature superconductors \cite{Lee2006,Keimer2015}, but similar questions arise in heavy-fermion compounds \cite{Stewart1984}, organic superconductors \cite{Lebed2008} and most recently twisted bilayer graphene \cite{Andrei2020}. However the problems are not limited to paired states of matter: the nature of charge carriers in the exotic normal phases of these strongly correlated systems is likewise debated and subject of active research.

Remarkably, the prevailing picture of quantum matter is one with more-or-less featureless charge carriers, with little or no rigid internal structure taken into account in calculations. Some theoretical approaches assume fractionalization of quantum numbers, which leads to rich and interesting physics, but short-range spatial fluctuations are rarely considered in detail. Among the reasons for this restriction is the difficulty to directly experimentally visualize spatial structures of quickly fluctuating charges. While ultracold atoms in optical lattices \cite{Tarruell2018,Bohrdt2021PWA} have taken very promising steps in this direction \cite{Chiu2019Science,Koepsell2019}, they too have not managed to fully visualize the internal structure of doped charge carriers yet. Another reason may be the lasting influence of Anderson's RVB theory of high-$T_c$ superconductivity \cite{Anderson1987}, which assumes point-like charge carriers moving in a surrounding spin-liquid -- in some sense the antithesis of any theory assuming spatial structures of charge carriers.

On the other hand, it is not for a lack of ideas which kinds of internal structures could emerge: Even before the discovery of high-$T_c$ superconductors, Bulaevksii et al. proposed the existence of string-like states with internal vibrational excitations \cite{Bulaevskii1968}, a view taken up by Brinkman and Rice to understand dynamical properties of charge carriers \cite{Brinkman1970}; in 1988, Trugman applied this idea to pairs of holes \cite{Trugman1988} and in the same year Shraiman and Siggia analyzed two-hole string states in greater detail \cite{Shraiman1988a}. Their conclusions at the time were mixed: while they found mechanisms supporting pairing, they also identified unfavorable effects such as frustration of the pair-kinetic energy due to the underlying Fermi-statistics of the holes. Today these works stand out as being among the few and first attempting a description starting from the strong coupling antiferromagnetic (AFM) Mott limit (large Hubbard-$U$). But it appears that the community then focused more on approaches inspired by the weak-coupling (small Hubbard-$U$) limit of the theoretical models \cite{Keimer2015}, which more naturally led to the magnon-exchange picture \cite{Dahm2009}. There, magnetic fluctuations provide the glue between point-like charge carriers, and the theoretical framework shares more similarities with the successful BCS theory of conventional superconductors. As almost all ideas in the field, this picture has been debated \cite{Anderson2007glue}.

Nevertheless, over time the idea that charge carriers have a pronounced spatial structure was reclaimed several times. In 1996, Laughlin and co-workers proposed a phenomenological parton theory of doped holes, including a confining linear string tension \cite{Beran1996,Laughlin1997}; different kinds of spatial strings, termed phase-strings, were introduced in 1996 by Weng and co-workers \cite{Weng1997,Sheng1996}, and their effect on pairing was recently analyzed \cite{Zhu2014}; signatures for the more traditional $S^z$-string fluctuations were reported in large-scale DMRG simulations by White and Affleck in 2001 \cite{White2001}; in 2007 Manousakis proposed a string-based interpretation of one-hole ARPES spectra and in the same work envisioned a pairing mechanism of holes constantly exchanging fluctuating strings \cite{Manousakis2007}; in 2013 exact numerical simulations  by Vidmar et al. in truncated bases, closely related to the string picture, have also revealed signatures for pairs with a rich spatial structure \cite{Vidmar2013}; already in 2000 and 2001 large-scale Monte Carlo simulations by Brunner et al.~\cite{Brunner2000} and Mishchenko et al.~\cite{Mishchenko2001} have revealed long-lived vibrational excitations of individual doped holes; and in the past few years, the present authors have added new evidence for the existence of long-lived rotational and vibrational string states of individual charge carriers \cite{Grusdt2018PRX,Bohrdt2020,Bohrdt2021PRL}; Vibrational peaks have also long been known to exist and recently further confirmed in linear-spin wave models of doped holes \cite{Kane1989,Liu1992,Diamantis2021,Nielsen2022arXiv}. Finally, Hubbard-Mott excitons formed by a bound pair of a doublon and a hole \cite{Alpichshev2017} have been proposed to have a rich internal structure \cite{Huang2020arXiv,Huang2022arXiv}.

In this article, we revisit the idea that mobile dopants, the charge carriers in doped AFM Mott insulators, can form bound states with a rich spatial structure. Specifically, we derive a semi-analytical theory of pairs of holes connected by a confining string which fluctuates only through the motion of the charges at its ends. Our approach is very similar in spirit to the much earlier work by Shraiman and Siggia mentioned above \cite{Shraiman1988a}; in fact we confirm several of their predictions and discuss them in the context of three decades worth of new results including from far advanced numerics. This includes one of their most exciting -- though often overlooked -- prediction that two fermionic holes can form infinitely heavy pairs with a flat-band dispersion at very low energies. In fact, we find that that these flat bands have $d$-wave character. This result sheds new light on the wealth of competing ordered states observed in cuprates and numerically found in the closely related $t-J$ and Fermi-Hubbard models \cite{Corboz2014,Huang2017,Qin2019PRX}.

Back-to-back with this article, we are publishing a separate work focusing on a numerical analysis of rotational two-hole spectra in the $t-J_z$ and $t-J$ models \cite{Bohrdt2022pairingPrep}. There we achieve full momentum resolution on extended four-leg cylinders, and compare our numerical results to the semi-analytic calculations performed in the present article. The focus of the present article is on the semi-analytical method itself, including its formal derivation. Moreover, the calculations we perform here are applicable in a larger class of models: as explained below, our main assumption is that two partons on a square lattice are connected by a rigid string $\Sigma$ which creates a memory of the parton's motion in the wavefunction. In the case of mobile holes doped into an AFM Mott insulator, the string directly encodes the prevalent spin-charge correlations. But similar situations can be found in lattice gauge theories in a strongly confining regime where the dynamics of the corresponding electric field strings is dominated by fast charge fluctuations.

The goal of our present work is primarily to understand the properties of pairs of mobile dopants, i.e. their spatial structure and energy spectrum, including their rotational quantum numbers and effective mass. Previously much emphasis has been on the question whether the dopants pair up; i.e. about the magnitude and sign of the binding energy
\begin{equation}
E_{\rm bdg} = 2 (E_1 - E_0 ) - (E_2 - E_0)
\label{eqEbdgDef}
\end{equation}
where $E_n$ is the ground state energy in the presence of $n$ dopants. While we also view this as an important issue, we argue that it is often not a well-suited question for a semi-analytical approach, since it strongly depends on details; e.g. addressing it requires precise knowledge of the one-hole energy. In this article we take the view that the structure of the paired state can be very different from the structure of one-dopant states. Our goal is to understand the former, and we leave the question of how and when excited (or ground) states of pairs can decay into individual single-dopant states to future analysis (only a brief discussion within our model will be provided). Moreover, we note that the question of pairing is not identical to a question about the existence superconductivity: instead of condensing into a superconductor, pairs of holes may also crystalize and form a pair-density wave at finite doping \cite{Hamidian2016}.

Nevertheless, we will address the origins of pairing within our semi-analytical framework. To this end we identify competing effects which tend to increase or decrease the binding energy. Our results shed new light on earlier predictions \cite{Shraiman1988a,Trugman1988}: (i) fermionic statistics of the dopants frustrates the pair's kinetic energy, which is unfavorable for pairing; (ii) the most detrimental effect for pairing comes from the hard-core property of two dopants, which cannot occupy the same site; Moreover, we reveal (iii) a geometric spinon-chargon repulsion in dimensions $d \geq 2$ which enhances the one-dopant energy $E_1$ \cite{Grusdt2019PRB} and thus favors pairing. In addition to all of these, further contributions stemming from low-energy spin-fluctuations in the background are expected, whose quantitative effects are more challenging to predict and will thus be left to future work to explore. Finally, we note that in specifically tailored settings our simplifying assumptions become quantitatively accurate: In Ref.~\cite{Bohrdt2021pairing} we proposed a mixed-dimensional bilayer model with strong rung singlets where we demonstrated strong string-based pairing of doped holes with a binding energy scaling as $E_{\rm bdg} \simeq t^{1/3} J^{2/3}$, when the hole tunneling $t \gg J$ exceeds the rung super-exchange $J$. Recently, in a closely related mixed-dimensional two-leg ladder, ultracold fermions directly observed strong hole binding \cite{Hirthe2022arXiv}, realizing a decades old toy model of pairing \cite{Schulz1999,Troyer1996}. The internal structure of hole pairs in these models can also be described by the theoretical model we develop here.

This article is organized as follows. In Sec.~\ref{secMicModel} we briefly discuss the microscopic models motivating our analysis. Sec.~\ref{secEffStringModel} constitutes the main body of our article: there we develop the effective string model to describe bound states of two mobile partons connected by a strongly confining string. Our focus is on two holes, but as we discuss in detail in Appendix \ref{SuppMatToyModel}, the formalism we develop is also applicable to pairs of more general partons, in particular spinon-chargon pairs \cite{Bohrdt2021PRL}. In the second main part of the article, Sec.~\ref{secResults}, we present results from our analytical formalism and discuss possible implications for general pairing mechanisms in doped AFMs. We close with a summary and outlook in Sec.~\ref{secSummary}.

\section{Microscopic models}
\label{secMicModel}
In this article we introduce and solve an effective theory describing bound states of holes. As described in detail in Sec.~\ref{secEffStringModel}, we will make approximations on the level of both the Hilbert space and the Hamiltonian. Nevertheless, our starting point are microscopic models of doped AFM Mott insulators to which, we argue, our results apply within some approximations. Critical minds should simply view these models as motivating our effective theory, although our numerical analysis in Ref.~\cite{Bohrdt2022pairingPrep} indicates remarkable similarities with the semi-analytical predictions derived here.

The system most closely related to our effective theory is constituted by the 2D $t-J_z$ model on a square lattice, with Hamiltonian
\begin{multline}
\H_{t-J_z} = - t ~ \hat{\mathcal{P}} \sum_{\ij} \sum_\sigma \l  \hat{c}^\dagger_{\vec{i},\sigma} \c_{\vec{j},\sigma} + \text{h.c.} \r  \hat{\mathcal{P}}  + \\
+ J_z \sum_{\ij} \hat{S}^z_{\vec{i}} \hat{S}^z_{\vec{j}} - \frac{J_z}{4} \sum_{\ij} \hat{n}_{\vec{i}} \hat{n}_{\vec{j}} .
\label{eqHtJz}
\end{multline}
Here $\c_{\vec{j},\sigma}$ defines the underlying particles, with spin-index $\sigma=\uparrow,\downarrow$ and density $\hat{n}_{\vec{j}} = \sum_\sigma \hat{c}^\dagger_{\vec{j},\sigma} \c_{\vec{j},\sigma}$; the tunneling amplitude $t$ describes hopping of these particles, and $\hat{\mathcal{P}}$ projects on a sector with a given total number of doublons, holes and singly-occupied sites. $J_z > 0$ is an AFM Ising coupling between the spins $\hat{S}_{\vec{j}}^z = \sum_\sigma (-1)^\sigma \hat{c}^\dagger_{\vec{j},\sigma} \c_{\vec{j},\sigma}$ which we assume to be antiferromagnetic throughout.

In principle the Hamiltonian in Eq.~\eqref{eqHtJz} can be defined with particles $\c_{\vec{j},\sigma}$ of any exchange statistics, bosonic or fermionic. While this makes no difference for zero and one mobile dopant, the statistics plays an important role if two dopants of the same type are considered. Since the fermionic case is more closely related to the celebrated Fermi-Hubbard model with its AFM ground state at half filling, it usually takes center stage. However, we find it instructive to consider the bosonic version -- without any direct connection to a Hubbard Hamiltonian -- as well. In fact, quantum simulators using ultracold atoms have been proposed which allow the realization of both cases in experiments \cite{Duan2003,Trotzky2008,Gorshkov2011tJ,GuardadoSanchez2017}.

\begin{figure}[t!]
\centering
\epsfig{file=FIG1.pdf, width=0.49\textwidth}
\caption{We work in an effective Hilbertspace consisting of pairs of dopants (red and blue) connected by a string $\Sigma$ on a square lattice. Every state $\ket{\vec{x}_1,\Sigma}$ avoiding double occupancies of any site with two dopants is associated with a unique state $\ket{\Psi(\vec{x}_1,\Sigma)}$ in the microscopic model. (a) A typical example with string length $\ell_\Sigma=3$. (b) Rare loop configurations leading to double-occupancies of dopants have no correspondence in the microscopic model.
}
\label{figStrings}
\end{figure}

Likewise, we can consider the AFM $t-J$ model in a 2D square lattice with arbitrary underlying statistics. Its Hamiltonian is given by
\begin{multline}
\H_{t-J} = - t ~ \hat{\mathcal{P}} \sum_{\ij} \sum_\sigma \l  \hat{c}^\dagger_{\vec{i},\sigma} \c_{\vec{j},\sigma} + \text{h.c.} \r \hat{\mathcal{P}} + \\
+ J \sum_{\ij} \hat{\vec{S}}_{\vec{i}} \cdot \hat{\vec{S}}_{\vec{j}} - \frac{J}{4} \sum_{\ij} \hat{n}_{\vec{i}} \hat{n}_{\vec{j}}.
\label{eqDefHtJ}
\end{multline}
Now $J>0$ is the strength of antiferromagnetic SU(2)-invariant Heisenberg interactions between spins $\hat{\vec{S}}_{\vec{j}} = \sum_{\alpha,\beta} \hat{c}^\dagger_{\vec{j},\alpha} \frac{1}{2} \vec{\sigma}_{\alpha \beta} \c_{\vec{j},\beta}$.

In both the $t-J$ and $t-J_z$ models, different types of dopants can be considered. The most often studied case, which is also our primary focus, constitutes pairs of two indistinguishable holes; owing to the particle-hole symmetry of the models, one can interchangeably consider two indistinguishable doublons however. A second, closely related case corresponds to Hubbard-Mott excitons \cite{Huang2020arXiv}, where pairs of doublons and holes can form. In this case, exchange statistics again plays no role on the level of the $t-J_{(z)}$ model, since the dopants define distinguishable conserved particles. Experimentally, all these situations can be addressed by state-of-the-art ultracold atom experiments \cite{Tarruell2018,Bohrdt2021PWA}.

\section{Effective string model}
\label{secEffStringModel}
In this section, we introduce an effective string model to describe tightly bound pairs of two holes. Our approach is motivated by considering two dopants moving in a N\'eel state, modeled by the 2D $t-J^z$ or $t-J$ model. To make analytical progress, we perform approximations on the effective Hilbert space (see Fig.~\ref{figStrings}) as well as the effective Hamiltonian. While we expect our approximate description to be most accurate for the $t-J^z$ model, it should also capture the essential physics of related models, such as the $t-J$ model, as long as charge fluctuations dominate, $t \gg J$, and they feature strong local AFM correlations at zero doping. In these cases, the geometric string approach, developed originally for single dopants, can be applied \cite{Bulaevskii1968,Brinkman1970,Shraiman1988a,Grusdt2018PRX,Chiu2019Science}.

In the subsequent sections, we will always talk about paired holes and consider different kinds of statistics. However, as discussed in Sec.~\ref{secMicModel}, these can interchangeably be considered as general types of dopants, or even more generally as two partons.

\subsection{String Hilbert space and effective Hamiltonian}
\label{subsecStringHilbertSpace}
In a perfect N\'eel state the motion of a hole leaves behind a string $\Sigma$ of displaced spins, which allows to associate distinct hole trajectories with orthogonal states $\ket{\Psi(\Sigma)}$ in the quantum many-body system (Fig.~\ref{figStrings}); here strings denote hole trajectories with all self-retracing components removed. Exceptions, where two different strings $\Sigma_1 \neq \Sigma_2$ correspond to identical many-body states $\ket{\Psi(\Sigma_1)} = \ket{\Psi(\Sigma_2)}$, are associated with so-called Trugman loops \cite{Trugman1988}. Since the number of Trugman loops is small compared to the exponentially growing number of string states (Tab.~\ref{tabHoleDbleOccpcy}), their effect is generally found to be small \cite{Trugman1988,Grusdt2018PRX}. Although in exceptional cases the small effect of loops can still dominate, e.g. for the very narrow center-of-mass dispersion of a single hole in the 2D $t-J^z$ model \cite{Trugman1988,Chernyshev1999}, we neglect such loops in the following and study only the dominant effects of string formation. Loop effects can be re-introduced perturbatively in the end \cite{Grusdt2018PRX}.

\subsubsection{Hilbert space}
While the set of string states $\{ \ket{\Psi(\Sigma)} \}$ defines an over-complete basis, we approximate our Hilbertspace and formally define a set of two-hole string states:
\begin{equation}
 \ket{\vec{x}_1, \Sigma}, \qquad \vec{x}_1 \in \mathbb{Z}^2, \quad \Sigma \in {\rm BL}.
 \label{eqTwoHoleStringBasis}
\end{equation}
Here $\vec{x}_1$ denotes the location of the first hole in the 2D square lattice, and $\Sigma$ is the string which connects $\vec{x}_1$ to $\vec{x}_2 = \vec{x}_1 + \vec{R}_\Sigma$ at its opposite end (Fig.~\ref{figStrings}). The strings $\Sigma$ can be represented by the sites of a Bethe lattice (BL), or Cayley tree, with coordination number $z=4$, see Fig.~\ref{figHeffHopping} (a). Similar to the construction of the celebrated Rokhsar-Kivelson quantum dimer model \cite{Rokhsar1988}, we postulate that the new basis states are orthonormal,
\begin{equation}
 \bra{\vec{x}_1', \Sigma'} \vec{x}_1, \Sigma \rangle = \delta_{\Sigma,\Sigma'} \delta_{\vec{x}_1,\vec{x}_1'}.
\end{equation}

Every new basis state is associated with a unique microscopic two-hole state $\ket{\Psi(\vec{x}_1,\Sigma)}$ in the original $t-J^z$ or $t-J$ model. Some states in the new model describe unphysical double occupancies with holes: in this case the associated state becomes $\ket{\Psi(\vec{x}_1,\Sigma)}=0$ (Fig.~\ref{figStrings}). The fraction of such strings is relatively small, however, and decreases with increasing string lengths (Tab.~\ref{tabHoleDbleOccpcy}).

\begin{table}
\begin{center}
\begin{tabular}{ |c|c|c|c| }
 \hline
$\ell_\Sigma$ &no. of states $d(\ell_\Sigma)$ & double-occupancies & Trugman loops \\
 \hline
$1$ & $4$ & 	$0 $ & $0$ \\
$2$ & $12$ & 	$0$ & $0$ \\
$3$ & $36$ & 	$0$ & $0$ \\
$4$ & $108$ & $7.4 \%$ & $0$ \\
$5$ & $324$ & $0$ & $0$ \\
$6$ & $972$ & $2.5 \%$ & $0$ \\
$7$ & $2,916$ & $0$ & $0.55 \%$ \\
$8$ & $8,748$ & $0.91 \%$ & $1.3 \%$ \\
\hline
\end{tabular}
\end{center}
\caption{\textbf{Imperfections of the string model.} The string basis with two distinguishable holes includes states $\ket{\vec{x}_1,\Sigma}$ corresponding to unphysical double-occupancies of holes in the associated microscopic states $\ket{\Psi(\vec{x}_1,\Sigma)}$. Their relative fraction of all string states $N(\ell_\Sigma)$ of a given length $\ell_\Sigma$ is indicated. The relative number of Trugman loop configurations \cite{Trugman1988} is also shown.
}
\label{tabHoleDbleOccpcy}
\end{table}

So far we work in first quantization and assign separate labels to the two holes. Later (see \ref{SubSecIndistHoles}) we will generalize our approach to situations with indistinguishable holes, with bosonic or fermionic statistics.

\subsubsection{Effective Hamiltonian}
Next we define the effective Hamiltonian $\H^{\rm eff}$ in the approximated Hilbert space of the string model. To this end we require that the matrix elements satisfy:
\begin{equation}
 \bra{\vec{x}_1',\Sigma'} \H^{\rm eff} \ket{\vec{x}_1,\Sigma} = \bra{\Psi(\vec{x}_1',\Sigma')} \H \ket{\Psi(\vec{x}_1,\Sigma)},
 \label{eqHeffMatrixElements}
\end{equation}
where $\H$ corresponds to the respective microscopic model Hamiltonian ($t-J^z$, $t-J$,...).

As a result of the microscopic nearest-neighbor (NN) hopping $t_2$ of hole $2$ at site $\vec{x}_2$ on the lattice, we obtain NN hopping on the Bethe lattice (Fig.~\ref{figHeffHopping} (b)):
\begin{equation}
\H^{\rm eff}_{t,2}= t_2 \sum_{\vec{x}_1} \ket{\vec{x}_1}\bra{\vec{x}_1} \otimes \sum_{\langle \Sigma',\Sigma \rangle} \ket{\Sigma'} \bra{\Sigma} + \text{h.c.}
\end{equation}
The NN hopping $t_1$ of hole 1 gives rise to a correlated NN tunneling of $\vec{x}_1$ and a simultaneous change of the first string segment of $\Sigma \to \Sigma^{(1)}(\vec{x}_1',\vec{x}_1,\Sigma)$:
\begin{equation}
\H^{\rm eff}_{t,1}= t_1 \sum_{\langle \vec{x}_1',\vec{x}_1 \rangle} \sum_{\Sigma} \ket{\vec{x}_1'} \bra{\vec{x}_1} \otimes \ket{\Sigma^{(1)}(\vec{x}_1',\vec{x}_1,\Sigma)} \bra{\Sigma} + \text{h.c.} .
\end{equation}
If the first string segment in $\Sigma$ points along $\vec{x}_1'-\vec{x}_1$, it is removed to obtain $\Sigma^{(1)}$ (Fig.~\ref{figHeffHopping} (c)); otherwise, the string is extended by adding a new string segment pointing along $\vec{x}_1'-\vec{x}_1$ at the beginning of $\Sigma$ to obtain $\Sigma^{(1)}$ (Fig.~\ref{figHeffHopping} (d)).

\begin{figure}[t!]
\centering
\epsfig{file=FIG2.pdf, width=0.45\textwidth}
\caption{The hopping part of the effective Hamiltonian $\H^{\rm eff}_t$ describes NN tunneling of holes 1 ($t_1$) and 2 ($t_2$) on the square lattice. The string $\Sigma$ from hole 1 to 2 changes accordingly. We illustrate a typical initial state (a), which is coupled to neighboring string states on the Bethe lattice by $t_2$ (b). The coupling $t_1$ creates a string state $\Sigma^{(1)}$ which is a further neighbor of $\Sigma$ on the Bethe lattice, either by re-tracing (b) or extending (d) the first string segment.
}
\label{figHeffHopping}
\end{figure}

Similarly, we obtain the potential terms $\H^{\rm eff}_{J}$ in the effective Hamiltonian. They do not change the positions $\vec{x}_{1,2}$ of the holes, and we neglect off-diagonal matrix elements \eqref{eqHeffMatrixElements} for which $\Sigma' \neq \Sigma$. Hence, formally we can write:
\begin{equation}
\H^{\rm eff}_{J} = \sum_{\vec{x}_1} \sum_\Sigma V_\Sigma ~ \ket{\vec{x}_1,\Sigma} \bra{\vec{x}_1,\Sigma},
\end{equation}
where $V_\Sigma$ is a function of $\Sigma$ on the Bethe lattice only. Note that Trugman loops correspond to local minima in the Bethe lattice potential $V_\Sigma$, which allows for a systematic tight-binding treatment of Trugman loops within our model so far, even when $t_{1,2} \gg J$ or $J_z$ \cite{Grusdt2018PRX}.

The complete effective Hamiltonian we consider is
\begin{equation}
\H^{\rm eff} = \H^{\rm eff}_{t,1} + \H^{\rm eff}_{t,2} + \H^{\rm eff}_{J}.
\label{eqHeff}
\end{equation}

\subsubsection{Linear string approximation}
\label{subsubsecLinStringThy}
Since the dimension of the string Hilbert space grows exponentially with the maximum string length $\ell_{\rm max}$,
\begin{equation}
d_{\rm tot}(\ell_{\rm max}) = \sum_{\ell_\Sigma=1}^{\ell_{\rm max}} \underbrace{4 \times 3^{\ell_\Sigma-1}}_{=d(\ell_\Sigma)},
\end{equation}
further approximations are required to make analytical progress. For a general string potential $V_\Sigma$, all states in the string Hilbert space are coupled to each other by the tunneling terms. Next, by simplifying the potential $V_\Sigma$, many symmetry sectors emerge which are only weakly coupled and can be described by a simpler effective Hamiltonian that will be derived.

Since we are mostly interested in the regime $t \gg J$, we expect that inhomogeneities of $V_\Sigma$ on the scale of $J$ play a sub-dominant role. Such fluctuations of $V_\Sigma$ from string to string (or site to site on the Bethe lattice) result from string-string interactions \cite{Wrzosek2021}, and appear like a weak disorder potential. We can include their effect on a mean-field level by averaging the potential over all strings of a given length:
\begin{equation}
 V(\ell) = \frac{1}{d(\ell)} \sum_{\Sigma: \ell_\Sigma=\ell} V_\Sigma.
 \label{eqVellDef}
\end{equation}
The resulting problem is highly symmetric since all branches on the Bethe lattice corresponding to the same string length are equivalent.

As a further approximation, we can estimate the string-length potential $V(\ell) \approx V_{\rm LST}(\ell)$ by considering only straight strings in Eq.~\eqref{eqVellDef}. Since string-string interactions are always attractive, this linear string theory (LST) estimate also defines an upper bound for the averaged potential:
\begin{equation}
 V(\ell) \leq V_{\rm LST}(\ell) = \frac{dE}{d\ell} \times (\ell_\Sigma-1) + g^{(0)}_{\rm cc} \delta_{\ell_\Sigma,1} + \mu_{\rm cc}.
 \label{eqStringPot}
\end{equation}
Here $dE / d\ell$ denotes the linear string tension, $g^{(0)}_{\rm cc}$ is a nearest-neighbor hole-hole interaction, and $\mu_{\rm cc}$ an overall energy shift. The overall energy of the two holes is measured relative to the undoped parent antiferromagnetic state.

In the case of a microscopic $t-J^z$ model, we obtain:
\begin{equation}
 \frac{dE}{d\ell} = J^z, \quad g^{(0)}_{\rm cc} = - \frac{J^z}{2}, \quad \mu_{\rm cc} = 4 J^z.
 \label{eqStringPottJz}
\end{equation}
More generally, we can derive the LST potential by applying the frozen spin approximation and expressing the potential in terms of local spin-spin correlations of the undoped parent antiferromagnet, see Refs.~\cite{Grusdt2018PRX,Grusdt2019PRB}. For a doped $J_1-J_2$ spin model on a square lattice, with NN hopping of the holes, this yields:
\begin{flalign}
 \frac{dE}{d\ell} &= 2 J_1  \l C_2 - C_1 \r + 2 J_2 \l C_1 + C_4 - 2 C_2 \r,  \\
 g^{(0)}_{\rm cc} &= J_1 C_1 - \frac{J}{4}, \label{eqdEdlJ1J2} \\
 \mu_{\rm cc} &= - 8 \l J_1 C_1 + J_2 C_2 - \frac{J}{4} \r.
\end{flalign}
Here $J_1$ and $J_2$ denote the NN and diagonal NNN Heisenberg couplings on the square lattice; $J$ denotes the strength of the local attraction $- J/4 \hat{n}_{\vec{i}} \hat{n}_{\vec{j}}$ in Eq.~\eqref{eqDefHtJ} and is treated as an independent parameter here. The correlators are $C_1 = \langle \hat{\vec{S}}_{\vec{i}} \hat{\vec{S}}_{\vec{i}+\vec{e}_x} \rangle$ (NN), $C_2 = \langle \hat{\vec{S}}_{\vec{i}} \hat{\vec{S}}_{\vec{i}+\vec{e}_x+\vec{e}_y} \rangle$ (NNN), $C_3 = \langle \hat{\vec{S}}_{\vec{i}} \hat{\vec{S}}_{\vec{i}+2 \vec{e}_x} \rangle$ and $C_4 = \langle \hat{\vec{S}}_{\vec{i}} \hat{\vec{S}}_{\vec{i}+2\vec{e}_x+\vec{e}_y} \rangle$; they depend on the ratio $J_1/J_2$ \cite{Richter2010}.

\begin{figure*}[t!]
\centering
\epsfig{file=FIG3.pdf, width=0.87\textwidth}
\caption{The rotational basis consists of superpositions of string states on the Bethe lattice, defined around a rotational center. There is (a) one $m_4$ quantum number around the central site, and (b) one $m_3$ quantum number per every node except the center. Angles $\nu \pi/2$ are measured relative to the $x$-axis and $\lambda \pi/2$ relative to the preceding string segment. Note that in (b) only one branch of the Bethe lattice, with fixed $\nu$, is shown.}
\label{figRotationalBasis}
\end{figure*}

\subsubsection{Momentum basis}
So war we have defined the string model in the real-space basis. Because of the overall translational symmetry of the model, $[\H^{\rm eff}, \hat{T}]=0$ where the operator $\hat{T}$ translates both holes and the string by a discrete lattice vector, we can also work in the momentum-space basis. We define the following total momentum states,
\begin{equation}
  \ket{\vec{k}, \Sigma} = \frac{1}{\sqrt{V}} \sum_{\vec{x}_1} e^{i \vec{k} \cdot \vec{x}_1} \ket{\vec{x}_1, \Sigma},
\end{equation}
where $V=L^2$ denotes the total area of the square lattice.

By the symmetry, the effective Hamiltonian is block-diagonal, with entries
\begin{equation}
 \H^{\rm eff}(\vec{k}) = \sum_{\Sigma,\Sigma'} \bra{\vec{k},\Sigma'} \H^{\rm eff} \ket{\vec{k},\Sigma} ~ \ket{\Sigma'}\bra{\Sigma},
 \label{eqHkeff}
\end{equation}
and we can calculate the bound state properties independently at different total momenta $\vec{k}$.

\subsection{Rotational excitations and truncated basis}
\label{subsecRotExct}
In the simplified string potential $V(\ell_\Sigma)$ all strings of a given length are equivalent. This symmetry is conserved by the hopping term $t_2$ of the second hole, which leads to delocalization on the Bethe lattice. Hence, as long as the first hole cannot tunnel, i.e. for $t_1=0$, our problem with an exponentially large Hilbert space can be mapped to a single particle problem on a semi-infinite one-dimensional lattice $\ell_\Sigma=1,2,...., \infty$ in a central potential $V(\ell_\Sigma)$, see \cite{Grusdt2018PRX}. This situation is of great relevance for describing spinon-chargon bound states in an effective string basis \cite{Grusdt2018PRX,Bohrdt2021arXiv}, as will be discussed further in \ref{subsubSpinonChargon}.

In general, the tunneling of the first hole, $\propto t_1$, will break the symmetry between the string states, because some couple to longer, others to shorter strings depending on their orientation. Nevertheless, we find it useful to work in a basis of string states which takes into account the equivalence of strings when $t_1=0$: As we show later, for $C_4$-rotation invariant total momenta $\vec{k}^{\rm C4IM}$ (C4IM), even $t_1 \neq 0$ keeps the symmetry intact. Away from the C4IM, $t_1$ introduces weak couplings between the symmetric eigenstates, which we will explicitly take into account.

The quantum numbers characterizing the symmetric model correspond to one discrete rotational eigenvalue for each node in the Bethe lattice \cite{Grusdt2018PRX}. At the first node, around the central site $\ell_\Sigma=0$ in the Bethe lattice, one obtains a discrete $C_4$ rotational symmetry with eigenvalues $\exp(i m_4 \pi/2)$ where $m_4=0,1,2,3$. The symmetric states with string length $\ell_\Sigma=1$ are
\begin{equation}
 \ket{\vec{k},\ell_\Sigma=1,m_4} = \frac{1}{2} \sum_{\nu=0}^3 e^{i m_4 \nu \pi/2} \ket{\vec{k},\ell_\Sigma=1,\nu},
\end{equation}
where $\nu \pi/2$ denotes the angle of the string segment $\Sigma$ with the $x$-axis, see Fig.~\ref{figRotationalBasis} (a).

Each higher node, around sites on the Bethe lattice corresponding to string length $\ell \geq 1$, is associated with a discrete $C_3$ rotational symmetry. The corresponding eigenvalues are $\exp(i m_3^{(\ell)} 2\pi/3)$, with $m_3^{(\ell)}=0,1,2$. The corresponding symmetric states are defined as
\begin{equation}
 \ket{\vec{k},\ell_\Sigma,m_3^{(\ell)}} = \frac{1}{\sqrt{3}} \sum_{\lambda=1}^3 e^{i m_3^{(\ell)} \lambda 2\pi/3} \ket{\vec{k},\ell_\Sigma-1,\lambda},
\end{equation}
where $\lambda \pi /2$ denotes the angle between the previous and last string segment, see Fig.~\ref{figRotationalBasis} (b).

The most general states are defined by an entire set of angular momentum quantum numbers,
\begin{equation}
\vec{m} = (m_4, m_3^{(1)}, m_3^{(2)},...).
\end{equation}
I.e. instead of the original basis states $\ket{\vec{x},\Sigma}$, we work in the basis: $\{ \ket{\vec{k}, \ell_\Sigma, \vec{m}_\Sigma} \}$.

As a final step, we simplify our problem further by working with a truncated set of basis states. We discard non-trivial $m_3$ rotational excitations on nodes higher than one, i.e. we only take into account the quantum number $m_3^{(1)}$ and set all $m_3^{(n)}=0$ for $n\geq 2$. The truncated basis consists of the states
\begin{equation}
 \{ \ket{\vec{k}, \ell_\Sigma, m_4, m_3^{(1)}} \} \qquad \text{truncated~basis}
 \label{eqTruncBasis}
\end{equation}
and allows us to work with very large cut-offs in the string length, see Fig.~\ref{figTruncatedBasis}.

The main motivation for discarding higher rotational states is their higher energy in the purely symmetric limits (at $\vec{k}^{\rm C4IM}$ or when $t_1=0$). By studying the importance of the $m_3^{(1)}$ states within our reduced basis, we obtain an estimate for the effect of higher rotational manifolds in general.

We note that for cases with distinguishable holes, or more generally distinguishable partons, one should place the heavier parton (corresponding to a weaker hopping element) in the center. Around this first hole, labeled by $n=1$, the rotational states can then be defined. In the case of spinon-chargon pairs in the context of the $t-J$ \cite{Bohrdt2021PRL} or $t-J_z$ model \cite{Grusdt2018PRX}, which we discuss in more detail below, this is typically the spinon since $t < J_{(z)}$ \cite{Beran1996}.

We emphasize that for $t_1=t_2$, by choosing to define rotational states around parton $n=1$, the symmetry between the two partons is explicitly broken in the truncated basis; we will refer to this ansatz as the \emph{asymmetric approximation}. Below, in Sec.~\ref{SubSecIndistHoles}, we will show how (anti-) symmetrization of both cases, with parton $n=1$ and $n=2$ in the center respectively, leads to a more accurate truncated basis in the case $t_1=t_2$, going beyond the asymmetric approximation.

\begin{figure}[t!]
\centering
\epsfig{file=FIG4.pdf, width=0.38\textwidth}
\caption{The truncated basis consists of rotational states with $m_4$ and $m_3^{(1)} \equiv m_3$ quantum numbers only, defined around the first hole (asymmetric approximation). Later we will symmetrize this basis to treat both chargons on an equal footing. Solid lines indicate non-zero matrix elements which preserve the rotational symmetries in a linear string model. The state with $\ell_\Sigma=0$ can be formally added, but should be removed to describe cases where two holes cannot occupy the same site.}
\label{figTruncatedBasis}
\end{figure}

\subsection{Distinguishable partons: Mott-Hubbard excitons}
\label{subsecDistHoles}
Now we are in a position to further simplify the effective string problem, still assuming distinguishable holes and working in the asymmetric approximation. The first step is to calculate all required matrix elements in Eq.~\eqref{eqHkeff} for the truncated basis from Eq.~\eqref{eqTruncBasis}.

The string potential is completely diagonal, and we simply get
\begin{multline}
\bra{\vec{k}, \ell_\Sigma', m_4',m_3'} \H^{\rm eff}_J \ket{\vec{k}, \ell_\Sigma, m_4,m_3} = \\
=  V(\ell_\Sigma) ~ \delta_{\ell_\Sigma',\ell_\Sigma} \delta_{m_4',m_4} \delta_{m_3',m_3}.
\end{multline}
The tunneling of the second hole, $\propto t_2$, is diagonal in the rotational quantum numbers by construction \cite{Grusdt2018PRX} but changes the string length by one unit. For $\ell_\Sigma \geq 2$ we obtain
\begin{multline}
\bra{\vec{k}, \ell_\Sigma', m_4',m_3'} \H^{\rm eff}_{t,2} \ket{\vec{k}, \ell_\Sigma, m_4,m_3} = \\
= \sqrt{3} t_2  \l \delta_{\ell_\Sigma',\ell_\Sigma-1} + \delta_{\ell_\Sigma',\ell_\Sigma+1} \r  \delta_{m_4',m_4} \delta_{m_3',m_3},
\end{multline}
where the pre-factor $\sqrt{3} = \sqrt{z-1}$ related to the coordination number $z=4$ of the square lattice appears. The states with $\ell_\Sigma=1$ only have an $m_4$ quantum number and couple to $\ell_\Sigma'=2$ states, i.e.
\begin{equation}
 \bra{\vec{k}, \ell_\Sigma', m_4',m_3'} \H^{\rm eff}_{t,2} \ket{\vec{k}, 1, m_4} =
 \sqrt{3} t_2  ~ \delta_{\ell_\Sigma',2}  ~ \delta_{m_4',m_4} ~ \delta_{m_3',0};
\end{equation}
note that $\ell_\Sigma=0$ states are not included here since we assume hard-core holes. Overall, we obtain for the hopping amplitude $J_-^{(2)}$ of the second hole from $\ell_\Sigma$ to $\ell_\Sigma-1$,
\begin{multline}
J_-^{(2)}(\vec{k};m_4', m_3'; m_4, m_3)  \equiv \\
\equiv  \bra{\vec{k}, \ell_\Sigma-1, m_4',m_3'} \H^{\rm eff}_{t,2} \ket{\vec{k}, \ell_\Sigma, m_4,m_3} = \\
= \sqrt{3} t_2 ~  \delta_{m_4',m_4} \delta_{m_3',m_3},
\label{eqJmin2}
\end{multline}
and $J_+^{(2)} = (J_-^{(2)})^*$ for the reverse process, $\ell_\Sigma \to \ell_\Sigma+1$.

More complicated matrix elements $\propto t_1$ are generated by the hopping of the first hole, which may change both angular momenta $m_4$ and $m_3$. We only need to calculate the matrix elements going from $\ell_\Sigma$ to $\ell_\Sigma' = \ell_\Sigma-1$, which reduces the calculational workload since the final state must have $m_3'=0$. The remaining non-zero matrix elements $J_+^{(1)} = (J_-^{(1)})^*$ describe the reverse process, going from $\ell_\Sigma$ to $\ell_\Sigma' = \ell_\Sigma+1$, and are directly obtained from the former by complex conjugation. We find from a detailed calculation
\begin{multline}
 	J_-^{(1)}(\vec{k};m_4', m_3'; m_4, m_3)  \equiv \\
\equiv  \bra{\vec{k}, \ell_\Sigma-1, m_4',m_3'} \H^{\rm eff}_{t,1} \ket{\vec{k}, \ell_\Sigma, m_4,m_3} = \\
=  \delta_{m_3',0} \frac{t_1}{4} \sum_{\nu=0}^3 e^{i \frac{\pi}{2} \nu (m_4 - m_4')} \chi_{m_4,m_3}(\nu,\vec{k}),
\end{multline}
where we worked in the string configuration basis to obtain
\begin{equation}
 \chi_{m_4,m_3}(\nu,\vec{k}) = \frac{e^{i \pi m_4}}{\sqrt{3}} \sum_{\nu'=1}^3  e^{ -i \nu' \frac{\pi}{2} m_4 + i \frac{2 \pi}{3} m_3 \nu' - i \vec{k} \cdot \vec{e}_{\nu-\nu'+2} }.
\end{equation}
In the last expression, we defined the unit vector $\vec{e}_\lambda$ for $\lambda=0,1,2,3 ~{\rm mod} ~4$ to point along $\lambda \pi/2$ relative to the $x$-axis, i.e. $\vec{e}_0 = \vec{e}_x$, $\vec{e}_1 = \vec{e}_y$, $\vec{e}_2 = -\vec{e}_x$ and $\vec{e}_3 = - \vec{e}_y$.

Summarizing, we obtain an effective Hamiltonian with the matrix elements calculated above,
\begin{multline}
    \H^{\rm eff}(\vec{k}) = \sum_{\ell_\Sigma}  \Biggl[ \sum_{\vec{m}}~ V(\ell_\Sigma) ~\ket{\ell_\Sigma,\vec{m}}\bra{\ell_\Sigma,\vec{m}} + \\
    + \sum_{\vec{m},\vec{m}'} \sum_{n=1,2}  \Bigl( J_-^{(n)}(\ell_\Sigma;\vec{k}; \vec{m}';\vec{m}) \\  \times \ket{\ell_\Sigma-1,\vec{m}'}\bra{\ell_\Sigma,\vec{m}} + \text{h.c.} \Bigr) \Biggr].
    \label{eqHeffDist}
\end{multline}
Here the summation over $\vec{m}$, $\vec{m}'$ is restricted to $m_4,m_3^{(1)}$ in our truncated basis; for later purposes we also formally included an $\ell_\Sigma$-dependence in $J^{(n)}_-$, although we found the expressions above to be independent of $\ell_\Sigma$. This Hamiltonian can be easily diagonalized using exact numerical techniques, with large cut-offs for the maximum string length $\ell_{\rm max}$.

\subsubsection{Soft-core holes and spinon-chargon pairs}
\label{subsubSpinonChargon}
Before proceeding, we note that the formalism introduced above for distinguishable holes is very general and can be extended to describe string-bound states of any two different partons. Examples of particular relevance include strongly bound holes in mixed-dimensional bilayers \cite{Bohrdt2021PWA,Bohrdt2021pairing,Hirthe2022arXiv} and one-hole spinon-chargon pairs \cite{Beran1996,Grusdt2018PRX,Bohrdt2020,Bohrdt2021PRL}. An accurate description of the latter is important, since their energy relative to the two-hole chargon-chargon string states determines whether tightly-bound pairs of holes are energetically favorable in an antiferromagnet.

In both examples we mentioned, the two partons may occupy the same site. While we already made an approximation and included states with two partons on one site in the effective string basis for string lengths $\ell_\Sigma >0$, we have not included the $\ell_\Sigma=0$ state with two partons on the same site and without a string. To treat more general parton models, the $\ell_\Sigma=0$ state can be easily included in our formalism, however. This leaves the form of the effective Hamiltonian Eq.~\eqref{eqHeffDist} unchanged, but the hopping elements now include couplings to $\ell_\Sigma=0$ and the potential term $V(\ell_\Sigma=0) \ket{\ell_\Sigma=0}\bra{\ell_\Sigma=0}$ has to be added as well.

Still assuming that both partons $n=1,2$ can only tunnel between NN sites on the physical square lattice, we obtain hopping elements from $\ell_\Sigma=1$ to $\ell_\Sigma=0$. For the first, central parton ($n=1$) we get
\begin{equation}
    J_-^{(1)}(\vec{k};m_4) |_{\ell_\Sigma=1 \rightarrow \ell_\Sigma=0} = \frac{t_1}{2} \sum_{\nu=0}^3 e^{i \nu \frac{\pi}{2} m_4} e^{- i \vec{k} \cdot \vec{e}_\nu}.
\end{equation}
Note that only one rotational index $m_4$ appears since the state $\ket{\vec{k},\ell=0}$ has no rotational quantum numbers, and the states $\ket{\vec{k},\ell=1,m_4}$ at $\ell=1$ have one $m_4$ number.

For the second, orbiting parton ($n=2$) we obtain
\begin{equation}
     J_-^{(2)}(\vec{k};m_4) |_{\ell_\Sigma=1 \rightarrow \ell_\Sigma=0} = 2 t_2 \delta_{m_4,0}.
\end{equation}
In contrast to the other hopping elements, see Eq.~\eqref{eqJmin2}, this matrix element has strength $2 = \sqrt{4}$ instead of $\sqrt{3}$. This is a consequence of the enhanced connectivity of the $\ell_\Sigma=0$ state, which couples to 4 instead of the usual 3 longer string states \cite{Bulaevskii1968,Grusdt2018PRX}.

Finally, for the spinon-chargon case in an antiferromagnet, it is natural to assume NNN hopping of the spinon due to spin-exchange processes in the Hamiltonian. An extension of our formalism to this case has been used in \cite{Bohrdt2021PRL}, and a self-contained derivation is provided in Appendix \ref{SuppMatToyModel} of this article (see also \cite{Bohrdt2021arXiv}).

\subsection{Indistinguishable holes (partons)}
\label{SubSecIndistHoles}
To use the effective parton theory for describing doped holes in the $t-J$ or $t-J_z$ model, we next consider the case of indistinguishable partons. I.e. the string states must be invariant, up to an overall sign, if holes $n=1$ and $n=2$ are exchanged. To treat both partons on equal footing, we need to go beyond the asymmetric approximation made earlier, where rotational basis states were defined around parton $n=1$.

The required (anti)-symmetrization procedure further complicates the use of the truncated basis. In this section we will show that a suitable generalization of the truncated basis allows to keep the rotational quantum numbers, and derive effective Hamiltonians separately for fermions and bosons. As a result, $m_4$ remains a good quantum number at C4IM.

Together, the obtained fermionic and bosonic eigenstates span the space of distinguishable partons. When $t_1=t_2$, any bosonic (fermionic) eigenstate we identify also constitutes an eigenstate of distinguishable partons -- however, the use of the (anti-) symmetrized truncated basis leads to an improved variational energy going beyond the asymmetric approximation from Sec.~\ref{subsecDistHoles}; there, the rotational basis was defined around just one of the two partons, breaking the symmetry implied by $t_1=t_2$.

The results of our approach will be presented in the subsequent section \ref{secResults}, which can be understood without following the technical details derived in the remainder of this section. However, subsection \ref{subsecMic2MacMapping} may be of interest, where we discuss the relation between effective bosonic / fermionic string states and the microscopic AFM $t-J_{(z)}$ models composed of bosons or fermions, respectively.

\subsubsection{First quantization formalism}
We start by defining the permutation operator which exchanges the labels of the holes,
\begin{equation}
    \P \ket{\vec{x}_1,\Sigma} = \ket{\vec{x}_2,\overline{\Sigma}},
    \label{eqDefP}
\end{equation}
where $\overline{\Sigma}$ describes the same string as $\Sigma$ but starting from the opposite end.

If we consider the case $t_1=t_2=t$ in our model, the Hamiltonian $\H^{\rm eff}$ in Eq.~\eqref{eqHeff} commutes with the permutation operator, $[\H^{\rm eff},\P]=0$. Hence every eigenstate belongs to one of two classes, either fermionic (with $\P$-eigenvalue $-1$) or bosonic (with $\P$-eigenvalue $+1$).

Making use of this symmetry is still challenging, however, due to the exponential size of the string Hilbert space. Moreover, the string states $\ket{\vec{k},\ell_\Sigma,\vec{m}}$ introduced above in general do not have well-defined exchange statistics, i.e. they are not eigenstates of $\P$. We will discuss below how this problem can be avoided and proper string-states can be defined.

\paragraph{String-length $\ell_\Sigma=1$ states.--}
The high-symmetry states at the shortest string length $\ell_\Sigma=1$ are of particular importance for defining quasiparticle weights of pairs \cite{Bohrdt2022pairingPrep} later on. At the C4IM $\vec{k}^{\rm C4IM}$ (for the considered square lattice these are $\vec{k}^{\rm C4IM} = \vec{0}$ and $\vec{k}^{\rm C4IM}=(\pi,\pi)$) one obtains
\begin{multline}
    \P ~ \ket{\vec{k}^{\rm C4IM}, \ell_\Sigma=1,m_4} = (-1)^{m_4} \\
   \times  e^{i (\vec{e}_x+\vec{e}_y)\cdot \vec{k}^{\rm C4IM} / 2 } ~ \ket{\vec{k}^{\rm C4IM}, \ell_\Sigma=1,m_4},
   \label{eqSpectralWeightsl1}
\end{multline}
i.e. these states are purely fermionic / bosonic respectively. Concretely, at $\vec{k}^{\rm C4IM} = \vec{0}$ the states $m_4=0,2$ ($m_4=1,3$) are bosonic (fermionic) respectively. Vice-versa, at $\vec{k}^{\rm C4IM} = (\pi,\pi)$ the states $m_4=0,2$ ($m_4=1,3$) are fermionic (bosonic).

It may seem counter-intuitive at first to have (spin-less) fermionic pairs of holes with $s$-wave pairing symmetry at $\vec{k}=(\pi,\pi)$. However, the center of mass momentum $\vec{k}=(\pi,\pi)$ effectively leads to an anti-symmetric wavefunction under exchange of the $A$ and $B$ sublattices. This ensures correct overall fermionic statistics even in the presence of an $s$-wave pairing symmetry, in the effective string model. We will discuss in Sec.~\ref{subsecRelationNeelPairs} how this relates to pairs in the microscopic N\'eel state, which breaks the discrete translational symmetry of the square lattice.

\paragraph{General string states.--}
For general string states, some additional insights can be gained at C4IM. For example, we find that the rotational ground states with $\vec{m}=\vec{0}$ satisfy: (i) all $\vec{k}=\vec{0}$ states are bosonic,
\begin{equation}
    \P \ket{\vec{k}=\vec{0} ,\ell_\Sigma, \vec{m}=\vec{0}} = \ket{\vec{k}=\vec{0} ,\ell_\Sigma, \vec{m}=\vec{0}};
\end{equation}
and, (ii), the statistics of states at $\vec{k}=\vec{\pi} = (\pi,\pi)$ alternates with the string length,
\begin{equation}
     \P \ket{\vec{k}=\vec{\pi} ,\ell_\Sigma, \vec{m}=\vec{0}} = (-1)^{\ell_\Sigma} \ket{\vec{k}=\vec{\pi} ,\ell_\Sigma, \vec{m}=\vec{0}}.
\end{equation}

\subsubsection{(Anti-) symmetrized truncated basis}
As explained above, the rotational basis states generally do \emph{not} have proper exchange statistics. To enforce the latter, we can explicitly (anti-) symmetrize the basis states by defining new basis states:
\begin{equation}
    \ket{\vec{k},\ell_\Sigma,\vec{m},\mu} = (1 + \mu \P) ~ \ket{\vec{k},\ell_\Sigma,\vec{m}}.
    \label{eqSymmBasis}
\end{equation}
Here $\mu =\pm 1$ denotes bosonic (fermionic) states when $\mu=+1$ ($\mu=-1$, respectively).
By restricting states in Eq.~\eqref{eqSymmBasis} to the lowest rotational quantum numbers $\vec{m}$, a new truncated basis with well-defined exchange statistics is obtained. However, the new states are no longer orthonormal, and we use a Gram-Schmidt procedure to construct an orthonormal basis (ONB) in the linear space spanned by states in Eq.~\eqref{eqSymmBasis}.

\paragraph{Gram-Schmidt procedure.--}
For each given set of good quantum numbers $\vec{k}$, $\ell_\Sigma$ and $\mu$ we use the standard Gram-Schmidt method to construct ONB states by mixing different rotational states $\vec{m}$. Denoting the new ONB states by $\tilde{\vec{m}}$, we have the general representation
\begin{equation}
    \ket{\vec{k},\ell_\Sigma,\tilde{\vec{m}},\mu} = \sum_{\vec{m}} \tilde{C}_{\vec{m},\tilde{\vec{m}}}(\vec{k},\ell_\Sigma,\mu) ~ \ket{\vec{k},\ell_\Sigma,\vec{m},\mu}
    \label{eqGramSchmidtDef}
\end{equation}
with coefficients determined by the matrix $\tilde{C}(\vec{k},\ell_\Sigma,\mu)$. Specifically, we choose
\begin{flalign*}
    \ket{\tilde{\vec{m}}&=\vec{0}} = \mathcal{N}_{0} \ket{\vec{m}=\vec{0}},\\
    .&..\\
    \ket{\tilde{\vec{m}}} &= \mathcal{N}_{\tilde{\vec{m}}} \Biggl( \ket{\vec{m}} - \sum_{\tilde{\vec{M}}<\tilde{\vec{m}}} \ket{\tilde{\vec{M}}} \bra{\tilde{\vec{M}}} \left. \! \vec{m} \right\rangle \Biggr),
\end{flalign*}
where we dropped all labels $\vec{k}$ and $\ell_\Sigma$ for simplicity and $\mathcal{N}_{\tilde{\vec{m}}}$ denotes normalization constants.

\paragraph{Calculating overlaps of string states.--}
To calculate overlaps between the old and new basis states, we use the following explicit expression for overlaps between the symmetrized states defined in Eq.~\eqref{eqSymmBasis}:
\begin{multline}
    \bra{\ell,m_4',m_3'} \left. \ell,m_4,m_3 \right\rangle = 2 \delta_{m_4',m_4} \delta_{m_3',m_3} + \frac{\mu}{4 \times 3^{\ell-1}} \\
    \times \sum_{\nu,\nu'=0}^4 e^{i \frac{\pi}{2}(\nu m_4 - \nu' m_4')} \sum_{\lambda,\lambda'=1}^3 e^{i \frac{2 \pi}{3} (\lambda m_3 - \lambda' m_3')} \\
    \times \Biggl( \sum_{\substack{\Sigma:~~ \ell_\Sigma=\ell \\ \varphi_\Sigma=\frac{\pi}{2} \nu',~ \vartheta_\Sigma=\frac{\pi}{2}\lambda' \\ \overline{\varphi}_\Sigma=\frac{\pi}{2} \nu,~ \overline{\vartheta}_\Sigma=\frac{\pi}{2}\lambda}} e^{i \vec{k} \cdot \vec{R}_\Sigma} ~~
    + \sum_{\substack{\Sigma:~~ \ell_\Sigma=\ell \\ \varphi_\Sigma=\frac{\pi}{2} \nu,~ \vartheta_\Sigma=\frac{\pi}{2}\lambda \\ \overline{\varphi}_\Sigma=\frac{\pi}{2} \nu',~ \overline{\vartheta}_\Sigma=\frac{\pi}{2}\lambda'}} e^{-i \vec{k} \cdot \vec{R}_\Sigma}
    \Biggr)
    \label{eqOvlpFormula}
\end{multline}
again dropping $\vec{k}$ and $\mu$ labels for simplicity; here we only included $m_4$ and $m_3^{(1)}$ rotational states. In the last line we sum over all string configurations $\Sigma$ whose length $\ell_\Sigma = \ell$ equals the length $\ell$ in the considered string sector, and with constraints on the last string segments (see Fig.~\ref{figStringOvlp}): $\varphi_\Sigma$ denotes the angle of the first string segment starting at the central hole ($n=1$) measured relative to the $x$-axis; $\vartheta_\Sigma$ denotes the angle of the second string segment starting at the central hole ($n=1$) measured relative to first string segment; $\vec{R}_\Sigma = \vec{x}_2 - \vec{x}_1$ denotes the vector in the 2D square lattice connecting the first hole at $\vec{x}_1$ to the second hole at $\vec{x}_2$; the angles $\overline{\varphi}_\Sigma$ and $\overline{\vartheta}_\Sigma$ are defined like $\varphi_\Sigma$ and $\vartheta_\Sigma$ but starting from the other hole ($n=2$), as shown in Fig.~\ref{figStringOvlp}.

\begin{figure}[t!]
\centering
\epsfig{file=FIG5.pdf, width=0.38\textwidth}
\caption{The orientations of the first two string segments starting from hole 1 in (a) (hole 2 in (b)) are labeled by the angle $\varphi_\Sigma$ ($\overline{\varphi}_\Sigma$) relative to the $x$-axis, and the angle $\vartheta_\Sigma$ ($\overline{\vartheta}_\Sigma$) measured relative to the preceding string segment.}
\label{figStringOvlp}
\end{figure}

For maximum string lengths up to around $\ell \simeq 13$ or slightly larger, the overlaps in Eq.~\eqref{eqOvlpFormula} can be calculated by an exact summation over all string states. For larger maximum string lengths, Metropolis Monte-Carlo sampling over string configurations can be used to obtain good estimates for the overlaps.

\subsubsection{Effective Hamiltonian}
To obtain the Hamiltonian in the Gram-Schmidt ONB basis $\ket{\vec{k},\ell_\Sigma,\tilde{\vec{m}}}$, we project the effective Hamiltonian to the symmetrized subspace,
\begin{equation}
    \tilde{\mathcal{H}} = \tilde{P} ~ \H^{\rm eff}~ \tilde{P}
\end{equation}
where $\tilde{P} = \sum_{\ell_\Sigma} \sum_{\tilde{\vec{m}}} \ket{\ell_\Sigma,\tilde{\vec{m}}} \bra{\ell_\Sigma,\tilde{\vec{m}}}$, and directly calculate its matrix elements $\bra{\ell_\Sigma',\tilde{\vec{m}}'}\tilde{\mathcal{H}} \ket{\ell_\Sigma,\tilde{\vec{m}}}$.

To this end, we first use Eq.~\eqref{eqGramSchmidtDef} to express the ONB vectors in terms of the non-ONB states. The action of the Hamiltonian on the latter, $\H^{\rm eff} \ket{\vec{k},\ell,\vec{m},\mu}$, is known from the case of distinguishable holes, noting that the (anti-) symmetrization $(1+\mu \P)$ is a linear operation. This leads to the processes $\propto J^{(n)}_{\pm}(\vec{k};\vec{m}',\vec{m})$ and $V(\ell_\Sigma)$.

Thus we obtain an expression for $\tilde{\mathcal{H}} \ket{\vec{k},\ell,\tilde{\vec{m}},\mu}$ as a sum over the non-ONB basis states $\ket{\vec{k},\ell,\vec{m},\mu}$ weighted by coefficients $\tilde{C}$ and couplings $J^{(n)}_\pm$ or $V(\ell_\Sigma)$. This allows us to express the non-ONB states in terms of ONB states again by
\begin{equation}
    \tilde{P} \ket{\vec{k},\ell_\Sigma,\vec{m},\mu} = \sum_{\tilde{\vec{m}}} C_{\tilde{\vec{m}},\vec{m}}(\vec{k},\ell_\Sigma,\mu) ~ \ket{\vec{k},\ell_\Sigma,\tilde{\vec{m}},\mu},
\end{equation}
with overlaps
\begin{equation}
    C_{\tilde{\vec{m}},\vec{m}}(\vec{k},\ell_\Sigma,\mu) = \bra{\vec{k},\ell_\Sigma,\tilde{\vec{m}}} \left. \vec{k},\ell_\Sigma,\vec{m} \right\rangle.
\end{equation}

Finally, this leads to the effective Hamiltonian in the ONB basis. It is most conveniently expressed by starting from the distinguishable-hole Hamiltonian in square matrix form:
\begin{equation}
    \uuline{H}= \uuline{V} + \l \uuline{J_-}  + \uuline{J_+} \r,
\end{equation}
where $\uuline{V}$ is a diagonal matrix with entries $V(\ell_\Sigma)\otimes 1_{\vec{m}}$. Moreover, $\uuline{J_{\pm}}$ are upper (lower) block-band matrices connecting states at $\ell_\Sigma$ and $\ell_\Sigma-1$ ($\ell_\Sigma+1$) respectively, and mixing $\vec{m}$ quantum numbers in the truncated basis; i.e. the matrix elements of $\uuline{J_{\pm}}$ are given by $J_\pm^{(1)}+J_\pm^{(2)}$. Defining block-diagonal matrices $\uuline{C}$ and $\uuline{\tilde{C}}$ with rectangular blocks $C_{\tilde{\vec{m}},\vec{m}}(\ell_\Sigma)$ and $\tilde{C}_{\vec{m},\tilde{\vec{m}}}(\ell_\Sigma)$ at each block of states with string length $\ell_\Sigma$, we obtain for indistinguishable holes in the ONB basis:
\begin{equation}
    \uuline{\tilde{H}} = \uuline{C}~  \uuline{H} ~\uuline{\tilde{C}}.
    \label{eqHONB}
\end{equation}
Note that through the matrices $\uuline{C}$ and $\uuline{\tilde{C}}$ the Hamiltonian $\uuline{\tilde{H}}$ depends on the particle statistics $\mu$; i.e. different Hamiltonians are obtained for bosonic, fermionic and distinguishable holes.

\subsubsection{Quasiparticle weights}
For the eigenstates $\ket{\psi_n(\vec{k},\mu)}$ of the Hamiltonian \eqref{eqHONB} in the ONB basis, we define the corresponding quasiparticle weights by their overlaps squared with (anti-) symmetrized states at string length $\ell_\Sigma=1$:
\begin{equation}
    Z_{n}(\vec{k},m_4,\mu) = |\bra{ \psi_n(\vec{k},\mu) } \left. \vec{k},1,m_4,\mu \right\rangle|^2.
    \label{eqZnkm4Def}
\end{equation}
This result can be expressed conveniently in terms of the matrices $\tilde{C}_{\vec{m},\tilde{\vec{m}}}$ and $C_{\tilde{\vec{m}},\vec{m}}$ and projectors $\uuline{P_{\ell=1}^{m_4}}$ to the string-length $\ell=1$ states with rotational quantum number $m_4$. We obtain:
\begin{equation}
    Z_{n}(\vec{k},m_4,\mu) = \bra{ \psi_n(\vec{k},\mu) } \uuline{C} ~\uuline{P_{\ell=1}^{m_4}} ~\uuline{C}^\dagger \ket{ \psi_n(\vec{k},\mu)}.
\end{equation}

\subsubsection{Exchange statistics in microscopic and effective models}
\label{subsecMic2MacMapping}
So far, our discussion of exchange statistics was only on the level of the effective string model, where particle exchange is defined through the application on the permutation operator $\hat{\mathcal{P}}$, see Eq.~\eqref{eqDefP}. On the other hand, the effective string states $\ket{\vec{x},\Sigma}$ can be directly related to two-hole states $\ket{\Psi(\vec{x},\Sigma)}$ in a classical N\'eel background composed of constituents $\c_{\vec{j},\sigma}$, see Secs.~\ref{secMicModel} and \ref{subsecStringHilbertSpace}. To this end, one starts from the classical N\'eel state $\ket{\rm N}$. Next, the first hole is created at site $\vec{x}$ and the second hole is created next to $\vec{x}$ at $\vec{R}_2=\vec{x}+\vec{\delta}_1$ along the direction of the first string segment in $\Sigma$, by applying $\c_{\vec{x}+\vec{\delta}_1,\sigma_2} \c_{\vec{x},\sigma_1} \ket{\rm N}$. Then the second hole is moved around, along the directions $\vec{\delta}_n$ of subsequent string segments in $\Sigma$, by applying hopping terms $\sum_\sigma \hat{c}^\dagger_{\vec{R}_n+\vec{\delta}_n,\sigma} \c_{\vec{R}_n}$, where $\vec{R}_{n+1} = \vec{R}_n+\delta_n$.

To understand how the exchange statistics $\mu$ in the effective string model relates to the statistics of the microscopic constituents $\c_{\vec{j},\sigma}$ of the N\'eel AFM, we start from a parton construction representing the underlying spins as $\c_{\vec{j},\sigma} = \hat{h}^\dagger_{\vec{j}} \hat{f}_{\vec{j},\sigma}$, subject to the local constraints $\sum_\sigma \hat{f}^\dagger_{\vec{j},\sigma} \hat{f}_{\vec{j},\sigma} + \hat{h}^\dagger_{\vec{j}}\hat{h}_{\vec{j}}=1$. One can choose different combinations of parton statistics $\mu_{f,h} = \pm 1$, as long as $\mu_c = \mu_f \mu_h$, where $\mu_c=+1$ ($\mu_c=-1$) if $\c_{\vec{j},\sigma}$ are bosons (fermions).

The simplest way to map effective string states $\ket{\vec{x},\Sigma}$ to microscopic two-hole states $\ket{\Psi(\vec{x},\Sigma)}$ is to choose bosonic spinons $\mu_f=+1$, i.e. $\hat{f}_{\vec{j},\sigma}$ are Schwinger bosons, and $\mu_h=\mu_c$. The spinons $\hat{f}_{\vec{j},\sigma}$ keep track of how the spin pattern is distorted by the hole motion creating the string $\Sigma$, and since we chose bosonic spinons the order of the spins in the background is irrelevant. Exchanging the two chargons $\hat{h}_{\vec{x}_1} \hat{h}_{\vec{x}_2} = \mu_h \hat{h}_{\vec{x}_2} \hat{h}_{\vec{x}_1}$ keeps track of additional minus signs in case $\mu_h=-1$. To reflect the resulting exchange statistics $\mu_h=\mu_c$ correctly, the statistics of the effective string states should equal $\mu = \mu_c$.

Making the opposite choice $\mu_f = -1$ and $\mu_h=-\mu_c$, i.e. choosing fermionic spinons $\hat{f}_{\vec{j},\sigma}$, leads to the same result, $\mu = \mu_c$. However, in this case one has to keep track of exchange signs $\mu_f$ picked up when spin operators in the background are exchanged, which requires additional book keeping.

\section{Results}
\label{secResults}
We compare the predictions by our string-based model of fermionic hole pairing to fully numerically obtained two-hole spectra (from matrix product states) on four-leg cylinders in Ref.~\cite{Bohrdt2022pairingPrep}. There we find good agreement for the $t-J_z$ and $t-J$ models, and we discuss how accurately our string model is able to describe the numerical observations. The main focus in the present article is on understanding the predictions of our effective model, the accuracy of the truncated bases we use, and the analytical insights that can be gained from our calculations.

\begin{figure}[t!]
\centering
\epsfig{file=FIG6.pdf, width=0.47\textwidth}
\caption{Two-hole eigenstates from the truncated string basis calculation. We compare fermionic (red) and bosonic (blue) holes, and show results from the asymmetric approximation (gray). The truncated basis used for the calculations includes all $m_4=0,..,3$ and $m_3^{(1)}=0,1,2$ sectors, and string lengths up to $\ell_{\rm max}=11$. We considered $t/J_z=3$ and a string potential for an Ising background.}
\label{figSpectratJz3Comb}
\end{figure}

\begin{figure*}[t!]
\centering
\epsfig{file=FIG7.pdf, width=0.9\textwidth}
\caption{String-length distributions for fermionic (red) and bosonic (blue) holes in (a) the ground state around $\vec{k} = \vec{0}$ and (b) in the lowest-energy state around $\vec{k} = (\pi,\pi)$. Predictions by the asymmetric approximation are also shown (gray). The rotational quantum numbers of all states can be extracted by analyzing spectral weights, and our results are indicated in the legend: s-wave ($m_4=0$), p-wave ($m_4=1$), d-wave ($m_4=2$). Note in (a) the bosonic prediction coincides with the asymmetric approximation. Throughout we used the truncated basis including all $m_4=0,..,3$ and $m_3^{(1)}=0,1,2$ sectors, and string lengths up to $\ell_{\rm max}=13$. We considered $t/J_z=3$ and a string potential for an Ising background.}
\label{figStringLengthsFlatAndDisp}
\end{figure*}

\subsection{Two-hole spectra}
We start by calculating the two-hole eigenstates along a high-symmetry cut through the Brillouin zone in Fig.~\ref{figSpectratJz3Comb}. We compare our results for different statistics of the holes: fermionic (relevant to pairing in the doped Hubbard model) and bosonic (as a theoretical reference); the case of distinguishable dopants (relevant to Hubbard-Mott excitons \cite{Huang2020arXiv}) corresponds to the combined bosonic and fermionic eigenstates. As a theoretical reference, we also show results from the less accurate asymmetric approximation, i.e. without (anti-) symmetrizing rotational states in the truncated basis. Other parameters are identical, we assume $t/J_z=3$ and calculate the string potential Eq.~\eqref{eqStringPot} for a $t-J_z$ model i.e. using Eq.~\eqref{eqStringPottJz} for all cases.

We observe several striking features, discussed in more detail shortly: (i) in all cases, the lowest energy state has zero momentum $\vec{k}=\vec{0}$; (ii) the ground state of distinguishable holes is bosonic and is captured by the asymmetric approximation; (iii) its energy is significantly below the lowest-energy fermionic state; (iv) the lowest-energy states are highly dispersive, with an effective mass scaling as $M_{\rm hh} \propto 1/t$; (v) in addition to several strongly dispersing bands, we observe numerous exactly flat bands; (vi) while their energy differs between different hole statistics, they are always present; (vii) in the fermionic case, the flat band constitutes the lowest-energy state over a significant portion of the Brillouin zone, except around $\vec{k}=\vec{0}$.

Some comments are in order: Regarding (i), we note that for the fermionic case, the dispersive $\vec{k}=0$ band is only slight lower in energy than the lowest flat band. This competition is found to be even more pronounced for smaller values of $t/J_z$ (not shown). For $t \ll J_z$ we find that the dispersive state becomes degenerate with the flat-band state at $\vec{k}=0$.

Regarding (ii), note that the lower variational energy at $\vec{k} \neq 0$ of the bosonic ground state demonstrates higher accuracy of the (anti-) symmetrized states as compared to the asymmetric approximation.

Regarding (iii), it has been pointed out previously by Trugman \cite{Trugman1988} that the fermionic sign effectively frustrates the hopping Hamiltonian of two fermionic holes connected by a string. We believe this explains the increase of their energy relative to the bosonic or distinguishable holes, as observed in Fig.~\ref{figSpectratJz3Comb}.

Regarding (iv), it has long been understood that the ability of one hole to follow the other may lead to a highly mobile bound state with bandwidth $\propto t$ in the $t-J_z$ model. An estimate for the effective mass, $M_{\rm hh}^{-1} = t \sqrt{3}$, was derived for $t \gg J_z$ from a string model in Ref.~\cite{Bohrdt2021pairing}. The inclusion of quantum fluctuations in the $t-J$ model is expected to cause additional polaronic dressing of these states and a corresponding mass enhancement.

Regarding (v), the existence of flat bands corresponding to $M_{\rm hh} \to \infty$ in an effective string model has been predicted by Shraiman and Siggia \cite{Shraiman1988a}, although they used a slightly different procedure to truncate their string basis. This suggests that flat bands of hole pairs are not an artifact but robust excitations of the system. Indeed, in our recent numerical work \cite{Bohrdt2022pairingPrep} we found strong evidence for the existence of flat bands of hole pairs in the $t-J_z$ model.

As in many flat-band systems, we believe that destructive quantum interference between different paths underlies the formation of self-localized flat-band states. Notably, we checked that this does not limit the string length of the pair: As shown in Fig.~\ref{figStringLengthsFlatAndDisp}, the string-length distribution of the dispersive bound state around $\vec{k}=0$ (a) is qualitatively similar -- with a broad peak around $\ell_\Sigma = 3$ -- to the flat-band bound state around $\vec{k}=(\pi,\pi)$. While different statistics of the holes and $m_4$ quantum numbers usually correspond to small differences in the string-length histograms, their overall shape is always similar with an average string-length $\langle \ell_\Sigma\rangle \approx 3-4$ for $t/J_z=3$.

Regarding (vii), we note that additional polaronic dressing by spin-waves \cite{Kane1989} or phonons in a solid is expected to lower the energy of the flat-band state further compared to the dispersive band, since the large recoil energy $\propto t$ associated with the dispersive band suppresses polaronic dressing. Moreover, strong interactions between the pairs may favor occupying the self-localized flat-band states, where localization costs no kinetic energy while the interaction energy can be minimized.

\begin{figure*}[t!]
\centering
\epsfig{file=FIG8.jpg, width=0.99\textwidth}
\caption{Two-hole spectra for distinguishable holes connected by a string, calculated using the asymmetric approximation, along a high-symmetry cut through the Brillouin zone of the square lattice. The plots are obtained from a spectral decomposition and using our semi-analytical theory with a truncated basis including all $m_4$ and $m_3^{(1)}$ sectors, up to a maximum string length $\ell_{\rm max}=11$. We considered $t/J_z=3$ and a string potential for an Ising background.
}
\label{figSpectratJz3Dist}
\end{figure*}

\begin{figure*}[t!]
\centering
\epsfig{file=FIG9.jpg, width=0.99\textwidth}
\caption{Two-hole spectra as in Fig.~\ref{figSpectratJz3Dist}, but for fermionic holes. Again, all $m_4$ and $m_3^{(1)}$ sectors were included and a maximum string length $\ell_{\rm max}=11$ was used. We considered $t/J_z=3$ and a string potential for an Ising background.
}
\label{figSpectratJz3Fermi}
\end{figure*}

\begin{figure*}[t!]
\centering
\epsfig{file=FIG10.jpg, width=0.99\textwidth}
\caption{Two-hole spectra as in Figs.~\ref{figSpectratJz3Dist}, \ref{figSpectratJz3Fermi}, but for bosonic holes. Again, all $m_4$ and $m_3^{(1)}$ sectors were included and a maximum string length $\ell_{\rm max}=11$ was used. We considered $t/J_z=3$ and a string potential for an Ising background.
}
\label{figSpectratJz3Bose}
\end{figure*}

\subsubsection{Spectral weights}
So far we have only calculated the energies of string-paired eigenstates in Fig.~\ref{figSpectratJz3Comb}. To reveal the nature of different states, we calculate their spectral weights $Z_n(\vec{k},m_4)$ for different hole statistics, i.e. their overlaps squared with string-length $\ell_\Sigma=1$ states of the same symmetries, see Eq.~\eqref{eqZnkm4Def}. For presentation purposes, the obtained spectral lines are broadened with a Gaussian of width $\sigma=J_z/5$, and the integrated spectral weight per peak equals the corresponding quasiparticle weight. Our results are shown as momentum cuts for different values of $m_4$ in Figs.~\ref{figSpectratJz3Dist}-\ref{figSpectratJz3Bose}.

Before we discuss our results, we emphasize that while the overlaps with string-length one states of a given momentum $\vec{k}$ and rotational quantum number $m_4$ are well-defined, the paired eigenstates only have a well-defined rotational eigenvalue for C4IM. For general $\vec{k} \notin {\rm C4IM}$, the momentum $\vec{k}$ explicitly breaks the $C_4$ symmetry and different rotational states can hybridize. Moreover, the A- and B- sublattice degree of freedom on the square lattice allows to realize $s$/ $d$-wave states of indistinguishable fermions ($p$ / $f$-wave states of indistinguishable bosons) at general $\vec{k}$, whose relation to microscopic states in a N\'eel AFM we discuss further in subsection \ref{subsecRelationNeelPairs}.

In Fig.~\ref{figSpectratJz3Dist} we show our results for the case of distinguishable holes, relevant e.g. to models of Mott excitons \cite{Huang2020arXiv} or mixD bilayer pairing \cite{Bohrdt2021pairing}. In the $s$-wave channel (a) we observe only dispersive bands and a collapse of spectral weight at high energies around $\vec{k}=(\pi,\pi)$. The latter effect was recently found numerically in a microscopic bilayer model \cite{Bohrdt2021pairing}. In the other channels $m_4=1,2,3$ shown in (b)-(d) we observe a mix of many flat and highly dispersive bands, where $d$-wave states are flat around $\vec{k}=(\pi,\pi)$ and dispersive around the corners of the Brillouin zone $\vec{k}=(0,\pi)$.

In Fig.~\ref{figSpectratJz3Fermi} we show the same plot for fermionic holes. A first striking difference is the complete suppression of spectral weight around $\vec{k}=0$ (for $s$- and $d$-wave pairs) and around $\vec{k}=(\pi,\pi)$ (for $p$- and $f$-wave pairs). This follows directly from the fixed and alternating statistics of the string-length $\ell_\Sigma=1$ states at those momenta, derived in Eq.~\eqref{eqSpectralWeightsl1}. Moreover, we find in (a) that the collapse of spectral weight at high energies in the $s$-wave channel around $\vec{k}=0$ remains a robust feature for fermionic holes. This has been observed in earlier exact diagonalization calculations for the $t-J$ model \cite{Dagotto1990} and confirmed by our recent numerical DMRG study \cite{Bohrdt2022pairingPrep} for the $t-J$ and $t-J_z$ models.

At low energies in Fig.~\ref{figSpectratJz3Fermi} we only find spectral weight in the $p$-, $d$- and $f$-wave channels, corresponding to the lowest-energy flat band. As we discuss further below, the prediction of a flat band with $d$-wave character around $\vec{k}=(\pi,\pi)$ is consistent with earlier exact diagonalization results reporting narrow $d$-wave and $p$-wave peaks at low energies \cite{Dagotto1990}, and further corroborated in \cite{Bohrdt2022pairingPrep}. In addition we predict a strongly dispersive band at low energies around $\vec{k}=\vec{0}$ which contributes only little spectral weight, however.

In Fig.~\ref{figSpectratJz3Bose} we show spectra for bosonic holes. As in the fermionic case, some regions have zero spectral weight owing to the nature of string-length one states, see Eq.~\eqref{eqSpectralWeightsl1}. In the bosonic case, as a consequence, no collapse of spectral weight at high energies can be observed. At low energies, the spectral weight is dominated by the lowest dispersive band. At higher energies, flat bands with $p$-, $d$- and $f$-wave character can be observed.

The spectrum for distinguishable holes, but with $t_1=t_2=t$ equal, can be obtained by the sum of the fermionic and bosonic spectrum. Qualitatively, this procedure matches our results in Fig.~\ref{figSpectratJz3Dist}; but note that the data in Fig.~\ref{figSpectratJz3Dist} is based on the less accurate asymmetric approximation, which leads to some quantitative deviations.

\subsubsection{Relation to pairs in a N\'eel-ordered state}
\label{subsecRelationNeelPairs}
To understand the relation of the two-hole spectra calculated in Figs.~\ref{figSpectratJz3Dist} - \ref{figSpectratJz3Bose} to studies of pairs in the $t-J_{(z)}$ or Hubbard models, see e.g. \cite{Dagotto1990}, note that we have so far worked in an effective parton basis. Namely, the spin background was assumed to define a featureless vacuum state and our starting point was the two-hole string basis $\ket{\vec{x}_1,\Sigma}$ defined in Eq.~\eqref{eqTwoHoleStringBasis}. Since a proper two-hole spectral function should connect the undoped state (our vacuum) to the paired eigenstates, the structure of the former matters.

To understand the effect of the spin background, we consider starting from a classical N\'eel state $\ket{\rm N}$, i.e. one of the symmetry-broken ground states of the 2D Ising model. We will assume that $\uparrow$ ($\downarrow$) spins occupy the A (B) sublattice and have bosonic or fermionic statistics characterized by $\mu$. This state breaks the lattice-translational symmetry, which leads to momenta well-defined only within the $2$-site magnetic Brillouin zone (MBZ). Nevertheless, within the effective string model with its one-site unit-cell, Umklapp scattering from outside to within the MBZ is not possible. Hence the two-hole bandstructure within the MBZ is obtained by simply folding the full dispersion into the MBZ.

In this process the rotational quantum numbers remain unchanged. As a result we can identify the following string-model states $\ket{\vec{k},m_4}$ for any statistics $\mu$
\begin{equation}
    \ket{\vec{k}=(\pi,\pi), m_4} \equiv 2 \ket{\vec{k}=0, m_4}_{\rm N}
    \label{eqNeelFoldingMBZ}
\end{equation}
with states $\ket{ \vec{k}, m_4 }_{\rm N}$ in the N\'eel background (proof below). I.e. for fermions, where only $m_4=0,2$ states exist at $\vec{k}=(\pi,\pi)$, we obtain the corresponding $s$- and $d$-wave peaks around $\vec{k} = 0$ in the MBZ; for bosons, where only $m_4=1,3$ states exist at $\vec{k}=(\pi,\pi)$, we obtain the corresponding $p$- and $f$-wave peaks around $\vec{k} = 0$ in the MBZ.

To show Eq.~\eqref{eqNeelFoldingMBZ}, we note that the two-hole string state corresponds to
\begin{equation}
    \ket{\vec{k},m_4} = \frac{1}{\sqrt{V}} \sum_{\vec{j}} e^{i \vec{k} \cdot \vec{j}}  \hat{h}^\dagger_{\vec{j}} \sum_{\vec{i}: \langle \vec{i}, \vec{j} \rangle} e^{i \frac{\pi}{2} m_4 \nu_{\langle \vec{i},\vec{j} \rangle}} \hat{h}^\dagger_{\vec{i}} \ket{\rm N},
\end{equation}
where $\frac{\pi}{2} \nu_{\langle \vec{i},\vec{j} \rangle} = {\rm arg}(\vec{i}-\vec{j})$ is the angle of $\vec{i}$ relative to $\vec{j}$. The hole creation operator acts on the N\'eel state as
\begin{flalign*}
    \hat{h}^\dagger_{\vec{j}} \ket{{\rm N}} &= \c_{\vec{j},\uparrow} \ket{{\rm N}} \qquad \vec{j} \in {\rm A}, \\
    \hat{h}^\dagger_{\vec{j}} \ket{{\rm N}} &= \c_{\vec{j},\downarrow} \ket{{\rm N}} \qquad \vec{j} \in {\rm B},
\end{flalign*}
and the state $\ket{\vec{k},m_4}$ becomes
\begin{multline}
    \ket{\vec{k},m_4} =  \frac{1}{\sqrt{V}} \sum_{\vec{j} \in {\rm A}} e^{i \vec{k} \cdot \vec{j}}  \c_{\vec{j},\uparrow} \sum_{\vec{i}: \langle \vec{i}, \vec{j} \rangle} e^{i \frac{\pi}{2} m_4 \nu_{\langle \vec{i},\vec{j} \rangle}} \c_{\vec{i},\downarrow} \ket{\rm N}+\\
    + \frac{\mu}{\sqrt{V}} \sum_{\vec{i} \in {\rm A}} \c_{\vec{i},\uparrow} \sum_{\vec{j}: \langle \vec{i}, \vec{j} \rangle} e^{i \frac{\pi}{2} m_4 (\nu_{\langle \vec{j},\vec{i} \rangle} + 2)} \c_{\vec{j},\downarrow} e^{i \vec{k} \cdot \vec{j}} \ket{\rm N}.
\end{multline}
In the second line we exchanged the order of the operators, which yields the exchange sign $\mu$, and used that $\nu_{\langle \vec{i},\vec{j} \rangle} = \nu_{\langle \vec{j},\vec{i} \rangle} + 2 ~{\rm mod} ~4$. At $\vec{k} = (\pi,\pi) \equiv \vec{\pi}$ we obtain
\begin{equation}
    \ket{\vec{\pi},m_4} = \frac{1-\mu e^{i \pi m_4}}{\sqrt{V}} \sum_{\vec{j}}  \c_{\vec{j},\uparrow} \sum_{\vec{i}: \langle \vec{i}, \vec{j} \rangle} e^{i \frac{\pi}{2} m_4 \nu_{\langle \vec{i},\vec{j} \rangle}} \c_{\vec{i},\downarrow} \ket{\rm N}
\end{equation}
which directly yields the result in  Eq.~\eqref{eqNeelFoldingMBZ}.

\subsection{String-based pairing mechanism}
Finally we discuss implications of our results for the binding energies $E_{\rm bdg}$ defined in Eq.~\eqref{eqEbdgDef}. Within the effective string model introduced in this article, we calculate the one-hole spinon-chargon energy $E_{1}$ and compare it to the energy $E_{2}$ of two holes bound by a string. Throughout we assume that the underlying Hamiltonian is of $t-J_z$ type, with pure Ising interactions in the background. As in the two-hole case we ignore (Trugman-) loop effects \cite{Trugman1988} or self-interactions of the strings \cite{Wrzosek2021}.

For more complicated microscopic models, such as the $SU(2)$-invariant $t-J$ model, corrections to the binding energy must be expected. While our quantitative predictions in this case are of limited use, they nevertheless allow us to identify relevant competing processes that tend to support or prevent pairs of holes from forming a tightly bound state.

\subsubsection{Binding energy in $t-J_z$ model}
In Fig.~\ref{figBindingEnergy} we show the negative binding energy $-E_{\rm bdg}$ in units of $t$, plotted as a function of $(J_z/t)^{2/3}$. When $-E_{\rm bdg}<0$, the bound state of two holes connected by a string is energetically favorable and we predict pairing. This is the case for bosonic, distinguishable and fermionic holes when $J_z$ is sufficiently large. Indeed, in the $t-J_z$ model, the asymptotic binding energy in the tight-binding limit $J_z \gg t$ is $E_{\rm bdg} = J_z/2$ (see discussion below). Within our approximations this value is closely approached for fermionic holes already at $J_z=t$, where bosonic and distinguishable holes are predicted to be significantly more strongly bound. As shown in Appendix~\ref{SuppMatMoreResults}, see Fig.~\ref{figEbndgLin}, all hole-types approach the asymptotic value $J_z/2$ when $J_z / t\gtrsim 10$.

Our effective theory clearly shows a significant increase of the binding energy for bosonic or distinguishable holes, as compared to the case with fermionic statistics. We attribute this to the frustrating effect of fermionic minus signs identified first by Trugman \cite{Trugman1988}, which are picked up when two indistinguishable fermionic holes are exchanged and the string is completely reversed.

Bosonic and distinguishable holes have identical binding energies, which is correctly predicted by the asymmetric approximation. This is due to the bosonic nature of bound states within the asymmetric approximation at the dispersion minimum around $\vec{k}=0$ which we discussed above. Because all rotational quantum numbers $\vec{m}$ are conserved within the linear string approximation introduced in Sec.~\ref{subsubsecLinStringThy}, we can solve the case with distinguishable holes at $\vec{k}=0$ for much larger cut-offs $\ell_{\rm max}$ within the asymmetric approximation (which is exact in this case) and test how strongly $E_{\rm bdg}$ depends on the maximum string length $\ell_{\rm max}$. To this end in Fig.~\ref{figBindingEnergy} we compare results formally including all $m_4$ and $m_3^{(1)}$ sectors with variable $\ell_{\rm max}$ up to $13$ to exact results at $\vec{m}=0$ with $\ell_{\rm max}=10^3$. We find that for $(J_z/t)^{2/3} > 0.2$, corresponding to $J_z/t > 0.09$ the finite-string length results with $\ell_{\rm max}=13$ are reasonably well converged.

For fermionic holes in the strong coupling regime $J_z \ll t$, we find that the binding energy approaches zero within our effective theory. In this regime larger values of $\ell_{\rm max}$ would be required to reach convergence, and thus we cannot draw a final conclusion whether holes are bound or unbound for $(J_z/t)^{2/3} > 0.2$ (i.e. $J_z/t > 0.09$).

\begin{figure}[t!]
\centering
\epsfig{file=FIG11.pdf, width=0.47\textwidth}
\caption{Binding energy from the effective theory, predicted for the $t-J_z$ model. We used the truncated basis including all $m_4=0,..,3$ and $m_3^{(1)}=0,1,2$ sectors, and string lengths up to $\ell_{\rm max}$ (as indicated in the plot). The string potential was calculated for an Ising background. Around $\vec{k}=0$ there is no coupling to rotational states for the indistinguishable holes and a linear string theory allows to work with much longer maximum string lenghts $\ell_{\rm max}=10^3$.
}
\label{figBindingEnergy}
\end{figure}

\begin{figure}[t!]
\centering
\epsfig{file=FIG12.pdf, width=0.49\textwidth}
\caption{Discussion of the underlying pairing mechanisms. (a) In the tight-binding limit $J_{(z)} \gg t$, kinetic contributions can be ignored. In this case the distortion of the N\'eel background can be minimized while also gaining maximal energy from the nearest-neighbor attraction term in the $t-J_{(z)}$ model. (b) In the strong coupling limit the hole's kinetic energy dominates, leading to a powerful pairing mechanism when one hole retraces the string of the other. This mechanism is supported by an effective spinon-chargon (sc) repulsion with a geometric origin \cite{Bohrdt2020} and suppressed by the hard-core chargon-chargon (cc) repulsion.
}
\label{figPairingMech}
\end{figure}

\subsubsection{Pairing mechanisms for holes}
For pairing to be energetically favorable one requires an intricate balance of different microscopic effects influencing the one- and two-hole ground state energies $E_{1,2}$ entering the expression for the binding energy \eqref{eqEbdgDef}. To shed more light on the underlying binding mechanism, we highlight various microscopic processes that either support or prevent pairing.

\emph{Tight-binding limit: $J_{(z)} \gg t$.--} In this case, hole motion can be ignored and the energetically most favorable locations of the holes can be determined, see Fig.~\ref{figPairingMech} (a). On one hand, the nearest-neighbor attraction $-J_{(z)}/4 \hat{n}_{\vec{i}} \hat{n}_{\vec{j}}$ in the microscopic models \eqref{eqHtJz}, \eqref{eqDefHtJ} favors tightly bound hole pairs. In addition, each hole breaks up antiferromagnetic bonds which contributed an energy $J_{(z)} C_1$ in the undoped ground state. Ignoring any back-action of the localized holes on their spin environment, one finds that neighboring holes feel an attractive binding energy $J_{(z)} (1/4 - C_1)$. In the Ising case this result is exact, yielding $E_{\rm bdg} = J_z / 2$ as claimed above.

\emph{Strong coupling limit: $J_z \ll t$.--} In this case, the kinetic energy of the holes plays an important role. As long as string formation is energetically favorable and the Nagaoka effect is suppressed \cite{White2001}, the leading order kinetic energy per hole is $E_{\rm kin}^{(0)} = - 2 t_{\rm eff}$. Here $t_{\rm eff} = \sqrt{z-1} t$ is the effective hopping constant between string states on the Bethe lattice, and $z$ is the coordination number of the underlying physical lattice \cite{Bulaevskii1968,Grusdt2018PRX}. This asymptotic zero-point contribution to the energy per hole is identical for spinon-chargon and chargon-chargon pairs, hence canceling exactly in the binding energy $E_{\rm bdg}$ \cite{Bohrdt2021pairing}.

Since the string potential is assumed to be linear in the string length within our theory, see Sec.~\ref{subsubsecLinStringThy}, the next-to-leading order contribution to the energy $E_n$ takes the universal form \cite{Bohrdt2021pairing}
\begin{equation}
 	E_n = - 2 t_{\rm eff} + \alpha_n (n t)^{1/3} J_{(z)}^{2/3} + \mathcal{O}(J_{(z)})
	\label{eqEnSC}
\end{equation}
for $n=1,2$-holes states. The factor of $n=2$ renormalizing the tunneling in the two-hole state reflects the reduced mass $m_{\rm red} = m/2 = 1/(2t)$ describing the relative motion of the hole pair attached to the string. This result can be directly obtained from Eq.~\eqref{eqHeffDist}; or see Ref.~\cite{Bohrdt2021pairing} for a simplified derivation. The pre-factors $\alpha_n > 0$ are non-universal constants, generally depending on $n$, which are determined by the details of the string potential.

From Eq.~\eqref{eqEnSC} we obtain a powerful binding mechanism if we assume $\alpha_1 \approx \alpha_2 = \alpha$. In this case,
\begin{equation}
 E_{\rm bdg} =  \left(  2 - 2^{1/3} \right) \alpha ~ t^{1/3} J_{(z)}^{2/3} > 0
 \label{eqEbdgScaling}
\end{equation}
asymptotically when $t \gg J_z$, see also Ref.~\cite{Bohrdt2021pairing}. This binding energy scales with a non-trivial power of $t^{1/3}$ and thus can easily exceed $J_{(z)}$ deep in the strong coupling limit. Indeed, for bosonic and distinguishable holes in Fig.~\ref{figBindingEnergy} we confirm the scaling predicted by Eq.~\eqref{eqEbdgScaling} for $J_z/t \to 0$, where $E_{\rm bdg} / t \propto (J_z/t)^{2/3}$ as expected. Plotted on a linear scale over $J_z/t$, see Fig.~\ref{figEbndgLin} in Appendix~\ref{SuppMatMoreResults}, we find a clear curvature of $E_{\rm bdg}/t$ deep in the strong coupling regime.

In general, the coefficients $\alpha_1$ and $\alpha_2$ are not identical. For example, as discussed above the fermionic statistics lead to an effective repulsion between holes, which can cause an increase of $\alpha_2$. Indeed, for fermionic holes at strong coupling due to finite-size effects we cannot identify a clear asymptotic behavior of the binding energy in Fig.~\ref{figBindingEnergy}. Over a significant range of values $J_z/t$ our results are consistent with $E_{\rm bdg} \propto J_z$, see also Fig.~\ref{figEbndgLin} in Appendix~\ref{SuppMatMoreResults}, which corresponds to $\alpha_1 / \alpha_2 \approx 2^{-2/3}$. In the following we discuss two further effects which influence the coefficients $\alpha_n$, as summarized in Fig.~\ref{figPairingMech} (b).

\emph{Hard-core repulsion of holes.--} A single hole forming a spinon-chargon bound state can realize the zero-length $\ell_\Sigma=0$ string state, corresponding to a spinon and a chargon on the same lattice site. In our effective description of two holes, we explicitly excluded such states to account for the hard-core nature of the holes. This effective chargon-chargon repulsion generally leads to a larger value of $\alpha_2 > \alpha_1$, suppressing the tendency to pairing.

In Ref.~\cite{Bohrdt2021pairing} we proposed another microscopic model with two layers, which allows to realize \cite{Hirthe2022arXiv} distinguishable hole pairs with opposite layer indices. In this system, the hard-core constraint of the holes is effectively removed and one can realize $\alpha_1=\alpha_2$. By comparison to numerical DMRG simulations we confirmed the strong pairing expected by Eq.~\eqref{eqEbdgScaling} in that model due to the hole's gain in kinetic energy \cite{Bohrdt2021pairing}. This demonstrates the importance of finding other mechanisms which allow to overcome the strong repulsion of hard-core holes, either by increasing $\alpha_1$ or decreasing $\alpha_2$.

\emph{Spinon-chargon repulsion in dimensions $d>1$.--} One such mechanism is the geometric spinon-chargon repulsion at strong coupling $J_{(z)} \ll t$ described in Ref.~\cite{Bohrdt2020}. This effect is due to a decreased zero-point kinetic energy around the string-length $\ell_\Sigma=0$ state, which leads to a localized repulsive interaction of strength $\propto t$. This in turn causes $\alpha_1$ to increase and approach $\alpha_2$, mimicking the hard-core repulsion of two holes and thus supporting pairing in the strong-coupling regime.

Quantitatively the spinon-chargon repulsion can be best understood by considering the ro-vibrational ground state with $\vec{m}=0$. As shown in Sec.~\ref{secEffStringModel} (see also Appendix \ref{SuppMatToyModel}), the latter can be described by a hopping problem on a semi-infinite lattice $\ell_\Sigma=0,1,2,...$ with an effective hopping strength $t_{\rm eff} = \sqrt{z-1} t$ in the bulk, and $t_0 = \sqrt{z} t$ between $\ket{\ell_\Sigma=0}$ and $\ket{\ell_\Sigma=1}$. To capture edge effects around $\ell_\Sigma=0$, this problem can be mapped to the even-parity sector of a hopping problem on an infinite one-dimensional lattice $\ell=...,-2,-1,0,1,2,...$ \cite{Bohrdt2020}, with tunneling strength $t_1 = \sqrt{z/2} t$ between $\ket{\ell= \pm 1}$ and $\ket{\ell=0}$ and $t_{\rm eff}$ otherwise. When $t \gg J_{(z)}$ this yields an effective repulsive interaction around $\ell=0$ with strength
\begin{equation}
 g_{\rm sc}(z) = \l \sqrt{z-1} - \sqrt{z/2} \r t,
\end{equation}
with $g_{\rm sc}(z) > 0$ for $z > 2$, i.e. in dimensions $d\geq 2$. In one dimension, the effect is absent, indicating a tendency to avoid pairing in $d=1$.

\emph{Spinon dynamics.--} Finally we note that spinon dynamics can also contribute to the energy of the one-hole spinon-chargon state. In the $t-J_z$ case the effective spinon hopping is due to Trugman loops \cite{Trugman1988} which lead to a negligibly small spinon kinetic energy a small fraction of $J_z$ \cite{Chernyshev1999,Grusdt2018PRX} for any ratio $t/J_z$. In the $t-J$ case, the spinon hopping can lead to a further reduction of the spinon-chargon energy on the order $J$, which provides another mechanism suppressing pairing for experimentally realistic ratios of $t/J$.

\section{Summary and Outlook}
\label{secSummary}
In summary, we have studied an effective model of a pair of mobile dopants bound together by a strongly confining string. Our work is motivated by the physics of holes moving in an antiferromagnet, which is believed to be at the heart of many strongly-correlated electron systems. While the model applies most directly to mobile dopants in an Ising antiferromagnet, we believe it also has relevance upon including spin-flip terms or in entirely different settings such as in a strongly confined regime of lattice gauge theories with dynamical matter where the strings correspond to gauge fields. We also studied the effect of exchange statistics of the charge carriers, which allowed us to extend our model to Hubbard-Mott excitons with distinguishable dopants (a doublon bound to a hole) or a theoretical scenario with bosonic spins featuring antiferromagnetic interactions.

To study the nature of the bound states predicted by the effective string model, we analyzed their ro-vibrational excitation spectra with full momentum resolution. This allowed us not only to reveal their binding energies, but more importantly gives access to the pair dispersion relations. Of the latter we revealed two types of bands: One, strongly dispersive bands, where one dopant retraces the string of the other giving rise to high-mobility of strongly bound pairs, even in situations where an isolated single dopant would feature a strongly renormalized mass. Second, flat bands, where the dopants still form a bound state with a strongly fluctuating average distance, but where destructive quantum interference effects suppress any center-of-mass motion of the pair. Such flat-band states of pairs have previously been mentioned \cite{Shraiman1988a} but, to our knowledge, never been analyzed in greater detail.

The main results of our work can be summarized as follows. For fermionic holes doped into a N\'eel state, most directly related to the problem of high-$T_c$ superconductivity in cuprates, we reveal a low-lying flat band corresponding to a pair with $d$-wave symmetry. Only in the direct vicinity of the $\Gamma$-point in the Brillouin zone we found a strongly dispersive paired state with $s$-wave symmetry. While the fate of these states upon including fully $SU(2)$-invariant spin-exchange terms remains to be clarified, the prospect of forming flat, or nearly-flat $d$-wave pairs at low energies suggests that they may be relevant for understanding competing phases featuring charge localization, such as stripe \cite{Kivelson2003} and pair density wave \cite{Agterberg2020} states. Thereby our analysis sheds new light on the question about the origin of superconductivity and possible connections to other systems, such as twisted bilayer graphene, believed to feature nearly flat bands.

For bosonic and distinguishable pairs of dopants, most relevant to Hubbard-Mott excitons, we found a stronger tendency towards pairing than in the fermionic case. We believe this is due to the frustrating effect of the fermionic exchange sign on the charge's kinetic energy \cite{Trugman1988}. In the spectra of bosonic and distinguishable dopants we revealed a strongly dispersive lowest-energy state with $s$-wave pairing symmetry, at energies significantly below the first flat-band bound state. We analyzed the underlying pairing mechanism in some detail, and predict a universal scaling of the binding energy $|E_{\rm bdg}| \propto t^{1/3} J^{2/3}$ for bosons and distinguishable charges in the strong-coupling regime. For fermions, in contrast, we revealed a reduced tendency towards pairing and due to finite-size effects the asymptotic behavior of $E_{\rm bdg}$ at strong couplings remains to be clarified.

Our work sets the stage for future extensions. In particular, it will be important to include the effects of quantum fluctuations, i.e. spin-flip terms, in the surrounding antiferromagnet. Starting from the effective string model developed here, the so-called generalized $1/S$-expansion technique \cite{Grusdt2018PRX} can be used to capture such effects as additional polaronic dressing of the bound states with magnons. For the strongly dispersive bands we found here, only relatively weak renormalization by magnon dressing can be expected to occur due to the large recoil energy associated with magnon emission from a pair. On the other hand, the flat band states with a localized center-of-mass are expected to be more strongly renormalized by magnons, and we anticipate that a weak dispersion can be induced. This picture is consistent with our recent DMRG results \cite{Bohrdt2022pairingPrep}.

Another noteworthy assumption we have made is to focus on tightly bound pairs of charges connected by a string. On the other hand, and as described in the text, individual charges can also form spinon-chargon bound states. By including magnon dressing of the latter, long-range effective interactions between these one-hole states can also be induced. These could lead to other types of bound states, with a character entirely different from the tightly-bound chargon-chargon pairs described in this article. Exploring such states, and their connection to the tightly-bound string states discussed here, will be a worthy future endeavor.

\section*{Acknowledgements}
We thank M. Hafezi, I. Bloch, M. Greiner, L. Vidmar, U. Schollwöck, T. Hilker, S. Hirthe and L. Homeier for fruitful discussions, and L. Homeier for critical reading of our manuscript. This research was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy -- EXC-2111 -- 390814868, by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programm (Grant Agreement no 948141  — ERC Starting Grant SimUcQuam), by the NSF through a grant for the Institute for Theoretical Atomic, Molecular, and Optical Physics at Harvard University and the Smithsonian Astrophysical Observatory. ED acknowledges support from the ARO grant number W911NF-20-1-0163.

\section*{References}