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

Yao Wang<EMAIL_ADDRESS>Department of Physics and Astronomy, Clemson
University, Clemson, South Carolina 29631, USA  Annabelle Bohrdt  Department
of Physics and Institute for Advanced Study, Technical University of Munich,
85748 Garching, Germany  Munich Center for Quantum Science and Technology
(MCQST), 80799 München, Germany  Joannis Koepsell  Max-Planck-Institut für
Quantenoptik, 85748 Garching, Germany  Munich Center for Quantum Science and
Technology (MCQST), 80799 München, Germany  Eugene Demler  Department of
Physics, Harvard University, Cambridge 02138, USA  Fabian Grusdt
<EMAIL_ADDRESS>Department of Physics and Arnold
Sommerfeld Center for Theoretical Physics (ASC), Ludwig-Maximilians-
Universität München, Theresienstr. 37, München D-80333, Germany  Munich Center
for Quantum Science and Technology (MCQST), 80799 München, Germany

###### Abstract

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.

###### pacs:

xxx

##  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 Dagotto (1994). The single-band (Fermi) Hubbard model has been
commonly considered as the minimal microscopic model for describing the
physics of high-$T_{c}$ curates Anderson (1987); Zhang and Rice (1988); Lee et
al. (2006). It correctly reproduces the antiferromagnetic (AFM) Mott
insulating state at half filling Singh and Goswami (2002), and is believed to
exhibit non-Fermi liquid properties Zheng et al. (2017); Huang et al. (2017);
Ponsioen et al. (2019); Kokalj (2017); Huang et al. (2019); Cha et al. (2020)
as well as $d$-wave pairing at finite doping Maier et al. (2005); Zheng and
Chan (2016); Ido et al. (2018); Jiang and Devereaux (2019); Chung et al.
(2002), though other lattice and material-specific factors beyond the Hubbard
model were also argued to be importantOhta et al. (1991); Lanzara et al.
(2001); Shen et al. (2004); Mishchenko and Nagaosa (2004); Cuk et al. (2004);
Reznik et al. (2006); Johnston et al. (2012); Peng et al. (2017); He et al.
(2018). Methods based on the mean-field 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 electronsDemler et al.
(2004); Scalapino (2012); Maier et al. (2016); LeBlanc et al. (2015).

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.Kane et al. (1989); Sachdev (1989); Liu and Manousakis (1991);
Martinez and Horsch (1991); Stephan and Horsch (1991); Manousakis (1991);
Dagotto et al. (1992); Liu and Manousakis (1992); Stephan and Horsch (1992);
Bała et al. (1995); Brunner et al. (2000) Fluctuations of these surrounding
spins, correlated with the dopants, should also play a crucial role in the
pairing of electrons or holes.Scalapino et al. (1986); Gros et al. (1987);
Kotliar and Liu (1988); Schrieffer et al. (1989); Scalapino (1995); Tsuei and
Kirtley (2000) This mechanism has been supported by multiple 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 cupratesGraf et al. (2007); He
et al. (2019); Chen et al. (2019); recently resonant inelastic x-ray
scattering experiments (RIXS) have also observed the persistent spin
fluctuations in a wide range of electron and hole dopingLe Tacon et al.
(2011); Dean et al. (2013a, b); Lee et al. (2014); Ishii et al. (2014).
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 Hasselmann et al. (2004).

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 CuO2 plane,
which induces charge impurities due to local chemical potential shifts. With
this substitution, nuclear magnetic resonance (NMR) and muon spin relaxation
($\mu$SR) experiments found that the magnetic moment is locally enhanced near
the Zn impurities Mahajan et al. (1994); Alloul et al. (1991); Xiao et al.
(1990); Mendels et al. (1994). At the same time, the substitution heavily
suppresses the $d$-wave superconductivity Xiao et al. (1990); Maeda et al.
(1990); Mendels et al. (1994); Fukuzumi et al. (1996); Nachumi et al. (1996);
Basov et al. (1998). While the mean-field theory can address this pair-
breaking phenomenon Lee (1993); Sun and Maki (1995); Salkola et al. (1996);
Hirschfeld and Putikka (1996); Franz et al. (1996); Tsuchiura et al. (1999),
it appears to provide an apparent contradiction to the spin-fluctuation-
mediated 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 Kivelson
et al. (1998); Emery et al. (1999); Kivelson et al. (2003). Therefore, to
address this question, one needs to understand the microscopic 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 Gross and Bloch (2017).
With the help of quantum gas microscopes, the direct measurement of
instantaneous correlations, especially those high-order correlations, becomes
possibleBakr et al. (2009); Sherson et al. (2010); Parsons et al. (2015);
Cheuk et al. (2015); Omran et al. (2015); Edge et al. (2015); Haller et al.
(2015); Mazurenko et al. (2017); Hilker et al. (2017); Chiu et al. (2019);
Schweigler et al. (2017). The full spatial resolution of individual lattice
sites has also been achieved. Furthermore, quantum simulators have control
over dimensionality Salomon et al. (2019) and can engineer tailored optical
potentials to precisely control model Hamiltonians Brown et al. (2019);
Nichols et al. (2019). 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 systemKoepsell et al. (2019). 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\small{-}J\small{-}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 fifth-order correlations of unpinned holes in
Ref. Bohrdt et al., 2020. 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 Koepsell et al. (2019);
Koepsell et al. (2020a), 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

$\displaystyle\mathcal{H}_{H}$ $\displaystyle=$ $\displaystyle-\sum_{{\bf
i},{\bf j},\sigma}\left(t_{\bf ij}\hat{c}^{\dagger}_{{\bf
j}\sigma}\hat{c}_{{\bf i}\sigma}+h.c.\right)$ (1) $\displaystyle+U\sum_{{\bf
i}}\left(\\!\hat{n}_{{\bf i}\uparrow}-\frac{1}{2}\right)\left(\hat{n}_{{\bf
i}\downarrow}-\frac{1}{2}\right),$

where $\hat{c}^{\dagger}_{{\bf i}\sigma}$ ($\hat{c}_{{\bf i}\sigma}$) and
$\hat{n}_{{\bf i}\sigma}$ denotes the creation (annihilation) and density
operator at site i of spin $\sigma$; $t_{\bf ij}$ is the hopping, restricted
here to nearest-neighbors (nn) $t_{\bf\langle ij\rangle}\\!=\\!t$. As
mentioned above, the single-band Hubbard model is believed to capture the
essential physics of strongly correlated materials such as cuprates Anderson
(1987); Zhang and Rice (1988) and can be precisely simulated by the cold-atom
experimentsHofstetter et al. (2002); Ho et al. (2009); Jördens et al. (2010);
Esslinger (2010); Mazurenko et al. (2017); Hilker et al. (2017); Tarruell and
Sanchez-Palencia (2018); Chiu et al. (2019).

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\sim 0.5t$), we also consider
the low-energy approximation of the Hubbard model near half-filling – the
$t\small{-}J$-like spin model. Through a $t/U$ expansion to the lowest order,
one can simplify the Hubbard model to the $t\small{-}J$ model Chao et al.
(1977, 1978); Belinicher and Chernyshev (1994); Belinicher et al. (1994).

$\mathcal{H}_{t\\!-\\!J}=-t\\!\sum_{\langle{\bf i},{\bf
j}\rangle,\sigma}\\!\left(\tilde{c}^{\dagger}_{{\bf j}\sigma}\tilde{c}_{{\bf
i}\sigma}\\!+\\!h.c.\right)\\!+\\!J\sum_{\langle{\bf i},{\bf
j}\rangle}\left(\textbf{S}_{\bf i}\\!\cdot\\!\textbf{S}_{\bf
j}-\frac{\hat{n}_{\textbf{i}}\hat{n}_{\textbf{j}}}{4}\right),$ (2)

where $\textbf{S}_{\bf i}\cdot\textbf{S}_{\bf j}\\!=\\!\hat{S}_{\bf
i}^{z}\hat{S}^{z}_{\bf j}\\!+\\!\frac{1}{2}\left(\hat{S}_{\bf
i}^{+}\hat{S}_{\bf j}^{-}+\hat{S}_{\bf i}^{-}\hat{S}_{\bf j}^{+}\right)$, with
$\hat{S}_{\bf i}^{z}\\!=\\!(\hat{n}_{{\bf i}\uparrow}\\!-\\!\hat{n}_{{\bf
i}\downarrow})/2$ and $\hat{S}_{\bf i}^{+}\\!=\\!(\hat{S}_{\bf
i}^{-})^{\dagger}\\!=\\!\tilde{c}^{\dagger}_{{\bf i}\uparrow}\tilde{c}_{{\bf
i}\downarrow}$. The constrained fermionic operators acting in the Hilbert
space without double occupancy are defined as $\tilde{c}^{{\dagger}}_{{\bf
i}\sigma}=\hat{c}_{{\bf i}\sigma}^{{\dagger}}(1-\hat{n}_{{\bf
i}\bar{\sigma}})$. Although sometimes negligible, there is a 3-site term at
the same lowest order ($\sim t^{2}/U$) that contributes to the dopant’s motion
Stephan and Horsch (1992); Jefferson et al. (1992); Bała et al. (1995); Spałek
(1988); von Szczepanski et al. (1990); Eskes and Oleś (1994); Eskes et al.
(1994); Belinicher and Chernyshev (1994); Belinicher et al. (1994); Belinicher
et al. (1996a, b); Psaltakis (1992); Eskes and Eder (1996). This defines the
$t\small{-}J\small{-}3s$ model with the Hamiltonian given by
$\mathcal{H}_{t\\!-\\!J\\!-\\!3s}=\mathcal{H}_{t\\!-\\!J}+\mathcal{H}_{3s}$:

$\mathcal{H}_{3s}=-\frac{J}{4}\\!\\!\sum_{\langle{\bf i},{\bf
j}\rangle,\langle{\bf i},{\bf j}^{\prime}\rangle\atop{\bf j}\neq{\bf
j}^{\prime},\sigma}\\!\\!\left(\tilde{c}^{\dagger}_{{\bf
j}^{\prime}\sigma}\tilde{n}_{{\bf i}\bar{\sigma}}\tilde{c}_{{\bf
j}\sigma}\\!+\\!\tilde{c}^{\dagger}_{{\bf
j}^{\prime}\sigma}\tilde{c}^{\dagger}_{{\bf i}\bar{\sigma}}\tilde{c}_{{\bf
i}\sigma}\tilde{c}_{{\bf j}\bar{\sigma}}\right).$ (3)

We use both the Hubbard and $t\small{-}J\small{-}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\small{-}J\small{-}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 (${\bf r}_{0}=\bf 0$), controlling the
mobility of the single hole. The Hamiltonian becomes:

$\mathcal{H}=\mathcal{H}_{0}+V\sum_{\sigma}\hat{n}_{{\bf 0}\sigma}\,$ (4)

where $\mathcal{H}_{0}$ is the Hubbard or $t\small{-}J\small{-}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 Koepsell et al. (2019). Hence, the high-order
correlations analyzed in this paper and Ref. Bohrdt et al., 2020 are directly
accessible to these experiments. Note that the strength of spin-exchange and
3-site terms in the $t\small{-}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|\leq U$. The difference between Hubbard and $t\small{-}J\small{-}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$\times$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

$\left\langle\hat{O}\right\rangle={\mathrm{Tr}}\left[\frac{e^{-{\mathcal{H}}/T}}{\mathcal{Z}}\hat{O}\right]\approx\sum_{n<n_{\rm
max}}\frac{e^{-E_{n}/T}}{\mathcal{Z}}\langle n|\hat{O}|n\rangle\,,$ (5)

in which $\mathcal{Z}$ is the partition function. The $n_{\rm max}$ sets the
numerical truncation of excited states, which satisfies $E_{n_{\rm
max}}-E_{0}\gg T$ for all temperatures considered in our work. In the results
presented in Secs. III and IV, we keep $n_{\rm max}\sim 650$ states for the
Hubbard model, giving reliable results up to $T\sim 0.4t$, whereas $n_{\rm
max}\sim 13,000$ states for the $t\small{-}J\small{-}3s$ model, giving
reliable results up to $T\sim t$.

In this paper, we employ the parallel Arnoldi and Paradeisos algorithm to
determine the equilibrium ground-state wavefunction and expectation
values.Lehoucq et al. (1998); Jia et al. (2018) 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\times 12$ cluster with cylindrical boundary
condition, to investigate possible artifacts caused by the finite system size.

##  III The Third-Order Spin-Hole Correlations

Figure 1:  Schematic showing the third-order spin-hole correlators for (a)
nearest spins ($|\mathbf{d}|=1$) and (b) next-nearest spins
($|\mathbf{d}|=\sqrt{2}$) .

To describe the spin correlations with respect to a doped hole, we consider
the following third-order correlation function:

$B(\textbf{r},\textbf{r}^{\prime},\mathbf{d})=4\left\langle\hat{n}^{h}_{\textbf{r}}\hat{S}^{z}_{\textbf{r}^{\prime}}\hat{S}^{z}_{\textbf{r}^{\prime}+\mathbf{d}}\right\rangle/\langle\hat{n}_{\textbf{r}}^{h}\rangle$
(6)

where the hole density operator
$\hat{n}_{\textbf{r}}^{h}=(1-\hat{n}_{\textbf{r}\uparrow})(1-\hat{n}_{\textbf{r}\downarrow})$.
In this paper, we focus on the $|\mathbf{d}|=1$ (neighboring spins) and
$|\mathbf{d}|=\sqrt{2}$ (diagonal spins), allowing the relative distance
between the hole and the spin “bond” $|\textbf{r}-\textbf{r}^{\prime}|$ to
vary (see Fig. 1).

Figure 2 gives an overview of the spatial $(\textbf{r}^{\prime}-\textbf{r})$
distribution of the $B(\textbf{r},\textbf{r}^{\prime},\mathbf{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 nearest-neighbor ($|\mathbf{d}|=1$)
correlations are always negative, while the diagonal ($|\mathbf{d}|=\sqrt{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(\textbf{r},\textbf{r}^{\prime},\mathbf{d})$ is weakened for shorter
distances $|\textbf{r}-\textbf{r}^{\prime}|$. This spatial distribution can be
understood in terms of the spin polaron Grusdt et al. (2019); Blomquist and
Carlström (2019), 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(\textbf{r},\textbf{r}^{\prime},\mathbf{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 $|\textbf{r}^{\prime}-\textbf{r}|$’s
asymptotically approach those of an undoped AFM system.

Figure 2:  (a) The spatial distribution of
$B(\textbf{r},\textbf{r}^{\prime},\mathbf{d})$ for $|\mathbf{d}|=1$ and
$|\mathbf{d}|=\sqrt{2}$ in a $n=1/16$ doped Hubbard model at $T=0$, with
pinning potential $V=0$, $V=2.5t$ and $V=5t$, respectively. The numbers in the
bottom-right corner of each panel mark the values of the correlators mentioned
in the main text. (b) The $V$-dependence of the differential third-order
correlators $\Delta B_{\rm nn}$ for $T=0.1t$ (blue), 0.2$t$ (green), 0.3$t$
(orange), and 0.4$t$ (red). The dots and solid lines denote the results
obtained for the Hubbard and $t\small{-}J\small{-}3s$ model, respectively. The
inset in (b) explains the definition of the differential correlators.  Figure
3:  (a) The spatial distribution of connected parts
$B^{c}(\textbf{r},\textbf{r}^{\prime},\mathbf{d})$ in a $n=1/16$ doped Hubbard
model at $T=0$, with pinning potential $V=0$, $V=2.5t$ and $V=5t$,
respectively. (b) The $V$-dependence of the connected third-order correlators
$B^{c}(\textbf{r},\textbf{r}^{\prime},\mathbf{d})$ for $T=0.1t$ (blue), 0.2$t$
(green), 0.3$t$ (orange), and 0.4$t$ (red). The dots and solid lines denote
the results obtained for the Hubbard and $t\small{-}J\small{-}3s$ model,
respectively. The inset in (b) explains the definition of the differential
correlators. (c) The temperature and pinning potential dependence of the
genuine third-order correlator
$B^{c}(\textbf{r},\textbf{r}+\hat{x},\mathbf{d})$ with
$|\mathbf{d}|=\sqrt{2}$.

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(\textbf{r},\textbf{r}^{\prime},\mathbf{d})$ between the hole (at r) and the
spins at $\textbf{r}^{\prime}$ and $\textbf{r}^{\prime}+\mathbf{d}$. As shown
in the right two panels of Fig. 2(a) and Fig. 2(b),
$B(\textbf{r},\textbf{r}^{\prime},\mathbf{d})$ becomes stronger for shorter
distances $|\textbf{r}-\textbf{r}^{\prime}|$ (i.e., for $|\mathbf{d}|=1$ it
becomes more negative, while for $|\mathbf{d}|=\sqrt{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 Koepsell et al. (2019). Indirectly, this phenomenon is also
consistent with the increased local moments measured in NMR and $\mu$SR
experiments in Zn-substituted cuprates Mahajan et al. (1994); Alloul et al.
(1991); Xiao et al. (1990); Mendels et al. (1994). 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) Essler et al. (1991); Seabra et al. (2014). 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 third-order correlators, we
define the difference between the nearest-neighbor correlators as [see the
illustration in the inset of Fig. 2(b)]

$\Delta B_{\rm
nn}=B(\textbf{r},\textbf{r}+2\hat{x},\hat{y})-B(\textbf{r},\textbf{r}+\hat{x},\hat{y}).$
(7)

Through the difference, the uniform (positive or negative) background in the
“bare” correlator $B(\textbf{r},\textbf{r}^{\prime},\mathbf{d})$ is removed.
Thus, we refer to the $\Delta B_{\rm nn}$ as differential third-order
correlators. As shown in Fig. 2(b), $\Delta B_{\rm nn}$ is always negative for
$V=0$, indicating the screening effect of the spin polaron; the $\Delta B_{\rm
nn}$ turns to positive for a finite $V$, indicating the onset of anti-
screeningdeg . As such, the crossover between screening and anti-screening
natures is clearly reflected by the differential third-order correlators,
while the original $B(\textbf{r},\textbf{r}^{\prime},\hat{y})$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.,
$\langle\hat{S}^{z}_{\textbf{r}}\hat{S}^{z}_{\textbf{r}+\hat{y}}\rangle$),
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 many-body physics. Such a genuine third-order
correlator can be more intuitively defined as the connected part of
$B(\textbf{r},\textbf{r}^{\prime},\mathbf{d})$

$B^{c}(\textbf{r},\textbf{r}^{\prime},\mathbf{d})=B(\textbf{r},\textbf{r}^{\prime},\mathbf{d})-\sum_{\textbf{r}^{\prime\prime}}\left\langle\hat{S}^{z}_{\textbf{r}^{\prime\prime}}\hat{S}^{z}_{\textbf{r}^{\prime\prime}+\mathbf{d}}\right\rangle/N,$
(8)

where $N$ denotes the number of lattice sites. As such, the
$B^{c}(\textbf{r},\textbf{r}^{\prime},\mathbf{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}(\textbf{r},\textbf{r}^{\prime},\mathbf{d})$’s for small
$|\textbf{r}-\textbf{r}^{\prime}|$ 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}(\textbf{r},\textbf{r}^{\prime},\mathbf{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\small{-}J\small{-}3s$ model. As shown in Figs. 3(b) and (c), the $\Delta
B_{\rm 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\small{-}J\small{-}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 $\mathcal{O}(t^{2}/U)$. On the large $V$ side,
$B^{c}(\textbf{r},\textbf{r}^{\prime},\mathbf{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_{\rm
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

Figure 4:  Schematic showing the fifth-order ring-spin correlators for (a)
$\mathbf{d}=(0,0)$ and (b) $\mathbf{d}=(1,1)$.  Figure 5:  (a) The
distribution of $C_{\lozenge}(\textbf{r},\textbf{d})$ for different distances
d, r in the center and $V=0$, calculated for a 6$\times$12 $t\small{-}J$ model
at zero temperature by DMRG. The weak asymmetry in the $x$-direction is caused
by the non-central location of the hole at $x=6$, while the system is
translational invariant along the $y$-direction due to the cylindrical
boundary condition. (b) The distribution of
$C_{\lozenge}(\textbf{r}=\textbf{r}_{0},\textbf{d})$ for different distances
d, strengths $V$ of the pinning potential in the center at $\textbf{r}_{0}$,
and temperatures $T$, calculated for the Hubbard model in a $4\times 4$
cluster. The central white circle denotes the position of the dopant hole
$\textbf{r}_{0}$ on the site of the pinning potential.

The third-order correlators $B(\textbf{r},\textbf{r}^{\prime},\mathbf{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 AFMGrusdt et al.
(2018, 2019). 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) ring-spin correlators with respect to a doped hole, which is
introduced recently in Ref. Bohrdt et al., 2020 as

$C_{\lozenge}(\textbf{r},\textbf{d})=2^{4}\left\langle\hat{n}^{\rm
h}_{\textbf{r}}\hat{S}^{z}_{\textbf{r}+\textbf{d}+\hat{x}}\hat{S}^{z}_{\textbf{r}+\textbf{d}+\hat{y}}\hat{S}^{z}_{\textbf{r}+\textbf{d}-\hat{x}}\hat{S}^{z}_{\textbf{r}+\textbf{d}-\hat{y}}\right\rangle/\langle\hat{n}_{\bf
r}^{h}\rangle.$ (9)

As illustrated in Fig. 4, this correlator reflects the ring-spin correlations
in a co-moving frame of the doped hole. By including non-zero displacements
$\bf d$ (between the ring-center and the hole), Eq. (9) generalizes the five-
point correlators introduced in Ref. Bohrdt et al., 2020.

Figure 6:  (a,b) The bare on-site ring-spin correlator
$C_{\lozenge}(\textbf{r}_{0},\textbf{d}=\mathbf{0})$ as a function of pinning
potential $V$ and temperature $T$, calculated using the (a)
$t\small{-}J\small{-}3s$ model and (b) Hubbard model. (c,d) Same as (a,b) but
for the connected part of the on-site ring-spin correlator
$C^{c}_{\lozenge}(\textbf{r}_{0},\textbf{d}=\mathbf{0})$ as a function of
pinning potential $V$ and temperature $T$.

We first investigate the dependence of $C_{\lozenge}(\textbf{r},\textbf{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_{\lozenge}(\textbf{r},\textbf{d})$ is only a function of d. As shown in
Fig. 5 (a) and the top-left panel of (b), $C_{\lozenge}$ is positive for any
$|\textbf{d}|>1$, but becomes sizable and negative for
$\textbf{d}=\mathbf{0}$. This behavior is reminiscent of a $\mathbb{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_{\lozenge}$ correlator, spins which are part of the
string originate from the opposite sub-lattice and contribute a negative
signBohrdt et al. (2020). Since $C_{\lozenge}(\textbf{r},\mathbf{0})$ involves
only one site belonging to the string, it is expected to be negative; whereas
$C_{\lozenge}(\textbf{r},{\bf d\neq\mathbf{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\small{-}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$\times$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 $\textbf{r}_{0}$. As shown in the upper panel of Fig. 5(b), the on-site
($\bf d=\mathbf{0}$) correlator $C_{\lozenge}(\textbf{r}_{0},\mathbf{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_{\lozenge}(\textbf{r},\mathbf{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_{\lozenge}(\textbf{r}_{0},\mathbf{0})$ remains negative until
$V\sim 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_{\lozenge}(\textbf{r}_{0},\mathbf{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_{\lozenge}(\textbf{r}_{0},\textbf{d})$s decrease with increasing
temperature, which occurs dramatically near $T\sim J$ and indicates the
thermal melting of the surrounding antiferromagnet. Given the presence of four
spins in the fifth-order correlator, its temperature dependence is more
evident than the third-order one presented in Fig. 3(c).

As such, the temperature dependence of this fifth-order 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. Bohrdt et al., 2020:

$C^{c}_{\lozenge}(\textbf{r},\textbf{d})=C_{\lozenge}(\textbf{r},\textbf{d})-2^{4}\biggl{[}\sum_{\mathbf{l}\notin(\mathbf{i},\mathbf{j},\mathbf{k})}\frac{\langle\hat{n}^{\rm
h}_{\mathbf{r}}\hat{S}^{z}_{\textbf{r}+\textbf{d}+\mathbf{i}}\hat{S}^{z}_{\textbf{r}+\textbf{d}+\mathbf{j}}\hat{S}^{z}_{\textbf{r}+\textbf{d}+\mathbf{k}}\rangle_{c}}{\langle\hat{n}^{\rm
h}_{\mathbf{r}}\rangle}\frac{\langle\hat{n}^{\rm
h}_{\textbf{r}}\hat{S}^{z}_{\textbf{r}+\textbf{d}+\mathbf{l}}\rangle}{{\langle\hat{n}^{\rm
h}_{\mathbf{r}}\rangle}}+\sum_{(\mathbf{i},\mathbf{j})\notin(\mathbf{k},\mathbf{l})}\frac{\langle\hat{n}^{\rm
h}_{\textbf{r}}\hat{S}^{z}_{\textbf{r}+\textbf{d}+\mathbf{i}}\hat{S}^{z}_{\textbf{r}+\textbf{d}+\mathbf{j}}\rangle_{c}}{\langle\hat{n}^{\rm
h}_{\mathbf{r}}\rangle}\frac{\langle\hat{n}^{\rm
h}_{\textbf{r}}\hat{S}^{z}_{\textbf{r}+\textbf{d}+\mathbf{k}}\hat{S}^{z}_{\textbf{r}+\textbf{d}+\mathbf{l}}\rangle_{c}}{\langle\hat{n}^{\rm
h}_{\mathbf{r}}\rangle}\\\
+\sum_{\mathbf{i}\neq\mathbf{j}\notin(\mathbf{k},\mathbf{l})}\frac{\langle\hat{n}^{\rm
h}_{\textbf{r}}\hat{S}^{z}_{\textbf{r}+\textbf{d}+\mathbf{i}}\rangle}{\langle\hat{n}^{\rm
h}_{\mathbf{r}}\rangle}\frac{\langle\hat{n}^{\rm
h}_{\textbf{r}}\hat{S}^{z}_{\textbf{r}+\textbf{d}+\mathbf{j}}\rangle}{\langle\hat{n}^{\rm
h}_{\mathbf{r}}\rangle}\frac{\langle\hat{n}^{\rm
h}_{\textbf{r}}\hat{S}^{z}_{\textbf{r}+\textbf{d}+\mathbf{k}}\hat{S}^{z}_{\textbf{r}+\textbf{d}+\mathbf{l}}\rangle_{c}}{\langle\hat{n}^{\rm
h}_{\mathbf{r}}\rangle}+\frac{\langle\hat{n}^{\rm
h}_{\textbf{r}}\hat{S}^{z}_{\textbf{r}+\textbf{d}+\hat{x}}\rangle}{\langle\hat{n}^{\rm
h}_{\mathbf{r}}\rangle}\frac{\langle\hat{n}^{\rm
h}_{\textbf{r}}\hat{S}^{z}_{\textbf{r}+\textbf{d}+\hat{y}}\rangle}{\langle\hat{n}^{\rm
h}_{\mathbf{r}}\rangle}\frac{\langle\hat{n}^{\rm
h}_{\textbf{r}}\hat{S}^{z}_{\textbf{r}+\textbf{d}-\hat{x}}\rangle}{\langle\hat{n}^{\rm
h}_{\mathbf{r}}\rangle}\frac{\langle\hat{n}^{\rm
h}_{\textbf{r}}\hat{S}^{z}_{\textbf{r}+\textbf{d}-\hat{y}}\rangle}{\langle\hat{n}^{\rm
h}_{\mathbf{r}}\rangle}\biggr{]},$ (10)

where the $\langle\cdots\rangle_{c}$ denotes the connected correlator with
lower-order terms subtracted.

Figures 6(c,d) show the dependence of the connected correlator
$C^{c}_{\lozenge}(\textbf{r},\mathbf{0})$ on the temperature and pinning
potential. We find that this correlator is sizably negative even for large
pinning potentials. Although the characteristic temperature, above which
$|C^{c}_{\lozenge}|$ drops dramatically, decreases slightly for the systems
with pinned holes, the absolute value of $C^{c}_{\lozenge}$ 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 lower-order correlators.

In the remainder of this section, we elucidate the physical origin and
properties of the genuine fifth-order ring-spin 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
$\langle\hat{S}^{z}\rangle=0$. Alternatively, in the former case, we could
consider an ensemble with definite $S^{z}=-1/2$, which gives the same results:
the higher-order correlators $C_{\lozenge}$ and $C_{\lozenge}^{c}$ each
involve an even number of spin operators $\hat{S}^{z}$. Moreover, the
underlying spin or Hubbard Hamiltonian is invariant under a global spin-flip
($S^{z}_{j}\to-S^{z}_{j}$). Hence, the correlators in the sectors with
definite $S^{z}=\pm 1/2$ are equivalent,

$C_{\lozenge}|_{S^{z}=1/2}=C_{\lozenge}|_{S^{z}=-1/2},\quad
C_{\lozenge}^{c}|_{S^{z}=1/2}=C_{\lozenge}^{c}|_{S^{z}=-1/2}.$ (11)

With the same argument, in the spin-balanced ensemble,
$\langle\hat{S}^{z}\rangle=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 solid-state system, both spin ensembles can
be addressed experimentally with ultracold atoms and with full spin and charge
resolutionBoll et al. (2016); Koepsell et al. (2020b). This ensemble can be
achieved by post-selecting experimental data with specific $S^{z}$ values.

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\to 0$. In the
resulting spin-balanced ensemble $\langle\hat{S}^{z}\rangle=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_{\lozenge}^{c}|_{\langle S^{z}\rangle=0}$
as a function of $t/J$ (see Fig. 1 in Bohrdt et al., 2020) 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
Grusdt et al. (2019).

Figure 7:  (a) Genuine fifth-order correlator $C^{c}_{\lozenge}$ for different
spin-ensembles, $\langle\hat{S}^{z}\rangle=0$ and $S^{z}=\pm 1/2$, at zero
temperature and without pinning potential. The results are obtained by the
DMRG simulations for a $t\small{-}J$ Hamiltonian on a $6\times 12$ cylinder.
(b) The bare fifth-order correlator $C_{\lozenge}$ and (c) the genuine
correlator $C^{c}_{\lozenge}$ 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\small{-}J\small{-}3s$ Hamiltonian.

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_{\lozenge}$ 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}_{\lozenge}$ 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}_{\lozenge}$ decreases in magnitude in the 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,
$\langle\hat{M}^{z}\rangle=\sum_{\textbf{r}}(-1)^{r_{x}+r_{y}}\langle\hat{S}^{z}_{\textbf{r}}\rangle\neq
0$, in the lattice. In a finite-size system, like the one we study, the latter
requires spin-imbalance 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 terms – e.g.
$\langle\hat{S}^{z}_{\textbf{r}}\rangle$ – 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\gg t$) and take into account only Ising couplings $\sim
J^{z}\hat{S}^{z}_{\mathbf{i}}\hat{S}^{z}_{\mathbf{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 $|{\rm N},+\rangle$ 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 $|{\rm
N},-\rangle$, obtained from $|{\rm N},+\rangle$ by flipping all spins. For
both Néel states, the bare ring-spin correlator
$C_{\lozenge}(\textbf{r},\mathbf{0})|_{S^{z}=1/2}=C_{\lozenge}(\textbf{r},\mathbf{0})|_{S^{z}=-1/2}$
is equal since all four spins surrounding the hole are aligned. Hence,
$C_{\lozenge}(\textbf{r},\mathbf{0})|_{\langle\hat{S}^{z}\rangle=0}=C_{\lozenge}(\textbf{r},\mathbf{0})|_{S^{z}=\pm
1/2}$ also takes the same value in the spin-balanced ensemble, defined by the
even statistical mixture $\\{|{\rm N},+\rangle,|{\rm N},-\rangle\\}$. Now we
turn to the connected correlators. For each of the two Néel states $|{\rm
N},\pm\rangle$ taken individually, $C_{\lozenge}^{c}|_{S^{z}=\pm 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 $\hat{S}^{z}$ operators. In contrast, in the spin-
balanced ensemble where $\langle\hat{S}^{z}\rangle=0$, only lower-order terms
with an even number of $\hat{S}^{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}_{\lozenge}(\textbf{r},\mathbf{0})|_{\langle\hat{S}^{z}\rangle=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\small{-}J\small{-}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 fifth-
order ring-spin correlators reflect the underlying spin polaron 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 spin-polaron, 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 high-order 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) fifth-order
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 solid-state 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.

## Acknowledgments

Y.W. acknowledges support from the National Science Foundation (NSF) award
DMR-2038011. A.B. and F.G. acknowledge support from Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s
Excellence Strategy – EXC-2111 – 390814868. J.K. acknowledges funding from
Hector Fellow Academy and support from the Max Planck Harvard Research Center
for Quantum Optics (MPHQ). E.D. acknowledges support from Harvard-MIT CUA, ARO
grant number W911NF-20-1-0163, the NSF award OAC-193471, and Harvard Quantum
Initiative. This research used resources of the National Energy Research
Scientific Computing Center (NERSC), a U.S. Department of Energy Office of
Science User Facility operated under Contract No. DE-AC02- 05CH11231.

Figure 8:  (a) Dependence of the differential third-order correlators $\Delta
B_{\rm nn}$ with bond lengths $|\mathbf{d}|=1$ as a function of the trapping
potential $V$ and temperature $T$, obtained for the $t\small{-}J\small{-}3s$
model (top) and the Hubbard model (bottom). (b) The same as (a) but for
$\Delta B_{\rm diag}$ with $|\mathbf{d}|=\sqrt{2}$. The insets in both (a) and
(b) explain the definition of the differential correlators.

##  Appendix A Details about the Differential Third-Order Correlator

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 $\Delta B_{\rm nn}$ in Fig. 8, following the same
manner of Fig. 3(c) in the main text.

In addition to the $\Delta B_{\rm nn}$ defined in Eq. 7, we further quantify
the difference between the diagonal (hole-spin) correlators as

$\Delta B_{\rm
diag}=B(\textbf{r},\textbf{r}+\hat{x},\hat{y}-\hat{x})-B(\textbf{r},\textbf{r}+2\hat{x},\hat{y}-\hat{x})\,.$
(12)

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 $\Delta B_{\rm nn}$ [Fig. 8(a)] and the
$\Delta B_{\rm diag}$ [Fig. 8(b)]: although the screening effect has been
dramatically suppressed for large pinning potentials following the similar
crossover, $\Delta B_{\rm 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 $\Delta B_{\rm
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}(\textbf{r},\textbf{r}^{\prime},\mathbf{d})$ over the differential
correlator $\Delta B_{\rm 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.

## References

  * Dagotto (1994) E. Dagotto, Rev. Mod. Phys. 66, 763 (1994).
  * Anderson (1987) P. W. Anderson, science 235, 1196 (1987).
  * Zhang and Rice (1988) F. Zhang and T. Rice, Phys. Rev. B 37, 3759 (1988).
  * Lee et al. (2006) P. A. Lee, N. Nagaosa, and X.-G. Wen, Rev. Mod. Phys. 78, 17 (2006).
  * Singh and Goswami (2002) A. Singh and P. Goswami, Phys. Rev. B 66, 092402 (2002).
  * Zheng et al. (2017) B.-X. Zheng, C.-M. Chung, P. Corboz, G. Ehlers, M.-P. Qin, R. M. Noack, H. Shi, S. R. White, S. Zhang, and G. K.-L. Chan, Science 358, 1155 (2017).
  * Huang et al. (2017) E. W. Huang, C. B. Mendl, S. Liu, S. Johnston, H.-C. Jiang, B. Moritz, and T. P. Devereaux, Science 358, 1161 (2017).
  * Ponsioen et al. (2019) B. Ponsioen, S. S. Chung, and P. Corboz, Phys. Rev. B 100, 195141 (2019).
  * Kokalj (2017) J. Kokalj, Phys. Rev. B 95, 041110 (2017).
  * Huang et al. (2019) E. W. Huang, R. Sheppard, B. Moritz, and T. P. Devereaux, Science 366, 987 (2019).
  * Cha et al. (2020) P. Cha, A. A. Patel, E. Gull, and E.-A. Kim, Phys. Rev. Research 2, 033434 (2020).
  * Maier et al. (2005) T. A. Maier, M. Jarrell, T. Schulthess, P. Kent, and J. White, Phys. Rev. Lett. 95, 237001 (2005).
  * Zheng and Chan (2016) B.-X. Zheng and G. K.-L. Chan, Phys. Rev. B 93, 035126 (2016).
  * Ido et al. (2018) K. Ido, T. Ohgoe, and M. Imada, Phys. Rev. B 97, 045138 (2018).
  * Jiang and Devereaux (2019) H.-C. Jiang and T. P. Devereaux, Science 365, 1424 (2019).
  * Chung et al. (2002) C.-M. Chung, M. Qin, S. Zhang, U. Schollwöck, and S. R. White, Phys. Rev. B 102, 041106 (2002).
  * Ohta et al. (1991) Y. Ohta, T. Tohyama, and S. Maekawa, Phys. Rev. B 43, 2968 (1991).
  * Lanzara et al. (2001) A. Lanzara, P. Bogdanov, X. Zhou, S. Kellar, D. Feng, E. Lu, T. Yoshida, H. Eisaki, A. Fujimori, K. Kishio, et al., Nature 412, 510 (2001).
  * Shen et al. (2004) K. Shen, F. Ronning, D. Lu, W. Lee, N. Ingle, W. Meevasana, F. Baumberger, A. Damascelli, N. Armitage, L. Miller, et al., Phys. Rev. Lett. 93, 267002 (2004).
  * Mishchenko and Nagaosa (2004) A. Mishchenko and N. Nagaosa, Phys. Rev. Lett. 93, 036402 (2004).
  * Cuk et al. (2004) T. Cuk, F. Baumberger, D. Lu, N. Ingle, X. Zhou, H. Eisaki, N. Kaneko, Z. Hussain, T. Devereaux, N. Nagaosa, et al., Phys. Rev. Lett. 93, 117003 (2004).
  * Reznik et al. (2006) D. Reznik, L. Pintschovius, M. Ito, S. Iikubo, M. Sato, H. Goka, M. Fujita, K. Yamada, G. Gu, and J. Tranquada, Nature 440, 1170 (2006).
  * Johnston et al. (2012) S. Johnston, I. Vishik, W. Lee, F. Schmitt, S. Uchida, K. Fujita, S. Ishida, N. Nagaosa, Z. Shen, and T. Devereaux, Phys. Rev. Lett. 108, 166404 (2012).
  * Peng et al. (2017) Y. Peng, G. Dellea, M. Minola, M. Conni, A. Amorese, D. Di Castro, G. De Luca, K. Kummer, M. Salluzzo, X. Sun, et al., Nat. Phys. 13, 1201 (2017).
  * He et al. (2018) Y. He, M. Hashimoto, D. Song, S.-D. Chen, J. He, I. Vishik, B. Moritz, D.-H. Lee, N. Nagaosa, J. Zaanen, et al., Science 362, 62 (2018).
  * Demler et al. (2004) E. Demler, W. Hanke, and S.-C. Zhang, Rev. Mod. Phys. 76, 909 (2004).
  * Scalapino (2012) D. J. Scalapino, Rev. Mod. Phys. 84, 1383 (2012).
  * Maier et al. (2016) T. A. Maier, P. Staar, V. Mishra, U. Chatterjee, J. Campuzano, and D. J. Scalapino, Nature communications 7, 1 (2016).
  * LeBlanc et al. (2015) J. P. LeBlanc, A. E. Antipov, F. Becca, I. W. Bulik, G. K.-L. Chan, C.-M. Chung, Y. Deng, M. Ferrero, T. M. Henderson, C. A. Jiménez-Hoyos, et al., Phys. Rev. X 5, 041041 (2015).
  * Kane et al. (1989) C. L. Kane, P. A. Lee, and N. Read, Phys. Rev. B 39, 6880 (1989).
  * Sachdev (1989) S. Sachdev, Phys. Rev. B 39, 12232 (1989).
  * Liu and Manousakis (1991) Z. Liu and E. Manousakis, Phys. Rev. B 44, 2414 (1991).
  * Martinez and Horsch (1991) G. Martinez and P. Horsch, Phys. Rev. B 44, 317 (1991).
  * Stephan and Horsch (1991) W. Stephan and P. Horsch, Phys. Rev. Lett. 66, 2258 (1991).
  * Manousakis (1991) E. Manousakis, Rev. Mod. Phys. 63, 1 (1991).
  * Dagotto et al. (1992) E. Dagotto, A. Moreo, F. Ortolani, D. Poilblanc, and J. Riera, Phys. Rev. B 45, 10741 (1992).
  * Liu and Manousakis (1992) Z. Liu and E. Manousakis, Phys. Rev. B 45, 2425 (1992).
  * Stephan and Horsch (1992) W. Stephan and P. Horsch, Int. J Mod. Phys. B 6, 589 (1992).
  * Bała et al. (1995) J. Bała, A. Oleś, and J. Zaanen, Phys. Rev. B 52, 4597 (1995).
  * Brunner et al. (2000) M. Brunner, F. F. Assaad, and A. Muramatsu, Phys. Rev. B 62, 15480 (2000).
  * Scalapino et al. (1986) D. Scalapino, E. Loh Jr, and J. Hirsch, Phys. Rev. B 34, 8190 (1986).
  * Gros et al. (1987) C. Gros, R. Joynt, and T. Rice, Zeitschrift für Physik B Condensed Matter 68, 425 (1987).
  * Kotliar and Liu (1988) G. Kotliar and J. Liu, Phys. Rev. B 38, 5142 (1988).
  * Schrieffer et al. (1989) J. Schrieffer, X. Wen, and S. Zhang, Phys. Rev. B 39, 11663 (1989).
  * Scalapino (1995) D. J. Scalapino, Physics Reports 250, 329 (1995).
  * Tsuei and Kirtley (2000) C. Tsuei and J. Kirtley, Rev. Mod. Phys. 72, 969 (2000).
  * Graf et al. (2007) J. Graf, G.-H. Gweon, K. McElroy, S. Zhou, C. Jozwiak, E. Rotenberg, A. Bill, T. Sasagawa, H. Eisaki, S. Uchida, et al., Phys. Rev. Lett. 98, 067004 (2007).
  * He et al. (2019) J. He, C. R. Rotundu, M. S. Scheurer, Y. He, M. Hashimoto, K.-J. Xu, Y. Wang, E. W. Huang, T. Jia, S. Chen, et al., Proc. Natl. Acad. Sci. U. S. A. 116, 3449 (2019).
  * Chen et al. (2019) S.-D. Chen, M. Hashimoto, Y. He, D. Song, K.-J. Xu, J.-F. He, T. P. Devereaux, H. Eisaki, D.-H. Lu, J. Zaanen, et al., Science 366, 1099 (2019).
  * Le Tacon et al. (2011) M. Le Tacon, G. Ghiringhelli, J. Chaloupka, M. M. Sala, V. Hinkov, M. Haverkort, M. Minola, M. Bakr, K. Zhou, S. Blanco-Canosa, et al., Nat. Phys. 7, 725 (2011).
  * Dean et al. (2013a) M. Dean, G. Dellea, R. Springell, F. Yakhou-Harris, K. Kummer, N. Brookes, X. Liu, Y. Sun, J. Strle, T. Schmitt, et al., Nat. Mater. 12, 1019 (2013a).
  * Dean et al. (2013b) M. Dean, A. James, R. Springell, X. Liu, C. Monney, K. Zhou, R. Konik, J. Wen, Z. Xu, G. Gu, et al., Phys. Rev. Lett. 110, 147001 (2013b).
  * Lee et al. (2014) W. Lee, J. Lee, E. Nowadnick, S. Gerber, W. Tabis, S. Huang, V. Strocov, E. Motoyama, G. Yu, B. Moritz, et al., Nat. Phys. (2014).
  * Ishii et al. (2014) K. Ishii, M. Fujita, T. Sasaki, M. Minola, G. Dellea, C. Mazzoli, K. Kummer, G. Ghiringhelli, L. Braicovich, T. Tohyama, et al., Nat. Commun. 5 (2014).
  * Hasselmann et al. (2004) N. Hasselmann, A. C. Neto, and C. M. Smith, Phys. Rev. B 69, 014424 (2004).
  * Mahajan et al. (1994) A. Mahajan, H. Alloul, G. Collin, and J. Marucco, Phys. Rev. Lett. 72, 3100 (1994).
  * Alloul et al. (1991) H. Alloul, P. Mendels, H. Casalta, J. Marucco, and J. Arabski, Phys. Rev. Lett. 67, 3140 (1991).
  * Xiao et al. (1990) G. Xiao, M. Z. Cieplak, J. Xiao, and C. Chien, Phys. Rev. B 42, 8752 (1990).
  * Mendels et al. (1994) P. Mendels, H. Alloul, J. Brewer, G. Morris, T. Duty, S. Johnston, E. Ansaldo, G. Collin, J. Marucco, C. Niedermayer, et al., Phys. Rev. B 49, 10035 (1994).
  * Maeda et al. (1990) A. Maeda, T. Yabe, S. Takebayashi, M. Hase, and K. Uchinokura, Phys. Rev. B 41, 4112 (1990).
  * Fukuzumi et al. (1996) Y. Fukuzumi, K. Mizuhashi, K. Takenaka, and S. Uchida, Phys. Rev. Lett. 76, 684 (1996).
  * Nachumi et al. (1996) B. Nachumi, A. Keren, K. Kojima, M. Larkin, G. Luke, J. Merrin, O. Tchernyshöv, Y. Uemura, N. Ichikawa, M. Goto, et al., Phys. Rev. Lett. 77, 5421 (1996).
  * Basov et al. (1998) D. Basov, B. Dabrowski, and T. Timusk, Phys. Rev. Lett. 81, 2132 (1998).
  * Lee (1993) P. A. Lee, Phys. Rev. Lett. 71, 1887 (1993).
  * Sun and Maki (1995) Y. Sun and K. Maki, Phys. Rev. B 51, 6059 (1995).
  * Salkola et al. (1996) M. Salkola, A. Balatsky, and D. Scalapino, Phys. Rev. Lett. 77, 1841 (1996).
  * Hirschfeld and Putikka (1996) P. Hirschfeld and W. Putikka, Phys. Rev. Lett. 77, 3909 (1996).
  * Franz et al. (1996) M. Franz, C. Kallin, and A. Berlinsky, Phys. Rev. B 54, R6897 (1996).
  * Tsuchiura et al. (1999) H. Tsuchiura, Y. Tanaka, M. Ogata, and S. Kashiwaya, Journal of the Physical Society of Japan 68, 2510 (1999).
  * Kivelson et al. (1998) S. A. Kivelson, E. Fradkin, and V. J. Emery, Nature 393, 550 (1998).
  * Emery et al. (1999) V. Emery, S. Kivelson, and J. Tranquada, Proc. Nat. Acad. Sci. 96, 8814 (1999).
  * Kivelson et al. (2003) S. A. Kivelson, I. P. Bindloss, E. Fradkin, V. Oganesyan, J. Tranquada, A. Kapitulnik, and C. Howald, Rev. Mod. Phys. 75, 1201 (2003).
  * Gross and Bloch (2017) C. Gross and I. Bloch, Science 357, 995 (2017).
  * Bakr et al. (2009) W. S. Bakr, J. I. Gillen, A. Peng, S. Fölling, and M. Greiner, Nature 462, 74 (2009).
  * Sherson et al. (2010) J. F. Sherson, C. Weitenberg, M. Endres, M. Cheneau, I. Bloch, and S. Kuhr, Nature 467, 68 (2010).
  * Parsons et al. (2015) M. F. Parsons, F. Huber, A. Mazurenko, C. S. Chiu, W. Setiawan, K. Wooley-Brown, S. Blatt, and M. Greiner, Phys. Rev. Lett. 114, 213002 (2015).
  * Cheuk et al. (2015) L. W. Cheuk, M. A. Nichols, M. Okan, T. Gersdorf, V. V. Ramasesh, W. S. Bakr, T. Lompe, and M. W. Zwierlein, Phys. Rev. Lett. 114, 193001 (2015).
  * Omran et al. (2015) A. Omran, M. Boll, T. A. Hilker, K. Kleinlein, G. Salomon, I. Bloch, and C. Gross, Phys. Rev. Lett. 115, 263001 (2015).
  * Edge et al. (2015) G. J. Edge, R. Anderson, D. Jervis, D. C. McKay, R. Day, S. Trotzky, and J. H. Thywissen, Phys. Rev. A 92, 063406 (2015).
  * Haller et al. (2015) E. Haller, J. Hudson, A. Kelly, D. A. Cotta, B. Peaudecerf, G. D. Bruce, and S. Kuhr, Nature Physics 11, 738 (2015).
  * Mazurenko et al. (2017) A. Mazurenko, C. S. Chiu, G. Ji, M. F. Parsons, M. Kanász-Nagy, R. Schmidt, F. Grusdt, E. Demler, D. Greif, and M. Greiner, Nature 545, 462 (2017).
  * Hilker et al. (2017) T. A. Hilker, G. Salomon, F. Grusdt, A. Omran, M. Boll, E. Demler, I. Bloch, and C. Gross, Science 357, 484 (2017).
  * Chiu et al. (2019) C. S. Chiu, G. Ji, A. Bohrdt, M. Xu, M. Knap, E. Demler, F. Grusdt, M. Greiner, and D. Greif, Science 365, 251 (2019).
  * Schweigler et al. (2017) T. Schweigler, V. Kasper, S. Erne, I. Mazets, B. Rauer, F. Cataldini, T. Langen, T. Gasenzer, J. Berges, and J. Schmiedmayer, Nature 545, 323 (2017).
  * Salomon et al. (2019) G. Salomon, J. Koepsell, J. Vijayan, T. A. Hilker, J. Nespolo, L. Pollet, I. Bloch, and C. Gross, Nature 565, 56 (2019).
  * Brown et al. (2019) P. T. Brown, D. Mitra, E. Guardado-Sanchez, R. Nourafkan, A. Reymbaut, C.-D. Hébert, S. Bergeron, A.-M. Tremblay, J. Kokalj, D. A. Huse, et al., Science 363, 379 (2019).
  * Nichols et al. (2019) M. A. Nichols, L. W. Cheuk, M. Okan, T. R. Hartke, E. Mendez, T. Senthil, E. Khatami, H. Zhang, and M. W. Zwierlein, Science 363, 383 (2019).
  * Koepsell et al. (2019) J. Koepsell, J. Vijayan, P. Sompet, F. Grusdt, T. A. Hilker, E. Demler, G. Salomon, I. Bloch, and C. Gross, Nature 572, 358 (2019).
  * Bohrdt et al. (2020) A. Bohrdt, Y. Wang, J. Koepsell, E. Demler, and F. Grusdt (2020).
  * Koepsell et al. (2020a) J. Koepsell, D. Bourgund, P. Sompet, S. Hirthe, A. Bohrdt, Y. Wang, F. Grusdt, E. Demler, G. Salomon, C. Gross, et al., arXiv:2009.04440 (2020a).
  * Hofstetter et al. (2002) W. Hofstetter, J. I. Cirac, P. Zoller, E. Demler, and M. Lukin, Phys. Rev. Lett. 89, 220407 (2002).
  * Ho et al. (2009) A. Ho, M. Cazalilla, and T. Giamarchi, Phys. Rev. A 79, 033620 (2009).
  * Jördens et al. (2010) R. Jördens, L. Tarruell, D. Greif, T. Uehlinger, N. Strohmaier, H. Moritz, T. Esslinger, L. De Leo, C. Kollath, A. Georges, et al., Phys. Rev. Lett. 104, 180401 (2010).
  * Esslinger (2010) T. Esslinger (2010).
  * Tarruell and Sanchez-Palencia (2018) L. Tarruell and L. Sanchez-Palencia, Comptes Rendus Physique 19, 365 (2018).
  * Chao et al. (1977) K. Chao, J. Spalek, and A. Oleś, J Phys C: Solid State Phys 10, L271 (1977).
  * Chao et al. (1978) K. Chao, J. Spałek, and A. Oleś, Phys. Rev. B 18, 3453 (1978).
  * Belinicher and Chernyshev (1994) V. Belinicher and A. Chernyshev, Phys. Rev. B 49, 9746 (1994).
  * Belinicher et al. (1994) V. Belinicher, A. Chernyshev, and L. Popovich, Lett. B 50, 13768 (1994).
  * Jefferson et al. (1992) J. H. Jefferson, H. Eskes, and L. F. Feiner, Phys. Rev. B 45, 7959 (1992).
  * Spałek (1988) J. Spałek, Phys. Rev. B 37, 533 (1988).
  * von Szczepanski et al. (1990) K. J. von Szczepanski, P. Horsch, W. Stephan, and M. Ziegler, Phys. Rev. B 41, 2017 (1990).
  * Eskes and Oleś (1994) H. Eskes and A. M. Oleś, Phys. Rev. Lett. 73, 1279 (1994).
  * Eskes et al. (1994) H. Eskes, A. M. Oleś, M. B. J. Meinders, and W. Stephan, Phys. Rev. B 50, 17980 (1994).
  * Belinicher et al. (1996a) V. Belinicher, A. Chernyshev, and V. Shubin, Phys. Rev. B 53, 335 (1996a).
  * Belinicher et al. (1996b) V. Belinicher, A. Chernyshev, and V. Shubin, Phys. Rev. B 54, 14914 (1996b).
  * Psaltakis (1992) G. C. Psaltakis, Phys. Rev. B 45, 539 (1992).
  * Eskes and Eder (1996) H. Eskes and R. Eder, Phys. Rev. B 54, R14226 (1996).
  * Lehoucq et al. (1998) R. B. Lehoucq, D. C. Sorensen, and C. Yang, _ARPACK Users’ Guide: Solution of Large-Scale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods_ (Siam, 1998).
  * Jia et al. (2018) C. Jia, Y. Wang, C. Mendl, B. Moritz, and T. Devereaux, Comput. Phys. Commun. 224, 81 (2018).
  * Grusdt et al. (2019) F. Grusdt, A. Bohrdt, and E. Demler, Phys. Rev. B 99, 224422 (2019).
  * Blomquist and Carlström (2019) E. Blomquist and J. Carlström, arXiv:1912.08825 (2019).
  * Essler et al. (1991) F. H. Essler, V. E. Korepin, and K. Schoutens, Phys. Rev. Lett. 67, 3848 (1991).
  * Seabra et al. (2014) L. Seabra, F. H. Essler, F. Pollmann, I. Schneider, and T. Veness, Phys. Rev. B 90, 245127 (2014).
  * (115) The fact that the correlation switches sign for any infinitesimal $V$ results from the extra hypercubic symmetry in a $4\times 4$ Hubbard model. This symmetry leads to a 6-fold degeneracy for the ground state without $V$, and is broken for any finite $V$. The critical $V$ should lie at a small finite value at the thermodynamic limit, similar to the presented results at finite temperature.
  * Grusdt et al. (2018) F. Grusdt, M. Kánasz-Nagy, A. Bohrdt, C. S. Chiu, G. Ji, M. Greiner, D. Greif, and E. Demler, Phys. Rev. X 8, 011046 (2018).
  * Boll et al. (2016) M. Boll, T. A. Hilker, G. Salomon, A. Omran, J. Nespolo, L. Pollet, I. Bloch, and C. Gross, Science 353, 1257 (2016).
  * Koepsell et al. (2020b) J. Koepsell, S. Hirthe, D. Bourgund, P. Sompet, J. Vijayan, G. Salomon, C. Gross, and I. Bloch, arXiv:2002.07577 (2020b).