Particle-hole asymmetry in the dynamical spin and charge responses of corner-shared 1D cuprates

Although many experiments imply that oxygen orbitals play an essential role in the high-temperature superconducting cuprates, their precise role in collective spin and charge excitations and superconductivity is not yet fully understood. Here, we study the doping-dependent dynamical spin and charge structure factors of single and multi-orbital (pd) models for doped one-dimensional corner-shared spin-chain cuprates using several numerically exact methods. In doing so, we determine the orbital composition of the collective spin and charge excitations of cuprates, with important implications for our understanding of these materials. For example, we observe a particle-hole asymmetry in the orbital-resolved charge excitations, which is directly relevant to resonant inelastic x-ray scattering experiments and not captured by the single-band Hubbard model. Our results imply that one must explicitly include the oxygen degrees of freedom in order to fully understand some experimental observations on cuprate materials. The underlying mechanism of the superconductivity in high-Tc cuprates and the role spin and charge excitations play has been under debate since the cuprates were first discovered. Here, the authors perform a series of calculations for multi- and single-orbital model of a 1D cuprate system revealing the importance of oxygen degrees of freedom and their relevance to resonant inelastic x-ray scattering experiments.

T he high-temperature (high-T c ) superconducting cuprates are governed by multiple intertwined orders, including unconventional superconductivity, pseudogap behavior, and various spin-and charge-orders 1 . Determining how well the single-band Hubbard model describes these orders and their associated dynamical fluctuations has become a frontier problem in condensed matter physics [2][3][4] . Nevertheless, significant progress has been made towards understanding the physics of the singleband Hubbard model, including its superconducting [3][4][5][6][7] and normal-state transport properties [8][9][10] , its pseudogap [11][12][13] , and stripe orders [14][15][16][17] , and its dynamical response functions 11,[18][19][20] . While it is now clear that the predictions of the single-band Hubbard model are consistent with many experimental observations in the cuprates, recent x-ray scattering and nuclear magnetic resonance experiments have implied that the oxygen 2p orbitals provide important contributions to the spin-and chargeorders [21][22][23] . These observations raise questions on the validity of the single-band Hubbard model for describing some properties of the cuprates, including the fluctuations of their intertwined orders.
Studying the ground and excited-state properties of the singleand multi-band Hubbard model on large two-dimensional (2D) lattices is challenging. For example, theoretical studies using dynamical mean-field theory and its cluster and diagrammatic extensions [24][25][26][27] have been limited to relatively small clusters. The study of the Hubbard model can be extended to larger lattices with determinant quantum Monte Carlo (DQMC) 9,17,18,28,29 but this method is more limited by the sign problem in comparison to embedded cluster methods. Some zero temperature techniques, such as the tensor network 30,31 and the path-constrained auxiliary-field quantum Monte Carlo 32 can access large lattices but have difficulty in calculating the single-and two-particle dynamical response functions probed by angle-resolved photoemission (ARPES), inelastic neutron scattering (INS), and resonant inelastic x-ray scattering (RIXS) experiments. Importantly, contrasting results obtained from the single-band model with multiorbital models adds another layer of difficulty since the inclusion of the O 2p orbitals significantly increases the complexity of the problem 29,[33][34][35] .
In this work, we study single-and multi-orbital models for quasi-one-dimensional (1D) corner-shared spin-chain cuprates like Sr 2 CuO 3 and the recently synthesized doped chains Ba 2−x Sr x CuO 3+δ 36,37 . Specifically, we compute and contrast the momentum-resolved dynamical spin and charge responses in these models using DQMC, density matrix renormalization group (DMRG), and exact diagonalization (ED, see Supplementary Note 1). By focusing on 1D cuprate models, we can perform reliable calculations for the single-and two-particle response functions for large system sizes with good momentum resolution, even in the multi-orbital case. These calculations are enabled by the fact that 1D systems generally have manageable fermion sign problems 38 and algorithmic advances in computing dynamical response functions using DMRG 39,40 .
We consider a four-orbital pd-model for corner-shared cuprates, which includes the Cu 3d x 2 Ày 2 and O 2p x/y orbitals near the Fermi level, (Fig. 1a). The Hamiltonian for the fourorbital pd-model, written in hole language, is given by Here, 〈 ⋯ 〉 denotes a sum over nearest neighbor orbitals; d y i;σ and p y j;γ;σ create a hole with spin σ ( = ↑, ↓) on the ith Cu 3d x 2 Ày 2 orbital and the jth O 2p γ (γ = x, ±y) orbital, respectively; ϵ d and ϵ p,γ are the onsite energies;n d i;σ andn p j;γ;σ are the number operators for the Cu 3d x 2 Ày 2 orbital and O 2p γ orbital, respectively; t ij pd and t jj 0 pp are the nearest-neighbor Cu−O and O−O hopping integrals, whose phase factors are drawn in Fig. 1a; U d and U p are the onsite Hubbard interactions on the Cu and O orbitals, respectively, and U pd is the nearest-neighbor Cu−O Coulomb repulsion; Finally, μ is the chemical potential, which is adjusted to control the hole density in our DQMC simulations. Throughout, we adopt parameters determined from density functional theory calculations and comparisons to experiments 53,54 . Specifically, we set (in units of eV) ϵ d = 0, ϵ p,x = 3, ϵ p,y = 3.5, |t (p,x)d | = 1.5, |t (p,y)d | = 1.8, |t pp | = 0.75, U d = 8, U p = 4, and U pd = 1. For these parameters, 60% of the holes are distributed on the Cu orbitals and 40% of the holes reside on the O orbitals at half filling (see Supplementary Note 2).
To isolate the influence of the oxygen degrees of freedom and assess the validity of a single-band effective model, we also studied a single-band t-t 0 Hubbard model with model parameters selected to reproduce results from our multi-orbital model (see Supplementary Note 3). To remain consistent with the pd-model, our single-band Hubbard model is also written in the hole language, and is given by Here, t i,j = t and t 0 are the nearest-and next-nearest-neighbor hopping integrals (we set all longer range hopping to zero). We adopt t = 0.5 eV based on the analysis in ref. 55. To determine the remaining parameters, we adjusted t 0 and U to fit the dynamical magnetic susceptibility of the pd-model and found that U = 5.32t = 2.66 eV and t 0 ¼ 0:06t ¼ 0:03 eV produce the best agreement between these two models (see Supplementary Note 3). All of our calculations are performed on N = 20 chains (for a total of 80 orbitals total in the multi-orbital case). We have checked that finite-size effects are minimal for chains of this length (see Supplementary Note 4). The simulation temperature was held at T = 0.0625 eV for both the single-and multi-orbital DQMC calculations. We also perform simulations for the same system size at T = 0 using the DMRG.
Electronic structure. There is a significant orbital overlap between the Cu and O orbitals in quasi-1D and 2D cuprate materials, which hybridizes the Cu 3d x 2 Ày 2 and O 2p x,y orbitals. In the non-interacting limit, these orbitals form bonding (pd), nonbonding (pd) 0 , and anti-bonding (pd) * bands, as shown in Fig. 1b−d. In the corner-shared cuprates, there is also a flat band, which originates from a nonbonding combination of 2p ±y orbitals that does not hybridize with the 3d x 2 Ày 2 orbital.
Throughout, we work in hole language where hni ¼ 1 corresponds to half-filling (i.e., 1 hole/Cu), and hni > 1 (<1) corresponds to hole (electron) doping. At half-filling, when the interactions are turned on, the bonding (pd) band is split into the lower and upper Hubbard band (LHB/UHB) as well as an additional Zhang-Rice band (ZRS) located above the Fermi level 56,57 . These bands can be easily resolved in the total interacting density of states (DOS) shown in Fig. 1e, and in both the single-particle spectral functions shown in Figs. 2 and 3. For our parameters, the energy of the upper Hubbard band is close to the O-derived bands, leading to a broad peak in the DOS centered around ω = 6 eV.
To study the electronic structure of the interacting pd-model, we plot the momentum-dependent spectral functions A(k, ω) in For simplicity, we begin by discussing the spectra at half-filling. Both DQMC and DMRG have a clear gap at the Fermi level E F , consistent with a Mott-insulating state. In the low energy region [−4,4] eV, the DMRG spectra also have footprints of spin-charge separation, consistent with ARPES measurements on SrCuO 2 58 . Besides, the DMRG results show a dispersing band between 6 and 8 eV with a significant amount of O 2p character and two flat bands at ω = 6 eV and ω = 7.2 eV. The flat band at lower energy is composed almost entirely of the O 2p ±y orbitals and corresponds to the flat band shown in the noninteracting band structure. The higher-energy flat band is mixed between Cu and O and corresponds to the UHB with additional weak Cu satellites at ω ≈ 10 eV. Note that the contribution of the Cu orbitals to the anti-bonding (pd) * band is reduced significantly in the interacting case compared to the noninteracting case. Due to a combination of thermal broadening and the Maximum Entropy method, these fine structures blend into a broad peak in the DQMC results in the same energy region. By considering the full pd-model, we can access the Cu and O components to the spinon and holon states in the low energy region, which have not been reported in the literature to our knowledge. Our DMRG results show that the Cu and O weights of the main structures are comparable over the entire Brillouin zone. We note, however, that the holon-shadow bands are dominated by the Cu component below the Fermi level near k = 0 but dominated by the O component above the Fermi level near k = ±π/a. Our DQMC results show a similar composition for the main structures, but the intensity of the shadow bands is too weak to be captured by our Maximum Entropy method.
The system undergoes a metal−insulator transition as the system is either electron-or hole-doped. In our pd-model, we see that spinon-antiholon branches are responsible for the spectral weight crossing the Fermi level, consistent with the results of the single-band Hubbard model 59 . The orbital components in both low and high-energy regions do not change much compared to the half-filled case.
The magnetic excitations. We now examine the collective magnetic excitations of the pd-model.   At low-temperature and hni ¼ 1, the corner-shared cuprates are charge-transfer insulators with antiferromagnetic correlations.  The system's elementary magnetic excitations are spinons in this limit, which must be created in pairs. The magnetic excitation spectrum is dominated by a two-spinon continuum, which has been explicitly observed in, for example, Sr 2 CuO 3 using INS 43 and RIXS 45,46 . The upper and lower boundaries of the twospinon continuum are given by ω þ ðqÞ ¼ πJj sinðqa=2Þj and ω À ðqÞ ¼ πJ 2 j sinðqaÞj, respectively, where J is the Cu−Cu antiferromagnetic superexchange energy. Our DQMC and DMRG results for both the pd- (Figs. 4a and 5a) and single-band (Figs. 4f and 5f) models reproduce this behavior. Specifically, we observe a continuum of magnetic excitations confined within ω ± (q) but with a maximum spectral weight ω m (q) [indicated by the cyan points] concentrated at energies near the lower boundary ω − (q). This distribution suggests that the low-energy magnetic excitations of the pd-model are close to the Heisenberg limit, where the low-energy physics is described by the 1D t-J model 55,60 . By assuming that the locations of the maximum spectral weight in Figs. 4a and 5a correspond to ω − (q), we can estimate J ¼ 2 π ω m = sinð π 2 Þ % 350 meV for the multi-orbital model. This value is a little smaller than the value J ≈ 4t 2 /U = 376 meV that one would obtain from the single-band Hubbard model in the large U limit. (For this estimate, we have neglected the presence of t 0 ¼ 0:06t, which would induce a negligible frustration to the system.) Upon doping hni ¼ 1 þ x, the location of the zero-energy mode q s shifts from π a to ð1 À jxjÞ π a , consistent with the prior demonstration that the deviation of the wave vector from π/a is proportional to the excess hole or electron concentration 60 . We note that using DQMC, at the largest doping, |x| = 0.2, the spin excitations at q s ¼ 0:8 π a appear to acquire a finite gap at T = 0.0625 eV but remain gapless at zero temperature. This behavior may be an artifact of the Maximum Entropy Method and/or the finite temperature effect of the simulation.
We observe a hardening of the spin excitation spectrum as the total hole concentration increases. For example, Fig. 6a shows the change of the energy of the maximum intensity I max ðqÞ of S(q, ω) as a function of the hole density at q = 0.4π/a. Here, the solid triangles and circles represent DQMC results, while the open triangles and circles represent DMRG results. Since the DQMC data can be broad due to thermal broadening and the use of the Maximum Entropy method, we also provide approximate error bars, which are estimated as the energy range over which Sðq; ωÞ ≥ 0:99I max ðqÞ. The two dashed lines are guides to the eye. For both the pd-and downfolded single-band models, we find that the energy of the peak increases with hni. This observation implies that the spin excitation energy hardens with hole-doping and softens with electron doping. We also note that the spin excitations obtained with DQMC at finite temperature are slightly higher in energy compared to the DMRG results, implying that the spin excitations shift to higher energies as the temperature increases.
The observed particle−hole asymmetry in the energy scale of the magnetic excitations is opposite to what is found in the twodimensional Hubbard model and experimental observations for the 2D cuprates 18,20 . Since we observe consistent behavior in our 1D single-and multi-orbital models, we attribute this difference to the value of t 0 used in the models. For example, ref. 18 showed that the single-band Hubbard model could account for the observed behavior in 2D, provided t 0 % À0:3t. Similarly, ref. 61 reported that the spin excitations in a 1D single-band Hubbard model harden in the electron-doped case and soften significantly in the hole-doped case after assuming t 0 ¼ À0:3t to remain consistent with the 2D cuprates. Later, ref. 62 showed that changing the sign of t 0 reverses this behavior in 1D, consistent with our current results. It is important to stress that the nextnearest-neighbor hopping process t 0 in the 1D case is actually analogous to the next-next-nearest-neighbor process t ″ in the 2D cuprates. With this in mind, it is then natural to expect that the value of t 0 in 1D will be very different from the one found in 2D with jt 0 1D j % jt 00 2D j < jt 0 2D j. We circumvent this issue by comparing the single-band description directly to the multi-orbital description, where explicitly determining the value of t 0 is not needed. To reproduce the multi-orbital behavior, we found that t 0 ¼ 0:06t > 0, in line with our expectations that jt 0 1D j < jt 0 2D j. In the strong coupling limit, the magnetic moment m = S(S + 1) = 0.75 in the single-band Hubbard model; however, this value will be reduced for finite U due to double occupancy and additional covalency effects in the pd-model 43 . To determine by how much, we calculated the local moment from the expectation value of the local spin operator m ¼ hðŜ total i Þ 2 i, whereŜ total i ¼Ŝ i andŜ total i ¼ ∑ αŜi;α measures the total spin in the unit cell for the single and multi-orbital models, respectively. The results are shown in Fig. 6b, c for the single-and multi-band cases, respectively. For both models, the DMRG and DQMC results are in excellent agreement and predict an average local moment m ≈ 0.6 at half-filling, which decreases as the system is doped with electrons or holes. This decrease is mostly symmetric in the single-band case and slightly asymmetric in the multi-band case. Besides, we note that the temperature used in our DQMC simulations is low enough so that the local magnetic moment is the same as that at zero temperature. These results imply that the effective moment in both models has a slight reduction due to double occupancy, consistent with previous INS studies 43,44,63 .
We can also estimate the effective size of the system's local moment using a sum rule that relates it to an integral over the spin structure factor m ¼ lim Here, N is the number of the unit cells, α, β are orbital indices, r x,α and r x,β are the x-components of the basis vectors, and the factor of three comes from the sum over the three spin components, which contribute equally because of the unbroken SU (2) symmetry. Since INS experiments typically only access the low energy region, in practice, one needs to cut off the integration to a finite value of ω c (= 3 eV in our case). Performing the integral for the single-band case, we find that the total sum rule is almost completely recovered. However, for the pd-model, this energy window only recovers ≈ 75−80% of the sum rule, depending on the doping. This discrepancy is not due to the use of the Maximum Entropy method, as we obtain very consistent results using DMRG and DQMC. Rather, it reflects that some of the magnetic moment has been transferred to higher energies in the pd-model. To verify this, we also computed S(q, ω) to high energies and confirmed that extending to ω c = 13 eV recovers the missing weight (see Supplementary Note 7). Reference 43 previously estimated the local moment in Sr 2 CuO 3 by integrating INS data for S(q, w) up to ω c = 1 eV. The authors were only able to recover 80% of the sum rule after accounting for covalency effects and attributed the missing weight to Debye−Waller effects. Our results demonstrate that the missing weight has instead been transferred to higher energies. Remarkably, the magnetic excitations of the full multi-orbital are well reproduced by our effective single-band model, even though ≈20% of the overall magnetic moment is transferred to higher energies in the pd-model. Next, we analyze the contribution of each orbital to the total magnetic moment. The inset of Fig. 6c plots the weight of spin excitations between two Cu orbitals, two O orbitals, and Cu and O orbitals at half-filling. Here, the intraorbital spin excitations between neighboring Cu and O orbitals are labeled as "Cu" and "O," respectively, while interorbital spin excitations between Cu and O are labeled as "Cu−O". Interestingly, we observe that the spin excitations on the Cu sites have the maximum weight (about 52%), and interorbital spin excitations between Cu and O account for sizable 39% of the total excitations.
The charge excitations. We now examine the charge excitations of the pd-model and compare them to the excitations predicted by the single-band model. Figure 7 plots DQMC results for the dynamical charge structure factor N(q, ω) obtained from both models, again at T = 0.0625 eV and for different hole densities. As was the case with our spin results, N(q, ω) represents the total charge response, and panels a−e show the excitation spectrum of the pd-model while panels f−j show the excitation spectrum of the single-band Hubbard model. Figure 8 plots DMRG results for the same models at zero temperature, following the same format.
The charge excitation spectrum of the pd-model can be divided into low-and high-energy sectors, with dividing line occurring at ω ≈ 5 eV, as shown in Figs. 7 and 8. The high-energy region corresponds to particle−hole excitations from the low-energy bands crossing E F to the high-energy oxygen-derived bands and UHB observed in the spectral function. The charge excitations in the low-energy region appear as a sharp cosine-like excitation that is gapped at half-filling and gapless in the doped systems. These low-energy features originate from scattering within holon and ZRS bands near the Fermi level.
At zero temperature, shown in Fig. 8, additional fine structure in the high-and low-energy regions of the spectral function develops. Here, we observe that the low-energy region consists of two distinct branches, one gapped and the other gapless. The gapless excitations are intraband scattering within the holon branch of the doped system and are notably absent in the spectra at half-filling. The gapped excitations are then scattering from the holon band to the remnant of the ZRS band (located at ω ∈ [−3,−1] eV in Fig. 3a and ω ∈ [1, 4] eV in Fig. 3c,  respectively).
When comparing our pd-model results to the single-band model, we focus on the low-energy section of the charge excitation spectrum because the high-energy excitations are understandably absent in the single-band model. Overall, we find that the low-energy charge excitations of the pd-model are qualitatively well described by the single-band model. The DMRG results show that both the single-band and pd-models have gapped and gapless charge excitations when hni ≠ 1. The gapped excitations of the single-band model originate from the scattering between the holon band and the upper Hubbard band, while these gapped excitations of the pd-model come from the scattering between the holon band and the remnant ZRS band. In the DQMC results, due to the broadening of the finite temperature and the Maximum Entropy method, the sharp gapped spectrum is replaced by a broad spectrum, connecting to the spectrum of the gapless excitation.
To better visualize the charge excitations of the DQMC results at low-and high-energy, we plot N(q, ω) of the pd-model in Fig. 9 for hni ¼ 0:9, hni ¼ 1:1, hni ¼ 0:8, and hni ¼ 1:2 at q = π/a. We also include the DMRG results for reference. We decomposed the DQMC results into three Gaussian functions to distinguish the low-and high-energy charge excitations. The centers of these three Gaussian functions coincide with the peak position of the DMRG results. The sum of these three Gaussian functions matches the original DQMC results very well. Besides, we observe an asymmetry in the intensity of the lowest energy peak between hole-and electron-doped regimes consistent with the asymmetric orbital content between Cu and O in the undoped ground state.
Even though the single-band model can describe the lowenergy total charge excitations, the orbital-resolved results show a fascinating behavior, which could help account for the particlehole asymmetry observed in RIXS experiments 64,65 . Figure 10 shows the orbital-resolved N γ;γ 0 ðq; ωÞ, which is evaluated from both DMRG and DQMC calculations. Red, cyan, and blue colors in Fig. 10

Conclusion
We have studied the dynamical spin and charge structure factors of a four-orbital pd-model relevant for one-dimensional cuprate spin chains using numerically exact DQMC, DMRG, and ED. We also compared the two-particle response functions of the pdmodel against those predicted by an effective downfolded singleband Hubbard model.
Our results show that the single-band Hubbard model can describe the low-energy total spin and charge excitations of the pd-model upon hole-or electron-doping. In the case of the magnetic excitations, we found that the collective excitations harden (soften) with hole (electron) doping. This asymmetry is captured by the single-band Hubbard model with a small positive hopping between next-nearest-neighbors. We find even richer physics in our orbital-resolved results in the four-orbital pdmodel. For example, we find that the low-energy spin excitations mainly consist of intraorbital Cu−Cu and interorbital Cu−O  Our work has important implications for numerical studies of competing and intertwined orders in one-and two-dimensional cuprates. For example, in the 2D superconducting cuprates, a clear particle-hole asymmetry has been identified in the excitations probed by RIXS experiments 64,65 . Here, a branch of collective modes has been observed in the electron-doped Nd 2−x Ce x CuO 4 , which has been challenging to observe in holedoped cuprates. This particle-hole asymmetry of the collective modes is consistent with the behaviors of the charge excitations on the Cu site observed here, suggesting that the explicit inclusion of the oxygen degrees of freedom may be required to capture this physics. It would, therefore, be interesting to extend this study to higher dimensions by considering multi-leg ladders on the route towards full 2D models. Finally, our results provide extensive predictions for the evolution of the spin and charge excitations in doped cuprate spin chains, which have been synthesized recently for a wide range of hole doping 37 .

Methods
Determinant quantum Monte Carlo. The details of the DQMC algorithm applied to the multi-orbital Hubbard models can be found in ref. 66. DQMC works in the grand canonical ensemble, where the expectation value of an observableÔ is given by hÔi ¼ Z À1 Tr½Ôe ÀβH , where Z ¼ Tr½e ÀβH is the partition function and β is the inverse temperature.
To study the model's excited state properties, we measured the imaginary-time dynamical magnetic χ s and charge χ c susceptibilities. They are given by whereŜ z q ¼ ∑ i;γ e iqr i;γŜ z i;γ andn z q ¼ ∑ i;γ e iqr i;γn i;γ . Here, r i,γ represents the position of the orbital γ. Besides, we also calculate the spin and charge responses between Cu (O) and O sites, where the operator on the O site is given bŷ O ¼Ô p x þÔ p y þÔ p Ày .
To examine the spectral properties, we then used the method of the maximum of entropy 67 to analytically continue the imaginary-time susceptibilities to the real frequency axis. Some of the authors have used the same analytic continuation method to study the 2D cuprates, and reasonable results were obtained 20 . The dynamical spin and charge structure factors are calculated by the fluctuationdissipation theorem, which simplifies to S γ;γ 0 ðq; ωÞ ¼ Imχ γ;γ 0 s ðq; ωÞ 1 À e Àβω ð8Þ and N γ;γ 0 ðq; ωÞ ¼ Imχ γ;γ 0 c ðq; ωÞ The primary drawback to DQMC is the Fermion sign problem 68 , which limits the range of accessible temperatures and Hubbard interactions. In general, we have found that the sign problem is alleviated in 1D systems 38 , and the smallest value of the sign we obtained in our simulations is about 0.78, much larger than the sign value of the 2D three-orbital pd-model 29 .
Density matrix renormalization group. The DMRG 69,70 calculations were carried out with the correction-vector method 39 using the Krylov decomposition 40 , as implemented in the DMRG++ code 71 . This approach requires real-space representations for the dynamical structure factors in Eqs. (4) and (5), which can be found in ref. 72. Here, we calculated the response functions for N = 20 unit cell long chains with open boundary conditions, corresponding to total system sizes of N and 4N + 1 orbitals for the single-and multi-orbital cases, respectively. We kept up to m = 1000 DMRG states to maintain a truncation error below 10 −7 and introduced a spectral broadening in the correction-vector approach fixed at η = 0.1 eV for both the single-and multi-band calculations. We note that the accuracy of our calculation deteriorates at high excitation energies even with such small truncation errors. Panels a−d show the dynamical charge structure factor N(q, ω) as a function of energy ω at q = π/a for fillings of a hni ¼ 0:9, b hni ¼ 1:1, c hni ¼ 0:8, and dn ¼ 1:2, respectively. The DQMC spectra (blue line) have been fit with a set of Gaussian distributions, whose energies correspond well with the main peaks observed in the DMRG data (black line). The parameters for the pd-model, which includes the full Cu and O basis, are identical to those used in Figs. 4

and 5.
This behavior is expected within our Krylov DMRG correction vector approach, as introduced in ref. 40. In the Krylov method, another source of error is given by the Lanczos error, which occurs in the tridiagonalization of the Hamiltonian for the construction of the correction vector c:v: j i¼ 1 ωÀHþiηÔ g:s: , whereÔ is the relevant operator. In the current work, we have imposed a maximum Lanczos error of 10 −8 allowing at most 200 iterations (as opposed to standard 10 −12 Lanczos precision of ground state computations) to avoid the proliferation of a larger number of Lanczos steps at high excitation energies and the consequent need for CPU consuming reorthogonalization of the Lanczos vectors.

Data availability
The data that support the findings of this study are available from the corresponding authors on reasonable request.