Computing Resonant Inelastic X-Ray Scattering Spectra Using The Density Matrix Renormalization Group Method.

We present a method for computing the resonant inelastic x-ray scattering (RIXS) spectra in one-dimensional systems using the density matrix renormalization group (DMRG) method. By using DMRG to address this problem, we shift the computational bottleneck from the memory requirements associated with exact diagonalization (ED) calculations to the computational time associated with the DMRG algorithm. This approach is then used to obtain RIXS spectra on cluster sizes well beyond state-of-the-art ED techniques. Using this new procedure, we compute the low-energy magnetic excitations observed in Cu L-edge RIXS for the challenging corner shared CuO4 chains, both for large multi-orbital clusters and downfolded t-J chains. We are able to directly compare results obtained from both models defined in clusters with identical momentum resolution. In the strong coupling limit, we find that the downfolded t-J model captures the main features of the magnetic excitations probed by RIXS only after a uniform scaling of the spectra is made.

Introduction

Resonant inelastic x-ray scattering (RIXS) has emerged as a powerful and versatile probe of elementary excitations in quantum materials[1,2]. One of the most commonly used approaches for computing RIXS spectra is small cluster exact diagonalization (ED)[3-21]. This approach is limited by the exponential growth of the Hilbert space, however, which restricts clusters to a relatively small size, thus limiting momentum resolution. For example, ED treatments of multi-orbital spin-chain systems such as the edge-shared CuGeO3 or corner shared Sr2CuO3 have been limited to no more than six CuO4 plaquettes[3,6,8,12,22], while studies carried out using downfolded singleband Hubbard (or t-J) chains have been limited to ~16–22 sites[4,10,19,20].

The density matrix renormalization group (DMRG) is the most powerful method for computing the ground state properties of strongly correlated materials in one dimension (1D)[23-25]. Within the DMRG framework, several efficient methods are available for computing dynamical correlation functions, including: time-dependent DMRG[26,27], which computes dynamical correlation functions in the time domain with a subsequent Fourier transform into frequency space[28]; correction-vector methods, which compute the dynamical correlator directly in frequency space[29-32]; continued fraction methods[33-35]; and Chebyshev polynomial expansion methods[36,37]. In this work, we present an efficient algorithm to compute the dynamical correlation function representing the RIXS scattering cross section with DMRG directly in frequency space. We then apply this approach to computing the Cu L-edge RIXS spectra of a quasi-1D corner-shared cuprate (e.g., Sr2CuO3, see Fig. 1b), a geometry that is challenging for ED calculations due to significant finite size effects[3,6,8]. We consider a multi-orbital Hubbard model that retains the Cu and O orbital degrees of freedom, as well as a downfolded t-J model. Using our DMRG-based approach, we access systems sizes beyond those accessible to ED, thus enabling us to directly compare the results obtained from the two models on large clusters with comparable momentum resolution.

Figure 1

(a) A sketch of the algorithm for computing the real space Kramers-Heisenberg formula [Eq. (3)] using the DMRG method at a fixed value of the energy loss Ω = ω. (b) A sketch of the multi-orbital pd-model describing the corner shared spin chain cuprates (e.g. Sr2CuO3).

The Kramers-Heisenberg Formalism

In a RIXS experiment, photons with energy ωin and momentum kin (ℏ = 1) scatter inelastically off of a sample, transferring momentum q = kout − kin and energy Ω = ωout − ωin to its elementary excitations. The resonant nature of the probe arises because ωin is tuned to match one of the elemental absorption edges, such that it promotes a core electron to an unoccupied level of the crystal.

The intensity of the RIXS process I(q, Ω) is given by the Kramers-Heisenberg formalism[1,2], withHere, E and E are the energies of the ground |g〉 and final states |f〉 of the system, respectively. The scattering amplitude F is defined aswhere is the dipole transition operator describing the core-hole excitation. In what follows, we consider the Cu L-edge (a Cu 2p → 3d transition). In this case, the dipole operator is defined as , where adds an electron to the valence band orbital , and destroys a spin σ electron (creates a hole) in a core 2p orbital on site j located at R. The prefactor is the matrix element of the dipole transition between the core 2p orbital and the valence orbital, , which we set to 1 for simplicity. Γ is the inverse core-hole lifetime, and , where describes the Coulomb interaction between the core hole and the valence electrons and is the many-body Hamiltonian of the system.

Under the assumption that the core-hole is completely localized, and only one Cu 2p orbitals is involved in the RIXS process, Eq. (2) simplifies towhere we have defined the local dipole-transition operator and , with .

Reformulation of the Problem for DMRG

The primary difficulty in evaluating Eq. (1) lies in computing the final states |f〉. This task is often accomplished using ED on small clusters meant to approximate the infinite system. Obtaining these same final states is usually impossible with DMRG, which targets only the ground state; however, we will show that to accomplish this task one can use the Lanczos method, which projects the state onto a Krylov space[38]. Some of the present authors introduced this alternative method to calculate the correction vectors for frequency-dependent correlation functions with DMRG[32].

We can formulate an efficient DMRG algorithm by expanding the square in Eq. (1), yielding a real space version of the Kramers-Heisenberg formula. To compact the notation, we define vectors . Using this definition, Eq. (1) can be written asHere, η is a broadening parameter, which plays the same role as the Gaussian or Lorentzian broadening introduced in ED treatments of the energy-conserving δ-function appearing in Eq. (1). Throughout this work, we set it to 75 meV, which is in the range of the typical energy resolution currently attainable in RIXS experiments[39]. Note that the vectors |α〉 must be computed for each value of ωin and Γ.

The X-ray absorption spectrum (XAS) can be computed using a similar formalism. Its intensity is given by

Finally, we note that we have removed the elastic line from all spectra shown in this work. The precise method for doing this is discussed in Supplementary Note IV.

Computational Procedure

The algorithm to compute the RIXS spectra using Eq. (3) is as follows (see also Fig. 1a):

Step 1: Compute the ground state |g〉 of using the standard ground state DMRG method. The vector |g〉 must be stored for later use.

Step 2: Restart from the ground state calculation, reading and then targeting the ground state vector calculated earlier and using a different Hamiltonian , where j = c is the center site of the chain. Construct the vector |α〉 at the center of the chain using the Krylov-space correction vector approach[32] where we have performed a Lanczos tridiagonalization with starting vector , and a subsequent diagonalization of the Hamiltonian , represented in its diagonal form Dch, in the Krylov basis. The vector |α〉 should also be stored for later use. Because the cluster is not periodic, the use of a central site here represents an approximation that will become exact in the thermodynamic limit. This central site “trick” was used for the first time in the application of time-dependent DMRG[26].

Step 3: Restart from previous run, now using a different Hamiltonian . Read and then target (in the DMRG sense) the ground state vector |g〉, as well as the vector |α〉 constructed in Step 2. For each site j, except for the center site considered in Step 2, construct the vector with a Lanczos tridiagonalization with starting vector , and a subsequent diagonalization of . This step of the algorithm requires a number of runs which is equal to the number of sites minus 1, i.e., L − 1. These can be run in parallel on a standard cluster machine, restarting from Step 2. Performing Step 2 and Step 3 in this sequence is crucial for having the vectors |α〉 and |α〉 in the same DMRG basis. The vector |α〉 should also be stored for later use.

Step 4: Restart using the original Hamiltonian . Read and then target the ground state vector |g〉, the vector |α〉, as well as the vector |α〉 constructed in Step 3. For a fixed Ω = ω, compute the correction vector of |α〉 using again the Krylov-space correction vector approach as with a Lanczos tridiagonalization (using as the seed) and a subsequent diagonalization of the Hamiltonian , with D being the diagonal form of in the Krylov basis. This is a crucial part of the algorithm, which amounts to computing the correction vector of a previously calculated correction vector |α〉. Execute this computation NΩ times for .

Step 5: Finally, compute the RIXS spectrum in real space (in I(Ω) we omit the spin indices γ, γ′, σ, σ′ in order to lighten the notation) and then Fourier transform the imaginary part to obtain the RIXS intensity.

Computational Complexity

The computational cost required for DMRG to compute the RIXS spectrum can be easily estimated, assuming that the ground state of the Hamiltonian has already been calculated. Let C2–3 be the computational cost (i.e., the number of hours) for a single run in Step 2 (1 run only) or Step 3 (L − 1 runs in total). Let C4 be the computational cost for a single run in Step 4. The total computational time needed to compute the RIXS spectrum is then CPUcost = C2–3L + C4LNΩ, where NΩ is the number of frequencies needed in a given interval of energy losses. As explained in the previous section, we use a center site “trick” to reduce the computational cost by a factor of the order of L (Eqs (3–8)). For the largest system size considered in this work (20 plaquettes in the CuO4 multi-orbital model at half-filling, using up to m = 1000 DMRG states), the typical values for CPUcost on a single core of a standard computer cluster are: C2–3 ~ 2 hours, while C4 ~ 2 − 24 hours. The computational cost C4 for Step 4 follows the typical performance profile of the Krylov-space approach found in ref.[32], where less CPU time is needed to compute the spectra at lower energy-losses. We also note that the calculation of each energy loss is trivially parallelizable. From these assumptions, we estimate the proposed method can compute the RIXS spectrum of a cluster as large as Cu20O61 in less than a day if enough cores are available. In this specific case, one single core run was needed for ground state calculation, 80 single core runs for Steps 2–3, and 800 single core runs for Step 4.

Numerical Results for the t-J Model

We first apply our approach to compute the RIXS spectrum of the 1D t-J model as an effective model for the antiferromagnetic corner-shared spin chain cuprate Sr2CuO3 (see Methods). Throughout this paper, we adopt open boundary conditions, work at half-filling, and set t = 0.3 eV for the nearest neighbor hopping and J = 0.25 eV for the antiferromagnetic exchange interaction. These values are typical for Sr2CuO3[10,20,21,40-46].

Before scaling up our DMRG calculations to large systems, we benchmarked our method by directly comparing our DMRG results to ED. The results for a L = 16 sites t-J chain are presented in Supplementary Note I. (We provide a similar comparison for a four-plaquette multi-orbital cluster in Supplementary Note II.) Our DMRG approach gives perfect agreement with the ED result for both the XAS and RIXS spectra, for the largest clusters we can access with ED. All of the DRMG simulations presented in this work used up to m = 1000 states, with a truncation error smaller than 10−6.

We now turn to results obtained on a L = 64 site chain, as shown in Fig. 2. Here, we present results for the spin-conserving (ΔS = 0) and non-spin-conserving (ΔS = 1) contributions to the total RIXS intensity. The ΔS = 0 contribution corresponds to the σ = σ′ and γ = γ′ terms in the Kramers-Heisenberg formula Eq. (3). In this case, only two configurations ( and , ) have to be explicitly calculated with DMRG, as the other two possible spin-conserving configurations contribute equally by symmetry. The remaining terms with σ ≠ σ′ and γ ≠ γ′ determine the non-spin-conserving ΔS = 1 contributions to the spectrum. In this case, only one configuration (, , , ) has been simulated with DMRG, as the flipped configuration (, , , ) contributes equally by symmetry. The remaining two possible non-spin-conserving configurations also give zero contribution to the RIXS spectrum by symmetry.

Figure 2

DMRG results for the RIXS intensity I(q, Ω) of a half-filled t-J chain. Results are shown for a L = 64 site chain, in the (a) ΔS = 0 and (b) ΔS = 1 channels. The remaining parameters are t = 0.3 eV, J = 0.25 eV, η = 75 meV and ωin = 0.1 eV (which corresponds to the resonance observed in the XAS).

In Fig. 2, the ΔS = 1 part of the RIXS spectrum shows a continuum of excitations resembling the two spinon continuum commonly observed in the dynamical spin structure factor S(q, ω) of one-dimensional spin-1/2 antiferromagnets[47-50]. The ΔS = 0 contribution in Fig. 2a shows two broad arcs with maxima at q = π/2a. Notice also the perfect cancellation of the RIXS signal at the zone boundary, which is in open boundary conditions. Our results agree with the ED results of refs[4],[19], but with much better momentum resolution. We find that the finite size effects of the magnetic excitations in the t-J model are mild; we observe only small differences between results obtained on L = 16 (shown in Supplementary Note I) and L = 64 site clusters.

Magnetic Excitations in the Multi-orbital pd-Model

In the strong coupling limit, the low-energy magnetic response of the spin-chain cuprates are believed to be effectively described by a single orbital Hubbard or t-J model[51,52]. According to this picture, holes predominantly occupy the Cu orbitals at half-filling, while the oxygens along the Cu-Cu direction provide a pathway for superexchange interactions between the nearest-neighbor Cu orbitals. Since our DMRG approach provides access to large cluster sizes, we now compute the RIXS spectrum of a more realistic multi-orbital model. Here, we consider the challenging corner-shared geometry, which suffers from slow convergence in the cluster size. To address this, we consider finite 1D CuO3 clusters, with open boundary conditions, as illustrated in Fig. 1b for the n = 4 case. The Hamiltonian is given in the Methods section. We evaluated the Cu L-edge RIXS intensity for this model as a function of n for up to n = 20 CuO4 plaquettes.

The RIXS spectra for spin-conserving (ΔS = 0) and non-spin-conserving contributions (ΔS = 1) calculated with our DMRG method are shown in Fig. 3. Similar to the t − J spectra, panels (a–f) in Fig. 3 show two broad arcs with maxima at ±π/2a. Here, we observe significant finite size effects in the RIXS spectra. Some of these effects are the result of our use of the “center-site approximation” in evaluating the Kramers-Heisenberg formula. For example, the downward dispersing low-energy peak centered at q = 0 seen in the smaller clusters is the result of this approximation. These features in the spectra can be minimized by carrying out calculations on larger clusters. Because of this, to observe well defined spectral features, we need to consider at least fourteen plaquettes. The pd model also shows that the low-energy ΔS = 1 part of the RIXS spectrum is characterized by a two-spinon-like continuum of excitations (panels (g–l) in Fig. 3).

Figure 3

The response I(q, Ω) obtained with DMRG for a multi-orbital pd model as a function of the number of CuO4 plaquettes. Panels (a–f) show the spin-conserving ΔS = 0 channel of the RIXS intensity, while panels (g–l) show the non-spin-conserving ΔS = 1 channel. Results are shown for 8 to 20 unit cells at half filling, computed at resonance with ωin = 2.5, Γ = 0.2, and V = 4.0 (in eV units).

Comparing the Multi-Orbital and Effective t-J Models

Over the past decade, there has been a considerable research effort dedicated to quantitatively understand the intensity of magnetic excitations probed by inelastic neutron scattering (INS)[44,49,53]. This effort is motivated by the desire to understand the relationship between the spectral weight of the dynamical spin response S(q, ω) and the superconducting transition temperature T of unconventional superconductors[54]. To this end, several studies have set out to determine whether the observed INS intensity can be accounted for by the Heisenberg model in low-dimensional strongly correlated cuprates. In this context, the highest degree of success has been achieved in quasi-1D materials, where accurate theoretical predictions for S(q, ω) are available[44,49]. Many of these studies find that the low-energy Heisenberg model can indeed account for the INS intensity, after including corrections due to effects such as the degree of covalency, its impact on the form factor, and Debye-Waller factors.

RIXS has also been applied to study magnetic excitations in many of the same materials[10,43,46]. It is therefore natural to ponder how covalency modifies the magnetic excitations as viewed by RIXS. In the limit of a short core-hole lifetime, or under constraints in the incoming and outgoing photon polarization, the RIXS intensity for single orbital Hubbard and t-J chains is well approximated by S(q, ω)[1,5,19,46]. However, to the best of our knowledge, no systematic comparison of the RIXS intensity, as computed by the Kramers-Heisenberg formalism, has been carried out for multi-orbital and downfolded Hamiltonians.

Figure 3 demonstrated that DMRG grants access to large system sizes. We are, therefore, in a position to make such a comparison for the multi-orbital spin-chain cuprates. Figure 4 compares the spectra computed on a L = 20 site t-J chain against those computed on a Cu20O61 cluster, such that the momentum resolution of the two clusters is the same. The parameters for the multi-orbital model are identical to those used in Fig. 3. To facilitate a meaningful comparison with the t-J model, we adopted t = 0.5 eV and J = 0.325 eV. These values are obtained by diagonalizing a Cu2O7 cluster (see methods). Note that we use the same value of the core hole potential V = 4 eV in both cases. In Supplementary Note III, we show results for a reduced value of V for the t-J model, which are very similar. To compare the two spectra, the results for the t-J model have been scaled by a factor of 0.26 such that the maximum intensity of the ΔS = 1 excitations is the same at the zone boundary. This factor presumably accounts for covalent factors and differences in how the core-hole interacts with the distribution of electrons in the intermediate state.

Figure 4

A comparison of the magnetic RIXS excitations computed using DMRG for a 20-site t-J chain (solid red line) and multi-orbital pd model (dashed blue line) with twenty unit cells at half filling. Results are shown for the (a) ΔS = 0 and (b) ΔS = 1 channels. The parameters for the t-J model are t = 0.5 eV, J = 0.325 eV, V = 4 eV, Γ = 0.2 eV, and ωin = 0.14 eV. The parameters for the multi-orbital model are given in the main text. The incident photon energy is ωin = 2.5 eV, the inverse core hole lifetime is Γ = 0.2 eV, and the core hole potential is V = 4 eV. The results for the t-J model have been scaled by a factor of 0.26 such that the maximum intensity of the ΔS = 1 excitations are the same at the zone boundary.

After we have rescaled the spectra, we find excellent overall agreement between the two calculations: the amplitude of the broad arcs for the magnetic excitations, both in the ΔS = 0 and in the ΔS = 1 channels of the RIXS spectra are well captured by the t-J model. There are, however, minor quantitative differences related to the spectral weight of the excitations appearing near q = π/2a in the ΔS = 0 channel. In this channel, it is natural that the CuO4 spectrum is more noisy than the t − J because we are using the same number of DMRG states for both models, but the Hilbert space of the 20-site CuO4 cluster is much larger than the 20-site t − J cluster. Overall, the t-J model concentrates the magnetic excitations at slightly lower values of the energy loss in the ΔS = 0 channel. This discrepancy might be compensated for by taking a different value of J; however, this would come at the expense of the agreement in the ΔS = 1 channel. These differences should be kept in mind when one calculates the low-energy magnetic RIXS spectra using an effective t-J or single-band Hubbard model. Nevertheless, our results suggest that in the strong coupling limit, the magnetic RIXS spectrum can be described qualitatively by the effective t-J model.

Figure 4 shows that the overall agreement between the full multi-orbital model and the t-J model is much better in the ΔS = 1 channel than in the ΔS = 0 channel. We have verified that this discrepancy is not linked to loss of accuracy in our computational algorithm in the ΔS = 0 channel. Rather, we can naively understand this difference by recalling the role of charge fluctuations in the two magnetic excitation pathways. The ΔS = 1 RIXS excitations are possible in a system with strong spin-orbit coupling in the Cu 2p orbitals, which allows the spin of the core-hole to flip in the intermediate state of the RIXS process[4,10,19]. The ΔS = 0 pathway, however, requires a double spin-flip between neighboring Cu spins in the final state[4,19]. At the Cu L-edge, such processes occur due to charge fluctuations between the neighboring Cu sites in the intermediate state. The multi-orbital model treats such charge fluctuations differently owing to the presence of the ligand oxygen orbitals. This difference accounts for the discrepancy between the two models in the ΔS = 0 channel. At the Cu L-edge, however, the strong core-hole potential suppresses this difference by repelling holes from the site where it was created resulting in only minor differences between the predictions of the two models.

Concluding Remarks

We have presented a novel DMRG approach to computing the RIXS spectra and benchmarked this method against traditional ED. Using our DMRG algorithm, we can compute the RIXS spectra on 1D clusters much larger than those accessible to state-of-the-art ED methods. Using the proposed technique, we modeled the magnetic excitations probed by RIXS at the Cu L-edge in 1D antiferromagnets on the largest cluster sizes to date. We found that both the full multi-orbital cluster and the effective t-J model provide comparable descriptions of the magnetic excitations in the ΔS = 1 channel, while there were slight differences in the ΔS = 0 channel. This discrepancy was explained by noting the different ways that magnetic excitations are created in these channels.

Finally, we note that the bottleneck to RIXS simulations using ED is the exponential growth of the Hilbert space. Our approach shifts the computational burden to the availability of CPUs, thus opening the door to calculations for much larger systems. For example, one can envision extending this approach to the 2D models currently under active study by the DMRG community. Indeed, our RIXS-DMRG method is not restricted to 1D systems and can be applied to a 2D lattice geometry in the same sense specified in ref.[55]. One can always map a finite 2D N × N lattice with short-range interactions and hoppings into a 1D lattice with N2 sites and long-range interactions and hoppings. Once such a mapping is obtained, for instance, by drawing a 1D path that scans through the lattice following a “snake”-like pattern, one can straightforwardly apply the conventional DMRG algorithm, and similarly our RIXS-DMRG approach. It is well known, however, that there is a cost to this simplification. Since the interactions have long-range character, the numerical simulations become computationally more expensive as more DMRG states are needed to achieve converged results (area-law entanglement growth). Instead, DMRG applications to 2D systems usually adopt a mapping to multi-leg cylinders[55]. These lattice structures consist of coupled 1D chains or “legs” (usually up to 12 legs for spin systems[55], and up to 4–6 legs for fermionic systems[56-58]). DMRG simulations are then performed using periodic boundary conditions along the short (“y”) direction and open boundary conditions along the long (“x”) direction. The area-law entanglement growth is still the main limitation of the DMRG algorithm in this case as the number of legs of the cylinder increases[56-58]. The computation of dynamical response functions within the DMRG framework on multi-leg cylinders is computationally even more challenging, but it has been shown to be feasible. Indeed, the computation of spin (S(q, ω)) and charge (N(q, ω)) dynamical response functions of 4-leg t-t′-J cylinders has recently appeared in ref.[59]. Our RIXS-DMRG algorithm can be carried out on similar multi-leg cylinders.

Methods

The multi-orbital pd-Hamiltonian describing the corner-shared spin-chains, given in the hole-picture, isHere, denotes a sum over nearest neighbor orbitals; () creates a spin σ hole on the ith Cu 3 orbital (the jth O 2p orbital, γ = x, ±y); ε and ε are the on-site energies; () is the number operator for the Cu 3 orbital (the jth O 2p orbital); and are the Cu-O and O-O overlap integrals, respectively (the ij and jj′ dependence only indicates the ± differences among hoppings, Fig. 1b); U and U are the onsite Hubbard repulsions of the Cu and O orbitals, respectively, and U is the nearest-neighbor Cu-O Hubbard repulsion. The phase convention for the overlap integrals is shown in Fig. 1b. In this work, we adopt (in units of eV) ε = 0, ε = 3, ε = 3.5, |t(| = 1.5 |t(| = 1.8, |t| = 0.75, U = 8, U = 4, and U = 1, following ref.[21].

In the limit of large U, one integrates out the oxygen degrees of freedom and maps Eq. (9) onto an effective spin-1/2 t-J Hamiltonian[52]Here, is the annihilation operator for a hole with spin σ at site i, under the constraint of no double occupancy, is the number operator, and S is the spin operator at site i.

To facilitate a direct comparison between the two models, one can extract the hopping t and exchange interaction J from an ED calculation of a two-plaquette Cu2O7 cluster with open boundary conditions[60]. Here, we obtain the hopping (t = 0.5 eV) by diagonalizing this cluster in the ()-hole sector, and setting 2t to be equal to the energy separation between the bonding and antibonding states of the Zhang-Rice singlet. Similarly, we can obtain the superexchange (J = 0.325 eV) by diagonalizing the cluster in the -hole sector, and setting the singlet-triplet splitting of the Cu (d9d9) configurations equal to J.

Data and Code availability

The data that support the findings of this study are available from the corresponding author upon request. In the Supplementary Note V we provide details about the code used to obtain the DMRG results.

Supplementary Material