Modeling L_2,3-edge X-ray absorption spectroscopy with linear response exact two-component relativistic time-dependent density functional theory

X-ray absorption spectroscopy (XAS) is an important tool that can give insight into local molecular geometry and electronic structure through the excitation of core electrons in molecular compounds. Advances in synchrotron technology with greatly improved temporal and spectroscopic resolution have made XAS a powerful tool for investigating the electronic and nuclear structure of molecules and condensed matter.^1–10 X-ray absorption is element specific, meaning that the absorption spectra for different elements are highly separated energetically. Because of this energetic separation, XAS has the unique ability to isolate the selected elements for physical study. For example, XAS has been especially successful in the characterization and study of metal complexes, including but not limited to charge transfer pathways,¹¹ oxidation states,¹² observation of spin crossover events,¹³ and understanding solvation effects.^14,15

The near edge X-ray absorption fine structure (NEXAFS), also known as X-ray absorption near edge structure (XANES), contains excitations to bound electronic states close to the ionization potential.¹⁶ In the extended X-ray absorption fine structure (EXAFS) region, continuum effects become more prominent due to ionization and Feshbach resonances.¹⁷ The NEXAFS region is difficult to capture without an accurate description of the electronic structure of the absorbing atom and its neighbors. Previous theoretical methods have been well studied in the nonrelativistic regime, but newer relativistic methods have not been characterized in great detail.^18–22 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 energy. For K-edge XAS (excitations from the 1s core level), neglecting relativistic corrections in ab initio calculations uniformly redshifts the spectra compared to experiments.^23–27 However, the relative peaks in the spectrum remain similar, so the overall characterization of K-edge spectra is commonly treated by using a nonrelativistic Hamiltonian and uniformly shifting the computed spectrum.

This is not the case for L-edge spectra, where the core orbitals that are being excited are in principal quantum number n = 2. Like K-edge XAS, L-edge XAS is also element specific, but L-edge spectra have finer linewidths due to longer core-hole lifetimes.¹⁶ The finer line width allows for a higher sensitivity of element specific characterization, making L-edge XAS a unique probe of molecular properties. Transition probabilities for L-edge spectra are electric dipole-allowed, so issues with origin-dependence need not be considered in computations.^24,28–31 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 2p_1/2 and 2p_3/2 sets, denoted as the L₂ and L₃ edge in XAS, respectively.

Theoretically, spin-orbit coupling terms fall out of the relativistic Dirac equation. As a result, relativistic Hamiltonians are a natural choice for the ab initio treatment of L-edge spectroscopy. We do note that it is possible to include spin-orbit coupling perturbatively in a nonrelativistic or scalar relativistic Hamiltonian, and methods of this type have been implemented and shown to provide a qualitative description for light elements, such as ZORA,^32–34 MCSCF,³⁵ RASPT2,³⁶ time-dependent density functional theory (TDDFT),³⁷ and restricted-open-shell configuration interaction with singles (ROCIS).^38,39 However, perturbative treatment of relativistic effects will break down for heavier elements.

Previously, variational treatments of spin-orbit coupling in the electronic Hamiltonian within the real-time time-dependent electronic structure framework have been used to model X-ray absorption spectra.^40,41 These approaches, without explicitly solving for the eigenvectors of excited states, have successfully produced XAS spectra in excellent agreement with experiments. Other alternatives to explicit eigensolvers have also been explored, including frequency dependent-response^19,42–48 and model order reduction.^49,50 When eigenvectors of electronic adiabats are needed for the interpretation of the chemical properties of XAS, the linear response formalism of the time-dependent Hartree-Fock (LR-TDHF) and time-dependent density functional theory (LR-TDDFT) with the restricted^25,51–54 or growing excitation window^26,54,55 techniques have been particularly successful for predicting K-edge XAS. However, for systems that exhibit a high degeneracy of excited states in the spectral region of interest, e.g., L-edge XAS of a metal complex, solving the relativistic LR-TDDFT equations is often subject to nontrivial convergence problems. To obtain the interior eigenspace of challenging spectroscopic systems, we have recently developed a hybrid method using the generalized preconditioned locally harmonic residual (GPLHR)⁵⁶ algorithm for solving the real-valued TDDFT equation.⁵⁷ 

In this work, we introduce the implementation of the GPLHR algorithm to solve the complex-valued relativistic two-component TDDFT. The aim of this work is to investigate and calibrate the ability of linear response formalism of exact two-component transformed (X2C) time-dependent Hartree-Fock (X2C-TDHF) and time-dependent density functional theory (X2C-TDDFT)^58,59 for modeling the L-edge X-ray absorption spectra with a focus on first and second row transition metal complexes.

Relativistic two-component methods begin with the four-component (4c) relativistic Dirac Hamiltonian for the electron,^60–62 

where V is the scalar potential, σ is the vector of Pauli matrices, p is the linear momentum operator, and I is the rank-two identity matrix. The speed of light and mass of the system is given by c and m, respectively. The eigenstates of this Hamiltonian can be partitioned into large (L) and small (S) components,

that can be subsequently partitioned into spin-up and spin-down components. Since there are both positive and negative energy solutions to the Dirac equation, a variety of algorithms have been devised that isolate the electronic (positive energy) solutions via a unitary transformation,

The Hamiltonians Ĥ⁺ and Ĥ⁻ have two-component (2c) eigenstates Ψ^2c and describe positive and negative energy states, respectively. Although there are many different techniques to form the transformation matrix $U$⁠, in this work, we use the exact two-component (X2C) method.^{58,59,63–75} In our implementation, we restrict the 4c → 2c transformation to the one-electron operator. Because of this approximation, an additional scaling factor for the spin-orbit coupling terms is included to account for the two electron terms in an approximate manner.⁷⁶ For Sec. II B, we assume that the 4c → 2c transformation has been carried out so that all matrix quantities are in the transformed two-component framework, as represented by the tilde.

The absorption spectrum for molecules can be obtained using TDDFT.^77,78 One of the most common approaches is to solve the TDDFT equations using the linear response formalism.^79–81 This amounts to solving an eigenvalue problem of the form

where ω is the resonant energy, and $X̃$ and $Ỹ$ are the transition amplitudes. The matrices $Ã$ and $B̃$ are given by

and

where a, b denote virtual MOs and i, j occupied MOs. The coupling matrix $K̃$ is given in the AO basis as

The functional second derivatives in Eq. (7) depend on the particular exchange and correlation functionals. In this work, we use a direct atomic-orbital based X2C transformation with the torque-free spin-density approach as described in Refs. 58, 59, and 74. We note that the extension of standard collinear functionals to 2-component noncollinear DFT is nontrivial, and refer readers to Refs. 74 and 82 for a recent review on noncollinear DFT methods. The solution of Eq. (4) can be performed in the X2C-transformed frame, since the transition energies E and vectors |XY⟩do not depend on the transformation matrix $U$⁠.⁵⁸ However, the evaluation of properties must be done in the nontransformed basis, as is required to avoid the so-called “picture change” error.^58,83

In general, the size of the matrix in Eq. (4) is too large to store in memory and so a variety of iterative methods have been used. Although there is extensive discussion in the literature on iterative eigensolvers, two specific aspects of this problem warrant further elaboration. First, for high energy excitations such as those in X-ray spectroscopy, some type of energy windowing is desirable in the iterative diagonalization of Eq. (4). Obtaining only a subset of the spectrum greatly reduces the computational cost by removing the need to converge all lower-energy roots. Energy-specific Davidson^27,55 and restricted excitation window^52–54,84 methods have proved to be fairly successful in this regard. However, recently some alternative methods such as GPLHR^56,57,85 and IVI^75,86 have shown great promise in tackling challenging interior eigenproblems in TDDFT.

Second, although no less important, is the need to use complex arithmetic. Typically, implementations of TDDFT in real arithmetic make use of the fact that A and B are Hermitian to simplify Eq. (4) to an equivalent problem of half the original dimension,⁷⁸ 

However, in the general complex case, A ≠ A^* and B ≠ B^* and Eq. (4) cannot be reduced to the form of Eq. (8). Despite this, the problem still possesses a block structure that can be taken advantage of to improve convergence. Most notably, while general non-Hermitian eigenproblems can have complex eigenvalues, solutions to Eq. (4) still have real eigenvalues. To eliminate intermediate complex eigenvalues, one can enforce the block structure of Eq. (4) on the subspace problem at each iteration. This is accomplished by incorporating matching trial vectors into the space. While in the real case for every positive energy solution (ω, |XY⟩), there is a negative energy solution (−ω, |YX⟩) for the complex non-Hermitian problem, for each (ω, |XY⟩), there is a corresponding “paired” solution (−ω^*, |Y^*X^*⟩). With both trial vectors included, we retain the same structure as Eq. (4) for the subspace problem at each iteration.

The full technical report and derivations of the GPLHR algorithm are already published,⁵⁶ so here we only provide an overview of the algorithm. The TDDFT equation in Eq. (4) is a generalized eigenvalue problem of the form

The GPLHR algorithm attempts to find the n interior eigenvectors with eigenvalues Ω closest to a specified value σ through a harmonic Rayleigh-Ritz procedure. This procedure is general for a complex non-Hermitian eigenvalue problem, which can in principle have complex eigenvalues, but in the case of LR-TDDFT, $H$ represents the matrix in Eq. (4). For physically meaningful solutions, $H$ will have all eigenvalues real; complex eigenvalues indicate an instability in the reference. However, intermediate complex eigenvalues can arise during the iterative diagonalization, and care must be taken to ensure these are handled properly.

From an orthonormalized initial guess of n (right) vectors $V$⁠, the set of vectors $Q$ is formed in the σ-shifted space,

After orthonormalization of $Q$⁠, the initial approximations for the eigenvalues $ω̃i$ are found by solving the generalized eigenproblem in the reduced space given by

This can be solved by generalized Schur decomposition (also known as QZ factorization) and is preferable over diagonalization when the matrices might become low-rank. From the Schur decomposition of Eq. (11), we obtain the triangular factors $R1$⁠, $R2$⁠, as well as Schur bases $YL$⁠, $YR$⁠, such that

The eigenvalues are given by the ratio of the diagonal elements of the triangular factors,

The Schur vectors $YL$ and $YR$ can be used to update $V$ and $Q$⁠, with $V→VYR$ and $Q→QYL$⁠.

Next, a subspace is generated. This trial subspace $Z=[V,W,S1,…,Sm,P]$ corresponds to the Krylov-Arnoldi sequence generated by the preconditioned residuals. Note that unlike the Davidson method which increases the subspace size at each iteration, the subspace is of fixed size determined by the integer parameter m. $P$ is an additional block of saved vectors and is not used on the first iteration. As is done in the Davidson algorithm, we approximate the preconditioner $T$ as the inverse of the difference between the approximate eigenvalue $ω̃$ and the diagonal elements of $H$⁠, which are given by the orbital energy differences. Additionally, GPLHR projects out components against $V$ and $Q$⁠. For the trial space $Z$⁠, the analogous set of vectors in the σ-shifted test space is formed,

The new reduced-dimensional eigenvalue problem $(U†HZ,U†MZ)$ is solved by Schur decomposition as before. Since the dimension of the eigenproblem in Eq. (15) is now larger than the number of roots we seek, the eigenvalue-eigenvector pairs are ordered by closeness to the shift value σ and can then be used to update $V$⁠, $Q$⁠, and $P$⁠. The convergence of eigenvectors and eigenvalues is evaluated to determine if additional iterations are necessary.

To better understand how well the DFT-based approach should be expected to perform generally, a variety of different density functional and basis set combinations are used to compute L_2,3 XAS for a set of transition metal complexes. This set includes metal centers in several rows of the periodic table to test the robustness of the method to capture different symmetry environments as well as varying strengths of bonding and relativistic effects. The basis sets used were the nonrelativistic Pople-type 6-311G(d),^87,88 the relativistically optimized double-ζ and triple-ζ Sapporo sets,⁸⁹ and the relativistically optimized aug-cc-pVTZ-dk.^90–93 The DFT functionals included B3LYP,^94,95 PBE0,^96,97 as well as HSE06,⁹⁸ CAM-B3LYP,⁹⁹ and M062X.¹⁰⁰ Results from Hartree-Fock (HF) are also given for comparison. All calculations were run using a locally modified version of the development version of the Gaussian suite of programs.¹⁰¹ For each complex, the geometry was optimized using the LANL2DZ^102–104 effective core potential and basis set with each density functional. The optimized geometries for the test set of complexes are shown in Fig. 1. Convergence of the linear response X2C-TDDFT problem was done using a combination of the energy-specific Davidson method^27,55 and GPLHR^56,57 as a hybrid method. It is also important to note that one of the key challenges in solving the linear response X2C-TDDFT problem for L-edge spectra is that there are a large number of roots required to obtain the spectrum. The degeneracies in magnetic microstates (2J + 1) lead to a large number excited states, even if they do not contribute oscillator strength to the overall spectrum.

The use of GPLHR has been shown to provide improved convergence in TDDFT problems, particularly where there is a dense manifold of states associated with high degeneracy.⁵⁷ For complex two-component TDDFT problems, the proper handling of complex intermediates is essential. We remark that large scale systems can introduce a high level of degeneracy of excited states which are spatially distinct in nature. For systems with nonzero total angular momentum, J > 0, an additional level of degeneracy arises from magnetic microstates. In this section, we only focus on the convergence behavior of a complex GPLHR interior eigensolver for obtaining excited states in the L-edge spectral region using complex X2C-TDDFT. We refer interested readers to Ref. 57 for additional discussion of the computational performance of this type of eigensolver.

In Fig. 2, we show the convergence profile and computational cost for computing excited states in the L-edge of SiCl₄ using the aug-ccpVTZ-dk basis and the CAM-B3LYP functional. The transitions obtained here are for the first 5 states above 103 eV and belong to the L₃-edge. Since the matrix-vector product (m-vec) is the time consuming step, the total number of m-vecs is used to measure the computational cost. As seen in the convergence profile, the energy-specific Davidson is less smooth and not monotonic. In the energy-specific Davidson algorithm, a fixed number of trial eigenvectors are selected for a given energy window. In cases where there is a high level of degeneracy, there can be a fluctuating selection of trial eigenvectors, resulting in an ill-conditioned convergence profile. While this hinders the Davidson algorithm, the problem appears less severe than high levels of spatial degeneracy, presumably because similar sets of orbitals are needed to describe the degeneracies due to different spin states. By contrast, in the GPLHR hybrid method, there is no hard threshold used to include or exclude certain eigenvalue-eigenvector pairs. As a result, the GPLHR hybrid exhibits a faster and smooth convergence.

For each molecular complex, several peaks were identified for comparison to experimental spectra,^34,105–108 resulting in a total of 40 transitions used in the analysis. However, since the number of states required to cover the X-ray spectrum is quite high and each experimental peak may be composed of multiple roots, instead of analyzing individual eigenvalue-eigenvector solutions, we obtain a Gaussian-broadened spectrum and identify peaks and shoulders using a spline fit and solving for the minima of the second derivative. The broadening parameters were chosen to match qualitatively with the experimental spectra and are given in the supplementary material. Note that we were unable to obtain the full spectra for certain combinations of DFT functionals and basis sets due to strong numerical instabilities in the noncollinear formalism. The errors in energy are computed using several methods: unshifted, a single uniform shift for both L₂ and L₃, and using an additional shift for L₂ relative to the shifted L₃. A single uniform shift is used to correct for self-interaction errors in DFT,¹⁰⁹ and the additional L₂-edge shift corrects for the approximate two-electron spin-orbit coupling.⁴¹ Figure 3 plots the shifted spectra of $[FeCl6]3−$ compared to experimental measurements.¹⁰⁵ 

The mean absolute error (MAE) and range of the excitation energies are given in Table I. The MAE is also shown in Fig. 4 graphically. Since the scale of excitations changes considerably, the data from PdCl₂ are not included in Table I. When only unshifted errors are considered, self-interaction errors dominate the MAE analysis among the DFT results. There is little difference in MAE between the use of the nonrelativistic and the relativistically optimized basis sets, nor is there a significant difference between different functionals. The large error for HF results is likely due to the lack of electron correlation.

After applying a uniform shift to align the L₃ edge with the experimental value, it is clear that basis sets optimized for relativistic calculations provide a better description of the spectrum. Unsurprisingly, using the triple-ζ Sapporo set performs better than the double-ζ. Most noticeably, the HSE06/Sapporo(TZ) and HSE06/aug-cc-pVTZ-dk levels of theory show a remarkably small MAE of 0.24 eV and 0.34 eV, respectively. These two levels of theory also show a small error range (1.53 eV and 1.56 eV, respectively). Table II shows the relative error in the computed L₂–L₃ splitting compared to experiments. All levels of theory, except for the HF/Sapporo(TZ) and HF/aug-cc-pVTZ-dk, underestimate the magnitude of spin-orbit coupling. It is clear that the description of spin-orbit coupling is greatly improved in aug-cc-pVTZ-dk and Sapporo(TZ) compared to 6-311G(d), most likely due to tight functions near the core. L₂–L₃ splittings from the range-separated HSE06 functional with relativistic-optimized basis sets are in the best agreement with experiments.

When an additional spin-orbit correction is included by separately aligning the L₂ edge with the experimental value relative to L₃, Table I shows that there is not a significant difference in quality between nonrelativistic and relativistic basis sets. This is because after all relativistic corrections are applied the overall shape of each edge is dominated by the description of the valence orbitals, which both relativistic and nonrelativistic basis sets are able to model with a similar quality. This error analysis also suggests that approximate two-electron spin-orbit treatment, such as the Boettger scaling factor used in this work,⁷⁶ should be sufficient for the calculations of L-edge spectra.

We continue the analysis by including heavier elements (PdCl₂) in the benchmark analysis, and no longer include the results with the 6-311G(d) basis. Since this test set contains molecules with metal centers on several periods, the L-edge energies vary from around 100 eV to over 3 keV. In order to appropriately quantify error for transition energies of different orders of magnitude, it is useful to examine the percent error,

In Table III, the percent errors across the entire test set are shown and also displayed graphically in Fig. 5. Mean percent errors (MPEs) and error ranges for unshifted spectra are significant, especially given the large L-edge excitation energy. After the uniform shift, both the mean percent error and error range are drastically improved, with CAM-B3LYP/aug-cc-pVTZ-dk and CAM-B3LYP/Sapporo(TZ) showing a mean percent error of only 0.08%. When an additional spin-orbit correction is applied, these errors are further improved. For the functionals with the PBE correlation (HSE06 and PBE0), the aug-cc-pVTZ-dk basis performs better than Sapporo(TZ). However, the reverse is true for B3LYP and CAM-B3LYP. There does not appear to be any significant difference in overall errors across the different functionals, although CAM-B3LYP had slightly better average error and one of the smallest ranges of error. Surprisingly, the range-corrected HSE06 does not perform better than PBE0 in this study, while CAM-B3LYP performed slightly better than B3LYP, which seems to suggest that range-corrected functionals do not consistently improve L-edge spectra.

Although obtaining correct excitation energies is important, the computed oscillator strengths must also be correct to match the experimental spectrum. Additionally, relative intensities such as the ratio between the L₂ and L₃ regions are known to contain information about oxidation state and covalency.^110–112 The branching ratio is defined as

where $IL2$ and $IL3$ are the integrated intensities of L₂ and L₃ edges, respectively. As plotted in Fig. 6, the computed branching ratios of $[FeCl6]3−$ reflect how accurately the method captures the peak intensities. The double-ζ basis consistently underestimates the experimental branching ratio, most likely due to less variational flexibility. While the triple-ζ basis sets perform more accurately on average, they are also subject to more variation with respect to the DFT functional.

We continue our discussion to assess the benefits and drawbacks of using TDDFT to model L-edge spectra. The L-edge X-ray absorption spectrum of $[FeCl6]3−$ has been previously studied in detail with experimental data and multiple different theoretical techniques, including the restricted active space second order perturbation (RASPT2) method,³⁶ restricted-open-shell configuration interaction with singles (ROCIS),³⁹ and the ligand field multiplet (LFM) method.¹⁰⁵ The experimental¹⁰⁵ and computed spectra using RASPT2,³⁶ ROCIS/DFT,³⁹ and X2C-TDDFT from this work for $[FeCl6]3−$ are compared in Fig. 7. Each spectrum is normalized and uniformly shifted to match the main experimental L₃ edge peak at 708.5 eV. All three methods are able to qualitatively reproduce the main features in the L-edge spectrum of $[FeCl6]3−$⁠. Both ROCIS/DFT and RASPT2 results show more peak features than that computed using X2C-TDDFT. However, the additional peak features from ROCIS/DFT seem to arise from overestimated intensities. The RASPT2 result more closely matches the experimental spectrum. Additionally, RASPT2 is the only method able to capture the peak at 713.5 eV, which is due to a “shake up” ligand-to-metal charge transfer excitation.³⁶ This is possible because RASPT2 includes electronic configurations with more than just single excitation character.

In order to investigate the role of multideterminantal effects in L-edge spectra further, we examine the electronic structure characteristics of VOCl₃. Figure 8 shows the experimental³⁴ and computed L-edge spectra using LFM¹¹³ and X2C-TDDFT for VOCl₃. Each spectrum is normalized and uniformly shifted to the central L₃ edge peak centered at 513 eV. The spectrum denoted as LFM contains the (2p⁵) (3d¹) electronic configuration on the metal center.¹¹³ Including additional doubly excited determinants with ligand-to-metal charge transfer character leads to a more accurate spectral prediction, denoted as LFM+CT.¹¹³ 

As shown in Fig. 8, the X2C-TDDFT result is similar to that computed using LFM+CT for the L₃-edge, with the LFM+CT spectrum agreeing better with the experiment. The LFM approach performs less satisfactorily compared to LFM+CT and X2C-TDDFT. This is likely due to the fact that the X2C-TDDFT result used a hybrid functional, B3LYP, while the spin-orbitals in LFM were computed using the local density approximation. The addition of ligand-to-metal charge transfer from the doubly excited determinants in LFM+CT seems to significantly improve the L-edge prediction. The L₂-edge spectrum computed from LFM+CT shows the best agreement with experiments. This is because the ligand-to-metal charge transfer shake-up states are more important for the electronic transitions in the L₂ region for VOCl₃.¹¹³ As a result, the X2C-TDDFT and LFM L₂ spectra are in less satisfactory agreement with experiments.

Generally, the relativistic X2C-TDDFT method exhibits excellent agreement with experiments for main L_2,3 peaks dominated by single electron transitions. However, within the linear response formalism, X2C-TDDFT is unable to resolve spectral features that need multideterminantal treatment.

In this paper, we have presented the LR-X2C-TDDFT approach using the hybrid GPLHR-Davidson diagonalization method for the computation of L_2,3-edge spectra. Several different density functionals and basis sets are used in the benchmarking and comparison. While X2C-TDDFT cannot model satellite “shake-up” peaks due to their doubly excited character, it accurately captures the single excitation features and the energetic splitting of the L₂ and L₃ peaks. The choice of basis set did not have a strong effect on the qualitative character of computed spectra, although a reasonable quantitative improvement is seen from using a triple-ζ rather than double-ζ basis. By contrast, the effect of the choice of exchange-correlation functional has a much stronger effect. Among standard GGAs and hybrid functionals such as B3LYP or PBE0, there was no significant difference with each performing well in these data. Using range-corrected functionals had only marginal impact on overall quality, with HSE06 slightly less accurate than PBE0, and CAM-B3LYP slightly better than B3LYP. From our test set of metal-centered compounds, the combination of CAM-B3LYP and a relativistically optimized triple-ζ basis set gives us the best result for predicting L-edge spectra.

See supplementary material for Cartesian coordinates of molecular systems, broadening parameters, computed branching ratios, and uniform shift values.

The development of the energy-specific TDDFT method was supported by the National Science Foundation (Grant No. CHE-1856210). The development of the relativistic two-component DFT/TDDFT method is funded by the U.S. Department of Energy (Grant No. DE-SC0006863). Computational L-edge spectroscopy was 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 the Argonne National Laboratory under Contract No. DE-AC02-06CH11357. Computations were facilitated through the use of advanced computational, storage, and networking infrastructure provided by the Hyak supercomputer system at the University of Washington, funded by the Student Technology Fee and the National Science Foundation (Grant No. MRI-1624430).