Quantum Many-Body Effects in X-Ray Spectra Efficiently Computed using a Basic Graph Algorithm

Quantum Many-Body Effects in X-Ray Spectra Efficiently Computed using a Basic Graph Algorithm

Yufeng Liang yufengliang@lbl.gov The Molecular Foundry, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA    David Prendergast The Molecular Foundry, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA

The growing interest in using x-ray spectroscopy for refined materials characterization calls for accurate electronic-structure theory to interpret x-ray near-edge fine structure. In this work, we propose an efficient and unified framework to describe all the many-electron processes in a Fermi liquid after a sudden perturbation (such as a core hole). This problem has been visited by the Mahan-Noziéres-De Dominicis (MND) theory, but it is intractable to implement various Feynman diagrams within first-principles calculations. Here, we adopt a non-diagrammatic approach and treat all the many-electron 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 linear-dependence of determinants describing different final states involved in the spectral calculations. An elementary graph algorithm, breadth-first 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 tree-structure of the many-body expansion, which mimics a path-finding process. We demonstrate that the determinantal approach is computationally inexpensive even for obtaining x-ray spectra of extended systems. Using Kohn-Sham orbitals from two self-consistent fields (ground and core-excited state) as input for constructing the determinants, the calculated x-ray spectra for a number of transition metal oxides are in good agreement with experiments. Many-electron aspects beyond the Bethe-Salpeter equation, as captured by this approach, are also discussed, such as shakeup excitations and many-body wave function overlap considered in Anderson’s orthogonality catastrophe.

I Introduction

There is a fast-growing interest in using first-principles computational methods to interpret x-ray 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 x-ray spectral fingerprints for given systems. Central to a first-principles spectroscopic theory is solving the dynamics of a many-electron Hamiltonian upon excitation of a core electron by an x-ray photon, for realistic systems ranging from molecules to solids, in an efficacious manner.

From a fundamental viewpoint, the approaches to tackle a many-body problem fall into two major categories. Quantum-field-theoretical methods wen2004quantum ; bruus2004many ; tsvelik2007quantum ; mattuck2012guide focus on describing the trajectories of a many-body system. Through computing the path integrals of all trajectories from one many-body state to another, one obtains the transition probability between the two. The field-theoretical approach has given rise to a set of powerful first-principles tools such as the and Bethe-Salpeter-Equation (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 many-body 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 coupled-cluster technique pople1987quadratic ; bartlett2007coupled ; booth2013towards , or exact diagonalization for solving strongly-correlated systems caffarel1994exact ; dagotto1994correlated . Currently these methods are mostly applied to systems with electrons limited by the exponential growth of the configuration space.

For x-ray excitations and associated spectra, we have witnessed the success of the constrained-occupancy density functional theory (SCF) taillefumier2002x ; prendergast2006x ; drisdell2013probing ; pascal2014x ; velasco2014structure ; ostrom2015probing ; drisdell2017determining , which approximates an x-ray excited state with one empty Kohn-Sham (KS) orbital in the final state. Recently, we highlighted the shortcomings in this single-particle (1p) approach for a class of transition metal oxides (TMOs) and indicated an lack of generation to higher-order excitations involving multiple electron-hole (e-h) pairs liang2017accurate . Driven by these deficiencies, we proposed a better many-body wavefunction ansatz that approximates the initial and final state with a single Slater determinant. The initial-state Slater determinant is constructed from the KS orbitals of the ground-state system, while the Slater determinant for a specific final-state is derived from the KS orbitals of the core-excited 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 x-ray near-edge structures in general? (b) is it practicable for calculations of extended systems, given the huge configuration space? (c) can this approach permit access to higher-order excitations and describe various many-body x-ray 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 many-body states to consider for configurations with double e-h pairs . However, only a small portion of these determinants have significant transition amplitudes, due to the spatially localized nature of core-level excitations, as can be tested by brute-force calculations. Motivated by this observation, we adopt a breadth-first 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 speed-up can often be achieved via heuristically pruning the search tree pearl1984heuristics ; zhou2006breadth . For the many-body configuration problem, we design the BFS to search for active “pathways” from the initial state to many excited-state configurations. Instead of directly accessing a large number of high-order configurations, the search algorithm first visits its ascendant configurations with fewer e-h 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 high-order configurations are generated. We will show that this tree-pruning technique can typically lead to at least 100-fold speed-up in the calculation of x-ray 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 Mahan-Noziéres-De Dominicis (MND) model mahan1967excitons ; nozieres1969singularities in which multiple electrons interact with a core hole. Hence, this approach can naturally incorporate all many-electron processes in the MND theory, which includes the direct and exchange diagrams as in the BSE rohlfing2000electron ; onida2002electronic ; vinson2011bethe , the zig-zag diagrams, and the diagrams with a core hole dressed by many e-h bubbles. While the BSE diagrams mainly describe e-h attraction, or excitonic effects, the zigzag or bubble diagrams describe higher order e-h excitations that lead to shakeup features brisk1975shake ; mcintyre1975x ; ohtaka1990theory ; enkvist1993c ; calandra2012k ; mahan2013many ; lemell2015real or many-body effects due to reduced wave-function overlap. A reduction in many-body 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 BSE-like equations and using a cumulant expansion kas2015real ; kas2016particle , would be required. Here, the determinant formalism, in conjunction with the first-principles KS orbitals, provides a efficient means to investigate all many-electron 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 () x-ray absorption spectra (XAS) of various TMOs and find this approach can faithfully reproduce the experimental x-ray 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 x-ray spectra can be challenging, and the conclusions often depend sensitively on intricate near-edge line shapes.

This work is organized as follows. Sec. II.1 revisits the many-body 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 many-electron 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 so-called -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 e-e 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 many-body aspects beyond the BSE as captured by this method will be discussed in Sec. III.4 and III.5, using the half-metal \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 Independent-electron model and diagrammatic approaches

We first revisit the conceptually simple MND model from the perspective of Feynman diagrams. The incorporation of first-principles calculations will be deferred to Sec. II.6. In the MND model nozieres1969singularities ; ohtaka1990theory ; mahan2013many , the electrons only interact with the core hole and electron-electron (e-e) interactions are neglected. Consider a supercell with one of the atoms replaced by its core-excited version. This is typically a good approximation to a core-excited 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 e-e interactions reads


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 two-body term in is the Coulomb interaction between the valence orbitals and the core level, as described by , in which the core-hole potential is defined by


where ’s are the 1p wave functions and is the (effective) Coulomb potential. The two-body interaction accounts for the electron scattering from orbital to due to the core-hole potential.

The x-ray photon field can be described by a current operator mahan2013many that promotes one core electron to a valence orbtial


The transition operator is the electric field polarization-projected position operator that couples the core level to valence orbitals: , in the limit of zero-momentum 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 independent-electron model was originally considered by the MND theory mahan1967excitons ; nozieres1969singularities ; mahan2013many using diagrammatic techniques. The time-evolution of the many-electron system after photon absorption is described by the Kubo current-current correlation function


where is the vertex that represents the absorption of a photon to create an e-h pair ( represents the opposite process). The x-ray absorption spectrum (XAS) is the spectral function of the photon self-energy in the frequency domain


In the following discussion, we focus on the e-h correlation function as defined in Eq. (4)


which includes all the many-electron processes in x-ray absorption.

Figure 1: Four distinct types of e-h processes in the second-order Feynman diagrams in the MND theory. There are exactly two Coulomb lines (at and ) in each diagram, as marked by vertical dashed lines.

We exemplify these many-electron processes by four types of second-order 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 e-h attraction as described by the ladder diagram Fig. 1 (a), and e-h exchange as described by the diagram in Fig. 1 (b). In these diagrams, there is only one e-h pair present at any time of the propagation. However, there are other diagrams with more e-h 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 e-h 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 e-h 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 e-h bubbles as shown in Fig. 1 (d). These e-h bubbles tend to reduce the many-body 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 near-edge structure of x-ray spectra. In essence, it is found that denominators in the BSE diagrams involve , which is roughly the energy required to create an electron-core-hole excitation, while the denominators in the zigzag or bubble diagrams involve an offset of , the energy required to create an additional (valence) e-h 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 e-h excitations.

In practical first-principles BSE calculations however, e-e 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 empty-bubble diagrams within the random-phase approximation rohlfing2000electron ; onida2002electronic ; vinson2011bethe , which, to some extent, describe the many-electron screening effects in x-ray excitations.

ii.2 An alternative MND solution based on many-body wave functions

In the last section we have discussed the diagrammatic approach, or many-body 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 initial-state Hamiltonian is simply For the final state, there is exactly one core hole, i.e., , and the final-state Hamiltonian also becomes quadratic


Within the quadratic forms, it is straightforward to construct the many-body wave functions of the initial- and final-state. The initial state is simply a Slater determinant that consists of valence electrons occupying the lowest-lying orbitals:


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 final-state XAS wave functions can be expressed in a similar manner, but using the eigenvectors of


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 final-state basis set in terms of the initial-state one:


where ’s are the transformation coefficients: .

With these expressions for and , the many-body transition matrix element for any one-body operator has been calculated in previous work stern1983many ; ohtaka1990theory ; liang2017accurate


in which the transition amplitude also takes a determinantal form


The row index goes over occupied final-state orbitals , the column index over the lowest-lying initial-state orbitals plus one empty orbital labeled by (This empty orbital is coupled to the core level with the one-body operator ). This determinantal form reflects how these electrons transit from the initial to final state in the x-ray 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 initial-final 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 final-state can be obtained by direct summation of 1p-orbital energies


where are taken from the diagonalized . A relative energy may also be defined for later discussion, where is the energy of the lowest-lying : .

For ease of calculation, previously we have also regrouped the final-state multi-electron 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) e-h 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 final-state are shown in Fig. 2 (schematics on the second row).

Figure 2: Definitions of the multi-electron configurations used in the initial(i)- and final(f)-state picture according to the convention in quantum chemistry. A final-state configuration (a single Slater determinant) at the order of can be hybridized from a number of initial-state configurations at multiple orders, as shown by the thick opaque downarrows, which illustrates the spirit of Eq. (9). The solid uparrows in a configuration indicate one possible multi-electronic pathway to access that configuration from the ground-state. The dashed uparrows show the other possible pathway to access the configuration.

ii.3 Interpretation of the final-state many-body approach from an initial-state perspective

In this section, we provide a comparison between the outlined determinant formalism and MBPT using Feynman diagrams. While the determinant formalism constructs many-electron states using both initial- and final-state orbitals, MBPT, such as BSE, relies on initial-state quantities only. To relate the two theories, we can express the MND many-electron final states in Eq. (9) using only the initial-state orbitals. Rewriting final-state operators according to a linear combination of the initial-state operators [Eq. (10)] and expressing the wave function in terms of :


Expanding the product of the operators and regrouping like terms,


The leading-order 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 e-h pair generated on the top of the electron-core-hole pair. This term takes into account the second-order many-electron processes: the valence e-h excitations induced by the core-hole 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 e-h-pair than the last, and more complicated shakeup processes with multiple e-h-pairs are included. A full schematic for the relation of one single final-state configuration (written as one Slater determinant using final-state orbitals) in terms of initial-state 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 final-state Hamiltonian is expanded over the single-e-h-pair space and the eigenvector coefficients (analogous to ) refer to this single-e-h 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 e-h pair involved.

By contrast, the determinant formalism does not restrict the number of e-h-pairs in the final-state configuration space. When is projected onto as in Eq. (15), a superposition of single, double, and high-order terms naturally arises, although only the leading-order coefficients are relevant for calculating matrix elements of one-body operator. In this way, the zig-zag and bubble diagrams, present within MND theory, which involve multiple e-h-pair 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 brute-force calculation is rarely used because the many-electron 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 half-filled system with orbitals and () electrons, even the group has configurations. Iterating the index of [Eq. (12)] over all empty initial-state 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 core-hole potential, higher-order terms such as are typically needed for testing convergence, which leads to a higher time complexity of . Such a brute-force calculation that scales up quickly with number of states is not very practical for realistic core-hole 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 x-ray 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 final-state configuration , obtaining requires calculating only one determinant. More specifically, we rewrite as


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 final-state 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 ground-state 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 low-rank 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


Assume has been calculated from scratch and is nonzero (assume is full-rank). If is replaced by a new vector , which can be considered as a rank-1 update, the updated determinant can be obtained by expanding in terms of


where is the expansion coefficient defined as


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 already-known .

Now if the last two lines of are replaced by two new row vectors and , the rank-2 updated determinant is


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 final-state orbital as:


Then the ground-state 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 :


Rewrite the above matrix multiplication in a compact form, we have , where . Then can be obtained easily via matrix inversion and multiplication


Note that is a matrix and we find it is typically invertible in practical calculations.

is of size . Its column indices map onto the lowest-lying 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 single-determinant 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 breadth-first 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 non-vanishing 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 bottom-up 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 non-vanishing only when at least one of these minors is non-vanishing. 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 non-vanishing minors. When proceeding to order, one can construct the determinant from the short list of non-vanishing 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 breadth-first search (BFS) algorithm to enumerate all important minors of . A configuration can be considered as a descendent of via creating one more e-h pair with the configuration. Through arranging according to this inheritance relation, a tree-like structure of the many-body expansion is formed, as illustrated in Fig. 3. The BFS algorithm visits this tree-like 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

1:initialize *
4:     for  do
5:         extract the indices of :
6:         for all  satisfying  do *
7:              if  and  then *
8:                  Obtain a composite index at order:
9:                  *
10:                  if  then
11:                       Add to
14:                   *                             
15:     for  do
16:         if  then Delete *               
17:     Calculate the spectral contribution from
19:until the spectrum converges
Algorithm 1 Breadth-First Search for Pathways

Below are further instructions on the lines marked by asterisks.

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

  2. is a threshold for small matrix elements. One can set , where and is a user-defined relative threshold.

  3. 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.

  4. Compare to , contains one more e-h 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 .

  5. 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.

  6. 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 user-defined 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 matrix-element 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 open-source PYTHON simulation package.

Figure 3: The search tree in the BFS algorithm for finding all nontrivial minors of . The digits in the bracket denote the configuration, e.g., means . The semi-transparent configurations are discarded in the search process so that they don’t spawn any child configuration.

Here we demonstrate the BFS algorithm with a toy model with orbitals and valence electrons. Suppose of the system is


The BFS algorithm for this example of is carried out as follows:

First, the non-zero 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 non-zero 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 single-digit 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 e-h pair formed by orbital and ; (2) the core electron is first promoted to orbital and then coupled with the e-h 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 higher-order 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 brute-forcely 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 first-principles calculations into the determinant formalism

In the above sections, we have demonstrated an efficient solution to the MND model using many-body wave functions for simulating x-ray transition amplitudes. However, in order to simulate reliable x-ray 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 many-electron excited states (Eq. (13)).

For the final state, we employ the standard SCF core-hole approach to obtain the KS orbitals and eigenenergies. The core-excited 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 electron-core-hole pair, the core-excited 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 trade-off, the electron is only placed onto the lowest unoccupied orbital (), which we have dubbed the excited-state core-hole (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 core-hole (FCH) approach, in which the ground-state 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 core-excited atom is replaced by a ground-state 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 core-hole calculation, we can compute the orbital overlap integral for computing the determinantal amplitudes. To reduce computational cost, we employ projector-augmented-wave (PAW) form of ultrasoft pseudopotentials kresse1999ultrasoft ; blochl1994projector ; taillefumier2002x to model electron-ion interactions. The excited atom potential has deeper energy levels and more contracted orbitals so its PAW construction differ from the ground-state atom. In Appendix A, we derive a formalism to calculate expectation values between of the ground-state system and of the core-excited system. Initial dipole matrix elements are also evaluated within this PAW formalism blochl1994projector ; taillefumier2002x .

As in typical DFT impurity calculations, some low-lying excited states in the core-hole approach could be bound to the core-excited atom, resembling mid-gap localized electronic states near an impurity. In this situation, the electronic structure is well described by using a single k-point (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 k-point and take the k-point-weighted 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 plane-wave 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 k-points, the overlap -matrix for every k-point can be computed quickly as in Appendix B.

After the XAS is calculated by the first-principles determinantal approach, the established formation-energy calculation can be adopted to align spectra for core-excited 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 e-e interaction terms in the MND theory, which results in a single-determinant solution to the many-body wave functions, we argue that this first-principles determinantal approach does not entirely neglect valence e-e interactions. The self-consistent-field (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 x-ray 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 Fermi-liquid picture. Finally, corrective DFT (DFT + or DFT with exact-exchange functionals) or the self-consistent approximation bruneval2006effect ; van2006quasiparticle ; kang2010enhanced ; sun2017x can also improve QP energies and wave functions to be used in the determinantal approach. We imagine the many-body 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 core-hole potential described with the chosen exchange-correlation functional.

ii.7 Comparison with the one-body SCF core-hole approach

Before the determinantal formalism, the many-body transition amplitudes in the SCF core-hole approach are often approximated with 1p matrix elements


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 many-electron system (excluding the electron-core-hole 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 core-level transition and the many-electron response are not entangled. This is also the so called sudden or frozen approximation.

We know that from the diagrammatic interpretation of the x-ray many-body processes in Fig. 1, the photon first decays into an initial-state e-h pair instantaneously, and then the other electrons see the core-hole 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 initial-state orbital, and the subsequent many-electron 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


where are the matrix elements on the last column of as in Eq. (12) and is the minor complementary to . Since in the one-body core-hole approach only the terms are summed, we limit our analysis here to the many-body terms and condense the configuration tuple into a single index: . Then the matrix elements can be written as


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 core-hole 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 open-shell system.

The first term in Eq. (27) is more directly relevant to our previous one-body approximation. Here, is the minor of , the transformation matrix without its column and row. It reflects the -electron many-body overlap between the initial and final state occupied orbitals and should reflect the extent to which the electron density is modified by the core-hole perturbation. Since does not depend on , we can relate it to the many-body prefactor that appears in the final-state rule of Eq. (25): . Using the completeness relation: , the first term in the expansion of Eq. (27) can be expressed as


If it happened that , then this expression would amount to the final state matrix element as defined in the one-body final-state 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 one-body approach when this is not the case.

It appears useful to focus on to reveal the role of hybridization in modulating near-edge spectral intensity. To quantify the contribution of the second term in Eq. (28), we introduce the projection spectrum


in which the single index sums over all empty final-state 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 transition-metal oxides

In this section, we discuss an important application of the determinantal approach to computing core-excited state transition amplitudes, that is, to predict the x-ray 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 one-body core-hole approach. It has been found for a number of TMOs, that the one-body approach systematically underestimates the intensity of near-edge 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 X-ray 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 density-of-states related to the shell more easily interpretable in terms of band theory or effective 1p states.

The angularly-averaged (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 near-edge part of the spectra, i.e., the spectral features below eV contain the most useful information for material characterization. For these TMOs, the near-edge 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 near-edge part, so that one can interpret the spectra on a first-principles basis. More specifically, we use the ratio of the intensity of the first (lowest-lying) peak to that of the second (unless otherwise specified) as a metric for the accuracy of different levels of approximation.




\sidesubfloat[(b)] \sidesubfloat[(c)]

Figure 4: (Color online) (a) XAS for the selected crystal structures obtained from experiments (black), one-body FCH approach (blue), and the many-electron determinant approach (red) introduced in this work. The XAS calculated with the configurations are shown by dashed orange curves. The energy axes for \ceNiO and \ceSiO_2 are relative. (b) Comparison of experimental peak intensity ratios compared with the ones predicted by the one-body (circles) and the many-electron (triangles) formalism. Each color represents the result for one system. The peak intensity ratio refers to the ratio of the lowest-energy maxima to the second of the spectrum, unless otherwise specified by the numbers in (a). The spectra are broadened to the best as compared with experimental broadenings. (c) Schematics showing how the one-body and the many-electron formalism treats x-ray excitations, using a metal--\ceO- molecular model in both the initial (i) and final (f) state. The one-body approach mainly relies on the single-particle (1p) matrix element and has skipped (red arrow) the dynamics of the many-electron charge relaxation, while the many-electron formalism considers the actual multiple-step (blue arrows) excitation process that involves all the electrons in the system.

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