From “weak” to “strong” hole conﬁnement in a Mott insulator

We study the problem of a single hole in an Ising antiferromagnet and, using the magnon expansion and analytical methods, determine the expansion coefﬁcients of its wave function in the magnon basis. In the 1D case, the hole is “weakly” conﬁned in a potential well and the magnon coefﬁcients decay exponentially in the absence of a string potential. This behavior is in sharp contrast to the 2D square lattice where the hole is “strongly” conﬁned by a string potential and the magnon coefﬁcients decay superexponentially. The latter is identiﬁed here to be a ﬁngerprint of the strings in doped antiferromagnets that can be recognized in the numerical or cold atom simulations of the 2D doped Hubbard model. Finally, we attribute the differences between the 1D and 2D cases to the magnon-magnon interactions being crucially important in a 1D spin system.


Introduction
A tendency towards particle delocalization is an ubiquitous phenomenon in quantum mechanics, for it is encoded in the Heisenberg uncertainty principle for momentum and position operators [1]. Perhaps one of its most iconic examples is the so-called particle tunneling under a finite potential barrier: even if the energy of the particle is below the potential amplitude, the probability to find a particle outside the potential well is finite. Yet, the particle is considered localized, for its wave function decays exponentially with an increasing distance from the potential well. As an important feature, the potential acts here only locally and does not increase with the distance. This example of electron localization by a potential well will be called "weak confinement" in what follows. The fact that a particle can delocalize beyond a potential barrier has tremendous implications for the electron wave functions typically found in crystals. It allows for an electron tunnelling under the potential barrier of the periodic potential formed by the ions, leading to the modulus of the electron wave function following the periodicity of the ionic potential [2]. This is the essence of the Bloch theorem and means that an electron in a typical crystal is completely delocalized over all ionic sites.
Nevertheless, localization of electrons in crystals is possible. For instance, this can happen in the Mott insulators-crystals for which strong electron-electron interactions determine electron localization [3,4]. The Mott localization is still not fully-understood and is an area of active research-both for an integer filling [5][6][7][8][9], as well as in doped systems [10][11][12][13][14][15][16][17][18][19][20]. This lack of a complete understanding of the problem is largely due to the fact that the most widelyused models describing the problem (e.g. two-dimensional (2D) Hubbard [21], t-J [22,23] or even the t-J z [4] models) cannot be solved exactly [3] and the wave function of an electron in a Mott insulator is not known in general.
Therefore, here we concentrate on perhaps the simplest, though still nontrivial and realistic, 1 problems of electron localization in the Mott insulator-the problem of the confinement of a particle by an effective potential that takes place when a single hole is added to the ordered ground state of the half-filled t-J z model [24][25][26][27][28][29][30][31][32][33][34][35][36][37][38]. Using an improved version of the recently developed magnon expansion (ME) method 2 [39][40][41][42][43] and analytic calculations, we unambiguously show that a single hole: (i) in 2D or higher dimensional models is "strongly" confined in the ground state and its wave-function coefficients decay superexponentially, i.e., much faster than in the textbook case of a single finite potential well, (ii) in a 1D chain experiences "weak" confinement, i.e., has the wave-function coefficients decaying exponentially, just as in a potential well. Interestingly, as we show below, these differences between the 1D and higher-dimensional cases can be easily understood in the magnon language as originating from the crucial role played by the magnon-magnon interactions in a 1D spin system. Alto-gether, this means that lowering dimensionality and adding interactions may in fact remove "strong confinement" in favor of "weak confinement" in a strongly correlated system.

Model and methods
The Hamiltonian of the t-J z model reads [4]: whereñ i = σc † iσc jσ . It describes the constrained hopping of fermions ∝ t, in the restricted Hilbert space without double occupancies, i.e.,c † iσ ≡ c † iσ (1−n iσ ), along the bond 〈i j〉, and the antiferromagnetic (AF) exchange ∝ J > 0 between the z-th components of the S = 1 2 spins of the localized c-electrons [44,45]. In what follows we study the properties of the eigenstates of the Hamiltonian (1) for the case of a single hole introduced into the half-filled limit (i.e., one electron localized at each site). As in the half-filled case, the ground state of the model (1) is an Ising antiferromagnet; the problem studied here is that of a propagation of a single hole in an Ising antiferromagnet.
The method used in the paper requires first to express the above model (1) in the magnon language. To this end, the magnons are introduced here by means of the Holstein-Primakoff transformation [46] which maps interacting spins on a boson problem, Here i is a magnon creation operator, and the sign ± alternates between the AF sublattices in the Ising antiferromagnet. At the same time, removing a spin generates a hole: for ↑-sublattice:c since removing an inverted spin also requires its realignment; complementary operations generate a hole at ↓-sublattice [24]. Crucially, the maximal number of bosons and holes has to be limited to maximally one at each site i (constraint C1). We note that in Eqs. (2) the aforementioned constraint C1 is implicitly imposed, allowing us to omit the "square root" multipliers on the right hand side of these definitions. As a result, we get an exact representation of the t-J z model in the magnon language for α = 1, This ensures that the terms ∝ J are present only on bonds without holes in the t-J z model (constraint C2). Most importantly, the last term in (4) is the only magnon-magnon interaction (at α = 1) in this Hamiltonian-a constant term active only if two magnons (i.e., inverted spins) are present on the bond 〈i j〉. It is precisely this term which is neglected in the well-known linear spin wave (LSW) theory (at α = 0).
Our primary method, the ME, is a novel, numerical technique of solving the polaronic models, in the Green's function formalism, through the expansion of the equations of motion [40]. It has been successfully applied to several polaronic problems [41][42][43], including spin, orbital, and spin-orbital polarons. The central idea is that the relevant processes are limited to the hole's neighborhood. Thus, one can perform the expansion in real space and apply a Hilbert space cutoff based on the size and spread of the bosonic cloud surrounding the hole. This greatly reduces the Hilbert space while all the states relevant to the dynamics, as well as the C1 and C2 constraints, are included.
On the practical level, the method consists in multiple applications of the Dyson equation where V is the interaction of the problem, and G 0 (ω) is the free single-particle Green's function operator corresponding to the exactly solvable, non-interacting single particle Hamiltonian. Taking the expectation value of the (full) single-particle Green's function operator G(ω) operator in a single particle state, usually the Bloch state and expanding the right-hand-side of Eq. (5) yields a single equation of motion (EOM) for the Green's function.
Since the matrix element of the G 0 (ω) operator is in principle known, the central problem of the expansion is evaluating the effect of the interaction V on the given single particle state. We assume here that this interaction is a hole-boson coupling, which either creates or destroys additional bosons in the system. For instance, acting with the kinetic coupling, as found in our problem, on the Bloch state |k〉 leads to a hole-magnon state with momentum k: Therefore, the expansion procedure introduces higher order Green's functions involving states with different hole-boson spatial configurations. These Green's functions are also unknown and are subject to the same expansion as before. Thus, the variational expansion is controlled by the maximal number of bosons created in the system-the process is continued until the EOM system closes at the desired level of expansion. Once generated, the EOM system can be solved numerically to yield the Green's function G(k, ω), as well as all the other generalized functions appearing in the expansion, which are crucial for recreating the wavefunction and the resulting P n string length probability distributions. We stress that the ME method is a numerically exact method of calculating the Green's function of a particular polaronic model on a particular lattice and for a given number of magnons. Thus, when applied to a hypercubic lattice and once the calculated observables cease to change with increasing number of magnons (i.e., the method "converges"), the obtained result contains all physical processes governing the propagation of a hole in the polaronic model under study; e.g. it can include the quite-often disregarded Trugman processes [26].

Results
The central object that is calculated in this paper is the probability distribution {P n } of observing n magnons in the wave function of a hole doped into the Ising antiferromagnet. This is achieved by calculating the coefficients of the expansion of the ground state wave function of the half-filled t-J z model (1) with a single added hole in the above-explained "magnon language" basis [47] using the ME method in the numerically converged case, which requires keeping up to ca. 100 bosons in the calculations.

1D AF Ising chain
The results, obtained for the "genuine" 1D t-J z model, i.e., when the magnon-magnon interactions are correctly included (α = 1), are presented in Fig. 1(a1). We observe that, irrespectively of the value of the spin exchange J, the probabilities {P n } always decay exponentially with the increasing number of magnons n. This result is further elucidated by the analytic calculations (see Appendix A) which are in perfect agreement with the above 1D numerical results, see Fig. 1(a1). Thus, in the case of the full 1D model (1) we obtain an exponential decay of the probabilities, with the decay length l given by the inverse of the logarithm of the ground state energy of the AF Ising chain with a single hole ( GS ); see Appendix A for the exact expressions for the parameters {A, l, GS }. For completeness, let us note that we could have obtained this exact result also in the spinon language [36]: in this case the ground state of the AF Ising chain with a single hole is most easily understood as a bound state of a mobile holon with a spinon, which arises as a result of the single hole being confined in a 1D potential well. Nevertheless, for a better comparison with the numerics as well as with the 2D case studied below we employ here the magnon language. Interestingly, especially when contrasted with the 2D studies below, the above exponential decay is not obtained when the magnon-magnon interactions are not correctly taken into account in the magnon language expression for the 1D t-J z model, i.e., when α ∈ [0, 1) in Eq. (4). Instead a qualitatively "faster" decay is then clearly observed in the numerical ME calculations, see Fig. 1(b1). Again the numerical behavior is exactly reproduced by the analytically derived values of the probabilities {P n }, see Fig. 1(b1) and Appendix A. In the quantitatively simplest case of α = 0 the expression for {P n } reads: where J is the Bessel function of the first kind and the explicit equations for the constants {B, GS } are given in Appendix A. The expressions for {P n } for any value of α ∈ [0, 1) are qualitatively similar but quantitatively more complex-they are explicitly given by Eqs. (26) in Appendix A. Thus, in the 1D case a first order quantum phase transition can be observed at α = 1. We point out that the decay given by Eq. (9) is well approximated by a Γ function, i.e., the decay is faster than that of an exponential and can be described as a superexponential. The superexponential character of the decay is also visible in the asympotic behavior of ln P n obtained for large n in both of the above-mentioned cases (cf. Appendix A for details): Hence, in the asymptotic limit the superexponential decay of ln P n can be understood as being described by a function decreasing faster with increasing argument than a linear function with a negative coefficient.

2D AF Ising model on a square lattice
The central result of this paper is that the superexponential decay of probability distribution {P n } is found for the hole wave function of the 2D Ising antiferromagnet-not only when the magnon-magnon interactions are neglected but also for the "genuine" t-J z model with all the interactions correctly included. While such a superexponential decay is already visible from the plots presenting the numerical ME results in Figs. 1(c1) and 1(d1), we have further confirmed this behavior by fitting the numerical results with the following approximate expressions of the probability distributions {P n } that is valid for α ∈ [0, 1], see Appendix A: where {z, } are fitting parameters (note that the constant C also depends on {z, }; for more details, including the explicit expression for C, see Appendix A). Again one can obtain the asymptotic behavior for large n of ln P n -which is the same as in Eq. (10), i.e., which further confirms that the 2D case is qualitatively similar to the 1D case without the magnon-magnon interactions. However, we note that a systematic error occurs in the 2D square lattice antiferromagnet when magnon-magnon interaction is neglected, see also below. The approximate expression for the probability distribution {P n } in a square lattice 2D model [Eq. (11)] is motivated by the analytically exact expression obtained in the 1D case without the magnon-magnon interactions, see Eq. (9). Next, the 1D result can be extended to the case of the Bethe lattice for a given arbitrary coordination number z and the respective ground state energy GS . Finally, to account for the fact that the t-J z model is studied on the 2D square lattice and not on the Bethe lattice, we allow for some variation of the coordination number z and the ground state energy GS and leave them as the fitting parametersz and , respectively.

Intuitive understanding: cartoons and effective potential
Let us now try to gain some intuitive understanding of the results presented above. It turns out that this can be rather easily achieved by looking at the cartoon figures, showing the hole propagation in an Ising antiferromagnet in all four studied cases, see Figs. 1(a2)-(d2). Let us first concentrate on probably the easiest case, i.e., on the 1D model without magnonmagnon interactions (α = 0 or the LSW case), cf. Fig. 1(a2). Here, a single hole hop to the nearest-neighbor site (i +2) generates a boson at site (i +1) and thus n i+1 = 1. Since after the previous step also n i = 1, the term J(n i+1 + n i )/2 in (4) increases the energy of this state by J once magnon interactions are absent. Therefore, not only the further the hole moves the more magnons are created, but also an energy cost ω 0 ≡ J is paid after creating each magnon according to the LSW theory, while there is no mechanism to reduce energy due to magnons arising on the hole's path. The well-known linear string potential found in the 2D case [24][25][26][27][28][29][31][32][33][34][35] is clearly observed here, see Fig. 1(a2). Crucially, such a phenomenon is absent once the magnon-magnon interactions are turned on in the 1D model [cf. Fig. 1(b2)]: in that case, even though the hole creates magnons at each step of its motion, there is no energy cost associated with this process as the −αJ n i+1 n i term in the t-J z Hamiltonian (4) cancels completely the energy cost in LSW approximation after creating all but the first magnon. This shows why the magnon-magnon interactions play such a unique role in the 1D case.
The above simple understanding changes to some extent in the 2D case. Here the magnonmagnon interactions are no longer qualitatively relevant, since, unlike in the 1D case, the energy associated with creating a magnon during hole motion cannot be canceled completely by the magnon-magnon interactions, see Fig. 1(c2)-(d2). This shows that the string-like picture [24][25][26][27][28][29][31][32][33][34][35] is valid in the 2D model even when the magnon interactions are included. However, our analysis and the comparison with the ME approach confirms that the string energy evaluated for the 2D square lattice is overestimated in the SCBA and could be corrected by including the effects of magnon-magnon interactions within the modified-SCBA method [31]. De facto, an unphysical part of string is generated on the hole path itself, pretty much the same as in the 1D model. It is removed when the magnon-magnon interactions are included.
The above discussion can be rationalized in the language of the hole being mobile in an effective potential. Then the "genuine" 1D case (i.e., with magnon-magnon interactions correctly included) corresponds to a hole effectively moving in a potential well. On the other hand, all other cases (and in particular the 2D case with the magnon-magnon interactions) can be approximated by a hole being mobile in a (discrete version) of a linear "string-like" potential. Thus, the exponential (superexponetial) decay of the wave functions coefficients can be associated with a hole moving in a potential well (a linear potential), respectively.
Although we just referred above to the discrete version of the linear potential, we would like to emphasize that there exists actually a difference in the asymptotic behavior of {P n } in the discrete and continuous versions of the linear potential. Whereas the former was discussed above [cf. Eq. (12)] and is of the −n ln n form, the asymptotic behavior of the logarithm of the Airy function (which is the ground state wave function of a hole in the continuous linear potential [25]) is different and is given by ∼ − 2 3 n 3/2 . Nevertheless, in both cases the logarithms of both asymptotes "decay faster" than a linear function and are thus understood as superexponential.

Optical lattice experiments
The recent optical lattice experiments [48][49][50], as well as numerical simulations of the 2D doped Hubbard model [51,52], reported on the histograms of the lengths of strings of mis-aligned spins in the ground state hole wave function. As this quantity can be reliably approximated by the aforementioned probability distribution {P n }, this means that a detailed study of the functional form of such histograms can be used to verify to what extent the hole is confined in the 2D Hubbard model.
More precisely, this work allows us to formulate a condition for the observation of a linear string potential acting on the hole in the 2D doped Hubbard (or t-J) models. First, we note that the case with the presence (absence) of the linear string potential is actually naturally defined in the 1D t-J z model with (without) interactions, respectively. 3 Second, as discussed in detail above, once the linear string potential is present (absent) in the latter model, the coefficients of the ground state wave function for a hole in the "magnon language" basis (i.e., {P n }) need to decay superexponentially (exponentially), respectively. Therefore, combining these two observations we can conclude that, if the above-mentioned histograms observed in the optical lattice experiments: (i) showed a superexponential decay with the growing string length l, then this would strongly indicate that a linear string potential indeed plays a dominant role in the hole motion in the 2D doped Hubbard model; (ii) showed merely an exponential decay, then the presence of a linear string potential would be ruled out.
Can we apply the above condition to verify the existence of the linear string potential in the recent experimental or large-scale numerical simulations of the Hubbard model? While the latest optical lattice experiments [48][49][50] still show too few data points to unambiguously conclude whether the observed dependence is exponential or superexponential, we suggest that future optical lattice experiments might be capable of delivering more conclusive data. Moreover, the probability distribution {P n } can in principle be easily calculated in the (future) large-scale numerical simulations of the doped t-J or Hubbard models.

Spectral functions
The above calculations show how the ground state properties of the t-J z model depend on the dimensionality of the problem as well as on the inclusion of the magnon-magnon interactions. This has implications primarily for the optical lattice experiments. However, also the excited states are affected in a somewhat similar manner and this has implications for the understanding of the spectral function A(k, ω) of a number of correlated compounds. Thus, we start by defining the spectral function A(k, ω) in the usual way, where G(k, ω) = 〈Φ 0 |c † kσ (ω − H + E 0 ) −1c kσ |Φ 0 〉 is the momentum-dependent interacting electron Green's function (the spin index σ can be suppressed here), and |Φ 0 〉 is the ground state of the model Eq. (1) in the half-filled limit (Ising antiferromagnet) with energy E 0 . As the spectral function is calculated using the ME method, we first need to express it in the magnon language-the constrained electron Green's function transforms then into the single-particle hole Green's function, G h (k, ω) = 〈Φ 0 |h k (ω − H + E 0 ) −1 h † k |Φ 0 〉. In 1D chain the spectral function A(ω) is always momentum-independent, but as the results obtained using the ME method show, its form depends crucially on the inclusion of the magnon-magnon interactions, see Fig. 2(a). If these interactions are included, which should always be done for a "genuine" representation of the t-J z model, it consists of a single dispersionless δ-like peak at low energy, which indicates a quasiparticle-like state with infinite mass (a bound state), accompanied by an incoherent spectrum at higher energies. On the other Figure 2: Spectral function A(ω) of a single hole in the 1D t-J z model with (α = 1) and without (α = 0) magnon-magnon interactions correctly included in the model calculated using: (a) the ME method, (b) the analytically exact or SCBA approach. Parameters: J = 0.4t, broadening δ = 0.05t; the bound state at the low-energy onset of A(ω) splits off from the continuum at larger J/t (lower δ). Note that in the SCBA calculations for α = 0 (α = 1) constraints C1, C2 are excluded (included), respectively.
hand, when this interaction is switched off in the LSW approximation [by putting α = 0 in (4)], the calculated 1D spectral function is ladder-like, suggesting that the string-like potential builds up [24,28,29,31], see Fig. 2(a). This striking difference between the 1D spectral function, with / without magnon-magnon interactions is also recovered using other methods: First, let us turn to the result obtained for the "genuine" representation of the t-J z model with the magnon interactions correctly included. In this case, taking advantage of an exact analytical result for the t-J z model obtained in the spinon language and using the continued fractions [36][37][38], we observe that the same spectrum is obtained as that of the full model in the ME method, see Fig. 2(b). This shows that the exact result is indeed fully recovered using the magnon language-provided that the magnon interactions are properly taken into account. However, we emphasize that such result is only obtained in the converged case; for the non-converged ME method, i.e., once only a few or tens of bosons are retained in the ME calculations, the spectra consist of several δ-like peaks (not shown).
Second, turning now to the approximate (LSW) representation of the t-J z model in the magnon language, we note that a ladder-like solution is also obtained when the widely-used self-consistent Born approximation (SCBA) method [24] is applied to the polaronic model with the magnon-magnon interactions switched off, see Fig. 2(b). Interestingly, this spectrum differs a bit with respect to the one obtained using the ME method, the reason for this being Figure 3: Spectral function A(k, ω) (13) of a single hole in the 2D t-J z model with (α = 1) and without (α = 0) magnon-magnon interactions correctly included in the model calculated using: (a) the ME method, (b) the SCBA approach following Ref. [31]. All parameters as in Fig. 2 and k = S ≡ (π/2, π/2). Note that in the SCBA calculations for α = 0 (α = 1) constraints C1, C2 excluded (included), respectively. that the above-mentioned constraints C1 and C2 are typically implicitly neglected in the SCBA method (but are automatically taken into account in the ME method). We note that the failure of the LSW approach can be understood by comparing with the well-known 1D result [36] obtained in the "spinon language", for we have shown that one spinon corresponds here to an infinite number of magnons.
A radically different situation occurs on the 2D square lattice, as evident from the spectral functions calculated within the ME, see Fig. 3(a). Here, irrespectively of whether the magnon-magnon interactions are included or not, the spectrum is ladder-like, suggesting that the string-like potential [24][25][26][27][28][29][31][32][33][34][35] develops always in the 2D model: not only in the t-J z model approximated in terms of a polaronic model without magnon interactions (i.e., the LSW approximation) but also in the "genuine" t-J z model, i.e., when magnon-magnon interactions are correctly included. We note that the distances between neighboring maxima found in the ladder spectrum are lower when magnon-magnon interactions are present, see Fig. 3(a), which indicates that the string potential grows then in a slower way with increasing distance. This is indeed confirmed by a slower superexponential decrease of the wave function coefficients {P n }, cf. Figs. 1(c1) and 1(d1), and can also be read off from the cartoon Figs. 1(c2) and 1(d2).
For a 2D square lattice the (converged) ME method is the numerically exact method and, since an analytically exact solution does not exist in this case, it can be used as a benchmark for the other more approximate methods. In particular, the 2D problem can be approximately solved by implementing the SCBA equations derived in Ref. [31]-this time not only without the magnon-magnon interactions and C1, C2 constraints excluded but also with the magnonmagnon interactions and the C1 and C2 constraints included, cf. Fig. 3(b). While the SCBA results qualitatively agree with the ME spectra (e.g. the SCBA also has lower distances between neighboring maxima in the ladder spectrum when magnon-magnon interactions are present), there exist two differences between these two results: (i) the higher energy peaks contain incoherent spectral weight in the ME method whereas they are of delta-like ("quasiparticle") character on the SCBA level; (ii) although the energy of the ground state in the ME method and in the SCBA method (for α = 1 and the "canonical" value of J = 0.4t) is basically the same at k = (π/2, π/2) point (E = −1.58t), there is a small difference between the two results at, e.g., k = (0, 0) point (δE = 0.05t, since according to the ME method the ground state energy reads then E = −1.63t; unshown), in agreement with Ref. [31] which suggests a slight variance between the SCBA and numerical methods once k = (π/2, π/2) and J = 0.4t. We attribute these differences to the important role played by the closed loops (Trugman loops [26]) in the 2D square lattice. 4 The latter, which are neglected by definition within the SCBA, lead to the hole propagating by "cutting the strings" and thus "disrupting" the string potential-which partially destroys the observed ladder-like spectrum by adding a small incoherent weight to the higher energy peaks.
The above results rather naturally lead to the following two questions: First, a relatively important role played by the Trugman processes in obtaining the observed spectrum in 2D means that the the spectral function is momentum-dependent and the hole can be thought as delocalised [26]. Is it then justified to think of the hole in the ground state as being "confined" (as discussed above)? This paradox can be resolved by realising that the Trugman processes almost do not affect the probability distribution {P n } which is of the superexponential type (see discussion in Sec. 3.2) and thus the hole in the ground state can indeed be well-described as being confined in a (discrete) linear potential. Second, one might ask if the physics related to the spectral functions studied above is observable. There exists a number of compounds which can possibly be modeled by the t-J z Hamiltonian-these are predominantly the 1D Ising-like antiferromagnetic chains of BaCo 2 V 2 O 8 [53][54][55] and SrCo 2 V 2 O 8 [56,57], the 2D ferromagnets with alternating orbital order found in K 2 CuF 4 [42,43] and Cs 2 AgF 4 [44,45], and maybe even partially the high-T c cuprates [58]. Thus, we expect that the angle resolved photoemission spectra (ARPES) obtained on these crystals should be qualitatively similar to the spectra shown in Figs. 2-3 (once the magnon-magnon interactions are included). This, however, could have already been predicted using the existing results reporting the spectral functions of the 1D [36][37][38] and 2D [24][25][26][27][28][29][30][31][32][33][34][35] t-J z model. A far more interesting consequence of the present spectral function study is merely theoretical: it shows that the qualitative differences between the 1D and 2D spectral functions originate solely in the different role played by the magnon-magnon interactions in the 1D and 2D models.

Conclusions
In this paper we investigated in detail the particle localization in a Mott insulator that takes place when a single hole effectively moves in a confining potential (hole confinement). The latter appears in the Mott insulating ground state with the Ising antiferromagnetic order and turns out to be of two kinds. In the 1D Ising antiferromagnet the confining potential acts only locally and does not increase with the distance. Hence, the hole is "weakly" confined-the coefficients decay exponentially-just as in the textbook example of a particle localized in a potential well. On the other hand, in the 2D (or higher-dimensional) Ising antiferromagnet, the hole is subject to a (discrete version of the) linear string potential and "strongly" confined, for its coefficients in the magnon language basis decay superexponentially.
The obtained results have two important consequences: First, on the pragmatic side, observation of a superexponential decay of the wave function coefficients in the magnon basis may serve as a fingerprint of the existence of a (long soughtafter) linear string potential in the doped 2D antiferromagnets. This can be performed by a detailed analysis of the recent (and future) optical lattice [48][49][50] or numerical simulations [51,52] of the 2D doped Hubbard model. Second, on the more abstract level, the use of the magnon language not only in 2D but also in 1D model allows us to understand the qualitative differences between the behavior of the hole in the 2D (or higher dimensions) and the 1D Ising antiferromagnets. These differences, which are not only observed in the decay of the hole wave function coefficients but also in the hole spectral functions, are attributed to the crucial role played by the magnon-magnon interactions in one dimension. In fact, in the 1D t-J z model with a single hole a first order quantum phase transition is observed when magnon-magnon interaction can no longer compensate the string potential felt by the mobile hole. 5 Last but not least, this study demonstrates how a particular 1D antiferromagnetic problem can be described in the "magnon language". Such an approach can be seen as being complementary to, e.g., the application of the "spinon language" to magnons in the 2D non-frustrated antiferromagnet [59][60][61][62].

A Derivation of the formulae for the probabilities {P n }
In what follows we derive explicit formulae for the probabilities {P n } of observing n magnons in the single hole ground state of the t-J z model [47] on the Bethe lattice. This result will be later applied to the case of the one-dimensional (1D) chain and to obtain an approximate formula for the P n on a two-dimensional (2D) square lattice.
We start by choosing the basis, B = {|n〉}, where n ∈ , and |n〉 denotes a normalized sum of states with n magnons in a chain connecting the hole with a site at which the hole was created. In this basis the matrix of the t-J z Hamiltonian H [Eq. (1)] for both interacting and non-interacting magnons becomes tridiagonal, Using the Gaussian elimination procedure for a finite matrix and taking the limit n → ∞, one can reduce the equation for the coefficients v n of the ground state vector v corresponding to the ground state energy GS , into a recurrence relation, A.1 1D with α = 1 For the case of the matrix of the 1D t-J z Hamiltonian H with the magnon-magnon interactions included, calculated with respect to the energy E 0 of the half-filled system, we have a 1 = J, a 2 = a 3 = . . . = 3 2 J, b 1 = −t 2, b 2 = b 3 = . . . = −t, so the recurrence relation simplifies to where c 1 = ( GS − J)/(t 2) and c 2 = c 1 / 2. Normalization of the vector v requires |v 0 | 2 to be equal to the residue of the Green's function at the quasiparticle energy GS = 3 2 J − 1 2 J 2 + 16t 2 [30,36], i.e., Finally, we obtain analytical expressions for the probability P n of finding n magnons in the single hole ground state, The above geometric progression can also be expressed in terms of the exponential function, where A = 2|v 0 | 2 and l −1 = 2 ln(|c 2 |).
It is straightforward to observe that the asymptotic behavior for large n is ln P n ∼ − n l .

A.2 1D with α < 1
For the case where the interactions between magnons are not correctly included (α < 1), the coefficients of the matrix become a 1 = J, a n>1 = (1 − α)(n − 2) + 3 2 J, b 1 = −t 2, b 2 = b 3 = . . . = −t. Therefore, the recurrence relation cannot be simplified but one can express the coefficients c n as continued fractions which, due to the linear dependence appearing on the diagonal, can be further expressed (cf. Ref. [29] for details) in terms of the Bessel functions of the first kind [63], The ground state energy GS must satisfy the relation, GS − J = Σ( GS ), with Similarly to the interacting case the normalization of the ground state vector is given by the value of the residue of the Green's function at the pole ω = GS , Here, neither GS nor |v 0 | 2 are given explicitly in the non-interacting case but given an expression for Σ(ω) one can calculate them numerically. Finally, using the recurrence relation for v n , most of the Bessel functions cancel out which leads to a simple formula for the probability P n of finding n magnons in the single hole ground state, We can also approximate the above expressions in terms of the Gamma functions. This is because for 0 < x < n + 1 the Bessel functions can be expanded, leading to J n (2x) x n Γ (n + 1) , (27) so for large n the string length probabilities can be approximated as where Finally, one can obtain the asymptotic behavior for large n of ln P n ∼ −2n ln n. This follows either from the asymptotic behavior of the Bessel function ln J z (2x) ∼ −z ln z or of the Gamma functions ln Γ (z + 1) ∼ z ln z.
Interestingly, the asymptotic behavior for large n of P n is more complex. We obtain: where A.3 2D with α ≤ 1 Before we jump to the case of a 2D square lattice we would like to generalize the above 1D approach to an equivalent problem on the Bethe lattice with coordination number z > 2. Once again, using the basis B, we end up with the tridiagonal form of the matrix of the Hamiltonian. For any α ≤ 1 the coefficients of the matrix are a 1 = zJ 2 , a n>1 = ( z 2 − α)(n − 2) + z − 1 2 J, b 1 = −t z, b 2 = b 3 = . . . = −t z − 1. The equivalent expressions for c n reads, The above result, with z = 4 and α = 1, is equal to the self-energy calculated using Eqs. (21)(22)(23) in Ref. [31]: one merely needs to substitute in Eqs. (21-23) → − 2J. This change is due to the differently defined zero energy level: in Ref. [31] the zero energy level corresponds to the Ising antiferromagnet with one hole whereas in the present paper the zero energy level corresponds to the Ising antiferromagnet.
Similarly to the 1D case the normalization of the ground state vector is given by the value of the residue of the Green's function at the pole ω = GS , Also this time, using the recurrence relation for v n , most of the Bessel functions cancel out which leads to a simple formula for the probability P n of finding n magnons in the single hole ground state on the Bethe lattice,