Higher-Order Spin-Hole Correlations around a Localized Charge Impurity

Analysis of higher-order correlation functions has become a powerful tool for investigating interacting many-body systems in quantum simulators, such as quantum gas microscopes. Experimental measurements of mixed spin-charge correlation functions in the 2D Hubbard have been used to study equilibrium properties of magnetic polarons and to identify coherent and incoherent regimes of their dynamics. In this paper we consider theoretically an extension of this technique to systems which use a pinning potential to reduce the mobility of a single dopant in the Mott insulating regime of the 2D Hubbard model. We find that localization of the dopant has a dramatic effect on its magnetic dressing. The connected third-order spin correlations are weakened in the case of a mobile hole but strengthened near an immobile hole. In the case of the fifth-order correlation function, we find that its bare value has opposite signs in cases of the mobile and of fully pinned dopant, whereas the connected part is similar for both cases. We study suppression of higher-order correlators by thermal fluctuations and demonstrate that they can be observed up to temperatures comparable to the spin-exchange energy $J$. We discuss implications of our results for understanding the interplay of spin and charge in doped Mott insulators.


I. INTRODUCTION
Understanding strongly correlated systems and their emergent phases of matter has been a central task in modern condensed matter physics. Cuprate superconductors have provided a particularly strong motivation for this line of inquiry because of both potential applications of high-temperature superconductivity and unusual thermodynamic and transport properties of these materials 1 . The single-band (Fermi) Hubbard model has been commonly considered as the minimal microscopic model for describing the physics of high-T c curates 2-4 . It correctly reproduces the antiferromagnetic (AFM) Mott insulating state at half filling 5 , and is believed to exhibit non-Fermi liquid properties [6][7][8][9][10][11] as well as d-wave pairing at finite doping [12][13][14][15][16] , though other lattice and materialspecific factors beyond the Hubbard model were also argued to be important [17][18][19][20][21][22][23][24][25] . Methods based on the meanfield approximation are not applicable for analyzing the Hubbard model, because of the presence of many competing instabilities and the importance of quantum antiferromagnetic fluctuations that underlie the emergent non-local attraction between electrons [26][27][28][29] .
Theoretically, the motion of a single hole in a doped Hubbard model takes the form of spin-polaron propagation (sometimes also referred to as magnetic polaron), where the dopant is dressed by a cloud of spin defects in the AFM background. [30][31][32][33][34][35][36][37][38][39][40] Fluctuations of these surrounding spins, correlated with the dopants, should also play a crucial role in the pairing of electrons or holes. [41][42][43][44][45][46] This mechanism has been supported by mul-tiple solid-state experiments. For example, photoemission (ARPES) experiments have revealed the presence of strong correlations, such as the "high-energy anomaly" and the non-quasiparticle features, even in overdoped cuprates [47][48][49] ; recently resonant inelastic x-ray scattering experiments (RIXS) have also observed the persistent spin fluctuations in a wide range of electron and hole doping [50][51][52][53][54] . Therefore, the correlations between spin and dopants generally exist in strongly correlated materials and may be crucial for the observed emergent phases, including also the spin-glass phase at low doping in disordered cuprate compounds 55 .
An additional line of inquiry into the nature of strongly correlated electron systems such as cuprates comes from using impurities. One of the most well-studied examples is the non-magnetic Zn-substitution into the CuO 2 plane, which induces charge impurities due to local chemical potential shifts. With this substitution, nuclear magnetic resonance (NMR) and muon spin relaxation (µSR) experiments found that the magnetic moment is locally enhanced near the Zn impurities [56][57][58][59] . At the same time, the substitution heavily suppresses the d-wave superconductivity [58][59][60][61][62][63] . While the mean-field theory can address this pair-breaking phenomenon [64][65][66][67][68][69] , it appears to provide an apparent contradiction to the spin-fluctuationmediated pairing mechanism: It remains unclear why superconductivity is suppressed while the local moment increases.
As an insightful perspective, we know that it is the dynamical instead of static spin correlations that enhance the d-wave superconductivity [70][71][72] . Therefore, to address this question, one needs to understand the microscopic arXiv:2101.00721v1 [cond-mat.str-el] 3 Jan 2021 correlation between the carrier (electron or hole) and the spin fluctuations near the impurity. This is important not only for cuprates but also generally for all correlated materials: As we will show in this paper, a local charge impurity has outsized effects on high-order spin correlations. However, such a correlation involves a dopant hole/electron and at least two neighboring spins, which cannot be directly accessed in existing solid-state spectral measurements. Fortunately, this task was recently achieved by a distinct approach -the ultracold-atom based quantum simulator 73 . With the help of quantum gas microscopes, the direct measurement of instantaneous correlations, especially those high-order correlations, becomes possible [74][75][76][77][78][79][80][81][82][83][84] . The full spatial resolution of individual lattice sites has also been achieved. Furthermore, quantum simulators have control over dimensionality 85 and can engineer tailored optical potentials to precisely control model Hamiltonians 86,87 . As an example, an optical tweezer can be used to engineer a localized potential for holes, which allows for continuous control of a local impurity in a clean system 88 . This paradigmatic scenario of a single disordered lattice site with a tunable on-site potential will be analyzed in detail throughout this work.
Here, to systematically understand the high-order spin-hole correlations in the doped Hubbard model and the impact of an impurity on the correlations, we present an exact diagonalization calculation of t−J−3s and Hubbard models. Compared to the well-known spin polaron dressing in the case of a mobile dopant in a Mott insulator, we find that the presence of a local pinning field of variable strength results in significantly different spin correlation patterns. Such a difference is reflected in both third-order and fifth-order spin-hole correlation functions, suggesting a crossover from a spin-polaron surrounded by weakened magnetic correlations to a geometric defect surrounded by enhanced magnetic fluctuations. This work also complements our analysis of the fifthorder correlations of unpinned holes in Ref. 89. Since the proposed correlation functions are composed of local spin and charge operators, they can be measured directly using the state-of-the-art quantum gas microscopy 88,90 , which can ultimately extend the quantitative conclusion to the thermodynamic limit. This paper is organized as follows. We first introduce the models and methods used in Sec. II. Then we investigate the third-order and the fifth-order correlation functions in Sec. III and Sec. IV, respectively. We conclude our study in Sec. V.

II. MODELS AND METHODS
As the simplest microscopic model depicting Mott physics in a strongly correlated electronic system, the 2D (Fermi) Hubbard model is described by the Hamiltonian whereĉ † iσ (ĉ iσ ) andn iσ denotes the creation (annihilation) and density operator at site i of spin σ; t ij is the hopping, restricted here to nearest-neighbors (nn) t ij = t. As mentioned above, the single-band Hubbard model is believed to capture the essential physics of strongly correlated materials such as cuprates 2,3 and can be precisely simulated by the cold-atom experiments [81][82][83][91][92][93][94][95] .
The Hubbard model's Hilbert space dimension is relatively large, limiting both the system size and the number of states accessible in exact numerical calculations. To investigate the temperature dependence in the regime accessible by cold-atom experiments (typically T ∼ 0.5t), we also consider the low-energy approximation of the Hubbard model near half-filling -the t−J-like spin model. Through a t/U expansion to the lowest order, one can simplify the Hubbard model to the t−J model [96][97][98][99] .
The constrained fermionic operators acting in the Hilbert space without double occupancy are defined asc † iσ =ĉ † iσ (1 −n iσ ). Although sometimes negligible, there is a 3-site term at the same lowest order (∼ t 2 /U ) that contributes to the dopant's motion 38,39,[98][99][100][101][102][103][104][105][106][107][108] . This defines the t−J−3s model with the Hamiltonian given by H t−J−3s = H t−J + H 3s : We use both the Hubbard and t−J−3s models to study high-order correlations, where the latter allows for access to higher temperatures. Comparisons between the two models allow identifying which effects are captured by effective spin models and do not require higher-order expansions in t/U . Unless otherwise specified, we use U = 8t for the Hubbard model and accordingly J = 0.5t for the t−J−3s model throughout this paper.
To test the correlations with respect to an immobile hole, we also consider an extra pinning potential in the origin (r 0 = 0), controlling the mobility of the single hole. The Hamiltonian becomes: where H 0 is the Hubbard or t−J−3s Hamiltonian, and V is the strength of a local pinning potential. This pinning potential can be realized experimentally by an optical tweezer in an ultracold atom system 88 . Hence, the high-order correlations analyzed in this paper and Ref. 89 are directly accessible to these experiments. Note that the strength of spin-exchange and 3-site terms in the t−J model involving a virtually doubly-occupied central site (which is perturbed by the added pinning potential) must be modified when |V | becomes comparable to the Hubbard interaction U . Throughout this work, we neglect such corrections and focus on the regime |V | ≤ U . The difference between Hubbard and t−J−3s models are quantitatively compared and discussed whenever relevant.
To resolve the wavefunction and the corresponding correlators over a range of low temperatures in a 2D D 4 symmetric system, we perform exact diagonalization calculations on a 4×4 cluster with periodic boundary conditions. Throughout this paper, we evaluate the expectation values of observables in a canonical ensemble. For convenience, we denote the expectation as the thermal average in which Z is the partition function. The n max sets the numerical truncation of excited states, which satisfies E nmax − E 0 T for all temperatures considered in our work. In the results presented in Secs. III and IV, we keep n max ∼ 650 states for the Hubbard model, giving reliable results up to T ∼ 0.4t, whereas n max ∼ 13, 000 states for the t−J−3s model, giving reliable results up to T ∼ t.
In this paper, we employ the parallel Arnoldi and Paradeisos algorithm to determine the equilibrium groundstate wavefunction and expectation values. 109,110 Unless otherwise indicated, we include all total-S z sectors in the thermal ensemble. Part of the results in Sec. IV are benchmarked using the density matrix renormalization group (DMRG) at zero temperature in a 6 × 12 cluster with cylindrical boundary condition, to investigate possible artifacts caused by the finite system size.

III. THE THIRD-ORDER SPIN-HOLE CORRELATIONS
To describe the spin correlations with respect to a doped hole, we consider the following third-order correlation function: where the hole density operatorn h r = (1 −n r↑ )(1 −n r↓ ). In this paper, we focus on the |d| = 1 (neighboring spins) and |d| = √ 2 (diagonal spins), allowing the relative distance between the hole and the spin "bond" |r − r | to vary (see Fig. 1). Figure 2 gives an overview of the spatial (r − r) distribution of the B(r, r , d) in a single-hole-doped Hubbard model at zero temperature, with r corresponding to the coordinate where the pinning potential V is applied (white dot). For different pinning potentials, the nearestneighbor (|d| = 1) correlations are always negative, while the diagonal (|d| = √ 2) correlations are mostly positive due to the strong AFM order. However, the spatial variation of these correlators indicates different underlying physics triggered by the pinning potential, as discussed below.
We first consider systems with a mobile hole, i.e., for V = 0. As shown in the left panel of Fig. 2(a), the (absolute value of the) third-order correlator B(r, r , d) is weakened for shorter distances |r − r |. This spatial distribution can be understood in terms of the spin polaron 111,112 , where a dopant's motion is dressed by a cloud of spin defects, forming a polaronic quasiparticle. As the system is translationally invariant for V = 0, the correlator B(r, r , d) depicts the concentration of spin defects in the co-moving frame of the mobile hole. Due to the small radius of the spin polaron (local screening of spin correlations), the amplitudes of correlators for large |r − r|'s asymptotically approach those of an undoped AFM system.
In contrast to the free hole, the spin correlations in the presence of a pinning potential V > 0 lead to significantly different patterns. With the pinning potential applied at site r, we examine the the third-order correlator B(r, r , d) between the hole (at r) and the spins at r and r + d. As shown in the right two panels of Fig. 2(a) and Fig. 2(b), B(r, r , d) becomes stronger for shorter distances |r − r | (i.e., for |d| = 1 it becomes more negative, while for |d| = √ 2 it becomes more positive). To distinguish from the screening of correlations in the spin-polaron picture, we denote this phenomenon as "anti-screening". A similar phenomena has recently been observed in experiments, albeit at higher temperatures, where it has been attributed to imperfections of the pinning potential 88 . Indirectly, this phenomenon is also consistent with the increased local moments measured in NMR and µSR experiments in Zn-substituted cuprates [56][57][58][59] . One can intuitively understand this phenomenon by observing the fact that spin correlations in systems with a lower coordination numbers (e.g., a 1D chain) are stronger than in the 2D plane with otherwise identical parameters due to the existence of reduced frustrations (an AFM state favors spin to form singlets with all its nearest neighbors) 113,114 . Thus, a system with a strong pinning potential tends to form a 1D edge at the boundary surrounding the impurity. Here the system can lower its energy by forming stronger singlet bonds closer to the impurity, while retaining the bulk spin order at further-away sites. In addition to the change in spatial distribution, the presence of a pinning potential also leads to an overall enhancement of all correlators (in terms of absolute values).
To quantify the different distributions of the thirdorder correlators, we define the difference between the nearest-neighbor correlators as [see the illustration in the inset of Fig. 2 Through the difference, the uniform (positive or negative) background in the "bare" correlator B(r, r , d) is removed. Thus, we refer to the ∆B nn as differential thirdorder correlators. As shown in Fig. 2(b), ∆B nn is always negative for V = 0, indicating the screening effect of the spin polaron; the ∆B nn turns to positive for a finite V , indicating the onset of anti-screening 115 . As such, the crossover between screening and anti-screening natures is clearly reflected by the differential third-order correlators, while the original B(r, r ,ŷ)s do not have a sign change. The reason is that the original ones contain a substantial contribution from the lower-order correlators (i.e., Ŝ z rŜ z r+ŷ ), which is sizably negative. A spatial difference with respect to the hole, however, cancels (or at least heavily reduces) this lower-order background.
The above observation indicates that a "genuine" correlator is required while describing the underlying manybody physics. Such a genuine third-order correlator can be more intuitively defined as the connected part of B(r, r , d) where N denotes the number of lattice sites. As such, the B c (r, r , d) reflects the net hole-spin-spin correlation on top of the AFM background, without the need to extract through a spatial difference. As shown in Fig. 3(a), the B c (r, r , d)'s for small |r − r | flip sign in the presence of the pinning potential, which allows for a better discrimination of the screening and anti-screening regimes.
Given that thermal fluctuations are expected to disrupt long-range ordering, the crossover between the screening and anti-screening situations should depend on both temperature and the pinning potential. To address this, we calculate the T -and V -dependence of B c (r, r , d). As mentioned in Sec. II, the Hilbert space size of the Hubbard model limits the accessible maximal temperature; therefore, we also calculate the above correlators for the t−J−3s model. As shown in Figs. 3(b) and (c), the ∆B nn characterizes a crossover from a screening (negative) to the anti-screening (positive) regime for any temperature T < 0.3t. The critical V value increases with temperature, because thermal fluctuations help with the hole's delocalization and, therefore, its mobility becomes larger for the same strength of pinning potential.
By comparing the Hubbard and t−J−3s models in Fig. 3(c), we note that the results obtained from both models agree quantitatively at the small V side, suggesting that the third-order spin-hole correlations originate from hole motion in an AFM background instead of any interactions beyond the O(t 2 /U ). On the large V side, B c (r, r , d) displays slightly different temperature dependence for two models, which can be attributed to the modified exchange coupling in the vicinity of the pinning site: Between this site and its nearest neighbors, the superexchange coupling becomes J eff = JU 2 /(U 2 − V 2 ) > J. At much higher temperatures (T > J), thermal fluctuations smear out the AFM spin correlations and thus suppress the third-order correlator.

IV. THE FIFTH-ORDER RING-SPIN CORRELATORS
The third-order correlators B(r, r , d) reflect the spin fluctuations around the hole. These correlators are expected to be weak for mobile holes because the underlying motion of the dopants admixes spins from different sub-lattices in the surrounding AFM 111,116 . To minimize the effect of averaging over different trajectories and provide deeper insights into the underlying spin-charge correlation, one may consider constructing higher-order correlators. To this end, we examine the (fifth-order) ringspin correlators with respect to a doped hole, which is (9) As illustrated in Fig. 4, this correlator reflects the ringspin correlations in a co-moving frame of the doped hole. By including non-zero displacements d (between the ringcenter and the hole), Eq. (9) generalizes the five-point correlators introduced in Ref. 89.
We first investigate the dependence of C ♦ (r, d) on the ring displacement d, and consider a mobile hole with no pinning potential (V = 0). In this case, the system has translational invariance and C ♦ (r, d) is only a function of d. As shown in Fig. 5 (a) and the top-left panel of (b), C ♦ is positive for any |d| > 1, but becomes sizable and negative for d = 0. This behavior is reminiscent of a Z 2 Gauss law expected in a string description of the single spin polaron: when the mobile hole forms a spin polaron, its motion creates a "string" of displaced spins; In the C ♦ correlator, spins which are part of the string originate from the opposite sub-lattice and contribute a negative sign 89 . Since C ♦ (r, 0) involves only one site belonging to the string, it is expected to be negative; whereas C ♦ (r, d = 0) has statistically more contributions from an even number of spins affected by a string, displaying a positive value. We find they reach a quantitative agreement by comparing the t−J model and the Hubbard model with different system sizes. This indicates that the observed correlations can be understood from the spin-exchange picture, while higher-order terms in the t/U expansion can be ignored. The consistency between the DMRG calculation in a sizeable system and the ED calculation in a 4×4 cluster further precludes the influence of system size and models.
We then consider systems with a finite pinning potential V in the center at site r 0 . As shown in the up-per panel of Fig. 5(b), the on-site (d = 0) correlator C ♦ (r 0 , 0) flips sign with a V -field. This is a direct consequence of a hole becoming immobile: when the pinning potential traps the dopant hole, the spin polaron breaks down, and the remaining system becomes a half-filled AFM system with a static defect. Since next-nearest spins are expected to be parallel in a pure Néel state, the C ♦ (r, 0) becomes positive in this immobile case. This difference induced by the hole's mobility further reflects the screening to anti-screening crossover, consistent with the observation using the third-order correlators in Sec. III.
The sensitivity to a pinning potential is affected by the thermal fluctuation at finite temperature. As shown in the lower panel of Fig. 5(b), the on-site correlator C ♦ (r 0 , 0) remains negative until V ∼ 2.5t at a higher temperature T = 0.4t. As mentioned in Sec. III, this is a consequence of the hole's delocalization enhanced by thermal fluctuations. To better visualize this crossover, we extract the C ♦ (r 0 , 0) and plot its T and V dependence in Figs. 6(a) and (b). As temperatures rise, the regime dominated by the spin-polaron screening (blue) extends to a larger pinning potential. We observe a similar crossover here, as in Fig. 3; note that the crossover's exact position may be different for different observables as it is not a phase transition. Except for this crossover, we find that all C ♦ (r 0 , d)s decrease with increasing temperature, which occurs dramatically near T ∼ J and indicates the thermal melting of the surrounding antiferromagnet. Given the presence of four spins in the fifthorder correlator, its temperature dependence is more evident than the third-order one presented in Fig. 3(c).
As such, the temperature dependence of this fifthorder correlator reflects the screening to anti-screening crossover similar to the third-order one. The difference is that this phenomenon has already been reflected in the "bare" fifth-order correlator, without spatial derivative or subtracting disconnected parts. In contrast, the genuine (connected) fifth-order correlator reflects a different level of information. To extract the genuine fifth-order correlators, we subtract off the disconnected pieces of the correlator following the same manner as Ref. 89: where the · · · c denotes the connected correlator with lower-order terms subtracted. Figures 6(c,d) show the dependence of the connected correlator C c ♦ (r, 0) on the temperature and pinning potential. We find that this correlator is sizably negative even for large pinning potentials. Although the char- acteristic temperature, above which |C c ♦ | drops dramatically, decreases slightly for the systems with pinned holes, the absolute value of C c ♦ does not change much with V . This observation indicates that the genuine fifth-order correlator is always present, significantly deviating from the classical Néel state regardless of whether the hole is mobile. Such an intrinsic nature is invisible in lowerorder correlators.
In the remainder of this section, we elucidate the physical origin and properties of the genuine fifth-order ringspin correlators. We discuss the role of different statistical ensembles on the fifth-order ring-spin correlators.
To this end, we compare a spin-imbalanced ensemble with definite total spin S z = 1/2 to a spin-balanced statistical mixture with Ŝ z = 0. Alternatively, in the former case, we could consider an ensemble with definite S z = −1/2, which gives the same results: the higherorder correlators C ♦ and C c ♦ each involve an even number of spin operatorsŜ z . Moreover, the underlying spin or Hubbard Hamiltonian is invariant under a global spinflip (S z j → −S z j ). Hence, the correlators in the sectors with definite S z = ±1/2 are equivalent, With the same argument, in the spin-balanced ensemble, Ŝ z = 0, all lower-order correlators involving an odd number of spin operators vanish. While an ensemble with definite S z = 1/2 or −1/2 is hard to realize in a solidstate system, both spin ensembles can be addressed experimentally with ultracold atoms and with full spin and charge resolution 117,118 . This ensemble can be achieved by post-selecting experimental data with specific S z val-ues.
The two spin ensembles correspond to different ways of taking the zero-temperature limit. The usual canonical ensemble converges to a balanced mixture of the two sectors S z = +1/2 and S z = −1/2 when T → 0. In the resulting spin-balanced ensemble Ŝ z = 0, and all but the second sum on the right-hand side of Eq. (10) vanish. In contrast, for the ground state with a definite S z = 1/2, all terms in Eq. (10) contribute.
Nevertheless, for a mobile hole we find that the calculated connected correlators do not change significantly if switched to the imbalanced ensemble: In Fig. 7(a), we plot C c ♦ | S z =0 as a function of t/J (see Fig. 1 in 89) and compare it to the spin imbalanced ensemble. Qualitatively, the same behavior is found for all considered values of t/J. The deviations are largest for small values of t/J, where we expect a smaller spin-polaron radius in the ground state 111 .
The situation becomes slightly different when the hole is trapped by a pinning potential. In Fig. 7 we compare connected and bare correlators for a spin-imbalanced ensemble with the fixed total S z = 1/2. The bare correlator C ♦ in Fig. 7(b) reflects a crossover similar to Fig. 6(a), with slightly stronger temperature dependence. This sensitivity to temperature results from the fact that the fixed total spin S z = 1/2 gives finite contributions to the odd terms in Eq. (10) as lower orders and thus induces more fluctuations when temperature increases. For the same reason, the connected fifth-order correlator C c ♦ also displays a slightly stronger sensitivity to the pinning potential V s [see Fig. 7 (b)].
Most significantly, we observe in Fig. 7(c) that the connected correlator C c ♦ decreases in magnitude in the FIG. 7: (a) Genuine fifth-order correlator C c ♦ for different spin-ensembles, Ŝ z = 0 and S z = ±1/2, at zero temperature and without pinning potential. The results are obtained by the DMRG simulations for a t−J Hamiltonian on a 6 × 12 cylinder. (b) The bare fifth-order correlator C ♦ and (c) the genuine correlator C c ♦ in a spin-imbalanced ensemble with fixed total S z = 1/2, as a function of pinning potential V and temperature T . We consider a t−J−3s Hamiltonian. spin-imbalanced ensemble for large V and small temperatures T ; in contrast, it remains sizably negative in this regime for the spin-balanced ensemble [see Fig. 6(c)]. We attribute this influence of the spin ensemble on the connected higher-order correlator to the formation of Néel order and a non-zero staggered magnetization, M z = r (−1) rx+ry Ŝ z r = 0, in the lattice. In a finite-size system, like the one we study, the latter requires spinimbalance and a pinning potential for the hole to avoid mixing of A-and B-sublattices in the co-moving frame of the hole where we define the higher-order correlators. When a non-zero staggered magnetization forms in the co-moving frame with the hole, the lower-order termse.g. Ŝ z r -provide sizable non-vanishing contributions to the bare ring-spin correlators. Since we observe that bare ring-spin correlators are only weakly affected by different spin ensembles, the genuine higher-order terms can differ significantly for different ensembles, as we find in Fig. 7(b).
Finally, to further illustrate the differences between spin-balanced and imbalanced ensembles, we consider a strongly simplified physical setting at zero temperature. We assume that the hole is fully pinned in the center of the system (V t) and take into account only Ising couplings ∼ J zŜz iŜ z j between the spins. The ground state of this toy model in the spin-imbalanced ensemble with S z = 1/2 is a classical Néel state |N, + where a missing down-spin at r defines the pinned hole. Similarly, the ground state in the spin-imbalanced ensemble with S z = −1/2 is the opposite Néel state |N, − , obtained from |N, + by flipping all spins. For both Néel states, the bare ring-spin correlator C ♦ (r, 0)| S z =1/2 = C ♦ (r, 0)| S z =−1/2 is equal since all four spins surrounding the hole are aligned. Hence, C ♦ (r, 0)| Ŝz =0 = C ♦ (r, 0)| S z =±1/2 also takes the same value in the spin-balanced ensemble, defined by the even statistical mixture {|N, + , |N, − }. Now we turn to the connected correlators. For each of the two Néel states |N, ± taken individually, C c ♦ | S z =±1/2 = 0 vanishes: By construction, the higher-order connected correlators are identically zero for product states. To obtain this result from Eq. (10), it is important to account for all lower-order contributions, including those with an odd number ofŜ z operators. In contrast, in the spin-balanced ensemble where Ŝ z = 0, only lower-order terms with an even number ofŜ z operators contribute in Eq. (10). This leads to a strongly modified connected correlator, with an overall negative value: indeed, counting all terms for the classical Néel states gives the estimate C c ♦ (r, 0)| Ŝz =0 = −2 in our toy model with Ising interactions.

V. DISCUSSIONS AND CONCLUSIONS
Using the third-order and fifth-order spin-hole correlators, we have analyzed spatial correlation of a doped hole and its surrounding spin fluctuations in the single-hole doped Hubbard and t−J−3s models. We particularly investigate the impact of impurity on these high-order correlations, mimicked by a localized pinning potential with varying strength V . Interestingly, we find that even an extremely local impurity imposes an outsized effect on the high-order correlations. With the increase of the pinning potential, we identified a crossover from screening to anti-screening of the spin fluctuations surrounding the hole: For a mobile hole, the weakened third-order correlators in proximity to the dopant and the negative fifthorder ring-spin correlators reflect the underlying spin po-laron formation -here the motion of the hole screens the spin defects; for an immobile hole trapped by the pinning potential, the third-order correlators in proximity to the dopant are enhanced and the bare fifth-order ring-spin correlators are positive -as a result of the anti-screened magnetic moment of the pinned hole. In this case, a geometric defect and a staggered field are induced in the AFM background. This microscopic transition from a spin-polaron to an anti-screened defect provides a route to understand the Zn-substitution experiments, where the magnetic moment is enhanced while superconductivity is suppressed. Due to the breaking of the spinpolaron, different holes are no longer paired by sharing the spin fluctuations. Instead, the locally enhanced spin fluctuations surrounding, the impurity-trapped hole, indicates a repulsion to other holes and therefore suppresses the coherence of Cooper pairs. This crossover in highorder correlators may also provide intuitions to impurity effects in other correlated materials.
We have also examined the impact of temperature on the considered spin-charge correlations. We find that thermal fluctuations lead to an extended spin-polaron regime with a weakly pinned hole, and identify the crossover to the pinned hole as a function of both temperature T and pinning potential V . In contrast to the crossover, we found genuine (i.e. connected) fifthorder ring-spin correlators carrying different underlying physics: they are sizably negative and robust against the pinning potential.
Based on these observations, we demonstrated that higher-order correlation functions provide a new perspective on quantum many-body systems with strong correlations. This is a particularly promising direction in the context of correlated fermionic systems such as the Fermi-Hubbard model, where traditional theoretical approaches are limited and a comprehensive physical picture is still lacking. The intrinsic properties of this type of quantum many-body systems may be invisible in traditional two-point correlations of solid-state measurements. The specific situation discussed here, where we analyzed the effect of a localized pinning potential on mobile dopants, is particularly relevant in the solidstate context for understanding the effects of disorder on strong spin-charge correlations. On the other hand, the higher-order correlators we propose to measure are directly accessible in state-of-the-art quantum gas microscopy experiments with ultracold atoms in optical lattices. The required temperatures to observe non-trivial effects have already been reached, and we expect that even lower temperatures will become accessible in the near future.
The differential third-order correlator has been briefly discussed in Sec. III. To reveal the impact of temperature and the pinning potential, we present its V -T dependence of ∆B nn in Fig. 8, following the same manner of Fig. 3(c) in the main text.
In addition to the ∆B nn defined in Eq. 7, we further quantify the difference between the diagonal (hole-spin) correlators as ∆B diag = B(r, r +x,ŷ −x) − B(r, r + 2x,ŷ −x) . (A1) The spatial relations of both differences are illustrated in the inset of Fig. 8. Note, to compensate for the negative sign induced by the AFM sublattices, we swapped the order while defining these two differences.
We notice a difference between the ∆B nn [ Fig. 8(a)] and the ∆B diag [ Fig. 8(b)]: although the screening effect has been dramatically suppressed for large pinning potentials following the similar crossover, ∆B diag does not exhibit an evident sign change. Since spin correlations are inhomogeneous when the system breaks translational symmetry regardless of the hole, the differential correlator ∆B diag has been affected by this background inhomogeneity, reducing the signal-to-noise ratio for describing the underlying properties of the spin-hole composite.
The advantage of the genuine correlator B c (r, r , d) over the differential correlator ∆B diag is reflected by comparing Fig. 3(c) and Fig. 8(b). Although both describe the diagonal spin correlations relative to the hole, the genuine correlator clearly separates the spin polaron (where it is sizeably negative) and the anti-screening regime (where it is sizeably positive). Especially in the large V limit, the genuine correlator characterizes a narrower crossover compared with Fig. 3(a), due to the exclusion of noise in lower-order correlators.