Exact dynamics of two holes in two-leg antiferromagnetic ladders

We study the motion of holes in a mixed-dimensional setup of an antiferromagnetic ladder, featuring nearest neighbor hopping $t$ along the ladders and Ising-type spin interactions along, $J_\parallel$, and across, $J_\perp$, the ladder. We determine exact solutions for the low-energy one- and two-hole eigenstates. The presence of the trans-leg spin coupling, $J_\perp$, leads to a linear confining potential between the holes. As a result, holes on separate legs feature a super-linear binding energy scaling as $(J_\perp / t)^{2/3}$ in the strongly correlated regime of $J_\perp,J_\parallel \leq t$. This behavior is linked to an emergent length scale $\lambda \propto (t/J_\perp)^{1/3}$, stemming from the linear confining potential, and which describes how the size of the two-hole molecular state diverges for $J_\perp,J_\parallel \ll t$. On the contrary, holes on the same leg unbind at sufficiently low spin couplings. This is a consequence of the altered short-range boundary condition for holes on the same leg, yielding an effective Pauli repulsion between them, limiting their kinetic energy and making binding unfavorable. Finally, we determine the exact nonequilibrium quench dynamics following the sudden immersion of initially localized nearest neigbhor holes. The dynamics is characterized by a crossover from an initial ballistic quantum walk to an aperiodic oscillatory motion around a finite average distance between the holes due to the confining potential between them. In the strongly correlated regime of low spin couplings, $J_\perp, J_\parallel \leq t$, we find this asymptotic distance to diverge as $t / J_\perp$, showing a much stronger scaling than the eigenstates. The predicted results should be amenable to state-of-the-art quantum simulation experiments using currently implemented experimental techniques.


I. INTRODUCTION
Quantum simulation experiments have matured to the level, where they push our understanding of many-body quantum dynamics and inspire new approximate theoretical tools [1][2][3][4][5] that allow us to explore the complex spatial structures arising in e.g. Fermi-Hubbard systems [6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22]. A major driver for this line of research is to better understand the microscopic origins of high-temperature superconductivity [23], which basic phenomenology may be explained by the interaction and ensuing pairing of dopants in Fermi-Hubbard systems [24][25][26]. Recent experiments [27] have for the first time successfully demonstrated that cold-atom simulators can achieve and probe the formation of such pairs in a particular kind of spin ladders. Whereas these experiments were still limited to rather small system sizes, they have shown a great promise in how we can understand these mechanisms from the bottom-up perspective, and the approximate theoretical description of this situation [28] suggests that the system supports a strong binding of the dopants, in contrast to the usual scenario in twodimensional square lattice geometries [29][30][31][32][33][34][35]. Importantly, spin ladders also arise in compound materials supporting unconventional superconductivity [36,37]. While these compounds are mainly probed in scattering experiments [38][39][40][41][42][43], cold-atom simulators give direct access to spatial correlations and nonequilibrium dynamics [10,17,20].
Inspired by this development, we analyze a situation, in which we may gain exact results for the binding and nonequilibrium dynamics of dopants in a mixed-dimensional Fermi-Hubbard system [Figs. 1(a)-1(b)]. The main theoretical difficulty in previous studies [28] has been to fully describe isotropic spin couplings, coming with the complication of an underlying order-disorder phase transition as the trans-leg spin coupling is increased [22,[44][45][46][47][48]. However, by restricting the spin interactions to the Ising type, the underlying spin lattice always supports a perfectly Neél-ordered ground state. Based FIG. 1. Mixed-dimensional t-J model featuring spin-1/2 particles on a two-leg ladder geometry with two holes on separate legs (a) or on the same leg (b). The spins can hop to nearest neighbor vacant sites along the ladder with amplitude −t, and have nearest neighbor Ising interactions J ∥ , J ⊥ along and across the ladder, respectively. (c) Average distance between two holes versus time τ , which initially sit at nearest neighbor sites. At short times, the holes blow apart ballistically as described by a quantum walk [red line], after which they oscillate around a well-defined long-time average. This is shown for J ⊥ /t = 0. on this simplification, we find that we can describe the lowenergy single-and two-hole eigenstates exactly in this case, whether they be on the same or separate legs [Figs. 1(a)-1(b)]. Furthermore, using the precise insights into the twohole eigenstates, we calculate the exact quench dynamics following the sudden creation of two holes as nearest neighbors. Here, Figs. 1(c)-1(d) show the result of holes on separate legs. We find that the dynamics can be divided into two characteristic regimes. First, the holes perform independent quantum walks, meaning that they blow apart ballistically. Second, as the holes diverge from each other, a confining string of overturned spins forms between them. Eventually, the holes are slowed down by this confinement and aperiodic oscillations in the strings, and thereby in the inter-hole distance, take place around a well-defined long-time average.
A major challenge in previous experiments with doped Fermi-Hubbard systems [14,21] has been to reach the strongly correlated regime, which is interesting both from the perspective of the physics of the cuprate materials supporting high-temperature superconductivity [23], and for understanding many-body phenomena outside the realm of perturbation theory. We note that this system is a natural experimental candidate for that, because the effective coupling strength between the holes is 4t/J ⊥ . Consequently, the crossover timescale from the quantum walk to the string oscillation behavior in Fig. 1(d) changes from a perturbative (t/J ⊥ ) 2 scaling to a strongly correlated scaling of t/J ⊥ already for J ⊥ ≲ 4t. Importantly, the crossover time is still quite moderate in terms of hopping times 1/t, and remains below 3/t for J ⊥ > t, which should make it possible to experimentally observe the departure from the quantum walk. While the mixed-dimensional property of this model has already been achieved experimentally [27], the ability to tune the spin interactions to the Ising limit can be facilitated by Rydberg-dressed atoms [49][50][51][52][53][54], polar molecules [55], and trapped ions [56]. This setup is, therefore, within reach for modern quantum simulation experiments.
The paper is organized as follows. In Sec. II, we set up the mixed-dimensional t-J model. In Sec. III, we determine the low-energy single-and two-hole eigenstates. In Sec. IV, we study the nonequilibrium quench dynamics of two holes, before we conclude in Sec. V. Throughout the paper, we set the reduced Planck constant, ℏ, and the lattice spacing to 1.

II. MODEL
We consider a system of spin-1/2 particles placed along a two-leg ladder, described by a mixed-dimensional t-J model with HamiltonianĤ =Ĥ t +Ĥ J [ Fig. 1(a)]. The spins σ =↑ , ↓ can hop to nearest neighbors along the legs µ = 1, 2, under the constraint that there is at most a single spin on each site. This is enforced by the modified particle operator c † j,µ,σ =ĉ † j,µ,σ (1 −n µ,j ), withn µ,j = σĉ † j,µ,σĉ j,µ,σ the local density operator. The nearest neighbor antiferromagnetic spin-spin coupling is assumed to be fully polarized in the z- with J ⊥ , J ∥ > 0. Such mixed-dimensional models [57,58] have recently been proposed to yield strong binding of holes through an emergent confining string potential of overturned spins [28], and was recently implemented successfully in the case of fully symmetric spin couplings [27]. The polarized Ising-type interaction explored here, enables us to derive exact results for low-energy single-and two-hole eigenstates, as well as the full nonequilibrium dynamics of two initially localized holes.

III. LOW-ENERGY EIGENSTATES
In the absence of holes, the polarized AFM coupling in Eq. (2) results in a perfect Neél ordered state, |AF⟩, for any values of J ∥ , J ⊥ > 0. For periodic boundary conditions of N spins along each of the two legs, this results in the ground state energy owing to a nearest neighbor spin bond energy of J ∥ /2 (J ⊥ /2) along (across) the ladder. This should be contrasted to the case of isotropic spin couplings, in which case there is a quantum phase transition between Neél order along the ladder and spin singlet formation along the rungs as J ⊥ /J ∥ is increased [22,[44][45][46][47][48]. Utilizing this simplification, we can find the single-hole and two-hole ground states. To have a more efficient description, we employ a Holstein-Primakoff transformation and describe the system in terms of holes,ĥ, and bosonic spin excitationsŝ. The latter operators are defined with respect to the antiferromagnetic ground stateŝ |AF⟩ = 0. The hopping Hamiltonian then reads [59] Here, F (ĥ,ŝ) = 1 −ĥ †ĥ −ŝ †ŝ ensures that there is at most a single spin excitation and a single hole on each site. The spin-coupling Hamiltonian likewise becomeŝ We emphasize that the spins can both be fermions and hardcore bosons. In fact, if the holes sit on separate legs, they are distinguishable by which leg they move in. If they move along the same leg, they are equivalently distinguishable by which one is to the left and which one is to the right. As a result, the statistics never come into play in what follows, only the hard-core constraint and the one-dimensional nature of the motion. The results, therefore, apply equally well to fermionic and hard-core bosonic spins, as one might expect from the general duality of fermions and imprenatable bosons in one dimension [60].

A. Single-hole eigenstates
Central to the analysis of a single hole is the insight that a single hole doped into the two-leg antiferromagnetic Ising ladder is localized. Due to inversion symmetry, the low-energy eigenstates may then be written as (assuming that the hole resides in leg 1) This describes that for hole positions d sites away from the central site, with total amplitude C (d) , a resulting string of overturned spins of length d appears. Taking the energy of a single stationary hole, E 0 + J ∥ + J ⊥ /2 as reference, and utilizing the Schrödinger equation,Ĥ |Ψ 1 ⟩ = E 1 |Ψ 1 ⟩, we obtain the equations of motion The lower equation applies for d ≥ 2. The motion of the hole d sites away leaves behind a single frustrated spin bond in leg 1, as well as d frustrated spin bonds across the ladder. This results in the linear string potential confining the hole around to its origin. The obtained equations of motion are identical in form to the recently obtained exact results in general Bethe lattices [61]. Utilizing the same techniques for solving the equations of motion in Eq. (7), without loss of generality we seek for a recursive structure of the amplitudes for d ≥ 1. Inserting this into Eq. (7), we obtain the selfconsistency equation for a yet to be determined function f 1 . As Eq. (10) is independent of the distance d, f 1 (E) is as well. The self-consistency equation (10) can, finally, be used to find a closed-form expression of f 1 (E) in terms of Bessel functions of the first kind, with Ω(E) = −2E/J ⊥ , similar to the results in Refs. [61][62][63]. Inserting f 1 (E 1 ) = f 1 (E − V 1 (2)) in the equation for d = 1 in Eq. (7), hereby, yields C (1) = √ 2t · f 1 (E − V 1 (1))C (0) . Inserting this into the equation for d = 0 in Eq. (7), then finally results in the equation for the single-hole energy, hereby defining the single-hole self-energy Σ 1 (E). Equation (12) actually supports a discrete series of eigenstates similar to a single hole in a Bethe lattice [61]. Here, however, our main focus is on the ground state as this is important to decipher whether two holes will bind or not. The recursive structure of the amplitudes along with the result in Eq. (11), thus, allows us to construct the full many-body eigenstate with Here, by normalizing the wave function, 1 = ⟨Ψ 1 |Ψ 1 ⟩.
At strong coupling, J ⊥ , J ∥ ≪ t, the hole spreads out more and more, resulting in a continuum limit. This yields the asymptotic single-hole energy with −a (0) ≃ −1.02 the first zero of the derivative of the Airy function, Ai ′ (x). In Fig. 2(b), we plot the single hole energy as a function J ⊥ for a few indicated values of J ∥ . We see good agreement between Eq. (14) and the full solution of Eq. (12) for J ∥ = 0.01t and J ⊥ ≪ t.

B. Two-hole eigenstates
We now focus on the low-energy two-hole eigenstates. We both consider holes moving on separate legs [ Fig. 3], as well as holes moving along the same leg [ Fig. 4]. For holes traveling on separate legs, the breaking of spin bonds within a leg can be completely avoided by starting from a configuration of spins in which the spins to the right of the holes is moved by one lattice site to the right. Hence, if the holes move alongside each other, the perfect Neél order is retained and no spin bonds are broken. The appropriate two-hole eigenstates are, therefore, In blue is shown ±4t cos(k/2) for reference. The spectrum, here, depends on both the trans-(J ⊥ ) and intra-leg (J ∥ ) spin couplings. For J ⊥ → 0 (right), a quasiparticle band appears below the continuum, when J ∥ ≥ 4t cos(k/2). For J ∥ = 2t, this corresponds to k ≥ 2π/3. delocalized along the ladder. We, thus, define the states for a linear distance d between the two holes. Here, we assume N sites in each leg and periodic boundary conditions. In this manner, d > 0 (d < 0) indicates that the hole in leg 2 has moved |d| sites to the right (left). The appearance of the string operator, l>jŝ † l,1 m>j+dŝ † m,2 , is due to the shift of the underlying AFM order by one lattice site to the right at j and j + d. These states have crystal momentum k ∈ (−π, π], as translating the holes and spin stringŝ h † j,1ĥ † j+d,2 l>jŝ † l,1 m>j+dŝ † m,2 |AF⟩ by one lattice site to the right results in an additional phase of e −ik . As no spin frustration within a leg occurs for this configuration of holes, the resulting low-energy eigenstates are independent of the intraleg spin coupling J ∥ .
For holes moving along the same leg, the most favorable configuration of the spins is now obtained by taking out two adjacent spins. Once again, if the holes move alongside each other, no spin bonds are broken. The delocalized states for a distance d between the holes in this case, therefore, becomes in which we see that a string of overturned spins forms between the two holes. Since two holes cannot sit on top of each other, let alone pass through one another, the distance is now always greater than one lattice site, d ≥ 1. The full two-hole eigenstate can, hereby, be written as where s =⊥, ∥ denotes whether the holes move on separate legs (⊥) or along the same leg (∥). Additionally, the band index n = 0, 1, 2, . . . specifies that a discrete series of twohole bands emerge, which will become essential for describing the nonequilibrium dynamics in Sec. IV. The normalization condition is 1 = ⟨Ψ Crucially, the hopping Hamiltonian only couples the states within Eqs. (15) and (16) in a well-defined hierachy. In particular, it couples holes moving on separate legs as follows: |Ψ 2s (k)⟩ for holes on separate legs (s =⊥) and on on the same leg (s =∥) leads to the equations of motion where the lower line in Eq. (19)  (k)⟩ via two pathways. For holes on separate legs, this corresponds to the hole in leg 2 hopping to the right with amplitude te −ik/2 , and the hole in leg 1 hopping to the left with amplitude te +ik/2 . For holes on the same leg, it similarly corresponds to the hole to the right to hop further to the right with amplitude te −ik/2 , and the hole to the left to hop further to the left with amplitude te +ik/2 . In any case, the associated quantum interference of these pathways leads to the total coupling of t(e −ik/2 + e +ik/2 ) = 2t cos(k/2), as was also recognized previously [28]. Here, we define the two-hole linear string potentials using the energy of two separate stationary holes, E 0 + J ⊥ + 2J ∥ , as reference. We emphasize that for holes traveling on separate legs, the string potential does not contain the spin coupling along the ladder J ∥ , because no intra-leg spin bond is broken in this case. For holes moving along the same leg, the intra-leg spin coupling J ∥ only appears at d = 1, as the two holes here share a frustrated intra-leg spin bond. Similar to the single-hole case, we propose the recursion relations Here, the upper line applies for both configurations of holes for d ≥ 0 (s =⊥) and d ≥ 1 (s =∥). The lower line is solely for holes on separate legs in the case of d ≤ 0. Inserting this into Eqs. (18) and (19) results in the self-consistency equation applying both to holes on separate legs and on the same leg. Analogous to the single hole case, we set . As Eq. (22) is independent of d and the hole configuration s =∥, ⊥, so is f 2 . Note that for holes on separate legs, this also means that f (d) = f (−d) , and that C (n) (k, −d) = C (n) (k, d) as one might expect from inversion symmetry of the system. Equation (22) has the exact same structure as Eq. (10) for the equivalent function f 1 in the single-hole case. As a result, we may simply replace t → 2t cos(k/2) to once again obtain a closed-form expres-sion in terms of Bessel functions of the first kind As for the single-hole case, we can use the recursion relations in Eq. (21) along with Eq. (23) to explicitly write the coefficient of the many-body wave function. For holes on separate legs, we get for |d| ≥ 1, using Ω(E − V 2⊥ (d)) = Ω(E) + |d| − 1. For holes on the same leg, we similarly get for d ≥ 2. Analogous to the single-hole case, Z (n) 2s (k) ] −1 is the residue of the two-hole Green's function for η = 0 + , Σ 2⊥ (k, ω + iη) = 8t 2 cos 2 (k/2) · f 2 (k, ω + iη), and Σ 2∥ (k, ω + iη) = 4t 2 cos 2 (k/2) · f 2 (k, ω + iη − J ⊥ /2). The poles of G 2s , hereby, determine the spectra for states of the forms in Eqs. (15) and (16). These are all states that have a nonzero overlap with finding holes at adjacent sites with no frustrated spin bonds. Importantly, this subfamily of states contains the two-body states with the lowest energy. The spectral functions, A 2s (k, ω) = −2ImG 2s (k, ω) are shown in Figs. 3 and 4 for a few indicated values of the spin couplings. From here, the discrete bands, E 2s (k), with n = 0, 1, 2, . . . are now apparent. In the limit of J ⊥ /t → 0 + , a continuum of states form between ±4t cos(k/2). Below this continuum of states, holes traveling on the same leg [ Fig. 4(b)], feature a well-defined quasiparticle state if J ∥ > 4t cos(k/2), in which case Eq. (24) may be solved to yield ending up at −J ∥ /2 at the Brillouin zone boundary, k = π.
For J ∥ > 4t [bottom right in Fig. 4(b)], this state appears for any k and a full quasiparticle band remains even for J ⊥ → 0.
For 0 < J ∥ < 4t, on the other hand, a quasiparticle state appears only for crystal momenta close enough to the boundary of the Brillouin zone [top right in Fig. 4(b)]. We note that at k = π for general J ⊥ > 0, there seems to be an equal spacing of the bands. In fact, inspecting Eqs. (18) and (19) we see that the two hopping pathways destructively interfere here, giving a vanishing total hopping amplitude, 2t cos(π/2) = 0. The states are, therefore, completely immobile and their energies are, consequently, determined by the string potentials where n ≥ 0 (n ≥ 1) in the upper (lower) line. This gives a spacing of J ⊥ /2 at k = π. We note that the overall structure of the spectra in Fig. 3 are similar to the isotropic spin coupling case in the regime of J ⊥ ≫ J ∥ , where the underlying spin lattice resides in a disordered regime of spin-singlets on each rung [28]. Finally, for the lowest energy two-hole state, k = 0 and n = 0, we find that the energies at strong coupling, J ⊥ , J ∥ ≪ t, behaves as E In Fig. 5(a) binding mechanism for holes on separate legs. For holes on the same leg, however, since the prefactor in front of this term is negative, 2a (0) − 2 1/3 a (1) ≃ −0.90 < 0, two holes on the same leg actually energetically prefer to unbind. Similar to recent cold-atom experiments with isotropic spin couplings [27], this difference can be understood from a Pauli repulsion effect. In fact, the hard-core constraint means that the boundary condition at d = 0 is altered from being soft for holes on separate legs to exactly zero for holes on the same leg. This results in the different prefactors of a (0) ≃ 1.02 and a (1) ≃ 2.34 in the two cases, which will become apparent when we investigate the spatial distribution of the holes below. We may note, however, that already for moderate values of the intraleg spin-coupling J ∥ , this unbinding is overcome and eventually reaches J ∥ /2 for J ∥ ≫ t. In fact, in the extreme limit of J ⊥ /t → 0 + Eq. (28) in combination with E 2∥ (0) = −4t for J ∥ /t < 4t, results in the positive binding energy shown with a black line in Fig. 5(b). Here, we use Eq. (28) for J ∥ < 4t and E 2∥ (0) = −4t for J ∥ ≥ 4t, as well as E 1 = J ∥ /2 − 2t 1 + (J ∥ /4t) 2 by solving Eqs. (10) and (12) for J ⊥ → 0 + . Hence, in this limit the binding energy  interpolates between two linear behaviors in the intra-leg spin coupling, from an initial J ∥ to J ∥ /2 behavior. This illustrates that two holes on the same leg bind unless both the trans-and intra-leg spin couplings are small. Furthermore, we stress once again that these results, including the unbinding mechanism for holes on the same leg, ensues regardless of the statistics of the spins and only depends on the hard-core constraint, as one should also expect for in a system with one-dimensional motion [60].
In this manner, we have given a detailed account of the lowenergy behavior for both intra-and trans-leg configurations. Holes on separate legs always bind with a super-linear scaling of t · (J ⊥ /t) 2/3 for J ⊥ , J ∥ ≪ t. For holes on the same leg, however, the hard-core constraint results in an energy cost proportional to t · (J ⊥ /t) 2/3 for low J ⊥ , J ∥ and leads to unbinding in this regime. However, for higher values of either spin coupling the holes will, once again, bind.
Whereas a determination of the two-hole binding energy is direct proof of their ability to bind, it is simultaneously notoriously difficult to measure in modern quantum simulation experiments with ultracold atoms in synthetic lattices, such as optical lattices and Rydberg arrays. The simple reason is that the required spectroscopy entails single atom detection in e.g. time of flight or rather advanced band-mapping techniques [16,64]. On the other hand, the combination of the lattice experiments and the development of quantum gas microscopy has enabled the direct and precise measurement of spatial correlations, and has successfully been employed to measure antiferromagnetic correlations in Fermi-Hubbard systems [7,9], as well as characterizing the spatial properties [14] and formation dynamics of magnetic polarons [21] in such systems. For two holes, the two-point hole-hole correlators provide such a spatial probe of their binding, as was also recently used in experiments [27]. In Eq. (32), the average is taken for the states |Ψ   s (k, d)| 2 /N , whereas the uniform spreading of the holes means that ⟨ĥ † µ,jĥ µ,j ⟩ k = 1/N , for both legs µ = 1, 2. In Figs. 6(a)-6(b), we plot these correlators as a function of d for several values of J ⊥ . For lower values of J ⊥ /t, the holes separate more and more from each other as one expects for a higher mobility. We note that already for J ⊥ = 3t, the probability of finding the holes as nearest neighbors has dropped to around 50%. For J ⊥ /t ≪ 1 -and J ∥ /t ≪ 1 in the intra-leg case -the relative wave functions of the holes, C with the effective inverse length scale λ(k) = [J ⊥ /(4t cos(k/2))] 1/3 , and the normalization constants A j . We compare to these continuum results and see remarkably good agreement away from d = 0 even for relatively large values of J ⊥ . Additionally, we show the average distance between the holes ⟨|d|⟩ sk = d |d|g (2) s (k, d)/N as a function of J ⊥ /t in Fig. 6(c) for the ground state at k = 0. This reveals the strong-correlation scaling for J ⊥ /t ≪ 1. Figure 6(c) shows that this asymptotic form is already accurate for J ⊥ /t ≤ 1. We attribute this to the fact that the effective interaction strength for two holes is 4t/J ⊥ , rather than just t/J ⊥ . For weak correlations, we similarly get ⟨d⟩ ∥k − 1 = ⟨|d|⟩ ⊥k /2 = 1/λ 6 (k) ∝ (t/J ⊥ ) 2 , becoming accurate for J ⊥ /t ≥ 10 in Fig. 6(c). Importantly, we emphasize that for holes on the same leg, the hard-core constraint g (2) ∥ (k, d = 0) = 0 results in a different boundary condition in the continuum limit. This change in the boundary condition alters the relative wave function from being on the form of Ai(λ(k)|d| − a (0) ) to Ai(λ(k)d − a (1) ), and changes the prefactor of the (J ⊥ /t) 2/3 term in the two-hole energy from a (0) ≃ 1.02 to a (1) ≃ 2.34. This also results in a significant qualitative change in the relative spatial distribution of the holes. For holes on separate legs, the holes are always most likely to be found as nearest neigbors, whereas this is not true at all for holes on the same leg. This leads to more distant holes in the intra-leg configuration and explains the extra factor of a (1) ≃ 2.34 in ⟨d⟩ ∥k . In Fig. 6, we focus on the ground state behavior, i.e. k = 0. We note, however, that as the momentum approaches the edge of the Brillouin zone, the correlator compresses more and more and eventually the holes only sit next to each other: g (2) ⊥ (k = π, 0) = N and g (2) ∥ (k = π, 1) = N . Equation (33) reveals that for holes on separate legs, the correlator at d = 0 scales with the inverse length scale λ(k) ∝ (J ⊥ /t) 1/3 . Since the binding energy scales with t(J ⊥ /t) 2/3 , we get that E b⊥ /J ⊥ ∝ 1/g (2) (0, 0) at strong coupling. More precisely, ⊥ (0, 0)/N with c ⊥ = 2 −4/3 (1 − 2 −2/3 ) for J ⊥ , J ∥ ≪ t. This is very valuable for quantum simulation experiments, as it provides an indirect probe of the binding energy. In fact, in Ref. [27] an approximate relation at finite temperatures between the binding energy and the two-point correlator was used in this regard. For the configuration with two holes on the same leg, Eq. (33) similarly gives g (2) ∥ (k, 1) ∝ λ 3 (k). The asymptotic relation between the binding energy and the correlation function, therefore, now takes on the form ∥ (0, 0)/N ] 1/3 with c ∥ = 2 −1/3 (a (0) − 2 −2/3 a (1) ). This asymptotic relationship indicates that g ∥ (0, 0) must be much smaller to observe an impact on the binding energy. To explore these behaviors further, we plot the binding energy versus g (2) in Fig. 7. For holes on separate legs, this reveals a monotonic relation between the binding energy and the g (2) correlator for nearest neighbor holes for any value of J ∥ , which indeed enables experiments to infer a binding energy from a measured g (2) function. In the case of holes on the same leg, however, the applicability of this approach may, however, depend quite crucially on the value of J ∥ . In fact, for J ∥ ≪ t, we see that only at extremely low value of g (2) ∥ (0, 0) does the binding energy start to change significantly, which will naturally make it much harder to infer a binding energy from a measured g (2) function.

IV. NONEQUILIBRIUM DYNAMICS
In this section, we investigate the nonequilibrium dynamics of two initially localized holes. Such a quench experiment is a natural choice for quantum simulation experiments, and have recently been considered for the motion of a hole in a Fermi-Hubbard background [21], in which they were able to see the crossover dynamics from an initial free ballistic motion of the hole, signatures of string oscillations, and finally to the ballistic motion of magnetic polaron quasiparticles at long times [5,67]. In the current setup, we investigate the situation where the holes are localized and start out as nearest neighbor, described by the wave functions for the separate-legs (⊥)and same-leg (∥) configurations using τ as the variable for time to distinguish it from the hopping amplitude t. In the second line, as well as the last expression of the third line, we utilize that the initial state is the superposition of all the crystal momentum states in Eqs. (15) and (16)  full dynamics, we calculate the overlap of the initial states with the two-hole eigenstates in Eq. (17) Since the eigenstates are delocalized over the entire lattice, there is an overall factor of 1/  2∥ (k)⟩, respectively. See also Eqs. (25) and (26). We note that it is crucial to take into account that states in all the bands n have an overlap with the initial state and must be taking into account. The nonequilibrium wave functions are then for the separate-legs configuration, and for holes on the same leg. To describe the two-hole dynamics more concisely, we use Eqs. (39) and (40) to compute the probability of finding the holes at a distance d as a function of time describing the relative wave function versus time. Figures  8(a)-8(f) shows the dynamics of these probability distributions for several indicated distances, d. At short times, the holes initially blow apart ballistically as described by the quantum walks derived in Appendix B. For holes on the same leg, lower line in Eq. (42), the hard-core property of the holes constrains their motion and slightly alters it from the quantum walk of independent holes on separate legs. On longer timescales, the distribution of the holes is seen to oscillate around the time-averaged distributions  Fig. 6, while the scaling in the strong correlation limit is changed from (t/J ⊥ ) 1/3 for the eigenstates to t/J ⊥ for the dynamics.
which denotes the steady state approximately reached on long timescales. We note, however, that because the two-hole spectra in Figs. 3 and 4 consists of a discrete set of coherent peaks for any nonzero value of the trans-leg spin coupling J ⊥ , the motion will generally be highly aperiodic and never settle at its long-time average. As a result, the system does not fully equilibrate. It does still, however, give the characterize distribution of the holes at long times. To understand this further, in Figs. 8(g)-8(h), we compare it to the probability distribution for the two holes in the ground state. We observe that the behavior of the steady state is markedly different from the ground state. First and foremost, the state will dynamically extend over much larger length scales than its ground state counterpart. This is challenging for a cold-atom simulation experiment, and may hinder the observation of the long-time dynamics. However, we stress that already over a few hopping timescales 1/t, does the dynamics start to deviate from the quantum walk.
To investigate this more quantitatively using a simple experimental probe, we compute the average distances as a function of time. Two exemplary results are shown in Figs. 9(a)-9(b). For times τ < 2/t, holes on the same leg will de-part slightly slower than holes on separate legs, simply because there are more configurations available for holes on separate legs in the very first hop. After that, the hard-core constraint leads to faster divergent motion for holes on the same leg, but the motion remains a ballistic quantum walk. When the holes cross their long-time average, ⟨|d|⟩ ⊥ , ⟨d⟩ ∥ , the motion starts to deviate significantly from the initial ballistic behavior and instead crosses over to oscillations around ⟨|d|⟩ ⊥ , ⟨d⟩ ∥ . We use this to define the dynamical regimes in Fig. 1(d). In fact, the interhole distance in the separate-legs configuration quickly evolves linearly in time, ⟨|d|⟩ (q.w.) ⊥ = 13/8(t · τ ). We, therefore, simply define the crossover timescale in Fig. 1(d) as the time τ at which ⟨|d|⟩ (q.w.) ⊥ ≃ 13/8(t·τ ) = ⟨|d|⟩ ⊥ . We, hereby, note that the crossover from the quantum walk to the string oscillation regime for say J ⊥ = 3t, happens already around τ ≃ 1/t. This should be a signifant help to see at least the onset of the oscillation regime in a cold-atom simulation [21]. Figure 9(c), finally, shows the long-time average distances as a function of the trans-leg spin coupling, J ⊥ . For the sameleg configuration, this is, furthermore, shown for indicated values of the intra-leg spin coupling, J ∥ . In Equation (45), we use that the probability to find the holes as nearest neighbors in a given eigenstate n, k is given by the quasiparticle residues |C (n) 2∥ (k). The expressions to the right in Eq. (45), hereby, reveal that the long-time averages of the nonequilibrium average distances are given by an appropriate mean value of the inter-hole average distances, ⟨|d|⟩ (n) ⊥k , ⟨|d|⟩ (n) ∥k , of the eigenstates. One could, therefore, naïvely expect these to scale in the same manner as the eigenstates with t/J ⊥ . For weak correlations, J ⊥ ≫ t, this is indeed the case, where we find that this distance is the same as for the ground states in Fig. 6(c) and, thus, vanishes as (t/J) 2 . For strong correlations J ⊥ ≪ t, however, we see that the distance between the holes reaches a universal t/J ⊥ -scaling. For the same-leg configuration, this also reguires J ∥ ≪ t. This same scaling was found for the motion of a single hole in antiferromagnetic Bethe lattice structures [61], and shows a remarkable qualitative difference to the equilibrium situation with eigenstates supporting only a much weaker (t/J) 1/3 scaling for the eigenstates, Fig. 6(c). This quantifies the qualitative picture drawn from Figs. 8(g)-8(h) that the quenched holes already for intermediate values of J ⊥ ∼ t spread out much more than one would expect from the spatial size of the ground state.
For the computation of the dynamics, we increase the system size N and the number of included bands n max until the results have converged. As a rule of thumb, this is achieved when the system size is a few times larger than the mean distance between the holes. For the most strongly correlated case of J ⊥ = 0.1t, we go up to N = 600 sites and n max = 88 bands. Utilizing the inversion symmetry of the system, we, hereby, need to resolve the energy and residue of N/2 · n max = 26400 states. This emphasizes that we need a very thorough understanding of the eigenspectrum to be able to simulate the quench dynamics in this manner.

V. CONCLUSIONS AND OUTLOOK
Inspired by the recent experimental realization of hole pairing in a cold-atom quantum simulator [27] of a mixeddimensional t-J model [28], we have investigated a simplified setup of Ising spin interactions. This allowed us to determine the exact low-energy single-and two-hole eigenstates. We used this to rigorously show that two holes on separate legs bind strongly to each other in the strongly correlated regime of J ⊥ , J ∥ ≪ t, in that it features a superlinear binding energy: Furthermore, we used this exact description to rigorously account for the nonequilibrium quench dynamics following two initially localized holes at adjacent sites. Similar dynamics has previously been investigated for a single hole in a square lattice geometry [21], whose analysis provided evidence of emergent dynamical regimes, describing the crossover from a quantum walk on short timescales to string oscillations at intermediate timescales and finally the ballistic motion of magnetic polaron quasiparticles at long times [5,67]. In the present mixed-dimensional setup, we found a similar dichotomy of the dynamics for two holes with two major differences. First, the holes are confined to each other, such that their distance remains finite. Second, the string oscillations in the present scenario have an infinite lifetime, and, therefore, persist indefinitely, hindering the long-time equilibration of the system.
These results pave the way for a precise comparison with state-of-the-art cold-atom quantum simulation experiments. There are three essential ingredients that makes the system interesting from this perspective. First, our mixed-dimensional model may be implemented both with fermions and hard-core bosons. Second, the effective interaction strength of 4t/J ⊥ means that the experiments can more easily access a strongly correlated regime already for J ⊥ ≲ 4t. Third, this is particularly relevant for the quench dynamics, where the crossover from the quantum walk to the string oscillations already happens around times of τ ≃ 1/t in this intermediate parameter regime. We, therefore, believe that it should be possible to experimentally access the crossover from the quantum walk to the confinement-induced oscillations.
Furthermore, such experiments also naturally lead to two interesting roads ahead. First, by systematically increasing the number of legs in the ladder, one can carefully analyze if the system supports the formation of stripes [68][69][70][71] inherent to the phenomenology of high-temperature superconductors. Such inquiries were investigated in Ref. [58] using quantum Monte Carlo calculations, in the special case where the transand intra-leg spin interactions are equal. We speculate that our methodology may lend exact insights into this scenario at zero temperature. Second, we believe that it is possible to use the present methodology also at nonzero temperatures. This would require a thorough analysis of the eigenstates as more and more spins are flipped. This would enable the exact determination of the nonequilibrium dynamics of holes at finite temperatures, and could be used to answer whether the holes will deconfine [72] from each other as a result of thermal spin fluctuations.

(B5)
From here, we may, then, calculate the probabilities to find the holes a distance d apart. Since we are not interested in the absolute position of the holes, j, we get