# Unravelling the Nature of Spin Excitations Disentangled from Charge
Contributions in a Doped Cuprate Superconductor

Wenliang Zhang<EMAIL_ADDRESS>Photon Science Division, Swiss Light
Source, Paul Scherrer Institut, CH-5232 Villigen PSI, Switzerland  Cliò
Efthimia Agrapidis  Institute of Theoretical Physics, Faculty of Physics,
University of Warsaw, Pasteura 5, PL-02093 Warsaw, Poland  Yi Tseng  Photon
Science Division, Swiss Light Source, Paul Scherrer Institut, CH-5232 Villigen
PSI, Switzerland  Teguh Citra Asmara  Photon Science Division, Swiss Light
Source, Paul Scherrer Institut, CH-5232 Villigen PSI, Switzerland  Eugenio
Paris  Photon Science Division, Swiss Light Source, Paul Scherrer Institut,
CH-5232 Villigen PSI, Switzerland  Vladimir N. Strocov  Photon Science
Division, Swiss Light Source, Paul Scherrer Institut, CH-5232 Villigen PSI,
Switzerland  Enrico Giannini  Department of Quantum Matter Physics, University
of Geneva, Quai Ernest-Ansermet 24, CH-1211 Genève 4, Switzerland  Satoshi
Nishimoto  Department of Physics, Technical University Dresden, 01069 Dresden,
Germany  Institute for Theoretical Solid State Physics, IFW Dresden, 01171
Dresden, Germany  Krzysztof Wohlfeld<EMAIL_ADDRESS>Institute
of Theoretical Physics, Faculty of Physics, University of Warsaw, Pasteura 5,
PL-02093 Warsaw, Poland  Thorsten Schmitt<EMAIL_ADDRESS>Photon
Science Division, Swiss Light Source, Paul Scherrer Institut, CH-5232 Villigen
PSI, Switzerland

###### Abstract

The nature of the spin excitations in superconducting cuprates is a key
question toward a unified understanding of the cuprate physics from long-range
antiferromagnetism to superconductivity. The intense spin excitations up to
the over-doped regime revealed by resonant inelastic X-ray scattering bring
new insights as well as questions like how to understand their persistence or
their relation to the collective excitations in ordered magnets (magnons).
Here, we study the evolution of the spin excitations upon hole-doping the
superconducting cuprate Bi2Sr2CaCu2O8+δ by disentangling the spin from the
charge excitations in the experimental cross section. We compare our
experimental results against density matrix renormalization group calculations
for a $t$-$J$-like model on a square lattice. Our results unambiguously
confirm the persistence of the spin excitations, which are closely connected
to the persistence of short-range magnetic correlations up to high doping.
This suggests that the spin excitations in hole-doped cuprates are related to
magnons—albeit short-ranged.

INTRODUCTION
Starting from antiferromagnetic Mott insulators, the cuprate high-temperature
superconductors go through various quantum states with the charge carrier
doping as the tuning parameter and form a universal doping-temperature phase
diagram. For the hole-doped cuprates, superconductivity emerges in the
intermediate doping regime in close proximity to the notorious pseudogap state
and the recently established charge density wave state [1]. Disentangling the
physics behind these intertwined states is a major challenge for constructing
a complete theory of the superconducting cuprates. Fundamentally, the
competition between the exchange energy of the localized spins and the kinetic
energy of the doped holes is believed to dominate the basic physics, and the
two energy scales are naturally regulated by the amount of doped holes [2].
While the holes tend to delocalize and turn the system into a Fermi liquid at
high doping level, electron correlations and local antiferromagnetic
correlations survive with modest doping in the superconducting regime and
coexist with the well-defined quasiparticles [2, 3, 4, 5]. Exactly how these
spin and charge degrees of freedom act and interact throughout the doping-
temperature phase diagram is therefore a crucial question towards the
formulation of a definitive theory of superconductivity in the doped cuprates.

Experimentally, the dynamics and interactions of the spin and charge degrees
of freedom are studied by assessing the momentum and energy dependence of the
elementary excitations across the phase diagram. As suggested in numerous
papers over the last 11 years [6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18,
19, 20, 21, 22, 23, 24, 25, 26, 27], resonant inelastic X-ray scattering
(RIXS) observes intense spin excitations in a large momentum space area around
the Brillouin zone center, which persist well into the over-doped region.
Crucially the spin excitations in doped samples disperse along the ($\uppi$,
0) direction similarly to the magnons in the antiferromagnetic phase with the
damping increasing moderately, although they rapidly become overdamped with
doping along the ($\uppi$, $\uppi$) direction and show almost nondispersive
profiles [14, 15, 17, 20, 21, 22, 23, 26, 27]. This suggests that in some
particular area of the Brillouin zone, that is mostly in the ($\uppi$, 0)
direction, the magnetic excitations completely ‘ignore’ the existence of a
critical value of the doping $\delta$ at which Fermi liquid behavior takes
over the correlated magnetism and, moreover, that these excitations resemble
the well-known magnons in undoped cuprates (hence their name—paramagnons).
Naturally, such results are very much counterintuitive for they lead to an
apparent paradox related to the small changes of the spin excitations in the
doped cuprates, despite the rapid collapse of the long-range magnetic order
upon doping and the dominant Fermi-liquid nature in the over-doped regime.
This has sparked an intensive discussion on whether the observed magnetic
excitations should indeed be viewed as paramagnons—or rather incoherent
particle-hole excitations with a spin-flip [14, 15, 24, 28, 29, 30].

To justify the nature of the spin excitations in cuprates as well as the
reason of their persistence upon doping, it is necessary to precisely evaluate
the momentum and doping evolution of the intrinsic spin excitations and
comprehensively compare to theoretical calculations. This is a difficult task
due to the experimental difficulties in extracting $S(\mathbf{q},\omega)$ from
RIXS spectra and also due to the problems in reliably calculating
$S(\mathbf{q},\omega)$. One of the major experimental obstacles comes from the
mixing of spin and charge excitations in the RIXS spectra of doped cuprates
[31, 32]. With increasing doped holes, one would expect the charge excitations
to be stronger, which will worsen the ‘mixing’ problem. This strongly hampers
the correct assignment of the spectral profile to solely spin-flip containing
excitations, and casts doubts on whether RIXS indeed observes the persistence
of the intrinsic spin excitations.

Here we report a systematic study on the momentum and doping evolution of the
disentangled intrinsic spin and charge excitations in superconducting
Bi2Sr2CaCu2O8+δ (Bi2212). By applying the azimuthal-dependent analysis based
on the distinct scattering tensors [33], we show that the low-energy
excitations of the Cu $L_{3}$-edge RIXS spectra can be well described by spin-
flip (spin) and non-spin-flip (charge) components for all studied doping
levels and momenta, which allows us to extract the intrinsic spectral weights
of the two components. We find that the spin spectral weight only slightly
increases (decreases) with doping at intermediate momentum q along the
($\uppi$, $\uppi$) [($\uppi$, 0)] direction, which unequivocally confirms the
persistence of the spin excitations in doped cuprates. We then compare the
above experimental results to state-of-the-art density matrix renormalization
group (DMRG) calculations of $t$-$J$-like models. The detailed comparison
reveals the key characteristics of the spin excitations in the doped cuprates:
On one hand, we unravel the crucial role of the longer-range hoppings, as the
second- and third-nearest-neighbor hopping are needed to fully reproduce the
experimental spin spectrum; on the other hand, we show that solely the short-
range magnetic correlations are enough to qualitatively reproduce the
persistence of the intensity of the spin excitations upon doping. Our results
thus suggest that RIXS indeed observes the persistent spin excitations in
hole-doped cuprates and that their paramagnetic nature can be understood as
stemming from localised spins with short-range correlations.

Fig. 1: Experimental setup and the azimuthal dependence of non-spin-flip and
spin-flip RIXS responses. a Schematic temperature-doping phase diagram for
Bi2212. The purple, green and cyan areas represent the antiferromagnetic,
pseudogap, and superconducting states, respectively. The solid circles
indicate the locations of the three measured samples (UD, OD1, OD2) in the
phase diagram. b The scattering geometry of the azimuthal dependent RIXS
experiment on a wedged sample. c – f The azimuthal $\phi$ dependence of the
non-spin-flip and spin-flip RIXS responses in cuprates with a $40^{\circ}$
wedge angle ($\theta_{\textrm{w}}$). c and d are the ‘bare’ $\phi$ dependence
on infinite thin samples, while e and f give the intensity evolution on bulk
samples including the self-absorption effect at an energy loss of -0.25 eV.

RESULTS
Disentangling the spin and charge excitations in the RIXS spectra
We studied the spin and charge excitations of three different doped Bi2212
samples from under-doped to over-doped regime. The three samples are labeled
as UD ($T_{\textrm{c}}$ = 73 K), OD1 ($T_{\textrm{c}}$ = 88 K) and OD2
($T_{\textrm{c}}$ = 65 K), as shown in Fig 1a. To disentangle the overlapping
spin and charge excitations in the Cu $L_{3}$-edge RIXS spectra, we applied
the recently proposed azimuthal dependent method [33], which resolves the two
kinds of excitations based on their distinct scattering tensors. In this
method, a sample on a wedged sample holder is rotated to change the
orientation of the photon polarization in the sample space (see Fig. 1b and
Methods section), which gives rise to different rotation dependences according
to the scattering tensors [33]. Fig. 1c,d show the azimuthal dependences of
the non-spin flip and spin-flip excitations at $\sigma$ and $\pi$ incident
polarizations with a 40∘ wedge angle. In the RIXS experiment, self-absorption
effects need to be considered, which depends on the scattering geometry
(azimuthal angles in this experiment) as well as the polarization and energy
of the photons. Fig. 1e,f show the azimuthal dependence of spin and charge
response after including the self-absorption effect at an energy loss
($E_{\textrm{f}}-E_{\textrm{i}}$) of -0.25 eV (see Supplementary Note 1). The
spin and charge responses show clearly different azimuthal dependences, which
allows to disentangle them in the low-energy RIXS spectra.

Fig. 2: Azimuthal dependent RIXS spectra and the decomposition into the spin
and charge contributions. a – b Azimuthal dependence of the low-energy RIXS
spectra of the OD2 sample, measured with $\sigma$ and $\pi$ incident
polarizations, respectively, and at an in-plane momentum transfer of
$\mathbf{q}$ = (0.33, 0). c – f Constant-energy-loss cuts of the azimuthal
dependent RIXS intensity at 0 eV (c and d) and -0.25 eV (e and f). The solid
green and red lines represent the decomposed spin and charge components
$w_{\mathrm{s(c)}}{\cdot}A_{\mathrm{s(c)}}(\phi)$, respectively, and the solid
dark lines are the sums of the two components. g \- h The RIXS energy-loss
spectra compared to the spin and charge components at grazing incidence
($\phi$ = 0∘) with $\sigma$ polarization and grazing emission ($\phi$ = 180∘)
with $\pi$ polarization. The error bars represent the combined errors
including the statistic errors and the errors in determining the zero-energy-
loss positions (See Methods).

Fig. 2a and b show the RIXS intensity map of OD2 sample at $\mathbf{q}$ =
(0.33, 0) as a function of azimuthal angle $\phi$ and energy loss $E$ at
$\sigma$ and $\pi$ incident polarization, respectively. As the energies of
$d$-$d$ excitations are situated well above 1 eV [14, 15, 34], it is natural
to assume that the low-energy excitations are mainly composed of the single-
spin-flip and non-spin-flip excitations involving the 3$d_{x^{2}-y^{2}}$
orbital. Although the double spin flip with a net spin change of zero
(bimagnon with $\Delta$$S$ = 0) is also allowed in RIXS process in cuprates
[35, 36, 25], earlier studies show that its spectral weight quickly diminishes
with hole-doping [37, 38, 39], and becomes negligible in Cu $L_{3}$ RIXS when
doping is beyond 0.08 [25]. In addition, the bimagnon intensity in Cu $L_{3}$
RIXS is maximal at the zone center, and becomes significantly weaker at large
momentum [25]. The low-energy RIXS intensity at a certain momentum
$\mathbf{q}$ is thus expressed as a linear combination of the spectral weights
of the single-spin-flip and non-spin-flip components modified by their
corresponding azimuthal dependence:

$I_{\mathrm{RIXS}}(E,\phi,\boldsymbol{\epsilon}_{\mathrm{i}})=w_{\mathrm{s}}(E){\cdot}A_{\mathrm{s}}(E,\phi,\boldsymbol{\epsilon}_{\mathrm{i}})+w_{\mathrm{c}}(E){\cdot}A_{\mathrm{c}}(E,\phi,\boldsymbol{\epsilon}_{\mathrm{i}})$
(1)

where $w_{\mathrm{s(c)}}(E)$ is the spectral weight of spin (charge)
components, and $A_{\mathrm{s(c)}}(E,\phi,\boldsymbol{\epsilon}_{\mathrm{i}})$
is the azimuthal dependence which is already known from the scattering tensor
and self-absorption effects, and $\boldsymbol{\epsilon}_{\mathrm{i}}$ is the
incident polarization. Fig. 2c-f show the $\phi$ dependence of RIXS intensity
at constant energy loss of 0 and -0.25 eV with the decomposed spin and charge
contributions $w_{\mathrm{s(c)}}{\cdot}A_{\mathrm{s(c)}}(\phi)$. As can be
seen, the RIXS azimuthal dependences can be well fitted by these two
components, which verifies that the above analysis correctly describes the
low-energy RIXS response in Bi2212. In addition, the quasi-elastic scattering
at 0 energy loss is dominated by charge-like $\phi$ dependence, consistent
with the charge nature of the quasi-elastic peak. In Fig. 2g,h, we compare the
RIXS spectra at two special geometries, grazing incidence ($\phi$ = 0∘) with
$\sigma$ polarization and grazing emission ($\phi$ = 180∘) with $\pi$
polarization, with the decomposed spin and charge components. The grazing
emission with $\pi$ polarization is usually used to measure the magnetic
excitations in cuprates, since the charge component is largely suppressed in
this geometry as shown in Fig. 2h. Nonetheless, the charge component is still
considerable, which could influence the correct evaluation of the profile and
intensity of magnetic excitations. It is therefore necessary to fully
disentangle the spin and charge components to precisely study their nature. We
note here that the obtained spectral functions $w_{\mathrm{s(c)}}(E)$ are
solely related to the properties of the studied samples, as the angle and
polarization related geometry factors and the self-absorption effect are
removed by the knowledge of
$A_{\mathrm{s(c)}}(E,\phi,\boldsymbol{\epsilon}_{\mathrm{i}})$. This allows
the direct and unambiguous comparison between $w_{\mathrm{s(c)}}(E)$ and
theoretical calculations based on different models, which could provide vital
knowledge to understand the spin and charge excitations in cuprates.

Fig. 3: Momentum and doping dependence of the decomposed spin and charge
spectral functions. a – g (i – o) display the spin (charge) spectral functions
for samples of different doping. The error bars indicate the combined errors
of the fitting, statistic, and errors in determining the zero-energy-loss
positions (See Methods). The solid lines in a – g indicate the fittings by a
damped harmonic oscillator model. h Energy integrated intensities of the spin
spectral functions. p Energy integrated intensities of the charge spectral
functions with the quasi-elastic peaks subtracted. The error bars in h and p
represent the integration errors assuming that the spectral errors have an
energy correlation defined by the resolution function (see Methods).

Fig. 3 presents the obtained spin and charge spectral functions
$w_{\mathrm{s(c)}}(E)$ of the three different doped samples at different
momenta $\mathbf{q}$. Fig. 3a-g show the decomposed spin spectral functions,
which all show a single peak with a damped profile. Error bars describe the
fitting errors in the disentanglement (see Methods). Fig. 3h plots the energy
integrated intensity $I_{\mathrm{s}}(\mathbf{q})$ of the spin spectral
functions. The $I_{s}(\mathbf{q})$ of all three samples show a similar
$\mathbf{q}$ dependence: $I_{\mathrm{s}}(\mathbf{q})$ monotonically increases
with increasing $\mathbf{q}$, and has a slightly larger intensity when
approaching large $\mathbf{q}$ along ($\uppi$, $\uppi$) direction than
($\uppi$, 0) direction. This $\mathbf{q}$ dependence is qualitatively
consistent with the results in a previous study which calibrate the geometry
influences by comparing to the INS results [27]. Fig. 3i-o show the decomposed
charge response at different $\mathbf{q}$. There are two main components in
the charge response: a quasi-elastic peak including the low-energy phonons
close to zero energy loss, and a broad peak around -0.4 eV which extends to
high energy loss. The quasi-elastic peak is enhanced at (0.13, 0.13), which is
due to the structure modulation at (0.125, 0.125) in Bi2212 samples. Fig. 3p
shows the integrated intensity of the broad peak after the quasi-elastic peak
is subtracted. In contrast to the spin response, this charge response shows
much stronger intensity increases along ($\uppi$, 0) direction while it
saturates around (0.25, 0.25) along ($\uppi$, $\uppi$) direction. The
difference further highlights the distinct nature of the two decomposed
components.

By comparing different doping levels, one can notice that the spin excitations
show different development along the ($\uppi$, 0) and ($\uppi$, $\uppi$)
directions: the total spectral weight increases with increasing doping at
intermediate q along ($\uppi$, $\uppi$), while it slightly decreases along
($\uppi$, 0) direction, as shown by $I_{\mathrm{s}}(\mathbf{q})$ in Fig. 3h.
On the other hand, the spectral weights of the decomposed charge excitations
simply increases with increasing doping along both directions, while the
increase is more remarkable along ($\uppi$, 0) direction than ($\uppi$,
$\uppi$) as shown in Fig. 3p. A previous RIXS study on single layer
(Bi,Pb)2(Sr,La)2CuO6+δ [26] also investigated the influence of doping on the
spin excitations using the grazing-emission and incident $\pi$-polarization
geometry which enhances the scattering contribution from the spin-flip
channel. In contrast, they found the intensity of spin excitations increases
with doping at small and intermediate $\mathbf{q}$ along both ($\uppi$, 0) and
($\uppi$, $\uppi$) directions, but crosses over to decrease at large
$\mathbf{q}$. These different results could originate from either the
differences between single and double layer cuprates, or a small residual
mixture with charge excitations in the grazing-emission and incident
$\pi$-polarization geometry. Note that the ($\uppi$, 0) direction will endure
more influences from the residual charge excitations as the charge excitations
are more intense and doping-dependent along the ($\uppi$, 0) direction (see
Fig. 3p). The distinct spin excitations response to the hole doping along
($\uppi$, 0) and ($\uppi$, $\uppi$) directions could provide important
information for the understanding of the spin dynamics in doped cuprates,
which can be obtained by comparing with state-of-the-art theoretical
calculations discussed in detail in the next part of the paper.

Fig. 4: Fitting parameters of the spin spectral functions in Fig. 3 by a
damped harmonic oscillator model. a Bare frequency $\omega_{0}$. b Propagation
frequency $\omega_{\mathrm{p}}$. c Amplitude $A$. d Damping $\gamma$. The
error bars represent the fitting errors.

For a more quantitative analysis of the spin response, we fit the spin
spectral functions $w_{\mathrm{s(c)}}(E)$ by a damped harmonic oscillator
(DHO) model convoluted with a resolution function (see Methods). As shown by
the solid lines in Fig. 3a–g, the results can be well fitted by the DHO model.
Fig. 4 presents the fitting results. The bare frequency $\omega_{0}$ is
similar for all three dopings, while the damping $\gamma$ increases with
increasing doping and has a much larger value along ($\uppi$, $\uppi$)
direction. This is consistent with previous studies on doped cuprates showing
that the magnetic excitations are much more damped along ($\uppi$, $\uppi$)
direction [14, 15, 17, 20, 21, 22, 23, 26, 27]. With the charge excitations
excluded, we can now rule out that the over-damped profile of the spin
excitations along ($\uppi$, $\uppi$) direction comes from an increase of
charge contributions with doping. We note that the fitted damping factors of
the two smallest momenta of OD2 sample are large and out of the main trend of
the momentum dependence. We attribute this to the decomposed shape of spin
excitations bearing more influences from the uncertainties in the azimuthal
fitting when the spin excitation peak is getting closer to the elastic peak at
small momenta. Fig. 4b shows the propagation frequency defined as
$\omega_{\mathrm{p}}=\sqrt{\omega_{0}^{2}-\gamma^{2}}$, where a zero value is
assigned when the system is over damped, i.e. $\omega_{0}<\gamma$. One can see
that the spin excitations in the over-doped OD2 sample are fully over damped
along the ($\uppi$, $\uppi$) direction. Fig. 4c shows the fitted proportional
amplitude $A$ to DHO model. It increases with increasing doping along
($\uppi$, $\uppi$) direction while it changes little along ($\uppi$, 0)
direction, which is a bit different from the integrated total spectral weight
$I_{\mathrm{s}}(\mathbf{q})$ shown in Fig. 3h. This is mostly due to
$I_{\mathrm{s}}(\mathbf{q})$ including both the effects from the proportional
amplitude $A$ and the damping $\gamma$, while $A$ excludes the effect of
damping $\gamma$ which suppresses the total intensity.

Unravelling the character of the spin excitations in doped cuprates
The model of choice to study the evolution of the spin excitations upon doping
the cuprates is the “celebrated” $t$-$J$-like model [40], defined by the
Hamiltonian on a 2D square lattice:

$\displaystyle H=$ $\displaystyle-t\sum_{\langle
i,j\rangle}(\tilde{c}_{i}^{\dagger}\tilde{c}_{j}+h.c.)-t^{\prime}\sum_{\langle\langle
i,j\rangle\rangle}(\tilde{c}_{i}^{\dagger}\tilde{c}_{j}+h.c.)$
$\displaystyle-t^{\prime\prime}\sum_{\langle\langle\langle
i,j\rangle\rangle\rangle}(\tilde{c}_{i}^{\dagger}\tilde{c}_{j}+h.c.)+J\sum_{\langle
i,j\rangle}\left(\mathbf{S}_{i}\cdot\mathbf{S}_{j}-\frac{1}{4}\tilde{n}_{i}\tilde{n}_{j}\right),$
(2)

where $\tilde{c}_{i}^{\dagger}$ operator creates an electron at site $i$ in
the constrained Hilbert space without double electron occupancies,
$\mathbf{S}_{i}$ is a spin-1/2 operator at site $i$ and $\tilde{n}_{i}$ is the
on-site electron density at site $i$:
$\tilde{n}_{i}=\tilde{c}_{i}^{\dagger}\tilde{c}_{i}$ . The model parameters
$t$, $t^{\prime}$ and $t^{\prime\prime}$ denote the hopping integrals between
first, second, and third neighbors, respectively, whereas $J$ is the
antiferromagnetic Heisenberg interaction between nearest neighbor spins. In
our calculations, we take a realistic and widely-accepted (cf. [41] or a quite
similar choice in [42]) choice of the values of the
$t$-$t^{\prime}$-$t^{\prime\prime}$-$J$ model parameters $t^{\prime}=-0.3t$,
$t^{\prime\prime}=0.15t$ and $J=0.4t$. These values slightly differ from the
doping-dependent values suggested for Bi2212 in [43], but we keep these more
standard values to make our study universal. Furthermore, to better understand
the role of the longer-range hoppings, we consider switching off the
$t^{\prime\prime}$ hopping (the $t$-$t^{\prime}$-$J$ model), both the
$t^{\prime\prime}$ and $t^{\prime}$ hoppings (the $t$-$J$ model) as well as
substantially increasing the value of $t^{\prime}$ in the $t$-$t^{\prime}$-$J$
model calculations. Finally, note that, while to describe electronic
properties of the cuprates probably the charge-transfer ($pd$) model would be
more appropriate, the spin excitations are believed to be well-described by
models with oxygens being integrated out [44, 31, 26, 44, 45, 46, 47, 48, 49].
The $t$-($t^{\prime}$-$t^{\prime\prime}$-)$J$ model on a square lattice has
been intensively studied as a host of the superconducting state for a long
time (for example, see [50, 51] and references therein), and a possible
existence of the superconducting phase in a wide range of hole doping has been
proposed [52].

Next our goal here is to compute the static spin structure factor
$S(\mathbf{q})$, typically defined as:

$\displaystyle S(\mathbf{q})$
$\displaystyle=\langle\mathbf{S}(-\mathbf{q})\cdot\mathbf{S}(\mathbf{q})\rangle=\frac{1}{N}\sum_{i,j}\langle\mathbf{S}_{i}\cdot\mathbf{S}_{j}\rangle\mathrm{e}^{i\mathbf{q}(\textbf{r}_{i}-\textbf{r}_{j})},$
(3)

where the $i$, $j$ indices run over all sites, $N$ is the number of sites and
$\textbf{r}_{i}$ denotes the position of the site in the lattice. This is done
by involving state-of-the-art DMRG calculations which are carried out on an
$N=6\times 6$ square lattice with open boundary conditions (OBC). We chose a
6x6 OBC cluster to investigate a wide range of parameters with high accuracy.
The use of OBC enables us to avoid an artificial enhancement of specific
period correlations, which frequently occurs in periodic and cylindrical
systems. However, due to the OBC, the charges tend to localize at the edges
when charge imbalance, i.e. holes, is introduced. To counterbalance this
effect, we introduce an edge factor $\lambda$ that multiplies the electronic
hopping parameters $t$, $t^{\prime}$ and $t^{\prime\prime}$ as well as the
spin exchange coupling $J$ acting on the sites on the perimeter of our
cluster, see Methods and Supplementary Note 6.

By using the real space approach we have full control over which contributions
to $S(\mathbf{q})$ we include in Equation 3. In fact, we can select the
distance $\ell$ at which the sum in Equation 3 is truncated. In particular, it
is possible to consider the different total averages for the different spin-
spin correlations $\langle\mathbf{S}\cdot\mathbf{S}\rangle_{\ell}$ with
$\ell=1,2,3,\dots$ defining the considered neighbors. The difference in
$S(\mathbf{q})$ between using the singular values of the spin-spin
correlations and the averages is minimal (see Supplementary Note 7). While
shielding us from accessing $S(\mathbf{q},\omega)$, this approach allows us to
thoroughly study the possible magnonic character of the persistent spin
excitations.

Fig. 5: Comparison between the static spin structure factor in experiment and
theory a The experimental proxy for the static spin structure factor, i.e. the
integrated intensity of the spin excitations as seen by RIXS on Bi2212, see
Fig. 3h. The points are normalised to the value of the intensity at
$\mathbf{q}=(0.13,0.13)$ for the OD2 doping level. The error bars represent
the integration errors assuming that the spectral errors have an energy
correlation defined by the resolution function (see Methods). b Static spin
structure factor $S(\mathbf{q})$ obtained using DMRG on a $6\times 6$ cluster
(see text for further details) for the $t$-$t^{\prime}$-$t^{\prime\prime}$-$J$
model and three different hole-doping levels. The results are normalised to
the value of $S(\mathbf{q})$ at $\mathbf{q}\simeq(0.13,0.13)$ for the doping
level $n=0.22$. Model parameters: $J=0.4t$, $t^{\prime}=-0.3t$,
$t^{\prime\prime}=0.15t$. c Schematic comparison between the face-value of the
experimental and numerical sequences of intensities: points refer to RIXS
experiments and lines to the DMRG numerical results. The sequence at
$\mathbf{q}=(0.28,0.28)$ momentum refer to the sequence seen in the numerical
results at $\mathbf{q}=(0.267,0.267)$, see text for further details. Only
$\mathbf{q}=(0.18,0)$ momentum cannot be reproduced by our theoretical
calculations.

We now discuss how our theoretical calculations compare to the experimental
results presented in Fig. 3. Our aim within the calculations is to reproduce
the switching in the sequence of intensities upon doping. Hence, our focus is
on this qualitative aspect of the experimental results, rather than on the
quantitative reproduction of the experimental data within our theoretical
calculations. The main results are shown in Fig. 5: whereas Fig. 5a presents
the experimental integrated intensity of spin spectral weights for the three
measured doping levels UD, OD1, and OD2 (see above), Fig. 5b shows the
calculated $S(\mathbf{q})$ with the $t$-$t^{\prime}$-$t^{\prime\prime}$-$J$
model in the restricted Brillouin zone kinematically available to the
experiments. Finally, we schematically compare the doping evolution of the
experimental and theoretical intensities at all experimental momenta in Fig.
5c. The crucial message here is that, overall, there is good qualitative
agreement between theory and experiment. First, the experimentally observed
small anisotropy between the ($\uppi$, $\uppi$) and ($\uppi$, 0)
directions—the ($\uppi$, $\uppi$) direction shows larger intensities than
($\uppi$, 0) at high $\mathbf{q}$—is also reproduced by our calculations,
although it is not as small as in the RIXS experiment. Second, at five crucial
momentum points the theoretical calculations give the same sequence of
intensities of the spin structure factor as a function of doping as the face
value of the experimental results (see five grey circles in Fig. 5b and Fig.
5c). In addition, the sequence of intensities at $\mathbf{q}=(0.28,0.28)$ can
be reproduced with a slightly smaller momentum (see dashed grey circle in Fig.
5b and dashed lines in Fig. 5c). This indicates that a small fine-tuning of
the model parameters might give a full agreement between theory and experiment
also at this momentum. Last but not least we note that the aforementioned
sequence of intensities, and therefore the agreement with the theoretical
results, is largely confirmed also when the experimental error bars are taken
into account (see Supplementary Note 2).

Besides the overall agreement, there are two important discrepancies between
the experiments and calculations. First, the experimentally observed sequence
of intensities at $\mathbf{q}=(0.18,0)$ is not reproduced by the calculations.
Although the calculations do produce a crossing to a decreasing sequence upon
doping at a relatively large momentum [${\bf q}>(0.28,0)$], which satisfies
the experimental results at $\mathbf{q}=(0.33,0)$, they fail to reproduce the
intensity sequence at the intermediate momentum $\mathbf{q}=(0.18,0)$. The
experimental results thus indicate this crossing should happen at a smaller
momentum [${\bf q}<(0.18,0)$]. Within the $t$-$J$-like model and the parameter
sets we considered in this work, the current choice with both $t^{\prime}$ and
$t^{\prime\prime}$ gives the best agreement (see discussion on Fig. 6 below).
To fully reconcile this discrepancy, one may need to further fine-tune the
Hamiltonian parameters, and in particular, consider their doping-dependence
[43].

Second, the slope of the momentum dependence of the integrated experimental
intensities is smaller compared with the calculated $S(\mathbf{q})$. The slope
extrapolates to a non-zero residual intensity at the Brillouin zone center
${\bf q}=(0,0)$, making it gap-like, which is not found in the numerical
results. The reason for this apparent disagreement is likely the inter-layer
interaction in this bilayer cuprate, which gives rise to the onset of the
optical and acoustic magnon branches [53, 54, 7, 34]. While the 2D
antiferromagnetic acoustic spin wave has a zero structure factor at the zone
center [${\bf q}=(0,0)$], the optical branch due to the inter-layer coupling
will show a non-zero intensity, thus leading to a ‘gap’ in the zone center. As
our calculations do not include the coupling between the layers, this gap
feature is absent. However, due to the inter-layer interaction being much
smaller than the intra-layer one, the two branches quickly merge as
$\mathbf{q}$ moves away from the zone center ($\gtrsim$ 0.1 r.l.u.) [53, 54,
7], and their total intensities will be dominated by the intra-layer
parameters. Therefore, the comparison between our calculations and the
experiments is in reasonable agreement apart from the zone center.

Fig. 6: Theoretical static spin structure factor and the longer-range
electronic hoppings. Static spin structure factor $S(\mathbf{q})$ obtained
using DMRG on a $6\times 6$ cluster (see text for further details) for three
different hole-doping levels and a the $t$-$J$ model, b the
$t$-$t^{\prime}$-$J$, and c the $t$-$t^{\prime}$-$t^{\prime\prime}$-$J$ model.
Model parameters as in Fig. 5. Note the enlarged momentum coverage w.r.t. Fig.
5 and different scales of $S(\mathbf{q})$ for the ($\uppi$, $\uppi$) and
($\uppi$, 0) directions of the Brillouin zone.

The agreement between the $t$-$t^{\prime}$-$t^{\prime\prime}$-$J$ model and
experiments can be appreciated even more after looking at Fig. 6, where we
present the comparison between the theoretical results of the three different
$t$-$J$-like models. Fig. 6a shows $S(\mathbf{q})$ calculated for the ‘bare’
$t$-$J$ model. Within this model, the sequence of intensities of the spin
structure factor as a function of doping does not change in the experimental
momentum range, unlike the results of the experiments. Fig. 6b shows the same
quantity as Fig. 6a, but for the $t$-$t^{\prime}$-$J$ model. In this case, the
calculated $S(\mathbf{q})$ in the ($\uppi$, $\uppi$) direction shows the same
qualitative behavior as in the experimental case. However, the spin structure
factor $S(\mathbf{q})$ in ($\uppi$, 0) direction keeps increasing with doping
in the studied momentum range, which is still inconsistent with the
experiments. A different parameter choice with a larger value of $t^{\prime}$
in the $t$-$t^{\prime}$-$J$ model also does not improve the agreement between
calculations and experiments (See Supplementary Figure 4 and Note 4). The
decreasing sequence of intensity upon doping in the ($\uppi$, 0) direction is
only achieved after including the third neighbor hopping $t^{\prime\prime}$ as
shown in Fig. 6c. However, it appears at relatively larger $\mathbf{q}$ than
the experimental results. Nevertheless, it suggests the importance of long
range hoppings in achieving better agreement between experiment and theory.
Further improvements may require fine-tuning of the parameters (see above).

Fig. 7: Theoretical static spin structure factor and the short-range magnetic
correlations. a The $6\times 6$ cluster used in the DMRG calculations and the
nearest neighbor (NN), next nearest neighbor (NNN) and third neighbor
($3^{rd}$N) spin bonds; the color scale applied to the sites shows the value
of the average spin-spin correlation for those particular neighbors. b Static
spin structure factor $S(\mathbf{q})$ obtained using DMRG on a $6\times 6$
cluster (see text for further details) for three different hole-doping levels
and for the $t$-$t^{\prime}$-$t^{\prime\prime}$-$J$ model (parameters as in
Fig. 5) with keeping only first, second and third neighbor spin-spin
correlations in the Fourier Transform of Eq. (5). Note the enlarged momentum
coverage w.r.t. Fig. 5 and different scales of $S(\mathbf{q})$ for the
($\uppi$, $\uppi$) and ($\uppi$, 0) directions of the Brillouin zone.

Due to our real space approach, we have full control over the contribution of
spin correlations of different range to the static spin structure factor. In
fact, we can cut the summation in the Fourier transform [cf. Eqs. (5-6) in the
Methods] to include only up to a certain number of neighbors. This analysis
allows us to investigate the least effective range of magnetic correlations
that can qualitatively produce the spin structure factor $S(\mathbf{q})$ upon
doping. In Fig. 7, we show the main results of this analysis on the
$t$-$t^{\prime}$-$t^{\prime\prime}$-$J$ model. Fig. 7a is a cartoon
description of the real space spin-spin correlations between one sample site
around the centre of the $6\times 6$ cluster and its neighbors up to third-
nearest ones. The color scale represents the value of such real space
correlations. In Fig. 7b, we plot the spin structure factor $S(\mathbf{q})$
calculated by including only up to third-neighbor spin-spin correlations.
Since the short-range correlations mostly account for the large-$\mathbf{q}$
dynamics, they become poor in sketching the small-$\mathbf{q}$ properties,
which leaves an artificial gap around $\mathbf{q}=(0,0)$. Except the
small-$\mathbf{q}$ region, where the gap appears, the results show an
ascending intensity as a function of doping in the intermediate $\mathbf{q}$
range in both the ($\uppi$, $\uppi$) and ($\uppi$, 0) directions and a
crossover to descending sequence in larger $\mathbf{q}$, which qualitatively
agree with the spin structure factor $S(\mathbf{q})$ calculated from the full
Fourier transform depicted in Fig. 6c. We also notice that almost all of the
spectral weight of $S(\mathbf{q})$ at ($\uppi$, $\uppi$) is already accounted
for once solely the short-range magnetic correlations up to the third
neighbors are taken into account, cf. Fig. 5b and Fig. 7c. This is in stark
contrast with the undoped case (see Supplementary Figure 3): in the latter
case, the spectral weight at ($\uppi$, $\uppi$) is strongly underestimated
when only the short-range magnetic correlations are taken into account. This
is due to the importance of long-range correlations in the ordered
antiferromagnet stabilized at half-filling. An in-depth discussion of these
properties, as well as a similar analysis with different range of correlations
for the $t$-$J$ and $t$-$t^{\prime}$-$J$ models can be found in the
Supplementary Note 3.

DISCUSSION
Our main experimental result is the unambiguous assessment of the momentum and
doping evolution of the disentangled intrinsic spin and charge excitations
measured by Cu $L_{3}$-edge RIXS in doped Bi2212 samples. The disentangled
spin responses show profiles that can be well fitted by a damped harmonic
oscillator model, although they become over damped along ($\uppi$, $\uppi$)
direction. The momentum and doping dependence of the spin and charge responses
are different, implying the distinct nature of the two responses. In addition,
the disentangled spin excitations well persist into the over doped sample,
which confirms that RIXS indeed observes the persistence of spin excitations
in a large part of the Brillouin zone in doped cuprates.

The obtained experimental results on the spin excitations are qualitatively
reproduced by numerical simulations of the $t$-$J$-like model. It turns out
that the bare $t$-$J$ model is not enough and instead this model has to be
supplemented by longer-range hoppings—which points out the decisive role of
such hoppings in reproducing the experimentally observed paramagnons in doped
cuprates. Furthermore, the extensive real space analysis shows that short-
range magnetic correlations are needed in order to cause the observed
persistence of spin excitations in doped cuprates, meaning they need to be
paramagnonic in nature. From that, two important consequences follow: On one
hand, within the class of cuprate models with localized spins, those without
any magnetic correlations at all seem not to be realistic for doped cuprates.
(We note in passing that the class of models with localized spins is not only
restricted to the studied $t$–$J$ models. The Hubbard (or charge transfer)
model description of the doped cuprates also possesses localized moments whose
value is not substantially reduced compared with the $t$–$J$ model case, since
the number of doubly occupied sites is substantially suppressed in the
Hubbard-like models with realistic values of the on-site Coulomb repulsion,
cf. Fig. 1 of [55] or Fig. 6 of [56].) On the other hand, this means that
longer-range magnetic correlations do not play a crucial role in the doped
cuprates. Altogether, this helps in resolving the paradox related to the
persistence of the spin excitations upon doping the cuprates—despite a rapid
collapse of the long-range magnetic correlations.

We also briefly comment on the charge excitations we obtained in the
disentanglement analysis. With the same model calculations as for the spin
structure factor, we can get the charge structure factor (see Supplementary
Figure 5), which can reproduce the increase in intensities as a function of
doping in both the ($\uppi$, 0) and ($\uppi$, $\uppi$) directions. However,
the intensity anisotropy between the ($\uppi$, 0) and ($\uppi$, $\uppi$) is
missing. This might signal the onset of the bimagnons in the charge channel of
the RIXS spectrum [44, 31] or suggest that the longer-range Coulomb
interactions are important [57].

On the theory side, there are three important implications of the results
presented in this paper: The first one is that this work shows that a recent
theoretical suggestion that the spin excitations are responsible for the
T-linear dependence of the electronic scattering in the Hubbard model [58]
might indeed become a realistic scenario for the cuprates. This follows from
the above-stated conclusion that, without any ambiguities, the collective spin
excitations persist in the doped cuprates. The second one follows from the
suggested crucial role played by the longer-range hoppings in the $t$-$J$
models. Such a result goes in line with, inter alia, recent works advocating
for the strong sensitivity of the phase diagram of the $t$-$J$ like models to
the value of the next-nearest neighbor hopping $t^{\prime}$ [59, 60] and thus
‘sweetens the bad news’ coming from the study suggesting the lack of
superconductivity in the ground state of the 2D Hubbard model without longer-
range hoppings [61]. The third point relates to the fact that, as stated
above, solely the short-range magnetic correlations are needed to explain the
persistence of the intensity of the paramagnons in doped cuprates, similar to
the recently established well-defined magnons in the random $t$-$J$ model up
to $33\%$ hole doping [62]. Thus, an interesting task for theory would be to
gain an intuitive understanding of why the short-range magnetic correlation
alone can lead to the lack of changes of the paramagnons along the ($\uppi$,
0) direction.

We close the paper by presenting the impact of the results presented above on
the issue of pairing mechanism in the cuprate superconductors. While there are
many competing theories describing this issue, one of the widely-spread
concepts suggests that the pairing is mediated by the magnetic excitations
close to ($\uppi$, $\uppi$) momentum [63]. More precisely, it was shown by
inter alia Nocera et al. [47] and Huang et al. [64] that the spin fluctuations
mediated pairing is mostly due to the low-energy spin excitations with
momentum ${\bf q}$ = ($\uppi$, $\uppi$) and that the ‘doping-persistent’ high-
energy magnetic excitations away from that momentum, i.e. the ones which are
observed by RIXS, are not important to pairing. This finding has been
corroborated by the experimental results of Meyers et al. [22] and Dean et al.
[8]. Interestingly, here we have shown that even the ${\bf q}$ = ($\uppi$,
$\uppi$) paramagnons are largely a result of the short-range magnetic
correlations. This is because the spectral weight in the magnetic response
close to the ${\bf q}$ = ($\uppi$, $\uppi$) momentum is very similar both in
the ‘full’ spin structure factor calculations as well as in the ones
containing solely the short-range magnetic correlations, see Fig. 6(c) versus
Fig 7. Note that this is in contrast to the undoped case, for which the peak
around ${\bf q}$ = ($\uppi$, $\uppi$) is known to be hugely sensitive to the
longer-range magnetic correlations (see Supplementary Figure 3). Since the
nature of the ($\uppi$, $\uppi$) paramagnons is central to any of the spin
fluctuations mediated superconducting pairing mechanism, this finding plays an
important role in our deeper understanding of the puzzle of superconductivity
in cuprates. Finally, as a side note on the issue of pairing, the present
study stresses the importance of the proper choice of the longer range
hoppings ($t^{\prime}$ and $t^{\prime\prime}$) in any realistic cuprate
modelling. Thus, we believe that any model trying to explain superconductivity
in the cuprates should include realistic values of these parameters as slight
variations might qualitatively affect the observed physics.

METHODS
RIXS Experiments and samples
The RIXS experiments were carried out with the SAXES spectrometer at the
ADRESS beamline of the Swiss Light Source at the Paul Scherrer Institut [65,
66]. The incident X-ray energy was set at the Cu $L_{3}$ resonance peak at
approximately 933 eV. The instrument resolution was determined by the elastic
peak measured on carbon tape, giving an overall energy resolution of $\sim$
100 meV full width half maximum (FWHM). The single-crystal samples of
Bi2Sr2CaCu2O8+δ were prepared by the floating zone method as described in Ref.
[67]. The samples were cleaved at base temperature to present fresh surfaces
before the measurements. All the data were collected at base temperature
$\sim$ 24 K. The momentum transfer $\mathbf{q}$ is denoted in reciprocal
lattice units (r. l. u.) using the pseudo-tetragonal unit cell with a = b =
3.82 Å.

Azimuthal dependent RIXS measurements
In two-dimensional superconducting cuprates, the single-spin-flip and non-
spin-flip excitations in the 3$d_{x^{2}-y^{2}}$ orbital as well as the other
$d$-$d$ excitations in the Cu $L_{3}$-edge RIXS spectra show distinct geometry
and polarization related characters, which are determined by their different
scattering tensors [35, 68, 69]. This allows to assess the nature of these
excitations by either resolving the polarizations of both the incident and
scattered photons [70] or by azimuthal dependent RIXS measurements [33]. In
Fig.1b, we show the scattering geometry and the sample rotation in the
azimuthal dependent experiments. The directions of the incident and emitted
x-rays are fixed through the scattering angle to 130∘, thus the total momentum
transfer is fixed at $\mathbf{q}$. Two linear polarizations ($\sigma$ and
$\uppi$) are used for the incident x-rays while the polarization of the
emitted light is not resolved. The plate-like sample is mounted on a wedged
sample holder (with wedge angle $\theta_{w}$ = 10∘, 20∘, 40∘ and 50∘ in the
experiments) to have a certain in-plane momentum transfer. The azimuthal
rotation axis is parallel to the total momentum transfer $\mathbf{q}$, so that
the projections of $\mathbf{q}$ in the sample reciprocal frame are unchanged
during rotation, while the projections of the photon polarization are
changing. This allows measuring the azimuthal dependence of the excitations at
fixed momentum in the sample momentum space. When rotating the sample, the
photon polarization will be continuously rotated in the sample space, and the
scattering tensors will then result in different rotation dependences for
different excitations.

Error estimations of the experimental data
The major errors of the RIXS spectra are the statistical errors of the photon
counting, which are expressed as the square roots of total photon counts.
Another error comes from the uncertainties in determining the zero-energy-loss
positions in the spectra, which is done in examining the positions of the
elastic peaks. Here we assume that this error is about $\pm 5$% of the
resolution, which is about $\pm 5$ meV. This error is converted to the error
of spectral intensity by multiplying the derivative of the spectrum. Other
random errors are accounted for by the fitting errors (95% confidence
interval) in the azimuth-dependent fitting. All the above errors are treated
as independent at each energy-loss point and summed to a total error by the
square root of their sum of squares. The errors of the integrated intensities
are obtained by assuming that the errors of the decomposed spectral functions
have an energy correlation defined by the resolution function. The error bars
in the fitting of the damped harmonic oscillator (Fig.4) are the fitting
errors.

Fitting by damped harmonic oscillator model
The formula of the damped harmonic oscillator (DHO) model used for the fitting
of the spin spectral funtions in Fig. 3 is:

$A{\cdot}\frac{\gamma\omega}{(\omega^{2}-\omega_{0}^{2})^{2}+4\gamma^{2}\omega^{2}}.$
(4)

A Gaussian resolution function with 100 meV FWHM is convoluted with the above
DHO model to fit the results. The fitting energy range is [-0.8, 0.4] eV for
all data except those with the smallest momenta ($\mathbf{q}$ = (0.09, 0) and
(0.065, 0.065)), which are fitted in [-0.65, 0.4] eV. This is to reduce the
influence of the long tail at high energy loss.

DMRG calculations
In numerical calculations with OBC cluster, a proper correction is often added
into the Hamiltonian to minimize the effects of missing terms at open edge. In
this study, we introduced the edge factor to uniformize the mobility of
charge.

The correct edge factor is calculated for each doping level and each different
model ($t$-$J$, $t$-$t^{\prime}$-$J$,
$t$-$t^{\prime}$-$t^{\prime\prime}$-$J$). We extrapolate it by computing the
dispersion $\delta=n_{\text{in}}-n_{\text{out}}$ where $n_{\text{in}}$ is the
averaged electron density taken over the sites which do not belong to the
edges and $n_{\text{out}}$ is the averaged electron density taken over the
sites which form the edge of the cluster. We computed the dispersion $\delta$
for different values of the edge factor $\lambda$ for each doping level and
model and take the final value of the edge factor $\lambda$ as that at which
the dispersion $\delta=0$. Nevertheless, the obtained values ($\lambda\sim
0.9-1.2$) are close enough to 1 to smoothly connect the inside and edge of the
cluster. Furthermore, by introducing this edge factor the Friedel oscillations
as well as the most important finite-size effect coming from using OBC can be
significantly reduced, as can be seen in the Supplementary Note 6.

We keep up to $m=7000$ states in the DMRG calculations, leading to an error
$\epsilon/N=10^{-6}$. We make sure that the local density in the system is
isotropic by using the edge factor $\lambda$ as described above and we compute
the real space spin-spin correlations
$\langle\mathbf{S}_{i}\cdot\mathbf{S}_{j}\rangle$ for all pairs $(i,j)$
labelling the system’s sites. This way, we can compute the spin static factor
$S(\mathbf{q})$ as:

$S(\mathbf{q})=\frac{1}{N}\sum_{i,j}\langle\mathbf{S}_{i}\cdot\mathbf{S}_{j}\rangle\cos(q_{x}(x_{i}-x_{j})+q_{y}(y_{i}-y_{j})).$
(5)

This method is only valid when the local $z$-component of the spins is small
enough, e.g. $S^{z}_{i}\leq 10^{-5}$ for all sites $i$. However, for certain
doping levels and models, it was not possible to reach this level of accuracy
for the local value of $S^{z}$. Therefore, we renormalized the $S^{z}S^{z}$
correlations and computed $S^{z}(\mathbf{q})$ as:

$\displaystyle S^{z}(\mathbf{q})$ $\displaystyle=\frac{1}{N}\sum_{i,j}(\langle
S^{z}_{i}\cdot S^{z}_{j}\rangle-\langle S^{z}_{i}\rangle\langle
S^{z}_{j}\rangle)$
$\displaystyle\times\cos(q_{x}(x_{i}-x_{j})+q_{y}(y_{i}-y_{j})).$ (6)

Due to the symmetries of the considered models,
$S(\mathbf{q})=3S^{z}(\mathbf{q})$. Equation 6 was used for the doping level
$n=0.22$ for both the $t$-$t^{\prime}$-$J$ and the
$t$-$t^{\prime}$-$t^{\prime\prime}$-$J$ models (see Fig. 6).

DATA AVAILABILITY
All data needed to evaluate the conclusions in the paper are present in the
paper and/or the Supplementary Materials. Additional raw experiment data to
reproduce the analysis as well as data and scripts to reproduce the
theoretical results are available in [71].

CODE AVAILABILITY
The code for the numerical calculations will be made available from the
corresponding authors upon reasonable request.

## References

  * [1] Keimer, B., Kivelson, S. A., Norman, M. R., Uchida, S. & Zaanen, J.  From quantum matter to high-temperature superconductivity in copper oxides.  _Nature_ 518, 179–186 (2015).
  * [2] Lee, P. A., Nagaosa, N. & Wen, X.-G.  Doping a Mott insulator: Physics of high-temperature superconductivity.  _Rev. Mod. Phys._ 78, 17–85 (2006).
  * [3] Xu, G. _et al._ Testing the itinerancy of spin dynamics in superconducting Bi2Sr2CaCu2O8+δ.  _Nat. Phys._ 5, 642–646 (2009).
  * [4] Damascelli, A., Hussain, Z. & Shen, Z.-X.  Angle-resolved photoemission studies of the cuprate superconductors.  _Rev. Mod. Phys._ 75, 473–541 (2003).
  * [5] Tranquada, J. M.  Spins, stripes, and superconductivity in hole-doped cuprates.  _AIP Conference Proceedings_ 1550, 114–187 (2013).
  * [6] Braicovich, L. _et al._ Magnetic excitations and phase separation in the underdoped ${\mathrm{La}}_{2-x}{\mathrm{Sr}}_{x}{\mathrm{CuO}}_{4}$ superconductor measured by resonant inelastic X-Ray scattering.  _Phys. Rev. Lett._ 104, 077002 (2010).
  * [7] Le Tacon, M. _et al._ Intense paramagnon excitations in a large family of high-temperature superconductors.  _Nat. Phys._ 7, 725–730 (2011).
  * [8] Dean, M. P. M. _et al._ Persistence of magnetic excitations in La2-xSrxCuO4 from the undoped insulator to the heavily overdoped non-superconducting metal.  _Nat. Mater._ 12, 1019–1023 (2013).
  * [9] Dean, M. P. M. _et al._ High-energy magnetic excitations in the cuprate superconductor ${\mathrm{Bi}}_{2}{\mathrm{Sr}}_{2}{\mathrm{CaCu}}_{2}{\mathrm{O}}_{8\mathrm{+}\delta}$: Towards a unified description of its electronic and magnetic degrees of freedom.  _Phys. Rev. Lett._ 110, 147001 (2013).
  * [10] Dean, M. P. M. _et al._ Magnetic excitations in stripe-ordered La1.875Ba0.125CuO4 studied using resonant inelastic x-ray scattering.  _Phys. Rev. B_ 88, 020403(R) (2013).
  * [11] Le Tacon, M. _et al._ Dispersive spin excitations in highly overdoped cuprates revealed by resonant inelastic X-Ray scattering.  _Phys. Rev. B_ 88, 020501(R) (2013).
  * [12] Ishii, K. _et al._ High-energy spin and charge excitations in electron-doped copper oxide superconductors.  _Nat. Commun._ 5, 3714 (2014).
  * [13] Lee, W. S. _et al._ Asymmetry of collective excitations in electron- and hole-doped cuprate superconductors.  _Nat. Phys._ 10, 883–889 (2014).
  * [14] Dean, M. P. M. _et al._ Itinerant effects and enhanced magnetic interactions in Bi-based multilayer cuprates.  _Phys. Rev. B_ 90, 220506(R) (2014).
  * [15] Guarise, M. _et al._ Anisotropic softening of magnetic excitations along the nodal direction in superconducting cuprates.  _Nat. Commun._ 5, 5760 (2014).
  * [16] Minola, M. _et al._ Collective nature of spin excitations in superconducting cuprates probed by resonant inelastic X-Ray scattering.  _Phys. Rev. Lett._ 114, 217003 (2015).
  * [17] Wakimoto, S. _et al._ High-energy magnetic excitations in overdoped ${\mathrm{La}}_{2-x}{\mathrm{Sr}}_{x}{\mathrm{CuO}}_{4}$ studied by neutron and resonant inelastic X-Ray scattering.  _Phys. Rev. B_ 91, 184513 (2015).
  * [18] Peng, Y. Y. _et al._ Magnetic excitations and phonons simultaneously studied by resonant inelastic X-Ray scattering in optimally doped ${\mathrm{Bi}}_{1.5}{\mathrm{Pb}}_{0.55}{\mathrm{Sr}}_{1.6}{\mathrm{La}}_{0.4}{\mathrm{CuO}}_{6+\delta}$.  _Phys. Rev. B_ 92, 064517 (2015).
  * [19] Ellis, D. S. _et al._ Correlation of the superconducting critical temperature with spin and orbital excitations in $({\mathrm{Ca}}_{x}{\mathrm{La}}_{1-x})({\mathrm{Ba}}_{1.75-x}\text{La}_{0.25+x}){\mathrm{Cu}}_{3}{\mathrm{O}}_{y}$ as measured by resonant inelastic X-Ray scattering.  _Phys. Rev. B_ 92, 104507 (2015).
  * [20] Huang, H. Y. _et al._ Raman and fluorescence characteristics of resonant inelastic X-Ray scattering from doped superconducting cuprates.  _Sci Rep_ 6, 19657 (2016).
  * [21] Monney, C. _et al._ Resonant inelastic X-Ray scattering study of the spin and charge excitations in the overdoped superconductor ${\mathrm{La}}_{1.77}{\mathrm{Sr}}_{0.23}{\mathrm{CuO}}_{4}$.  _Phys. Rev. B_ 93, 075103 (2016).
  * [22] Meyers, D. _et al._ Doping dependence of the magnetic excitations in ${\mathrm{La}}_{2-x}{\mathrm{Sr}}_{x}{\mathrm{CuO}}_{4}$.  _Phys. Rev. B_ 95, 075139 (2017).
  * [23] Ivashko, O. _et al._ Damped spin excitations in a doped cuprate superconductor with orbital hybridization.  _Phys. Rev. B_ 95, 214508 (2017).
  * [24] Minola, M. _et al._ Crossover from collective to incoherent spin excitations in superconducting cuprates probed by detuned resonant inelastic X-Ray scattering.  _Phys. Rev. Lett._ 119, 097001 (2017).
  * [25] Chaix, L. _et al._ Resonant inelastic X-Ray scattering studies of magnons and bimagnons in the lightly doped cuprate ${\mathrm{La}}_{2-x}{\mathrm{Sr}}_{x}{\mathrm{CuO}}_{4}$.  _Phys. Rev. B_ 97, 155144 (2018).
  * [26] Peng, Y. Y. _et al._ Dispersion, damping, and intensity of spin excitations in the monolayer ${(\text{Bi,Pb})}_{2}{(\text{Sr,La})}_{2}{\mathrm{CuO}}_{6+\delta}$ cuprate superconductor family.  _Phys. Rev. B_ 98, 144507 (2018).
  * [27] Robarts, H. C. _et al._ Anisotropic damping and wave vector dependent susceptibility of the spin fluctuations in ${\mathrm{La}}_{2-x}{\mathrm{Sr}}_{x}{\mathrm{CuO}}_{4}$ studied by resonant inelastic X-Ray scattering.  _Phys. Rev. B_ 100, 214510 (2019).
  * [28] James, A. J. A., Konik, R. M. & Rice, T. M.  Magnetic response in the underdoped cuprates.  _Phys. Rev. B_ 86, 100508(R) (2012).
  * [29] Benjamin, D., Klich, I. & Demler, E.  Single-band model of resonant inelastic X-Ray scattering by quasiparticles in high-${T}_{c}$ cuprate superconductors.  _Phys. Rev. Lett._ 112, 247002 (2014).
  * [30] Kanász-Nagy, M., Shi, Y., Klich, I. & Demler, E. A.  Resonant inelastic X-Ray scattering as a probe of band structure effects in cuprates.  _Phys. Rev. B_ 94, 165127 (2016).
  * [31] Jia, C., Wohlfeld, K., Wang, Y., Moritz, B. & Devereaux, T. P.  Using RIXS to uncover elementary charge and spin excitations.  _Phys. Rev. X_ 6, 021020 (2016).
  * [32] Tsutsui, K. & Tohyama, T.  Incident-energy-dependent spectral weight of resonant inelastic X-Ray scattering in doped cuprates.  _Phys. Rev. B_ 94, 085144 (2016).
  * [33] Kang, M. _et al._ Resolving the nature of electronic excitations in resonant inelastic X-Ray scattering.  _Phys. Rev. B_ 99, 045105 (2019).
  * [34] Peng, Y. Y. _et al._ Influence of apical oxygen on the extent of in-plane exchange interaction in cuprate superconductors.  _Nat. Phys._ 13, 1201–1206 (2017).
  * [35] Ament, L. J. P., van Veenendaal, M., Devereaux, T. P., Hill, J. P. & van den Brink, J.  Resonant inelastic X-Ray scattering studies of elementary excitations.  _Rev. Mod. Phys._ 83, 705–767 (2011).
  * [36] Bisogni, V. _et al._ Bimagnon studies in cuprates with resonant inelastic X-Ray scattering at the O ${K}$ edge. I. Assessment on La2CuO4 and comparison with the excitation at Cu ${L}_{3}$ and Cu ${K}$ edges.  _Phys. Rev. B_ 85, 214527 (2012).
  * [37] Devereaux, T. P. & Hackl, R.  Inelastic light scattering from correlated electrons.  _Rev. Mod. Phys._ 79, 175–233 (2007).
  * [38] Sugai, S., Suzuki, H., Takayanagi, Y., Hosokawa, T. & Hayamizu, N.  Carrier-density-dependent momentum shift of the coherent peak and the LO phonon mode in p-type high-${T}_{c}$ superconductors.  _Phys. Rev. B_ 68, 184504 (2003).
  * [39] Sugai, S. _et al._ Superconducting pairing and the pseudogap in the nematic dynamical stripe phase of La2-xSrxCuO4.  _J. Phys.: Condens. Matter_ 25, 475701 (2013).
  * [40] Chao, K. A., Spałek, J. & Oleś, A. M.  Canonical perturbation expansion of the Hubbard model.  _Phys. Rev. B_ 18, 3453–3464 (1978).
  * [41] Delannoy, J.-Y. P., Gingras, M. J. P., Holdsworth, P. C. W. & Tremblay, A.-M. S.  Low-energy theory of the $t$-$t{{}^{\prime}}$-$t{{}^{\prime\prime}}$-${U}$ hubbard model at half-filling: Interaction strengths in cuprate superconductors and an effective spin-only description of ${\mathrm{La}}_{2}{\mathrm{CuO}}_{4}$.  _Phys. Rev. B_ 79, 235130 (2009).
  * [42] Sénéchal, D. & Tremblay, A.-M. S.  Hot spots and pseudogaps for hole- and electron-doped high-temperature superconductors.  _Phys. Rev. Lett._ 92, 126401 (2004).
  * [43] Drozdov, I. K. _et al._ Phase diagram of Bi2Sr2CaCu2O8+δ revisited.  _Nat. Commun._ 9, 5210 (2018).
  * [44] Jia, C. J. _et al._ Persistent spin excitations in doped antiferromagnets revealed by resonant inelastic light scattering.  _Nat. Commun._ 5, 3314 (2014).
  * [45] Eder, R., Ohta, Y. & Maekawa, S.  Anomalous spin and charge dynamics of the $t$-${J}$ model at low doping.  _Phys. Rev. Lett._ 74, 5124–5127 (1995).
  * [46] Tohyama, T., Horsch, P. & Maekawa, S.  Spin and charge dynamics of the $t$-${J}$ model.  _Phys. Rev. Lett._ 74, 980–983 (1995).
  * [47] Nocera, A., Patel, N. D., Dagotto, E. & Alvarez, G.  Signatures of pairing in the magnetic excitation spectrum of strongly correlated two-leg ladders.  _Phys. Rev. B_ 96, 205120 (2017).
  * [48] Tohyama, T., Mori, M. & Sota, S.  Dynamical density matrix renormalization group study of spin and charge excitations in the four-leg $t\text{$-$}{t}^{{}^{\prime}}\text{$-$}{J}$ ladder.  _Phys. Rev. B_ 97, 235137 (2018).
  * [49] Pärschke, E. M. _et al._ Numerical investigation of spin excitations in a doped spin chain.  _Phys. Rev. B_ 99, 205102 (2019).
  * [50] Jiang, S., Scalapino, D. J. & White, S. R.  Ground-state phase diagram of the $t-t^{\prime}-{J}$ model.  _Proc. Natl. Acad. Sci. U.S.A._ 118, e2109978118 (2021).
  * [51] Gong, S., Zhu, W. & Sheng, D. N.  Robust $d$-wave superconductivity in the square-lattice $t\text{$-$}{J}$ model.  _Phys. Rev. Lett._ 127, 097003 (2021).
  * [52] Jiang, H.-C., Weng, Z.-Y. & Kivelson, S. A.  Superconductivity in the doped $\mathit{t}-\mathit{J}$ model: Results for four-leg cylinders.  _Phys. Rev. B_ 98, 140505 (2018).
  * [53] Reznik, D. _et al._ Direct observation of optical magnons in YBa2Cu3O6.2.  _Phys. Rev. B_ 53, R14741(R) (1996).
  * [54] Hayden, S. M., Aeppli, G., Perring, T. G., Mook, H. A. & Doğan, F.  High-frequency spin waves in YBa2Cu3O6.15.  _Phys. Rev. B_ 54, R6905(R) (1996).
  * [55] Moreo, A., Scalapino, D. J., Sugar, R. L., White, S. R. & Bickers, N. E.  Numerical study of the two-dimensional Hubbard model for various band fillings.  _Phys. Rev. B_ 41, 2313–2320 (1990).
  * [56] Kung, Y. F. _et al._ Doping evolution of spin and charge excitations in the Hubbard model.  _Phys. Rev. B_ 92, 195108 (2015).
  * [57] Fidrysiak, M. & Spałek, J.  Unified theory of spin and charge excitations in high-${T}_{\textrm{c}}$ cuprate superconductors: A quantitative comparison with experiment and interpretation.  _Phys. Rev. B_ 104, L020510 (2021).
  * [58] Wú, W., Wang, X. & Tremblay, A.-M.  Non-Fermi liquid phase and linear-in-temperature scattering rate in overdoped two-dimensional Hubbard model.  _Proc. Natl. Acad. Sci. U.S.A._ 119, e2115819119 (2022).
  * [59] Jiang, H.-C. & Devereaux, T. P.  Superconductivity in the doped Hubbard model and its interplay with next-nearest hopping $t^{\prime}$.  _Science_ 365, 1424–1428 (2019).
  * [60] Jiang, Y.-F., Zaanen, J., Devereaux, T. P. & Jiang, H.-C.  Ground state phase diagram of the doped hubbard model on the four-leg cylinder.  _Phys. Rev. Research_ 2, 033073 (2020).
  * [61] Qin, M. _et al._ Absence of superconductivity in the pure two-dimensional Hubbard model.  _Phys. Rev. X_ 10, 031016 (2020).
  * [62] Kumar, A., Sachdev, S. & Tripathi, V.  Quasiparticle metamorphosis in the random $t\text{$-$}{J}$ model.  _Phys. Rev. B_ 106, L081120 (2022).
  * [63] Scalapino, D. J.  A common thread: The pairing interaction for unconventional superconductors.  _Rev. Mod. Phys._ 84, 1383–1417 (2012).
  * [64] Huang, E. W., Scalapino, D. J., Maier, T. A., Moritz, B. & Devereaux, T. P.  Decrease of $d$-wave pairing strength in spite of the persistence of magnetic excitations in the overdoped Hubbard model.  _Phys. Rev. B_ 96, 020503 (2017).
  * [65] Ghiringhelli, G. _et al._ SAXES, a high resolution spectrometer for resonant X-Ray emission in the 400-1600ev energy range.  _Rev. Sci. Instrum._ 77, 113108 (2006).
  * [66] Strocov, V. N. _et al._ High-resolution soft X-ray beamline ADRESS at the Swiss Light Source for resonant inelastic X-ray scattering and angle-resolved photoelectron spectroscopies.  _J. Synchrotron Rad._ 17, 631–643 (2010).
  * [67] Weyeneth, S., Schneider, T. & Giannini, E.  Evidence for Kosterlitz-Thouless and three-dimensional ${XY}$ critical behavior in ${\text{Bi}}_{2}{\text{Sr}}_{2}{\text{CaCu}}_{2}{\text{O}}_{8+\delta}$.  _Phys. Rev. B_ 79, 214504 (2009).
  * [68] Ament, L. J. P., Ghiringhelli, G., Sala, M. M., Braicovich, L. & van den Brink, J.  Theoretical demonstration of how the dispersion of magnetic excitations in cuprate compounds can be determined using resonant inelastic X-Ray scattering.  _Phys. Rev. Lett._ 103, 117003 (2009).
  * [69] Haverkort, M. W.  Theory of resonant inelastic X-Ray scattering by collective magnetic excitations.  _Phys. Rev. Lett._ 105, 167404 (2010).
  * [70] Fumagalli, R. _et al._ Polarization-resolved Cu ${L}_{3}$-edge resonant inelastic X-Ray scattering of orbital and spin excitations in ${\mathrm{NdBa}}_{2}{\mathrm{Cu}}_{3}{\mathrm{O}}_{7-\delta}$.  _Phys. Rev. B_ 99, 134517 (2019).
  * [71] Zhang, W. _et al._ Unravelling the nature of the spin excitations disentangled from the charge contributions in a doped cuprate superconductor.  _zenodo.7286412_ (2022).
  * [72] Bricogne, G. _et al._ X-ray absorption, refraction and resonant scattering tensors in selenated protein crystals: implications for data collection strategies in macromolecular crystallography.  _J. Appl. Cryst._ 38, 168–182 (2005).
  * [73] Nakajima, R., Stöhr, J. & Idzerda, Y. U.  Electron-yield saturation effects in ${L}$-edge X-ray magnetic circular dichroism spectra of Fe, Co, and Ni.  _Phys. Rev. B_ 59, 6421–6429 (1999).
  * [74] Henneken, H., Scholze, F. & Ulm, G.  Lack of proportionality of total electron yield and soft X-ray absorption coefficient.  _J. Appl. Phys._ 87, 257–268 (2000).
  * [75] Ruosi, A. _et al._ Electron sampling depth and saturation effects in perovskite films investigated by soft X-ray absorption spectroscopy.  _Phys. Rev. B_ 90, 125120 (2014).
  * [76] Ronay, M. _et al._ A new correlation for ${T_{\text{c}}}$ from Cu 2p absorption.  _Solid State Communications_ 77, 699–702 (1991).
  * [77] Nücker, N. _et al._ Site-specific and doping-dependent electronic structure of ${\mathrm{YBa}}_{2}$${\mathrm{Cu}}_{3}$${\mathrm{O}}_{\mathit{x}}$ probed by O 1s and Cu 2p X-ray-absorption spectroscopy.  _Phys. Rev. B_ 51, 8529–8542 (1995).
  * [78] Chabot-Couture, G. _et al._ Polarization dependence and symmetry analysis in indirect ${K}$-edge RIXS.  _Phys. Rev. B_ 82, 035113 (2010).
  * [79] Achkar, A. J. _et al._ Bulk sensitive X-ray absorption spectroscopy free of self-absorption effects.  _Phys. Rev. B_ 83, 081106(R) (2011).

Acknowledgements The experiments were performed at the ADRESS beamline of the
Swiss Light Source at the Paul Scherrer Institut (PSI). The experimental work
at PSI is supported by the Swiss National Science Foundation through project
no. 200021$\\_$178867, and the Sinergia network Mott Physics Beyond the
Heisenberg Model (MPBH) (SNSF Research Grants CRSII2$\\_$160765/1 and
CRSII2$\\_$141962). The theoretical work is supported by the Narodowe Centrum
Nauki (NCN, Poland) under Project Nos. 2016/22/E/ST3/00560 (C.E.A. and K.W.),
2021/40/C/ST3/00177 (C.E.A.) and 2016/23/B/ST3/00839 (K.W.). We acknowledge
support by the Interdisciplinary Centre for Mathematical and Computational
Modeling (ICM), University of Warsaw (UW), within grant no G81-4. S.N.
acknowledges support from SFB 1143 project A05 (project-id 247310070) of the
Deutsche Forschungsgemeinschaft. T.C.A. acknowledges funding from the European
Union’s Horizon 2020 research and innovation programme under the Marie
Skłodowska-Curie grant agreement No. 701647 (PSI-FELLOW-II-3i program). We
acknowledge the valuable discussions with Jonathan Pelliciari, Mingu Kang,
Riccardo Comin and Ekaterina Pärschke. We kindly thank U. Nitzsche for
technical support. For the purpose of Open Access, the author has applied a
CC-BY public copyright licence to any Author Accepted Manuscript (AAM) version
arising from this submission

Competing Interests: The authors declare that they have no competing
interests.

Author Contributions: W.Z. and C.E.A. contributed equally to this work. W.Z.
and T.S. conceived the project. K.W. conceived the theory approach. W.Z.,
Y.T., T.C.A., E.P. and T.S. performed the RIXS experiments with the support of
V.S. C.E.A. performed the DMRG calculations with the support of S.N. E.G. grew
and characterized the single crystals. W.Z. and C.E.A. analyzed the data in
discussion with T.S., S.N., and K.W. W.Z., C.E.A., K.W., and T.S. wrote the
manuscript with the input from all authors. T.S. and K.W. coordinated the
research.

Supplementary Note 1: Absorption coefficients and self-absorption correction

Self-absorption is always present in soft X-ray RIXS measurements on bulk
samples. It will largely modify the RIXS intensity and has to be considered in
azimuthal dependent measurement. To derive its contribution, one needs the
knowledge of the linear absorption coefficient $\mu(E,\boldsymbol{\epsilon})$
of the studied sample, which depends on the energy ($E$) of the X-rays as well
as the orientation of the polarization vector ($\boldsymbol{\epsilon}$) in the
sample space. For a single crystal, the absorption coefficient can be
expressed in a tensor form, which is constrained by the point group symmetry
of the crystal. The Bi2212 sample has a tetragonal structure, leading to a
diagonal absorption tensor with only two independent elements, $f_{aa}(E)$ and
$f_{cc}(E)$, which correspond to the absorption coefficients with polarization
vector in the sample $a$-$b$ plane and along the out-of-plane $c$ axis,
respectively [72]. The absorption coefficient with arbitrary polarization
vector is,

$\mu(E,\boldsymbol{\epsilon})=\boldsymbol{\epsilon}^{T}\cdot\begin{pmatrix}f_{aa}(E)&0&0\\\
0&f_{aa}(E)&0\\\ 0&0&f_{cc}(E)\end{pmatrix}\cdot\boldsymbol{\epsilon}.$ (S1)

In this study, we use X-ray absorption spectroscopy (XAS) measured by total
electron yield (TEY) to evaluate the absorption tensor elements $f_{aa}(E)$
and $f_{cc}(E)$ of the Bi2212 samples around the Cu $L_{3}$-edge. The inset of
Fig. S1 displays the geometry of the measurements. The in-plane element
$f_{aa}(E)$ is probed by $\sigma$ polarized X-rays, while the out-of-plane
element $f_{cc}(E)$ is probed at grazing incidence and $\pi$ polarization. As
the fully grazing incidence is not possible, here we use 20∘ incident angle,
which gives a result of
${\sin}^{2}20^{\circ}{\cdot}f_{aa}(E)+{\cos}^{2}20^{\circ}{\cdot}f_{cc}(E)$.
The TEY signal is usually proportional to the absorption coefficient given
that the electron escape depth $L$ is much smaller than the photon penetration
depth $\lambda{\cdot}{\sin}\alpha$, where $\lambda$=1/$\mu$ and $\alpha$ is
the incident angle [73, 74, 75]. Although $L\ll\lambda$ is usually true in
cuprates [75], ${\sin}\alpha$ will reduce the photon penetration depth when
$\alpha$ is small at grazing incidence, which will reduce the proportionality
especially at the resonance where $\mu$ is large. Such a saturation effect can
be described by the following relation [73, 74, 75]:

$\text{TEY}(E,\alpha,\boldsymbol{\epsilon})\propto\frac{ML}{\lambda(E,\boldsymbol{\epsilon})\cdot\sin\alpha}\cdot\frac{1}{1+L/(\lambda(E,\boldsymbol{\epsilon})\cdot\sin\alpha)}.$
(S2)

where $M$ is a material dependent factor and the second fraction on the right
side expresses the saturation effect. One can see that this saturation factor
will distort the shape of TEY when (1)
$L/[\lambda(E,\boldsymbol{\epsilon})\cdot\sin\alpha]\sim$ 1 and (2) the
$\lambda(E,\boldsymbol{\epsilon})$ is strongly energy dependent such as around
the absorption edge. To address the possible saturation effect in our samples,
we measured the TEY at several different incident angles at $\sigma$
polarization, as shown in Supplementary Figure 1(a) – (c). The saturation
effect is relatively small with only a slight intensity reduction on the
resonance peak in the UD and OD1 samples but becomes more obvious in the OD2
sample. Supplementary Figure 1(d) shows the incident angle dependence of
$\text{TEY}\cdot\sin\alpha$ at the resonance peak and pre-edge background. The
solid lines are the fitting by a function of $A/(1+1/(r\cdot\sin\alpha))$, in
which $r$ gives an estimate of $\lambda/L$. The fitted values of $r$ are 22,
40, and 11 for the UD, OD1, and OD2 samples, respectively. We note that this
is just a rough estimate due to the limited points of measured angles.
Nevertheless, these values qualitatively agree with the previous measurement
on a YBCO film with $r$$\sim$20 [75]. The results suggest that the TEY at
normal incidence is nearly proportional to the absorption coefficient. So we
take the TEY at normal incidence and $\sigma$ polarization as the in-plane
absorption tensor element $f_{aa}(E)$. For the TEY at $\alpha=20^{\circ}$ and
$\pi$ incidence, some saturation effects exist. However, as the resonance peak
amplitude is relatively small ($\sim 15\%$) compared to the pre-edge
background, the saturation factor is much less energy dependent compared to
the $\sigma$ polarization. Thus, we can use an overall energy-independent
correction factor determined by the change of the pre-edge background from
$\alpha=90^{\circ}$ to $\alpha=20^{\circ}$ at $\sigma$ polarization to correct
the spectra. By subtracting the contribution of
${\sin}^{2}20^{\circ}{\cdot}f_{aa}(E)$, one can then get the out-of-plane
element $f_{cc}(E)$. Supplementary Figure 1(e) shows the obtained $f_{aa}(E)$
and $f_{cc}(E)$ of all three samples, which are normalized to make the pre-
edge background the same. As can be seen, there are almost no resonance peaks
in $f_{cc}(E)$, which confirms the in-plane 3$d_{x^{2}-y^{2}}$ character of
the samples. $f_{aa}(E)$ and $f_{cc}(E)$ are almost the same for the three
samples, except for a small intensity reduction in the resonance peak of
$f_{aa}(E)$ in OD2 sample. As the doped holes mainly go to the oxygen ligands,
this main peak should be mostly unaffected while the shoulder at a bit higher
energy ($\sim$934.2 eV) increases with the hole doping [76, 77]. The reduction
is thus most likely due to some residual saturation effects in normal
incidence in the OD2 sample due to its smaller $\lambda/L$ ratio at resonance.
In our analysis, we use the same set of $f_{aa}(E)$ and $f_{cc}(E)$ from the
OD1 sample for all three samples.

The self-absorption effect in RIXS is described by the following formula [78,
79, 33]:

$I_{\mathrm{t}}^{\mathrm{exp}}(E,\phi,\boldsymbol{\epsilon}_{\mathrm{i}})=\sum_{\boldsymbol{\epsilon}_{\mathrm{f}}}\frac{I_{\mathrm{t}}(E,\phi,\boldsymbol{\epsilon}_{\mathrm{i}},\boldsymbol{\epsilon}_{\mathrm{f}})}{\mu_{\mathrm{i}}(E_{\mathrm{i}},\phi,\boldsymbol{\epsilon}_{\mathrm{i}})+\mu_{\mathrm{f}}(E_{\mathrm{f}},\phi,\boldsymbol{\epsilon}_{\mathrm{f}})\cdot\frac{-\widehat{\mathbf{k}}_{\mathrm{i}}\cdot\widehat{\mathbf{n}}_{(}\phi)}{\widehat{\mathbf{k}}_{\mathrm{f}}\cdot\widehat{\mathbf{n}}_{(}\phi)}}.$
(S3)

Here the
$I_{\mathrm{t}}^{\mathrm{exp}}(E,\phi,\boldsymbol{\epsilon}_{\mathrm{i}})$is
the experimental ($\mathrm{exp}$) RIXS intensity for a certain type
($\mathrm{t}$) of excitation which can be spin, charge or $d$-$d$ excitations,
and
$I_{\mathrm{t}}(E,\phi,\boldsymbol{\epsilon}_{\mathrm{i}},\boldsymbol{\epsilon}_{\mathrm{f}})$
is the intrinsic RIXS intensity with $E=E_{\mathrm{f}}-E_{\mathrm{i}}$ the
energy loss. The denominator describes the self-absorption effect, where
$\mu_{\mathrm{i(f)}}$ is the absorption coefficient of the incident (emitting)
photons, $\widehat{\mathbf{k}}_{\mathrm{i(f)}}$ is the unit vector of the
photon propagation direction, and $\widehat{\mathbf{n}}_{(}\phi)$ is the unit
vector normal to the sample surface. As the final polarization
($\boldsymbol{\epsilon}_{\mathrm{f}}$) is not resolved in this experiment, the
measured intensity is a sum of all the possible
$\boldsymbol{\epsilon}_{\mathrm{f}}$ ($\sigma$ and $\pi$ polarizations). The
azimuthal $\phi$ dependence of
$I_{\mathrm{t}}(E,\phi,\boldsymbol{\epsilon}_{\mathrm{i}},\boldsymbol{\epsilon}_{\mathrm{f}})$
can be calculated based on the scattering tensor of the studied excitation and
a rotation matrix which links the laboratory coordinates and sample
coordinates, as already demonstrated in reference [33]. Fig. 1c –d in the main
text show the $\phi$ dependence of
$\sum_{\boldsymbol{\epsilon}_{\mathrm{f}}}I_{\mathrm{t}}(E,\phi,\boldsymbol{\epsilon}_{\mathrm{i}},\boldsymbol{\epsilon}_{\mathrm{f}})$
for the non-spin-flip and spin-flip excitations. With the knowledge of the
absorption tensor of our samples,
$\mu_{\mathrm{i(f)}}(E_{\mathrm{i(f)}},\phi,\boldsymbol{\epsilon}_{\mathrm{i(f)}})$
can be determined precisely, and we can therefore calculate the self-
absorption effect and include it in the final azimuthal dependence
$I_{\mathrm{t}}^{\mathrm{exp}}(E,\phi,\boldsymbol{\epsilon}_{\mathrm{i}})$ for
a certain type of excitation. The results at an energy loss of -0.25 eV are
shown in Fig. 1e and f for the non-spin-flip and spin-flip excitations. The
final RIXS spectrum is a sum of all possible types of excitations
$\sum_{\mathrm{t}}I_{\mathrm{t}}^{\mathrm{exp}}(E,\phi,\boldsymbol{\epsilon}_{\mathrm{i}})$,
and $I_{\mathrm{t}}^{\mathrm{exp}}(E,\phi,\boldsymbol{\epsilon}_{\mathrm{i}})$
can be expressed as a product of the intrinsic spectral weight of the
excitation $w_{\mathrm{t}}(E)$ and the azimuthal dependent geometry factor
$A_{\mathrm{t}}(E,\phi,\boldsymbol{\epsilon}_{\mathrm{i}})$. In the low-energy
part of the spectra in cuprates, with the spin-flip and non-spin-flip
excitations related to 3$d_{x^{2}-y^{2}}$ orbital dominating, the RIXS
intensity is then explained by eq. (1), which can be decomposed into the two
components by their azimuthal dependence.

Supplementary Note 2: Consequences of including the error bars for the doping-
dependence of the sequence of experimental spin structure intensities.

In order to show the consequences of including the error bars for the doping-
dependence of the sequence of experimental spin structure intensities, we
start by presenting the error bars of the experimental results for each
available momentum and doping level, see Supplementary Table 1.

We note that we define the two spin structure intensities as different if the
error bars do not overlap, i.e., their differences are twice the widths of the
error bars. Thus, the doping development is clear for five momenta points
[$(0.33,0)$, $(0,18,0)$, $(0.09,0)$, $(0.13,0.13)$ and $(0.23,0.23)$] when the
error bars are taken into account. Although at some momenta [$(0.33,0)$ and
$(0.18,0.18)$ $(0.23,0.23)$], the error bars of OD1 and OD2 are overlapped,
they are quite well separated with UD ones, which means there is still clear
doping changes from UD to OD2. The intensities of different doping at (0.065,
0.065) and (0.28, 0.28) are quite similar within the error bars, but they are
consistent with the small difference in the calculated $S({\bf q})$.

Supplementary Table 1: Error bars of experimental spin structure factors (integrated spin excitation intensities obtained from RIXS) given in the same units as the integrated spin excitation intensities of Fig. 5a of the main text.  | Momentum
---|---
Doping |  (0.33, 0) |  (0.18, 0) |  (0.09, 0) |  (0.065, 0.065) |  (0.13, 0.13) |  (0.23, 0.23) |  (0.28, 0.28)
UD |  0.0305 |  0.0254 |  0.0259 |  0.0273 |  0.0302 |  0.0287 |  0.0302
OD1 |  0.0247 |  0.0264 |  |  |  0.0269 |  0.0190 |  0.0294
OD2 |  0.0254 |  0.0280 |  0.0267 |  0.0304 |  0.0271 |  0.0256 |  0.0295

Supplementary Note 3: Detailed discussions of the influence of the short-range
magnetic correlations and longer-range electronic hoppings on the static spin
structure factor

In what follows we study the impact of including solely a restricted number of
short-range spin-spin correlations in the Fourier transform defining the
static spin structure factor $S({\bf q})$, see Eq. (5) in the main text of the
paper. The main results are presented in Supplementary Figure 2 which shows
$S({\bf q})$ calculated in the three distinct $t$-$J$–like models considered
in this paper and considering up to the third neighbor spin-spin correlations
in the Fourier transform. (In the main text, we discuss the case of the
$t\text{-}t^{\prime}\text{-}t^{\prime\prime}\text{-}J$ model considering up to
the third neighbor spin-spin correlations.)

First of all, let us look at the results shown in the first column of
Supplementary Figure 2, i.e. once solely spin-spin correlations up to first-
neighbors are included (by definition this includes also on site spin-spin
correlations) when calculating $S({\bf q})$. Except for small quantitative
differences, the three considered versions of the $t$-$J$ model show the same
behavior, meaning the effect of the longer range hopping $t^{\prime}$ and
$t^{\prime\prime}$ is minimal when the cut-off in the Fourier transform is on
the nearest-neighbor correlations. This can be understood when one realizes
that first neighbors correlations are always relatively large and
antiferromagnetic and the subtle effects induced by the longer range hopping
terms are minimal. Moreover, due to the fact that the considered correlations
are antiferromagnetic (i.e. negative), a negative spectral weight sets in for
small ${\bf q}$ in $S(\mathbf{q})$. We would like to underline that this
result is not linked to the instability of the ground state, as also the
ground state in these calculations is still the exact ground state of each
considered model; instead, this is due to the approximate calculations of
$S({\bf q})$, namely the respective cuts in the Fourier series (see above). As
a side remark let us note that when one considers up to first neighbors in
$S(\mathbf{q})$, its behavior is qualitatively the same in the nodal and anti-
nodal directions of the Brillouin zone. In particular, the damping of the peak
at $(\pi,\pi)$ is recovered, but the same behavior happens at $(\pi,0)$.

Next, we move on to the second column of Supplementary Figure 2, where we
consider up to second-neighbor spin-spin correlations in the Fourier transform
for $S({\bf q})$. Now the two directions in the Brilloin zone, $(\pi,\pi)$ and
$(\pi,0)$, show different behaviors. The damping of the $(\pi,\pi)$ peak is
more prominent and the intensity at this point is increased for all three
models. No negative weight is present. The main features present in the
$S(\mathbf{q})$ are already present for the $t\text{-}J$ model. A distinct
behavior is now visible when longer range hoppings are included: in the
$(\pi,0)$ direction the different doping lines do not cross nor overlap, while
they do so once the “bare” $t\text{-}J$ model is considered. In the
$(\pi,\pi)$ direction, both the $t\text{-}t^{\prime}\text{-}J$ and
$t\text{-}t^{\prime}\text{-}t^{\prime\prime}\text{-}J$ models already show the
onset of a line crossing around $(\pi/2,\pi/2)$. However, as previously
stated, no crossing is present in the $(\pi,0)$ direction, meaning we need to
include further neighbor spin correlations in the Fourier transform to recover
this behavior.

Finally, we focus on the third column of Supplementary Figure 2, where up to
third-neighbor spin-spin correlations are included in the Fourier transform.
When considering the simpler $t\text{-}J$ model, there is again an onset of
negative weight, which has disappeared in the previous case. This comes from
the inclusion of third-neighbor spin-spin correlations, which are
antiferromagnetic, i.e., negative. Being quite strong, they are not
compensated by the ferromagnetic, i.e., positive, second-neighbors
correlations, which leads effectively to a negative weight. On the other hand,
the intensity of the $(\pi,\pi)$ peak is now doubled compared to the one seen
when only first-neighbor correlations are included. Qualitatively, there is
not a big difference if compared to the previous case where up to second-
neighbors correlations are being considered. If longer range hoppings are
included, the negative weight problem is solved. Let us now focus on the
$t\text{-}t^{\prime}\text{-}J$ case: in the $(\pi,\pi)$ direction, the
intensity of the peak has increased and it is comparable with that of the full
$S(\mathbf{q})$. Furthermore, the crossing at $(\pi/2,\pi/2)$ is present. If
we now look at the $(\pi,0)$ direction, we see the onset of a crossing close
to $(\pi,0)$. Compared to the full $S(\mathbf{q})$ case [see Fig. 6b], the
0.11 doping level is now “meeting” the other doping level lines at $(\pi,0)$.
The last case is the one which includes also the $t^{\prime\prime}$ hopping as
already shown in the main text. The behavior in the $(\pi,\pi)$ direction is
similar to that seen in the $t\text{-}t^{\prime}\text{-}J$ model just
discussed, therefore we will only examine what happens in the $(\pi,0)$
direction. Here the crossing close to $(\pi,0)$ is clearly visible and
includes all three doping levels, as seen in experiments and in the full
$S(\mathbf{q})$.

Lastly, we would like to comment on the presence of a gap at $(0,0)$ momentum
transfer observed in the results for both the second and third column of
Supplementary Figure 2: it is clear that when one considers the long-range
distances in real spaces, this translates to small $\mathbf{q}$ values in
$q$-space. Therefore, the “wrong” behavior seen at small $\mathbf{q}$ is
expected to improve more and more when further neighbors are included.

In conclusion, in order to understand the behavior of the integrated intensity
$S(\mathbf{q})$, it is important to include longer range hopping up to third
neighbors, but only shorter range spin-spin correlations are needed to recover
the main properties seen in the experimental data.

To show that the above analysis is only possible for the doped system, we
carried out a similar analysis for the undoped (half-filled) system. We
present results in Supplementary Figure 3. For a better comparison, we set the
$S(\mathbf{q})$ range to be the same for all panels. This cuts part of the
line in Supplementary Figure 3a, where we are showing the results for the
Fourier transform carried out on up to the first-neighbor correlation. For
small $\mathbf{q}$, the static spin structure factor turns out to be negative
and is thus cut out of the presented figure. Looking at the progression of
$S(\mathbf{q})$ when further neighbors are included in the Fourier transform,
we observe that the full spin structure factor $S(\mathbf{q})$ is gradually
better reproduced. Nevertheless, even when considering only up to third
neighbors in the pure antiferromagnetic case we cannot obtain a good
quantitative agreement with the full $S(\mathbf{q})$. This is because the
spectral weight of the peak at $(\pi,\pi)$ is still only at about half of its
value when up to third neighbors are included.

Supplementary Note 4: Effect of larger $t^{\prime}$ on the spin static
structure factor

To check the importance of long-range hopping $t^{\prime\prime}$ in our model,
we consider again the $t$-$t^{\prime}$-$J$ model, but we set
$t^{\prime}=-0.5t$ rather than $t^{\prime}=-0.3t$. While this value is too
large for cuprates, it allows us to understand the trend of the spin static
structure factor when $t^{\prime}$ is changed. We plot our results in
Supplementary Figure 4 and show them only in the Brillouin zone available for
the experiments. The first thing we notice is that the agreement in the
$(\pi,\pi)$ direction is lost. Indeed, the $0.22$ hole-doping line loses
intensity at a rather small $\mathbf{q}$, so that the sequence of intensities
at (0.28, 0.28) is completely different between theory and experiments, in
contrast to the results shown in Fig. 5b of the main text. When looking at the
$(\pi,0)$ direction, we see that the 0.22 hole-doping line does cross the
other two intensities, becoming the lowest line, which does not happen for
$t^{\prime}=-0.3t$. However, the disagreement at $\mathbf{q}=(0.18,0)$
persists also in this model. Hence, when compared to the option of including
longer range hopping $t^{\prime\prime}$, it seems that the latter is to be
preferred as the lowering of the intensity of the 0.22 hole-doping case is
achieved for larger $\mathbf{q}$ without compromising the agreement in the
$(\pi,\pi)$ direction.

Supplementary Note 5: Static charge structure factor

We calculate the static charge structure factor

$N(\mathbf{q})=\langle
n(-\mathbf{q})n(\mathbf{q})\rangle=\frac{1}{N}\sum_{i,j}(\langle
n_{i}n_{j}\rangle-\langle n_{i}\rangle\langle
n_{j}\rangle)\mathrm{e}^{i\mathbf{q}(\textbf{r}_{i}-\textbf{r}_{j})}.$ (S4)

for the $t$-$t^{\prime}$-$t^{\prime\prime}$-$J$ model on a $6\times 6$ cluster
using DMRG. This is done in a similar manner as the static spin structure
factor calculated elsewhere in this manuscript; in particular, the applied
edge factors for the different doping levels are the same as those applied to
compute the static spin structure factor $S(\mathbf{q})$.

We plot the charge structure factor using DMRG in the part of the Brillouin
zone available to the RIXS experiments and compare these results to the
experimentally extrapolated integrated intensities in the non-spin-flip RIXS
sector (Supplementary Figure 4). The sequence of intensities at the probed
momentum points is qualitatively well reproduced in the DMRG results for
$N(\mathbf{q})$. We believe that some differences, especially the much larger
relative intensity along the antinodal directions observed in the experiment,
can either be a result of other contribution besides the charge ones in the
non-spin-flip channel in the experiments (in particular, the bimagnon
contributions [44, 31]) or show that the longer-range Coulomb interactions are
important here [57].

Supplementary Note 6: Benchmarking the use of the edge factor in the numerical
method

In order to perform calculations on an open cluster, we have introduced the
edge factor $\lambda$ as described in the main text. To prove that this method
leads to correct results, we calculated the static spin structure factor
$S(\mathbf{q})$ on the half-filled $6\times 6$ cluster. The results are shown
in Supplementary Figure 3 and fully agree with the textbook behavior of the
Heisenberg model on a square lattice at zero temperature [a dominant peak at
$(\pi,\pi$)].

Moreover, as mentioned in the main text, using open clusters could provide
much better results than periodic clusters if the open edge is managed
properly. To confirm this statement, we compare the static spin structure
factor $S(q)$ between short and long open chains in the 1D $t$-$J$ chain with
$J/t=0.4$ and $n=2/3$, since the result for a large 2D system is not
available. In Supplementary Figure 5, the spin static structure factors for
$L=6$ and $L=60$ open chains are compared. The result for $L=60$ is expected
to be almost identical to taking the thermodynamic limit. For the $L=6$ open
cluster, small local potentials $0.15$ and $-0.3$ have been added on the first
and second edge sites, respectively, as edge factors used to achieve a uniform
density distribution, i.e., $<n_{i}>=2/3$ for all $i$, $i$ labeling the sites
in the chain. We find a good agreement between $S(q)$ with $L=6$ and $L=60$
except for the peak height around $q=2\pi/3$. This discrepancy is caused
mainly by an enhanced finite-size effect due to strong quantum fluctuations in
1D. Thus, such discrepancy will be much smaller in 2D systems. The structure
factor for an $L=6$ periodic chain is also plotted. It is difficult to compare
it with the thermodynamic limit result, because only discrete momenta are
allowed. Besides, the artificial enhancement (suppression) of $S(q)$ at
$q=2\pi/3$ (at $q\neq 2\pi/3$) is clearly seen. Consequently, we can suggest
that the use of an (small) open cluster is a reasonable way to capture the
overall features of $S(q)$, unless $S(q)$ shows a very complex behavior.

Supplementary Note 7: Averaged compared to non-averaged correlations

In the main text as well as in Supplementary Note 3, the results presented are
based on the Fourier transform of the averaged values of the different
neighbor spin-spin correlations. To support our choice of showing results
based on averaged correlations, Supplementary Figure 6 compares the spin
static structure factor $S(\mathbf{q})$ calculated using the Fourier transform
with keeping up to third neighbor spin-spin correlations for the
$t$-$t^{\prime}$-$t^{\prime\prime}$-$J$ model—with averaged (over the whole
cluster) and non-averaged values of the spin-spin correlations. The
differences between the two latter cases are almost invisible, suggesting that
our choice of considering only averaged correlations does not introduce any
additional approximations.

Supplementary Figure 1: X-ray absorption spectra and the extracted absorption
tensor elements $f_{aa}(E)$ and $f_{cc}(E)$ for the three samples. (a) -(c)
The incident-angle dependent XAS of UD, OD1 and OD2 samples, respectively. (d)
The intensity of $\text{TEY}\cdot\sin\alpha$ of the resonance peak and pre-
edge background as a function of incident angle $\alpha$. (e) The extracted
absorption tensor elements $f_{aa}(E)$ and $f_{cc}(E)$ of UD, OD1 and OD2
sample, respectively. Inset: geometry and photon polarizations of the
measurements. Supplementary Figure 2: Detailed contributions of the short-
range magnetic correlations and the longer-range electronic hoppings to the
theoretical static spin structure factor. Static spin structure factor
$S(\mathbf{q})$ obtained using DMRG on a $6\times 6$ cluster (see text for
further details) for the $t$-$J$ (top panels), the $t$-$t^{\prime}$-$J$
(middle panels), and $t$-$t^{\prime}$-$t^{\prime\prime}$-$J$ model (bottom
panels). In all cases only a restricted number of short-range spin-spin
correlations is nonzero in the Fourier Transform of Eq. 5: solely first
neighbors (left panels); solely first and second neighbors (middle panels);
and solely first, second, and third neighbors (right panels). Model parameters
as in Fig. 5. Note the enlarged momentum coverage w.r.t. Fig. 5 and different
scales of $S(\mathbf{q})$ for the nodal and anti-nodal directions of the
Brillouin zone. Supplementary Figure 3: Theoretical static spin structure
factor at half-filling. Static spin structure factor S(q) obtained using DMRG
on a $6\times 6$ cluster for the half-filled $t$–$J$ model and (as in all
other cases) including the edge factor $\lambda$ for a square lattice and
keeping only up to first (a panel), second (b panel) and third (c panel)
neighbor spin-spin correlations in the Fourier Transform. d Static spin
structure factor $S(\mathbf{q})$ obtained using DMRG on a $6\times 6$ cluster
for the half-filled $t$–$J$ model and (as in all other cases) including the
edge factor $\lambda$ for a square lattice. Note the enlarged momentum
coverage w.r.t. Fig. 5 and the same scale in the nodal and anti-nodal
direction as well as the agreement with the textbook result of $S(\mathbf{q})$
for the Heisenberg model. Supplementary Figure 4: a The experimental proxy for
the static spin structure factor, i.e. the integrated intensity of the spin
excitations as seen by RIXS on Bi2212. The points are normalised to the value
of the intensity at $\mathbf{q}=(0.13,0.13)$ for the OD2 doping level. b
Static spin structure factor $S(\mathbf{q})$ obtained using DMRG on a $6\times
6$ cluster (see text for further details) for the $t$-$t^{\prime}$-$J$ model
with $t^{\prime}=-0.5t$, $J=-0.4t$. The results are normalised to the value of
$S(\mathbf{q})$ at $\mathbf{q}\simeq(0.13,0.13)$ for the doping level
$n=0.22$. Supplementary Figure 5: a The experimental proxy for the static
charge structure factor, i.e. the integrated intensity of the charge
excitations as seen by RIXS on Bi2212, see Fig. 3p in the main text. The
points are normalised to the value of the intensity at
$\mathbf{q}=(0.13,0.13)$ for the OD2 doping level. b Static charge structure
factor $N(\mathbf{q})$ obtained using DMRG on a $6\times 6$ cluster (see text
for further details) for the $t$-$t^{\prime}$-$t^{\prime\prime}$-$J$ model and
three different hole-doping levels. The results are normalised to the value of
$N(\mathbf{q})$ at $\mathbf{q}\simeq(0.13,0.13)$ for the doping level
$n=0.22$. Model parameters: $J=0.4t$, $t^{\prime}=-0.3t$,
$t^{\prime\prime}=0.15t$. Supplementary Figure 6: Usage of OBC vs PBC
condition. We plot results for $S(q)$ for the $t$-$J$ model on a chain with
doping level $n=2/3$. Lines show results for OBC on chains of size $L=60$ and
$L=6$, edge factors are applied on the smaller cluster. The overall agreement
in the evolution of $S(q)$ is incredibily good. Open squares show results at
the allowed momenta for a closed chain with $L=6$. Overall, the $S(q)$ values
do not compare well with those of the $L=60$ cluster. Supplementary Figure 7:
Theoretical static spin structure factor including averaged and non-averaged
short-range magnetic correlations. Static spin structure factor
$S(\mathbf{q})$ obtained using DMRG on a $6\times 6$ cluster for the doped
$t$-$t^{\prime}$-$t^{\prime\prime}$-$J$ model with only up to third neighbor
spin-spin correlations included in the Fourier transform defining
$S(\mathbf{q})$. Top (bottom) panels shows results obtained using an averaged
(non-averaged) values of the short magnetic correlations, respectively. Model
parameters as in Fig. 5. Note the enlarged momentum coverage w.r.t. Fig. 5 and
different scales of $S(\mathbf{q})$ for the nodal and anti-nodal directions of
the Brillouin zone.