Taming Convergence in the Determinant Approach for X-Ray Excitation Spectra

TL;DR

This paper introduces theorems based on SVD analysis to estimate the significance of many-electron configurations in x-ray spectra simulations, improving convergence understanding in the determinant formalism.

Contribution

The work provides a theoretical framework to quickly estimate the contribution of higher-order configurations, aiding convergence analysis in determinant-based x-ray spectra calculations.

Findings

Theorems based on SVD quantify configuration contributions.

Application to metallic systems up to fifth excitation order.

Improved understanding of convergence behavior in determinant formalism.

Abstract

A determinant formalism in combination with \emph{ab initio} calculations proposed recently has paved a new way for simulating and interpreting x-ray excitation spectra. The new method systematically takes into account many-electron effects in the Mahan-Nozi\'eres-De Dominicis (MND) theory, including core-level excitonic effects, the Fermi-edge singularity, shakeup excitations, and wavefunction overlap effects such as the orthogonality catastrophe, all within a universal framework using many-electron configurations. A heuristic search algorithm was introduced to search for the configurations that are important for defining spectral lineshapes, instead of enumerating them in a brute-force way. The algorithm has proven to be efficient for calculating \ce{O} $K$ edges of transition metal oxides, which converge at the second excitation order (denoted as $f^{(n)}$ with $n = 2$), i.e., the final-state configurations with two \emph{e-h} pairs (with one hole being the core hole). However, it remains unknown how the determinant calculations converge for general cases and at which excitation order $n$ one should stop the determinant calculation. Even with the heuristic algorithm, the number of many-electron configurations still grows exponentially with the excitation order $n$. In this work, we prove two theorems that can indicate the order of magnitude of the contribution of the $f^{(n)}$ configurations, so that one can estimate their contribution very quickly without actually calculating their amplitudes. The two theorems are based on the singular-value decomposition (SVD) analysis, a method that is widely used to quantify entanglement between two quantum many-body systems. We examine the $K$ edges of several metallic systems with the determinant formalism up to $f^{(5)}$ to illustrate the usefulness of the theorems.

Tables

Table 1. Table 1: The largest five singular values of the zeta matrices of the studies systems, and their accumulative products.

System	$s^1$	$s^2$	$s^3$	$s^4$	$s^5$	$s^1{s}^2$	$s^1{s}^2{s}^3$	$s^1{s}^2{s}^3{s}^4$	$s^1{s}^2{s}^3{s}^4{s}^5$
1D Chain	9.9683	1.0287	0.2810	0.0592	0.0094	10.2545	2.8820	0.1707	0.0016
Li	478.8850	1.9354	0.9774	0.6528	0.4657	926.8103	905.8236	591.2829	275.3625
\ceTiO2	7092.6499	0.5008	0.4417	0.3118	0.2006	3552.0365	1569.0875	489.1961	98.1572
\ceRuO3	267.8267	1.5267	0.5323	0.3790	0.2721	408.8852	217.6645	82.4916	22.4428
Graphite	7.8617	0.5691	0.3137	0.2325	0.1913	4.4744	1.4035	0.3264	0.0624
\ceCu	1076.8833	1.7903	0.4386	0.4087	0.3919.	1927.8910	845.5359	345.5765	135.4185
\ceCrO2$\uparrow$	37.0903	1.5178	0.6071	0.4586	0.3631	56.2962	34.1768	15.6735	5.6917
\ceCrO2$\downarrow$	310.3019	0.3667	0.3556	0.3157	0.1634	113.7812	40.4625	12.7736	2.0867
\ceLiCoO2	113.3031	0.5227	0.4932	0.3257	0.2632	59.2231	29.2065	9.5120	2.5036

Equations

$$\begin{array}{r}⟨{\Psi }_f|\mathcal{O}|{\Psi }_i⟩=\sum _c(A_c^f{)}^{\ast }⟨{\psi }_c|o|{\psi }_h⟩\end{array}$$

$$\begin{array}{r}{A}_c^f=\det \left[\begin{array}{ccccc}{\xi }_{f_1,1}& {\xi }_{f_1,2}& \cdots & {\xi }_{f_1,N}& {\xi }_{f_1,c}\\ {\xi }_{f_2,1}& {\xi }_{f_2,2}& \cdots & {\xi }_{f_2,N}& {\xi }_{f_2,c}\\ ⋮& & \ddots & & ⋮\\ {\xi }_{f_{N+1},1}& {\xi }_{f_{N+1},2}& \cdots & {\xi }_{f_{N+1},N}& {\xi }_{f_{N+1},c}\end{array}\right]\end{array}$$

$$\begin{array}{rl}{A}^f& \equiv ⟨{\Psi }_i|\mathcal{O}|{\Psi }_f⟩=\det {\mathbf{A}}^f\\ {\mathbf{A}}^f& =\left[\begin{array}{ccccc}{\xi }_{f_1,1}& {\xi }_{f_1,2}& \cdots & {\xi }_{f_1,N}& {\sum }_c{\xi }_{f_1,c}{w}_c^{\ast }\\ {\xi }_{f_2,1}& {\xi }_{f_2,2}& \cdots & {\xi }_{f_2,N}& {\sum }_c{\xi }_{f_2,c}{w}_c^{\ast }\\ ⋮& & \ddots & & ⋮\\ {\xi }_{f_{N+1},1}& {\xi }_{f_{N+1},2}& \cdots & {\xi }_{f_{N+1},N}& {\sum }_c{\xi }_{f_{N+1},c}{w}_c^{\ast }\end{array}\right]\end{array}$$

$$\begin{array}{r}{a}_i=\left[\begin{array}{cccc}{\xi }_{i,1}& \cdots & {\xi }_{i,N}& {\sum }_c{\xi }_{i,c}{w}_c^{\ast }\end{array}\right]\end{array}$$

$$\begin{array}{r}\left[\begin{array}{c}{a}_{N+1}\\ a_{N+2}\\ ⋮\\ a_M\end{array}\right]=\left[\begin{array}{cccc}0& 0& \cdots & 1\\ {\zeta }_{N+2,1}& {\zeta }_{N+2,2}& \cdots & {\zeta }_{N+2,N+1}\\ ⋮& ⋮& & ⋮\\ {\zeta }_{M,1}& {\zeta }_{M,2}& \cdots & {\zeta }_{M,N+1}\end{array}\right]\left[\begin{array}{c}{a}_1\\ a_2\\ ⋮\\ a_{N+1}\end{array}\right]\end{array}$$

$$\begin{array}{rl}{A}^c& ={\zeta }_{c,N}{A}^{\text{ref}}\\ A^{c;c_1,v_1}& =\det \left[\begin{array}{cc}{\zeta }_{c,v_1}& {\zeta }_{c,N}\\ {\zeta }_{c_1,v_1}& {\zeta }_{c_1,N}\end{array}\right]A^{\text{ref}}\\ A^{c;c_1,v_1;c_2,v_2}& =\det \left[\begin{array}{ccc}{\zeta }_{c,v_2}& {\zeta }_{c,v_1}& {\zeta }_{c,N}\\ {\zeta }_{c_1,v_2}& {\zeta }_{c_1,v_1}& {\zeta }_{c_1,N}\\ {\zeta }_{c_2,v_2}& {\zeta }_{c_2,v_1}& {\zeta }_{c_2,N}\end{array}\right]A^{\text{ref}}\end{array}$$

$$\begin{array}{r}{\zeta }_{ij}\sim a_i{b}_j\end{array}$$

$$\begin{array}{r}{\zeta }_{ij}=\sum _k{s}^k{a}_i^k{b}_j^k\end{array}$$

$$\begin{array}{r}D=\sum _{k_1<k_2<\cdots <k_n}{s}^{k_1}{s}^{k_2}\cdots s^{k_n}{D}_{[k_{\mu }]}^a{D}_{[k_{\mu }]}^b\end{array}$$

$$\begin{array}{r}D=\sum _{l_1{l}_2\cdots l_n=1}^n{\epsilon }_{l_1{l}_2\cdots l_n}{\zeta }_{1l_1}{\zeta }_{2l_2}\cdots {\zeta }_{n{l}_n}\end{array}$$

$$\displaystyle\begin{split}D&=\sum^{n}_{l_{1}l_{2}\cdots l_{n}=1}\varepsilon_{l_{1}l_{2}\cdots l_{n}}\big{(}\sum_{k_{1}}s^{k_{1}}a^{k_{1}}_{1}b^{k_{1}}_{l_{1}}\big{)}\cdots\big{(}\sum_{k_{n}}s^{k_{n}}a^{k_{n}}_{n}b^{k_{n}}_{l_{n}}\big{)}\\ &=\sum^{n}_{l_{1}l_{2}\cdots l_{n}=1}\varepsilon_{l_{1}l_{2}\cdots l_{n}}\\ &\times\bigg{[}\sum_{k_{1}k_{2}\cdots k_{n}}s^{k_{1}}s^{k_{2}}\cdots s^{k_{n}}a^{k_{1}}_{1}b^{k_{1}}_{l_{1}}\cdots a^{k_{n}}_{n}b^{k_{n}}_{l_{n}}\bigg{]}\\ &=\sum_{k_{1}k_{2}\cdots k_{n}}s^{k_{1}}s^{k_{2}}\cdots s^{k_{n}}a^{k_{1}}_{1}\cdots a^{k_{n}}_{n}\\ &\times\bigg{[}\sum^{n}_{l_{1}l_{2}\cdots l_{n}=1}\varepsilon_{l_{1}l_{2}\cdots l_{n}}b^{k_{1}}_{l_{1}}\cdots b^{k_{n}}_{l_{n}}\bigg{]}\end{split}$$

$$\begin{array}{r}{D}_{[k_{\mu }]}^b\equiv \sum _{l_1{l}_2\cdots l_n=1}^n{\epsilon }_{l_1{l}_2\cdots l_n}{b}_{l_1}^{k_1}\cdots b_{l_n}^{k_n}\end{array}$$

$$\begin{array}{rl}D& =\sum _{k_1{k}_2\cdots k_n}{s}^{k_1}{s}^{k_2}\cdots s^{k_n}{a}_1^{k_1}\cdots a_n^{k_n}{D}_{[k_{\mu }]}^b\\ & =\sum _{k_1<k_2<\cdots <k_n}{s}^{k_1}{s}^{k_2}\cdots s^{k_n}\\ & \times \sum _{l_1{l}_2\cdots l_n}{a}_1^{l_1}\cdots a_n^{l_n}{\epsilon }_{l_1{l}_2\cdots l_n}{D}_{[k_{\mu }]}^b\end{array}$$

$$\begin{array}{r}D=\sum _{k_1<k_2<\cdots <k_n}{s}^{k_1}{s}^{k_2}\cdots s^{k_n}{D}_{[k_{\mu }]}^a{D}_{[k_{\mu }]}^b\end{array}$$

$$\begin{array}{r}|D^{n\times n}|\le \sum _{k_1<k_2<\cdots <k_n}{s}^{k_1}{s}^{k_2}\cdots s^{k_n}\end{array}$$

$$\begin{array}{r}{D}_{[k_{\mu }]}^a=\left|\begin{array}{cccc}{a}_1^1& a_2^1& \cdots & a_n^1\\ a_1^1& a_2^2& \cdots & a_n^2\\ ⋮& ⋮& \ddots & ⋮\\ a_1^n& a_2^n& \cdots & a_n^n\end{array}\right|\end{array}$$

$$\begin{array}{r}|{\mathbf{a}}_{1\times n}^k|\le 1\end{array}$$

$$\displaystyle\begin{split}|D^{n\times n}|&\leq\sum_{k_{1}<k_{2}<\cdots<k_{n}}\big{|}s^{k_{1}}s^{k_{2}}\cdots s^{k_{n}}D^{a}_{[k_{\mu}]}D^{b}_{[k_{\mu}]}\big{|}\\ &\leq\sum_{k_{1}<k_{2}<\cdots<k_{n}}s^{k_{1}}s^{k_{2}}\cdots s^{k_{n}}\end{split}$$

$$\displaystyle\begin{split}H=-t\sum^{N}_{j=1}\big{(}c^{\dagger}_{j+1}c_{j}+c^{\dagger}_{j}c_{j+1}\big{)}\end{split}$$

$$\begin{array}{r}|k\sigma ⟩=\frac{1}{\sqrt{N}}\sum _{j=0}^N{e}^{ikj}|j\sigma ⟩\end{array}$$

$$\begin{array}{r}⟨k\sigma |j=0,\sigma ⟩=\frac{1}{\sqrt{N}}\end{array}$$

$$\begin{array}{rl}{f}^{(1)}:& s^1\\ f^{(2)}:& s^1{s}^2\\ f^{(3)}:& s^1{s}^2{s}^3\\ \cdots & \end{array}$$

Full text

I Introduction

X-ray spectroscopy has become increasingly important for providing insights into many problems in materials characterization at microscopic scale de2008core ; fadley2010x ; singh2010advanced ; ament2011resonant ; van2016x , especially in recent times when this area is propelled by the development of light sources enabled by free-electron lasers emma2010first ; schlotter2012soft ; lemke2013femtosecond . We consider resonant x-ray excitations, where a core electron of a specific type of atoms is promoted into orbitals localized on or near that atomic site, revealing a wealth of information of local chemical environment and electronic structure. The near-edge part of a x-ray spectrum, i.e., a few eVs above onset, is of particular interest, where most intriguing physics and chemical processes are related to de2008core . The interpretation of x-ray spectra, however, is often a nontrivial task that requires accurate first-principles modeling of both atomic and electronic structure of interest and their associated spectra prendergast2006x ; cabaret2010first ; gougoussis2009first . While the structural properties of a wide range of materials can be mostly captured by density-functional theory (DFT) parr1980density ; dudarev1998electron ; bickelhaupt2000kohn ; wang2006oxidation ; hautier2012accuracy , predicting x-ray excited-state spectra in a reliable and efficient manner presents a greater theoretical challenge.

There has been a broad spectrum of theoretical approaches for simulating x-ray excitation spectra. At one extreme, exact diagonalization chen1991electronic ; hybertsen1992model ; haverkort2012multiplet ; kourtis2012exact ; uldry2012systematic ; chen2013doping ; jia2016using ; lu2017nonperturbative has been applied to rigorously solve the many-electron Hamiltonian. This method represents the most accurate solutions and is most suitable for localized x-ray excitations that occur within a few atomic orbitals, due to the exponential growth of the many-electron Hilbert space. Wannier down-folding is typically required for reducing the size of the Hamiltonian haverkort2012multiplet ; uldry2012systematic and the Coulomb interaction is usually simplified as an on-site Hubbard $U$ term. Similar methodology is also employed in the quantum chemistry community within the configuration interaction sherrill1999configuration ; tubman2016deterministic and other post-Hartree-Fock methods bartlett2007coupled ; szalay2011multiconfiguration , in which realistic Coulomb interactions and superposition of many-electron configurations are considered explicitly in the calculation. Likewise, increasing the size of the excited-state calculation is hindered by the exponential barrier associated with the size of the Hilbert space, although placing restrictions on the active orbital space can mitigate the problem to some extent roemelt2013excited ; pinjari2014restricted ; maganas2017restricted . Due to its rigor, this class of methods can be used for system with strong electron correlation, but due to its computational inefficiency, it is limited to specific cases such as small molecules roemelt2013excited ; maganas2017restricted or clusters lu2017nonperturbative and $3d$ metal $L$ edges (dominate d by localized atomic multiplet effects) haverkort2012multiplet ; jia2016using so far.

At the other extreme, excited-state spectra of condensed phases can be routinely obtained using ab initio methods based on DFT and many-body perturbation theory (MBPT), where electron correlation is treated with less rigor but higher efficiency. There are two representative methods. One is the Delta-self-consistent-field ($\Delta$SCF) approach nyberg1999core ; taillefumier2002x ; prendergast2006x ; liang2017accurate that treats the core hole as an external potential, and maps a many-electron excited state to an empty orbital. Many-electron response is taken into account by DFT electronic relaxation, and then the transition matrix elements are calculated using the Fermi’s Golden rule. In previous formalisms, however, only one-body orbitals are employed for final states, which does not account for proper time ordering of the many-electron processes in x-ray excitations liang2017accurate ; liang2018quantum , leading to possible failures in predicting the intensity of an absorption edge. The other method is the core-level Bethe-Salpeter Equation (BSE) shirley1998ab ; olovsson2009all ; vinson2011bethe within MBPT, that utilizes DFT ground state as a zero-order approximation and incorporate many-electron correlations in a perturbative manner. A subset of Feynman diagrams generated by the direct and exchange kernels are included to account for e-h interactions (excitonic effects). and hence the correct time-ordering is retained in this approach. Formulations akin to the BSE are also adopted in some versions of time-dependent DFT (TDDFT) besley2009time ; besley2010time ; lopata2012linear ; zhang2012core , in which exchange-correlation kernels play the roles of e-h interaction kernel, and real-time evolution could be involved lopata2012linear . There are many successful application of the $\Delta$SCF core-hole approach nyberg1999core ; taillefumier2002x ; prendergast2006x ; pascal2014finite ; velasco2014structure and the core-level BSE vinson2011bethe ; liang2017accurate ; zhan2017enabling in extended systems, and even in molecules nyberg1999core where the Coulomb interaction is not well screened.

I.1 Motivations and Advantages of the Determinant Formalism

Following the philosophy of the post DFT methods, we recently proposed a determinant formalism liang2017accurate ; liang2018quantum for simulating x-ray excitation spectra based on the one-body core-hole approach. The determinant formalism has three main advantages: First, it provides an exact solution to all the many-electron effects considered within the Mahan-Noziéres-De Dominicis (MND) model mahan1967excitons ; nozieres1969singularities ; ohtaka1990theory ; mahan2013many , which is beyond the scope of e-h attraction in the BSE. These effects include the power-law edge singularity as considered by Mahan using Feynman diagrams mahan2013many , and the many-body wavefunction overlap effect considered by Anderson in the orthogonality catastrophe anderson1967pw . In a less dramatic manner, the latter often manifests as shakeup effects kim1975x ; stern1983many ; calandra2012k , but not edge-rounding effects. The shake-up effects are the fragmentation of an absorption feature to higher energy, which would not be described by the BSE.

Secondly, the determinant formalism adopts many-electron configurations in the calculation, and hence provides a conceptually simple picture for understanding x-ray excited states. In the x-ray final-state system, a many-body state is simply a single determinant (single reference) of the occupied Kohn-Sham orbitals. Each excited state is now mapped to its composite orbitals. We denote the group of singles as $f^{(1)}$, which comprises a core hole and an electron; the groups of doubles as $f^{(2)}$, which comprises an additional valence e-h pair; and so forth. The $n$ as in $f^{(n)}$ can be understood as the order of excitations. This enables a straightforward assignment of absorption features to excited states, making the interpretation of the spectrum simple.

Thirdly, an alternative solution for solving the MND Model using DFT orbitals as input was introduced in Ref. wessely2005ab ; wessely2006dynamical ; wessely2007dynamical ; ovcharenko2013specific , and good agreement has been achieved in some carbon systems and metals. In this method, a time-dependent matrix integral equation needs to be solved and matrix inversion is required for each point on the discretized time axis. Then a Fourier transformation is performed to obtain the spectrum. The determinant approach we proposed only requires a one-shot matrix inversion and no real-time evolution is involved. It is physically more intuitive and computationally less complex.

Finally, although many-electron configurations are employed for condensed-phase systems with hundreds to thousands of electrons in a supercell, the determinant calculation is not at all intractable. We have developed a heuristic search algorithm for finding the many-electron configurations that are important for determining the x-ray spectra liang2018quantum . For transition metal oxides (TMOs), it is found the x-ray absorption spectra (XAS) can be well defined by just $10^{5}$ configurations up to the $f^{(2)}$ order (not their superposition and no diagonalization of many-body Hamiltonian is needed), which are inexpensive calculations given the current computational capability.

It is, however, unknown yet how expensive the determinant calculations are in other systems apart from TMOs. It remains unclear if $f^{(n)}$ configurations with $n>2$ are important for shaping the x-ray excitation spectra. Even with the heuristic search algorithm liang2018quantum that truncates the number of many-electron configurations, the number of meaningful $f^{(n)}$ configurations still tends to diverge exponentially with respect to the shakeup order $n$. If one can estimate the spectral contributions from the $f^{(n)}$ configurations before actually calculating them, it will be of great benefit for saving computational resources. In this work, we propose a useful criterion for estimating the contributions of $f^{(n)}$, which helps one to decide at what $n$ one should stop the determinant calculation. It is based on an singular-value decomposition (SVD) analysis of the $\zeta$-matrix, the auxiliary matrix used to obtain the determinant transition amplitude. An SVD analysis of the $\zeta$-matrix reflects how the one-body basis set is rotated due to the perturbation of the core hole, and how the final-state occupied manifold is entangled with the initial state (typically the many-electron ground state). Besides TMOs, we have chosen several metallic systems such as \ceLi and \ceCu metal and performed exhaustive calculations up to $n=5$ to test the convergence criterion in this work.

I.2 Review of the Determinant Formalism

The central formula for calculating the x-ray absorption amplitudes from the ground-state (initial state) $|\Psi_{i}\rangle$ to a specific final state $|\Psi_{f}\rangle$ is liang2017accurate ; liang2018quantum

[TABLE]

in which the transition amplitude also takes a determinant form

[TABLE]

and $\xi_{ij}=\langle\psi_{j}|\tilde{\psi}_{i}\rangle$, where $|\tilde{\psi}_{i}\rangle$ is one final-state orbital and $|\psi_{j}\rangle$ is one initial-state orbital. Alternatively, $A^{f}_{c}$ can be regarded as a “Slater determinant”of a set of final-state orbitals $(\tilde{\psi}_{f_{1}},\tilde{\psi}_{f_{2}},\cdots,\tilde{\psi}_{f_{N+1}})$, expanded over the lowest $N+1$ initial-state orbitals $(\psi_{1},\psi_{2},\cdots,\psi_{N+1})$ (rather than $(\bm{r}_{1},\bm{r}_{2},\cdots,\bm{r}_{N+1})$). The first step we use to simplify the calculation is to move the summation over $c$ into the determinant expression:

[TABLE]

There are an enormous number of combinations of $(f_{1},f_{2},\cdots,f_{N+1})$, representing possible final-state configurations, but it is only meaningful to visit a small subset of this space. The determinant amplitude $A^{f}$ is a significant number only when the most of $(f_{1},f_{2},\cdots,f_{N+1})$ overlaps significantly with the lowest-energy configuration $(1,2,\cdots,N+1)$, because in most realistic materials, the core hole only slightly perturbs the system and the orthogonality catastrophe does not occur. Therefore, $(f_{1},f_{2},\cdots,f_{N+1})$ may only differ from $(1,2,\cdots,N+1)$ by a few indices.

To simplify the notation, we may denote $(f_{1}=1,f_{2}=2,\cdots,f_{N}=N,f_{N+1}=c)$ as a single, or $f^{(1)}$ configuration, $(f_{1}=1,f_{2}=2,\cdots,f_{v_{1}-1}=v_{1}-1,f_{v_{1}}=v_{1}+1,\cdots,f_{N-1}=N,f_{N}=c,f_{N+1}=c_{1})$ as a double, or $f^{(2)}$ configuration, and so forth.

Each line of $\bm{A}^{f}$ can be denoted as:

[TABLE]

Next we can introduce the $\zeta$-matrix to calculate all the $A^{f}$’s in an iterative manner. The $\zeta$-matrix is the linear transformation from $a_{i}$’s of the occupied orbitals ($i\leq N$ plus $a_{N+1}$) to the $a_{i}$’s of the empty orbitals (i ¿ N)

[TABLE]

Rewriting the above matrix multiplication in a compact form, we have $\bm{A}^{\text{new}}=\bm{\zeta}\bm{A}^{\text{ref}}$, where $(\bm{\zeta})_{ij}=\zeta_{N+i,j}$. Thus the $\zeta$-matrix can be obtained from $\bm{\zeta}=\bm{A}^{\text{new}}(\bm{A}^{\text{ref}})^{-1}$.

With the auxiliary $\zeta$ matrix, we can quickly calculate the determinants for many excited-state configurations without repeatedly using the $\mathcal{O}(N^{3})$ determinant algorithm. Instead, an $\mathcal{O}(1)$ algorithm can be used. Once the determinant for the ground state, $A^{\text{ref}}$, is obtained, the determinant of each excited state can be computed by multiplying $A^{\text{ref}}$ by a pre-factor. The pre-factor is a small determinant that composes the matrix elements of the $\zeta$ matrix. For example, the amplitude of a single ($f^{(1)}$), double ($f^{(2)}$), and triple ($f^{(3)}$) configuration can be obtained respectively as

[TABLE]

These are the essential formulae to obtain the amplitudes for the excited-state configuration with the $\mathcal{O}(1)$ updating algorithm. A $f^{\text{(n)}}$ configuration corresponds to a $n\times n$ sub-determinant of the $\zeta$ matrix.

Directly enumerating all such sub-determinants will be a computationally expensive task. Also, it is not necessary to do so because the $\zeta$ matrix could be a sparse matrix. In a previous work liang2018quantum , we have proposed a heuristic search algorithm for quickly finding significant sub-determinants of a sparse matrix. Our general search algorithm is not merely specific to the $\zeta$ matrix for x-ray spectroscopic problems.

What was overlooked in the previous work, however, is that the $\zeta$ matrix for x-ray excitations does have structures. It can be seen that the $\zeta{}$ matrices for \ceSiO2, \ceTiO2, \ceCrO2 all display a vertical (horizontal) stripe pattern (Fig. 8 (d) of Ref. liang2018quantum ). This stripe pattern can also be seen in other chosen examples, which are discussed later in this work. These stripe patterns imply a further simplification of the $\zeta$ matrix that it can be approximately expressed as the Kronecker product of two vectors:

[TABLE]

If $\zeta_{ij}\sim a_{i}b_{j}$ strictly holds, then all the $n\times n$ sub-determinants for $n>1$ will be exactly zero and only the $f^{(1)}$ amplitudes are non-vanishing. It is the deviation of $\zeta_{ij}$ from $a_{i}b_{j}$ that determines the size of higher-order terms $f^{(n)}$ ($n>1$). If one can expand $\zeta_{ij}$ into just a few terms, then it is highly probable that the size of higher-order terms can be quickly estimated. In this regard, a singular-value decomposition (SVD) of the $\zeta$ matrix is most relevant for this problem. SVD has been widely used to analyze the entanglement structure between two quantum many-body systems vidal2007classical ; li2008entanglement ; schollwock2011density .

In the following discussion, we first provide and prove two theorems that will give upper bounds on the size of the $f^{(n)}$ terms, for a specific $n$, using SVD analysis for the $\zeta$ matrix. The bounds will enable one to determine the contribution of the $f^{(n)}$ terms to the x-ray excitation spectrum, without explicitly calculating these terms. This will save a substantial amount of computational cost and help one obtain a meaningful x-ray excitation spectrum faster. Then we apply the theorems to several small band-gap and metallic systems, in which higher-order terms $f^{(n)}$ ($n>1$) are expected to contribute to the spectrum significantly. It is, however, found that in none of the chosen systems, the contribution from $f^{(n)}$ ($n>2$) can significantly alter the spectral lineshapes (more precisely, the peak intensity ratios). In other words, the spectra have already taken shape at the order of $n=2$.

II Results and Discussion

II.1 Two theorems about sub-determinants

Theorem 1**.**

Let $D$ be the determinant of an $n\times n$ sub-matrix that spans over rows $i_{1},i_{2},\cdots,i_{n}$ and columns $j_{1},j_{2},\cdots,j_{n}$ of an$N\times M$ matrix $\bm{\zeta}$. Suppose the singular-value decomposition (SVD) of $\bm{\zeta}$ is

[TABLE]

where $\{s^{k}\}$ are the singular values of $\bm{\zeta}$, and $a^{k}_{i}$ ($b^{k}_{j}$) is a normalized vector for a given $k$. Then the determinant $D$ can be expanded as the summation

[TABLE]

in which $D^{a}_{[k_{\mu}]}$ ($D^{b}_{[k_{\mu}]}$) is the determinant of the sub-matrix that spans over rows $i_{1},i_{2},\cdots,i_{n}$ and columns $k_{1},k_{2},\cdots,k_{n}$ of the matrix $a^{k}_{i}$ ($b^{k}_{j}$).

Proof.

Without loss of generality, we may assume $i_{1}=j_{1}=1$, $i_{2}=j_{2}=2$, $\cdots$, $i_{n}=j_{n}=n$. According to the definition of determinant, $D$ can be expanded using the Levi-Civita symbol:

[TABLE]

Inserting the SVD expression of the matrix element $\zeta_{ij}$ as in Eq. [3] (examples of $\bm{\zeta}$ can be found in Fig. 3),

[TABLE]

Note that the inner summation with respect to $l_{\nu}$ for a specific $\{k_{\mu}\}$ ($\mu,\nu=1,2,\cdots,n$) gives rise to a determinant, which can be denoted as $D^{b}_{[k_{\mu}]}$:

[TABLE]

This determinant corresponds to the sub-matrix $b^{k}_{l}$ formed by row $1,2,\cdots,n$ and column $k_{1},k_{2},\cdots,k_{n}$.

In the outer summation of Eq. [11], each index $k_{\mu}$ can range from $1$ to the number of singular values. $D^{b}_{[k_{\mu}]}$ is non-zero only when $k_{1},k_{2},\cdots,k_{n}$ are not equal to each other, thus placing constraints on values of $k_{\mu}$ in the outer summation.

Next, we may consider the case where the $n$ values of $k_{\mu}$ are taken from $1,2,\cdots,n$, without loss of generality. Note that there is no ordering presumed in the outer summation of $k_{\mu}$, therefore summing over all $k_{\mu}$’s will generate $n!$ permutations of $1,2,\cdots,n$. For these $n!$ permutations, the corresponding $D^{b}_{[k_{\mu}]}$ will have the same absolute value, and its $\pm$ sign depends on whether the permutation is odd or even, due to the nature of determinant. If we enforce ordering $k_{1}<k_{2}<\cdots<k_{n}$, we may use the Levi-Civita symbol to represent the sign due to permutation.

[TABLE]

where the tuple $l_{1}l_{2}\cdots l_{n}$ iterate over all permutation of $k_{1}<k_{2}<\cdots<k_{n}$. Therefore,

[TABLE]

∎

Theorem 2**.**

The absolute value of the determinant of an $n\times n$ sub-matrix of $\bm{\zeta}$ is bound by the product of its largest $n$ singular values provided by SVD

[TABLE]

where the singular values satisfy $s^{1}\geq s^{2}\geq\cdots\geq s^{n}>0$.

Proof.

We may start by proving the absolute value of the determinant $D^{a}_{[k_{\mu}]}$ ($D^{b}_{[k_{\mu}]}$) is bound by $1$. Without loss of generality, we may again assume $i_{1}=1,i_{2}=2,\cdots,i_{n}=n$ and $k_{1}=1,k_{2}=2,\cdots,k_{n}=n$. According to its definition:

[TABLE]

Each row of the above determinant is a vector $\bm{a}^{k}_{1\times n}=(a^{k}_{1},a^{k}_{2},\cdots,a^{k}_{n})$, which is a truncation of a full $1\times N$ vector $\bm{a}^{k}=(a^{i}_{1},a^{i}_{2},\cdots,a^{i}_{N})$. $N$ is the first dimension of $\bm{\zeta}$ and $n\leq N$. Since the matrix $a^{k}_{i}$ is obtained from SVD of $\bm{\zeta}$, each of its row vector is normalized to 1: $|\bm{a}^{k}|=1$, and therefore:

[TABLE]

The geometric meaning of the determinant $D^{a}_{[k_{\mu}]}$ is the volume of a parallelepiped spanned by $n$ vectors $\bm{a}^{k}$, where $k=1,2,\cdots,n$. Since the length of each of its edge $|\bm{a}^{k}_{1\times n}|\leq 1$, then the volume of the parallelepiped will be no larger than 1, and thus $|D^{a}_{[k_{\mu}]}|\leq 1$ ($|D^{b}_{[k_{\mu}]}|\leq 1$).

Using the conclusion of theorem 1,

[TABLE]

Note that all the singular values of $\zeta$ are non-negative. ∎

II.2 An analysis of $\zeta$-matrices with singular-value decomposition

We analyze several representative examples here to illustrate the usefulness of the two theorems for $\bm{\zeta}$ matrices: a 1D single-atom chain at half filling, the \ceC $K$ edge of graphite, the \ceCu $K$ edge copper, the \ceLi $K$ edge of lithium metal, and the \ceO $K$ edge of rutile \ceTiO2, \ceCrO2, \ceRuO3, and \ceLiCoO2. The crystal structures of \ceRuO3, and \ceLiCoO2 are shown in Fig. 1. Both systems are layered structures that will allow lithium insertion/removal, and are being studied as cathode prototypes of rechargeable batteries. The chosen systems are gapless except for \ceTiO2 and \ceLiCoO2, which have a DFT-PBE band gap of $2.1$ and $1.8$ eV (with a Hubbard $U$ value of $3.3$ eV on the \ceCo atom) respectively.

The tight-binding model for the 1D chain reads:

[TABLE]

where a periodic boundary condition is employed and each site has double occupancy. The above model is considered as the initial-state Hamiltonian $H_{i}$. The core-hole potential is assumed to act on a single site (the site at $j=0$): $V_{c}=\Delta Vc^{\dagger}_{0}c_{0}$, thus the final-state Hamiltonian is $H_{f}=H_{i}+V_{c}$. It can be solved that the 1p wavefunctions for $H_{i}$ are:

[TABLE]

where $k=0,\frac{1}{N},\cdots,\frac{N-1}{N}$ and $\sigma=\uparrow,\downarrow$. We then define the 1p XAS matrix element as:

[TABLE]

In the actual calculation, we choose the number of sites $N=200$, the number of electrons $N_{e}=200$ (half-filled), $t=1$, a perturbation potential of $\Delta V=-100$ at the excited site to simulate the core hole effect (set to 100 for exaggeration). We find the determinant spectrum barely changes after $\Delta V<-100$.

First, we plot the $\zeta$ matrices for the XAS (\ceLi, \ceC, \ceO, and \ceCu $K$ edges) of six chosen systems in order to exemplify the aforementioned stripe pattern. Once we have completed the $\Delta$SCF calculation for both the initial and final state, the $\zeta$ matrices can be calculated using Eq. [5]. Several $\zeta$ matrices are shown in Fig. 2. Although the chosen systems are vastly different in terms of crystal symmetry and bonding nature, it can be seen that all the plotted $\zeta$ matrices display a cross-strip pattern that runs both vertically and horizontally. This cross-strip pattern suggested that the $\zeta$ matrix for x-ray excitations can be universally expressed in terms of Kronecker product of two vectors: $\zeta_{ij}\sim a_{i}b_{j}$, with a few small residual terms. Thus it becomes natural to introduce an SVD analysis to the $\zeta$ matrix, and apply the two theorems proved above to determine the spectral convergence with respect to shakeup order.

The twenty most significant singular values of the $\zeta$ matrices of the chosen systems are shown in Fig. 3 (the singular values are plotted logarithmically). The first 5 largest singular values (denoted as $s^{i}$) are also tabulated in Table. 1. It can immediately be seen for all the $\zeta$ matrices that only a few singular values are larger than $1$, especially for graphite, \ceTiO2, \ceLiCoO2 and the insulating (spin-down) channel of \ceCrO2 with only 1 $s^{i}>1.0$. The latter three are systems with a large band gap. The other systems only have two singular values larger than $1$. Note that in \ceLi, $s^{3}=0.977$ and $s^{4}=0.653$, which can be rounded to 1. In all the cases, $s^{i}<0.1$ for $i>10$, which are orders of magnitude smaller than $s^{1}$. This means only a few singular values $s^{i}$ are relevant for the analysis of shakeup orders.

With these numerical results and the properties of $\{s^{i}\}$, we can revisit the two theorems and discuss the contributions from different shake-up order $f^{(n)}$ to the x-ray spectra. Because the vast majority of $s^{i}$ are small, we could take the product of the leading $s^{i}$ to form an analysis of order of approximation. As indicated in the Theorem 1, the $n\times n$ sub-determinants corresponding to the $f^{(n)}$ terms are associated with the coefficients $s^{k_{1}}s^{k_{2}}\cdots s^{k_{n}}$. We may thus take the largest $n$ singular values $s^{1}$ and calculate their product to estimate the contribution of $f^{(n)}$:

[TABLE]

The accumulative product $s^{1}s^{2}\cdots s^{n}$ of some leading-order singular values are also shown in Table. 1. For \ceLi, \ceRuO3, \ceCu, \ceCrO2$\uparrow$, $s^{1}s^{2}$ is significantly larger than $s^{1}$ because $s^{2}>1.5$. This suggests the contribution of the $f^{(2)}$ terms is of the same order of magnitude of $f^{(1)}$. An extreme case is \ceLi in which even $s^{1}s^{2}s^{3}s^{4}\approx 591$ is larger to $s^{1}\approx 479$. This requires one to go beyond $f^{(1)}$ in the XAS determinant calculation for these systems.

II.3 XAS calculated by the determinant formalism

To test how good this empirical estimate is, we perform determinant calculations for XAS of the chosen systems, at least at the order of $f^{(2)}$. The obtained spectra decomposed by different shake-up orders are shown in Fig. 4 (odd rows). We find it is indeed true that the contribution from the $f^{(2)}$ terms is as significant as $f^{(1)}$ in \ceLi, \ceRuO3, \ceCu, \ceCrO2$\uparrow$ (as studied in Ref. liang2017accurate ; liang2018quantum ). In every case, the contribution from $f^{(2)}$ constitutes more than $40\%$ of the entire computed spectrum. In particular for \ceLi, a trend of convergence is only seen after including the $f^{(4)}$ terms. The contributions from $f^{(1)}$, $f^{(2)}$, $f^{(3)}$, $f^{(4)}$, and $f^{(5)}$ to the full spectrum ($f^{(1)}+f^{(2)}+f^{(3)}+f^{(4)}+f^{(5)}$) are $29.7\%$, $51.3\%$, $14.8\%$, $3.7\%$, and $0.5\%$. This is consistent with the fact that $s^{1}s^{2}\approx 927$, $s^{1}s^{2}s^{3}\approx 906$, and $s^{1}s^{2}s^{3}s^{4}\approx 591$ are comparable to $s^{1}\approx 479$ in this system. But note that the contribution $f^{(n)}$ is not entirely proportional to the accumulative product $s^{1}s^{2}\cdots s^{n}$. It is observed that $f^{(n)}$ decays more quickly than $s^{1}s^{2}\cdots s^{n}$ as $n$ increases. This is because we haven’t taken into account the decomposed determinants $D^{a}_{[k_{\mu}]}$ and $D^{b}_{[k_{\mu}]}$ in this empirical estimate, and $|D^{a}_{[k_{\mu}]}|,|D^{b}_{[k_{\mu}]}|<1$ holds strictly as proved in Theorem 2.

In other investigated systems such as graphite, \ceLiCoO2, \ceTiO2, and \ceCrO2$\downarrow$, the contribution from the $f^{(2)}$ terms is noticeably less significant, constituting less than $30\%$ of the full spectrum. Their corresponding second largest singular values $s^{2}<0.6$, suggesting that $s^{1}s^{2}<s^{1}$ and hence $f^{(2)}$ would not be as important as the $f^{(1)}$ terms. One severe deviation from this estimate would be the tight-binding 1D chain. The corresponding $s^{1}s^{2}\approx 10.25$, which is comparable to $s^{1}\approx 9.97$. However, the $f^{(2)}$ contribution is tiny, about $6.1\%$ of the full ($f^{(1)}+f^{(2)}$) spectrum. Again, this is because the decomposed determinants $D^{a}_{[k_{\mu}]}$ and $D^{b}_{[k_{\mu}]}$ are missing. In this regard, it would be better to view the accumulative product $s^{1}s^{2}\cdots s^{n}$ as the upper bound of the $f^{(n)}$ contribution. In other words, if $s^{1}s^{2}\cdots s^{n}$ is significantly small compared to the leading order terms ($s^{1}$ and $s^{1}s^{2}$), then it is already safe to neglect the $f^{(n)}$ terms.

To define the smallness of the cumulative product, we may introduce the ratio $\eta_{n}$ such that: $\eta_{n}\max\limits_{m}\{s^{1}s^{2}\cdots s^{m}\}=s^{1}s^{2}\cdots s^{n}$, where $\max\limits_{m}\{s^{1}s^{2}\cdots s^{m}\}$ is the maximum of the cumulative product. According to the examples studied in this work, we find one could safely neglect the $f^{(n)}$ configurations with $\eta_{n}<0.5$.

Although we have gone beyond $f^{(2)}$ configurations in the determinant calculation, we still find the spectra for all the cases studied in this work can be well defined by just the configurations up to $f^{(2)}$. This applies even to \ceLi, in which up to $f^{(5)}$ configurations are included. The $f^{(2)}$ configurations slightly increase intensities of the two absorption humps at $61$ and $64$ eV compared to the absorption peak at $56$ eV. And the spectra beyond $f^{(2)}$ (summed from $f^{(1)}$ to $f^{(n)}$ for $n=3,4,5$) can basically be reproduced by scaling up the $f^{(1)}+f^{(2)}$ spectrum using an energy-independent factor.

Here, we take a closer examination of the intensity contribution from $f^{(2)}$ configurations (as shown by black dotted curves in Fig. 4). The $f^{(2)}$ contributions in the 1D chain, graphite, and \ceLiCoO2 (the first row of Fig. 4) are different from the ones as in \ceLi, \ceRuO3, and \ceCu (the third row of Fig. 4). In the former group, the $f^{(2)}$ contribution does not have as many absorption features compared to $f^{(1)}$. In particular, the $f^{(2)}$ contribution does not have a well-defined peak at the absorption edge as manifests in $f^{(1)}$. In the latter group, the $f^{(2)}$ contribution mimics the $f^{(1)}$ contribution in that they have similar absorption structures (peaks). All the absorption peaks that appear in $f^{(1)}$ also appear in $f^{(2)}$, including the peaks at the absorption onset. A primary reason for the $f^{(2)}$ contribution mimicking the $f^{(1)}$ one in \ceLi, \ceRuO3, and \ceCu is due to low-energy e-h pair production in these metallic systems. There are many orbitals of similar energies near the Fermi level in the final-state system (the supercell that contains an impurity for modeling core-excited states). More speficially, there are more than 12 orbitals within an energy window of 0.3 eV at the $\Gamma$-point of the supercell, as found by the $\Delta$SCF calculations. Some of these orbitals are occupied and some are not. The majority of the $f^{(2)}$ configurations in these systems can be understood as a transition from the core level to one empty final-state orbital $c$, coupled with a low-energy e-h pair (as labeled by $c_{1}$ and $v_{1}$), whose excitation energy is negligibly small. Therefore the spectral contribution from $f^{(2)}$ mimics the one from $f^{(1)}$, which is only defined by empty final-state orbitals, in a one-body manner. This partially explains why the one-body final-state method can often reproduce the XAS of many metallic systems adequately, even though the physics of MND theory is completely missing.

The $f^{(2)}$ contributions in the 1D chain, graphite, and \ceLiCoO2 are more complex to analyze, which involves not just low-energy e-h pairs near the Fermi level / band gap. In graphite, for instance, many $f^{(2)}$ configurations that contribute modestly to XAS are composed of a core-level transition to $c$ which is coupled to an e-h pair with $c_{1}$ being a low-lying empty orbital ($\pi^{*}$), and $v_{1}$ going over the continuum of occupied orbitals that span more than $8$ eV ($\pi$ and $\sigma$ continuum). So the transition energy should be added from the two e-h pairs: $E=(\varepsilon_{c}-\varepsilon_{h})+(\varepsilon_{c_{1}}-\varepsilon_{v_{1}})$, in which both $\varepsilon_{c}-\varepsilon_{h}$ and $\varepsilon_{c_{1}}-\varepsilon_{v_{1}}$ can vary across a wide energy range. This explains why the $f^{(2)}$ spectral contribution in graphite is smeared out without well-defined peaks. Similar analysis applies to the $f^{(2)}$ continuum between $530$ and $534$ eV in \ceLiCoO2, which involves transitions from the core level and some valence orbitals $v_{1}$ (mainly \ceO $2p$) around $4$ eV below the valence band maximum (VBM) to two low-lying empty orbitals $c$ and $c_{1}$. This $f^{(2)}$ continuum fills out the energy gap between $530$ and $534$ eV, and explains why there is no spectral gap between near-edge peaks ( $<6$ eV from onset) and the high-energy humps ($>6$ eV from onset) in the \ceO $K$ edge of TMOs, although there is no single-body orbital within this gap. Such is true for the energy gap near $6$ eV in \ceRuO3.

It should be noted that configuration interaction is absent from the above analysis of $f^{(2)}$ configurations. Interaction within $f^{(2)}$ may introduce excitonic effects between $c_{1}$ and $v_{1}$ and plasmon excitations. So far there is no electron-plasmon (plasmaron) coupling present in the determinant approach. Configuration interaction between $f^{(1)}$ and $f^{(2)}$ may remix the spectral contribution from the two sets of configurations and modify the shakeup effects. So far the ratio of the intensities of the near-edge peaks to the high-energy humps in TMOs is still too high in the determinant calculation up to $f^{(2)}$, as compared with experiments. Introducing configuration interaction may resolve this problem, which is beyond the scope of this work.

II.4 Initial-state, final-state, and the determinant XAS

Lastly, we rationalize the determinant calculation by comparing the determinant spectra with the initial- and final-state spectra (even rows in Fig. 4).

For the 1D tight-binding chain, the MND effects manifest as an asymmetric, power-law singularity that diverges at the absorption edge. Note that this singularity has already been reproduced at the $f^{(1)}$ level. Thanks to the final-state orbitals in the determinant formalism, one does not need to go over many orders of Feynman diagrams expanded in the initial-state orbitals to produce the edge singularity. The determinant spectrum does not resemble to the one-body final-state spectrum, which looks like a ”half dome”, making the inclusion of MND effects essential to the spectral calculation. The singularity at $2$ eV is the van Hove singularity due to the 1D band edge.

Similar MND effects also manifest in the XAS (polarization vector is $45\deg$ off-plane) of graphite. After the correction of the determinant approach, the first-peak (around $286$ eV, due to $\pi^{*}$) intensity is significantly magnified compared with the final-state XAS. The corrected intensity ratio of the first peak to the second (around $293$ eV) is $1.25$ ($0.65$ in the final-state spectrum), which is in good agreement with $1.35$ in a previous experiment weser2010induced . The spectral plateau between $288$ and $292$ eV (due to the constant joint DOS in 2D systems) is also tilted upward at the low-energy end, due to the MND effects. Such MND effects in graphite were also obtained from first-principles using a more complex approach based on Green’s function, which involves energy and time integral, and Fourier transformation wessely2005ab ; wessely2006dynamical ; ovcharenko2013specific . Here, the determinant approach provides equivalent spectra by only a one-shot matrix inversion and a heuristic search of computationally accessible configuration space.

Significant intensity correction of the peak at onset is also observed in \ceLiCoO2, \ceLi, and \ceRuO3. In \ceLiCoO2, the intensity ratio of the peak at $528$ eV to the peak at $532$ eV as found by the determinant approach is $2:1$, which is in good agreement with a previous measurement juhin2010angular . In the same work, the intensity ratio by the one-body final-state approach with GGA + U is $1:1$, which is consistent with our calculation. In \ceRuO3, the intensity ratios of two near-edge peaks are reversed after correction, which was also reported in our previous work on 3d TMOs liang2017accurate ; liang2018quantum .

The Li $K$ edge calculated from the determinant formalism (at the level of $f^{(4)}$) for \ceLi metal is in good agreement with a previous experiment miedema2014in_situ . Thermal vibration of the lattice could further modify the spectrum pascal2014finite , and may broaden the peaks at 61 and 63 eV, making it in closers agreement with the experiment. However, thermal effects will not be discussed in this work, because we focus on the convergence with respect to shakeup orders.

The core-hole attraction effect is most significant in the 1D Chain, graphite, and \ceLi, although they are gapless systems. After the inclusion of the core hole, the initial-state spectra dramatically redshift. However, the core-hole effect only causes the initial-state spectra to redshift rigidly in \ceLiCoO2, \ceRuO3, and \ceCu. The redshifts of the lowest-energy peaks in \ceLiCoO2 and \ceRuO3 are $1.28$ and $0.34$ eV respectively, which are not negligibly small, although the core hole is at the \ceO site where the near-edge orbitals are mainly composed of TM $d$ orbitals. This explains why initial-state spectra are sometimes good approximation to XAS for \ceO $K$ edges in TMOs, as discussed in Ref. de2008core .

III Conclusions

In summary, we have introduced two theorems for regulating the convergence of the determinant calculation, using a SVD analysis over the auxiliary $\zeta$-matrix. The convergence with respect to the excitation order $n$ depends on the number of the significant singular values of the $\zeta$ matrix. We show that the cumulative product of the singular values $s^{1}s^{2}\cdots s^{n}$ can be used as an effective estimate for the $f^{(n)}$ contribution. It is found empirically that it is safe to neglect the $f^{(n)}$ contribution and higher-order when $s^{1}s^{2}\cdots s^{n}<0.5\max\limits_{m}\{s^{1}s^{2}\cdots s^{m}\}$. However, satisfactory determinant spectra have been achieved at the order of $f^{(2)}$ for all the examined cases (in this work and the TMOs as in Ref. liang2018quantum ).

Acknowledgement

Theoretical and computational work was performed by Y. L. and D. P. at The Molecular Foundry, which is supported by the Office of Science, Office of Basic Energy Sciences, of the United States Department of Energy under Contact No. DE-AC02-05CH11231. Computations were performed with the computing resources at the National Energy Research Scientific Computing Center (NERSC).

Appendix A Computational Details

The single-body energies and orbitals for both the initial- and final-state systems are obtained with DFT $\Delta$ SCF calculations, as described in Ref. liang2017accurate ; liang2018quantum . Modified pseudopotentials (by changing $1s^{2}$ to $1s^{1}$) are generated for simulating the core-hole potential of $K$ edges. Supercell dimensions are chosen to be approximately $10$Å that is sufficient to minimize spurious periodic interactions among the core-hole impurities. The $\Delta$ SCF calculations for TMOs are performed using the DFT+$U$ theory dudarev1998electron with the $U$ value adopted from Ref. wang2006oxidation . A $5\times 5\times 5$ $\bm{k}$-grid is used to sample the BZ of the supercell, which is essential for high-energy scattering states and metallic systems. The DFT part of the calculations is performed using a local repository of the ShirleyXAS code, which is available at the David Prendergast’s group at the Molecular Foundry.

The determinant calculations are performed using an open source software package, MBXASPY, which is available at https://github.com/yufengliang/mbxaspy.