\section{Introduction}
Field theoretic approaches to quantum spin models in lattices, such as the Heisenberg antiferromagnet (AFM), provide very successful descriptions of these paradigmatic quantum many-body systems \cite{Haldane1983a} and have led to a thorough understanding of their various quantum phase transitions in different dimensions \cite{Sachdev2011}. Key to their success is the underlying hypothesis that the coarse-grained fields on long length-scales feature similar behavior as the microscopic local magnetic moments underlying the spin model. More formally, a simple renormalization-group (RG) procedure yielding the effective low-energy field theory does not change the particle-content of the analyzed fields. However, in dimensions larger than one and with mobile dopants included, this approach has not been able to explain the rich phase diagram of high-temperature superconductors so far.

In this article, we take a different perspective and explore emergent structures, at low- to intermediate energies, in a doped quantum spin chain. The zoo of constituents we find defies a naive field-theoretic description: we identify emergent parton structures of spinons and chargons, forming mesonic bound states with a rich spectrum of ro-vibrational internal excitations. Moreover, these mesons interact with collective magnon excitations in the surrounding spin system, which leads to polaronic dressing on the one hand and, more exotically, to long-lived meson-magnon bound states. In a phenomenological field-theoretic model, each of these constituents should be described by a separate quantum field, with mutual interactions between all of them. Describing how these new fields emerge at intermediate length- or energy-scales in a thorough RG procedure is a challenging task, even for the simple toy model we consider. Hence we focus on a microscopic description of the individual emergent bound states and analyze their characterizing properties, such as their dispersion relations, zero-point energies, and mutual interactions. To this end, we apply the powerful theoretical tools developed for the description of Bose polarons \cite{Chevy2010,Devreese2020,Grusdt2015Varenna,Rath2013,Shi2018}.

Concretely, we consider doped one-dimensional spin chains. When featuring ${\rm SU}(2)$ invariance, these systems display spin-charge separation \cite{Kim1996,Sing2003,Ogata1990,Hilker2017,Vijayan2020}: the collective excitations of the spin-chain are fractionalized spinons, which co-exist with free chargons. In this limit, non-trivial bound states of the constituents are absent \cite{Ogata1990} and bosonization techniques provide a powerful field-theoretic description of the doped system in terms of Luttinger liquids \cite{Giamarchi2004}. As we demonstrate, the situation changes drastically when a staggered magnetic field is included, breaking the ${\rm SU}(2)$ symmetry, see Fig.~\ref{Fig1} a): Now spinons and chargons are confined \cite{Giamarchi2004,Borla2020,Kebric2021}, the undoped ground state has gapped collective magnon excitations, and spin-charge separation breaks down. Despite this confinement, the situation is far from trivial: As we will show, doped holes in this system host a zoo of excitations reflecting their rich internal structure, and their interaction with gapped magnons can lead to even more complicated multi-parton bound states, see Fig.~\ref{Fig2}.

In several regards, our 1D model is motivated by the physics of mobile holes in a ${\rm SU}(2)$-invariant 2D Hubbard model. In contrast to the 1D case with ${\rm SU}(2)$ symmetry, the ground state of the two-dimensional (2D) Heisenberg model has long-range magnetic order and gapless spin-$1$ magnon excitations. This effect is mimicked by the external staggered magnetic field we consider in our model, which introduces magnetic order and leads to similar spin-$1$ magnon modes, although with a non-vanishing gap. There is strong evidence that a doped hole in the 2D AFM features a rich internal meson structure \cite{Beran1996}, with discrete vibrational \cite{Brunner2000,Mishchenko2001,Bohrdt2020PRB} and rotational excitations \cite{Grusdt2018,Bohrdt2021}. In our 1D model, we reveal similar structures and develop an effective strong-coupling description.

\begin{figure}
    \centering
    \includegraphics[scale=0.5]{FIG1.png}
    \caption{Particle zoo in a doped spin chain: We consider a doped mobile hole in a spin-chain subject to a staggered magnetic field. a) Ignoring transverse spin fluctuations, the staggered field leads to a confining force between spinons and chargons connected by a string of overturned spins (top row), which leads to meson formation. A similar situation is realized in a mixed-dimensional model, where a strong gradient prevents hole motion along the direction of the gradient (bottom row). b) Transverse spin couplings give rise to Holstein-Primakoff (HP) magnon fluctuations in the surrounding spin background. The latter interact with the spinon, which is surrounded by the strongly fluctuating but tightly bound chargon cloud. All constituents making up the zoo of excitations are summarized in c), where we also indicate the background Ising fields $\tau^z_j$, affected by the hole motion, around which we expand in the generalized $1/S$ approximation employed here.}
    \label{Fig1}
\end{figure}

Remarkably, the coupling of mesons to collective magnon excitations remains poorly understood in 1D and 2D, in particular around zero momentum where we show that the competition of mesons and magnons is most pronounced. In the present article we fill this gap and apply a powerful theoretical framework, the so-called generalized $1/S$ expansion \cite{Grusdt2018}, to describe the coupling of mesons to magnons in a systematic manner. As a result, we obtain an effective polaron Hamiltonian describing the dressing, or even binding, of a spinon-chargon meson with additional magnon excitations. Our work paves the way for similar studies in doped 2D Mott insulators, and may lead to a better understanding of the charge carriers and their interactions with magnons in underdoped copper oxides. In particular we expect that our formalism will be useful for understanding transport measurements involving magnetic polarons, such as the long-time spreading dynamics of a hole reported in~\cite{Hubig2020,Bohrdt2020Dyn,Ji2021,Nielsen2022arXiv}.

Experimentally, the model we consider can be realized with ultracold atoms in optical lattices, which have recently made significant advances in studying doped quantum magnets \cite{Bohrdt2021review}. On a mean-field level, our model moreover maps to a doped mixed-dimensional $t-J$ model \cite{Grusdt2018SciPost}, which can be realized by subjecting a Fermi-Hubbard system to a strong tilt along one of the lattice directions \cite{Hirthe2022arXiv}. Ultracold atom realizations allow to measure spectra \cite{Stewart2008,Feld2011,Bohrdt2018,Brown2019} like the ones we calculate here to identify the emergent zoo of excitations; moreover they can directly visualize string patterns \cite{Endres2011,Hilker2017,Chiu2019} or the dressing cloud of magnetic polarons in configuration-space \cite{Koepsell2019}, making them ideal platforms to explore the emergent structures we predict on intermediate length scales.

\section{Model and main results}
In this article, we study a simple but rich one-dimensional model of a doped AFM. Our starting point is an ${\rm SU}(2)$-invariant Heisenberg spin chain. An additional staggered magnetic field of strength $\pm h$ on alternating sites along the $z$-direction breaks the ${\rm SU}(2)$ symmetry, introduces long-range magnetic correlations, and leads to collective magnon excitations with a tunable gap controlled by $|h|$. To describe mobile holes doped into this model, we use a $t-J$ Hamiltonian:
\begin{eqnarray}\label{model}
    \nonumber \H = -t \sum_{j,\sigma} \P \l \hat{c}^\dagger_{j+1,\sigma}\c_{j,\sigma} + \text{h.c.} \r \P\\
    + J \sum_{j} \hat{\vec{S}}_{j+1} \cdot \hat{\vec{S}}_j - h \sum_j \l -1\r^{j} \hat{S}^z_j.
\end{eqnarray}
Since we will only consider a single doped hole in this article, we dropped the nearest-neighbor interaction $-J/4 ~ \hat{n}_{j+1} \hat{n}_j$ typically included in the $t-J$ model \cite{Auerbach1994}. A similar model, including phonons, has been studied in Ref.~\cite{Kogoj2014}.

\subsection{Lattice gauge Hamiltonian}
For later purposes, we find it convenient to write the Hamiltonian as a sum of two separate parts: (i) a $t-J_z$ part which conserves each individual spin in the so-called squeezed space \cite{Ogata1990,Kruis2004a,Hilker2017} obtained by removing holes from the chain:
\begin{multline}
    \H_{t-J_z} = -t \sum_{j,\sigma} \P \l \hat{c}^\dagger_{j+1,\sigma}\c_{j,\sigma} + \text{h.c.} \r \P\\
    + J_z \sum_{j} \hat{S}^z_{j+1} \hat{S}^z_j - h \sum_j \l -1\r^{j} \hat{S}^z_j.
    \label{EqHIsing}
\end{multline}
To keep our analytical formalism later on general, we introduced the coupling $J_z$, which is simply $J_z=J$ for the original model in Eq.~\eqref{model}.

Remarkably, the Hamiltonian in Eq.~\eqref{EqHIsing} is exactly equivalent to a $\mathbb{Z}_2$ lattice gauge theory (LGT), as shown in Refs.~\cite{Grusdt2020PRL,Kebric2021}. In this mapping, chargons (i.e. spin-less holes) and spinons (i.e. Ising domain walls) carry $\mathbb{Z}_2$ gauge charges and are connected by a $\mathbb{Z}_2$ electric string $\tau^x_{\langle i,j \rangle}$. The staggered field $\pm h$ leads to a term $h \sum_{\langle i,j \rangle} \hat{\tau}^x_{\langle i,j \rangle}$ in the $\mathbb{Z}_2$ gauge invariant Hamiltonian. The latter has been shown to cause spinon-chargon confinement for any infinitesimal $h \neq 0$ \cite{Borla2020,Kebric2021}. This $\mathbb{Z}_2$ LGT formalism forms the basis for our mesonic description of a doped hole.

In addition, the full Hamiltonian in Eq.~\eqref{model} includes (ii) transverse spin fluctuations,
\begin{equation}
    \H = \H_{t-J_z} + \H_{J_\perp},
\end{equation}
where we find it most convenient to write
\begin{equation}
    \H_{J_\perp} = \frac{J_\perp}{2} \sum_{j} \l \hat{S}^+_{j+1} \hat{S}^-_{j} + \text{h.c.} \r.
    \label{eqHJperp}
\end{equation}
Again we introduced the more general coupling strength $J_\perp$ in this term, although for our original model in Eq.~\eqref{model} $J_\perp = J$.
Later on, we will include such transverse spin fluctuations on top of a N\'eel ordered ground state distorted by the hole motion by introducing Holstein-Primakoff bosons (magnons), see Fig.~\ref{Fig1} b).

Finally, we note that in the limit $h / J_\perp \to \infty$, the transverse fluctuations $\H_{J_\perp}$ can always be treated perturbatively, independent of the ratios $J_z/J_\perp$ or $t/J_\perp$. To lowest order, only the $t-J_z$ part of the Hamiltonian, Eq.~\eqref{EqHIsing}, remains and it follows that the model has an emergent $\mathbb{Z}_2$ gauge structure for large values of $h$.

\subsection{Main results: particle zoo in the spin chain}
The separation of the Hamiltonian in two components lends a natural understanding of our results. Our main goal is to understand the ground and excited states of a mobile dopant in the spin chain. As described in detail below, we find that the $\mathbb{Z}_2$ gauge structure of the $t-J_z$ part of the Hamiltonian, or equivalently (in our model) the string-picture of magnetic polarons \cite{Bulaevski1968,Trugman1988,Grusdt2018}, introduces parton constituents, namely spinons and chargons, which are confined by the linear string tension generated by the staggered field $h$. The resulting mesonic spinon-chargon bound state has a rich internal structure constituted by inversion-even and inversion-odd vibrational modes of the $\mathbb{Z}_2$ electric string, or equivalently the string of overturned Ising spins, connecting the spinon and the chargon. We probe these states directly in spectra calculated by time-dependent matrix product states (td-MPS), see Sec.~\ref{secDMRGresults}, and compare to an effective strong-coupling description that we develop here, see Sec.~\ref{secStrongCplngThy}.

The transverse couplings introduced by $\H_{J_\perp}$ lead to vacuum fluctuations of magnons in the absence of a doped hole. This effect can be captured by a simple linear spin-wave expansion around the classical N\'eel state, which we achieve by a Holstein-Primakoff approximation. In the vicinity of the meson, the distortion of the N\'eel background caused by the spinon-chargon pair introduces additional couplings to magnons which give rise to additional rich physics: On one hand, they lead to polaronic dressing and weak mass renormalization of the meson around the dispersion minimum at momentum $k=\pi/2$. This is shown for the lowest-energy mesonic state (solid blue line) in Fig.~\ref{Fig2}.

More dramatically, the interactions with magnons can give rise to meson-magnon bound states. Since the magnon itself can be viewed as a bound state of two confined spinons, this state constitutes an emergent tetra-parton composite. As demonstrated in Fig.~\ref{Fig2}, for sufficiently large values of $h$ our effective model of the meson-magnon coupling predicts a low-lying meson-magnon bound state at relatively low energies below the meson-magnon scattering continuum. This should be contrasted with the higher excitation energies of ro-vibrational internal meson modes. We confirm our prediction of meson-magnon bound states in td-MPS calculations of one-hole spectra in a sector with total spin $S_z=3/2$, see Sec.~\ref{secDMRGresults}.

\begin{figure}
    \centering
    \includegraphics[width=0.5\textwidth]{FIG2.pdf}
    \caption{Polaronic bands in the presence of meson-magnon interactions at low energies: The overall ground state at momentum $k=\pi/2$ is realized by a weakly dressed meson (solid blue line). Before the broad meson-magnon continuum is reached at higher energies (wide red band), we predict a weakly dispersing meson-magnon bound state (dark red line), corresponding to a tetra-parton configuration. The black lines indicate the bare meson dispersion (solid) and edges of the meson-magnon continuum (dashed) in the absence of meson-magnon interactions, respectively. At higher energies (not shown) we find ro-vibrational internal meson excitations. Calculations were performed using the strong-coupling generalized $1/S$ approximation introduced in the text; we chose parameters $h=0.6 J$ and $t=5 J$.}
    \label{Fig2}
\end{figure}

Finally, meson-magnon interactions can have a pronounced effect on the quasiparticle dispersion of the dressed meson around momentum $k=0$. In this region of momentum space, the bare meson dispersion approaches the meson-magnon scattering continuum most closely, as indicated by the dashed and solid black lines in Fig.~\ref{Fig2}. Without meson-magnon interactions and for sufficiently weak fields $h\ll t,J$, we find that they can even cross, leading to a decaying bare meson state inside the meson-magnon scattering continuum. However, in Sec.~\ref{SubSecAvdMgnDcy} we analyze our effective meson-magnon Hamiltonian and find indications that meson-magnon interactions in the 1D chain are strong enough to avoid such quasiparticle decay \cite{Verresen2018}. Namely, the meson and magnon bands repel and an isolated quasiparticle band of the mesonic magnetic polaron survives even around $k=0$. This prediction is further supported by td-MPS simulations at small fields where the effect is most pronounced.

Methodologically, we deviate from the standard approach typically used to describe magnetic polaron formation in an AFM \cite{Kane1989,Sachdev1989,Martinez1991}. As mentioned above, we first take into account how the mobile hole distorts the N\'eel background with pure Ising interactions. This allows us to make a direct connection to the $\mathbb{Z}_2$ LGT and identify the parton content of the meson. Moreover, we can relatively easily capture the competition between the tunneling term $t$ and the linear string tension $\propto h$, to all orders in $t/h$. This is achieved within a strong-coupling theory. Next we introduce generalized Holstein-Primakoff bosons (loosely speaking, magnons) by expanding around the \emph{already distorted} N\'eel state (we refer to this approach as the generalized $1/S$ approximation \cite{Grusdt2018}). This yields additional couplings of the meson to the magnons; importantly, the strength of these couplings is only of order $J$, and a fraction of $t$ for some further corrections we identify. Hence, perturbative or simple variational approaches are sufficient to capture the additional meson-magnon interactions. This should be contrasted with the traditional $1/S$ approximation \cite{Kane1989,Sachdev1989,Martinez1991} where the hole hopping $t$ itself leads to magnon creation: as a result, the effective Hamiltonian is strongly coupled when $t>h$ and direct analytical insights are harder to obtain.

\subsection{Possible experimental realizations}
Experimentally, the model in Eq.~\eqref{model} we study can be realized in different ultracold atom setups. We propose to use ultracold fermionic Lithium or Potassium atoms which have very successfully explored the ${\rm SU}(2)$-invariant 2D Fermi-Hubbard model \cite{Bohrdt2021review}. The main obstacle in these systems is to implement the staggered magnetic field, which requires local addressability on the scale of an optical wavelength, see e.g.~\cite{Weitenberg2011}, and sizable magnetic moments in order to distinguish different spin states, in a regime close to an atomic Feshbach resonance to realize super-exchange couplings.

A first option is to use Potassium atoms in a quantum gas microscope \cite{Cheuk2015} which have a sizable magnetic moment \cite{Nichols2018}, allowing for a local modulation of the magnetic field. A second option is to work in a mixed-dimensional setting where tunneling is strongly suppressed by strong gradients along all but one lattice direction \cite{Grusdt2018SciPost}. Moreover, we assume that nearest-neighbor AFM Ising couplings between all spins are present, which dominate over the weak super-exchange couplings along the gradient directions. This can be realized in an optical lattice by adding Rydberg dressing \cite{GuardadoSanchez2017}. When doping only the central chain with one hole and keeping all neighboring chains at half filling, the surrounding spin chains can generate an effective staggered field term $\pm h$ if they are sufficiently cold. Here we assumed, in a mean-field spirit, that the wavefunctions of the different chains approximately factorize. Similarly, in mixed-dimensional settings with $SU(2)$ invariant spin-exchange interactions \cite{Grusdt2018SciPost,Hirthe2022arXiv} we expect a ground state with broken $SU(2)$ symmetry in qualitatively very similar physics.

Finally, we note that the 1D model in Eq.~\eqref{model} can be equally realized with bosons as long as one ensures to have AFM Heisenberg couplings between the spins \cite{Trotzky2008,Dimitrova2019}. The statistics of the dopants is irrelevant, as can be shown by a Jordan-Wigner transformation. Hence the model in Eq.~\eqref{model} can also be simulated in qubit arrays or digital quantum computers \cite{Arute2020}, without the need to incorporate fermionic statistics.

\section{Numerical DMRG spectra}
\label{secDMRGresults}
In this section we present our numerical results, largely based on td-MPS simulations \cite{Paeckel2019}, which support our main findings about the structure and interactions of doped holes in the 1D spin chain with a staggered field. We already compare our numerical results to predictions by the semi-analytical strong-coupling meson-magnon theory introduced in the subsequent sections. This theory provides a unified understanding of all our key numerical observations.

Detailed descriptions of the numerical td-MPS simulations we performed can be found in Refs.~\cite{Bohrdt2020PRB,Bohrdt2021}; our algorithm builds upon the earlier works \cite{Kjaell2013,Zaletel2015,Verresen2018spec}. To ensure proper convergence of the MPS calculations, we performed the same convergence checks, in time and bond-dimension, as described in \cite{Bohrdt2020PRB,Bohrdt2021}.

\subsection{Ground state: Dressed hole}
\label{subsecDMRGgs}
In Fig.~\ref{StandARPES} we start by showing the standard one-hole angle-resolved photoemission spectrum (ARPES), defined by
\begin{equation}
    S(k,\omega) = -\frac{1}{\pi} \text{Im}{\int_0^\infty dt \frac{e^{i (\omega-E_0) t}}{\sqrt{L}} \sum_j e^{i k j} \mathcal{G}_j(t)},
    \label{eqAkwDef}
\end{equation}
with the Green's function
\begin{equation}
    \mathcal{G}_j(t) = \sum_\sigma \bra{\Psi_0} \hat{c}^\dagger_{j,\sigma} e^{- i \H t} \c_{0,\sigma} \ket{\Psi_0}
    \label{eqGjStandard}
\end{equation}
where $\ket{\Psi_0}$ is the ground state with energy $E_0$ of $\H$ with zero holes.

\begin{figure}
    \centering

    \includegraphics[width=0.5\textwidth]{FIG3.pdf}
    \caption{The standard one-hole ARPES spectrum reveals a pronounced quasiparticle peak at the lowest energy. The dispersion minimum is located at $k=\pi/2$, as predicted by our semi-analytical theory (solid blue line). Here we consider $h=1.0J$ and $t=5J$; the color scale is in a.u..}
    \label{StandARPES}
\end{figure}

In the spectrum, we observe a pronounced quasiparticle peak at low energy which corresponds to the magnetic polaron. The comparison with our semi-analytical theory shows that it is located around the expected energy, and shows the same dispersion relation with a minimum at $k=\pi/2$. At higher energies the spectrum is relatively featureless for the considered value of $h/J=1.0$ in Fig.~\ref{StandARPES}. As we show next, additional features becomes visible for larger values of $h$.

\subsection{Ro-vibrational excitations: Mesonic states}
\label{subsecDMRGmeson}
Now we calculate a rotational variant of the ARPES spectrum, where spinon-chargon excitations with odd $(\xi=-1)$ and even ($\xi=+1$) inversion symmetry can be detected. It is defined as in Eq.~\eqref{eqAkwDef} but using the rotational Green's function \cite{Bohrdt2021}
\begin{equation}
    \mathcal{G}^{\rm rot}_{j,\xi}(t) = \sum_\sigma \bra{\Psi_0} \hat{c}^\dagger_{j,\sigma} \hat{X}_{j,\xi}^\dagger e^{- i \H t} \hat{X}_{0,\xi} \c_{0,\sigma} \ket{\Psi_0}
\end{equation}
where $\hat{X}_{j,\xi}=\sum_\sigma \hat{c}^\dagger_{j,\sigma} ( \c_{j+1,\sigma}+\xi \c_{j-1,\sigma})$ creates an additional excitation of the spinon-chargon pair.

In Fig.~\ref{RotARPES} we show our results for $h=4J$. In both parity sectors $\xi=\pm 1$ we observe pronounced vibrational peaks, which correspond to vibrational modes of the spinon-chargon string. The absence of even (odd) peaks in the odd (even) spectrum indicates that the parity $\xi$ is a good emergent quantum number at all momenta, not only at $k=0,\pi/2$ where the system is strictly inversion symmetric. This is a direct indication for the existence of an internal meson structure \cite{Bohrdt2021}.

In Fig.~\ref{RotARPES} we also compare the peak positions observed in td-MPS with predictions by our strong-coupling theory. The observed peaks in our full numerical spectra are in excellent agreement with our semi-analytical predictions. In the latter, for simplicity, we neglected corrections from magnon-dressing which are weak at large values of $h$. Nevertheless, note that significant charge fluctuations are present since we consider $t>h$ in the figure.

\begin{figure}
    \centering

    \includegraphics[width=0.47\textwidth]{FIG4a.pdf}\\
    \includegraphics[width=0.47\textwidth]{FIG4b.pdf}
    \caption{The rotational one-hole ARPES spectrum reveals a series of long-lived vibrational excitations with even ($\xi=+1$, top) and odd ($\xi=-1$, bottom) parity. We compare the td-MPS spectra with bare meson resonances calculated from our strong-coupling theory (gray solid lines: $\xi=+1$ even; gray dashed lines: $\xi=-1$ odd). Here we consider $h=4.0J$ and $t=5J$; the color scale is in a.u..}
    \label{RotARPES}
\end{figure}

\subsection{Meson-magnon bound states}
\label{ARPESBound}
Next we show that even more complex excitations can arise when the mesonic hole interacts with its spin environment. Specifically, the meson can form a stable bound state with a magnon excitation. To demonstrate the existence of such bound states, we first consider an even more involved type of spectral function. To obtain spectral weight in the sector with one hole and one extra magnon, we create an excitation with total spin $S^z=3/2$ by flipping a spin next to the hole. This corresponds to working with the meson-magnon Green's function
\begin{equation}
    \mathcal{G}^{\rm mes-mag}_j(t) = \sum_\sigma \bra{\Psi_0} \hat{c}^\dagger_{j,\sigma} ~\hat{\Xi}^\dagger_{j}~ e^{- i \H t} ~ \hat{\Xi}_{0}~ \c_{0,\sigma} \ket{\Psi_0},
\end{equation}
where $\hat{\Xi}_{j} = \sum_{\delta=\pm 1}\hat{c}^\dagger_{j+\delta,\overline{\sigma}} \c_{j+\delta,\sigma}$ flips an additional spin and $\overline{\uparrow}=\downarrow$ ($\overline{\downarrow}=\uparrow$).

\begin{figure}
    \centering

    \includegraphics[width=0.5\textwidth]{FIG5.pdf}
    \caption{The spin-flip one-hole ARPES spectrum probes the sector with total spin $S^z=3/2$. It reveals a long-lived weakly dispersing meson-magnon bound state. We compare the td-MPS spectrum with our semi-analytical prediction for the meson-magnon bound state (solid red line). Parameters are $h=1.0J$ and $t=5J$; the color scale is in a.u..}
    \label{StagFlipARPES}
\end{figure}

The resulting spin-flip one-hole ARPES spectrum is shown in Fig.~\ref{StagFlipARPES}. There we observe a low-lying pronounced quasiparticle peak featuring a weakly dispersing band. Comparison to our semi-analytical prediction in the one-hole plus one magnon sector yields good qualitative agreement up to a small overall energy shift a fraction of $J$. Hence we interpret the observed feature as a stable meson-magnon bound state.

To further analyze the robustness of the meson-magnon bound state, we need to check whether it lies energetically below the meson-magnon scattering continuum. This is the case for the bound state predicted by our semi-analytical theory: Indeed, in Fig.~\ref{Fig2} we observe an isolated one-magnon excited state (lowest red band) between the mesonic ground state (blue) and the meson-magnon continuum (filled red band). To test this scenario in our fully numerical DMRG simulations, we calculate the meson-magnon binding energy, which is defined as follows:
\begin{equation}
    E_{\rm mm} = (E_{{\rm 1h},3/2} - E_{{\rm 0h},0}) - (E_{{\rm 1h},1/2} + E_{{\rm 0h},1} - 2 E_{{\rm 0h},0}).
\end{equation}
Here $E_{n{\rm h},s}$ denotes the ground state energy in the sector with $n$ holes and total spin $S^z=s$. If $E_{\rm mm}<0$, the meson-magnon state is located below the scattering continuum and forms a stable bound state.

In Fig.~\ref{MesMagBndgEngy} we plot the numerically obtained meson-magnon binding energy for various field strengths $h/J$, at fixed $t/J=5$. We find consistently that  $E_{\rm mm}<0$ beyond a critical field strength $h_c = 0.3(1) J$, confirming the existence of a stable bound state as anticipated from the spin-flip ARPES spectrum. In the figure, we also compare our results to the semi-analytical theory (solid red) and an effective theory valid at large $h \gg J$ (see Appendix \ref{ApdxBoundStateLargeh}). For small $h$, our semi-analytical theory is in good agreement with the numerics. At larger values of $h$ we observe deviations, which can be attributed to some simplifying approximations we made, see Sec.~\ref{secSolSCthy} for a detailed discussion. The large-$h$ theory provides good qualitative agreement everywhere.

\begin{figure}
    \centering
    \includegraphics[width=0.5\textwidth]{FIG6.pdf}
    \caption{The meson-magnon binding energy evaluated from DMRG simulations of the ground state with and without an additional hole and magnon. We used an additional magnetic field of strength $h_{\rm edge}=2 J$ at the boundaries of the chain to avoid boundary effects.}
    \label{MesMagBndgEngy}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=0.5\textwidth]{FIG7.pdf}
    \caption{Avoided quasiparticle decay. a) The meson dispersion in the limit $t \gg h,J$ approaches its strong coupling shape $J_\perp^* \cos(2k)$. b) The magnon dispersion in the limit $h \ll J$ approaches the 1D spinon-continuum, with a lower edge at $(J \pi/2) |\sin(q)|$. c) Without interactions, the bare meson state ($m$) would enter the meson-magnon ($mm$) continuum when $t \gg J \gg h$. d) In the presence of sufficiently strong meson-magnon interactions, the mesonic quasiparticle band remains stable for all momenta.}
    \label{FigAvdQPdecay}
\end{figure}

\subsection{Avoided magnon decay in a weak field}
\label{secAvdMagDcySmry}
Finally, to study the effect of meson-magnon interactions around zero total system momentum $k=0$, we return to the standard ARPES spectrum, i.e. the Green's function in Eq.~\eqref{eqGjStandard}. However now we consider a parameter regime where $t \gg J > h$. Beyond a critical value $t>t_c(J,h)$ depending on $h$ and $J$, the bare spinon-chargon dispersion is predicted by our strong coupling theory to enter the magnon continuum in the absence of meson-magnon interactions. Let us begin by specifying what is meant by this.

We consider two states, both at the same total system momentum $k$. In the first -- the bare meson state -- all momentum is carried by the spinon-chargon pair. From the strong-coupling theory we predict its energy to be
\begin{equation}
    \varepsilon_{\rm sc}(k) = J_\perp^* \cos(2k),
\end{equation}
with a renormalized tunneling $J_\perp^*$, see Fig.~\ref{FigAvdQPdecay} a). In Sec.~\ref{subsecDMRGmeson} we already confirmed that the meson dispersion takes this general shape. The second, competing state we consider contains an additional magnon excitation with momentum $q$ and energy $\omega_q$, see Fig.~\ref{FigAvdQPdecay} b). To obtain the same total momentum $k$, the meson carries momentum $k-q$ and the total energy of the state is:
\begin{equation}
    E_{k,q}=\varepsilon_{\rm sc}(k-q) + \omega_q.
\end{equation}
This defines the meson-magnon continuum.

Since for $h \neq 0$ the magnon spectrum $\omega_q > 0$ is gapped, the lowest energy state at $k=\pi/2$ always corresponds to a single spinon-chargon pair. However, at $k=0$ the situation is much more interesting. In particular, the meson state with $k=0$ and energy $\varepsilon_{\rm m}=\varepsilon_{\rm sc}(k)=J_\perp^*$ is very competitive with the meson-magnon state at $q=\pi/2$ which has energy $\varepsilon_{\rm mm}=\omega_{\pi/2}-J_\perp^* = E_{k=0,q=\pi/2}$. Indeed, for very small values of $h \to 0$ we can show analytically using our strong-coupling theory that $J_\perp^* \to J$. Moreover, in this limit the magnon dispersion approaches the well-known spinon continuum which has a lower-edge at $\varepsilon_{\rm s}(q) = J \pi/2 |\sin(q)|$ \cite{Giamarchi2004}. Hence
\begin{equation}
    \varepsilon_{\rm mm} - \varepsilon_{\rm m} =  \frac{\pi}{2} J - 2 J = -0.429... J < 0,
\end{equation}
and we conclude that, in the absence of meson-magnon interactions, the meson-magnon state has lower energy; i.e. the bare meson enters the meson-magnon continuum in this limit, see Fig.~\ref{FigAvdQPdecay} c).

\begin{figure}[t!]
    \centering
    \includegraphics[width=0.5\textwidth]{FIG8.pdf}
    \caption{Avoided quasiparticle decay in the weak-field large-tunneling limit, seen in ARPES spectra obtained by td-MPS simulations. We observe a pronounced quasiparticle peak at the lowest energies for all momenta $k$, even around $k=0$ where the meson comes closest in energy to the meson-magnon continuum. We show td-MPS results for $t=5 J$ and $h=0.01 J$; the color scale is in a.u..}
    \label{FigAvoidedMagDcytdMPS}
\end{figure}

Before we proceed, we note that the same result is obtained when the meson plus two-magnon continuum is considered. This case becomes relevant if the meson can only couple to pairs of magnons. Making the same considerations as above does not change the outcome, as one can see by placing the second magnon in the $q=0$ state whose energy $\omega_{q=0} \to 0$ as $h \to 0$.

At first glance, the scenario we find appears reminiscent of supersonic polarons emitting Cherenkov phonons \cite{Seetharam2021}; But we completely ignored interactions between mesons and magnons so far, which can destroy the quasiparticle. However, recently it has been shown that sufficiently strong couplings of a quasiparticle to an excitation continuum can, on the contrary, stabilize the quasiparticle band and lead to a complete avoidance of quasiparticle decay \cite{Verresen2018}. This scenario is sketched in Fig.~\ref{FigAvdQPdecay} d). As we show below in Sec.~\ref{SubSecAvdMgnDcy}, our effective theory indicates that meson-magnon couplings in the doped spin chain Eq.~\eqref{model} become strong enough at long wavelengths to cause an avoided quasiparticle decay.

In Fig.~\ref{FigAvoidedMagDcytdMPS} we confirm this prediction by td-MPS simulations of the ARPES spectrum. In fact, all ARPES spectra we considered showed a clearly defined quasiparticle band at low energies, for all momenta. In Fig.~\ref{FigAvoidedMagDcytdMPS} we consider the most extreme regime where $t \gg J \gg h$ and our argument above predicts the meson quasiparticle band to enter the magnon continuum in the absence of interactions.
Specifically we assumed $t=5 J$ and $h=0.01 J$.\\

\section{Strong-coupling theory\\ of doped holes}
\label{secStrongCplngThy}
In this section we discuss the generalized $1/S$ expansion for the 1D $t-J$ model in a staggered field, Eq.~\eqref{model}. It combines a parton theory for individual holes doped into the spin system \cite{Grusdt2018} with linear spin wave theory (LST) above the magnetically ordered spin background. The main achievement of this approach is to combine advantages of both methods: We keep the clear physical picture afforded by the parton theory while including a back-action of partons on their spin-environment.

Concretely the idea of the method is to include small quantum fluctuations about a classical magnetic state of Ising spins in the lattice. In a key distinction from earlier approaches \cite{Kane1989,Sachdev1989,Martinez1991}, we allow the Ising configuration around which we expand to be displaced by the quantum motion of the doped hole. This method can be formalized by using a generalization of the $1/S$ expansion in the length $S$ of the considered spins \cite{Grusdt2018}. As usual, we send $S \to 1/2$ in the end to obtain predictions for the spin-$1/2$ model in Eq.~\eqref{model}.

\subsection{Generalized 1/S expansion}
We begin the discussion by using a Schwinger boson representation of the (fermionic) $t-J$ model. We introduce a spin-less fermionic chargon $\hat{h}^\dagger_j$ and a Schwinger boson $\hat{\beta}^\dagger_{j,\sigma}$ and write the original fermion operators as $\hat{c}^\dagger_{j,\sigma}=\hat{h}_j \hat{\beta}^\dagger_{j,\sigma}$ for $S=1/2$, see e.g. \cite{Auerbach1994}. For general values of spin $S$ of the underlying fermions, the physical Hilbert space we consider is realized by states satisfying
\begin{equation}\label{single-occ-constraint}
    \sum_{\sigma} \hat{\beta}^\dagger_{j,\sigma}\hat{\beta}_{j,\sigma}= 2S \l 1- \hat{h}^\dagger_j \hat{h}_j\r.
\end{equation}
This constraint ensures that there is either one vacancy and no spin, or one spin and no vacancy localized at the specified site $j$. We emphasize that for $S\neq 1/2$ this constraint -- which is at the heart of the generalized $1/S$ expansion -- is different from the condition $\sum_\sigma \hat{\beta}^\dagger_{j,\sigma}\hat{\beta}_{j,\sigma} + \hat{h}^\dagger_j \hat{h}_j = 2S$ more commonly used in the conventional $1/S$ expansions \cite{Kane1989}. Using the new constraint in Eq.~\eqref{single-occ-constraint}  has the advantage to treat spin and charge as mutually exclusive degrees of freedom per lattice site, but leads to a highly non-linear chargon hopping term correlated with the surrounding spins when the $t-J$ model is expressed in terms of $\hat{h}_j$ and $\hat{\beta}_{j,\sigma}$.

Our strategy to deal with the complicated constraint in Eq.~\eqref{single-occ-constraint} in the following is two-fold. First, on the level of the Ising part of the Hamiltonian Eq.~\eqref{EqHIsing} we can keep track of the full constraint Eq.~\eqref{single-occ-constraint} due to the classical nature of the Ising spins. Second, we make a Holstein-Primakoff-approximation around the Ising configuration to take into account transverse spin fluctuations from Eq.~\eqref{eqHJperp}.

\subsubsection{Zero Doping}
At zero doping, the ground state of the model in Eq.~\eqref{model} has long-range AFM correlations along the $z$-direction for any non-zero $h$, which breaks the ${\rm SU}(2)$ symmetry of the bare $t-J$ model explicitly. The corresponding low-energy magnon excitations can be described using linear spin-wave theory, which is equivalent to a first order expansion of the model in powers of $1/S$ \cite{Auerbach1994}. To lowest order in $1/S$ one obtains a classical configuration $\hat{S}^z_j=\tau_j^z/2$ where the Ising variables $\tau_j^z= \l-1 \r^{j}$ describe a N\'eel configuration.

Up to first order in $1/S$ the linear spin wave theory corresponds to the Holstein-Primakoff (HP) approximation, where the spin operators are represented as
\begin{eqnarray}\label{HPapprox}
    \label{HPapproxSz}
    \hat{S}^z_j= \tau_j^z\l S - \hat{a}^\dagger_j\a_j \r,\\
    \label{HPapproxSpm}
    \hat{S}^{\tau_j^z}= \sqrt{2S}\a_j, \quad \hat{S}^{-\tau_j^z}= \sqrt{2S}\hat{a}^\dagger_j.
\end{eqnarray}
Note that in principle this expansion can be performed for arbitrary configurations of the classical Ising field $\tau_j^z$. The bosonic operators $\a_j$ are related to the Schwinger bosons by $\hat{\beta}_{j,-\tau^z_j}=\a_j$ and $\hat{\beta}_{j,\tau^z_j} = \sqrt{S}$, i.e. bosons $\hat{\beta}_{j,\sigma}$ with $\sigma=\tau^z_j$ condense and to leading order fluctuations of the condensate fraction are ignored.

Magnon excitations in the undoped AFM are obtained by setting $\tau^z_j=\l -1 \r^{j} $ and inserting Eqs.~\eqref{HPapproxSz} and \eqref{HPapproxSpm}  in the Heisenberg Hamiltonian. This results in the well-known free spin-wave Hamiltonian,
\begin{equation}
    \H_{\text{mag}}^{(0)}= \sum_q \omega_q \hat{b}^\dagger_q \b_q
\end{equation}
where the sum is taken over lattice momenta $q\in [-\pi,\pi]$.\\
The spin-wave dispersion is given by
\begin{equation}
    \omega_q = J_z \sqrt{\l 1+ h/J_z \r² - \l J_\perp/J_z\r²\cos²{q} }
\end{equation}
where we allowed for anisotropic interactions $J_z$ ($J_\perp$) along $z$ ($xy$) direction in spin space; for our model in Eq.~\eqref{model}, $J_\perp=J_z=J$. The Bogoliubov operators $\b_q$ are related to the HP bosons $\a_j$ by a Fourier- and Bogoliubov transformation.

\subsubsection{Single Hole Doping - Hilbert Space}
To describe the properties of a single hole doped into a N\'eel state, we apply the generalized $1/S$ expansion outlined above and extrapolate our result to the case $S=1/2$ in the end. In order to include distortions of the N\'eel state by the chargon, we work with the constraint \eqref{single-occ-constraint} and promote the Ising field $\tau^z_{j}$, around which we perform the linear spin-wave expansion later on, to a dynamical field: i.e. $\tau^z_{j}$ depends explicitly on the instantaneous configuration of the spinon, string of displaced spins, and chargon. We will make this precise in the following, by constructing a complete set of low-energy basis states for the partons.

Before we proceed, we mention already that changes in the Ising fields will lead to corresponding changes in the magnon terms resulting from the HP approximation, see Eq.~\eqref{HPapproxSz}. The idea is to include HP bosons on all bonds of the lattice, as done in the zero doping case, and describe how each term in the Hamiltonian leads to parton and / or magnon processes in the effective model.

To leading order in $1/S$ we can ignore magnons completely and only take into account changes in the Ising fields $\tau^z_{j}$ induced by the chargon motion. This is equivalent to solving only the $t-J_z$ part, Eq.~\eqref{EqHIsing}, in our model, which we now do for one doped hole. To this end we construct a set of low-energy basis states with one hole: we start from a N\'eel state, i.e. $\tau^z_{j}=\left(-1\right)^{j}$, and remove a fermion from some lattice site $j^s$. In accordance with the spin representation \eqref{HPapproxSz} and the constraint \eqref{single-occ-constraint}, we enlarge the allowed values of $\tau^z$ and set $\tau^z_{j^s}=0$. Next, we construct all relevant basis states by applying the hopping part $\H_t$ in Eq.~\eqref{model}, still ignoring magnons. In 1D these states can be labeled by the position $j^s$, where the hole has been initially created, and the position $j^h$ reached by the chargon,
\begin{equation}\label{parton-Hilbert-space}
    \{\ket{j^h,j^s}\}, \quad j^h= j^s + \Sigma,\quad \Sigma\in \mathbb{Z}.
\end{equation}

When the chargon moves it displaces all spins along its path $\Sigma$ by one lattice site, which changes the Ising fields $\tau^z_j$ on the corresponding lattice sites. Thereby it distorts the N\'eel pattern. However, since $\hat{S}^z_j = S\tau^z_j$ in the absence of magnons, each displacement can be associated with a potential energy cost $\propto h$. Thus, to leading order in $1/S$, the problem is described by a single hole moving in a classical spin background $\tau^z_j$ with AFM Ising interactions in the staggered field $\pm h$.

The site $j^s$ where the hole has been initially created carries a surplus of spin and corresponds to a domain wall of two nearest neighbor aligned spins in the $t-J_z$ model. The domain wall can be interpreted as the charge-neutral spinon which carries a spin $\sigma$ opposite to the spin of the removed fermion. In addition to the already introduced fermionic chargon operator $\hat{h}_j$, we define a new bosonic spinon operator $\hat{s}_{j,\sigma}$, which will become useful when mapping the $t-J$ Hamiltonian to the parton basis, Eq.~\eqref{parton-Hilbert-space}. To this end we make the following identification for the parton basis,
\begin{equation}\label{2nd-parton-Hilbert-space}
    \ket{j^h,j^s}\equiv \sd_{j^s,\sigma}\hat{h}^\dagger_{j^h}\ket{0}.
\end{equation}

As described above, every spinon-chargon configuration given by Eq.~\eqref{2nd-parton-Hilbert-space} is uniquely related to a configuration of the Ising fields $\tau^z_j$: at the position of the chargon, $\tau^z_{j^h}=0$, and along the string $\Sigma$ of displaced spins connecting the chargon and spinon, the Ising fields have a reversed sign as compared to the original N\'eel state, $\tau^z_j= -(-1)^{j}$ for $j=j^s+ \delta \Sigma $ where $\delta \Sigma= 1,...,\Sigma$. The spinon corresponds to the domain wall located at the beginning of the string, formed by two aligned spins.

Since we are working in a subspace with only one spinon and chargon, we can omit the spin index $\sigma$ at the spinon operator. Once the hole is created in the spin chain, the spin of the spinon is specified by the sub-lattice index of the removed fermion, and within our approximation this sub-lattice index of the spinon cannot change.

Before we proceed to construct the effective parton Hamiltonian, we note that our construction above is equivalent to assuming that $\Sigma$ describes a $\mathbb{Z}_2$ electric string connecting a pair of $\mathbb{Z}_2$-charged spinon and chargon in a $\mathbb{Z}_2$ LGT; see Refs.~\cite{Grusdt2020PRL,Kebric2021} for a detailed discussion of this general mapping.

\subsubsection{Parton Hamiltonian}
Mapping our model \eqref{model} to the parton basis in Eq.~\eqref{parton-Hilbert-space} and introducing HP magnons yields an effective Hamiltonian of the form
\begin{equation}\label{partonHam}
    \H= -L \varepsilon_0 +\H_{\rm h}+\H_{\rm s}+\H_{\rm sh}+\H_{\text{mag}},
\end{equation}
where $\varepsilon_0= S^2 J_z+Sh $ denotes the classical Néel ground state energy per lattice site; the total number of sites in the system is denoted by $L$.

The next two terms describe free chargon and spinon terms, namely
\begin{equation}\label{freechargonH}
    \H_{\rm h}= \varepsilon_0^{\rm h}\sum_{j}\hat{h}^\dagger_j\hat{h}_j+ t \sum_{j} \l\hat{h}^\dagger_{j+1}\hat{h}_j+\text{h.c.} \r
\end{equation}
and
\begin{equation}\label{freespinonH}
    \H_{\rm s}= \varepsilon_0^{\rm s}\sum_{j}\sd_j\hat{s}_j+ \frac{J_\perp}{2} \sum_{j} \l \sd_{j+2}\hat{s}_j+\text{h.c.}\r.
\end{equation}
Here, $\varepsilon_0^{\rm h}= 2S^2 J_z+Sh$ and $\varepsilon_0^{\rm s}= S^2 J_z$ are the rest energies of the chargon and spinon in the $t-J_z$ model with the staggered field, respectively. We emphasize that we use second-quantized operators for the spinon and the chargon for notational convenience, keeping in mind that our derivation is valid for a single spinon-chargon pair.

An important caveat is that the spinon tunneling term, second term in Eq.~\eqref{freespinonH}, is only valid for the case $S=1/2$. It appears to first order in $1/S$ and its derivation is given in Appendix \ref{ApdxSpnDyn}.  The spinon dynamics results similarly as in the case of the 1D $t-J$ model without an external magnetic field where genuine spin-charge separation occurs \cite{Giamarchi2004}; in the present case the flip-flop terms $\propto J_\perp \hat{S}^+_{j+1}\hat{S}^-_j$ acting on bonds adjacent to the domain wall let the spinon move by two lattice sites. In our formalism, the action of these terms on other bonds away from the spinon creates HP bosons, as will be shown below. After this subsection we will consider the limit $S=1/2$, as the spinon tunneling is only valid in this case.

Next, the term $\H_{\rm sh}$ describes spinon-chargon interactions and consists of two terms specified below: a density-density interaction and a kinetic interaction,
\begin{equation}\label{spinon-chargon-int}
    \H_{\rm sh}=\hat{\mathcal{H}}_{\rm sh}^{\rm dd}+\hat{\mathcal{H}}_{ \rm sh}^{\text{kin}}.
\end{equation}
The first spinon-chargon interaction describes the linear confining potential stemming from the string of displaced spins in the $t-J_z$ Hamiltonian \eqref{EqHIsing},
\begin{equation}
    \hat{\mathcal{H}}_{\rm sh}^{dd}= \sum_{i,j} ~ \hat{h}^\dagger_i \hat{h}_i ~ \sd_j \hat{s}_j ~ V_{\rm sh}(|i-j|)
\end{equation}
with
\begin{equation}\label{linear-pot}
    V_{\rm sh}(\ell)= |h| \ell -S^2 J_z \delta_{\ell,0}, \quad \ell\ge0.
\end{equation}
The first term in \eqref{linear-pot} accounts for the spins residing on the energetically unfavorable sub-lattice along the string $\Sigma$ connecting the spinon and chargon. This leads to a string tension $\propto |h|$ which grows linearly with the length of the string $\ell=|\Sigma|$. The second term in Eq.~\eqref{linear-pot} describes a point-like attraction between the spinon and chargon. To understand its origin, note that the rest energies for the partons are different when they occupy the same lattice site compared to being on different lattice sites. The difference between these rest energies is accounted for by the point-like parton attraction.

The kinetic spinon-chargon term in Eq.~\eqref{spinon-chargon-int} describes how the spinon dynamics is constrained by the chargon. It only appears for the case $S=1/2$ as the spinon tunneling term,
\begin{equation}
    \hat{\mathcal{H}}_{\rm sh}^{\rm kin}= -\dfrac{J_\perp}{2} \sum_{j} \hat{h}^\dagger_{j+1}\hat{h}_{j+1}\left(\sd_{j+2}\hat{s}_j+\text{h.c.}\right).
    \label{eqHshKin}
\end{equation}
This term describes how the presence of a chargon disables the flip-flop term $J_\perp\hat{S}^+\hat{S}^-$ across it, which was originally assumed to give rise to spinon dynamics in Eq.~\eqref{freespinonH}. The minus sign in Eq.~\eqref{eqHshKin} subtracts this contribution.

The last term in the parton Hamiltonian \eqref{partonHam} summarizes all magnon contributions,
\begin{equation}\label{mag-contributions}
    \hat{\mathcal{H}}_{\text{mag}}= \sum_q \omega_q\l \hat{b}^\dagger_q\b_q+1/2\r + \hat{\mathcal{H}}_{\text{mag}}^{\text{int}}
\end{equation}
where the first term results from diagonalization of our model \eqref{model} in the case of zero doping. $\H_{\text{mag}}^{\text{int}}$ describes magnon-parton interactions -- to be specified below --  which involve HP magnons to quadratic order and have to be included due to the effects of the partons on the underlying N\'eel state. The full expressions are derived in detail in Appendix \ref{ApdxMagnonCont}.

\subsection{Strong coupling approximation}
A full solution of the spinon-chargon problem in the presence of magnons is not possible. To simplify our analytical formalism further, we now introduce a strong coupling (SC) theory of the meson formed by the spinon and the chargon. The meson, in turn, interacts with the bath of low-energy magnon excitations. We will demonstrate that magnons lead to polaronic dressing of the meson, or even to the formation of a meson-magnon bound state. The strong-coupling approach is valid for $t\gg J,h$.
It is based on the separation of time scales between the chargon motion and the spin degrees of freedom, i.e. the spinon dynamics and the magnon creation and annihilation processes. Namely, the chargon dynamics takes place on a much shorter time scale in this limit.

\subsubsection{Born-Oppenheimer Approximation - Meson Operator}
At strong couplings, $t\gg J,h$, the fast chargon can adiabatically follow the slow spinon dynamics and magnon creation and annihilation processes. Thus, the chargon motion can be treated first while fixing the spinon position. Moreover, we assume that the meson state is weakly affected by magnons and neglect spin fluctuations at first. This allows us to solve the parton part of the Hamiltonian, \eqref{partonHam}, i.e. excluding the magnon contribution $\H_{\rm mag}$, by making a Born-Oppenheimer ansatz for the meson state,
\begin{equation}\label{Born-Oppenheimer-ansatz}
    \ket{\psi_{\rm sh}^{(n,\xi)}(k)} = \dfrac{1}{\sqrt{L}} \sum_{j^s} e^{-i k j^s} \sd_{j^s}\ket{0} \otimes \ket{\psi_h^{(n\xi)}(j^s)}
\end{equation}

This wavefunction corresponds to plane waves where $k$ denotes the total momentum of the meson, which is carried by the heavy spinon. For a fixed spinon position $j^s$, the chargon wavefunction is given by $\ket{\psi_h^{(n\xi)}(j^s)}$. It only depends on the distance from the chargon to the spinon, and can be characterized by two quantum numbers $n \in \mathbb{Z}_{>0}$ and $\xi=\pm 1$. They correspond to vibrational and rotational states of the chargon. The derivation of the meson energies and states can be found in Appendix \ref{ApdxMesonStatesSC}.

The ansatz in Eq.~\eqref{Born-Oppenheimer-ansatz} leads to a meson dispersion of the form
\begin{equation}\label{meson-dispersion}
    \varepsilon_k^{(n\xi)} = E_h^{(n\xi)} + J_\perp^{(n\xi)} \cos(2k)
\end{equation}
which describes a tight-binding dispersion for the meson with a two-site tunneling term.
Here, the meson hopping amplitude $J_\perp^{(n\xi)}$ is related to the spin exchange coupling $J_\perp$ by a Franck-Condon factor,
\begin{equation}
    \label{FCoverlap}
    J_\perp^{(n\xi)}/J_\perp = \bra{\psi_h^{(n\xi)}(j^s+2)} 1- \hat{h}^\dagger_{j+1}\hat{h}_j \ket{\psi_h^{(n\xi)}(j^s)}.
\end{equation}
The energy offset $E_h^{(n\xi)}$ corresponds to the chargon eigenenergy characterized by the quantum numbers $n$ and $\xi$. As shown in Appendix \ref{ApdxMesonStatesSC}, for $t \gg h$ it scales as,
\begin{equation}
    E_h^{(n\xi)} \approx E_0 -2t + a_{\rm sh}^{(n\xi)} t^{1/3}h^{2/3} + \mathcal{O}(h),
\end{equation}
with numerical coefficients $a_{\rm sh}^{(n\xi)}$ related to the Airy-function \cite{Bulaevski1968}.

The center-of-mass momentum of the mesonic bound state is carried by the heavy spinon and the binding of the spinon to the fluctuating chargon leads to a renormalization of the hopping amplitude, $J^{(n\xi)}\le J$, see Fig.~\ref{Fig2}. Because the Franck-Condon factor $J^{(n\xi)}/J$ is independent of $k$ at strong couplings, the shape of the meson dispersion is identical to that of the spinon up to an overall rescaling.

Formally, we can define meson operators $\hat{f}^\dagger _{j^s,n\xi}$ to describe the bound state,
\begin{equation}
\label{mesonop}
    \hat{f}^\dagger _{j^s,n\xi} \ket{0} = \sd_{j^s}\ket{0}_s \otimes \ket{\psi_h^{(n\xi)}(j^s)}.
\end{equation}
This will help us to describe interactions between the meson and the bath of magnons in the next step. Note that in this description, the chargon quantum numbers $(n,\xi)$ take the role of band indices of the meson.

\begin{table*}[t]
\renewcommand{\arraystretch}{1.5}
    \centering
    \begin{tabular}{p{2.2cm}||p{5.5cm}|p{7.5cm}||p{1.2cm}}
        Theory & free variables & alternative representation & $\H$ \\
        \hline \hline
        $t-J$ model & $\c_{j,\sigma}$ (spin-$1/2$ fermions) & $\c_{j,\sigma} =  \hat{\beta}_{j,\sigma} \hat{h}^\dagger_{j}$ (Schwinger boson, slave fermion) & Eq.~\eqref{model}\\
        \hline
        parton theory &  $\hat{\tau}^z_j$ (Ising variables), $\a_j$ (HP magnon) & $\hat{s}_{j,\sigma}$ (spinon), $\hat{h}_j$ (chargon), $\b_q$ (Bogoliubov magnon)& Eq.~\eqref{partonHam} \\
        \hline
        polaron theory & $\hat{f}_{k,n\xi}$ meson, $\b_q$ (Bogoliubov magnon) & $\qquad \qquad \qquad$ -- & Eq.~\eqref{PolHam}
    \end{tabular}
    \caption{Overview of the different levels of approximation used in our semi-analytical theory, as described in the text.}
    \label{tabOverview}
\end{table*}

\subsubsection{Polaron Hamiltonian}
\label{SecPolHam}
Next we include magnon processes. As in the case of spinon dynamics, the time scales associated with magnon processes also correspond to longer times compared to the chargon motion. It is thus legitimate to extend the SC ansatz to the magnon contributions and treat magnon terms on a mean-field level. We assume that the meson is formed even in the presence of magnons, and neglect any back-action of the latter on the underlying SC meson wavefunction, such as a possible weak renormalization of the linear string tension.

Within the SC theory, the effective magnon Hamiltonian is obtained by averaging the parton Hamiltonian over the SC spinon-chargon wavefunction Eq.~\eqref{Born-Oppenheimer-ansatz}. Working in second quantized notation with the meson operators $\hat{f}^\dagger_{j^s,n\xi}$ we obtain:
\begin{equation}
    \label{PolHamProjection}
    \hat{\mathcal{H}}_{\rm pol}= \sum_{\substack{j^s,n\xi \\ j^{s\prime},n^\prime\xi^\prime}} \bra{\psi_{sh}^{(n^\prime\xi^\prime)}(j^{s\prime})}{\hat{\mathcal{H}}}\ket{\psi_{sh}^{(n\xi)}(j^s)} \hat{f}^\dagger_{j^{s\prime},n^\prime\xi^\prime}\hat{f}_{j^s,n\xi}.
\end{equation}
The resulting effective polaron Hamiltonian can be decomposed into different contributions,
\begin{equation}\label{PolHam}
    \hat{\mathcal{H}}_{\text{pol}}= \hat{\mathcal{H}}_{\text{mes}}^{(0)}+\hat{\mathcal{H}}_{\rm mag}^{(0)}+\hat{\mathcal{H}}_{\rm pol}^{\rm int}.
\end{equation}

The free meson and magnon Hamiltonians are
\begin{subequations}
\begin{eqnarray}
    \hat{\mathcal{H}}_{\text{mes}}^{(0)}&=& \sum_{k,n\xi}\varepsilon^{(n\xi)}_{k, \rm eff}\hat{f}^\dagger_{k,n\xi}\hat{f}_{k,n\xi}, \label{eqHmesFree}\\
    \hat{\mathcal{H}}_{\text{mag}}^{(0)}&=&\sum_q \omega_q \hat{b}^\dagger_q\b_q.
\end{eqnarray}
 \end{subequations}
Note that the free meson dispersion gets weakly renormalized by zero-point contributions of the magnons resulting from the Bogoliubov transformation. However, the analytic form of Eq.~\eqref{meson-dispersion} remains unchanged: $\epsilon^{(n\xi)}_{k} \to \epsilon^{(n\xi)}_{k, \rm eff} = E_{h,\rm eff}^{(n\xi)}+ J_{\perp, \rm eff}^{(n\xi)} \cos(2k)$. We do not include  full expressions for the small corrections here but refer the reader to Appendix \ref{ApdxMesMagInt} for details.

The meson-magnon interactions in Eq.~\eqref{PolHam} take the compact form
\begin{eqnarray}
    \label{PolInt}
    \nonumber\hat{\mathcal{H}}_{\rm pol}^{\rm int}=- \dfrac{1}{2L}\sum_{n\xi}\sum_{k,pq} \left[ \mathcal{V}^{(n\xi)}_{k,pq} ~ \hat{f}^\dagger_{k-p+q,n\xi}\hat{f}_{k,n\xi}~ \hat{b}^\dagger_{p}\b_q \right.\\
    \left. + \l \mathcal{W}^{(n\xi)}_{k,pq}~ \hat{f}^\dagger_{k-p-q,n\xi}\hat{f}_{k,n\xi}~ \hat{b}^\dagger_{p}\hat{b}^\dagger_q +\text{h.c.} \r \right],
\end{eqnarray}
where we ignore band-changing collisions, i.e. we only consider terms diagonal in $(n,\xi)$. This is justified in the SC limit by the large separation of energy scales.
The Bogoliubov operators $\b_q$ are related to the HP bosons $\a_j$ by a Fourier- and Bogoliubov transformation where the last transformation diagonalizes the free magnon terms. This leads to explicit expressions for the two couplings $\mathcal{V}^{(n\xi)}_{k,pq}$ and $\mathcal{W}^{(n\xi)}_{k,pq}$ describing normal and anomalous magnon terms, as derived in Appendix \ref{ApdxMesMagInt}.

This effective strong coupling Hamiltonian describes a polaron model of the mesonic impurity coupled to the bath of quantum spin fluctuations, represented by collective magnon excitations.

\subsection{Overview of the approach}
Above we described how the original $t-J$ model with a staggered field, Eq.~\eqref{model}, first leads to an effective parton model which subsequently maps to the simplified polaron Hamiltonian in Eq.~\eqref{PolHam}. In Tab.~\ref{tabOverview} we provide a summary of the involved fields and Hamiltonians appearing in the various stages of approximations. Next we will solve the effective polaron theory.

\section{Solutions of the strong-\\coupling theory}
\label{secSolSCthy}
In this section we apply and compare different methods to solve the effective polaron Hamiltonian from Eq.~\eqref{PolHam}. These are based on known analytical polaron techniques which have been successfully used to study Bose polaron problems in the past\cite{Chevy2010,Devreese2020,Grusdt2015Varenna,Rath2013,Shi2018}.

\subsection{LLP + Gaussian approach: Meson-magnon binding}
The above effective meson-magnon theory does not couple states where the total momentum of the meson plus magnon excitations is changed. Thus the total momentum is a conserved quantity which reflects the underlying translational invariance of the system. Below we make this conservation of the total momentum explicit by applying a Lee-Low-Pines (LLP) transformation \cite{Lee1953} to the Hamiltonian \eqref{PolHam}. To solve the resulting Hamiltonian $\tilde{\mathcal{H}} = \hat{U}^\dagger_{\text{LLP}}\hat{\mathcal{H}}\hat{U}_\text{LLP}$ in the LLP frame, we simplify it further by expanding to quadratic order in magnons. This allows us to solve it explicitly using multi-mode Gaussian states of magnons.

\begin{figure*}[t]
    \centering

    \includegraphics[width=19cm]{FIGwide1.pdf}
    \caption{Average densities of bare HP bosons $\langle \hat{a}^\dagger_j \a_j \rangle$ (red), Bogoliubov bosons  $\langle \hat{b}^\dagger_j \b_j \rangle$ (magnons, blue) and the chargon $\langle \hat{h}^\dagger_j \hat{h}_j \rangle$ (gray) for total momenta $k=0$, $k=\pi/4$ and $k=\pi/2$. The staggered magnetic field is chosen as $h=0.45J$ and the chargon tunneling as $t=5J$. At this particular magnetic field value our theoretical treatment shows that at a total momentum $k=0$ (left panel) the tetra-parton bound state of the meson with a magnon comes close in energy to the dressed meson state, while for $k=\pi/2$ (right panel) they are further apart form each other, see Fig.~\ref{Fig2}. The formation of the bound state is also visible here in the quasi-particle densities: While in the right panel the HP boson density is strongly suppressed around the spinon, the left panel shows that local spin flips start to accumulate around the meson.}
    \label{Densities}
\end{figure*}

\subsubsection{LLP Transformation}
The LLP transformation shifts the entire magnon state into the frame co-moving with the center-of-mass of the impurity -- in our case the mesonic bound state. It is represented by the unitary transformation,
\begin{equation}
\hat{U}_\text{LLP}= e^{-i \hat{X}_\text{mes}\hat{Q}_{b}},
\end{equation}
where we introduced the meson position operator, $\hat{X}_{\text{mes}}=\sum_{j^s,n\xi}j^s \hat{f}^\dagger_{j^s,n\xi}\hat{f}_{j^s,n\xi}$ and the total magnon momentum operator, $\hat{Q}_{b}= \sum_q q\hat{b}^\dagger_q\b_q$.

Now we apply the unitary LLP transformation to the effective polaron Hamiltonian \eqref{PolHam}. This is established by determining how the meson and magnon operators transform, namely
\begin{eqnarray}
    \hat{U}^\dagger_\text{LLP} \hat{f}^\dagger_{k,n\xi}\hat{U}_\text{LLP}= \hat{f}^\dagger_{k + \hat{Q}_b,n\xi} \\   \hat{U}^\dagger_\text{LLP}\hat{b}^\dagger_q\hat{U}_\text{LLP}= e^{iq\hat{X}_{\text{mes}}}\hat{b}^\dagger_q.
\end{eqnarray}
Insertion of these relations into $\hat{U}^\dagger_{\text{LLP}}\hat{\mathcal{H}}\hat{U}_\text{LLP}=\tilde{\mathcal{H}}$ yields a Hamiltonian which is block-diagonal in the total system momentum $K$. Thus we get a Hamiltonian of the form,
\begin{equation}
    \tilde{\mathcal{H}}= \sum_{K,n\xi}\hat{f}^\dagger_{K,n\xi}\hat{f}_{K,n\xi}\otimes \hat{\mathcal{H}}_a (K),
\end{equation}
where the term $\hat{\mathcal{H}}_a (K)$ depends only on magnon operators $\hat{b}^\dagger_j$. I.e. we eliminated the impurity degree of freedom from the problem. However, the transformed magnon Hamiltonian $\hat{\mathcal{H}}_a (K)$ now includes non-linearities in the magnon operators.

Specifically, the LLP transformation leads to a shift of the meson momenta by the total magnon momentum operator. By applying the LLP transformation, the following replacement occurs in the free meson term Eq.~\eqref{eqHmesFree},
\begin{equation}
    \hat{U}^\dagger_\text{LLP} \H_{\rm mes}^{(0)} \hat{U}_\text{LLP}= \sum_{K,n\xi} J_{\perp,\rm eff}^{(n\xi)} \cos(2K-2\hat{Q}_b) \hat{f}^\dagger_{K,n\xi} \hat{f}_{K,n\xi},
\end{equation}
i.e. effectively $\cos(2k)\to \cos(2K-2\hat{Q}_b)$. Expressing the cosine in terms of exponentials, we observe strong magnon non-linearities corresponding to the factors of the form $e^{-i2\hat{Q}_b}$ in the Hamiltonian.

\subsubsection{Linearization and Gaussian states}
To deal with the non-linearities encountered during the LLP transformation, we use the approximation
\begin{equation}
    e^{-i2\hat{Q}_{b}}\approx \openone + \sum_{q}\left(e^{-i2q}-1\right)\hat{b}^\dagger_q\b_q.
    \label{eqLinearization}
\end{equation}
This equation holds exactly within a subspace of no more than one magnon excitation. Hence it is similar in spirit to the HP approximation we made earlier, which also restricted us to consider low magnon densities only. Both approximations rely on the assumption that meson-magnon interactions are sufficiently weak to work at low excitation densities.

The resulting LLP Hamiltonian $\tilde{\mathcal{H}}$ following this truncation yields a quadratic Hamiltonian in the bosonic magnon operators $\hat{b}^\dagger_q$, making it exactly solvable via a multi-mode Bogoliubov transformation \cite{Xiao2009}. Within the LLP frame, our approach is Gaussian, in the sense that all non-Gaussian contributions are ignored. Nevertheless, we emphasize that the overall meson-magnon wavefunction includes non-Gaussian correlations, since we applied the non-Gaussian LLP transformation first.

Our ansatz is similar in spirit to the more general class of non-Gaussian variational states discussed in Ref.~\cite{Shi2020}. Indeed, a more sophisticated approach would be to avoid the linearization in Eq.~\eqref{eqLinearization} and fully solve the variational problem in the LLP frame, $\langle \tilde{\mathcal{H}}(K)\rangle \stackrel{!}{=} {\rm min}$, in the class of Gaussian states. Here we chose the simpler linearization method because it admits more direct analytical insights and immediately yields a full magnon excitation spectrum, but extensions beyond this simple limit constitute an interesting future direction.

In the following we solve the meson-magnon polaron problem as described, and discuss the properties of our solution in more detail.

\subsubsection{Magnon Distribution}
To gain better understanding of how magnons lead to polaronic dressing of the meson, we first plot the average number of bare HP bosons $n_a(j) = \langle \hat{a}^\dagger_j \a_j \rangle$ and Bogoliubov bosons $n_b(j) = \langle \hat{b}^\dagger_j \b_j \rangle$ in Fig. \ref{Densities}. Since we work in the LLP frame $j$ corresponds to the distance to the spinon, which is located at the core of the meson and thus at the origin of the LLP frame. Additionally, we plot the local distribution of the chargon in its ro-vibrational ground state, which indicates the extension of the mesonic bound state and the chargon cloud around the spinon. Far away from the meson, the number of HP bosons approaches a constant value (vacuum fluctuations) which is a consequence of the spin-flip terms $\hat{S}^+_{j+1} \hat{S}^-_j$. The asymptotic value of $n_a(|j|\gg 1)$ can be calculated straightforwardly from linear spin-wave theory.

In the ground state at $k=\pi/2$, Fig.~\ref{Densities} right, we observe a suppression of the number of HP bosons, i.e. fewer local spin flips, in the vicinity of the spinon. The suppression of $n_a(j)$ close to the spinon is a direct consequence of the formation of the geometric string. In the region where the chargon delocalizes, i.e. the spinon-chargon string fluctuates, the magnetization of the spin background is reduced and therefore we expect suppressed quantum fluctuations around the spinon. Note that this suppression of magnon fluctuations is dictated by the chargon distribution.

The substantial \emph{reduction} of quantum fluctuations around the spinon also provides an a-posteriori justification for the use of the lowest order Holstein-Primakoff approximation in Eqs.~\eqref{HPapproxSz}, \eqref{HPapproxSpm}. This is in contrast to the conventional $1/S$-expansion\cite{Kane1989,Sachdev1989,SchmittRink1988}, where non-linear terms in the magnons must be included to prevent excessive densities of excitations at strong couplings. Indeed, in the conventional $1/S$-expansion where magnons are defined relative to the undoped N\'eel state, a local \emph{enhancement} of bare HP magnon fluctuations is obtained around the mobile hole.

\subsubsection{Ground State Energy}
\begin{figure}
    \centering
    \includegraphics[width=0.42\textwidth]{FIG9.pdf}
    \caption{Ground state energy in our Gaussian LLP approach (solid red) compared to the numerical DMRG data (symbols).

    We also provide the ground state (G.S., dashed line) energy for the bare meson without including magnon couplings, to show the influence of spin fluctuations on the energy. The theory shows good agreement to the numerical data. We used  parameters $J_\perp=J_z=J=1$ and assumed that the meson carries a total momentum of $k=\pi/2$ in our calculations.}
    \label{GS}
\end{figure}

To acquire a better knowledge of how magnon excitations influence the dressed meson we show in Fig. \ref{GS} the ground state energy of one hole for various values of the staggered magnetic field $h/t$. We compare the ground state energy computed from our SC theory to direct numerical DMRG calculations of the ground state energy. The results show good agreement between our SC theory and the numerical data for all chosen magnetic field values. Note that we also display the ground state energy for the bare meson - neglecting magnon excitations - which follows from a pure parton theory for the doped hole. It shows good agreement to the numerics, especially for large values of $h$ where magnon dressing only leads to a weak renormalization of the meson dispersion, see Fig.~\ref{Fig2}. Thus the bare parton theory for the hole already provides a good description of the ground state properties of the meson.

\subsubsection{Polaron Spectrum}
\label{SecPolSpec}
As already elaborated, we have simplified the effective meson-magnon Hamiltonian~\eqref{PolHam} to a quadratic Hamiltonian in the Bogoliubov operators $\hat{b}^\dagger_q, \b_p$,
\begin{eqnarray}
    \label{H_aK}
    \nonumber \H_a (K) = E_{h,\rm eff}^{(n\xi)} + J_{\perp, \text{eff}}^{(n\xi)}\cos(2K)+ \sum_q \Omega^{(n\xi)}_{K,q} \hat{b}^\dagger_q \b_q\\
    - \dfrac{1}{2L} \sum_{pq} \left[\mathcal{V}^{(n\xi)}_{K,pq} \hat{b}^\dagger_p\b_q + \l \mathcal{W}^{(n\xi)}_{K,pq} \hat{b}^\dagger_p\hat{b}^\dagger_q + \text{h.c.}\r  \right]
\end{eqnarray}
by transforming the system into the LLP-frame and linearization -- corresponding to the frame co-moving with the meson -- which eliminates the meson degree of freedom. The total system momentum, $K$, is conserved in this frame and thus labels different total momentum blocks of the Hamiltonian.
The linearization~\eqref{eqLinearization} results in a renormalized free magnon dispersion at fixed total momentum $K$:
\begin{equation}
    \Omega_{K,q}^{(n\xi)} = \omega_q + J_{\perp,\text{eff}}^{(n\xi)} \l \cos(2K-2q) - \cos(2K) \r.
\end{equation}

To solve the quadratic Hamiltonian in Eq.~\eqref{H_aK}, we use a  multi-mode Bogoliubov transformation\cite{Xiao2009}. In the following we will sketch the idea and mention peculiarities arising with this multi-mode technique. It starts by introducing transformed Bogoliubov operators,
\begin{eqnarray}
    \label{MultiModeTrafo}
    \d_l = \sum_q \left[U_{l,q}\b_q -V_{l,q}\hat{b}^\dagger_q\right],\\
    \hat{d}^\dagger_l = \sum_q \left[U_{l,q}\hat{b}^\dagger_q - V_{l,q}\b_q \right],
\end{eqnarray}
where $U,V$ are real $L\times L$ matrices; $L$ is the number of lattice sites which equals the number of momentum modes.

Next one searches for $U$ and $V$ such that the following equation,
\begin{equation}
    \comm{\H_a (K)}{\hat{d}^\dagger_l} = w_l(K) \hat{d}^\dagger_l,
    \label{eqHaKDiagonalization}
\end{equation}
is fulfilled. Up to a constant energy shift, this ensures that the Hamiltonian takes the diagonal form $\H_d(K)= \sum_l w_l(K) \hat{d}^\dagger_l \d_l + {\rm const.}$  with dispersion $w_l(K)$. Eq.~\eqref{eqHaKDiagonalization} leads us to an eigenvalue equation of the form
\begin{equation}
    \label{MultiEig}
    \begin{pmatrix} A & B \\ -B & -A \end{pmatrix}
    \begin{pmatrix} U_l\\ V_l\end{pmatrix}
    = w_l(K) \begin{pmatrix} U_l\\ V_l \end{pmatrix}.
\end{equation}
$U_l$($V_l$) denotes the $l$-th column of the matrix $U$($V$) and $w_l(k)$ then correspond to the eigenenergies of the resulting matrix in the left hand site of Eq. \eqref{MultiEig}. Note that this matrix is not Hermitian and thus may lead to complex eigenvalues $w_l(K)$. If this is the case the Hamiltonian \eqref{H_aK} is not diagonalizable\cite{Xiao2009}.

The matrices $A$ and $B$ have the components
\begin{eqnarray*}
A_{pq} &=& \Omega_{K,q}^{(n\xi)} \delta_{pq} -\dfrac{1}{2L} \mathcal{V}^{(n\xi)}_{K,pq}\\
B_{pq} &=& -\dfrac{1}{L} \mathcal{W}^{(n\xi)}_{K,pq}.\\
\end{eqnarray*}
The eigenvalue equation \eqref{MultiEig} is solved numerically by exact diagonalization. Finally, inserting the transformation \eqref{MultiModeTrafo} and rearranging terms such that we can use the eigenvalue equation \eqref{MultiEig}, brings us to the diagonalized form of the Hamiltonian\cite{Xiao2009} \eqref{H_aK},
\begin{multline}
    \H_a(K) = E_{h,\rm eff}^{(n\xi)} + J_{\perp, \text{eff}}^{(n\xi)}\cos(2K) \\+ \sum_l w_l(K) \hat{d}^\dagger_l\d_l + \dfrac{1}{2} \l \sum_l w_l(K) - \text{tr}(A)\r.
\end{multline}
The second term in the second line corresponds to the zero-point energy of the $\hat{d}^\dagger$ bosons.

We can further simplify the Hamiltonian and bring it to the form
\begin{multline}
    \label{PolHLLP}
    \H_a(K) =  E_{h}^{(n\xi)} + \Gamma^{(n\xi)} + J_\perp^{(n\xi)} \cos(2K) \\ + \Delta^{(n\xi)}(K) + \sum_l w_l(K) \hat{d}^\dagger_l\d_l.
\end{multline}
$\Gamma^{(n\xi)}$ is the sum of all contributions coming from the vacuum fluctuations, not depending on the total momentum $K$. In addition we introduced $\Delta^{(n\xi)}(K)$ which sums up all the remaining $K$-dependent terms.
The first four terms of the polaron Hamiltonian~\eqref{PolHLLP} describe the polaron ground state band, i.e. the meson dressed by virtual magnon excitations, see the blue curves in Figs. \ref{Fig2} and \ref{4PolBands}. Each higher band (red) in the mentioned figures corresponds to the creation of one multi-mode boson $\hat{d}^\dagger_l$ with a specific $l$ value. In order to get a better understanding of the resulting polaron spectrum we included curves for the bare meson, $E^{(n\xi)}_h + J_\perp^{(n\xi)}\cos(2K)$ (black solid curve) and the lower and upper edges (black dashed curves) of the non-interacting meson-magnon continuum, defined by $E_h^{n\xi}+J_\perp^{(n\xi)}\cos(2P) + \omega_{K-P}$ with $P=\pi/2$ and $P=0$, respectively.

From inspection of Fig.~\ref{4PolBands} we see that the bandwidth of the meson-magnon continuum (red region) is reduced compared the non-interacting case, due to the interaction with the mesonic impurity, for the shown magnetic field values $h$. We note that for large $h\gg J_z$ the interacting continuum agrees with the non-interacting one. We also notice that the meson bandwidth gets smaller when lowering the magnetic field $h$, while for large $h\gg J_z$ it is equal to the one of the bare meson. Thus for growing magnetic fields $h$ the renormalization of the meson band and the meson-magnon continuum becomes less pronounced. The reason is that the energy gap for creating magnon excitations is growing with $h$. From the polaron spectrum in Fig.~\ref{4PolBands} it is clear that the multi-mode boson operators $\hat{d}^\dagger_l$ at some $l$ resemble a magnon excitation $\hat{b}^\dagger_{K-q}$ with the meson carrying momentum $K$.

\begin{figure}
    \centering

    \includegraphics[width=0.47\textwidth]{FIG10.pdf}
    \caption{Magnetic polaron bands plotted for various magnetic field values with the meson in its ro-vibrational ground state, $n=1,\xi=+1$. The chosen field values are a) $h=0.2J$, b) $h=0.3J$, c) $h=0.4J$ and d) $h=0.5J$. The black solid curve shows the bare meson band without magnon influence and the black dashed lines show the lower and upper edges of the non-interacting meson-magnon scattering continuum; blue lines correspond to the ground state magnetic polaron band; the dark red line corresponds to a meson-magnon bound state. We used system sizes of $L=51$ in our exact diagonalization of the effective polaron Hamiltonian.}
    \label{4PolBands}
\end{figure}

\subsubsection{Meson Magnon Bound State}
\label{SubSecBoundState}
Our theory contains further structure beyond a mere dressing of the meson. We predict a stable bound state of the meson with one local spin flip excitation, i.e. one HP boson. In Fig.~\ref{4PolBands} we see that one magnon band emerges below the meson-magnon scattering continuum but above the polaron ground state band for the chosen values of the staggered magnetic field $h$ (dark red line in Fig.~\ref{4PolBands}). We interpret this isolated band as a meson-magnon bound state as it lies energetically below the scattering continuum. Even for very large values of $h$ (not shown in Fig.~\ref{4PolBands}) we predict the existence of this stable bound state, which is further shown by an analytical calculation in appendix \ref{ApdxBoundStateLargeh}.

Another hint supporting the claim that a HP magnon may bind to the meson was already found in the density plot in Fig.~\ref{Densities}: In the region where the bound state band crosses the polaron ground state band, we see that local spin flips start to accumulate around the meson. This can be treated more formally by making a simplified variational Chevy-type ansatz, as we show explicitly in Sec.~\ref{SubSecChevy1} below.

On several occasions in this work we referred to the bound state as a tetra-parton bound state. The reasoning is that the spin-$1$ magnon excitation can be interpreted as a confined state of two spinons, developing an internal structure itself in the weak-field limit $h \to 0$ \cite{Piazza2015}. Indeed, an interesting future extension for the small-$h$ limit would be to construct an effective spinon-magnon model describing how the magnon resonates in and out of the spinon continuum \cite{Shao2017}, and adding couplings to the chargons as done in the present work. We expect this could improve results for the magnon dispersion, which is poorly represented by our simple linear-spin wave theory when $h/J \to 0$: we checked numerically by calculating the zero-doping dynamical spin structure factor that while the magnon gap is captured well by linear spin-wave theory, the overall shape of the magnon dispersion resembles more closely the spinon dispersion.

\subsection{Chevy approach}
\label{SubSecAvdMgnDcy}
In this section we will discuss an alternative solution of the effective polaron theory, based on Chevy's variational polaron wavefunction \cite{Chevy}. We distinguish two scenarios, depending on the number of magnon excitations we allow in the expansion. For simplicity, we omit the band indices $n,\xi$ of the meson from now.

\subsubsection{Two-magnon state and avoided magnon decay}
\label{SubSecChevy2}
Inspired by the so-called Chevy ansatz, originally introduced for an imbalanced Fermi gas \cite{Chevy}, we expand the ground state of the polaron Hamiltonian \eqref{PolHam} up to two magnon excitations above the bare meson. Taking into account the conservation of the total system momentum, we make the following variational ansatz,
\begin{equation}
    \ket{\psi_k}=  \sqrt{Z_k} \hat{f}^\dagger_{k} \ket{0} + \sum_{pq} \alpha_{k,pq} \hat{f}^\dagger_{k-p-q} \hat{b}^\dagger_p\hat{b}^\dagger_q\ket{0},
    \label{eq2MagChevy}
\end{equation}
where $\sqrt{Z}_k$ and $\alpha_{k,pq}$ are variational parameters satisfying the normalization condition $Z_k + 2 \sum_{pq} |\alpha_{k,pq}|^2=1$.

Minimizing the functional $\mathcal{L} := \bra{\psi} \H_{\text{pol}} - E \ket{\psi}$ with respect to the variational parameters $\sqrt{Z_k}$ and $\alpha_{k,pq}$, leads to the following coupled equations:
\begin{eqnarray}
    \label{sqrtZ}
    \l \varepsilon_k - E\r \sqrt{Z_k} = \dfrac{1}{L} \sum_{pq} \mathcal{W}_{k,pq}^* \alpha_{k,pq},\\
    \label{alpha}
    \sum_{q^\prime}\l \uuline{G}^{-1}_{kp,E}\r_{qq^\prime}\alpha_{k,pq^\prime} = -\dfrac{\sqrt{Z_k}}{2L}\mathcal{W}_{k,pq} ,
\end{eqnarray}
where the $L \times L$ matrix $\uuline{G}^{-1}_{kp,E}$ has the following components:
\begin{multline}
\label{GreensInv}
\l \uuline{G}^{-1}_{kp,E}\r_{qq^\prime} = \l E - \varepsilon_{k-p-q^\prime,\text{eff}} - \omega_p -\omega_{q^\prime} \r \delta_{qq^\prime}\\+ \dfrac{1}{L} \mathcal{V}_{k-p-q^\prime,qq^\prime}.
\end{multline}
The matrix $\uuline{G}_{kp,E}$ can be understood as the free retarted Green's function for two excited magnon excitations propagating along with the meson.

Solving the two above relations for the variational parameters, $\sqrt{Z_k}$ and $\alpha_{k,pq}$, by using a matrix inversion in Eq.~\eqref{alpha}, we get a variational energy of the form
\begin{equation}
    \label{VarEn}
    E_k = \varepsilon_{k,\text{eff}} + \Sigma_k(E_k)
\end{equation}
where we defined the self-energy,
\begin{equation}
    \Sigma_k(E_k) =  \dfrac{1}{2L^2} \sum_{pqq^\prime}\mathcal{W}^*_{k,pq} \l \uuline{G}_{kp,E_k}\r_{qq^\prime}\mathcal{W}_{k,pq^\prime}.
\end{equation}
It is important to note that $\l\uuline{G}_{kp,E_k}\r_{qq^\prime}$ are the components of the inverse of the matrix formed by the elements in Eq.~\eqref{GreensInv}.
Eq.~\eqref{VarEn} has to be solved self-consistently to get the variational energy $E_k$. The result is non-perturbative: it corresponds to re-summation over all diagrams describing two-magnon excitations \cite{Grusdt2015Varenna}. We will present our numerical results below.

Now we use the Chevy approach to argue that the meson at $k=0$ is stable for small magnetic fields $h$ and magnon decay as described in Sec.~\ref{secAvdMagDcySmry} is prevented by meson-magnon interactions. We already elaborated in previous chapters that the meson-magnon scattering continuum and the dressed meson band are repelling each other in this parameter regime. Now we will follow closely the arguments by Verresen et al.~\cite{Verresen2018}: they determined a threshold value where quasi-particle decay of a particle coupled into a continuum of states is prevented by strong interactions.

The condition for the existence of a stable state below the continuum is given by~\cite{Verresen2018}
\begin{equation}
    \dfrac{\varepsilon_{k,\text{eff}} - \omega^{-}_{k}}{|\Sigma_k(\omega^-_{k}+0^-)|} < 1   .
\label{eqCondStbleQP}
\end{equation}
Here $\omega^{-}_{k}$ denotes the lower edge of the meson plus two-magnon scattering continuum which is defined by $\omega^{-}_{k} = \min_{pq} (\varepsilon_{k-p-q,\text{eff}} + \omega_p +\omega_q )$. When the condition \eqref{eqCondStbleQP} is fulfilled, we expect an avoided band crossing between the two-magnon continuum and the meson band.

There are two cases to distinguish in the analysis of the condition~\eqref{eqCondStbleQP}:
\begin{itemize}
    \item When $\varepsilon_{k,\text{eff}} < \omega^{-}_{k}$ the condition~\eqref{eqCondStbleQP} is trivially fulfilled for any interaction strength. In this case the bare meson dispersion $\varepsilon_{k,\text{eff}}$ does not cross the two-magnon scattering continuum.
    \item For $\varepsilon_{k,\text{eff}} > \omega^{-}_{k}$, the integrated meson-magnon interactions must be strong enough, such that the self-energy in the denominator, $|\Sigma_k(\omega^-_{k}+0^-)|$, dominates. At a given interaction strength this can also be achieved by a sufficiently high density of states at low energies. This case corresponds to a non-trivial avoided quasi-particle decay.
\end{itemize}

In the following we analyze the self-energy obtained from the two-magnon Chevy ansatz introduced above. It turns out that the condition~\eqref{eqCondStbleQP} is trivially fulfilled for all values of the staggered magnetic field $h>0$. The reason is that the renormalization of the meson tunneling amplitude, $J_{\perp}^{(n\xi)} \to J_{\perp,\text{eff}}^{(n\xi)}$ by magnon zero-point contributions is strong enough to prevent the meson band from reaching the meson-magnon continuum. We emphasize that this effect results from meson-magnon interactions, which lead to the described renormalization, see Sec.~\ref{SecPolHam}.

In Fig.~\ref{figFCeff} we compare the bare meson tunneling amplitude $J_{\perp}^{(n\xi)}$, neglecting magnon fluctuations, with $J_{\perp,{\rm eff}}^{(n\xi)}$ renormalized by magnon zero-point contributions, assuming the  meson is in its ro-vibrational ground state $n=1,\xi=1$. For small values of the staggered magnetic field $h\lesssim 0.1$, the renormalized tunneling amplitude drops drastically. This is explained by the fact that the zero-point contributions of the magnons start to diverge for $h\to 0$ as we get a large magnon accumulation on all lattice sites of the chain at such low magnetic field values. In this regime, the linear spin-wave theory loses its applicability.

Finally, it is interesting to ask whether the condition~\eqref{eqCondStbleQP} would be fulfilled even if the bare meson band manages to leak into the meson-magnon scattering continuum. To this end, we repeat our above analysis and ignore the renormalization of the meson tunneling amplitude; i.e. we work with $\varepsilon_{k}$ instead of $\varepsilon_{k,\text{eff}}$ in Eq.~\eqref{eqCondStbleQP}. Moreover, we re-scale the magnon dispersion obtained from linear spin-wave theory by a numerical factor $\lambda$, $\omega_q \to \lambda \omega_q$, such that the analytically known spinon bandwidth $J \pi/2$ is correctly captured in the limit $h\to0$.

As shown in Fig.~\ref{figCond} we find that the resulting meson is non-trivially stable at $k=0$ due to meson-magnon interactions. There we plot the left hand side of Eq.~\eqref{eqCondStbleQP} under our simplifying assumptions, for different values of the staggered magnetic field $h$. For $h \lesssim 0.02 J$ non-trivial stabilization is found.

\begin{figure}
    \centering
    \includegraphics[width=0.5 \textwidth]{FIG11.pdf}
    \caption{Effective meson tunneling amplitude $J_{\perp, \text{eff}}^{(n\xi)}$ for $n=1,\xi=1$, which includes a reduction of the spinon tunneling due to magnon vacuum fluctuations compared to the meson tunneling amplitude $J_{\perp}^{(n\xi)}$ for $n=1,\xi=1$. The magnon zero-point contributions stabilize the meson band and prevent a decay into the two-magnon scattering continuum.}
    \label{figFCeff}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=0.5 \textwidth]{FIG12.pdf}
    \caption{Verification of the stability of the meson. We show the left hand side (l.h.s.) of Eq.~\eqref{eqCondStbleQP}, ignoring the renormalization of magnon tunneling due to vacuum magnon fluctuations as described in the text. For values smaller than zero the condition is trivially fulfilled and if the l.h.s. of the condition is between zero and one it is non-trivially fulfilled. We assume the meson to be in its ro-vibrational ground state $n=1,\xi=1$ and the chosen parameters are $t=5, J = 1$.}
    \label{figCond}
\end{figure}

\begin{figure}
    \centering

    \includegraphics[width=0.45\textwidth]{FIG13.png}
    \caption{Comparison of the results of all of our methods at low energies for $h=0.6J$ and $t=5J$: The upper row shows the semi-analytical results, the polaron bands from our LLP treatment with a Gaussian linearization (a) and the curves for two-magnon and one-magnon Chevy ansatz (b). The middle panel c) is the standard one-hole ARPES spectrum, with our theoretical curves plotted on top using the same color scheme as in a) and b). The lower figure d) shows the spin-flip ARPES spectrum with total spin-$3/2$, see \ref{ARPESBound}. In c) and d), for the Gaussian LLP approach we only included the curves corresponding to the polaron ground state and the meson-magnon bound state. The shapes of the theoretical curves and the ARPES dispersion are in good agreement while there is a small shift in total energy. This has already been seen in the ground state energy, see Fig. \ref{GS}.}
    \label{FigComparison}
\end{figure}

\subsubsection{One-magnon state and meson-magnon bound state}
\label{SubSecChevy1}
Instead of the two-magnon ansatz Eq.~\eqref{eq2MagChevy}, mixing even-magnon number states, we can also make an odd-magnon number Chevy ansatz. The lowest order, all one-magnon states should be considered. Since the next order is rather involved, including three magnons, we restrict our discussion to the simplest one-magnon states in the following.

In order to get a variational energy for one-magnon states, we project the full polaron Hamiltonian Eq.~\eqref{PolHam} onto the basis
\begin{equation}
    \ket{q}_k = \hat{f}^\dagger_{k-q}\hat{b}^\dagger_q \ket{0}
\end{equation}
which includes exactly one magnon excitation with variable momentum $q$. The state is constructed such that the total conserved momentum is $k$. We then diagonalize the projected Hamiltonian, defined by its matrix elements
\begin{eqnarray}
    \label{OneMagChevyHam}
    \nonumber \hat{H}_{pq} &=& {}_k\bra{p} \H_{\rm pol} \ket{q}_k\\
    &=& \l \varepsilon_{k-q,\text{eff}} + \omega_q \r \delta_{pq} - \dfrac{1}{2L} \mathcal{V}_{k-q,pq}.
\end{eqnarray}
Since our model Hamiltonian Eq.~\eqref{PolHam} only couples states of equal magnon-number parity, the one- and two-magnon Chevy wavefunctions must be considered independently and cannot couple.

In Fig.~\ref{FigComparison} b) we show the variational one-magnon eigenenergies (dark red curves) obtained by diagonalizing the Hamiltonian matrix in Eq.~\eqref{OneMagChevyHam}. In the same figure we include the result for the variational ground state energy of the two-magnon Chevy ansatz, \eqref{eq2MagChevy} obtained by solving Eq.~\eqref{VarEn} self-consistently (black curve). The obtained spectrum compares qualitatively well to those from the LLP + Gaussian approach, see Fig.~\ref{FigComparison} a). In particular the Chevy approach correctly predicts the broad meson-magnon continuum at high energies, as well as a stable tetra-parton bound state between the meson and the meson-magnon continuum.

\subsection{Comparison of Methods}
In Fig. \ref{FigComparison} we compare the numerical ARPES spectra with our two theoretical methods. We already discussed the polaron spectrum obtained by a Gaussian LLP approach in Sec. \ref{SecPolSpec} which is again shown here in panel a). The variational energies obtained from the two-magnon and one-magnon Chevy ansatz are shown in panel b). We put these semi-analytical curves on top of the numerically obtained standard one-hole ARPES spectrum [panel c)] where the ground state can be compared, and the spin-flip ARPES spectrum probing the sector with $S^z_{\rm tot}=3/2$ [panel d)] where the meson-magnon bound state can be compared.

In our discussion of the ground state energy we observed that there is a small deviation between the numerical and theoretical result, see Fig. \ref{GS}. This small deviation is also seen here in the comparison in Fig. \ref{FigComparison} c). The predicted shapes of the meson dispersion agree very well with each other, and with the full numerical result. We note that this agreement holds for all magnetic field values $h$ that we considered.

For the spin-$3/2$ spin-flip ARPES spectrum the situation is less clear. On the one hand, we find good qualitative agreement in that our numerical td-MPS simulations predict a pronounced quasiparticle peak at low energies, for all considered values of $h$. For the value of $h=0.6 J$ shown in Fig.~\ref{FigComparison} we find remarkable quantitative agreement with our 1-magnon Chevy prediction, which is also relatively close to the LLP prediction. However, for other values of the staggered field $h/J$, we observed larger quantitative differences.

\section{Summary and Outlook}
We have investigated the problem of a single hole doped into a spin chain with an external staggered Zeeman field and studied the interplay of the doped hole with quantum fluctuations of the spin background. In this simplistic setting, we found a remarkably rich zoo of quasiparticle excitations arising from the interplay of charge fluctuations and quantum magnetism. Our predictions can be tested in state of the art quantum gas microscopy experiments.

Based on a parton construction for the hole, capturing the Ising-limit of our model, we developed a simple semi-analytically solvable theory describing the hole as having an internal structure composed of two constituents: a spin-less chargon and a charge-neutral spinon. The hole is a confined mesonic bound state of these constituents, similar to mesons in high-energy physics \cite{Greensite2003}, and their binding potential is linear in nature. Formally, this setting is directly related to a $\mathbb{Z}_2$ lattice gauge theory \cite{Grusdt2020PRL,Kebric2021,Borla2020}. Similar to previous proposals for an analogous problem in 2D~\cite{Bohrdt2020PRB,Bohrdt2021} we find stable long-lived ro-vibrational excitations of the doped hole. Their spectroscopic measurement would constitute compelling evidence that mobile holes have a rich internal structure and can be understood as spinon-chargon bound states.

The main theoretical advancement of our work was the systematic inclusion of quantum fluctuations in the parton description. We used a generalized $1/S$ expansion technique\cite{Grusdt2018} to treat quantum spin fluctuations -- resulting from transverse spin couplings $\propto J_\perp$ in the $t-J$ model -- around the N\'eel state distorted by dominant charge fluctuations. This allowed us to derive an effective polaronic model, describing how a pre-formed meson interacts weakly with additional magnon excitations. As a main result of this approach, we were able to predict an additional stable meson-magnon bound state, i.e. a tetra-parton state of a chargon bound to three spinons. We confirmed this prediction by numerical DMRG simulations and an exact perturbative analysis in the strong-field limit, $h \gg t,J$. Another significant insight obtained by our method concerns the stability of the mesonic quasiparticle when its energy approaches the meson-magnon continuum around zero momentum. Namely, we found evidence that meson-magnon interactions are sufficiently strong to stabilize the quasiparticle peak for arbitrarily weak confining fields $h$.

Our theoretical formalism paves the way for many future extensions. For example, the simplification of the problem to a well-known and weakly coupled Bose-polaron model allows to study far-from equilibrium dynamics \cite{Zvonarev2007,Bruderer2008,Grusdt2014BO,Lampo2017bosepolaronas} of a mobile hole \cite{Ji2021}, going to non-zero temperatures and much longer times than accessible by the more accurate tensor-network methods \cite{Bohrdt2020Dyn,Hubig2020}. Moreover, our approach can be generalized to higher dimensions, where exact numerical methods become significantly more challenging. The dressing of strongly paired states of holes \cite{bohrdt2021strong} can also be investigated. Another promising direction would be the study of magnon excitations in mixed-dimensional systems \cite{Grusdt2018SciPost,Grusdt2020PRL} at finite doping. Finally, the microscopic connection we establish to an underlying $\mathbb{Z}_2$ lattice gauge theory may be more general, suggesting a new route how emergent gauge structures can arise in strongly correlated quantum matter.

The model Hamiltonian we considered in one dimension, namely a $t-J$ model in a staggered Zeeman field, also constitutes an interesting platform for future studies. Its close connections to other interesting models on one hand, such as the 1D $t-J_z$ or $t-J$ models or the 2D $t-J$ model which also has long-range magnetic correlations at zero doping, and its direct experimental realizability in ultracold atoms on the other hand make it an appealing system to study. In this article we limited our discussion to a single doped hole, but extensions to finite doping are straightforward. For example, it will be interesting to search for pairing or charge order at finite doping and investigate the role played by couplings to magnon excitations in the spin background. Exploring the connection to an underlying lattice gauge structure, and a possible breakdown of the meson picture with doping will also be worthwhile endeavours.

\section{Acknowledgements}
The authors benefited particularly from discussions with E. Demler and L. Pollet, and gratefully acknowledge discussions with U. Schollwöck, I. Bloch, F. Pollmann, M. Knap, R. Verresen, T. Shi, P. Bermees, Z. Zhu, P. Wrzosek, K. Wohlfeld and G. Uhrig. This research was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy -- EXC-2111 -- 390814868, from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programm (Grant Agreement no 948141) — ERC Starting Grant SimUcQuam, and by the NSF through a grant for the Institute for Theoretical Atomic, Molecular, and Optical Physics at Harvard University and the Smithsonian Astrophysical Observatory.