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.


Supplementary Note I: Benchmarks on a 16-site t-J chain
In this note, we compare the results of our DMRG method against the spectrum obtained from Lanczos ED. Supplementary Figure 1 directly compares the results from the two methods applied to a L = 16 site t-J chain, where our DMRG approach gives perfect agreement with the ED results for both the XAS and RIXS spectra. Here, we have assumed parameter values typical for a Cu L-edge measurement performed on Sr 2 CuO 3 with t = 0.3, J = 0.25, an inverse core hole lifetime Γ = 0.3 eV, and a core-hole repulsion V C = 2.0 eV. The two methods give a resonant absorption peak in the XAS for an incident energy ω in = 0.1 eV.

Supplementary Note III: Comparison of the Two Models
In the main text, we compared results for the magnetic RIXS spectra of the t-J and multi-orbital pd model, each with 20 unit cells. Here, Supplementary Figure 3 presents a similar comparison but for a different value of the core-hole potential used for the t-J model. The parameters for the multi-orbital pd 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 VC = 4 eV. The results for the t-J model have been scaled by a factor of 0.30 such that the maximum intensity of the ∆S = 1 excitations is the same at the zone boundary.

Supplementary Note IV: Removing the Elastic Line From the DMRG Calculations
We first rewrite Eq. (3) of the main text to explicitly indicate the center site c The ∆S = 0 contribution computed with DMRG is then given by the expression whereP = 1 − |g g| projects out the ground-state contribution. The expectation values α i,σ |D i,σ |g (and their hermitian conjugates) are calculated in Step 3 of the algorithm, and used in Step 4. Here, the contribution to the elastic peak of the spectra is removed by the subtraction of α i,σ |D i,σ |g . The ∆S = 1 contribution of the RIXS spectrum is given by In this case, the elastic contribution is absent because [Ĥ, Supplementary Note V: DMRG++ For the numerical results shown in this manuscript, PsimagLite version 2.02 and DMRG++ version 5.03 were used. Below, we will describe the input files needed to reproduce the data in Fig. 2 of the main text, corresponding to the RIXS spectrum of a t − J chain. For convenience, we include the ground state input file, inputGS.inp, valid for a t − J chain with a system size L = 8, using only m = 100 DMRG states. This input could be executed quickly on a standard laptop as Once the ground state |g has been obtained, the Step 2 of the algorithm described in Fig. 1 of the main text proceeds computing the vector |α c,↑ for c = L/2−1 restarting from the ground state calculations. The input input_Step2.inp to accomplish this job is This input could be executed quickly on a standard laptop as Above, the core-hole interaction strength V C = 2.0 enters as a local chemical potential shift at site c = L/2 − 1 (see below the potentialV flag). We have also used a core-hole linewidth Γ = 0.2 (CorrectionVectorEta) and incident energy ω IN = 0.1 (CorrectionVectorOmega). The solver option CorrectionVectorTargetting with CorrectionVectorAlgorithm=Krylov flag computes the |α c,↑ vector in Eq. (5) of the main text using the Krylov approach. Finally, notice that the operator matrixD c,↑ to be applied to the ground state enters directly at the end the input file. The main algorithm now bifurcates at Step 3. In fact, one needs to compute the vector |α j,↑ (for some site j = 5, for instance) in order to get the ∆S = 1 contribution to the RIXS spectrum, while |α j,↓ is needed to get the ∆S = 0 part. We here report only the input input_Step3_j=5.inp needed to get |α j,↑ .  This input could be executed quickly on a standard laptop as Above, as opposed to the Step 2 input, the core-hole interaction strength V C = 2.0 enters as a local chemical potential shift at site j (see below the potentialV flag). As in Step 2, we have used the same core-hole linewidth Γ = 0.2 (CorrectionVectorEta) and incident energy ω IN = 0.1 (CorrectionVectorOmega). The solver option TargetingRixsStatic with CorrectionVectorAlgorithm=Krylov flag computes the |α j,↑ vector in Eq. (6) of the main text using the Krylov approach and keeps track (in the DMRG sense) of the previously calculated |α c,↑ in Step 2. Even in this case, the operator matrixD j,↑ to be applied to the ground state enters directly at the end the input file. As shown in Fig. (1) of the main text, in the Step 3 of the main algorithm, contributions coming from runs with inputs input_Step3_j=x.inp for all the sites x of the chain are needed.
Finally, we describe a generic input file needed in Step 4 to compute the real space RIXS correlator for ∆S = 1. Once the vectors |α c,↑ and |α j,↑ are loaded restarting from Step 3, the following input input_Step3_j=5_om=4.inp computes the correction-vector |x c,↓,↑ and, when properly executed as ./dmrg -f input_Step3_j=5_om=4.inp ':cddn.txt' the correlator I j,c (Ω = ω 4 ) = α j,↑ |D j,↓ |x c,↓,↑ . Notice that the energy loss Ω = ω 4 has been computed. Also, the option ':cddn.txt' means that the operatorD j,↓ is being applied properly to get I j,c (Ω = ω 4 ). The input input_Step3_j=5_om=4.inp reads The solver option TargetingRixsDynamic with CorrectionVectorAlgorithm=Krylov flag computes the |x c,↓,↑ vector in Eq. (7) of the main text using the Krylov approach and by reading and keeping track of the previously calculated |α c,↑ in Step 2 and |α j,↑ in Step 3. Analogously to Step 3, the operator matrixD † c,↓ to be applied to |α c,↑ enters directly at the end the input file. We have used the same a peak broadening η = 0.075 (CorrectionVectorEta) and computed the correction-vector |x c,↓,↑ for the energy loss Ω = ω 4 = 0.4 (CorrectionVectorOmega). Finally, as shown in Fig. (1) of the main text, in the Step 4 of the main algorithm, contributions coming from runs with inputs input_Step4_j=x_om=y.inp for all the sites x of the chain and for all the energy losses y are needed.
In the last part of this section, we introduce the input file inputCuO.inp for the realistic multiorbital pd model introduced in the main text. Below, we consider a Cu 4 O 13 cluster, with open boundary conditions and use the same parameter values introduced in the main text. potentialV 34 3 3.5 0 3.5 3 3.5 0 3.5 3 3.5 0 3.5 3 3.5 0 3.5 3 3 3.5 0 3.5 3 3.5 0 3.5 3 3.5 0 3.5 3 3.5 0 3.5 3 InfiniteLoopKeptStates=128 FiniteLoops 5 7 1000 0 -15 1000 0 15 1000 0 -15 1000 0 15 1000 0 TargetElectronsUp=2  TargetElectronsDown=2  Threads=1 OutputFile=dataGS.txt SolverOptions=twositedmrg,MatrixVectorKron,useSvd Version=version TruncationTolerance=1e-9 For the hopping matrix, a snake-like one dimensional geometry has been used for the cluster. Also, it is worth mentioning that the matrix contains the repulsion term between the copper and oxygen orbitals, with interaction strength U pd = 1 eV . Finally, given the detailed discussion given above for the t-J model, inputs for the computation of the RIXS spectra at the Cu-L edge for the multiorbital pd model can be obtained straightforwardly.