Particle-hole asymmetry in the dynamical spin and charge structure factors of the corner-shared one-dimensional cuprates

The collective spin and charge excitations of doped cuprates and their relationship to superconductivity are not yet fully understood, particularly in the case of the charge excitations. Here, we study the doping-dependent dynamical spin and charge structure factors of single and multi-orbital models for the one-dimensional corner shared spin-chain cuprates using several numerically exact methods. We find that the singleband Hubbard model can describe the spin and charge excitations of the $pd$-model in the low-energy region, including the particle-hole asymmetry in the spin response. However, our results also reveal that the weight of the interorbital spin excitations between Cu and O orbitals is comparable to the weight of the spin excitations between two Cu orbitals. This finding elucidates the microscopic nature of the spin excitations in the 1D cuprates and sheds light on the spin properties of other oxides. Importantly, we find a particle-hole asymmetry in the orbital-resolved charge excitations, which cannot be described by the singleband Hubbard model and is relevant to resonant inelastic x-ray scattering experiments. Our results imply that the explicit inclusion of the oxygen degrees of freedom may be required to understand experimental observations.


Introduction
The high-temperature (high-T c ) superconducting cuprates are governed by competition between multiple intertwined orders including unconventional superconductivity, pseudogap behaviour, and various spin-and charge-orders 1 . Determining how well these competing orders and their associated dynamical fluctuations are described by the singleband Hubbard model 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-17 , and its dynamical response functions 11,[18][19][20] . While it is now clear that the predictions of the singleband Hubbard model are consistent with many experimental observations in the cuprates, recent x-ray scattering and nuclear magnetic resonance experiments have inferred that the oxygen 2p orbitals provide important contributions to the spin-and charge-orders 21-23 . These observations raise questions on the validity of the singleband 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 single-and multi-band Hubbard model on large twodimensional (2D) lattices is challenging. For example, theoretical studies using the dynamic 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-resolve photoemission (ARPES), inelastic neutron scattering (INS), and resonant inelastic x-ray scattering (RIXS) experiments. Importantly, contrasting results obtained from the singleband model with multi-orbital models adds another layer of difficulty, since the inclusion of the O 2p orbitals significantly increases the complexity of the problem 29, 33-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 . 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 are able to 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 37 and algorithmic advances in computing dynamical response functions using DMRG 38, 39 .

Corner-Shared Spin Chains
The active orbitals in spin-chain cuprates like Sr 2 CuO 3 and Ba 2−x Sr x CuO 3+δ are located in their CuO 4 plaquettes, which are arranged in a corner-shared geometry, as shown in Fig. 1a. Sr 2 CuO 3 has been extensively studied both experimentally 40-45 and using singleband models like the Heisenberg model, the t − J model 40-43, 45, 46 , and the extended Hubbard model 47,48 . These models generally provide an excellent description of the magnetic excitations -in this case a multi-spinon continuum -observed in INS 40-43 and RIXS 44,45 experiments. The spin and charge excitations of doped samples have received comparatively less attention because doping electrons or holes in Sr 2 CuO 3 has proved to be challenging [49][50][51] . However, the Ba 2−x Sr x CuO 3+δ system offers new possibilities in this regard 36 .
To determine the influence of the orbital degrees of freedom on the spin and charge dynamics of the cuprates, 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). To isolate the influence of the oxygen degrees of freedom, and assess the validity of a singleband effective model, we also studied a singleband t-t Hubbard model with model parameters selected to reproduce results from our multi-orbital model (see Methods and Supplementary Note 2).

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 Figs. 1b -1d. In the corner-shared cuprates, there is also a flat band, which originates from a nonbonding combination of 2p ±y orbitals that do not hybridize with the 3d x 2 −y 2 orbital.
Throughout, we work in the hole-language where n = 1 corresponds to half-filling (i.e. 1 hole/Cu), and n > 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 52,53 . These bands can be easily resolved in both the single-particle spectral functions shown in Figs. 2b and 3b, and in the total interacting density of states (DOS) shown in Fig. 1e. 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 Figs. 2 and 3 for hole concentrations n = 1.1, 1, and 0.9. Figures 2 and 3 show results obtained from finite-temperature DQMC and zero temperature DMRG calculations, respectively. From bottom to top, the momentum k increases from −π/a to π/a in each panel. The black dashed line represent the Fermi level E F . To highlight the orbital content of the spectral features, we indicate the Cu d x 2 −y 2 orbital of the spectral function with red color and the sum of the spectral weights of the O p x and O p ±y components in cyan. The zero-temperature spectral function shown in Fig. 3 exhibits sharper features and hence richer details compare to the DQMC results.
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 54 . 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 use of 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 singleband Hubbard model 55 . The orbital components in both low and high energy regions do not change much compared to the half filling case.

The Magnetic Excitations
We now examine the collective magnetic excitations of the pd-model. Figure 4 summarizes our finite-temperature (T = 0.0625 eV) DQMC results for the dynamical spin structure factor S(q, ω) for doping levels spanning from n = 0.8 − 1.2. (Additional DQMC data for the doping levels not shown here are provided in Supplementary Note 3.) Figure 5 shows comparable results at zero temperature obtained using DMRG. In both Figs. 4 and 5, panels a-e show the spectra obtained from the full pd-model while panels f-j show the corresponding results obtained from a singleband Hubbard model. In the latter case, we use t = 0.5 eV t = 0.06t, U = 2.66 eV to best fit to the pd-model results (see Methods and Supplementary Note 2 for a more detailed discussion). To compare to the singleband model, we present the total spin response here. We stress that both the cluster sizes and simulation temperatures are the same for the two models, and one can directly compare the results shown on both sides of each figure.
At low-temperature and n = 1, the corner-shared cuprates are charge transfer insulators with antiferromagnetic correlations. In this limit, the system's elementary magnetic excitations are spinons, 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 42 and RIXS 44,45 . The upper and lower boundaries of the two-spinon continuum are given by ω + (q) = πJ| sin(qa/2)| and ω − (q) = πJ 2 | sin(qa)|, respectively, where J is the Cu-Cu antiferromagnetic superexchange energy. Our DQMC and DMRG results for both the pd-  Fig. 5f] models reflect this behaviour. 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 deviate only mildly from the Heisenberg limit, where the low-energy physics is described by the 1D t-J model 56,57 . By assuming that the locations of the maximum spectral weight in Fig. 4a and Fig. 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.06t, which would induce a negligible frustration to the system.) Upon doping n = 1 + x, the location of the zero-energy mode q s shifts from π a to (1 − |x|) π a , consistent with the prior demonstration that the deviation of the wave vector from π/a is proportional to hole and electron doping 56 . 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 behaviour 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 singleband models, we find that the energy of the peak increases with n . This observation implies that the spin excitation energy hardens with hole-doping and softens with electron doping. This asymmetric behaviour is opposite to what is observed in the two-dimensional Hubbard model and experimental observations for the 2D cuprates 18, 20 . Since we observe consistent behaviour in our 1D single-and multi-orbital models, we attribute this difference to the dimensionality of the respective systems. 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 3/20 higher energies as the temperature increases.
We can estimate the effective size of the system's total local moment m from a sum rule that relates it to an integral over the spin structure where N is the number of the unit cells, 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, we cut off the integration to ω c = 2 eV here. In the limit of strong interactions, the magnetic moment m = S(S + 1) = 0.75 in the singleband Hubbard model; however, this value will be reduced for finite U due to double occupancy and additional covalency effects in the pd-model 42 . To determine by how much, we evaluated Eq. (1) by integrating our numerical data over binding energies in the interval between 0 and ω c = 2 eV, as shown in Fig. 6b. We note that the DQMC results are a little larger than the DMRG results, which is attributed to a broadening effect of the Maximum Entropy method and finite temperature. The magnetic moment is about 0.536 at half-filling and zero temperature for the Hubbard model, and it decreases somewhat symmetrically as the system is doped with electrons or holes. These results imply that the effective singleband model has some nonzero reduction due to double occupancy, consistent with previous INS studies 42,43,58 . For the pd-model, we observe that the magnetic moment is about 0.17 smaller compared to the singleband model, both at zero and finite temperature. This result implies that the hybridization with oxygen further reduces the total magnetic moment in cuprate chains. Next, we analyze the contribution of each orbital to the total magnetic moment. Figures 6c and 6d plot the weight of spin excitations between two Cu orbitals, two O orbitals, and Cu and O orbitals at half filling calculated using the DQMC and DMRG data. Here, the intraorbital spin excitations between neighbouring 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 55%), and interorbital spin excitations between Cu and O account for 38% of the total excitations, which is nonnegligible.
It is remarkable that the magnetic excitations of the full multi-orbital pd-model are well reproduced by our effective singleband model, apart from the overall magnetic moment per unit cell predicted by the two models. This difference is also reflected in the overall intensity of the spin structure factors.

The Charge Excitations
We now examine the charge excitations of the pd-model and compare them to the excitations predicted by the singleband model. Fig. 7 plots DQMC results for the dynamic 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 singleband Hubbard model. Fig. 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 Fig. 7 and Fig. 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 the 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 singleband model, we focus on the low-energy section of the charge excitation spectrum because the high-energy excitations are understandably absent in the singleband model. Overall, we find that the low-energy charge excitations of the pd-model are qualitatively well described by the singleband model. The DMRG results show that both the singleband model and the pd-model have gapped and gapless charge excitations whenn = 1. The gapped excitations of the singleband model originates from the scattering between the holon band and the upper Hubbard band, while this 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 regions, we plot N(q, ω) of the pd-model in Fig. 9 for n = 0.9, n = 1.1, n = 0.8, and n = 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 center of these three Gaussian functions coincide with the peak position of the DMRG results. The summation 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 low-energy total charge excitations can be described by the singleband model, the orbitalresolved results show a very interesting behaviour, which could help account for the particle-hole asymmetry observed in RIXS experiments 59,60 . Figure. 10 shows the orbital-resolved N γ,γ (q, ω), which is evaluated from both DMRG and DQMC calculations. Red, cyan, and blue colors in Fig. 10

Discussion
We have studied the dynamic 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 pd-model against those predicted by an effective downfolded singleband Hubbard model.
Our results show that the singleband 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 singleband Hubbard model with a small positive hopping between next-nearest neighbors. We find even richer physics in our orbital-resolved results in the four band pd-model. For example, we find that the low-energy spin excitations mainly consist of intraorbital Cu-Cu and interorbital Cu-O components, and we distinguish the dynamical spin behaviors on each orbital, including the effects of electron-electron interactions beyond first principle approaches 42  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 59,60 . Here, a branch of collective modes has been observed in the electron-doped Nd 2−x Ce x CuO 4 , which has been difficult to observe in hole-doped cuprates. This particle-hole asymmetry of the collective modes is consistent with the behaviours of the charge excitations on the Cu site observed here, which suggests 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 route towards full 2D models.

Models
The Hamiltonian for the four-orbital pd-model describing corner-shared cuprates [see Fig. 1 (a)], written in hole language, is given by Here, · · · denotes a sum over nearest neighbor orbitals; d † i,σ and p † j,γ,σ creates a hole with spin σ (=↑, ↓) on the i th Cu 3d x 2 −y 2 orbital and the j th 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 i j pd and t j j pp are the nearest-neighbor Cu-O and O-O hopping integrals, whose phase factors are drawn in Fig. 1(a); 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.
We map the low-energy spin properties of the four-orbital pd-model to a singleband t − t Hubbard model. To remain consistent with the pd-model, our singleband Hubbard model is also written in the hole language, and is given by Here, t i, j = t and t 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. 57. To determine the remaining parameters, we adjusted t and U to fit the dynamic magnetic susceptibility of the pd-model and found that U = 2.66 eV and t = 0.06t produce the best agreement between these two models (see Supplementary Note 2).

Determinant Quantum Monte Carlo
The details of the DQMC algorithm applied to the multi-orbital Hubbard models can be found in Ref. 63 To study the model's excited state properties, we measured the imaginary-time dynamic magnetic χ s and charge χ c susceptibilities. They are given by and Here, γ (γ ) is the orbital index, andn γ q,σ andŜ γ,z q are the Fourier transforms of the local density and spin-z operators.
To compare to the singleband model, we calculate the total spin and charge responses, which are given by and χ c (q, τ) = n z q (τ)n z −q (0) , 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 byÔ =Ô p x +Ô p y +Ô p−y .
To examine the spectral properties, we then used the method of the maximum of entropy 64 to analytically continue the imaginary-time susceptibilities to the real frequency axis. The same analytic continuation method has been used by some of the 6/20 authors to study the 2D cuprates, and reasonable results were obtained 20 . The dynamical spin and charge structure factors are calculated by the fluctuation-dissipation theorem, which simplifies to S(q, ω) = Imχ s (q, ω) 1 − e −β ω (8) and The primary drawback to DQMC is the Fermion sign problem 65 , which limits the range of accessible temperatures and Hubbard interactions. In general, we have found that the sign problem is alleviated in 1D systems 37 , 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 .
All of our DQMC calculations are performed on N = 20 chains (for a total of 40 orbitals total in the multi-orbital case). The simulation temperature was held at T = 0.0625 eV for both the single-and multi-orbital calculations.

Density Matrix Renormalization Group
The DMRG 66, 67 calculations were carried out with the correction-vector method 38 using the Krylov decomposition 39 , as implemented in the DMRG++ code 68 . This approach requires real-space representations for the dynamical structure factors in Eqs. (4) and (5), which can be found in Ref. 69. Here, we calculated the response functions for N = 20 unit cell long chains and open boundary conditions, which corresponds 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.