Quantum ManyBody Effects in XRay Spectra Efficiently Computed using a Basic Graph Algorithm
Abstract
The growing interest in using xray spectroscopy for refined materials characterization calls for accurate electronicstructure theory to interpret xray nearedge fine structure. In this work, we propose an efficient and unified framework to describe all the manyelectron processes in a Fermi liquid after a sudden perturbation (such as a core hole). This problem has been visited by the MahanNoziéresDe Dominicis (MND) theory, but it is intractable to implement various Feynman diagrams within firstprinciples calculations. Here, we adopt a nondiagrammatic approach and treat all the manyelectron processes in the MND theory on an equal footing. Starting from a recently introduced determinant formalism [Phys. Rev. Lett. 118, 096402 (2017)], we exploit the lineardependence of determinants describing different final states involved in the spectral calculations. An elementary graph algorithm, breadthfirst search, can be used to quickly identify the important determinants for shaping the spectrum, which avoids the need to evaluate a great number of vanishingly small terms. This search algorithm is performed over the treestructure of the manybody expansion, which mimics a pathfinding process. We demonstrate that the determinantal approach is computationally inexpensive even for obtaining xray spectra of extended systems. Using KohnSham orbitals from two selfconsistent fields (ground and coreexcited state) as input for constructing the determinants, the calculated xray spectra for a number of transition metal oxides are in good agreement with experiments. Manyelectron aspects beyond the BetheSalpeter equation, as captured by this approach, are also discussed, such as shakeup excitations and manybody wave function overlap considered in Anderson’s orthogonality catastrophe.
I Introduction
There is a fastgrowing interest in using firstprinciples computational methods to interpret xray spectroscopies for characterizations of materials and thereby enhance our basic understanding of electronic structure haerle2001sp ; nilsson2004chemical ; wernet2004structure ; prendergast2006x ; debeer2010calibration ; woicik2007ferroelectric ; rehr2009ab ; zhao2011visualizing ; tan2012unraveling ; liu2012phase ; drisdell2013probing ; eelbo2013adatoms ; t2013compensation ; velasco2014structure ; pascal2014x ; mcdonald2015cooperative ; wernet2015orbital ; lu2016quantitative ; drisdell2017determining ; liu2017highly . Fulfilling this task requires a reliable prediction of possible atomic structures that could lead to the observed spectra, and more challengingly, a generic theory that can predict accurate xray spectral fingerprints for given systems. Central to a firstprinciples spectroscopic theory is solving the dynamics of a manyelectron Hamiltonian upon excitation of a core electron by an xray photon, for realistic systems ranging from molecules to solids, in an efficacious manner.
From a fundamental viewpoint, the approaches to tackle a manybody problem fall into two major categories. Quantumfieldtheoretical methods wen2004quantum ; bruus2004many ; tsvelik2007quantum ; mattuck2012guide focus on describing the trajectories of a manybody system. Through computing the path integrals of all trajectories from one manybody state to another, one obtains the transition probability between the two. The fieldtheoretical approach has given rise to a set of powerful firstprinciples tools such as the and BetheSalpeterEquation (BSE) method hybertsen1986electron ; rohlfing2000electron ; onida2002electronic ; vinson2011bethe . In current implementations of these methods, only a finite set of diagrams are incorporated, due to the daunting complexity of evaluating all them. The other category of approaches focuses on the description of manybody wave functions based on Slater determinants knowles1984new ; szabo2012modern ; shao2006advances . This leads to methods that are used prevalently in quantum chemistry such as the full configuration interaction (FCI) approach and the coupledcluster technique pople1987quadratic ; bartlett2007coupled ; booth2013towards , or exact diagonalization for solving stronglycorrelated systems caffarel1994exact ; dagotto1994correlated . Currently these methods are mostly applied to systems with electrons limited by the exponential growth of the configuration space.
For xray excitations and associated spectra, we have witnessed the success of the constrainedoccupancy density functional theory (SCF) taillefumier2002x ; prendergast2006x ; drisdell2013probing ; pascal2014x ; velasco2014structure ; ostrom2015probing ; drisdell2017determining , which approximates an xray excited state with one empty KohnSham (KS) orbital in the final state. Recently, we highlighted the shortcomings in this singleparticle (1p) approach for a class of transition metal oxides (TMOs) and indicated an lack of generation to higherorder excitations involving multiple electronhole (eh) pairs liang2017accurate . Driven by these deficiencies, we proposed a better manybody wavefunction ansatz that approximates the initial and final state with a single Slater determinant. The initialstate Slater determinant is constructed from the KS orbitals of the groundstate system, while the Slater determinant for a specific finalstate is derived from the KS orbitals of the coreexcited system. Within this approximation, the transition amplitude can also be expressed as a determinant anderson1967infrared ; dow1980solution ; stern1983many ; ohtaka1990theory ; liang2017accurate comprising transformation coefficients between the two KS basis sets. We find this determinant approach can rectify the deficiency of the 1p SCF approach for a few TMOs liang2017accurate . It is natural to ask: (a) does this formalism provide a good approximation for xray nearedge structures in general? (b) is it practicable for calculations of extended systems, given the huge configuration space? (c) can this approach permit access to higherorder excitations and describe various manybody xray spectral features beyond the BSE?
In this work, we answer these questions by demonstrating an efficient yet simple approach to explore the large configuration space in the determinant formalism. A crucial first step is to relate similar determinants to one another via exterior algebra yokonuma1992tensor ; winitzki2010linear and then evaluate them via updates, rather than from scratch. Even so there are still to manybody states to consider for configurations with double eh pairs . However, only a small portion of these determinants have significant transition amplitudes, due to the spatially localized nature of corelevel excitations, as can be tested by bruteforce calculations. Motivated by this observation, we adopt a breadthfirst search (BFS) algorithm knuth1998art ; cormen2009introduction to look for nontrivial configurations rather than exhausting the entire configuration space.
The BFS algorithm is a basic algorithm for traversing a tree structure, finding the shortest path skiena1990dijkstra ; zeng2009finding , solving a maze moore1959shortest , and other combinatorial search problems. Although the BFS algorithm cannot guarantee answers within a polynomial time, substantial speedup can often be achieved via heuristically pruning the search tree pearl1984heuristics ; zhou2006breadth . For the manybody configuration problem, we design the BFS to search for active “pathways” from the initial state to many excitedstate configurations. Instead of directly accessing a large number of highorder configurations, the search algorithm first visits its ascendant configurations with fewer eh pairs. If multiple pathways to an ascendant configuration interfere destructively and result in a small transition amplitude, the search algorithm will discard the configuration before more highorder configurations are generated. We will show that this treepruning technique can typically lead to at least 100fold speedup in the calculation of xray spectra. This search algorithm is generic and can be generalized to any kind of sudden perturbation.
The determinant formalism is an exact solution to the MahanNoziéresDe Dominicis (MND) model mahan1967excitons ; nozieres1969singularities in which multiple electrons interact with a core hole. Hence, this approach can naturally incorporate all manyelectron processes in the MND theory, which includes the direct and exchange diagrams as in the BSE rohlfing2000electron ; onida2002electronic ; vinson2011bethe , the zigzag diagrams, and the diagrams with a core hole dressed by many eh bubbles. While the BSE diagrams mainly describe eh attraction, or excitonic effects, the zigzag or bubble diagrams describe higher order eh excitations that lead to shakeup features brisk1975shake ; mcintyre1975x ; ohtaka1990theory ; enkvist1993c ; calandra2012k ; mahan2013many ; lemell2015real or manybody effects due to reduced wavefunction overlap. A reduction in manybody wave function overlap is the origin of the Anderson orthogonality catastrophe anderson1967infrared ; mahan2013many . If one were to include all of these effects using the diagrammatic approach, a comprehensive set of techniques, such as solving BSElike equations and using a cumulant expansion kas2015real ; kas2016particle , would be required. Here, the determinant formalism, in conjunction with the firstprinciples KS orbitals, provides a efficient means to investigate all manyelectron effects within the MND model rigorously, for a wide energy range, within a simple unified framework.
This new determinant formalism has already shown great practicality to address realistic problems in materials characterization. We systematically study the \ceO edge () xray absorption spectra (XAS) of various TMOs and find this approach can faithfully reproduce the experimental xray line shapes for most of the investigated systems. This can be immediately applied to study various energy conversion and storage systems involving oxides yabuuchi2011detailed ; hu2013origin ; suntivich2011design ; lin2016metal ; luo2016charge ; strasser2010lattice ; matsukawa2014enhancing ; lebens2016direct ; de2016mapping , where the interpretation of xray spectra can be challenging, and the conclusions often depend sensitively on intricate nearedge line shapes.
This work is organized as follows. Sec. II.1 revisits the manybody effects captured by the MND theory in terms of Feynman diagrams. Sec. II.2 and II.3 provides a solution to the MND model from the perspective of manyelectron wave functions and introduces the determinant formalism. Sec. II.4 introduces exterior algebra to elucidate the linear dependence of the determinants that is encoded in the socalled matrix, followed by a BFS algorithm for an efficient evaluation in Sec. II.5. Sec. II.6 discusses how to combine this algorithm with DFT simulations and its validity in the presence of ee interactions. The simulated XAS of a variety of oxides are shown in Sec. III, with analysis of spectra obtained from different level of approximations. The manybody aspects beyond the BSE as captured by this method will be discussed in Sec. III.4 and III.5, using the halfmetal \ceCrO2 as an example. Finally, the numerical details and efficiency of the newly introduced algorithm are analyzed in Sec. IV.
Ii Theoretical Models and Methods
ii.1 Independentelectron model and diagrammatic approaches
We first revisit the conceptually simple MND model from the perspective of Feynman diagrams. The incorporation of firstprinciples calculations will be deferred to Sec. II.6. In the MND model nozieres1969singularities ; ohtaka1990theory ; mahan2013many , the electrons only interact with the core hole and electronelectron (ee) interactions are neglected. Consider a supercell with one of the atoms replaced by its coreexcited version. This is typically a good approximation to a coreexcited system at low photon flux. Assume there are valence electrons in its ground state and there is only one core level. The MND Hamiltonian without ee interactions reads
(1) 
where the diagonal part is composed of the valence orbitals ( iterates over both occupied and empty valence orbitals) and the core level (). and are electron and hole creation operators respectively. The only twobody term in is the Coulomb interaction between the valence orbitals and the core level, as described by , in which the corehole potential is defined by
(2) 
where ’s are the 1p wave functions and is the (effective) Coulomb potential. The twobody interaction accounts for the electron scattering from orbital to due to the corehole potential.
The xray photon field can be described by a current operator mahan2013many that promotes one core electron to a valence orbtial
(3) 
The transition operator is the electric field polarizationprojected position operator that couples the core level to valence orbitals: , in the limit of zeromomentum transfer and within the dipole approximation prendergast2006x ; de2008core . In principle, the transition operator can be any other local sudden perturbation, not necessarily limited to a core hole.
The independentelectron model was originally considered by the MND theory mahan1967excitons ; nozieres1969singularities ; mahan2013many using diagrammatic techniques. The timeevolution of the manyelectron system after photon absorption is described by the Kubo currentcurrent correlation function
(4) 
where is the vertex that represents the absorption of a photon to create an eh pair ( represents the opposite process). The xray absorption spectrum (XAS) is the spectral function of the photon selfenergy in the frequency domain
(5) 
In the following discussion, we focus on the eh correlation function as defined in Eq. (4)
(6) 
which includes all the manyelectron processes in xray absorption.
We exemplify these manyelectron processes by four types of secondorder Feynman diagram of , as shown in Fig. 1. The time axis runs from left to right and the Coulomb lines are vertical due to the neglect of dynamical effects in the Coulomb interaction . The BSE captures two kinds of processes: direct eh attraction as described by the ladder diagram Fig. 1 (a), and eh exchange as described by the diagram in Fig. 1 (b). In these diagrams, there is only one eh pair present at any time of the propagation. However, there are other diagrams with more eh pairs present at a time, e.g., the zigzag diagram in Fig. 1 (c). The corresponding process involves a core hole causing the ground state to decay into a valence eh pair ( and ) at . At a later time , the core hole assists the newly generated valence hole () to recombine with incoming electron (), leaving an outgoing electron () and the core hole. Lastly, it is also possible that the valence eh pair ( and ) generated earlier does not correlate with the incoming electron at all and simply annihilates at a later time . This leads to a bubble diagram with a freely propagating electron and a core hole dressed by eh bubbles as shown in Fig. 1 (d). These eh bubbles tend to reduce the manybody wave function overlap and are the causes for the Anderson orthogonality catastrophe anderson1967infrared .
The MND theory mahan1967excitons ; nozieres1969singularities ; mahan2013many systematically studies and estimates the impact of these diagrams on the nearedge structure of xray spectra. In essence, it is found that denominators in the BSE diagrams involve , which is roughly the energy required to create an electroncorehole excitation, while the denominators in the zigzag or bubble diagrams involve an offset of , the energy required to create an additional (valence) eh pair. This means the zigzag processes or the bubble diagrams can become significant in a metallic system where can be vanishingly small, or if the photon energy is sufficiently high to be in resonance with double eh excitations.
In practical firstprinciples BSE calculations however, ee interactions are taken into account and the bare Coulomb interactions in the direct (Fig. 1 (a) ) diagram are replaced by screened Coulomb interactions. The screened Coulomb interactions are typically modeled with the emptybubble diagrams within the randomphase approximation rohlfing2000electron ; onida2002electronic ; vinson2011bethe , which, to some extent, describe the manyelectron screening effects in xray excitations.
ii.2 An alternative MND solution based on manybody wave functions
In the last section we have discussed the diagrammatic approach, or manybody perturbation theory (MBPT), for solving the MND model [Eq. (1)]. However, this Hamiltonian is essentially quadratic and exactly solvable. For the initial state, no core hole is excited and and hence the initialstate Hamiltonian is simply For the final state, there is exactly one core hole, i.e., , and the finalstate Hamiltonian also becomes quadratic
(7) 
Within the quadratic forms, it is straightforward to construct the manybody wave functions of the initial and finalstate. The initial state is simply a Slater determinant that consists of valence electrons occupying the lowestlying orbitals:
(8) 
where goes over all the occupied valence orbitals, annihilates the core hole (fills the core level with one electron), and is the null state with no electrons. The finalstate XAS wave functions can be expressed in a similar manner, but using the eigenvectors of
(9) 
where the index is a tuple: , which denotes the valence orbitals that the electrons will occupy in the final state. (with tilde) correspond to the eigenvectors of so that . To apply the Fermi’s Golden rule, one needs to work within the same basis set. We express the finalstate basis set in terms of the initialstate one:
(10) 
where ’s are the transformation coefficients: .
With these expressions for and , the manybody transition matrix element for any onebody operator has been calculated in previous work stern1983many ; ohtaka1990theory ; liang2017accurate
(11) 
in which the transition amplitude also takes a determinantal form
(12) 
The row index goes over occupied finalstate orbitals , the column index over the lowestlying initialstate orbitals plus one empty orbital labeled by (This empty orbital is coupled to the core level with the onebody operator ). This determinantal form reflects how these electrons transit from the initial to final state in the xray excitation process. All the possible electronic pathways are taken into account by the transformation matrix in . The transition amplitude of an individual electron is quantified by the matrix elements, i.e., the initialfinal orbital overlap . The interference of these pathways is lumped into a determinant due to the fermionic nature of electrons.
For the quadratic , the energy of a finalstate can be obtained by direct summation of 1porbital energies
(13) 
where are taken from the diagonalized . A relative energy may also be defined for later discussion, where is the energy of the lowestlying : .
For ease of calculation, previously we have also regrouped the finalstate multielectron configurations according to the convention in quantum chemistry dow1980solution ; szabo2012modern ; bartlett2007coupled . The configuration with is dubbed as a single or a configuration because it has one electron(core)hole pair. A shorthand notation for an configuration can be employed, using to denote the the orbital of the excited valence electron. with and is dubbed as a double or configuration because it has one extra (valence) eh pair as defined by the electron (hole) index (). The shorthand notation for is . This definition can be extended to higher orders such as triples and so forth. For unique indexing, we require and in a index. Examples of finalstate are shown in Fig. 2 (schematics on the second row).
ii.3 Interpretation of the finalstate manybody approach from an initialstate perspective
In this section, we provide a comparison between the outlined determinant formalism and MBPT using Feynman diagrams. While the determinant formalism constructs manyelectron states using both initial and finalstate orbitals, MBPT, such as BSE, relies on initialstate quantities only. To relate the two theories, we can express the MND manyelectron final states in Eq. (9) using only the initialstate orbitals. Rewriting finalstate operators according to a linear combination of the initialstate operators [Eq. (10)] and expressing the wave function in terms of :
(14) 
Expanding the product of the operators and regrouping like terms,
(15) 
The leadingorder term comprises linear combinations of single electron(core)hole pairs, because there are creation operators and destruction operators in Eq. (14), leaving at least one creation operator for an unoccupied state. For this term, out of indices are chosen from so that ’s can cancel with ’s. There are such permutations, and reordering the fermionic operators gives rise to the determinantal form of the coefficients, as previously stated in Eq. (12).
The next term in Eq. (15) is a double term , which has one additional valence eh pair generated on the top of the electroncorehole pair. This term takes into account the secondorder manyelectron processes: the valence eh excitations induced by the corehole potential, which are also known as the shakeup excitations brisk1975shake ; mcintyre1975x ; ohtaka1990theory ; enkvist1993c ; calandra2012k ; mahan2013many ; lemell2015real , because an additional amount of energy is required to create these valence excitations. As the series expansion proceeds, each term will have one more valence ehpair than the last, and more complicated shakeup processes with multiple ehpairs are included. A full schematic for the relation of one single finalstate configuration (written as one Slater determinant using finalstate orbitals) in terms of initialstate configurations is shown in Fig. 2.
Within MBPT, the configuration series in Eq. (15) is typically truncated, and the coefficients are solved by expanding the Hamiltonian over the restricted configuration space and solving the eigenvalue problem. In the BSE, for instance, the finalstate Hamiltonian is expanded over the singleehpair space and the eigenvector coefficients (analogous to ) refer to this singleeh basis. In some sense, this approximation corresponds to the ladder and exchange diagrams: at any point in time of the propagation, there is only one eh pair involved.
By contrast, the determinant formalism does not restrict the number of ehpairs in the finalstate configuration space. When is projected onto as in Eq. (15), a superposition of single, double, and highorder terms naturally arises, although only the leadingorder coefficients are relevant for calculating matrix elements of onebody operator. In this way, the zigzag and bubble diagrams, present within MND theory, which involve multiple ehpair generation, are automatically incorporated.
ii.4 Efficient evaluation of determinantal transition amplitudes
The above determinantal formalism provides an alternative solution to the MND model in Eq. (1) without using diagrammatic approaches. If a sufficient number of final states are included, one may expect the determinantal method to give the spectrum as solved from the MND model. However, an bruteforce calculation is rarely used because the manyelectron configuration space grows factorially with the number of electrons. It does not seem to be practical to compute the large number of determinants that would represent all configurations.
For a halffilled system with orbitals and () electrons, even the group has configurations. Iterating the index of [Eq. (12)] over all empty initialstate orbitals multiplies the time complexity by a factor of . Calculating the determinant for each configuration requires a computational cost of . With all the three factors combined, obtaining the determinants for all of the configurations gives rise to a time complexity of . For metallic systems where the fermi surfaces are susceptible to the corehole potential, higherorder terms such as are typically needed for testing convergence, which leads to a higher time complexity of . Such a bruteforce calculation that scales up quickly with number of states is not very practical for realistic corehole calculations in which there could easily be to orbitals.
In this section, we introduce an efficient algorithm at much lower computational cost to access the determinants that are important for determining the xray spectrum. The determinant calculation needs to be performed only once for a given configuration, and subsequently the determinants for other configurations can be derived from it. More importantly, a BFS algorithm is employed to identify the important determinants above a specified threshold, largely reducing the number of configurations to be visited.
An apparent first step is to move the summation over in Eq. (11) into the definition of the transition amplitude coefficient, so that for each finalstate configuration , obtaining requires calculating only one determinant. More specifically, we rewrite as
(16) 
where . The summation in the column of can be calculated first before obtaining the determinant. This reduces the overall time complexity by a factor of .
Secondly, when considering transitions to various finalstate configurations, the determinants of interest in fact have many common rows so one can make use of the multilinearity of determinants to speed up the calculations significantly. For example, the tuple for a double configuration only differs from the groundstate one by 3 indices, meaning their corresponding determinants only differ by 3 rows. This observation motivates us to choose the determinant for the ground state as a reference, and evaluate other determinants for excited states via a lowrank updating technique.
To demonstrate this technique, it is most transparent to express the determinant in terms of the wedge (exterior) product yokonuma1992tensor ; winitzki2010linear of its row/column vectors. The wedge product is anticommutative and has similar algebra to the Fermionic operators. Suppose an arbitrary matrix has row/column vectors , its determinant can be expressed as
(17) 
Assume has been calculated from scratch and is nonzero (assume is fullrank). If is replaced by a new vector , which can be considered as a rank1 update, the updated determinant can be obtained by expanding in terms of
(18) 
where is the expansion coefficient defined as
(19) 
can be obtained via the matrix inversion of : . When multiplied by , only survives in the summation because . Then the new determinant is simply the product of an expansion coefficient and the alreadyknown .
Now if the last two lines of are replaced by two new row vectors and , the rank2 updated determinant is
(20) 
The minus sign arises from the anticommutative property of the wedge product: . Thus the new determinant is the product of a determinant composed of the expansion coefficients and . The above procedure can be carried out to more general situations where more row/column vectors are replaced. This remove the need to calculate the new determinant from scratch using the algorithm. For a rank update, one only needs to compute the product of the reference determinant and a small determinant containing , at the cost of .
In the context of the determinantal formalism as in Eq. (16), we define the row vector corresponding to the finalstate orbital as:
(21) 
Then the groundstate reference determinant can be expressed as . To access the determinants for excited states via this updating method, we formally introduce the auxiliary matrix () for a system with M orbitals and N valence electrons (), which is the transformation matrix from to :
(22) 
Rewrite the above matrix multiplication in a compact form, we have , where . Then can be obtained easily via matrix inversion and multiplication
(23) 
Note that is a matrix and we find it is typically invertible in practical calculations.
is of size . Its column indices map onto the lowestlying orbitals while the row indices map onto the to orbitals. A determinant can be obtained from the product of , a minor of , and an overall sign due to permutation of rows (trivial to consider in the singledeterminant case). The last column of this minor must be taken from the last (the ) column of , because there are electron orbitals and hole orbitals in a configuration, and the extra one electron can be viewed as removed from the hole on the orbital. The minor reflects the interference effect of pathways to access the configuration via permuting empty orbitals. The rows (columns) of the minor indicate the electrons (holes) that are excited in the given configuration: the minor formed by rows () and columns () corresponds to the configuration .
ii.5 Pruning the configuration space using the breadthfirst search algorithm
With this updating technique, we can access the determinants of many configurations without repeatedly carrying out the full determinant calculation for each. However, the number of configurations still grows exponentially with the order . Even the group grows rapidly as , and a system with orbitals may have configurations. The problem now becomes how to efficiently find all the significant minors of at all orders. Enumerating all of these minors will definitely be a hard problem that can not be solved within a polynomial time, and the question is whether it is necessary to visit all of them. In fact, we find that for the systems studied in this work (introduced in Sec. III.1) is sparse, with its nonvanishing elements concentrated in some regions, as will be shown in Sec. IV. A more efficient algorithm should be possible given the sparsity of .
To make the best use of the sparsity of , we investigate its minor determinants in a bottomup and recursive manner. According to the Laplace (cofactor) expansion, an determinant can be expanded into a weighted sum of minors of size . The determinant is nonvanishing only when at least one of these minors is nonvanishing. Physically, this means that a transition to an configuration is only probable when at least one of its parent configurations is probable, otherwise the transition to that configuration is forbidden. Assume that is sparse, and one can keep a short list of nonvanishing minors. When proceeding to order, one can construct the determinant from the short list of nonvanishing minors instead of exhaustively listing all of them.
This recursive construction of minors from the minors leads us to an ultimate improvement to the efficiency of the determinantal approach. We employ the breadthfirst search (BFS) algorithm to enumerate all important minors of . A configuration can be considered as a descendent of via creating one more eh pair with the configuration. Through arranging according to this inheritance relation, a treelike structure of the manybody expansion is formed, as illustrated in Fig. 3. The BFS algorithm visits this treelike structure in ascending order of . Note that a configuration can be accessed from its multiple parents via different pathways. If these pathways to the configuration interfere destructively such that the transition amplitude is vanishingly small, the BFS algorithm will discard this configuration, hence reducing the search space for the next order. Here is the detailed algorithm
Below are further instructions on the lines marked by asterisks.

of can be simply taken from the nonzero matrix elements on the last column of .

is a threshold for small matrix elements. One can set , where and is a userdefined relative threshold.

The determinant of is constructed via a Laplace expansion along its first column. ensures the chosen matrix element is always on the first column of the determinant.

Compare to , contains one more eh pair labeled by and . Because we require the ordering of and for unique indexing, the new index must obey the same order. The new sequence can be obtained by this procedure: first place the new in front of the sequence of that is already increasingly sorted, and then shift to the right by swapping indices till the whole sequence is also sorted. Define to be the number of swaps performed for deciding signs. can be obtained simply by placing at the end of .

is the cofactor of the Laplace expansion of a determinant. where is the proper position for inserting into , as defined above. At the end of the loop, there are at most contributions to the total amplitude of a specific configuration, corresponding to the transition amplitudes of different pathways from its parent configuration.

is a threshold for removing state with small oscillator strengths. Similar to , can be set to , where is the maximal oscillator strength and is a userdefined relative threshold. can be chosen to be the maximal intensity within the group which typically have the strongest oscillator strengths among all groups. can be related to the previously defined relative matrixelement threshold . If the contribution from a small were not added to , its intensity would be . Replacing with will lead to an error of . Therefore, choosing a such that can guarantee error in intensities smaller than . In practice, one can lower till convergence is achieved.
The detailed implementation of this search algorithm can be found at Ref. mbxaspy within an opensource PYTHON simulation package.
Here we demonstrate the BFS algorithm with a toy model with orbitals and valence electrons. Suppose of the system is
(24) 
The BFS algorithm for this example of is carried out as follows:
First, the nonzero configurations are initialized, (5), (6), and (8), whose determinants are simply the matrix elements: , , and , respectively. These configurations are considered as the roots of the BFS trees, as is shown in Fig. 3.
Next, the configurations are constructed based on the obtained configurations. Take the configuration in for example. There are 5 nonzero matrix elements that are to the left of and are not on the same row as , which are , , , , and . Paired up with these matrix elements, the configuration spawns 5 configurations: , , , , and (comma omitted due to the singledigit indices). Likewise, and spawn 6 and 4 configurations respectively.
Both and give rise to and the contributions from the configurations are merged: A(638) = . The two possible pathways are: (1) the core electron is first promoted to orbital and then coupled with the eh pair formed by orbital and ; (2) the core electron is first promoted to orbital and then coupled with the eh pair formed by orbital and . If is vanishingly small ( happens to be close with ) due to the destructive interference of the two pathways, Then will be removed from the list because it cannot contribute to the transition amplitude of any higherorder configuration. When the search process for is completed, nontrivial configurations are found.
Proceeding to the third order, the 13 configurations spawn configurations. Paring with and with both lead to , whose determinant is . Paring () with leads to (the other two pathways are forbidden because ). If and are small numbers such that their product is smaller than the specified threshold , then will be removed from . The above process can be repeated until all new determinants are small enough or no new determinants can be found.
If one bruteforcely enumerates all possible determinants, there are and determinants to examine for the above matrix. By contrast, the BFS algorithm only visits the nontrivial determinants and only and determinants are computed.
ii.6 Incorporating firstprinciples calculations into the determinant formalism
In the above sections, we have demonstrated an efficient solution to the MND model using manybody wave functions for simulating xray transition amplitudes. However, in order to simulate reliable xray spectra without fitting parameters from experimetns, we still need accurate approximations to the initial and final states and their energies. To this end, we rely on DFT calculations to obtain the KS eigenstate energies (for ) and wavefunctions (for both and ) as input for constructing the transformation matrix (Eq. (12)) and computing the energies of manyelectron excited states (Eq. (13)).
For the final state, we employ the standard SCF corehole approach to obtain the KS orbitals and eigenenergies. The coreexcited atom is treated as an isolated impurity embedded in the pristine system, and typical supercell settings for finite england2011hydration and extended drisdell2013probing ; velasco2014structure ; pascal2014x ; liang2017accurate systems can be employed. To simulate a electroncorehole pair, the coreexcited atom is modeled by a modified pseudopotential with a core hole, and an electron is added to the supercell system and constrained to one specific empty orbital. In principle, a SCF iteration needs to be performed for each case of constraint occupancy (for all ) which may lead to an expensive computational cost. As a tradeoff, the electron is only placed onto the lowest unoccupied orbital (), which we have dubbed the excitedstate corehole (XCH) method. After the SCF calculation is done, the KS equation with a converged charge density is used for .
Another important variation of the XCH method is the full corehole (FCH) approach, in which the groundstate occupation without the additional electron is used. The advantage of FCH is that it does not bias towards the lowest excited state and treat all of them on equal footing.
For the initial state, the same supercell as in the final state is used except that the coreexcited atom is replaced by a groundstate atom, using the occupation . A standard DFT calculation can be done to obtain the KS orbitals .
With the KS orbitals and obtained from the SCF corehole calculation, we can compute the orbital overlap integral for computing the determinantal amplitudes. To reduce computational cost, we employ projectoraugmentedwave (PAW) form of ultrasoft pseudopotentials kresse1999ultrasoft ; blochl1994projector ; taillefumier2002x to model electronion interactions. The excited atom potential has deeper energy levels and more contracted orbitals so its PAW construction differ from the groundstate atom. In Appendix A, we derive a formalism to calculate expectation values between of the groundstate system and of the coreexcited system. Initial dipole matrix elements are also evaluated within this PAW formalism blochl1994projector ; taillefumier2002x .
As in typical DFT impurity calculations, some lowlying excited states in the corehole approach could be bound to the coreexcited atom, resembling midgap localized electronic states near an impurity. In this situation, the electronic structure is well described by using a single kpoint (the point) to sample the Brillouin zone (BZ). However, for the purpose of spectral simulations, which include delocalized scattering states well above the band edges, we find that employing point sampling is necessary to improve the accuracy of the calculated line shape. Therefore, we perform the determinantal calculation individually for each kpoint and take the kpointweighted average spectrum as the final spectrum. The band structure and orbitals are interpolated accurately and efficiently using an optimal basis set proposed by Shirley shirley1996optimal ; prendergast2009bloch , whose size is much smaller than a planewave basis. In the Shirley construction, the periodic parts of Bloch wave functions across the first BZ are represented using a common basis which spans the entire band structure. Because there is only one optimal basis to represent the Bloch states for all kpoints, the overlap matrix for every kpoint can be computed quickly as in Appendix B.
After the XAS is calculated by the firstprinciples determinantal approach, the established formationenergy calculation can be adopted to align spectra for coreexcited atoms in different chemical contexts, using the XCH method to determine the excitation energy of the first transition england2011hydration ; jiang2013experimental .
Although there is no valence ee interaction terms in the MND theory, which results in a singledeterminant solution to the manybody wave functions, we argue that this firstprinciples determinantal approach does not entirely neglect valence ee interactions. The selfconsistentfield (SCF) procedure in the DFT updates the total charge density and KS orbitals simultaneously, and hence takes into account some degree of valence electronic screening. That said, the SCF approach should lead to a more realistic equilibrium total charge density for both the ground state and xray excited states, whereas the charge density in MBPT is only treated perturbatively. Moreover, a more accurate charge density may lead to a better approximation to quasiparticle (QP) wave functions. In fact, KS orbitals based on a converged SCF are often employed in MBPT to construct the Green’s functions and compute optical oscillator strengths hybertsen1986electron ; rohlfing2000electron ; onida2002electronic , which is typically a good approximation within the Fermiliquid picture. Finally, corrective DFT (DFT + or DFT with exactexchange functionals) or the selfconsistent approximation bruneval2006effect ; van2006quasiparticle ; kang2010enhanced ; sun2017x can also improve QP energies and wave functions to be used in the determinantal approach. We imagine the manybody effects captured in this framework can be described by the bolded version of the MND diagrams in Fig. 1 in which all the Green’s functions become dressed, and the bare Coulomb lines are replaced by the screened corehole potential described with the chosen exchangecorrelation functional.
ii.7 Comparison with the onebody SCF corehole approach
Before the determinantal formalism, the manybody transition amplitudes in the SCF corehole approach are often approximated with 1p matrix elements
(25) 
where the core orbital is in the initial state while the electron orbital is in the final states, both of which can be taken from DFT calculations. represents the response of rest of the manyelectron system (excluding the electroncorehole pair) due to the core hole, and it is normally assumed to be a constant for ease of calculations. This 1p form of the matrix element implies that: (a) the transition from the initial core level to the final electron orbital occurs instantaneously with the response of many other electrons in the system, with no particular time ordering; (b) the corelevel transition and the manyelectron response are not entangled. This is also the so called sudden or frozen approximation.
We know that from the diagrammatic interpretation of the xray manybody processes in Fig. 1, the photon first decays into an initialstate eh pair instantaneously, and then the other electrons see the corehole potential and begin to relax over a finite period of time. This physical reality can also be seen in the determinant formalism, in which the core hole is only coupled to an initialstate orbital, and the subsequent manyelectron response is described by the determinantal amplitude. So the question is why the simpler 1p matrix element in the frozen approximation still works for a good number of systems in the past.
In this section, we approach this question theoretically by relating the determinantal amplitude to the 1p matrix element. To do this, we first express the determinantal amplitude in terms of is minors (wavefunction overlaps of electron systems, such as ) by Laplace expansion along its last column
(26) 
where are the matrix elements on the last column of as in Eq. (12) and is the minor complementary to . Since in the onebody corehole approach only the terms are summed, we limit our analysis here to the manybody terms and condense the configuration tuple into a single index: . Then the matrix elements can be written as
(27) 
First, for systems with significant band gaps (insulators and semiconductors), we could expect that the overlap of the occupied final state orbitals with the unoccupied initial state orbitals could be quite small. For many orbitals unaffected by the localized corehole perturbation, for example, we might expect the final state occupied orbitals to closely resemble their initial state counterparts, which would render identically zero by orthogonality. Therefore, the sum over in Eq. (27) may only be significant in cases where the transformation matrix indicates mixing of unoccupied initial state character into the occupied final state orbitals, which might easily be the case for orbitals close to the Fermi level in a metal or otherwise openshell system.
The first term in Eq. (27) is more directly relevant to our previous onebody approximation. Here, is the minor of , the transformation matrix without its column and row. It reflects the electron manybody overlap between the initial and final state occupied orbitals and should reflect the extent to which the electron density is modified by the corehole perturbation. Since does not depend on , we can relate it to the manybody prefactor that appears in the finalstate rule of Eq. (25): . Using the completeness relation: , the first term in the expansion of Eq. (27) can be expressed as
(28) 
If it happened that , then this expression would amount to the final state matrix element as defined in the onebody finalstate rule (Eq. (25)). By the same arguments made above, for systems with limited mixing of orbital character across a significant band gap, then we might easily expect orthogonality (zero overlap) between occupied initial state and unoccupied final state orbitals. By the same token, we should be wary of limitations in the onebody approach when this is not the case.
It appears useful to focus on to reveal the role of hybridization in modulating nearedge spectral intensity. To quantify the contribution of the second term in Eq. (28), we introduce the projection spectrum
(29) 
in which the single index sums over all empty finalstate orbitals, and . The matrix element is nothing but Eq. (28) or the first term in Eq. (27) with . However, it is easier to calculate Eq. (28) because summation over all empty orbitals is avoided.
Iii Results and Discussion
iii.1 Applications to transitionmetal oxides
In this section, we discuss an important application of the determinantal approach to computing coreexcited state transition amplitudes, that is, to predict the xray absorption spectra (XAS) for transition metal oxides (TMOs). This is also our original motivation for proposing the determinantal approach liang2017accurate , which can be used to overcome the deficiency of the onebody corehole approach. It has been found for a number of TMOs, that the onebody approach systematically underestimates the intensity of nearedge features at the O edge that correspond to orbitals with hybridization between oxygen character and TM character. This underestimation can prevent reliable interpretations of the Xray absorption spectra for this important class of materials.
We use the newly developed determinantal approach to predict the XAS for eight TMOs: the rutile phase of \ceTiO2, \ceVO2 ( K), and \ceCrO2, the corundum \ceFe2O3, the perovskite \ceSrTiO3, \ceNiO, and \ceCuO. \ceSiO2 is also chosen for a comparative study. Their experimental XAS are extracted from Refs. yan2009oxygen ; koethe2006transfer ; stagarescu2000orbital ; shen2014surface ; zhu2012bonding ; lin2014hierarchically ; jiang2013experimental ; dudarev1998electron ; ma1992soft . The chosen TMOs cover a wide range of electronic and magnetic properties and therefore they are used as benchmark materials for the determinantal approach.
The \ceO edges are investigated here, i.e., the transitions from the \ceO level to shells. For TMOs, the \ceO orbitals are covalently hybridized with the transition metal orbitals, and hence the \ceO edge spectra can serve as an informative and sensitive probe for the electron physics de1989oxygen ; yabuuchi2011detailed ; hu2013origin ; suntivich2011design ; lin2016metal ; luo2016charge ; strasser2010lattice ; matsukawa2014enhancing ; lebens2016direct ; de2016mapping . Moreover, unlike transition metal edges (to transitions), in which atomic multiplet effects split spectral features into many closely space linesde2008core , the \ceO edges can provide a picture of the electronic densityofstates related to the shell more easily interpretable in terms of band theory or effective 1p states.
The angularlyaveraged (except in \ceCrO2, where the polarization is perpendicular to the hard axis) \ceO edge spectra for the chosen TMOs are shown in Fig. 4 (a). The very nearedge part of the spectra, i.e., the spectral features below eV contain the most useful information for material characterization. For these TMOs, the nearedge spectral fine structure exhibits two main peaks corresponding to the splitting of the orbitals into a and an manifold in the (quasi)octahedral crystal field. Our goal is to produce reliably all the spectral features, especially the very nearedge part, so that one can interpret the spectra on a firstprinciples basis. More specifically, we use the ratio of the intensity of the first (lowestlying) peak to that of the second (unless otherwise specified) as a metric for the accuracy of different levels of approximation.
[figure]subcapbesideposition=top
We first calculate the XAS for the chosen compounds using the conventional 1p FCH approach taillefumier2002x ; prendergast2006x ; liang2017accurate described above. A modified pseudopotential generated with the configuration is used for the