Ab initio methods for L-edge x-ray absorption spectroscopy

X-ray absorption spectroscopy (XAS) has become a critical tool to study molecular complexes and nanoclusters due to its ability to provide insight into local molecular geometry and electronic structure through the excitation of core electrons. While XAS measurements have been made since the discovery of x rays at the turn of the 20th century, advances in synchrotron technology and free electron lasers provide greatly improved resolution, both energetically and temporally.^1–10 Much of the power of XAS comes from the fact that it is element specific, meaning that the absorption spectra for different elements are highly separated energetically and are naturally resolved. With this unique ability to selectively probe elements within a compound, XAS is a key method in the investigation of the local electronic structure. For example, XAS has been successful in the characterization and study of metal complexes, including but not limited to charge transfer pathways,¹¹ oxidation states,¹² and spin crossover events.¹³ Recently, XAS has also been used to analyze the effect of the local environment around dopants or other defect centers in nanocrystals.^10,14,15

Since much of the power of XAS relies on the ability to correlate spectral features to particular electronic transitions, the interpretation of XAS is strongly aided by accurate electronic structure calculations. Most of the theoretical work to treat XAS from an ab initio perspective has focused on the K-edge region of XAS (excitations from the 1s core level), which is much more straightforward to compute than the L- or M-edges, which correspond to excitations from the n = 2 and n = 3 principle quantum numbers, respectively. However, L-edge XAS has several characteristics that make it particularly desirable. Like K-edge XAS, L-edge XAS is also element specific, but L-edge spectra tend to be more intense with formally allowed transitions and have finer linewidths due to longer core-hole lifetimes.^16,17 The finer linewidth allows for a higher sensitivity of element specific characterization, and rich spectral features, such as peak-splittings arising from spin–orbit coupling, can be resolved. Despite these experimental advantages, calculations of L-edge XAS are more challenging due to the fact that relativistic corrections and spin–orbit couplings play a larger role in the prediction of L-edge spectra.

Formally, relativistic effects are important to properly describe x-ray absorption because core electrons move at a significant percentage of the speed of light, causing core orbitals to contract and lower their energies. For K-edge XAS, not accounting for relativity in ab initio calculations uniformly redshifts the spectra compared to the experiment, but does not significantly alter the relative positions or intensities of spectral features.^18–22 Thus, for the most part, relativistic effects are often ignored or treated at the scalar-relativistic level in the ab initio treatment of K-edge XAS, and calculations are done using commonly implemented electron structure methods such as TDDFT,^20,22–28 coupled-cluster theory,^22,29–31 algebraic diagrammatic construction (ADC),^32–34 or transition potential theory.³⁵ This is not the case for L-edge spectra, where the excited core electrons are starting from the n = 2 principal quantum number atomic orbitals. Effects such as spin–orbit coupling cannot be accounted for by a uniform shift of the spectrum, but rather must be explicitly considered in the calculation. One way to include spin–orbit coupling is through perturbation theory in a configuration interaction picture. Multiconfigurational methods, such as the restricted active space (RAS) approach, have been applied to model both XAS by including the core orbitals in the active space.^36,37 A similar scheme has also been used to model L-edge spectra by combining restricted open-shell configuration interaction singles (ROCIS) and density functional theory.^38,39 On the other hand, more recently, both four- and two-component relativistic TDDFT methods, which include spin–orbit coupling variationally, have been developed to model L-edge XAS.^40–42 

The focus of this review is on recent work in electronic structure method development for L-edge XAS and is organized in three main sections. First, in Sec. II we describe the main physical processes underlying L-edge XAS in more detail. Next, in Sec. III a variety of the available theoretical approaches are discussed, starting with a brief overview of relativistic effects and the incorporation of spin–orbit coupling in traditional ab initio electronic structure theory methods. This is then followed by an overview of the variety of single-reference and multi-reference approaches that have been used to computationally treat L-edge XAS. Then, in Sec. IV we discuss several different numerical approaches to compute a spectrum, including time-domain approaches and iterative diagonalization techniques for interior eigenvalue problems, as well as damped frequency-dependent response and model order reduction. We conclude by presenting an outlook on L-edge XAS calculations in Sec. V along with future theoretical directions that we think would be beneficial for the x-ray community.

Core electron excitations can involve multiple physical processes due to the high energy nature of the experiment. For one, x-ray absorption spectroscopy measures both core electrons that are excited into bound valence states or ejected into the continuum. Additionally, multi-electron excitation processes are common, where core electrons are excited into either the continuum or a bound state in concert with a valence level charge transfer excitation, known as “shake-up” or “shake-down” transitions. A qualitative schematic of these processes is presented in Fig. 1. While these types of transitions are not a comprehensive description of all possible x-ray absorption processes, they represent the most common types that will be relevant in this review. We refer the interested reader to a more general review of the wide variety of x-ray methods from Norman and Dreuw.⁴³ 

There are two main energetic regimes in x-ray absorption spectroscopy, widely known as X-ray Absorption Near Edge Structure (XANES) and Extended X-ray Absorption Fine Structure (EXAFS).¹⁶ XANES typically refers to core excitations, which are excited into bound state orbitals, and EXAFS typically refers to core electrons excited into the continuum, which scatter off of neighboring atoms. In this topical review, we will focus on the theory of XANES, which is pertinent in examining the electronic structure of small molecules and clusters.

As mentioned previously, L-edge XAS offers a number of unique advantages over the K-edge. Specifically, L_2,3-edge spectra are of interest due to the core orbital being part of the 2p manifold. First, L_2,3-edge spectra are typically more intense than K-edge, due to it being a dipole-allowed $2p→d$ electronic transition for metal complexes with partially filled d manifolds. Additionally, the lower-energy nature of L-edge has a higher resolution due to a longer core-hole lifetime.¹⁶ In L-edge XAS, the 2s and 2p orbitals are not only contracted by relativistic effects, but the 2p orbitals are split in energy by spin–orbit coupling into 2$p1/2$ and 2$p3/2$ sets, denoted as the L₂ and L₃ edges in XAS, respectively. This phenomenon results in a rich structure of peak splittings of which the relative peak positions and intensities can be analyzed. These properties have made L-edge XAS a particularly good experimental probe of metal complex electronic structure in the d orbital manifold and has been used to resolve metal-ligand covalency.⁴⁴ 

As an example, the experimental L-edge spectrum of [FeCl₆]³⁻ is shown in Fig. 2 from Wasinger et al.⁴⁴ There are two main features in an L-edge spectrum, the L₃ and L₂ edge, which for [FeCl₆]³⁻ are centered approximately at 708 and 721 eV, respectively. These two separate edges arise from the spin–orbit splitting of the 2p manifold into the 2$p3/2$ and 2$p1/2$ manifolds. The branching ratio, defined by integration of the absorption signals in the two separate edges, $IL3$/$IL2$⁠, respectively, is statistically 2:1 for L-edge spectra, since there are four $p3/2$ states and two $p1/2$ states. However, this statistical branching ratio is typically only observed for single atom systems, since the covalency between the metal center and ligands and the symmetry of electronic wavefunction in complexes modify the electronic transition energies and amplitudes.^45–47 Within each of the L₂ and L₃ regions, there are additional subfeatures that correspond to transitions to the different sets of orbitals in the d-manifold. Additionally, in the [FeCl₆]³⁻ complex, a satellite peak is present at around 714 eV that is distinct from either of the two main absorption edges. This peak corresponds to a predominantly two-electron “shake-up” transition where a core 2p orbital in the L₃ manifold is excited along with a ligand-to-metal charge transfer excitation.^37,44 As we have highlighted so far, extracting chemical information from L-edge x-ray absorption spectra is indeed a complicated task, and an overview of the theoretical methods used to understand these spectroscopic signals are presented in Sec. III.

While there are a multitude of theoretical treatments for ab initio prediction of XAS, including muffin-tin potentials and real-space Green's functions,^48–50 we will focus only on electronic structure theory models using Gaussian basis functions with the linear combination of atomic orbitals (LCAO) approach. This is one of the most accurate and common choices for analyzing XANES in molecules and small finite clusters from first principles.

Although relativistic corrections are either ignored or treated at the scalar-relativistic level for K-edge features, since they largely amount to a uniform shift in the energy of the predicted spectrum and have little effect on the relative peak positions or intensities, this is not the case for excitations from the L-edge and M-edge core orbitals. In non-relativistic quantum chemistry methods, spin–orbit coupling is not included variationally in the Hamiltonian. As a result, the computed XANES spectrum for anything past K-edge will often be qualitatively incorrect for excited state calculations, although surprisingly good agreement can sometimes be achieved as seen in L-edge calculations of TiCl₄ by Stener et al.²³ 

In the weak coupling limit, perturbative spin–orbit coupling can be introduced into the non-relativistic Hamiltonian. For example, the Breit–Pauli spin–orbit operator,^51–56 

corresponding to the single- and two-electron spin–orbit effects, can be used to couple spin eigenstates obtained from a non-relativistic calculation. In Eq. (1), $si$ is the spin vector for each electron, $pi$ is the momentum operator, and index I refers to an atomic center. In practical calculations, the two-electron operator is usually neglected,^55,56 parameterized into the one-electron term,⁵⁷ or approximated using an atomic mean-field approach.^58–60 The Breit–Pauli spin–orbit operator can only be used perturbatively to couple spin-pure states because it is unbounded from below.^51,52 Diagonalization of the spin–orbit-perturbed Hamiltonian, through the interaction matrix elements,^{36–39,61,62} results in a new set of states that include spin–orbit contributions. In principle, if the full configuration interaction scheme is used, the resulting states are exact. The main advantage of the perturbative treatment of the spin–orbit coupling is that it can be carried out in the non-relativistic one-component computational framework. In this case, it is still generally necessary to variationally include scalar relativistic corrections to account for the effect of core orbital contraction.

While the spin–orbit interaction can be included perturbatively, it can also be included from first principles through explicitly relativistic approaches. Theoretically, spin–orbit coupling terms naturally arise from the relativistic Dirac equation. As a result, relativistic Hamiltonians are a natural choice for the ab initio treatment of L-edge spectroscopy. The four-component restrict-kinetic-balanced modified Dirac–Hartree–Fock equation for a molecular system is given in matrix form as^{51,52,63–68}

where c is the speed of light; ϵ is the orbital energy; and $T, V, S$ are the two-component matrices of the non-relativistic integrals for kinetic energy, potential energy, and the orbital overlap, respectively. For example,

where T is the N × N kinetic energy integral for N spatial basis functions. $C$ is the molecular orbital coefficient matrix, which has large and small components denoted by subscript L and S, respectively, and superscripts $+,−$ represent the positive energy electronic blocks and the negative energy positronic blocks, respectively. The W matrix can be expanded as

where $W0=p·Vp$ and $Wx,y,z=p×Vp$ give rise to scalar relativistic effects and spin–orbit couplings through the action of momentum operator p on the molecular potential. Inherently, this relation solves for both electronic and positronic solutions using the restricted-kinetic-balance condition,^63–67 where the definition of positronic and electronic basis functions are intrinsically coupled. Further details of four-component Dirac theory can be found in Refs. 51 and 52.

While the four-component Dirac equation can be used directly, there exist multiple different decoupling schemes that can separate the electronic and positronic blocks so that one can work with a smaller two-component wavefunction for electrons only. In the most general sense, decoupling these blocks corresponds to a unitary transformation of the Dirac Hamiltonian $H$ into a block diagonal framework. Mathematically, one seeks to find $U$ such that

The two-component block corresponding to electronic solutions of the four-component Hamiltonian can be calculated as

once the decoupling unitary matrix $U$ is known. A variety of approaches to determine $U$ and perform this transformation have been explored, such as the exact two-component method (X2C),^{41,42,69–76} Douglas-Kroll-Hess (DKH),^77–80 the zeroth-order regular approximation (ZORA),^68,81–83 and the Barysz-Sadlej-Snijders (BSS) method.^84,85 Most commonly the potential is assumed to only include the nuclear-electron attraction, with multi-electron effects only approximately included, using either a scaling factor⁵⁷ or atomic mean-field spin–orbit integrals (AMFI).⁸⁶ After this transformation, the electronic part of the Hamiltonian can be solved in a two-component spinor framework and has the same mathematical form as generalized Hartree–Fock (GHF), though we emphasize that from a theory perspective wavefunctions such as X2C-HF are not the same as GHF. As a result, this transformation from the four-component to two-component framework must also be taken into account when calculating properties to avoid “picture-change” error.

In nonrelativistic restricted Hartree–Fock (RHF) and unrestricted Hartree–Fock (UHF), the wavefunctions are one-component, in that each orbital is either spin-free as in RHF or the spin manifolds (α and β) are uncoupled (UHF) where they both preserve S_zspin symmetry. In two- and four-component Dirac–Hartree–Fock, spin symmetry constraints are lifted, where spinor orbitals are linear combinations of the α and β spin manifolds. Further, variational relativistic methods also require the relaxation of using real-value arithmetic, instead allowing the orbitals to have complex coefficients and the resulting wavefunction to take on complex values. A graphical representation of the relations between these methods are shown in Fig. 3.

While the perturbative approach within the one-component framework has the advantage of low computational cost, two- or four-component methods can variationally include the spin–orbit coupling at the molecular orbital level without assuming the weak perturbation limit. In addition, a variational reference can be directly incorporated into existing excited state electronic structure methods to provide a computationally simple workflow.

Single reference ab initio approaches are well established for simulating excited electronic states in molecular systems. The ground state of many molecular systems can be well described by a single Slater determinant of molecular orbitals computed with either Hartree–Fock or Kohn–Sham density functional theory (DFT). The excited state transitions are then solved for by calculating the response of the system to a perturbing electromagnetic field either in the frequency or the time domain.

Linear response time-dependent density functional theory (LR-TDDFT) is a well-established and computationally cheap single reference method for calculating the absorption spectra of molecular systems.^87–90 LR-TDDFT is commonly used to solve for lower-energy eigenstates that correspond to UV-Vis spectroscopic signatures. The computation of high-energy core to valence x-ray transitions in the XANES region can also be computed using the same theoretical foundation in an energy specific regime.^{20,21,23–26} The working equations for both UV-Vis and XAS for obtaining singly excited eigenstates can be found by computing the non-Hermitian eigenvalue problem

where the indices i, j and a, b are occupied and virtual orbital indices, respectively; ε is the orbital energy; $(ia|jb)$ are the two-electron repulsion integrals in Mulliken notation; $(ia|fxc|bj)$ is the DFT exchange-correlation contribution; and α is the Hartree–Fock exchange parameter when $0≤α<1$⁠. The equation reduces to time-dependent Hartree–Fock (TDHF) when α = 1. More details on the derivation of TDDFT and TDHF can be found in Refs. 90 and 91. In the approximation where $B=0$⁠, the problem is reduced to configuration interaction singles (CIS) (if α = 1) or the Tamm–Dancoff Approximation (TDA)^92,93 for DFT. Typically, these approximations are not as accurate, since the B matrix captures some ground state correlation effects when the reference state is approximate,⁹¹ which is true for single reference TDHF/TDDFT.

The dimension of the matrix depends on the number of occupied and virtual spin-orbitals, $2(Nocc×Nvir)$⁠. While all singly excited electronic transitions, from UV-Vis to XAS, are contained in the full matrix decomposition of Eq. (7), practical calculations of medium to large molecular systems require iterative eigensolvers for selected subspaces, since this matrix is typically too large to perform a full diagonalization. Therefore, in order to solve for selected XAS transitions, interior eigensolvers must be used to efficiently extract these transitions. A discussion of the algorithms used to solve this subspace problem for XAS transitions will be discussed further in Sec. IV.

Using TDDFT to solve for K-edge XAS is straightforward, since a non-relativistic reference can be used for computing the transitions. This is due to the fact that all of the electronic transitions start from core 1s orbitals that have zero orbital angular momentum. For L-edge, this is not the case, so spin–orbit coupling must be properly treated, as described previously in Sec. III A. The consequence of broken S_zsymmetry when using relativistic methods means that spins are free to rotate between spin-up and spin-down in three-dimensional space. This relaxation produces the added complication that spins can now be non-collinear, which requires the underlying electronic structure to adapt to this non-collinear framework and is typically not accounted for in traditional DFT method development. Recent work has expanded existing DFT functionals into the non-collinear regime, including two-component methods where variational relativistic effects are incorporated in the DFT reference determinants.^94–96 While this approach has had success, non-collinear relativistic TDDFT is still an active research field.

Instead of using non-collinear DFT, a variant of singly excited DFT methods, known as DFT/ROCIS has also been successful in predicting L-edge spectral features.^39,97 DFT/ROCIS uses Kohn–Sham DFT orbitals, but then performs configuration interaction singles on this reference with perturbative spin–orbit coupling. Unlike a two-component relativistic approach, this allows a standard one-component implementation to be used, which reduces computational cost. However, with a mixed configuration interaction and DFT method, some additional corrective parameters are required to obtain good agreement and prevent overcounting correlation. The uniformly shifted DFT/ROCIS spectrum is shown in comparison to experiment in Fig. 4. While qualitative features such as the relative intensities of the L₃ and L₂ edges are captured well, as expected, DFT/ROCIS (and relativistic TDDFT) both cannot capture the shakeup peak at 714 eV that has multi-electron excitation character. However, these methods are useful as a computationally efficient way to examine single electron transitions in L-edge XAS.

Using coupled-cluster theory, excited state energies and properties can be computed from several formalisms using a single reference state. One established method is known as equation-of-motion (EOM) coupled cluster. A more detailed review of this method can be found in Refs. 98 and 99, but we present a small overview of the method below.

In single reference coupled-cluster theory, an exponential cluster operator expands the Hartree–Fock ground state, $|Φ0⟩$⁠,

In Eq. (10), the coupled-cluster wavefunction is $|ΨCC⟩$⁠, and the generalized cluster operator is defined as $T̂=∑ktkτ̂k$⁠, where t_k is the cluster amplitude and $τ̂k$ is an N-body excitation operator. For example, truncating at double excitations known as CCSD is common, where the cluster operator is explicitly $T̂CCSD=∑iatiaâa†âi+14∑ijabtijabâa†âb†âjâi$⁠.

The EOM formalism^98,100 for solving excited states, involves the linear expansion of the coupled-cluster ground state represented by the operator $R$⁠:

From this relation, a non-Hermitian eigenvalue problem arises where we solve for the eigenvectors that represent excited state wavefunctions. In the EOM-CC method, all types of electronic excitations are treated on the same footing from the ground state reference. In order to model core excitations, one must solve for the interior high energy subsection of the whole eigenvalue problem, analogous to LR-TDDFT.

A multitude of different variants of excited state coupled-cluster have been developed to efficiently tackle the specific problem of solving for core electronic excitations. A straightforward approach comes from the energy-specific EOM-CC method, where high energy roots are solved in frequency space using a modified iterative solver that does not add any additional approximations.²² Alternatively, Nascimento and DePrince have applied a time-domain dipole autocorrelation variant of EOM-CCSD to calculate the K-edge XANES for several small molecules, including carbon, nitrogen, and oxygen.³¹ Extending this method to include a relativistic reference using the X2C framework has also been developed by Koulias et al.¹⁰¹ 

Another proposed method is using the core-valence separation (CVS) approximation to EOM-CC, in which an intermediate step projects out all eigenvectors that do not contain amplitudes from the core orbitals of interest during an iterative eigendecomposition.¹⁰² This provides a practical speed up for calculating XANES spectra since states that are not of interest are eliminated from the calculation. Extensions to this method, including a frozen core approach with an analytic alternative to the projection step, have also been explored.¹⁰³ Other related methods include the algebraic diagrammatic construction scheme (ADC). Typically this is truncated to second order, which is known as ADC(2). By applying the core-valence separation (CVS) approximation, Wenzel and co-workers have calculated K-edge spectra.^32,33

Additionally, recent work has explored computing L-edge spectra using CVS-EOM-CC and a perturbative Breit–Pauli spin–orbit Hamiltonian for a select group of small molecules.¹⁰⁴ In Fig. 5 the L-edge spectrum for SiH₄ using this method is presented alongside the experimental¹⁰⁵ spectrum. The L₃ and L₂ curves (102–104.5 eV) are captured qualitatively, although the broad shoulder peak near 102.5 eV is missing. Additionally, the Rydberg states from 105 eV onward are not easily compared, but this is expected due to the localized nature of using Gaussian type orbitals. This same system was calculated using the complex-polarization propagator method with relativistic TDDFT,¹⁰⁶ producing a spectrum that has similar qualitative results to the CVS-EOM-CC method. The main reason for the discrepancy with experiment is thought to be from vibrational broadening and cannot be remedied by introducing better theoretical descriptions of electronic correlation.

While dynamic electronic correlation and multi-configurational excited states are important for describing L-edge spectra, the problem of static correlation, where the ground state is highly multi-reference can also play a large role in the accuracy of computing L-edge spectra. The computation of L-edge spectra frequently involves open-shell systems and systems involving degenerate ground states. In such cases, in addition to accounting for relativistic effects, computation of the L-edge spectrum requires a description of static correlation in order to provide even a qualitatively correct electronic transitions. This may be achieved via use of multi-configurational wavefunctions, which also offer a convenient way to incorporate spin–orbit coupling perturbatively through a state interaction scheme.

One of the most common ways to build a multi-configurational reference ground state is through the complete active space (CAS) method. In CAS methods, the wavefunction $|Ψ⟩$ is written as a full interaction expansion of configurations (Slater determinants or configuration state functions) formed from a subset of the orthonormal molecular orbitals (MOs), known as the active space:

where $Ndet$ is the total number of configurations in the expansion. Minimization of the energy with respect to only the configuration interaction (CI) coefficients is referred to as CASCI. However, since the initial orbitals used to define the active space may not be optimal for the CAS wavefunction and for the state of interest, it is most common to re-optimize the orbitals self-consistently, such that the energy is stationary with respect to both the CI coefficients as well as the MO coefficients. This is the CASSCF procedure. Within CASSCF, in order to describe the core → valence excitations seen in L-edge spectra, both core and valence orbitals must be included in the active space. However, because CAS employs full configuration interaction, the scaling is poor as the size of the active space is increased, rapidly becoming intractable for the computation of L-edge spectra.

Generalizations such as the restricted active space (RAS) reduce the size of the space where full CI is used by partitioning the active space into different subspaces and imposing constraints on the number of holes and electrons in those subspaces in order to limit the number of configurations in the CI expansion.^107,108 The use of the restricted active space formalism naturally lends itself to the computation of x-ray spectra, as a limited number of excitations can be made from specified core orbitals in one subspace (referred to as RAS1) to the valence orbitals in RAS2, in which a full-CI is performed. This allows the simultaneous description of static correlation and the excitation processes inherent in x-ray spectroscopy. By selection of the appropriate orbital spaces, this approach has been used in the prediction of various x-ray spectra of transition metals, including K-edge, L-edge, and resonance inelastic scattering.^37,109–114

For the computation of L-edge XANES, spin–orbit coupling can either be included variationally during the orbital optimization procedure via a two-component or four-component formalism, allowing the orbitals to vary due to the presence of spin–orbit coupling, or it can be included as an a posteriori perturbative correction. To date, the perturbative RAS state interaction (RASSI) method^115,116 has been the most widely applied approach to include spin–orbit coupling for the study of L-edge spectra using RAS methods. In the state interaction method, RASSCF states of different spin multiplicity are converged before being allowed to interact via a spin–orbit coupling Hamiltonian, i.e., spin–orbit coupling is treated perturbatively after RASSCF convergence. At this stage, it is also common to include a RASPT2 correction to the energy of the states to account for dynamic correlation. This RAS approach has been used in studies of x-ray spectra of various (lighter) transition metal complexes, such as $FeCN64−, FeCl63−, FeCl4−$⁠, and extended iron porphyrin complexes.¹¹² Indicative of the applicability of this approach, it has been successful in reproducing experimental spectra for both high-spin (e.g., $FeCl63−$⁠) and low-spin (e.g., $FeCN6$⁠) ground states.

The form and performance of the approach can best be summarized using an example. By looking at the MO diagram for $FeCl63−$ (see Fig. 6), the simplest choice of RAS subspaces can be clearly identified. The L-edge XANES includes excitation of $t1u$ core (2p-like) orbitals and varied occupation of the valence $t2g$ and e_g orbitals. The $t1u$ core orbitals can be included in RAS1 and the $t2g$ and e_g orbitals in RAS2 where a full CI expansion is performed. By allowing a single excitation from RAS1 to RAS2, configurations are included in the CI expansion that have a single excitation from the RAS1 orbitals to the RAS2 orbitals, corresponding to the L-edge XANES. As noted above in Sec. II, the $FeCl63−$ XANES includes a shakeup peak at approximately 714 eV, corresponding to a coupled excitation of a $t1u$ electron and a σ electron. This peak can be described by including the σ bonding orbitals in the RAS2 space. The studies of Lundberg and co-workers³⁷ show that additional dynamic correlation can be included via inclusion of the $t2g$(π) and $t2g′$ orbitals. The computed L-edge XANES is shown in Fig. 7 and compared to that computed using X2C-TDDFT and the DFT/ROCIS methods.

The computed spectrum highlights the approach's proper description of shakeup transitions that require multi-electron excitations, which are not captured through other approaches like TDDFT. As with other methods used for the computation of L-edge spectra, such as TDDFT, a large number of states must be computed. However, due to a large CI expansion and inclusion of multi-electron effects, the number of states in RAS calculations is further increased. Since the cost of a CASSCF or RASSCF calculation severely limits the size of the active spaces, one approach that has been used to allow for correlation with a larger number of orbitals is by treating them with a truncated CI on the CASSCF wavefunction. This type of multireference configuration interaction (MRCI) allows for a better treatment of both static and dynamic correlation. For example, Sassi and co-workers used a MR-CISD method with perturbative spin–orbit coupling through the state interaction method to capture the L_2,3-edge for several Fe clusters with different multiplicities and oxidation states.¹¹⁷ An analysis of the spectra as shown in Fig. 8 was able to determine the relative abundance of these individual Fe centers in a magnetite sample and agreed well with the ideal stoichiometric ratio. These results agreed well with previous calculations from Ikeno et al.¹¹⁸ using a simpler CASCI-type approach as well as earlier atomic multiplet calculations.⁴⁶ 

We note that other methods such as atomic multiplet or ligand field multiplet calculations^46,119–122 are in a similar spirit to CAS and RAS methods, but use experimentally determined parameters (such as spin–orbit coupling, orbital covalency, etc.). These calculations model the atomic electronic structure of the metal center and treat molecular environment effects as a perturbation. Nevertheless, these simple models have been used to build up qualitative understanding of the important features and the physical phenomena underlying them.^44,123 However, given the reliance on determining parameters, they are less a tool of ab initio prediction than they are a tool of experimental analysis.

In formal theory, solving for the core-level excitations present in x-ray spectroscopies is no different than the valence excitations that appear in the UV-visible region of the electromagnetic spectrum. However, in practice there are several challenges that need to be addressed in algorithmic implementation. First, in iterative diagonalization it is necessary to selectively converge high energy excited states. Second, while eigenvector-based approaches are desirable for their interpretability of the orbitals and states involved, they are subject to non-trivial convergence problems and often need to converge a large number of states to cover the energy region of interest. In particular, since the number of states in a typical XANES calculation can be very large and lead to a dense eigenspectrum, alternatives to diagonalization of an appropriate Hamiltonian or response matrix are also possible and may offer better performance, such as when there are many (nearly) equivalent atoms being excited. One approach is to use the real-time time-dependent electronic structure framework, which enables the entire XANES region to be simulated at once.^40,41 Other alternatives to explicit eigensolvers or time propagation have also been explored, including frequency dependent-response^{29,30,106,124–128} and model order reduction.^28,129,130

By far, the most common approach to solve the preceding methods is to use an eigenvalue problem approach. That is, we solve a generalized eigenvalue problem of the form

where H is the Hamiltonian or response function to be diagonalized for eigenvectors $|xj⟩$ and eigenvalues λ_j, and M is a metric for normalization. Since the dimension of the matrices to be diagonalized are usually much larger than would be able to fit in memory, most implementations employ iterative diagonalization to only solve for the states of interest. For this, it is only required to be able to calculate the action of the matrix H on some generic vector $|b⟩$⁠. That is, instead of needing to build and store all of H, it is only necessary to form the matrix-vector product $H|b⟩$⁠. The Davidson algorithm is a common variant of Krylov subspace methods used to solve large eigenvalue problems in electronic structure theory.¹³¹ The basic idea of the Davidson algorithm is to work in a significantly reduced dimensional search space and find the best approximate (right) eigenvectors $xj$ as a linear combination of basis vectors $|bi⟩$⁠,

where the coefficients c_ij are found by diagonalization of the subspace problem. The only modification that needs to be made is to solve for large eigenvalues in the interior of the spectrum rather than only the lowest several roots. This “energy specific” method^22,25 involves selecting an energy threshold and only selecting the approximate eigenvectors with an energy greater than the threshold. Unfortunately, this method is not particularly efficient for dense manifolds of states as seen in XANES of large systems, since the convergence improves rather slowly as the search space expands. Alternatively, some methods propose eliminating unwanted excitations altogether by a restricted excitation window,^{20,23,24,26,132} growing excitation window,^21,25,26 or core-valence separation.^102,103,112

The use of restricted excitation window or core-valence separation methods has the advantage that it can generally alleviate the slow convergence issues in iterative diagonalizations, though depending on the implementation, its use may preclude the ability to capture processes such as shakeup transitions that include both core → virtual excitations as well as higher-lying valence → virtual excitations. Additionally, these approximations do not necessarily yield the same eigenvalue, eigenvector pairs that would be obtained from directly diagonalizing the full original matrix, though in many cases the difference is small. Several recent alternative approaches to specifically treat the interior eigenvalue problem are the generalized preconditioned locally harmonic residual (GPLHR)^133–135 method and iterative vector interaction (iVI).^136,137 These eigensolvers have been shown to be much more robust in locating the eigenvectors in the high-energy x-ray region, but without making any formal approximations. Despite the improved convergence properties, these methods still require significant computational cost to calculate the hundreds or even thousands of states in the desired region.

In the time-domain approach^138,139 (sometimes referred to as “real-time” to distinguish it from the response-based variants of methods such as linear response TDDFT), the time-dependent Schrödinger or Dirac equation is propagated explicitly. That is, one numerically solves

where $Ĥ$ is the Hamiltonian. This is often done using a symplectic integration scheme such as modified-midpoint unitary transformation (MMUT)^140–142 or 2nd-order Magnus^26,143 to ensure the conservation of energy and the number of electrons. Since the initial wavefunction is a stationary state, in order to get a spectrum, the system must first be perturbed. This is most commonly done by exciting all dipole-allowed electronic transitions with a δ-pulse along some direction q. The time-dependent Hamiltonian within the electric-dipole approximation is then given by

where $H0$ is the field free Hamiltonian, $κ(t)$ is the field strength, and $⟨rq⟩$ denotes the dipole integrals along direction q. Practically, the δ-pulse perturbation corresponds to a step function lasting for only the initial time step of width $Δt$⁠,

At each point in time, the expectation value of the dipole moment is then evaluated, $μ(t)$⁠. Following a Fourier transform to obtain $μ(ω)$ the (isotropic) dipole strength function is given by

One advantage of the real-time approach is the ability to efficiently calculate larger regions of the absorption spectra in a single simulation.^138,139 This is especially useful for the XANES region, which usually has a high density of states in the spectrum over tens of eV. Thus instead of computing hundreds or even thousands of eigenvectors to cover the spectral region of interest, a short time simulation provides the same spectra. Recent work has also shown that Padé approximants can also be used and produce comparable spectra with far shorter simulation time.^41,144 This is especially advantageous for XAS as the high energy requires a short time step. Kadek et al. used real-time propagation of the four-component Dirac–Kohn–Sham matrix to capture the L_2,3-edge region for SF₆, including a dipole-weighted transition analysis of the relevant spinors.⁴⁰ Recently the real-time X2C-TDDFT method was used by Kasper et al. to model the L_2,3-edge region for several molecular complexes. As seen in Fig. 9, the results are generally quite good, although there is some sensitivity to the choice of DFT functional and basis set.

Unlike the frequency or time-based approaches discussed previously, the damped response problem is inherently energy specific. In this method, one directly evaluates the response of the system (the output) to a given perturbation at frequency ω (the input).^145,146 The absorption cross section $σ(ω)$ for light at a given frequency ω is proportional to the trace of the dynamic polarizability $α(ω)$⁠,

where

and $η>0$ is a small damping parameter. The use of the complex $ω̃$ allows convergence near the poles where there are resonant excitations that formally diverge. The evaluation of the tensor $α(ω̃)$ can be performed by linear solvers according to

where H is the given response Hamiltonian, S is the overlap, and d is the vector of dipole operators, all in the molecular orbital basis. This is because $Tr[dTx(ω̃)]=α(ω̃)$⁠. A variety of Hamiltonians have been used to obtain the complex polarization propagator (CPP), including CPP-CC^30,124,125 and CPP-SCF,^146–149 and can be used for high-energy excitations.^{29,106,126–128,150} Now Eq. (22) can be evaluated for a given set of frequencies ω_i to obtain the spectrum over the range of values. This results in spectra that are similar to those obtained through the real-time approach, as there is not explicit calculation of the poles that corresponds to discrete transitions. For example, Frannson et al.¹⁰⁶ used CCP-DFT to compute the L_2,3-edge at both the four-component level of theory as well as X2C, as shown in Fig. 10.

In general, this procedure would require a large number of evaluations of Eq. (22) to cover a large energy range with sufficient resolution, but with techniques such as model-order reduction (MOR) the number of required evaluations can be made much smaller.²⁸ This essentially provides an interpolation scheme to evaluate a minimal number of linear systems to converge the spectrum to within some tolerance, since more frequencies are required to be evaluated near resonances with much fewer far away from resonance. Since all the linear system solves are independent problems for a given frequency ω, this approach is particularly good for its ease of parallelization. For larger systems where there are many nearly degenerate transitions such as the water clusters considered in Ref. 28 where the spectrum will be more dense and convergence of iterative diagonalization is harder, an approach using damped response and MOR is expected to perform well.

From the computational and theoretical standpoint, L-edge XAS presents some unique challenges over describing excitations in the UV-visible region, mainly due to the need to include relativistic corrections (scalar relativity and spin–orbit coupling), to compute the interior eigenspectrum, and to handle the dense manifold of excited states.

Throughout this topical review, a variety of different electronic structure methods and practical algorithmic strategies for computing molecular L-edge x-ray absorption spectra have been summarized. While new theory is still being developed, there are now many options for simulating L-edge XAS that each have their own strengths and weaknesses. For one, TDDFT and ROCIS/DFT are computationally cheap and provide good qualitative results for most single electronic excitations. However, TDDFT and ROCIS/DFT cannot capture shakeup excitations and other multi-electron effects that appear in x-ray spectroscopy. In order to capture multi-electron excitation effects, a higher level of theory such as coupled cluster, ADC, or CAS/RAS must be used. These are much more expensive than single excitation methods and currently are only practical for small complexes, though this may change with future breakthroughs.

Another important choice when calculating L-edge spectra is choosing between solving the problem in the frequency domain or the time domain. For example, TDDFT is typically solved as an eigenvalue problem in the frequency domain where eigenvalues are excitation energies and the eigenvectors represent the approximate excited state. Using these vectors, their energies, and their corresponding absorption strength, one can build up a spectrum with simple Gaussian or Lorentzian broadening. However, for systems with multiple atoms of the same type the number of states necessary quickly grows into the thousands, leading to high cost and convergence problems. Alternatively, time-domain methods that explicitly propagate the electronic density matrix in time can obtain a full spectrum in a single calculation instead of solving for many eigenvectors. This method can be advantageous if one desires only the shape of the spectrum, rather than a careful analysis of specific transitions. In contrast to either of these methods, the problem can also be recast as a damped response using the complex polarization propagator. This allows for explicit calculation of the damped response for a set of frequencies. Additionally, newer techniques such as MOR can reproduce similar results with great efficiency and are well-suited to parallelization across large supercomputers. For some systems, it might be advantageous to use a non-eigenvector-based method to first obtain a spectrum and then tailor an eigenvector-based approach to examine the states in a given feature of interest.

J.M.K., T.F.S., and A.J.J. contributed equally to this work.

X.L. acknowledges support from the U.S. Department of Energy, Office of Science, Basic Energy Sciences, in the Computational and Theoretical Chemistry program under Award No. DE-SC0006863 for the development of the first-principles electronic dynamics and relativistic TDDFT. The development of the relativistic multi-reference method is supported by the U.S. Department of Energy in the Heavy-Element Chemistry program (Grant No. DE-SC0021100). The development of the linear response TDDFT method for computational spectroscopy was supported by the National Science Foundation (Grant No. CHE-1856210). The development of the open source software package is supported by the National Science Foundation (Grant No. OAC-1663636). Computational L-edge spectroscopy applications are partially supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences and Biosciences, through Argonne National Laboratory under Contract No. DE-AC02–06CH11357. The relativistic coupled-cluster work was supported by the Computational Chemical Sciences (CCS) Program of the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Chemical Sciences, Geosciences and Biosciences Division in the Center for Scalable and Predictive methods for Excitations and Correlated phenomena (SPEC) at the Pacific Northwest National Laboratory. The model order reduction method was supported by IDREAM (Interfacial Dynamics in Radioactive Environments and Materials), an Energy Frontier Research Center funded by the U.S. Department of Energy (DOE), Office of Science, Basic Energy Sciences (BES).

The data that support the findings of this study are available from the corresponding author upon reasonable request.