Parton theory of magnetic polarons: Mesonic resonances and signatures in dynamics

When a mobile hole is moving in an anti-ferromagnet it distorts the surrounding Neel order and forms a magnetic polaron. Such interplay between hole motion and anti-ferromagnetism is believed to be at the heart of high-Tc superconductivity in cuprates. We study a single hole described by the t-Jz model with Ising interactions between the spins in 2D. This situation can be experimentally realized in quantum gas microscopes. When the hole hopping is much larger than couplings between the spins, we find strong evidence that magnetic polarons can be understood as bound states of two partons, a spinon and a holon carrying spin and charge quantum numbers respectively. We introduce a microscopic parton description which is benchmarked by comparison with results from advanced numerical simulations. Using this parton theory, we predict a series of excited states that are invisible in the spectral function and correspond to rotational excitations of the spinon-holon pair. This is reminiscent of mesonic resonances observed in high-energy physics, which can be understood as rotating quark antiquark pairs. We also apply the strong coupling parton theory to study far-from equilibrium dynamics of magnetic polarons observable in current experiments with ultracold atoms. Our work supports earlier ideas that partons in a confining phase of matter represent a useful paradigm in condensed-matter physics and in the context of high-Tc superconductivity. While direct observations of spinons and holons in real space are impossible in traditional solid-state experiments, quantum gas microscopes provide a new experimental toolbox. We show that, using this platform, direct observations of partons in and out-of equilibrium are possible. Extensions of our approach to the t-J model are also discussed. Our predictions in this case are relevant to current experiments with quantum gas microscopes for ultracold atoms.


I. INTRODUCTION
Understanding the dynamics of charge carriers in strongly correlated materials constitutes an important prerequisite for formulating an effective theory of hightemperature superconductivity. It is generally assumed that the Fermi Hubbard model provides an accurate microscopic starting point for a theoretical description of cuprates [1][2][3]. At strong couplings, this model can be mapped to the t − J model, which describes the motion of holes (hopping t) inside a strongly correlated bath of spins with strong anti-ferromagnetic (AFM) Heisenberg couplings (strength J). To grasp the essence of the complicated t − J Hamiltonian, theorists have also studied closely related variants, most prominently the t − J z model for which Heisenberg couplings are replaced with Ising interactions (J z ) between the spins [2,4].
While the microscopic t − J and t − J z models are easy to formulate, understanding the properties of their ground states is extremely challenging. As a result, most theoretical studies so far have relied on large-scale numerical calculations [5] or effective field theories which often cannot capture microscopic details [6][7][8][9][10][11][12][13]. Even the problem of a single hole propagating in a state with Néel order [4,9,10,[14][15][16][17][18][19][20][21][22][23][24][25][26][27][28], see Fig. 1 (a), is so difficult in general, that heavy numerical methods are required for its solution. This is true in particular at strong couplings t J, J z , where the tunneling rate t of the hole exceeds the couplings J, J z between the spins. The strong coupling regime is also relevant for high-temperature cuprate superconductors for which typically t/J ≈ 3 [3]. While several theoretical approaches have been developed, which are reliable in the weak-tointermediate coupling regime t J, J z , to date there exist only a few theories describing the strong coupling limit [19] and simple variational wavefunctions in this regime are rare. Even calculations of qualitative ground state properties of a hole in an anti-ferromagnet, such as the renormalized dispersion relation, require advanced theoretical techniques. These include effective model Hamiltonians [17,19], fully self-consistent Green's function methods [25], non-trivial variational wave functions [20,23,26] or sophisticated numerical methods such as Monte-Carlo [26,28] and DMRG [27,29] calculations. The difficulties in understanding the single-hole problem add to the challenges faced by theorists trying to unravel the mechanisms of high-T c superconductivity.
Here we study the problem of a single hole moving in an anti-ferromagnet from a different perspective, focusing on the t − J z model for simplicity. In contrast to most earlier works, we consider the strong coupling regime, t J z . Starting from first principles, we derive a microscopic parton theory of magnetic polarons. This approach not only provides new conceptual insights to the physics of magnetic polarons, but it also enables semianalytical derivations of their properties. We benchmark . This distortion carries a fractional spin S z = 1/2. We introduce two types of partons, spinons and holons, describing the spin and charge quantum numbers respectively. They are confined by a string of displaced spins connecting them (blue), similar to mesons which can be understood as bound states of confined quark antiquark pairs. (b) In the strong coupling limit the slow dynamics of the spinon can be decoupled from the fast dynamics of the holon. This leads to coherent motion of the spinon on the original square lattice, which is independent of the holon motion. The latter takes place on a fractal Bethe lattice which is co-moving with the spinon. The holon states on the Bethe lattice obey a discreteĈ4 symmetry that corresponds to rotations of the string configuration around the spinon position. (c) As a result, discrete rotational and vibrational excited states of the magnetic polaron can be formed. They are characterized by quantum numbers nm4..., where e iπm 4 /2 denotes the eigenvalue ofĈ4 and n labels vibrational excitations. In the figure we used labels S, P, D, F for m4 = 0, 1, 2, 3 and indicated the multiplicity of the degenerate states by numbers in circles. Different colors correspond to the different types of ro-vibrational excitations. For the calculation we used S = 1/2 and linear string theory (LST, described later in the text) with strings of length up to max = 100.
our calculations by comparison to the most advanced numerical simulations known in literature. Notably, our approach is not limited to low energies but provides an approximate description of the entire energy spectrum. This allows us, for example, to calculate magnetic polaron dynamics far from equilibrium. Note that in the extreme limit when J z = 0, Nagaoka has shown that the ground state of this model has long-range ferromagnetic order [14]. We will work in a regime where this effect does not yet play a role, see Ref. [27] for a discussion.

A. Partons and the t − Jz model
Partons have been introduced in high energy physics to describe hadrons [30]. Arguably, the most well known example of partons is provided by quarks. In quantum chromodynamics (QCD), the quark model elegantly explains mesons (baryons) as composite objects consisting of two (three) valence quarks. On the other hand, individual quarks have never been observed in nature, and this has been attributed to the strong confining force between a pair of quarks mediated by gauge fields [31]. Even though there is little doubt that quarks are truly confined and can never be separated at large distances, a strict mathematical proof is still lacking, and the quark confinement problem is still attracting considerable attention in high-energy physics, see for example Ref. [32].
To understand how the physics of holes moving in a spin environment with strong AFM correlations is connected to the quark confinement problem, consider removing a spin from a two-dimensional Néel state. When the hole moves around, it distorts the order of the surrounding spins. In the strong coupling regime, t J, J z , these spins have little time to react and the hole can distort a large number of AFM bonds. Assuming for the moment that the hole motion is restricted to a straight line, as illustrated in Fig. 1 (a), we notice that a string of displaced spins is formed. At one end, we find a domain wall of two aligned spins, and the hole is located on the opposite end. By analyzing their quantum numbers, we note that the domain wall corresponds to a spinon -it carries half a spin and no charge -whereas the hole becomes a holon -it carries charge but no spin. The spinon and holon are the partons of our model. Because a longer string costs proportionally more energy, the spinon can never be separated from the holon. This is reminiscent of quark confinement.
Partons also play a role in various phenomena of condensed matter physics. A prominent example is the fractional quantum Hall effect [33,34], where electrons form a strongly correlated liquid with elementary excitations (the partons) which carry a quantized fraction of the electron's charge [35][36][37]. This situation is very different from the case of magnetic polarons which we consider here, because the fractional quasiparticles of the quantum Hall effect are to a good approximation non-interacting and can be easily separated. Similar fractionalization has also been observed in one-dimensional spin chains [38][39][40][41][42][43][44], where holes decay into pairs of independent holons and spinons as a direct manifestation of spin-charge separation. Unlike in the situation described by Fig. 1 (a), forming a string costs no energy in one dimension and spinons and holons are deconfined in this case.
Confined phases of partons are less common in con-densed matter physics. It has first been pointed out by Béran et al. [45] that this is indeed a plausible scenario in the context of high-temperature superconductivity, and the t − J model in particular. In Ref. [45] theoretical calculations of the dispersion relation and the optical conductivity of magnetic polarons were analyzed, and it was concluded that their observations can be well explained by a parton theory of confined spinons and holons. A microscopic description of those partons has not yet been provided, although several models with confined spinonholon pairs have been studied [11,12,46]. The most prominent feature of partons that has previously been discussed in the context of magnetic polarons is the existence of a set of resonances in the single-hole spectral function [21,25,28,45,[47][48][49] which can be measured by angle-resolved photoemission spectroscopy (ARPES), see e.g. Ref. [50] for a discussion. Such longlived states in the spectrum can be understood as vibrational excitations of the string created by the motion of a hole in a Néel state [15,16,18,47,49]. In the parton theory they correspond to vibrational excitations of the spinon-holon pair [45], where in a semi-classical picture the string length is oscillating in time.
In this paper we present additional evidence for the existence of confined partons in the two-dimensional t − J z model at strong coupling. Using the microscopic parton theory, we show that besides the known vibrational states an even larger number of rotational excitations of magnetic polarons exist. This leads to a complete analogy with mesons in high-energy physics which we discuss next (in I B). The rotational excitations of magnetic polarons have not been discussed before, partly because they are invisible in traditional ARPES spectra. Quantum gas microscopy [51,52] represents a new paradigm for studying the t − J z model, and we discuss below (in I C) how it enables not only measurements of rotational excitations, but also direct observations of the constituent partons in current experiments with ultracold atoms.

B. Rotational excitations of parton pairs
Mesons can be understood as bound states of two quarks and thus are most closely related to the magnetic polarons studied in this paper. The success of the quark model in QCD goes far beyond an explanation of the simplest mesons, including for example pions (π) and kaons (K). Collider experiments that have been carried out over many decades have identified an ever growing zoo of particles. Within the quark model, many of the observed heavier mesons can be understood as excited states of the fundamental mesons. Aside from the total spin s, heavier mesons can be characterized by the orbital angular momentum of the quark antiquark pair [53] as well as the principle quantum number n describing their vibrational excitations. In Table I   Examples of meson resonances corresponding to rotational ( ) and vibrational (n) excitations of the quark antiquark pair (qq). This list is incomplete and the data was taken from Ref. [54]. The numbers in brackets denote the mass of the excited meson state in units of M eV /c 2 .
the fundamental pion (kaon) state π (K), many rotational states with = 1, 2, 3 (P, D, F ) can be constructed [53,54] which have been observed experimentally [55]. By changing n to two, the excited states π(1300) and K(1460) can be constructed. Due to the deep theoretical understanding of quarks, all these mesons are considered as composites instead of new fundamental particles. Similarly, rotationally and vibrationally excited states of magnetic polarons can be constructed in the t − J z model. They can be classified by the angular momentum (rotational) and radial (vibrational) quantum numbers and n of the spinon-holon pair, as well as the spin σ of the spinon. An important difference to mesons is that we consider a lattice model where the usual angular momentum is not conserved. However, there still exist discrete rotational symmetries which can be defined in the spinon reference frame.
To understand this, let us consider the Hilbert space H p of the parton theory introduced in this paper. As illustrated in Fig. 1 (b), it can be described as a direct product of the space of spinon positions on the square lattice H s and the space of string configurations H Σ emerging from the spinon, The latter is equivalent to the Hilbert space of a single particle hopping on a fractal Bethe lattice with four bonds emerging from each site. We recognize a discrete four-fold rotational symmetryĈ 4 of the parton theory, where the string configurations are cyclically permuted around the spinon position, see Fig. 1 (b). Therefore we can construct excited magnetic polaron states with eigenvalues e iπm4/2 ofĈ 4 with m 4 = 1, 2, 3 (P, D, F ), which are analogous to the rotational excitations of mesons. In the ground state, m 4 = 0. Similarly, there exist threefold permutation symmetriesP 3 in the parton theory corresponding to cyclic permutations of the string configuration around sites one lattice constant away from the spinon. This leads to a second quantum number m 3 = 0, 1, 2 required to classify all eigenstates. In Fig. 1 (c) we calculate the excitation energies of rotational and vibrational states in the magnetic polaron spectrum. We applied the linear string approximation, where self-interactions of the string connecting spinon and holon are neglected. It will be shown in the main part of this paper, that this description is justified for the low lying excited states of magnetic polarons in the t − J z model. In analogy with the meson resonances listed in Table I, we have labeled the eigenstates in Fig. 1 (c) by their ro-vibrational quantum numbers. Similar to the pion, the magnetic polaron corresponds to the ground state 1S. The lowest excited states are given by 1P, 1D, 1F . In contrast to their high-energy analogues [32] these states are degenerate which will be shown to be due to lattice effects. Depending on the ratio of J z /t, the next higher states correspond to a vibrational excitation (2S), or a second rotational excitation (states 1m 4 m 3 for m 3 = P, D and m 4 = S, P, D, F ).
The parton theory of magnetic polarons provides an approximate description over a wide range of energies at low doping. This makes it an excellent starting point for studying the transition to the pseudogap phase observed in cuprates at higher hole concentration [3,50,56], because the effective Hilbert space is not truncated to describe a putative low-energy sector of the theory. We discuss extensions of the parton theory to finite doping and beyond the t − J z model in a forthcoming work.

C. Quantum gas microscopy of the t − Jz model
Experimental studies of individual holes are challenging in traditional solid state systems. Ultracold atoms in optical lattices provide a promising alternative platform for realizing this scenario and investigating microscopic properties of individual holes in a state with Néel order. The toolbox of atomic physics offers unprecedented coherent control over individual particles. In addition, many powerful methods have been developed to probe these systems, including the ability to measure correlation functions [57][58][59][60], non-local string order parameters [44] and spin-charge correlations [60] on a single-site level. Moreover bosonic [51,52] and fermionic [61][62][63][64][65][66][67] quantum gas microscopes offer the ability to realize arbitrary shapes of the optical potential down to length scales of a single site as set by the optical wavelength [68,69].
These capabilities have recently led to the first realization of an anti-ferromagnet with Néel order across a finite system of ultracold fermions with SU (2) symmetry at temperatures below the spin-exchange energy scale J [70]. In another experiment, canted anti-ferromagnetic states have been realized at finite magnetization [67] and the closely related attractive Fermi Hubbard model has been investigated [71]. In systems of this type individual holes can be readily realized in a controlled setting and studied experimentally.
The t − J z model was long considered as a mere toy model, closely mimicking some of the essential features known from more accurate model Hamiltonians relevant The t − Jz model can be implemented for ultracold atoms with two internal (pseudo-) spin states in a quantum gas microscope by starting from a Mott insulator (a). Ising interactions can be realized by including Rydberg dressing for one of the two spin states, as demonstrated experimentally in Refs. [73,74]. Mobile holes can be doped into the system by removing atoms from the Mott insulator. Their statistics is determined by the underlying particles forming the Mott state. In principle the method works in arbitrary dimensions, although many Rydberg states have anisotropic interactions which can break the lattice symmetries.
in the study of high-temperature superconductivity. By using ultracold atoms one can go beyond this paradigm, and realize the t − J z model experimentally.
One approach suggested in Ref. [72] is to use ultracold polar molecules which introduce anisotropic and longrange dipole-dipole couplings between their internal spin states. The flexibility to adjust the anisotropy also allows to realize Ising couplings as required for the t−J z model. When the molecules are placed in an optical lattice, the ratio t/J z between tunneling and spin-spin couplings can be tuned over a wide range.
Other methods to implement Ising couplings between spins involve trapped ions [75,76] or arrays of Rydberg atoms [73,74,77]. We now discuss the second option in more detail, because it allows a direct implementation of the t − J z model in a quantum gas microscope with single-particle and single-site resolution.
As illustrated in Fig. 2 (a), one can start from a Mott insulating state of fermions or bosons with two internal (pseudo-) spin states. Strong Ising interactions between the spins can be realized by Rydberg dressing only one of the two states, see Fig. 2 (b). This situation can be described by the effective Hamiltonian [73,78,79] where nearest neighbor AFM Ising couplings J z are dominant. By doping the Mott insulator with holes, this allows to implement an effective t − J z Hamiltonian with tunable coupling strengths.
The statistics of the holes in the resulting t − J z model are determined by the statistics of the underlying particles forming the Mott insulator. For the study of a single magnetic polaron in this paper, quantum statistics play no role. At finite doping magnetic polarons start to interact and their statistics become important. Studying the effects of quantum statistics on the resulting many-body states is an interesting future direction.

Direct signatures of strings and partons
Quantum gas microscopes provide new capabilities for the direct detection of the partons constituting magnetic polarons, as well as the string of displaced spins connecting them. The possibility to perform measurements of the instantaneous quantum mechanical wavefunction directly in real space allows one to detect non-local order parameters [44,81] and is ideally suited to unravel the physics underlying magnetic polarons. Now we provide a brief summary of the most important signatures of the parton theory which can be directly accessed in quantum gas microscopes and will be discussed in this paper. The following considerations apply to a regime of temperatures T < J z where the local anti-ferromagnetic correlations are close to their zero-temperature values, although most of the phenomenology is expected to be qualitatively similar at higher temperatures [82].
Ro-vibrational excitations.-As mentioned in Sec. I A, the ro-vibrational excitations of magnetic polarons provide direct signatures for the parton nature of magnetic polarons. Their energies can be directly measured: Vibrational states are visible in ARPES spectra, which can also be performed in a quantum gas microscope [83]. Later we also discuss alternative spectroscopic methods based on magnetic polaron dynamics in a weakly driven system, which enable direct measurements of the rotational resonances.
Direct imaging of strings and partons.-In a quantum gas microscope, the instantaneous spin configuration around the hole can be directly imaged [60]. Up to loop configurations, this allows one to directly observe spinons, holons and strings and extract the full counting statistics of the string length for example. We show in this paper that this method works extremely accurately in the case of the t − J z model.
Scaling relations at strong couplings.-When t J z , the motion of the holon relative to the spinon can be described by an effective one-dimensional Schrödinger equation with a linear confining potential at low energies. This leads to an emergent scaling symmetry which allows to relate solutions at different ratios J z /t by a simple re-scaling of lengths [15]: x → λ 1/3 x when J z → λJ z . In ultracold atom setups the ratio t/J z can be controlled and a wide parameter range can be simulated. For doping with a single hole this allows to observe the emergent scaling symmetry by showing a data collapse after rescaling all lengths as described above. The scaling symmetry also applies for the expectation values of potential and kinetic energies in the ground state. In the strong coupling regime, t J z , to leading order both depend linearly on t 1/3 J 2/3 z [15]. By simultaneously imaging spin and hole configurations [60], the potential energy can be directly measured using a quantum gas microscope and the linear dependence on (J z /t) 2/3 can be checked.

Far-from-equilibrium experiments
Ultracold atoms allow a study of far-from equilibrium dynamics of magnetic polarons [82,[84][85][86][87][88][89]. For example, a hole can be pinned in a Néel state and suddenly released. This creates a highly excited state with kinetic energy of the order t J z , which is quickly transferred to spin excitations [85,86,90]. Such dynamics can be directly observed in a quantum gas microscope. It has been suggested that this mechanism is responsible for the fast energy transfer observed in pump-probe experiments on cuprates [88], but the coupling to phonons in solids complicates a direct comparison between theory and experiment.
As a second example, external fields (i.e. forces acting on the hole) can be applied and the resulting transport of a hole through the Néel state can be studied [85]. As will be shown, this allows direct measurements of the mesonic excited states of the magnetic polaron, analogous to the case of polarons in a Bose-Einstein condensate [91,92].

D. Magnetic polaron dynamics
To study dynamics of magnetic polarons in this paper, we use the strong coupling parton theory to derive an effective Hamiltonian for the spinon and holon which describes the dynamics of a hole in the AFM environment. By convoluting the probability densities for the holon and the spinon, we obtain the density distribution of the hole, which can be directly measured experimentally. Even though the spinon dynamics are slow compared to the holon motion at strong coupling, they determine the hole distribution at long times because the holon is bound to the spinon. We consider different nonequilibrium situations which can all be realized in current experiments with ultracold atoms.
Benchmark.-We benchmark our parton theory by comparing to time-dependent quantum Monte Carlo calculations [82,89] of the hole dynamics in the twodimensional t − J z model. To this end we study the far-from equilibrium dynamics of a hole which is initialized in the system by removing the central spin from the Néel state. A brief summary of the quantum Monte We benchmark the strong coupling parton theory of magnetic polarons by comparing to numerically exact timedependent quantum Monte Carlo (QMC) calculations. A quantum quench in the t − Jz model with spin S = 1/2 is considered, where initially a single hole is created in a Néel state on a square lattice by removing the spin on the central site. We calculate the root mean square distance, see Eq. (55), of the hole from the origin (a) and the return probability of the hole (b). For all values of Jz/t and for times accessible by the time-dependent Monte Carlo approach we obtain excellent agreement with calculations based on the linear string theory (LST) described in this paper. Non-linear string theory (NLST) predicts only small corrections and agrees with the QMC results slightly better at intermediate times. For details on the QMC method, see Appendix A. We compare results from NLST and LST at longer times in Appendix B.
Carlo method used for solving this problem can be found in Appendix A.
In Fig. 3 we show our results for the return probability n h (0) and the extent of the hole wave function x 2 h , for times accessible with our quantum Monte Carlo method. Heren h (x) is the density operator of the hole at site x, andx h is the position operator of the hole in first quantization. For all considered values of J z we obtain excellent agreement of the strong coupling parton theory with numerically exact Monte Carlo results.
Pre-spin-charge separation.-For the largest considered value of J z = t we observe a pronounced slow-down of the hole expansion in Fig. 3. This is due to the restoring force mediated by the string which connects spinon and holon. However, the expansion does not stop completely [4,18]. Instead it becomes dominated by slow spinon dynamics at longer times, as we show by an explicit calculation in Sec. V D and Fig. 22.
For smaller values of J z /t, the hole expansion slows down at later times before it becomes dominated by spinon dynamics. In the strong coupling regime, t J z , we obtain a large separation of spinon and holon time scales. This can be understood as a pre-cursor of spincharge separation: although the holon is bound to the spinon, it explores its Hilbert space defined by the Bethe lattice independently of the spinon dynamics. As illustrated in Fig. 1 (b) the entire holon Hilbert space is comoving with the spinon. This is a direct indicator for the parton nature of magnetic polarons.
At short-to-intermediate times the separation of spinon and holon energies gives rise to universal holon dynamics. Indeed, the expansion observed in Fig. 3 at strong coupling t J z is similar to the case of hole propagation in a spin environment at infinite temperature [89]. In that case an approximate mapping to the holon motion on the Bethe lattice is possible too [82].
Coherent spinon dynamics.-We consider a situation starting from a spinon-holon pair in its ro-vibrational ground state. In contrast to the far-from equilibrium dynamics discussed above, the holon is initially distributed over the Bethe lattice in this case. We still start from a state where the spinon is localized in the center of the system. In this case there exist no holon dynamics on the Bethe lattice within the strong coupling approximation, and a measurement of the hole distribution allows to directly observe coherent spinon dynamics.
We also present an adiabatic preparation scheme for the initial state described above, where the magnetic polaron is in its ro-vibrational ground state. The scheme can be implemented in experiments with ultracold atoms. The general strategy is to first localize the hole on a given lattice site by a strong pinning potential. By slowly lowering the strength of this potential, the ro-vibrational ground state can be prepared with large fidelity, as we demonstrate in Fig. 19. Details are discussed in Sec. V B.
Spectroscopy of ro-vibrational excitations.-To test the strong coupling parton theory experimentally, we suggest measuring the energies of rotational and vibrational eigenstates of the spinon-holon bound state directly. In Sec. V C we demonstrate that rotational states can be excited by applying a weak force to the system. As before, we start from the ro-vibrational ground state of the magnetic polaron. The force induces oscillations of the density distribution of the hole, which can be directly measured in a quantum gas microscope. We demonstrate that the frequency of such oscillations is given very accurately by the energy of the first excited state, which has a non-trivial rotational quantum number. This is another indicator for the parton nature of magnetic polarons.
Vibrational excitations can be directly observed in the spectral function [28]. In Sec. V E we briefly explain its properties in the strong coupling regime. Possible measurements with ultracold atoms are also discussed.

E. Outline
This paper is organized as follows. In Sec. II we introduce the microscopic parton theory describing holes in an anti-ferromagnet in the strong coupling regime, starting from first principles. We solve the effective Hamiltonian in Sec. III and derive the rotational and vibrational excited states of magnetic polarons. Direct signatures in the string-length distribution and its measurement in a quantum gas microscope are also discussed. Sec. IV is devoted to a discussion of the effective spinon dispersion relation. We derive and benchmark a semi-analytical tight-binding approach to describe the effects of Trugman loops, the fundamental processes underlying spinon dynamics in the t − J z model. In Sec. V we apply the strong coupling parton theory to solve different problems involving magnetic polaron dynamics, which can be realized in experiments with ultracold atoms. Extensions of the parton theory to the t − J model are discussed in Sec. VI. We close with a summary and by giving an outlook in Sec. VII.

II. MICROSCOPIC PARTON THEORY OF MAGNETIC POLARONS
In this section we introduce the strong coupling parton theory of holes in the t − J z model, which builds upon earlier work on the string picture of magnetic polarons [15,16,18,49]. After introducing the model in II A and explaining our formalism in II B we derive the parton construction in Sec. II C.
A. The t − Jz model For a single hole, j,σĉ † j,σĉ j,σ = N − 1, where N is the number of lattice sites, the t − J z Hamiltonian can be written aŝ Hereĉ † j,σ creates a boson or a fermion with spin σ =↑, ↓ on site j andP projects onto the subspace without double occupancies.
i,j denotes a sum over all bonds i, j between neighboring sites i and j, where every bond is counted once. The spin operators are defined byŜ z j = σ,τĉ † j,σ σ z σ,τĉj,τ /2. The second line of Eq. (3) describes the hopping of the hole with amplitude t and the first line corresponds to Ising interactions between the spins. See e.g. Ref. [4] and references therein for previous studies of the t−J z model. Experimental implementations of the t − J z Hamiltonian were discussed in Sec. I C 1. From now on we consider the case whenĉ † j,σ describes a fermion for concreteness, but as long as a single hole is considered the physics is identical if bosons were chosen.

Schwinger-boson representation and constraint
For our discussion of the t − J z model (3) we find it convenient to choose a parameterization in terms of Schwinger-bosonsb jσ and spinless fermionic holon operatorsĥ i satisfying the following constraint, In the original Hamiltonian from Eq. (3), the length of the spins is S = 1/2 and Eq. (4) is equivalent to the condition of no double occupancy of lattice sites. More generally, for S ≥ 1/2, Eq. (4) ensures that a given lattice site is either occupied by a single holon and no Schwingerbosons or no holon but exactly 2S Schwinger-bosons. From now on we will consider more general models with arbitrary integer or half-integer values of S. This approach is similar to the usual 1/S-expansion of the t−J and related Hamiltonians, see e.g. Ref. [93], except that in the latter a different constraint is used: σb † jσb jσ + h † jĥ j = 2S. For S = 1/2 the two constraints are identical, and we discuss in Appendix C how they differ for larger S > 1/2. In both cases, the spin operators arê In terms of holons and Schwinger-bosons, the second term in the t − J z Hamiltonian becomeŝ whereF † ij (S) involves only Schwinger-bosons and will be explained shortly. Note that here, in contrast to Eq. (3), the hopping rate t comes with a positive sign because holon creation corresponds to fermion annihilation,ĉ † jĉ i →ĥ jĥ † i = −ĥ † iĥ j . The termF † ij (S) describes the re-ordering of the spins during the hopping of the holon. This is required by the Schwinger-boson constraint in Eq. (4). In particular, F ij (S) describes how the spin state on site j is moved to site i while the hole is hopping from i to j. In this work, we are interested in the case S = 1/2, whereF ij (S) becomes [93] see Fig. 4 (a). Because the constraint Eq. (4) is fulfilled, the projectorsP from Eq. (3) can be dropped in the Schwinger-boson representation.

B. Generalized 1/S expansion
In this section we introduce a generalization of our model to large spins, S > 1/2. Technically we only achieve a re-formulation of the original Hamiltonian, and in the case of the t − J z model no real progress is made. However, as discussed in more detail in Appendix C, the formalism developed in this section can be straightforwardly generalized to include quantum fluctuations. In such cases, the leading order generalized 1/S expansion represents a t − J z model again, and the parton theory which we develop in Sec. II C can be applied here, too. This establishes our parton construction as a valuable starting point for analyzing a larger class of models. In a forthcoming work we apply the generalized 1/S expansion introduced below to situations including quantum fluctuations in the spin environment described by a general XXZ Hamiltonian [94,95].
In the conventional 1/S expansion, the holon motion is described by the operatorF ij (S) ≡F ij (1/2) from Eq. (8) and S only enters in the corresponding Schwinger-boson constraint, σb † jσb jσ +ĥ † jĥ j = 2S. As we explain in detail in Appendix C 1, this approach cannot capture strong distortions of the local Néel order parameter, or the local staggered magnetization, defined bŷ Here (−1) j denotes the sublattice parity, which is +1 (−1) for j from the A (B) sublattice. Within the conventional extension of the t − J z model to large values of S, the motion of the holon from site i to j is accompanied by changes of the spins S z i and S z j by ±1/2, see Fig. 4  (b). As a result, the sign of the local Néel order parameter Ω j cannot change when S 1/2 is large, unless the holon performs multiple loops.
To avoid these problems of the conventional 1/S expansion, and to ensure that the generalized Schwingerboson constraint Eq. (4) is satisfied byĤ t in Eq. (7), we replaceF ij (1/2) from Eq. (8) by a new termF ij (S). As we explain in detail in Appendix C 2, this generalized holon-hopping operatorF ij (S) describes a transfer of the entire spin state from site j to i, see also Fig. 4 (c).

Formalism: Ising variables and the distortion field
Now we introduce some additional formalism which is useful for the formulation of the microscopic parton theory. Our discussion is kept general and applies to arbitrary values of the spin length S.
Zero doping.-To describe the orientation of the local spin of length S, we introduce an Ising variableτ z j on the sites of the square lattice, For S = 1/2 the Ising variableτ z j = 2Ŝ z j is identical to the local magnetization. For S > 1/2 the situation is different becauseτ z j = ±1 can still only take two possible values, whereasŜ z j = −S, −S + 1, ..., S. The classical Néel state with AFM ordering along the z-direction corresponds to a configuration whereτ z j = +1 on the A-sublattice, andτ z j = −1 on the B-sublattice. To take the different signs into account, we define another Ising variableτ z j describing the staggered magnetization, The Néel state corresponds to the configuration [96] τ z j ≡ 1 for all j.
Doping.-Now we consider a systems with one hole. Our general goal in this section is to construct a complete set of one-hole basis states. This can be done starting from the classical Néel state by first removing a spin of length S on site j and next allowing for distortions of the surrounding spins. In this process, we assign the value ofτ z j = 1 to the lattice site where the hole was created. Note that this value is associated with a sublattice index of the hole and it reflects the spin σ which was initially removed when creating the hole.
Distortions.-The holon motion, described byĤ t , introduces distortions into the classical Néel state. Using the staggered Ising variableτ z j they correspond to sites withτ z j = −1. We will now show that the t − J z Hamiltonian can be expressed entirely in terms of the product defined on links,σ which will be referred to as the distortion field. On bonds withσ z i,j = 1 (respectivelyσ z i,j = −1) the spins are anti-aligned (aligned), see Fig. 1 (a) for an illustration.
Effective Hamiltonian.-The termĤ J in the t − J z Hamiltonian Eq. (3) can be re-formulated aŝ The first term corresponds to the ground state energy of the undistorted Néel state, where N is the number of lattice sites. The second term describes the energy cost of creating distortions. In Appendix C we explain how quantum fluctuations can be included within the generalized 1/S expansion. The termĤ t in the t−J z Hamiltonian Eq. (3) can also be formulated in terms of the distortion field. Consider the motion of the holon from site i to j. This corresponds to a movement of the spin on site j to site i, which changes the distortion fieldσ z l on links l = r, k including sites r = i and r = j. Such changes depend on the original orientation of the involved spins on sites k and can be described by the operatorsσ x l . We obtain the expression Becauseσ x i,j andσ z i,j are not commuting, the distortion fieldσ z i,j begins to fluctuate in the presence of the mobile hole with t = 0. The effect of the hole hopping term (15) is illustrated in Fig. 5.
By combining Eqs. (14) and (15) we obtain an alternative formulation of the single-hole t − J z model at S = 1/2. Even when S > 1/2, the leading order result in the generalized 1/S expansion is a similar effective S = 1/2 Hamiltonian, formulated in terms of the same distortion field. This method goes beyond the conventional 1/S expansion, where the distortion fieldσ z i,j ≡ 1 is kept fixed on all bonds and does not fluctuate. To describe the distortion of the Néel state introduced by the holon motion in the conventional 1/S expansion, one has FIG. 5. When a hole is moving in a spin state with AFM order, it distorts the Néel order (green) and creates a string of displaced spins (blue), (a)-(c). In the formalism used for the generalized 1/S expansion this is described by the coupling of the holon to the distortion field σ z i,j on the bonds. The end of the string (red) can be associated with a parton, the spinon, which carries the spin S of the magnetic polaron. Its charge is carried by another parton, the holon (gray).
to resort to magnon fluctuations on top of the undisturbed Néel state, which only represent sub-leading corrections in the generalized 1/S formalism. This property of the generalized 1/S expansion makes it much more amenable for an analytical description of the strong coupling regime where t J and the Néel state can be substantially distorted even by a single hole.
C. The spinon-holon picture and string theory So far we have formulated the t − J z model using two fields, whose interplay determines the physics of magnetic polarons: The holon operatorĥ j and the distortion field σ i,j on the bonds. By introducing magnonsâ j , more general models with quantum fluctuations can also be considered, but such terms are absent in the t − J z case. Our goal in this section is to replace the distortion field σ i,j by a simpler description of the magnetic polaron, which is achieved by introducing partons.

Spinons and holons
The Hamiltonian (15) describing the motion of the holon in the distorted Néel state determined byσ i,j , is highly non-linear. To gain further insights, we study more closely how the distortion fieldσ z i,j is modified by the holon motion. In particular we will argue that it carries a well-defined spin quantum number.
Quantum numbers.-Let us start from the classical Néel state and create a hole by removing the spin on the central site of the lattice. This changes the total charge Q and spin S z of the system by ∆Q = −1 and ∆S z = ±S, where the sign of ∆S z depends on the sublattice index of the central site. When the hole is moving, both S z and Q are conserved and we conclude that the magnetic polaron (mp) carries spin S z mp = ±S and charge Q mp = −1. There exists no true spin-charge separation for a single hole in the 2D Néel state [15,18,28], i.e. the spin degree of freedom of the magnetic polaron cannot completely separate from the charge. We can understand this for the case S = 1/2, where the magnetic polaron carries fractional spin S z mp = ±1/2. Because the elementary spin-wave excitations of the 2D anti-ferromagnet carry spin S z = ±1, see e.g. Ref. [93], they cannot change the fractional part S z mp mod 1 of the magnetic polaron's spin, which is therefore bound to the charge. This is in contrast to the 1D case, where fractional spinon excitations exist in the spin chain even at zero doping [38] and the hole separates into independent spinon and holon quasiparticles [39][40][41][42][43][44].
Partons.-Now we show that the magnetic polaron in a 2D Néel state can be understood as a bound state of two partons, the holonĥ j carrying charge and a spinon s i carrying spin. In the strong coupling regime, t J z , we predict a mesoscopic precursor of spin-charge separation: While the spinon and holon are always bound to each other, their separation can become rather large compared to the lattice constant, and they can be observed as two separate objects. Their bound state can be described efficiently by starting from two partons with an attractive interaction between them, similar to quarks forming a meson. When J z = 0 spinon-holon pairs can be completely separated [82,89].
In contrast to the usual slave-fermion (or slave-boson) approach [97], we will not define the holon and spinon by breaking up the original fermionsĉ j,σ on site j. Instead we notice that the spin quantum number S z mp of the magnetic polaron is carried by the distortion fieldσ z i,j . The latter determines the distribution of the spin on the square lattice. We have already introduced the spin-less holon operatorsĥ j in Eq. (4) by using a Schwinger-boson representation of the t − J z model.
Spin and charge distribution.-To understand how well the spin S z mp = ±1/2 of the magnetic polaron is localized on the square lattice, we study the motion of the holon described by Eq. (15). When the hole is moving it leaves behind a string of displaced spins [15,18] The stringtheory becomes exact in 1D. This picture is equivalent to the squeezed-space description of spin-charge separation [42], where the holon occupies the bonds between neighboring sites of a spin chain, see also Ref. [83]. −1) emerge, unless the hole returns to the origin. This corresponds to a surplus of spin on this site relative to the original Néel state, and we identify it with the location x s of the spinon. The spin σ of the spinon is opposite for the two different sublattices. Our definition of spinons can be considered as a direct generalization of domain walls in the one-dimensional Ising model, see Fig. 6 (a). In both cases, the fractional spin carried by the spinon is not strictly localized on one lattice site but extends over a small region around the assigned spinon position. To demonstrate this explicitly, we calculate the average magnetization 2 Ŝ z (x − x s ) in the spinon frame in Fig. 7 (a). We use the t − J z model in the strong coupling regime with t = 10J z . We observe that the checkerboard pattern of the Néel order parameter is completely retained, except for a spin-flip in the center. This shows that the spin of the magnetic polaron S z mp is localized around the spinon position, i.e. at the end of the string defined by the holon trajectory.
This result should be contrasted to the magnetization calculated in the holon frame, 2 Ŝ z (x − x h ) where x h is the holon position, shown in Fig. 7 (b). In that case, the spin of the magnetic polaron is distributed over a wide area around the holon. The AFM checkerboard structure is almost completely suppressed because it is favorable for the holon to delocalize equally over both sublattices at this large value of t/J z = 10. Similarly, the charge distribution n h (x − x s ) covers an extended area around the spinon, see Fig. 7 (c).
Formal definition of spinons.-Formally we add an additional labelŝ † i,σ |0 to the quantum states. Here i denotes the site of the spinon as defined above, and the spin index σ depends on the sublattice index of site i and will be suppressed in the following. When holon trajectories are included which are not straight but return to the origin, this label is not always unique, see Ref. [18] or Sec. IV A. Thus, by adding the new spinon label to the wave function, we obtain an over-complete basis. We will deal with this issue later in Sec. IV and argue that the use of the over-complete basis is a useful approach.
We also note that the spinon label basically denotes the site where we initialize the hole and let it move around to construct the basis of the model. In a previous work by Manousakis [49], this site has been referred to as the "birth" site without drawing a connection to the magnetization (the spin) of the magnetic polaron localized around this site, as shown in Fig. 7 (a).

String theory: an over-complete basis
Our goal in this section is to describe the distortion fieldσ i,j by a conceptually simpler string on the square lattice. This idea goes back to the works by Bulaevskii et al. [15] and Brinkman and Rice [16], as well as works by Trugman [18] and more recently by Manousakis [49]. In Refs. [18,49] a set of variational states was introduced, based on the intuition that the holon leaves behind a The string description can be obtained by replacing this basis by a closely related, but conceptually simpler set of basis states. When the holon propagates in the Néel state, starting from the spinon position, it modifies the distortion field σ z i,j differently depending on the trajectory Σ it takes. Here we use the convention that trajectories Σ are defined only up to self-retracing components, in contrast to paths which contain the complete information where the holon went. The holon motion thus creates a memory of its trajectory in the spin environment.
Given a trajectory Σ and the spinon position, we can easily determine the corresponding distortion field σ z i,j (Σ, x s ). In the following we will assume that for all relevant quantum states, the opposite is also true. Namely, that given the distortion of the Néel stateσ z i,j , we can reconstruct the trajectory Σ defined up to selfretracing components, as well as the spinon position. We will show that this is an excellent approximation. Using a quantum gas microscope this one-to-one correspondence can be used for accurate measurements of holon trajectories Σ in the Néel state by imaging instantaneous spin and hole configurations. We analyze the efficiency of this mapping in detail in Sec. III D.
In some cases our assumption is strictly correct, for example in the one-dimensional Ising model. In that case the spinon corresponds to a domain wall in the antiferromagnet, see Fig. 6 (a). When its location is known, as well as the distance of the holon from the spinon (i.e. the trajectory Σ), the spin configurationσ z i,j can be re-constructed, see Fig. 6 (b). A second example, which is experimentally relevant for ultracold atoms, involves a model where the hole can only propagate along one dimension inside a fully two-dimensional spin system [98].
For the fully two-dimensional magnetic polaron problem, there exist sets of different trajectories Σ which give rise to the same spin configurationσ z i,j . Trugman has shown [18] that the leading-order cases correspond to situations where the holon performs two steps less than two complete loops around an enclosing area, see Fig. 12 (a). Choosing a single plaquette, this requires a minimum of six hops of the holon before two states cannot be distinguished by the corresponding holon trajectory anymore.
By performing six steps, a large number of states can be reached in principle: In the first step, starting from the spinon, there are four possible directions which the holon can choose, followed by three possibilities for each of the next five steps. This makes a total of 4 × 3 5 = 972 states over which the holon tends to delocalize in order to minimize its kinetic energy. In contrast, there are only eight distinct Trugman loops involving six steps which lead to spin configurations that cannot be uniquely assigned to a simple holon trajectory.
In the following we will use an over-complete set of basis states, labeled by the spinon position and the holon Σ will be referred to as the string which connects the spinon and the holon. We emphasize that the string Σ is always defined only up to self-retracing components; i.e. two paths p 1,2 taken by the holon correspond to the same string Σ if p 1 can be obtained from p 2 by eliminating self-retracing components. The distortion fieldσ z i,j (Σ) is uniquely determined by the string configuration and no longer appears as a label of the basis states. Note that two inequivalent states in the over-complete basis can be identified with the same physical state if their holon positions as well as the corresponding distortion fieldsσ z i,j (Σ) coincide. Geometrically, the over-complete space [18] of all strings Σ starting from one given spinon position corresponds to the fractal Bethe lattice. The latter is identical to the tree defined by all possible holon trajectories without self-retracing components. In Fig. 1 (b) this correspondence is illustrated for a simple trajectory taken by the holon. When r = r(Σ) denotes the site on the Bethe lattice defined with the spinon in its origin andĥ † (r) creates the holon in this state, we can formally write the basis states aŝ Linear string theory.-Next we derive the effective Hamiltonian of the system using the new basis states (18). From Eq. (15) we obtain an effective hopping term of the holon on the Bethe lattice, where r, s ∈ BL denotes neighboring sites on the Bethe lattice. This reflects the fact that the system keeps a memory of the holon trajectory. The spin Hamiltonian Eq. (C4) without magnons can be analyzed by first considering straight strings. Their energy increases linearly with their length with a coefficient 4J z S 2 . To obtain the correct energy of the distorted state, we have to include the zero-point energies 4J z S 2 of the holon and 2J z S 2 of the spinon, both measured relative to the energy E cl 0 = −2N J z S 2 of the classical Néel state. Because the state with zero string length = 0 has energy 4J z S 2 , which is 2J z S 2 smaller than the sum of holon and spinon zero-point energies, we obtain a point-like spinon-holon attraction.
The resulting Hamiltonian readŝ (20) Here (r) denotes the length of the string defined by site r on the Bethe lattice andĥ † (0) creates a string with length zero. Because we neglect self-interactions of the string which can arise for configurations where the string is not a straight line, we refer to Eq. (20) as linear string theory (LST). Note that we have written Eq. (20) in second quantization for convenience. It should be noted however, that the spinon and the holon can only exist together and a state with only one of them is not a well defined physical state. Later we will include additional terms describing spinon dynamics. Non-linear string theory.-When the length of the string is sufficiently large, it can start to interact with itself. For example, when a string winds around a loop or crosses its own path, the energy of the resulting state becomes smaller than the value 4J z S 2 used in LST. We can easily extend the effective Hamiltonian from Eq. (20) by taking into account self-interactions of the string. If H J denotes the potential energy of the spin configuration, determined from Eq. (C4) by usingσ z i,j (Σ), we can formally write: The sum in Eq. (21), denoted with a prime, has to be performed over all string configurations Σ which do not include Trugman loops. This is required to avoid getting a highly degenerate ground state manifold, because Trugman loop configurations correspond to strings with zero potential energy. Such states are parametrized by different spinon positions in our over-complete basis from Eq. (18). As will be discussed in detail in Sec. IV, Trugman loops give rise to spinon dynamics, i.e. they induce changes of the spinon position. Their kinetic energy lifts the large degeneracy in the ground state of the potential energy operator.
If, on the other hand, we remove the spinon label from the basis states in Eq. (18) and exclude all string configurations with loops leading to double-counting of physical states in the basis, Eq. (21) corresponds to an exact representation of the single-hole t − J z model. By removing only the shortest Trugman loops and states with zero potential energy while allowing for spinon dynamics, one obtains a good truncated basis for solving the t − J z model.
Using exact numerical diagonalization, the spectrum of H NLST can be easily obtained. Because of the potential energy cost of creating long strings, the holon and the spinon are always bound, see for example Fig. 7 (b) and (c). In Sec. III we will discuss their excitation spectrum and the resulting different magnetic polaron states.

Parton confinement and relation to lattice gauge theory
Finally we comment on the definition of partons in our work and in the context of lattice gauge theories. In the latter case, one usually defines a gauge field on the links of the lattice which couples to the charges carried by the partons. Hence the partons interact via the gauge field, and the question whether they are confined or not becomes a question about the gauge fluctuations [31,99].
In our parton construction so far, we have not specified the gauge field, and partons interact via the string connecting them. Because the string is defined by displaced spins, it can be directly measured, see also Sec. III D, and thus represents a gauge-invariant quantity. An interesting question, which we devote to future research however, is whether a lattice gauge theory can be constructed which has a gauge-invariant field strength corresponding to the string. This would allow to establish even more direct analogies between partons in high-energy physics and holes in the t − J or t − J z model.

D. Strong coupling wavefunction
Because of the single-occupancy constraint enforcing either one spin or one hole per lattice site, see Eq. (4), the spin and charge sectors are strongly correlated in the original t − J z Hamiltonian. Even when a separation of timescales exists, as provided by the condition t J z , no strong-coupling expansion is known for conventional approaches developed to describe magnetic polarons, e.g. for the usual 1/S expansion [17]. This is in contrast to conventional polaron problems with density-density interactions, where strong coupling approximations can provide important analytical insights [100][101][102][103].
In the effective parton theory, the holon motion can be described by a single particle hopping on the fractal Bethe lattice. This already builds strong correlations between the holon and the surrounding spins into the formalism. Because the characteristic spinon and holon time scales are given by 1/J z and 1/t respectively, the magnetic polaron can be described within the Born-Oppenheimer approximation at strong couplings. This corresponds to using an ansatz wavefunction of the form We can solve the fast holon dynamics for a static spinon, and derive an effective low-energy Hamiltonian for the spinon dressed by the holon afterwards. We will make use of this strong-coupling approach throughout the following sections.

III. STRING EXCITATIONS
New insights about the magnetic polaron can be obtained from the simplified LST Hamiltonian in Eq. (20) by making use of its symmetries. Because the potential energy grows linearly with the distance between holon and spinon, they are strongly bound, i.e. the spinon and holon form a confined pair. The bound state can be calculated easily by mapping the LST to an effective one-dimensional problem, see Ref. [15]. After providing a brief review of this mapping, we generalize it to calculate the full excitation spectrum of magnetic polarons including rotational states. We check the validity of the effective LST by comparing our results to numerical calculations using NLST.
A. Mapping LST to one dimension: a brief review The Schrödinger equation for the holon moving between the sites of the Bethe lattice can be written in compact form as [104] Here the linear string potential is given by V = 2J z S 2 (2 − δ ,0 ), E denotes the energy, and z = 4 is the coordination number of the square lattice. In general the wave function ψ(Σ) depends on the index Σ ∈ BL corresponding to a site on the Bethe lattice, or equivalently a string Σ. A useful parameterization of Σ ∈ BL is provided by specifying the length of the string as well as angular coordinates s = s 1 , ..., s with values s 1 = 1...z and s j = 1...z − 1 for j > 1. This formalism is used in Eq. (23) and illustrated in Fig. 8 (a). In Eq. (23) only the dependence on s = s is shown explicitly. Because we started from the LST, the potential V is independent of s. The normalization condition is given by where the sum includes all sites of the Bethe lattice. The simplest symmetric wave functions ψ(Σ) only depend on and are independent of s. We will first consider this case, which realizes the rotational ground state of the magnetic polaron. It is useful to re-parametrize the wave function ψ by writing The normalization for the new wave function φ is given by the usual condition, ∞ =0 |φ | 2 = 1, corresponding to a single particle in a semi-infinite one-dimensional system with lattice sites labeled by .
The Schrödinger equation (23) for the 1D holon wave function φ becomes [15], Away from the origin = 0, the effective hopping constant t * in the 1D model is given by [15,16] The tunneling rate between = 0 and 1, on the other hand, is given by 2t * / √ z − 1 = 2t. Before we move on, we consider the continuum limit of the effective 1D model where φ → φ(x) and x ≥ 0 becomes a continuous variable, see Ref. [15]. This is a valid description in the strong coupling limit, where t J z . For simplicity we will ignore deviations of V from the purely linear form at = 0, as well as the renormalization of the tunneling t * → 2t from site = 0 to 1. As a result one obtains the Schrödinger equation [15] − where the effective mass is m * = 1/2t * , and the confining potential is given by By simultaneous rescaling of lengths, x → λ 1/3 x, and the potential J z → λJ z , one can show that the eigenenergies E in the continuum limit are given by [15,19,47] for some numerical coefficients a n . It has been shown in Refs. [15,19] that they are related to the eigenvalues of an Airy equation. The scaling of the magnetic polaron energy like t 1/3 J 2/3 z is considered a key indicator for the string picture. It has been confirmed in different numerical works for a wide range of couplings [21,25,28,48,105], both in the t − J and the t − J z models. Diagrammatic Monte Carlo calculations by Mishchenko et al. [28] have moreover confirmed for the t−J model that the energy −2 √ 3t is asymptotically approached when J → 0. However, for extremely small J/t on the order of 0.03 it is expected [27] that the ground state forms a ferromagnetic polaron [14] with ferromagnetic correlations developing inside a finite disc around the hole. In this regime Eq. (32) is no longer valid.
Using ultracold atoms in a quantum gas microscope the universal scaling of the polaron energy can be directly probed when J z /t is varied and for temperatures T < J. To this end the super-exchange energy Ĥ J can be directly measured by imaging the spins around the hole. Note that Ĥ J has the same universal scaling with t 1/3 J 2/3 z as the ground state energy at strong couplings.
The excited states of the effective 1D Schrödinger equation (31) correspond to vibrational resonances of the meson formed by the spinon-holon pair, labeled by the vibrational quantum number n. In a semi-classical picture, they can be understood as states where the string length is oscillating in time. Now we generalize the mapping to a 1D problem for rotationally excited states.

B. Rotational string excitations in LST
Within LST the entire spectrum of the magnetic polaron can easily be derived by making use of the symmetries of the holon Hamiltonian on the Bethe lattice. Around the central site, where = 0, we obtain a C 4 symmetry. The C 4 -rotation operator has eigenvalues e iπm4/2 with m 4 = 0, 1, 2, 3 and the eigenfunctions depend on the first angular variable s 1 in the following way: e iπm4s1/2 . So far we assumed that the wave function ψ only depends on the length of the string , which corresponds to an eigenvalue of C 4 which is m 4 = 0.
In addition, every node of the Bethe lattice at > 0 is associated with a P 3 permutation symmetry. The P 3 -permutation operator has eigenvalues e i2πm3/3 with m 3 = 0, 1, 2 and the eigenfunctions depend on the jth angular variable s j , j > 1, in the following way: e i2πm3sj /3 . The symmetric wave function ψ discussed in Sec. III A so far had m 3 = 0 for all nodes.
where we introduced labels (n, m 4 ) denoting the vibrational and the first rotational quantum numbers. Because the dependence of ψ (n,m4) 1,s1 on the first angular variable s 1 is determined by the value of m 4 as explained above, the first term in Eq. (33) becomes Because of the Kronecker delta function δ m4,0 on the right hand side, we see that for m 4 = 0 Eq. (33) becomes V 0 ψ = e i π 2 m4s1 λ φ (n,m4) , s 1 = 0, 1, 2, 3.
Here λ was defined in Eq. (25) and the radial part φ (n,m4) is the solution of the Schrödinger equation (27) - For = 0 we introduced a large centrifugal barrier, preventing the holon from occupying the same site as the spinon. This takes into account the effect of the Kronecker-delta function in Eq. (34), without the need to explicitly deal with the rotational variable s 1 in the wavefunction. Note that the effective 1D Schrödinger equation is independent of m 4 when m 4 = 0, and the same is true for the resulting eigenenergies E (n,m4) .

Higher rotationally excited states
Higher rotationally excited states with non-trivial P 3 quantum numbers m 3 = 0 at some node can be determined in a similar way. Let us consider m 3 = 0 at a node corresponding to a string of length P . The Schrödinger equation at this node reads t 3 s P +1 =1 ψ (n,m3) P +1,s P +1 + tψ (n,m3) Again the first term is only non-zero when m 3 = 0. This can be seen from the dependence of ψ (n,m3) P +1,s P +1 ∝ e i 2π 3 m3s P +1 on the angular variable s P +1 , which yields a Kronecker delta δ m3,0 when summed over s P+1 = 1, 2, 3.
For m 3 = 0, we obtain two sets of independent eigenequations. The first involves only strings of length ≤ P . If it has a non-trivial solution with ψ (n,m3) P = 0 and energy E (n,m3) , there is a second eigenequation involving strings of length ≥ P . In general the second equation cannot be satisfied, because the energy E (n,m3) is already fixed. The trivial choice ψ (n,m3) = 0 for > P does not represent a solution because there exists a nonvanishing coupling to ψ of the rotationally excited string is determined by the Schrödinger equation (27) - (29) for the potential V → V P3 , where In this case there exists an even more extended centrifugal barrier than for the C 4 rotational excitations. It excludes all string configurations of length ≤ P . Note that the effective 1D Schrödinger equation is independent of m 3 when m 3 = 0, and the same is true for the resulting eigenenergies E (n,m3) . Note that the rotationally excited states of the P 3 operator are highly degenerate. Because the wave function vanishes for all ≤ P , there exist 4 × 3 P−1 decoupled sectors on the Bethe lattice where ψ (n,m3) (Σ) = 0. Together with the two choices m 3 = 1 and 2, the total degeneracy D 3 ( P ) becomes C. Comparison to NLST and scaling laws In Fig. 1 (c) and Fig. 9 we show results for the eigenenergies from LST and NLST respectively. As indicated in  9. Excitation spectrum of magnetic polarons in the strong coupling parton theory as a function of Jz/t. Note that the ground state energy was subtracted, which scales as J 2/3 z t 1/3 . We used S = 1/2 and NLST with strings of length up to max = 8. A similar calculation is presented in Fig. 1 using LST with max = 100. The degeneracies of the lowest excited manifolds of states are indicated in circles. Dark blue lines correspond to vibrationally excited states without rotational excitations, i.e. m3 = m4 = 0. Orange lines correspond to purely rotationally excited states with at least one m3 or m4 non-vanishing. No calculations were performed for higher energies (shaded area). The finite gap predicted for small Jz/t by NLST is a finite-size effect caused by the maximal string length max = 8 assumed in the calculations. the figures, we have confirmed for the lowest lying excited states that LST correctly predicts the degeneracies of the low-lying manifolds of states obtained from the more accurate NLST. While these degeneracies are exact for LST, the self-interactions of the string included in the NLST open small gaps between some of the states and lift the degeneracies. In general we find good qualitative agreement between LST and NLST.

Excitation energies
For the energy of the first excited state with rotational quantum number m 4 > 0 and m 3 = 0, we find excellent quantitative agreement between LST and NLST. At small J z /t we note that NLST predicts a larger energy than LST which appears to saturate at a non-zero value when J z /t → 0. This is a finite-size effect caused by the restricted Hilbert space with maximum string length max = 8. Aside from this effect, we obtain the following scaling behavior for all rotationally excited states without radial (i.e. vibrational) excitations (n = 1). In contrast, the excited states with vibrational excitations (n > 1) show a scaling behavior as expected on general grounds from the effective onedimensional Schrödinger equation, see Eq. (32).

String-length distribution
In Fig. 10 we calculate the distribution function p of string lengths for t/J z = 10 well in the strong coupling regime. The comparison between results from NLST with a maximum string length max = 8, and LST calculations with max = 100 shows excellent quantitative agreement for p . Only for the highest excited state considered, with the largest mean string length, we observe some discrepancies at large values of around the cut-off max used in the NLST.
We confirm a strong suppression of p for ≤ P in the rotationally excited states due to the centrifugal barrier. We note that even within the NLST the C 4 rotational symmetry and the P 3 permutation symmetry at = 1 are strictly conserved. The good quantitative agreement for p should be contrasted to the predicted energies, where larger deviations are observed between LST and NLST in the strong coupling regime.

D. String reconstruction in a quantum gas microscope
As explained in Sec. I C the spin configuration around the holon in a t − J z model can be directly accessed [60]. Measurements of this general type are routinely performed in quantum gas microscopy. For example, they have been used to measure the full counting statistics of the staggered magnetization in a Heisenberg AFM, see Ref. [70], and non-local signatures of spin-charge separation, see Ref. [44]. These capabilities should allow imaging of the string attached to the holon, and extract the full distribution function of the string length. This makes quantum gas microscopes ideally suited for direct observations of the different excited states of magnetic polarons. Now we will assume that both spin states and the density can be simultaneously imaged [60].
Knowledge of the spin configuration enables the determination of the distortion fieldσ z i,j . In order to identify the string attached to a holon in a single shot, we now introduce a local measure for the distortion of the Néel state at a given site i, Here, i,j denotes a sum over all bonds i, j to neighboring sites j, where both sites i and j are occupied by spins. Therefore, the happiness θ i assumes integer values between 0 and 4, where 0 (4) corresponds to aligned (anti-aligned) spins on all adjacent bonds.
Since the AFM spin order is maintained along a string without loops, the distortion fieldσ z i,j = 1 is unity on bonds i, j that belong to the string and do not include the holon position. Since spins on the string are displaced with respect to the surrounding AFM, it holdsσ z i,j = −1 if site i is occupied by a spin and belongs to the Retrieving strings in the t−Jz model. We start from a perfect Néel state and create a spinon in the center. Within LST, the distribution of spin patterns is determined by the wavefunction of the holon on the Bethe lattice. The inset shows two snapshots for string lengths = 3 and = 7 where the holon position is marked red. We generated pictures by sampling the three different string length distributions plotted in the figure and choosing random directions of the string. By analyzing the happiness for every spin pattern as described in the text, we retrieve the string for every snapshot. We assume that both spin states and the density can be simultaneously measured. The distribution functions of the lengths of retrieved strings compare well with the sampled distributions, even for long strings.
string and site j is not part of the string. Therefore, the happiness θ i on sites i belonging to the string corresponds to the number of neighboring spins that are also part of the string. For sites j outside of the string, θ j = 4 − N s i , where N s i is the number of neighboring sites that belong to the string. By analyzing the happiness according to these rules, we can start from the holon position and reconstruct the attached string.
The scheme described above allows to directly observe spinons and strings in realizations of the t − J z model with quantum gas microscopes. To mimic this situation, we start from a perfect Néel state and initiate a hole in the center of the system. For a given string length > 0, we first move it randomly to one out of its four neighboring sites. If > 1, we randomly choose one out of the three sites which the holon did not visit in the previous step. This step is repeated until a string of length is created. Thereby, loops are allowed in the string, but the hole cannot retrace its previous path. In Fig. 11, we sample strings according to a given string length distribution function p and retrieve them from the resulting images by examining the values of θ i as described above. Comparison of the sampled distribution, taken from LST, with the retrieved string length distribution shows that even long strings can be efficiently retrieved in the images. Corrections by loops only lead to a small over-counting (under-counting) of short (long) strings on the level of a few percent.

IV. SPINON DISPERSION: TIGHT-BINDING DESCRIPTION OF TRUGMAN LOOPS
In the previous section we were concerned with the properties of the holon on the Bethe lattice. Motivated by the strong-coupling ansatz from Eq. (22) the spinon was treated as completely static so far, placed in the origin of the Bethe lattice. It has been pointed out by Trugman [18] that the magnetic polaron in the t − J z model can move freely through the entire lattice by performing closed loops around plaquettes of the square lattice which restore the Néel order, see Fig. 12 (a). By using the string theory in an over-complete Hilbert space we have not included the effects of such loops so far. Now we show that a conventional tight-binding calculation allows one to include the effects of Trugman loops in our formalism. They give rise to an additional term in the effective Hamiltonian describing spinon dynamics, Here i, j d denotes a pair of two next-nearest neighbor sites i and j on opposite ends of the diagonal across a plaquette in the square lattice; every such bond is counted once. The new term is denoted byĤ T because it describes how Trugman loops contribute to the spinon dynamics. Below we will derive an analytic equation for calculating t T and compare our predictions with exact numerical calculations of the spinon dispersion relation.

A. Trugman loops in the spinon-holon picture
In the spinon-holon theory we use the over-complete basis introduced in Eq. (17). To study corrections in our model introduced by this over-completeness, let us first consider only the LST approximation, see Eq. (20), where the holon is localized around the spinon and the latter has no dynamics. In this case, the over-completeness of the basis leads to very small errors.
To see this, note that the first physical state σ z i,j which can be identified with two inequivalent basis states |x s , Σ , |x s , Σ involves at least six string segments. Such states can be obtained from the so-called Trugman loops [18]: Starting with a holon at spinon position x s and performing one-and-a-half loops around a plaquette, one obtains a state |x s , Σ with string length (Σ) = 6 which corresponds to the same physical state as the string length zero state |x s , Σ with (Σ ) = 0, where x s is a diagonal next-nearest neighbor of x s . This is illustrated in Fig. 12 (a), where after performing the Trugman loop the surrounding Néel state is not distorted. When the confinement is tight, J z t, the wavefunction φ decays exponentially with the string length, and physical states which are represented more than once in the over-complete basis are very weakly occupied. Even when t J z -but still assuming the LST Hamiltonian  FIG. 12. Trugman loops in the spinon-holon picture: (a) When the holon performs one-and-a-half loops around one plaquette (arrows), the spin configuration of the Néel state is restored [18]. This process can be associated with two different configurations of the spinon and the string in the over-complete Hilbert space of the string theory. In the first case the spinon is located at (0, 0) and the string is oriented counterclockwise (blue). In the second case the spinon is located at (1, 1) and the string is oriented clockwise (pink). Trugman loops correspond to a correlated pair tunneling of the spinon and the holon diagonally across the plaquette, for which the holon has to overcome the potential energy barrier shown in (b). At strong couplings, t Jz, this barrier is shallow. The resulting effective spinon hopping element can be calculated by a tight-binding approach, where the difference in potential energiesĤNLST −ĤLST between the full non-linear and linear string theory are treated perturbatively (c). This is similar to the tight-binding calculation of the hopping element of an electron between two atoms. Like the string theory the latter relies on using an over-complete basis, with one copy of the entire Hilbert space for each atom.
-the fraction of physical states with multiple representations in the over-complete basis is small, see discussion above Eq. (17). The reason is essentially that specific loop configurations need to be realized out of exponentially many possible string configurations. Now we consider the effect of non-linear corrections to the string energy, as described in Eq. (21). We can loosely distinguish between two types of corrections: (i) for strings without loops, attractive self-interactions between parallel string segments can lower the string energy, and (ii) for strings with Trugman loops, the potential energy can vanish completely. While the first effect (i) only leads to small quantitative corrections of the spinon-holon energy, the second effect leads to a large degeneracy within the over-complete basis and needs to be treated more carefully.
When two states |x s , Σ and |x s , Σ in the overcomplete basis have zero potential energy and their holon positions coincide, their distortion fields are identical, σ z i,j (Σ) = σ z i,j (Σ ) and they correspond to the same physical state. By this identification we notice that the effect of Trugman loops is to introduce spinon dynamics: consider starting from a string length zero state around a spinon at x s . By holon hopping a final state in the over-complete basis can be reached which can be identified with another string length zero state but around a different spinon position x s . Our goal in this section is to start from degenerate eigenstates of the LST and describe quantitatively how the corrections by NLST, H NLST −Ĥ LST , introduce spinon dynamics. This can be achieved by a conventional tight-binding approach.
In the case when t < J z the Trugman loop process is strongly suppressed because it corresponds to a 6th order effect in t and the holon has to overcome an energy barrier of height 12J z S 2 , see Ref. [18] and Fig. 12 (b). Although these perturbative arguments based on an expansion in t/J z no longer work when t J z , we will show that Trugman loop processes only lead to small corrections t T J z , a small fraction of J z . The key advantage of using tight-binding theory is that the potential energy ∼ J z rather than the holon hopping ∼ t is treated perturbatively.

B. Tight-binding theory of Trugman loops
To explain how tight-binding theory allows us to take into account the over-completeness of the basis used in LST, see Sec. II C 2, we draw an analogy with conventional tight-binding calculations for Bloch bands. We start by a brief review of the tight-binding approach and explain how, implicitly, use is being made of an overcomplete basis. See e.g. Ref. [106] for an extended discussion of conventional tight-binding calculations.

Tight-binding theory in a periodic potential
Consider a quantum particle moving in a periodic potential W (x) = W (x + a). We assume that the lattice W (x) has a deep minimum at x j within every unit cell [ja, (j + 1)a], where the particle can be localized. This is the case for example if W (x) = jW (x − x j ) corresponds to a sum of many atomic potentialsW (x − x j ) created by nuclei located at positions x j . Similarly, the NLST potential is periodic on the Bethe lattice, although the geometry is much more complicated in that case.
The idea behind tight-binding theory is to solve the problem of a single atomic potentialW (x−x 0 ) first. This yields a solutionw(x − x 0 ) localized around x 0 . Then one assumes that the orbitalw(x − x 0 ) is a good approximation for the correct Wannier function w(x − x 0 ) defined for the full potential W (x), see Fig. 12 (c). The tight-binding orbitalw(x − x 0 ) defined by the potential W (x − x 0 ) around x 0 is similar to the holon state defined by the LST HamiltonianĤ LST around a given spinon position x s . The potential energy mismatch W (x) −W (x) can now be treated as a perturbation. Most importantly, it induces transitions between neighboring orbitalsw(x − x j ) andw(x − x j±1 ). This leads to a nearest neighbor tightbinding hopping element For this perturbative treatment to be valid, it is sufficient to have a small spatial overlap between the two neighboring orbitals, The energy difference W (x − x j ) −W (x − x j±1 ) on the other hand can be sizable. In practice, the most common reason for a small wavefunction overlap ν 1 is a high potential barrier which has to be overcome by the particle in order to tunnel between two lattice sites. In this case, once the barrier becomes too shallow, ν becomes sizable and the tightbinding approach breaks down. As another example, consider two superconducting quantum dots separated by a dirty metal. Even without a large energy barrier, the wavefunction of the superconducting order parameter decays exponentially outside the quantum dot due to disorder [107,108]. This leads to an exponentially small wavefunction overlap ν 1 and justifies a tight-binding treatment.
Similar to the case of the string theory of magnetic polarons, the effective Hilbert space used in conventional tight-binding calculations is over-complete. To see this, note that the tight-binding Wannier functionw(x − x j ) corresponding to lattice site j is defined on the space H j of complex functions mapping R → C with site j in the center. Assuming that every such Wannier function is defined in its own copy of this Hilbertspace H j , we obtain -by definition -that the resulting tight-binding wavefunctions are mutually orthogonal. But in reality, the physical Hilbert space H phys consists of just one copy of the space of all functions mapping R → C, and after identifying all states in H j with a physical state in H phys , the resulting tight-binding Wannier functions in H phys are no longer orthogonal in general: i.e. ν = 0. As long as ν 1, they still qualify as good approximations for the true, mutually orthogonal Wannier functions w(x − x j ), which justifies the tight-binding approximation.

Tight-binding theory in the spinon-holon picture
We can now draw an analogy for a single hole in the t − J z model. The true physical Hilbert space corresponds to all spin configurations and holon positions, see Eq. (16). The over-complete Hilbert space is defined by copies of the Bethe lattice, each centered around a spinon positioned on the square lattice. As a result of Trugman loops, certain string configurations can be associated with two spinon positions, see Fig. 12 (c).
We start by describing the tight-binding formalism for the Trugman loop hopping elements. When t J z , large energy barriers strongly suppress spinon hopping, see Fig. 12 (b). As a result, the overlap of the holon state corresponding to a spinon at x s with the holon state bound to a spinon at x s + r is small, and the tight-binding approximation is valid. Here r denotes a vector connecting sites from the same sublattice and n is the vibrational quantum number of the holon state. Note that in the definition of ν in Eq. (45), two states from the overcomplete basis are assumed to have unit overlap if they correspond to the same physical state. The holon wave functions on the Bethe lattice for different spinon positions x s are the equivalent of the tightbinding orbitalsw(x − x j ) corresponding to lattice sites x j discussed above. We will show below that the condition Eq. (45) remains true even when the barrier becomes shallow, J z t, and that even in this regime the following tight-binding calculation is accurate.
Transitions between states associated with different spinon positions x s and x s + r are induced by the potential energy mismatchĤ NLST −Ĥ LST , see Fig. 12 (b). This takes the role of W (x)−W (x) from the conventional tight-binding calculation. Accordingly, the tight-binding spinon hopping element is given by t T (n) = ψ n (r s + e x + e y )|Ĥ NLST −Ĥ LST |ψ n (r s ) . (46) Here we have set r = e x + e y which corresponds to the simplest Trugman loop process inducing spinon hopping diagonally across the plaquette. A more precise expression is obtained by noting that for long stringŝ H NLST −Ĥ LST also leads to longer-range spinon hopping, which should be counted separately. Below we will discuss in detail which contributions we include to obtain an accurate expression, see Eq. (50), for the diagonal Trugman loop spinon hopping t T from Eq. (42).
To understand the meaning of Eq. (46), consider the simplest configuration |x s , Σ for which (Ĥ NLST − H LST )|x s , Σ = 0. It has a string length of = 3 and can be reached by performing three quarters of a loop around a plaquette, starting with either a spinon in the lower left corner and going counter-clockwise, or a spinon in the top right corner and going clockwise. This is illustrated in Fig. 12 (a) and (c). Now we consider the case of a shallow potential barrier J z t and show that the overlap ν in Eq. (45) is exponentially suppressed, even for vibrationally highly excited states with a large average string length. The overlap ν can be written as a sum over terms of the form ψ * ψ 6− , where the wavefunction ψ defined on the Bethe lattice is related to the effective one-dimensional wavefunction φ by Eq. (25) with a typical extent given by L 0 . For simplicity we approximate |φ | 2 ≈ 1/L 0 for < L 0 and φ = 0 otherwise. Following Bulaevskii et al. [15] we equate the average kinetic and potential energies in the effective Schrödinger equation (31), 1/2m * L 2 0 = 4J z S 2 L 0 , to obtain the following estimate for the typical length of the string connecting the spinon to the holon. In the second step we used S = 1/2 and z = 4 for the 2D anti-ferromagnet. Because the Bethe lattice has a fractal structure with z − 1 new branches emerging from each node, the amplitude |ψ | 2 is reduced by the additional exponential factor λ 2 as compared to |φ | 2 , see Eq. (25). It holds |ψ | 2 ≈ (z − 1) 1− /(4L 0 ), from which we obtain using the estimate for L 0 from above. In terms of the original t−J z model, the exponential suppression of |ψ | 2 with can be understood as a consequence of the exponentially many spin configurations that can be realized by the motion of the hole. From Eq. (48) we see that already for = 3 we obtain |ψ 3 | 2 ≈ 0.023 (J z /t) 1/3 1 and thus ν 1. Even for excited string states the exponential suppression of the amplitude |ψ | 2 by λ 2 on the Bethe lattice leads to |ψ 3 | 2 ≈ 0.028|φ 3 | 2 1. Therefore the perturbative treatment of Trugman loops is justified for the entire range of parameters t/J z . We demonstrate this by an exact numerical calculation for a closely related toy model which is presented in Appendix D.
Next we provide additional physical intuition why the tight-binding approach works even in the regime when the potential barrier is shallow, J z t. In this limit, the kinetic energy of the holon is minimized by delocalizing the holon symmetrically over all possible z − 1 directions at each node on the Bethe lattice. This leads to the exponential decay of the holon wavefunction ψ with , which is equal for all directions on the Bethe lattice. It also yields a zero-point energy of −2 √ z − 1t when J z = 0. In order to minimize its energy further, the holon could make use of the non-linearity of the string potential described byĤ NLST : By occupying preferentially the directions on the Bethe lattice corresponding to Trugman loops, the average potential energy is lowered. This mechanism is very ineffective, because it requires localizing the holon on the particular set of states corresponding to the Trugman loops. This costs kinetic energy of order t: The average kinetic zero-point energy of a holon moving along one fixed direction is given by −2t instead of −2 √ z − 1t. Therefore, in the limit t J z , the kinetic term in the Hamiltonian dictates a symmetric distribution of the holon over all possible directions, which allows to treat perturbatively the effect of the reduced string potential for the very few specific directions defined by the Trugman loop string configurations. C. Application: Spinon hopping elements in the t − Jz model Next we apply the tight-binding theory of Trugman loops introduced above and derive a closed expression for the effective spinon hopping elements using Eq. (46).

Contributing string configurations
For every plaquette there exist two Trugman loops contributing to t T : a clockwise and a counter-clockwise one. In Fig. 12 (b) we show the energies V NLST for states along one of these loops. For strings of length 1 ≥ 3, measured from the first spinon position x s , they differ from the expression V = 4J z S 2 assumed in LST by for 1 = (3, 4, 5, 6). These configurations overlap with strings of length 2 = 6 − 1 , measured from the second spinon position x s ± e x ± e y on the opposite side of the plaquette in the square lattice. Using Eq. (46) these states contribute an amount 2 In addition there are overlaps for longer strings with 1 + 2 > 6. Let us start from a configuration with strings of lengths The only missing set of strings with 1 + 2 > 6 consists of cases where the first string (with length 1 measured from the first spinon) includes the position of the second spinon on the Bethe lattice. This situation corresponds to (0) 1 = 6 and (0) 2 = 0 using the notation from the previous paragraph. In this case there exist (z − 1) n configurations of this type, for most of which one obtains a mismatch of potential energies of δV 7 = −24J z S 2 . Along certain directions the magnitude of the energy mismatch between LST and NLST is larger, and this corresponds to a case where a second Trugman loop follows the first one. Because it gives rise to longer-range spinon hopping, we will neglect these contributions in the following and take into account the constant energy mismatch of δV 7 = δV 6 + 2J z S 2 in the tight-binding calculation of t T . This class of strings thus contributes 2δV 7 ∞ n=1 (z − 1) n λ n φ * n λ 6+n φ 6+n to t T .

Trugman loop hopping elements
Summarizing, we obtain the following expression for the tight-binding Trugman loop hopping element within LST. Note that the exponential factor (z − 1) n is canceled by powers of (z −1) appearing in the coefficients λ relating the wavefunction on the Bethe lattice to the effective one-dimensional wavefunction φ . In Fig. 14 we used Eq. (50) to calculate the diagonal spinon hopping elements t T for the lowest vibrational string states n = 1, ..., 5. For small J z /t 1 we find the largest hoppings in units of J z . For the ro-vibrational ground state of the string, n = 1, the Trugman loop hopping t T < 0 is always negative. This is due to the fact that we can choose φ > 0 in the ground state, such that all overlaps in Eq. (50) contribute a negative amount to t T . For vibrationally excited states φ has nodes where the wave function changes sign, thus reducing the Trugman loop hopping element. This effect explains the oscillatory behavior of t T (J z ) observed as a function of J z for vibrationally excited states.

Comparison to exact diagonalization
To test the strong coupling description of magnetic polarons in the t − J z model, we compare our prediction for the dispersion relation to exact numerical simulations in small systems. First we study the shape of the dispersion in a 4-by-4 lattice, followed by a discussion of the magnetic polaron bandwidth W which has been calculated in larger systems in Refs. [4,109].
In the parton theory of magnetic polarons, we use a product ansatz for the strong coupling wave function where the holon is bound to the spinon, see Eq. (22). As a result the momentum k of the magnetic polaron is entirely carried by the spinon; Its dispersion relation is determined from the effective spinon dispersion which is obtained from Eq. (42): In Fig. 15 our results for the magnetic polaron dispersion relation are shown for t/J z = 0.2. We note that the ground state energy of the magnetic polaron on the 4 × 4 torus is rather well reproduced by LST. The minimum of the dispersion relation is predicted at (0, 0) and (π, π) by LST. These two states are exactly degenerate in an infinite system because the Néel state breaks the sublattice symmetry. In the exact numerical calculation, the eigenstates with the two lowest energies are also located at momenta (0, 0) and (π, π). We note that the numerical results show an energy difference ∼ 0.25J z between (0, 0) and (π, π). This is a finite-size effect, indicating that the 4 × 4 lattice is not large enough to break the discrete translational symmetry of the t − J z model and form a Néel state.
The shape of the dispersion in the 4×4 system deviates somewhat from the expectation in the infinite system. The predictions from LST are compared to exact diagonalization (ED) for a small system with periodic boundary conditions. The solid black line is the LST prediction for an infinite system with a single hole. The dashed line is obtained by using the spinon dispersions for a 4 × 4 system, where the spinon hopping elements tT are calculated in the infinite system. Note that an overall energy shift was added to the dashed line to simplify a direct comparison of the predicted bandwidths.
In particular, the energies at (π, 0) and (±π/2, ±π/2) are equal in this case. It is understood that this is a consequence of a particular symmetry of the 4 × 4 system [105]. We can also easily understand this from the tightbinding theory for spinon hopping: On a 4 × 4 cylinder there exists an additional Trugman loop leading to next nearest neighbor hopping along the lattice direction. In this case the holon moves around the torus one-and-a-half times in a straight manner. As a result we have to add an additional term t T [cos(2k x ) + cos(2k y )] to the spinon dispersion relation in Eq. (51), which leads to an exact degeneracy of (π, 0) and (±π/2, ±π/2). This dispersion relation is plotted as a dashed line on top of the data in Fig. 15. Note that we added a small overall energy shift to the dashed line to simplify a comparison with the shape of the exact dispersion relation.
Comparison of the exact bandwidth with the LST calculation shows sizable quantitative deviations for the small 4 × 4 lattice. We expect that this is mostly due to finite size effects, as indicated for example by the absence of a degeneracy between momenta (0, 0) and (π, π) as discussed above. We found that this effect becomes even more dramatic for smaller values of J z /t. For example, at J z /t = 0.1 the energy difference between (0, 0) and (π, π) is comparable to the entire bandwidth.
To study finite-size effects more systematically, we compare the bandwidth W obtained from our tightbinding calculation to exact numerical results in larger systems obtained by Poilblanc et al. [109] and Chernyshev and Leung [4]. From Eq. (51) we see that the bandwidth can be directly related to the Trugman loop hop- ping, W = 8|t T |. In Fig. 16 we find excellent agreement between our tight-binding calculation and exact numerics in the regime where finite-size effects are small, for t 2J z . For larger values of t/J z , the numerical results begin to depend more sensitively on system size N . For the largest systems the behavior is even nonmonotonic with N . The tight-binding result is closest to the data obtained for the largest system of 32 sites solved by Chernyshev and Leung [4].
As a function of t/J z , we also observe non-monotonic behavior. For small t/J z the Trugman loop hopping is exponentially suppressed due to the presence of a large energy barrier. The amplitude |t T | (in units of t) reaches a maximum around t/J z ≈ 6. For larger t/J z the strings become very long, leading to a saturation of t T at a fraction of J z and thus a slow decay of |t T |/t. The numerics for finite-size systems shows a second increase of the bandwidth beyond some critical hopping t > t c . For the 4×4 system this value corresponds to the point where we observed that the bandwidth becomes comparable to the energy splitting between momenta (0, 0) and (π, π). This signals that the translational symmetry is not broken, which is a strong indication for a finite-size effect. Indeed, the critical value t c /J z quickly increases for larger system sizes. This is expected from the dependence of the average string length on t/J z .

Contributions from longer loops
Longer strings can generate higher-order loops, which gives rise to further-neighbor spinon hopping processes preserving the sublattice index. To a good approximation, they can be neglected because they are exponentially suppressed by the string length. Consider for example the simplest Trugman loop which renormalizes the next-neighbor hopping of the spinon linearly along the lattice direction. It involves a string of length = 10, four units longer than for the simplest Trugman loop. This already reduces the corresponding tight-binding element by a factor of 3 −4 ≈ 0.012. Similarly, there are four loops involving two plaquettes which renormalize the diagonal spinon hopping. They can change the value of t T on the level of 4 × 3 −4 ≈ 5%.

V. DYNAMICAL PROPERTIES OF MAGNETIC POLARONS
In this section we apply the strong coupling parton theory of magnetic polarons introduced in the previous section. We consider ultracold atom setups and calculate the dynamics of a single hole in the AFM for various experimentally relevant situations. We discuss how such experiments allow to test the parton theory of magnetic polarons. In the following, it will be assumed that the temperature T J is well below J where corrections by thermal fluctuations are negligible.

A. Effective Hamiltonian
The basic assumption in the strong coupling parton theory is that the holon can adiabatically follow the spinon dynamics. Together they form a magnetic polaron, which can be described by an operatorf † j,n creating a magnetic polaron at site j. Formally we can writê whereŝ † j creates a spinon on site j of the square lattice; |ψ (n) BL denotes the wave function of the spinon-holon bound state on the Bethe lattice with ro-vibrational quantum number n.
The effective Hamiltonian of the magnetic polaron in the t − J z model reads, The energy E(n) of the spinon-holon bound state, as well as the Trugman loop hopping element t T (n) (Eq. (46)) depend on the quantum number n, see Figs. 9 and 14.
Inter-band transitions, where n is changed by a hopping process of the spinon, can be safely neglected in the strong coupling limit.

B. Ballistic propagation: spinon dynamics
We begin by studying coherent spinon dynamics. As an initial state we consider a spinon localized in the origin of the square lattice, and the holon in its ro-vibrational ground state, A protocol how this state can be prepared experimentally is presented below.

Spinon expansion
The dynamics of the magnetic polaron can be characterized by the root-mean square (rms) radii of the density distributionsn s,h (x) of the spinon and the holon respectively, Both quantities can be measured using quantum gas microscopes. In such experiments the holon density is obtained directly by averaging over sufficiently many measurements of the hole position. By imaging the spin configuration [60] the position of the spinon as well as the string configuration can also be estimated, see discussion at the end of Sec. III C. In Fig. 17 the rms radii are shown as a function of the evolution time. Because the spinon is fully localized initially, it expands ballistically and its rms radius grows linearly in time. The initial holon extent, in contrast, is determined by the spinon-holon wave function on the Bethe lattice. When the radius of the spinon density distributionn s (x) becomes comparable to that of the holon distribution, the hole densityn h (x) also starts to expand ballistically with the same velocity v s as the spinon distribution. In the t − J z model the velocity v s is directly proportional to the Trugman loop hopping element t T of the spinon, v s ≈ 3|t T |, see Fig. 18. Experimentally this relation allows a direct measurement of the Trugman loop hopping element.

Preparation of the initial state
Experimentally the initial state (54) can be prepared by first pinning the hole on the central site of the lattice using a localized potential of strength g. When g is decreased slowly compared to the energy gap to the first vibrationally excited state, the holon adiabatically follows its ground state on the Bethe lattice.
When the potential is lowered quickly compared to the spinon hopping, during a time τ 1/|t T |, we can assume that the spinon remains localized on the central site during the quench. Because of the symmetry around the initial spinon position, the rotational quantum numbers of the holon eigenstate on the Bethe lattice do not change. They remain trivial as in the fully localized initial state. Thus, the central trapping potential can only couple different vibrationally excited states of the holon on the Bethe lattice during the adiabatic dynamics. In the limit when g = 0 these eigenstates are separated by an energy gap ∆E ∼ t 1/3 J 2/3 z J z . Hence one can adiabatically prepare the ro-vibrational ground state by choosing 1/|t T | τ 1/∆E. During this time τ required for the preparation scheme, the spinon essentially remains localized.
In Fig. 19 we calculate the overlap of the initially localized holon state with the vibrational eigenstates of the magnetic polaron in the absence of a pinning potential (g = 0) and using LST. As shown in the inset of the figure, we start from g = −5t and decrease the magnitude of the potential linearly over a time τ = 3/t. This is sufficiently fast to neglect spinon dynamics and, as explained above, this also justifies to ignore rotationally excited eigenstates. Fig. 19 demonstrates that even without optimizing the adiabatic preparation scheme, large overlaps with the ro- Adiabatic preparation of the ro-vibrational ground state of the holon. We start from a hole which is localized in the center by a pinning potential of strength g. The latter is decreased over a time τ = 5/t (inset) and the overlap of the holon state with the vibrational eigenstates of the magnetic polaron are calculated. We used LST for S = 1/2 and neglected spinon dynamics. This is justified for the considered value Jz/t = 0.2 in the strong coupling regime.
vibrational ground state of the holon (around 80% in this case) can be readily achieved using this method.

C. Bloch oscillation spectroscopy of the magnetic polaron spectrum
In the last section we explained how the slow dynamics of the hole density distribution allows to experimentally measure spinon properties (specifically t T of the vibrational ground state). Now we present a scheme allowing to study the rotational excitations of the holon experimentally and measure the corresponding excitation energies discussed in Sec. III.
As in the previous section, our starting point is a magnetic polaron in its ro-vibrational ground state. For simplicity we assume that the spinon is localized initially. As shown in Sec. V B 2 this state can be adiabatically prepared in experiments with a quantum gas microscope.
To populate the first excited state, which has a nontrivial rotational quantum number, we apply a force to the fermions,Ĥ This breaks the rotational symmetry around the initial spinon position. In a system at half filling with a single hole, the force corresponds to a potentialĤ F = −F ·x h acting on the hole at positionx h . A similar situation has been considered before by Mierzejewski et al. in Ref. [85], where interesting transport properties of the magnetic polaron have been predicted in the regime of a strong force. The term Eq. (56) in the Hamiltonian can be easily realized for ultracold atoms using a magnetic field or an optical potential gradient. Here we study the opposite limit of a weak force, F J z . The frequency of Bloch oscillations is given by ω B = aF , where a is the lattice constant. We assume that ω B is small compared to the gap ∆E to the first ro-vibrational excited state of the holon, We consider a situation where the force is suddenly switched on at time t 0 . In the subsequent time evolution, the population of the first excited state begins to oscillate with a high frequency given by ω = ∆E (we set = 1). This also manifests in oscillations of the density distribution of the holon, with the same frequency ω. Because the density distribution of the hole can be directly measured, this allows to extract the excitation energy of the first rotationally excited state of the magnetic polaron experimentally.
In Fig. 20 (a) we calculate the rms radius x 2 h 1/2 along the direction x of the applied force F = F e x , as a function of time. For various values of J z /t we observe clear, although not perfectly harmonic, oscillations. To extract their frequency and amplitude we performed the following fit to our numerical results, The fit parameters are the offset x 0 , the amplitude δx, the frequency ω and phase φ 0 , and the decay rate γ of the oscillations. In Fig. 20 (b) we plot the frequencies ω extracted from fits to the holon distribution for times t 0 , ..., t 0 + 300/t. These data points strongly depend on J z /t and are in excellent agreement with the energy gap to the first excited state of the magnetic polaron. We checked that the deviations from the exact value of the excitation energy ∆E are of the order of the Bloch oscillation frequency, ∆E − ω ≈ ω B . The decay rates γ ω B were negligible for the case F a = 0.05 considered in Fig. 20.
In Fig. 20 (c) we show the amplitude δx and the offset x 0 as a function of J z /t. We observe that the experimentally most interesting regime corresponds to cases where J z is not much larger than the Bloch oscillation frequency ω B . Here the amplitude is sizable and the signal-to-noise ratio required in an experiment is not too small. In Fig. 20 we have kept the force F constant. In order to obtain larger amplitudes of the oscillations it can also be tuned to larger values when J z /t is large.

Effect of spinon dynamics
Our calculations so far were performed using NLST in the spinon frame, see Fig. 20. Now we discuss the effects of spinon dynamics on the density distribution of the hole. First we note that the force also acts on the spinon: When the spinon moves from site i at position x i to site j at x j , the holon can follow adiabatically. The center of mass of the holon distribution is always given by the spinon position. When it changes from x i to x j , this leads to an overall energy shift F · (x i − x j ) which corresponds to a force F acting on the spinon.
In the t − J z model the spinon dynamics are slow, with a velocity set by the Trugman loop hopping element t T J z . This allows to realize a regime where the force F t T /a is large compared to t T . In this case the spinon is localized in a Wannier Stark state and it cannot move along the direction e x of the force. At the same time the force can be chosen to be small compared to the excitation energy ∼ J z of the first rotational state. This is necessary in order to observe coherent oscillations, see Eq. (57). Finally, the motion of the spinon in the direction e y orthogonal to the force has no effect on the rms radius of the hole distribution in x-direction which is calculated in Fig. 20 (a).
If the spinon hopping is comparable to the force, as ex- pected for isotropic Heisenberg interactions between the spins, the spinon will undergo Bloch oscillations. When Eq. (57) is satisfied, their frequency ω B is slow compared to the oscillations ω of the holon distribution. The resulting density of the hole in the laboratory frame, which can be directly measured experimentally, is a convolution of the spinon density and the holon distribution function. Thus by performing a Fourier analysis of the hole distribution, both frequency components ω B and ω can be extracted and analyzed separately.
D. Far-from equilibrium dynamics: releasing localized holes So far our analysis was restricted to situations close to equilibrium, where mostly the ro-vibrational ground state of the magnetic polaron was populated. Now we consider the opposite limit where the hole is initially localized on the central site and the deep pinning potential is suddenly switched off at time t 0 . This corresponds to the initial state In this case many different vibrational states n of the magnetic polaron are populated (with amplitudes ψ n ).
Recall thatĥ † (0) creates a holon in the center of the Bethe lattice defined in the spinon frame, see Sec. II C 2.
From the symmetry of the initial state underĈ 4 rotations of the holon around the spinon, it follows that no rotational excitations are created. While we focus on the zero temperature case here, the same problem was discussed in the opposite limit of infinite temperature and for J z = 0 in Refs. [82,89]. Here the vibrational quantum number n acts like an effective band index. Note that the rotational symmetry of the Bethe lattice around the spinon position is unbroken. In particular, we calculate the density of the hole n h (x) in the laboratory frame which is given by Hereĥ † j is the holon operator in the laboratory frame corresponding to site j on the square lattice. To evaluate the last expression we need the matrix elements which we calculated using Monte Carlo sampling over states on the Bethe lattice. Numerical results.-In Fig. 21 (a) we show the density distribution of the hole after the quench. It can be understood as a convolution of the spinon and the holon wavefunctions, whose density distributions are calculated in Fig. 21 (b) and (c). Note however that interference terms between terms from different bands n = n , included in Eq. (61), are taken into account in Fig. 21 (a).
In Fig. 22, we analyze the rms radius of the hole density distribution in the laboratory frame, for the same situation as in Fig. 21. We consider a large value of J z /t = 0.5 but checked that the qualitative behavior is identical at smaller values of J z /t. Our theoretical approach was benchmarked in Fig. 3. There we compared our results to time-dependent quantum Monte Carlo calculations at short-to-intermediate times and found excellent agreement with predictions by the LST.
As a main result of Fig. 22, we observe a clear separation of spinon and holon dynamics. At short times, the hole distribution is entirely determined by fast holon dynamics on the Bethe lattice, with a characteristic time scale ∼ 1/t, see also Fig. 3. Additionally, for short times up to ∼ 20/t the rms radii x 2 h 1/2 and (x h −x s ) 2 1/2 calculated in the lab and spinon frames respectively almost coincide in Fig. 22.
At long times, we observe a ballistic expansion of the hole distribution. This behavior is similar to the coherent spinon dynamics discussed in Sec. V B, see Fig. 17. Similar to that case, we observe that the rms radius of the hole is equal, up to a constant offset, to the rms radius of the spinon distribution x 2 s 1/2 at long times. The characteristic velocity is determined by the Trugman loop hopping element t T of the spinon. It is much smaller than the bare hole hopping t, explaining the large separation of time scales observed in Fig. 22.
The separation of characteristic spinon and holon timescales is a hallmark of the parton theory of magnetic polarons. It can be understood as a precursor of true spin-charges separation: on short length scales (short times), the holon behaves almost as if it was free, resembling the physics of asymptotic freedom known from quarks in high-energy physics. Only on long length scales (long times) the dynamics of magnetic polarons become dominated by the slow spinon, and it becomes apparent that the spinon and the holon are truly confined.
Another key indicator for the string theory of magnetic polarons is the scaling of the energy spacings between low-energy vibrational excitations with the nontrivial power-law t 1/3 J 2/3 z , see Ref. [15] and Sec. III. It is tempting to search for the same universal power-law in the far-from equilibrium dynamics considered here. To this end, we analyzed the position t 1 in time of the first pronounced maximum of the rms radius of the hole. For example, in Fig. 22 it is located around t 1 ≈ 3.5/t. By repeating this procedure for various ratios J z /t we found a trivial power-law t 1 ∝ J z .
To understand why we did not obtain the non-trivial power-law t 1/3 J 2/3 z we note that in the initial state the hole was localized on a single lattice site. The continuum approximation leading to the Schrödinger equation (31) is thus not justified and lattice effects become important. In particular this violates the scale invariance of the Schrödinger equation which gives rise to the characteristic J 2/3 z power-law. If on the other hand one starts from a low-energy state where the holon is distributed over several lattice sites, universal charge dynamics can be observed. Such situations can be realized experimentally by adiabatically releasing the hole from a pinning potential, see Sec. V B 2. A similar situation has been studied by Golez et al. in Ref. [110] where a different quench was considered starting from the ground state of the magnetic polaron. By rescaling the time axes by a factor of t 1/3 J 2/3 z the authors of Ref. [110] demonstrated a collapse of their data at various values of J z /t to a single universal curve for sufficiently short times.

E. Spectroscopy of magnetic polarons
The most direct way to confirm the t 1/3 J 2/3 z power-law describing the energy of vibrationally excited states is to measure the spectral function A(ω, k) of the magnetic polaron. Self-consistent Greens function calculations by Liu and Manousakis [48] as well as diagrammatic Monte Carlo calculations by Mishchenko et al. [28] assuming isotropic Heisenberg couplings between the spins have been performed within the t − J model. They showed strong evidence that the spectrum consists of several broad peaks above the magnetic polaron ground state, which have an energy spacing scaling like t 1/3 J 2/3 .
The spectral function can be measured in solid state systems using ARPES. Similarly, for ultracold fermions radio-frequency [111,112], Bragg [113] and lattice modulation [114] spectroscopy measurements have been carried out. In quantum gas microscopes the spectrum can also be measured with full momentum and energy resolution [83] by modulating the tunneling rate of ultracold atoms into an empty probe system.
Here we briefly discuss the properties of the magnetic polaron spectrum in the strong coupling regime where t J z . Because the dispersion of the magnetic polaron is determined by the small Trugman loop hopping t T , the momentum dependence of the spectrum can be neglected, i.e. A(ω, k) ≈ A(ω). This allows for a purely We have calculated the spectral weights Zn and energies En of all eigenstates n and broadened them by hand for clearer presentation. LST predicts a series of approximately equally spaced peaks, which evolve into a broad continuum if stringinteractions are included in the NLST. The gap from the magnetic polaron peak at the lowest energy to the first excited states with appreciable spectral weight scales with the nontrivial power-law t 1/3 J 2/3 z characteristic of the string theory. local measurement of the spectral function, which significantly simplifies implementations in a quantum gas microscope [83,115].
In general, a series of peaks is expected which correspond to vibrationally excited states of the string. As shown in Fig. 23 and previously in Ref. [19] this is indeed what LST predicts. In the NLST substantial corrections are present at high energies which lead to a broadening of the spectral lines. At low energies we observe a gap both within LST and NLST which scales like t 1/3 J 2/3 z . In a measurement of the spectral function a fermion is removed from the system and the spinon and the holon are always created at the same site. As a consequence, the initial state of the holon on the Bethe lattice is rotationally invariant around the spinon position. This gives rise to a selection rule and leads to vanishing spectral weight of string states with rotational excitations. Thus the latter cannot be observed in A(ω). This explains why the spectral gap in Fig. 23 corresponds to the gap to vibrational excitations, in contrast to the rotational states observed in the Bloch oscillation spectroscopy of Sec. V C. In that case the rotational symmetry was explicitly broken by the applied force.

VI. EXTENSIONS TO THE t − J MODEL
So far we have formulated the strong coupling parton description of magnetic polarons specifically for the t − J z model. As an important outcome of our work, the approach can be generalized to models including quantum fluctuations of the surrounding spin system. Most importantly, we can apply it to describe the dynamics of holes in the t − J model, which is believed to play a central role in the understanding of high-temperature superconductivity. In the following we summarize how such generalizations work and present first results of this method relevant for current experiments with quantum gas microscopes. A more detailed discussion will be provided in forthcoming papers.
Our work on the t − J z model enables two possible extensions to the t − J model. The first uses the generalized 1/S expansion, see Sec. II B and Appendix C. It is based on the idea to start from a classical Néel state which is distorted by the holon motion as in the case of the t − J z model. Transverse couplings between the spins are then included by applying linear spin wave theory around the distorted classical state. Even in the presence of quantum fluctuations this approach allows to introduce spinons and strings in the theory. The properties of the resulting mesons, for example the string tension, are renormalized in this case due to polaronic dressing by magnon fluctuations.
Here we discuss a second extension, which is based on an analogy with the squeezed space picture of the one-dimensional t − J model [43,44,116]. Instead of labeling the basis states by the eigenvalues of the spin operatorsŜ z j on lattice sites j, we keep track of the holon trajectory Σ defined on the Bethe lattice. The motion of the hole changes the original positions of the spins, and their quantum state is defined in a pure spin system without doping on the original square lattice; we identify the latter with the two-dimensional analogue of squeezed space. Similar to Eq. (1) in the case of the t − J z model, the Hilbert space used in the parton theory of a single hole in the t − J model becomes The first two terms H s and H Σ describe the geometric string connecting the spinon at one end with the holon at the other end, i.e. the meson degrees of freedom. The last term H mag includes quantum fluctuations in the surrounding spin system -e.g. magnons in the case of a quantum Néel state. At strong couplings, when t J, we can make a product ansatz to describe the magnetic polaron in the t − J model, |ψ mag.pol. = |ψ spinon ⊗ |ψ holon ⊗ |ψ spins sq .
The last term describes the underlying spins in squeezed space. For a single hole in the t − J model we choose it to be the ground state of a 2D Heisenberg AFM on a square lattice. The strong coupling ansatz is a direct generalization of Eq. (22) for the case of the t−J z model, where |ψ spins sq is given by a classical Néel state. As a first step, we confirm that spinons and holons are confined in the t − J model. We use the strong coupling wavefunction Eq. (63) as a variational ansatz and note that geometric strings in general correspond to highly excited states of the surrounding spin system. Two spins which are nearest neighbors (NN) in squeezed space can become next nearest neighbors (NNN) in real space after the holon has moved through the system, and vice-versa.
To estimate the energy E Σ of the resulting state, we consider straight holon trajectories for simplicity, for which E Σ ∝ is proportional to the length of the geometric string. The string tension obtained within this LST approximation is given by Here C 1 and C 2 denote the NN and NNN spin-spin correlation functions in squeezed space. In the case of the t − J z model C 1 = −C 2 = −S 2 , and by setting J = J z we obtain the string tension 4J z S 2 which we previously used for the LST, see Eq. (20). Because the ground state expectation values C 1,2 can be determined numerically for the 2D Heisenberg AFM, we can calculate the LST string tension for the t − J model: dE Σ /d ≈ 1.1J. For S = 1/2 the string tension for the t − J model in units of J is about 10% larger than for the t − J z model in units of J z . This can be understood by noting that the isotropic term JŜ i ·Ŝ j in the t − J model includes couplings in x and y-direction in spin space, in addition to the Ising coupling relevant in a weakly doped t − J model. We compare results by exact diagonalization (ED) of one hole in a 4 × 4 system with periodic boundary conditions to predictions from LST based on the strong-coupling wavefunction Eq. (63). The correlators are evaluated only for states where both sites 0 and d are occupied by spins, as indicated by the notation ... occ. In the ED calculations we set S z tot = −1/2 and simulated the sector with total conserved momentum of k = (π/2, π/2).
to the t − J z model. As a result, we expect that spinons and holons are confined in the t − J model with typical string lengths of comparable size in both models.
Next we check our assumption in Eq. (63) that the motion of the hole only has a negligible effect on the spins in squeezed space. To this end we compare the typical time scales for spin-exchange processes, τ fluc = 1/J, with the typical time τ h = 0 /t it takes the holon to cover a distance given by the average string length 0 . If τ h < τ fluc , quantum fluctuations do not have sufficient time to adapt to the new configuration of the geometric string and we expect that spin correlations in squeezed space are only weakly modified by the presence of the hole. Indeed, in the entire strong coupling regime t J we find that τ h τ fluc . When t/J = 3, as in the case relevant to high-temperature superconductors [3], we obtain τ h /τ fluc ≈ 0.4. This suggests that the strong coupling wavefunction (63) provides a good description of magnetic polarons in the t − J model.
As an application of the strong coupling meson theory, we use it to calculate local spin-spin correlations in the presence of a hole. They are directly accessible using quantum gas microscopes [60] and can be used to test our predictions experimentally. In Fig. 24 we compare exact numerical simulations of the t − J model for a single hole with results from LST. Remarkably, the strong coupling meson theory can explain quantitatively the dependence of the spin correlations on the ratio J/t, caused entirely by doping. It merely combines the holon motion on the Bethe lattice with a set of spin correlators in a pure spin system without doping: Essentially NN spin correlations in real space acquire contributions from NNN, or even more distant, correlators in squeezed space when the holon motion is taken into account. This leads to the observed decrease (increase) of NNN (NN) correlations.
We find good agreement of LST with ED simulations in Fig. 24 for NN correlators in the regime J 0.1t. For correlators between sites which are further apart, there is still perfect qualitative agreement, although some quantitative differences are visible. Only at very small values of J/t 0.1 we find strong deviations of LST from ED results, which can be associated with a transition to a state with strong ferromagnetic NN correlations reminiscent of the Nagaoka effect [14]. We expect that this transition is strongly influenced by the finite system size.
In a forthcoming work, we show how spinon dynamics can be included and more properties of magnetic polarons in the t − J model can be understood from a strong-coupling meson description. The most dramatic difference as compared to the t − J z model is that spinon motion is possible not only via Trugman loops, but also by direct spin-exchange processes. This readily explains why some properties of holes in the t − J and t − J z models -including for example the linear dependence of their energy on t 1/3 J 2/3 [28]-are identical, whereas the dispersion relations are entirely different.
Finally, we note that the strong coupling ansatz in Eq. (63) can be generalized to finite temperatures by using a product of density matrices, as described for one dimension in Ref. [44]. In this case the correlators C 1 and C 2 determining the string tension in Eq. (64) are calculated at finite temperature. The resulting weaker correlations lead to a reduced string tension, and thus to longer strings. As shown in Ref. [82], even at infinite temperature and for J t, the Hilbert space H Σ defined by geometric strings provides a qualitatively accurate description of the single-hole t − J model.

VII. SUMMARY AND OUTLOOK
When holes are doped into an anti-ferromagnet at low concentration, they form quasiparticles called magnetic polarons. In this paper we have introduced a strong coupling parton theory for magnetic polarons in the t − J z model. Starting from first principles, we introduced a parton construction where the magnetic polaron is described as a bound state of a spinon and a holon. The spinon carries its fractional spin S = 1/2 and the holon its charge Q = −1 quantum numbers. The two partons are connected by a string of displaced spins. This description combines earlier theoretical approaches to magnetic polarons, see in particular Refs. [15,18,19,49]. Our formalism can be generalized to Hamiltonians including quantum fluctuations, e.g. the isotropic t − J model.
Ultracold atoms provide a new opportunity to investigate individual magnetic polarons experimentally, on a microscopic level. In this paper, we have discussed several realistic setups to measure their properties and study their dynamics far from equilibrium. We start from a sin-gle hole which can be pinned by a tightly focused laser beam. To investigate magnetic polarons in equilibrium using quantum gas microscopes, we proposed an adiabatic preparation scheme where the pinning potential is adiabatically released. If the laser potential localizing the hole is suddenly switched off, on the other hand, interesting far-from equilibrium dynamics can be studied.
Using the capabilities of quantum gas microscopy, direct observations of the constituent partons forming the magnetic polaron are possible. We analyzed this problem and showed that the full distribution function p of string lengths connecting spinons and holons can be accurately measured.
As a hallmark of the parton theory, we point out the existence of rotational as well as vibrational excited states of magnetic polarons. Similar to mesons in highenergy physics, which are understood as bound states of two confined quarks, these excited states give rise to resonances with a well-defined energy. Creating the first vibrational excitation costs an energy ∼ t 1/3 J 2/3 z , whereas the first rotational excitation only costs an energy ∼ J z . While vibrational excitations of magnetic polarons have been discussed previously in the context of angle-resolved photoemission spectroscopy (ARPES), rotational resonances are invisible in ARPES spectra due to selection rules. To observe rotational states of magnetic polarons we propose to use Bloch oscillation spectroscopy: In this method a weak force is applied to the hole, which drives transitions to the first rotationally excited state at energy ∆E. This causes oscillations of the hole density distribution at frequency ∆E, which can be measured in experiments with ultracold atoms.
We have applied the strong coupling parton theory to study dynamical properties of holes inside an antiferromagnet. Our calculations are based on a strongcoupling product ansatz for describing the spinon and holon parts of the wavefunctions. In contrast to models with transverse couplings between the spins, the dynamics of spinons in the t − J z model is generated entirely by Trugman loop trajectories of the holon which restore the surrounding Néel order. To calculate the resulting hopping strength of the spinon we developed a tight-binding formalism which is valid for arbitrary values of J z /t. Our predictions are in excellent agreement with exact calculations for the largest system sizes which are numerically feasible [4,109].
We have considered a far-from equilibrium situation, where a single hole is created in a Néel state by removing the central spin. As a benchmark, we have performed time-dependent quantum Monte Carlo calculations for the t − J z model at short-to-intermediate times. We compared them to our strong coupling parton theory and obtained very good quantitative agreement for the accessible times. At short times, we predict universal holon expansion. At intermediate times ∼ 1/J z , we observe a saturation of the string length connecting the spinon and the holon. At much longer times, the ballistic expansion of the spinon becomes dominant.
The strong coupling parton description of magnetic polarons can be extended to study a larger class of problems. In a first approach, we will include magnon corrections describing quantum fluctuations around the classical Néel state. This allows us to study, for example, the dispersion relation of magnetic polarons in the usual t − J model with Heisenberg interactions between the spins. In addition, spin-hole correlation functions can be calculated, which are directly accessible in experiments with quantum gas microscopes [60]. The dynamics of the magnon cloud forming around the spinon-holon bound state can also be studied. In a second approach, we generalize the concept of squeezed space known from one dimension [116] and calculate properties of the magnetic polaron in the t − J model using a variational wavefunction. In this paper we have presented first results of this approach, which can be tested using current experiments with fermionic quantum gas microscopes.
The couplings to phonons present in real solids can also be included in our theory using the Landau-Pekar wavefunction [101]. Far-from equilibrium dynamics can then be studied, where the spinon becomes correlated with the magnons and the phonons. This can be of particular interest to describe recent experiments where cuprates have been driven far-from equilibrium by a short laser pulse [117]. Néel state by removing the central spin, as in Sec. V D. For short times we obtained excellent agreement between LST and NLST, which we benchmarked by comparing to time-dependent Monte Carlo calculations in Fig. 3.
In Fig. 25 (a) we show the rms radius of the hole density distribution for longer times. Note that LST calculations include the slow spinon dynamics, whereas the NLST results are performed without taking into account the motion of the spinon. For J z = t we find good agreement between NLST and LST. At long times the ballistic spinon expansion can be observed in the LST result, which is not included in NLST. The length of the string in the long-time limit is the same in LST and NLST in this case.
For J z /t = 0.5 we observe quantitative deviations between our results from NLST and LST in Fig. 25 (a). In particular, the LST predicts stronger oscillations at long times. In the NLST many quantum states with slightly different energies are included. This is expected to cause dephasing, which explains the observed suppression of coherent string dynamics.
For J z /t = 0.2 we observe large deviations between LST and NLST. The rms radius predicted by NLST suddenly saturates at intermediate times, while it continues to grow in the LST. This is an artifact of the NLST result, where the maximum string length used in the calculations was max = 10. For J z /t = 0.2 longer strings need to be included. Our LST calculations use strings up to a length max = 30, which is sufficient to obtain fully converged results. Note that LST can easily go to even longer values of the string length.
In Fig. 25 (b) we also compare predictions for the probability of finding the spinon and the holon on the same site. As before we obtain excellent agreement for short times. For longer times, the LST predicts more oscillations, but the qualitative behavior is the same in both theories. In this appendix we introduce the formalism for the generalized 1/S expansion and explain how it can be used to include quantum fluctuations in the strong coupling parton description of magnetic polarons.

Shortcoming of the conventional 1/S expansion
In the conventional 1/S expansion of magnetic polaron problems [17,19,48], S only enters in the corresponding Schwinger-boson constraint, σb † jσb jσ +ĥ † jĥ j = 2S. In this case the holon hopping is described by Eq. (7) with the operatorF ij (S) ≡F ij (1/2) from Eq. (8), for all values S. Note that the conventional Schwinger-boson constraint is respected by the resulting operatorĤ t .
Let us consider the effect of such holon motion on the local Néel order parameterΩ j , see Eq. (9). Within the conventional extension of the t − J z model to large values of S, the motion of the holon from site i to j is accompanied by changes of the spins S z i and S z j by ±1/2, see Fig. 4 (b). As a result, the sign of the local Néel order parameter Ω j cannot change when S 1/2 is large, unless the holon performs multiple loops.
We conclude that the conventional 1/S expansion cannot capture correctly the destruction of the Néel order parameter by the holon motion, at least when S 1/2 is large. This effect is particularly pronounced for a hole inside a one dimensional spin chain with S = 1/2. In this case a single hole introduces a domain wall into the system, where the local Néel order parameter changes signa direct manifestation of spin-charge separation [43,44]. The conventional 1/S expansion cannot capture this effect because it only allows for small local changes of the underlying Néel order parameter.
The shortcoming of the conventional 1/S expansion explained above is also reflected by the fact that the new Schwinger-boson constraint in Eq. (4) is not satisfied for S > 1/2. This means that the holon can move without fully distorting the underlying Néel order in this case. When the generalized constraint from Eq. (4) is enforced, the local Néel order can be destroyed by a single hole hopping event. This is an important feature of the t − J z model at S = 1/2, and it is fully captured -even for large values of S -by the generalized model which we introduce now.

Generalized holon hopping for large S
In order to ensure that the generalized Schwingerboson constraint Eq. (4) is satisfied byĤ t from Eq. (7), the operatorF ij (S) in Eq. (7) has to depend explicitly on the value of S. We will refer to the new termF ij (S) as the generalized holon-hopping operator.
To constructF ij (S) we recall the physical origin of nearest-neighbor hole hopping in the original model Eq. (3): The fermion initially occupying site j moves to site i without changing its spin state, see Fig. 4 (c). For S = 1/2 this process is accurately described within the Schwinger-boson language by the operatorF ij (1/2) from Eq. (8), which changes the z-component of the spins on sites i and j by ±1/2 respectively.
For larger S > 1/2 the entire spin state is transferred from site j to i by the operatorF ij (S). Because expressing this process in terms of Schwinger-bosons leads to a highly non-linear and cumbersome expression, we define it using the eigenstates |n l ↑ , n l ↓ ofb † l,↑b l,↑ (eigenvalue n l ↑ ) andb † l,↓b l,↓ (eigenvalue n l ↓ ) on site l = i, j now: |0 j , 0 j , n i ↑ , n i ↓ n j ↑ , n j ↓ , 0 i , 0 i |. (C1) Note thatĥ † jĥ iFij (S) respects the constraint in Eq. (4), because the operatorF ij (S) only acts on the subspace with zero Schwinger-bosons on site i and 2S Schwingerbosons on site j. expansion for describing next-nearest neighbor hopping processes. I.e. we treat the string as unbroken, but create magnon excitationsâ i when the holon is hopping to a next-nearest neighbor site. This reflects the distortion of the physical spin state: On the one hand the spin configuration for S = 1/2 is correctly described within the Holstein-Primakoff approximation, see Eq. (C2). On the other hand, the energy gain from the new spin configuration is taken into account. the wavefunction would need to be localized in one particular direction of the Bethe lattice, along the Trugman loop trajectory. This is energetically unfavorable because it only leads to a potential energy gain ∼ J z but costs localization energy ∼ t: Localizing the particle in one particular direction only leads to a zero-point energy −2t.