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

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 Bi$_2$Sr$_2$CaCu$_2$O$_{8+\delta}$ 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 * wenliang.zhang@psi.ch † krzysztof.wohlfeld@fuw.edu.pl ‡ thorsten.schmitt@psi.ch elementary excitations across the phase diagram. As suggested in numerous papers over the last 11 years , 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 (π, 0) direction similarly to the magnons in the antiferromagnetic phase with the damping increasing moderately, although they rapidly become overdamped with doping along the (π, π) direction and show almost nondispersive profiles [14, 15, 17, 20-23, 26, 27]. This suggests that in some particular area of the Brillouin zone, that is mostly in the (π, 0) direction, the magnetic excitations completely 'ignore' the existence of a critical value of the doping δ 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(q, ω) from RIXS spectra and also due to the problems in reliably calculating S(q, ω). 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 spinflip 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 Bi 2 Sr 2 CaCu 2 O 8+δ (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 (π, π) [(π, 0)] direction, which unequivocally confirms the persistence of the spin excitations in doped cuprates. We then compare the above experimental results to state-ofthe-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 secondand 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.

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 overdoped regime. The three samples are labeled as UD (T c = 73 K), OD1 (T c = 88 K) and OD2 (T 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  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 σ and π 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 f − E 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. 2a and b show the RIXS intensity map of OD2 sample at q = (0.33, 0) as a function of azimuthal angle φ and energy loss E at σ and π incident polarization, . 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 s(c) ·A s(c) (φ), 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 (φ = 0 • ) with σ polarization and grazing emission (φ = 180 • ) with π polarization. The error bars represent the combined errors including the statistic errors and the errors in determining the zero-energy-loss positions (See Methods).
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 3d x 2 −y 2 orbital. Although the double spin flip with a net spin change of zero (bimagnon with ∆S = 0) is also allowed in RIXS process in cuprates [25,35,36], 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 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: is the spectral weight of spin (charge) components, and A s(c) (E, φ, i ) is the azimuthal dependence which is already known from the scattering tensor and self-absorption effects, and i is the incident polarization. Fig. 2c-f show the φ dependence of RIXS intensity at constant energy loss of 0 and -0.25 eV with the decomposed spin and charge contributions w s(c) ·A s(c) (φ). 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 φ dependence, consistent with the charge nature of the quasielastic peak. In Fig. 2g,h, we compare the RIXS spectra at two special geometries, grazing incidence (φ = 0 • ) with σ polarization and grazing emission (φ = 180 • ) with π polarization, with the decomposed spin and charge components. The grazing emission with π 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 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 s(c) (E, φ, i ). This allows the direct and unambiguous comparison between w 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 presents the obtained spin and charge spectral functions w s(c) (E) of the three different doped samples at different momenta 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 s (q) of the spin spectral functions. The I s (q) of all three samples show a similar q dependence: I s (q) monotonically increases with increasing q, and has a slightly larger intensity when approaching large q along (π, π) direction than (π, 0) direction. This q dependence is qualitatively consistent with the results in a previous study which calibrate the geometry influences by comparing to the INS results [27]. 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 (π, 0) direction while it saturates around (0.25, 0.25) along (π, π) 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 (π, 0) and (π, π) directions: the total spectral weight increases with increasing doping at intermediate q along (π, π), while it slightly decreases along (π, 0) direction, as shown by I s (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 (π, 0) direction than (π, π) as shown in Fig. 3p. A previous RIXS study on single layer (Bi,Pb) 2 (Sr,La) 2 CuO 6+δ [26] also investigated the influence of doping on the spin excitations using the grazing-emission and incident π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 q along both (π, 0) and (π, π) directions, but crosses over to decrease at large 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 π-polarization geometry. Note that the (π, 0) direction will endure more influences from the residual charge excitations as the charge excitations are more intense and doping-dependent along the (π, 0) direction (see Fig. 3p). The distinct spin excitations response to the hole doping along (π, 0) and (π, π) directions could provide important information for the understanding of the spin dynamics in doped cuprates, which can be obtained by comparing with state-of-theart theoretical calculations discussed in detail in the next part of the paper.
For a more quantitative analysis of the spin response, we fit the spin spectral functions w 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 ω 0 is similar for all three dopings, while the damping γ increases with increasing doping and has a much larger value along (π, π) direction. This is consistent with previous studies on doped cuprates showing that the magnetic excitations are much more damped along (π, π) direction [14, 15, 17, 20-23, 26, 27]. With the charge excitations excluded, we can now rule out that the over-damped profile of the spin excitations along (π, π) 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 ω p = ω 2 0 − γ 2 , where a zero value is assigned when the system is over damped, i.e. ω 0 < γ. One can see that the spin excitations in the over-doped OD2 sample are fully over damped along the (π, π) direction. Fig. 4c shows the fitted proportional amplitude A to DHO model. It increases with increasing doping along (π, π) direction while it changes little along (π, 0) direction, which is a bit different from the integrated total spectral weight I s (q) shown in Fig. 3h. This is mostly due to I s (q) including both the effects from the proportional amplitude A and the damping γ, while A excludes the effect of damping γ 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: wherec † i operator creates an electron at site i in the constrained Hilbert space without double electron occupancies, S i is a spin-1/2 operator at site i andñ i is the on-site electron density at site i:ñ i =c † ic i . The model parameters t, t and t 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 -t -J model parameters t = −0.3t, t = 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 hopping (the tt -J model), both the t and t hoppings (the t-J model) as well as substantially increasing the value of t in the t-t -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 [26,31,[44][45][46][47][48][49]. The t-(t -t -)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(q), typically defined as: where the i, j indices run over all sites, N is the number of sites and 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×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 λ that multiplies the electronic hopping parameters t, t and t 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(q) we include in Equation 3. In fact, we can select the distance 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 S·S with = 1, 2, 3, . . . defining the considered neighbors. The difference in S(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(q, ω), this approach allows us to thoroughly study the possible magnonic character of the persistent spin excitations.
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 calculated S(q) with the t-t -t -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 (π, π) and (π, 0) directions-the (π, π) direction shows larger intensities than (π, 0) at high 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 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 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 [q > (0.28, 0)], which satisfies the experimental results at q = (0.33, 0), they fail to reproduce the intensity sequence at the intermediate momentum q = (0.18, 0). The experimental results thus indicate this crossing should happen at a smaller momentum [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 and t 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 dopingdependence [43].
Second, the slope of the momentum dependence of the integrated experimental intensities is smaller compared with the calculated S(q). The slope extrapolates to a non-zero residual intensity at the Brillouin zone center 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 [7,34,53,54]. While the 2D antiferromagnetic acoustic spin wave has a zero structure factor at the zone center [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 intralayer one, the two branches quickly merge as q moves away from the zone center ( 0.1 r.l.u.) [7,53,54], 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.
The agreement between the t-t -t -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.  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 -J model. In this case, the calculated S(q) in the (π, π) direction shows the same qualitative behavior as in the experimental case. However, the spin structure factor S(q) in (π, 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 in the t-t -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 (π, 0) direction is only achieved after including the third neighbor hopping t as shown in Fig. 6c. However, it appears at relatively larger 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).
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(q) upon doping. In Fig. 7, we show the main results of this analysis on the t-t -t -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 × 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(q) calculated by including only up to third-neighbor spin-spin correlations. Since the short-range correlations mostly account for the large-q dynamics, they become poor in sketching the small-q properties, which leaves an artificial gap around q = (0, 0). Except the small-q region, where the gap appears, the results show an ascending intensity as a function of doping in the intermediate q range in both the (π, π) and (π, 0) directions and a crossover to descending sequence in larger q, which qualitatively agree with the spin structure factor S(q) calculated from the full Fourier transform depicted in Fig. 6c. We also notice that almost all of the spectral weight of S(q) at (π, π) 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 (π, π) 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 -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 (π, π) 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 cor-relations 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 (π, 0) and (π, π) directions. However, the intensity anisotropy between the (π, 0) and (π, π) is missing. This might signal the onset of the bimagnons in the charge channel of the RIXS spectrum [31,44] 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 [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 (π, 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 (π, π) 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 q = (π, π) 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 q = (π, π) paramagnons are largely a result of the short-range magnetic correlations. This is because the spectral weight in the magnetic response close to the q = (π, π) 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 q = (π, π) is known to be hugely sensitive to the longer-range magnetic correlations (see Supplementary  Figure 3). Since the nature of the (π, π) 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 and t ) 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.

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 ∼ 100 meV full width half maximum (FWHM). The single-crystal samples of Bi 2 Sr 2 CaCu 2 O 8+δ 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 ∼ 24 K. The momentum transfer 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 3d 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 q. Two linear polarizations (σ and π) 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 θ 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 q, so that the projections of 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 ±5% of the resolution, which is about ±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 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 (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 -J, t-t -t -J). We extrapolate it by computing the dispersion δ = n in − n out where n in is the averaged electron density taken over the sites which do not belong to the edges and n out is the averaged electron density taken over the sites which form the edge of the cluster. We computed the dispersion δ for different values of the edge factor λ for each doping level and model and take the final value of the edge factor λ as that at which the dispersion δ = 0. Nevertheless, the obtained values (λ ∼ 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 /N = 10 −6 . We make sure that the local density in the system is isotropic by using the edge factor λ as described above and we compute the real space spin-spin correlations S i ·S j for all pairs (i, j) labelling the system's sites. This way, we can compute the spin static factor S(q) as: This method is only valid when the local z-component of the spins is small enough, e.g. S z i ≤ 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 (q) as: Due to the symmetries of the considered models, S(q) = 3S z (q). Equation 6 was used for the doping level n = 0.22 for both the t-t -J and the t-t -t -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. CRSII2 141962). The theoretical work is supported by the Narodowe Centrum Nauki (

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 µ(E, ) of the studied sample, which depends on the energy (E) of the X-rays as well as the orientation of the polarization vector ( ) 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, (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 σ polarized X-rays, while the out-of-plane element f cc (E) is probed at grazing incidence and π polarization. As the fully grazing incidence is not possible, here we use 20 • incident angle, which gives a result of sin 2 20 • ·f aa (E) + cos 2 20 • ·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 λ·sinα, where λ=1/µ and α is the incident angle [73][74][75]. Although L λ is usually true in cuprates [75], sinα will reduce the photon penetration depth when α is small at grazing incidence, which will reduce the proportionality especially at the resonance where µ is large. Such a saturation effect can be described by the following relation [73][74][75]: .
(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/[λ(E, )·sin α] ∼ 1 and (2) the λ(E, ) 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 σ 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 TEY·sin α at the resonance peak and pre-edge background. The solid lines are the fitting by a function of A/(1 + 1/(r · sin α)), in which r gives an estimate of λ/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∼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 σ polarization as the inplane absorption tensor element f aa (E). For the TEY at α = 20 • and π incidence, some saturation effects exist. However, as the resonance peak amplitude is relatively small (∼ 15%) compared to the pre-edge background, the saturation factor is much less energy dependent compared to the σ polarization. Thus, we can use an overall energy-independent correction factor determined by the change of the pre-edge background from α = 90 • to α = 20 • at σ polarization to correct the spectra. By subtracting the contribution of sin 2 20 • ·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 3d 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 (∼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 λ/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.
(S3) Here the I exp t (E, φ, i )is the experimental (exp) RIXS intensity for a certain type (t) of excitation which can be spin, charge or d-d excitations, and I t (E, φ, i , f ) is the intrinsic RIXS intensity with E = E f − E i the energy loss. The denominator describes the self-absorption effect, where µ i(f) is the absorption coefficient of the incident (emitting) photons, k i(f) is the unit vector of the photon propagation direction, and n ( φ) is the unit vector normal to the sample surface. As the final polarization ( f ) is not resolved in this experiment, the measured intensity is a sum of all the possible f (σ and π polarizations). The azimuthal φ dependence of I t (E, φ, i , 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 φ dependence of f I t (E, φ, i , f ) for the non-spin-flip and spin-flip excitations. With the knowledge of the absorption tensor of our samples, µ i(f) (E i(f) , φ, i(f) ) can be determined precisely, and we can therefore calculate the self-absorption effect and include it in the final azimuthal dependence I exp t (E, φ, 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 t I exp t (E, φ, i ), and I exp t (E, φ, i ) can be expressed as a product of the intrinsic spectral weight of the excitation w t (E) and the azimuthal dependent geometry factor A t (E, φ, i ). In the low-energy part of the spectra in cuprates, with the spin-flip and non-spin-flip excitations related to 3d 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  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(q), see Eq. 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(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 and t is minimal when the cut-off in the Fourier transform is on the nearestneighbor 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 q in S(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(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(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 (π, π) is recovered, but the same behavior happens at (π, 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(q). Now the two directions in the Brilloin zone, (π, π) and (π, 0), show different behaviors. The damping of the (π, π) 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(q) are already present for the t-J model. A distinct behavior is now visible when longer range hoppings are included: in the (π, 0) direction the different doping lines do not cross nor overlap, while they do so once the "bare" t-J model is considered. In the (π, π) direction, both the t-t -J and t-t -t -J models already show the onset of a line crossing around (π/2, π/2). However, as previously stated, no crossing is present in the (π, 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-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 (π, π) peak is now doubled compared to the one seen when only firstneighbor 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-t -J case: in the (π, π) direction, the intensity of the peak has increased and it is comparable with that of the full S(q). Furthermore, the crossing at (π/2, π/2) is present. If we now look at the (π, 0) direction, we see the onset of a crossing close to (π, 0). Compared to the full S(q) case [see Fig. 6b], the 0.11 doping level is now "meeting" the other doping level lines at (π, 0). The last case is the one which includes also the t hopping as already shown in the main text. The behavior in the (π, π) direction is similar to that seen in the t-t -J model just discussed, therefore we will only examine what happens in the (π, 0) direction. Here the crossing close to (π, 0) is clearly visible and includes all three doping levels, as seen in experiments and in the full S(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 longrange distances in real spaces, this translates to small q values in q-space. Therefore, the "wrong" behavior seen at small 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(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(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 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(q) when further neighbors are included in the Fourier transform, we observe that the full spin structure factor S(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(q). This is because the spectral weight of the peak at (π, π) is still only at about half of its value when up to third neighbors are included.

Supplementary Note 4: Effect of larger t on the spin static structure factor
To check the importance of long-range hopping t in our model, we consider again the t-t -J model, but we set t = −0.5t rather than t = −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 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 (π, π) direction is lost. Indeed, the 0.22 hole-doping line loses intensity at a rather small 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 (π, 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 = −0.3t. However, the disagreement at q = (0. 18,0) persists also in this model. Hence, when compared to the option of including longer range hopping t , 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 q without compromising the agreement in the (π, π) direction.

Supplementary Note 5: Static charge structure factor
We calculate the static charge structure factor N (q) = n(−q)n(q) = 1 N i,j ( n i n j − n i n j )e iq(ri−rj ) .
(S4) for the t-t -t -J model on a 6 × 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(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 (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 [31,44]) 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 λ as described in the main text. To prove that this method leads to correct results, we calculated the static spin structure factor S(q) on the half-filled 6 × 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 (π, π)].
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π/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π/3 (at q = 2π/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(q) calculated using the Fourier transform with keeping up to third neighbor spin-spin correlations for the t-t -t -J modelwith 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. 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(q) obtained using DMRG on a 6 × 6 cluster (see text for further details) for the t-J (top panels), the t-t -J (middle panels), and t-t -t -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(q) for the nodal and anti-nodal directions of the Brillouin zone.   Figure 3. Theoretical static spin structure factor at half-filling. Static spin structure factor S(q) obtained using DMRG on a 6 × 6 cluster for the half-filled t-J model and (as in all other cases) including the edge factor λ 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(q) obtained using DMRG on a 6 × 6 cluster for the half-filled t-J model and (as in all other cases) including the edge factor λ for a square lattice. Note the enlarged momentum coverage w.r.t.