In this paper, we present a method to compute the x-ray absorption near-edge structure (XANES) spectra of solid-state transition metal oxides using real-time time-dependent density functional theory, including spin–orbit coupling effects. This was performed on bulk-mimicking anatase titania (TiO2) clusters, which allows for the use of hybrid functionals and atom-centered all electron basis sets. Furthermore, this method was employed to calculate the shifts in the XANES spectra of the Ti L-edge in the presence of applied electric fields to understand how external fields can modify the electronic structure, and how this can be probed using x-ray absorption spectroscopy. Specifically, the onset of t2g peaks in the Ti L-edge was observed to red shift and the eg peaks were observed to blue shift with increasing fields, attributed to changes in the hybridization of the conduction band (3d) orbitals.
Transition metal oxides possess a wide range of optical,1 electrical,2 and magnetic properties3–5 that stem from the character and occupations of the d-orbitals.6–9 These properties can be tailored using external stimuli including electric,10 optical,11 and chemical fields12 due to the changes in the energy/hybridization of the transition metal d-orbitals.13,14 Subsequently, many modern technologies such as photovoltaics,15,16 solid-state display panels,17 and non-volatile memory devices18,19 that demand charge mobility and reversibility4 are currently based on transition metal oxides. Recent experiments, for example, have shown that the electronic structure of transition metal oxides can be modulated using electric fields induced by surface ligands, which can be used to tune the optoelectronic properties.12,20 For all of these applications, understanding the relationship between applied fields and the electronic structure is critical to elucidate the physical mechanisms and ultimately design new functional inorganic materials. Experimental approaches to this include UV–Vis absorption spectroscopy,21,22 x-ray photoemission/absorption spectroscopy,20,23,24 and magnetic measurements using the superconducting quantum interference device (SQUID).22,25
Among these, x-ray absorption near-edge structure (XANES) spectroscopy has emerged as a powerful tool for probing the electronic structure of transition metal oxides26–28 due to its atomic specificity and ability to capture subtle changes in the unoccupied electronic states (conduction band)29 that result from the changes in the lattice geometry,30 oxidation state,31 band spacing,32 and band populations.33 XANES has been applied to measure changes due to weak fields, such as the effect of surface ligands on the d-orbitals of TiO2,12 as well as strong field processes such as band-tunneling and transient metallization.34 Interpreting XANES spectra, however, can be quite challenging, which necessitates first-principles simulations for relating the observed spectra to the underlying electronic structure and/or dynamics.
There are numerous methods for computing XANES spectra for molecules and solids.35–37 Semi-empirical methods such as crystal-field multiplet (CFM)38 and charge-transfer multiplet (CTM)39,40 can give transition metal spectra that match with experiments quite well12,41 but require choosing empirical crystal field parameters. Alternatively, one can use ground-state-based first-principles (ΔSCF) methods to compute a spectrum directly from the transition between the core and valence states. These ΔSCF approaches require some description of the core hole, which is typically done by constraining the occupancy.42–47 For the electronic structure method, density functional theory (DFT)43–45,48–51 is often used due to a good tradeoff between accuracy and efficiency. Multiconfigurational methods are also applied to XANES, such as complete/restricted active space SCF (CAS/RASSCF),46,52 multireference configuration interaction (MRCI),53 and multireference coupled-cluster (MRCC).54 These are especially successful for partially occupied degenerate ground states and partially filled orbitals, but care must be taken in choosing the active space. While SCF-based methods naturally capture relaxation, they may suffer from variational collapse.55 Additionally, they may require modifications to explicitly enforce orthogonality.56
Excited-state methods, which do not suffer from these limitations, are also used to compute XANES spectra. Single-determinant excited-state methods such as static exchange (STEX),57,58 linear response (LR),59–61 and real-time (RT)59,62–70 time-dependent density functional theory (TDDFT) have been quite successful but may fail for double excitations and multiplet effects and often give inaccurate absolute energies due to incomplete core-hole relaxation.59 These problems can often be remedied somewhat by using post Hartree–Fock (HF) methods, such as equation-of-motion coupled-cluster (EOM-CC)54,71 and real-time EOM-CC,72,73 which are systematically improvable but with significantly increased computational cost. Green’s function (GW) approaches, such as multiple-scattering with SCF potentials,74,75 Bethe–Salpeter-Equation (BSE typically with the GW approximation),76–78 and algebraic diagrammatic construction (ADC),79 capture relaxation well but typically require transitions to be calculated separately, which can make them inconvenient for broadband spectroscopy.80
In particular, XANES calculations of transition metal oxides typically use some form of periodic boundary conditions with either grids or a planewave basis, primarily using DFT,48 multiple-scattering,74,75 and GW/BSE.76–78 These give a good description of the band structure and reliable spectra. For DFT, however, these basis sets are usually limited to local density (LDA) and generalized gradient (GGA) approximations or Hubbard-correct versions of these functionals (LDA/GGA+U).81–83 On the other hand, low concentration of dopants and defects can be challenging, and for practical reasons, these basis sets preclude the wavefunction-based methods such as hybrid DFT or post-HF techniques. As an alternative, finite simulations can be used, where a cluster is embedded chemically/electrostatically to emulate the bulk.62,67,84 This enables the use of all-electron atom-centered basis sets, which is convenient for inner-shell spectra and allows for efficient evaluation of exchange integrals. However, they can struggle to properly represent long range interactions and can suffer from unphysical finite size effects. Some examples of previous XANES simulations using cluster models include TiO2 using DFT85 and multiple-scattering,86 CaF2 using ROCIS,85 and iron oxides using MRCI.53 Finally, relativistic effects can be significant in transition metal oxides as the spin–orbit (SO) coupling is often on the order of the peak splitting in a XANES spectrum.87 In TiO2, for example, the SO splitting in the 2p orbitals is roughly 6 eV,12,13 meaning one cannot simply do a separate LIII and LII simulation as for larger Z elements.59
TDDFT-based methods are convenient for valence properties of condensed-phase systems88–92 due to the good tradeoff between reasonable computational cost, predictability, and ability to be extended to include SO coupling.64,67,87,93 Although less common than the valence, TDDFT has also been applied to solid-state XANES including silica,62 alkaline-earth oxides,94 and titania.95 In this context, atom-centered basis sets and cluster models are especially useful as they obviate the preparation of core holes or transition potentials and allow for efficient evaluation of hybrid functionals, which have been shown to give improved bandgaps.96,97 Additionally, hybrids have been shown to give better XANES spectra in molecules vs LDA/GGAs, with 50% exact exchange giving the most accurate absolute energies.98 LR-TDDFT can be challenging due to the convergence issues of the large number of roots required for the calculation of the XANES spectra.99 Real-time methods, where the density matrix or orbitals are propagated in time, are well-suited to solid-state XANES simulations as they yield the entire spectrum from valence-to-core transitions without having root convergence issues.100,101 Another advantage of using a real-time method is the ability to simulate non-linear63,66 and time-resolved processes, such as transient absorption spectroscopy.102,103 The methods developed in this paper were done with this in mind.
This paper demonstrates the use of RT-TDDFT simulations with the inclusion of SO coupling to elucidate the effect of external fields on the electronic structure (XANES spectra) of a prototypical transition metal oxide system, i.e., anatase TiO2.104,105 First, we develop bulk-mimicking clusters for anatase and validate SO-RT-TDDFT for computing the Ti L-edge XANES spectra. Next, we apply a range of static fields and calculate the field-modified XANES spectra to better understand how observed changes in spectra can be related to the changes in the Ti d-orbital energy landscape. Ultimately, this work has implications in solar cells,15 (bio)sensing,106,107 tunable displays,108 and capacitors,109 etc., where application of external fields can be used to modify the optoelectronic properties.
All calculations were performed using a customized development version of NWChem110 with Gaussian-type orbital (GTO) basis sets and relativistic effects described using zeroth order regular approximation (ZORA).111 The B3LYP exchange-correlation functional was used for the calculations in this study as it has performed well in previous studies for ground and excited-state properties of TiO2.96,112 All the geometry optimizations and similar calculations, along with the XANES calculations, were performed with the following basis sets:113 Ti atoms were given a Def2-TZVP basis set, O atoms were given a Def2-SVP basis sets, and the pseudo-H atoms were given a 6-31G basis set. In this work, we use bulk-mimicking clusters, which is a well-developed approach for weak-field vertical valence and core-level excited states in non-metallic materials where the excitations are localized in space.114 Herein, TixOy clusters were developed using a covalent embedding procedure.84,89 The experimental bulk anatase TiO2 structure was cut to yield a Ti centered finite cluster, which was then “chemically passivated” with pseudo-hydrogen atoms at the boundaries (Fig. S1). This is done by replacing the outermost atoms in the cluster by appropriately charged H atoms. In anatase TiO2, Ti forms a distorted octahedron with the six surrounding O atoms, while O atoms form a trigonal planar structure with the three surrounding Ti atoms. Based on the formal oxidation states of Ti (+4) and O (−2), each Ti atom will share an effective charge of +2/3 and each O atom will possess an effective charge of −2/3 to their neighboring atoms. The boundary pseudo-hydrogen atoms, which replace the outermost Ti atoms, will therefore have an effective charge of +2/3. These charges are independent of the applied electric field. Additionally, by varying the O–H bond lengths, one can control how much electron density is donated or withdrawn from the cluster. The O–H bond distance was chosen as 1.0 Å based on previous work.89 The clusters were then geometry optimized with the interior atoms allowed to move, while pinning the boundary O atoms and fixing the O–H bond lengths and H–O–H angles.
Next, the XANES spectra were computed using the NWChem real-time TDDFT module65 using a two-component SO approach, similar to previous relativistic RT-TDDFT implementations.61,67,93 A manuscript detailing the validation and technical aspects of this two-component implementation is in preparation.115 In this approach, the single particle density matrix is propagated in time after broadband pulse excitation. The equation of motion in the von Neumann representation is given by
where F′(t) and P′(t) are the Fock and density matrices in the canonical basis.65 The details regarding the integration of the equation of motion and calculation of the time-advanced Fock matrix are given in Ref. 116. To save computation time, our propagator is constructed via the exponential midpoint of the extrapolated (future) Fock matrix without self-consistent iteration,
This approximation, which was tested for a few spectra, is valid due to the relatively short time steps, which are required to capture x-ray frequency spectra. The time step used for these calculations is ∆t = 7.3 × 10−4 fs (0.03 a.u.), and the total time of propagation is 9.7 fs (400 a.u.). The time step was chosen to be short enough to adequately capture the frequency according to the energy range of the XANES spectra. To compute the spectra, the system was excited using a delta-function (broadband) electric field
For every simulation, the field amplitude was taken to be κ = 0.0001 a.u, the center of the pulse was to = 20 a.u = 0.48 fs, and the width was w = 0.024 fs. denotes the polarization. This electric field is coupled into the Fock matrix in the atomic orbital basis through an external potential (V) by its product with the transition dipole matrix (D), under the assumption of the uniform electric field across the system,
This approximation is valid for relatively low x-ray frequencies studied here. For higher energies, however, quadrupole and higher terms may need to be considered.
A Padé accelerated method of Fourier transform analysis116 was employed to convert the spectra from time-domain to frequency-domain using only the contributions from the Ti 2p core spin-orbitals. This converges much more rapidly with simulation time than conventional Fourier transform of the dipole moment and also allows one to only include contributions from the L-edge in the spectrum. We checked convergence with simulation time and observed that the spectra were converged for signals longer than 8.5 fs (350 a.u.). An alternative acceleration method is to use a time-correlation function with the energy of the core-level factored out to allow for a larger time step.117
To compute the spectra, the time-dependent dipole moment is first written as a sum of occupied-virtual pair dipoles,
where i = 1, …, mocc are the virtual orbitals and a = mocc+1, …, m are the occupied orbitals. These are computed by projecting the density matrix onto the ground state molecular orbitals,
where C′(0) is the eigenvector matrix of the ground state (in the canonical basis). Now the dipole polarizability for each dipole contribution () is obtained by taking Padé approximant to the Fourier transform for each dipole contribution signal separately,
Here, and the dd subscripts denote the on-diagonal part of the polarizability tensor, e.g., xx means x dipole resulting from x polarized kick. The dipole strength function is then computed from the polarizability as
Finally, to better match experimental spectra that have lifetimes that generally decrease with increasing energy above the edge, we apply energy-dependent broadening to our spectrum separately for both LIII and LII edges. This is achieved by fitting the spectrum to a number of Lorentzian curves and then broadening those curves according to τ (E), which takes the form of an exponentially decreasing core-hole lifetime given by
where α has dimensions of inverse energy, Eo stands for edge energy and is equal to 458.1 eV for the LIII edge, and τ varies from 70 s−1 to 60 s−1. For the LII edge, Eo = 463.4 eV and τ varies from 47 s−1 to 25 s−1. Another option for energy-dependent broadening is to add a small imaginary potential to the Fock matrix with the values chosen phenomenologically, e.g., exponential in the eigenvalues in the basis of the Kohn Sham orbitals.118 This results in peaks with increasing widths as you go higher above the edge. However, since the operator is non-Hermitian, strong applied fields will cause significant ionization, which can give unphysical spectra.
Before computing the effect of electric fields on the XANES spectra, we first validate our approach for the case of anatase TiO2 without an applied field. To confirm convergence of results with the cluster size, we checked the optical gaps and orbital character for a Ti9O38H60 (107 atom) and a smaller Ti3O14H24 (41 atom) cluster (Fig. 1) carved from the large one. For convenience, we computed the optical gaps using LR-TDDFT with 10 roots. In principle, RT-TDDFT, could also be used but would require longer simulation times. The optical gap of the 107 atom cluster was computed to be 3.1 eV with the “valence band” dominated by O 2p and the “conduction band” dominated by Ti 3d orbitals. The smaller 41 atom cluster had molecular orbitals similar to that of the large cluster with an optical gap of 3.7 eV. This is an overestimate of the experimental value of 3.2 eV,119 likely due to the quantum confinement effects observed in smaller clusters. Although the L-edge spectra should not necessarily depend on the specific value of the optical gap, the value and character of the gap serve as indicators that the cluster is semiconductor-like. Based on these results, the smaller 41 atom cluster adequately mimics bulk anatase and is thus employed for all subsequent calculations. For the XANES calculations, we use our spin–orbit (SO)-RT-TDDFT version of NWChem since the LR code in NWChem does not have SO coupling.
Figure 2 shows the resulting Ti L-edge XANES spectrum of the bulk-mimicking Ti3O14H24 anatase cluster. This cluster has 370 basis functions and took approximately two days on using 80 processors to complete the time propagation. To match the experiment, we include energy-dependent broadening and then shift our XANES spectrum by +9.33 eV to account for core-hole relaxation effects that are inadequately captured by TDDFT.120,121 Overall, it is observed that there is good agreement with the experimental spectrum, including the value of crystal field splitting energy (10 Dq). The two-component SO-RT-TDDFT captures both the LIII (Ti 2p3/2 → Ti 3d) and LII (Ti 2p1/2 → Ti 3d) edges, as well as the energy splitting between them (∼6 eV). The peak splitting of the eg peak in the LIII edge is attributed to and orbitals due to the deviation of the Ti from Oh symmetry in anatase.122 However, in the simulated spectra, the peak intensities of these eg peaks are reversed, likely due to finite size effects. While RT-TDDFT would be tractable for the 9-Ti (107) atom cluster, we found that the 3-Ti (41) atom cluster already gave adequate agreement with experiment, and thus, we did not perform the larger calculation due to computational cost. It would require 3000 processors for the 9-Ti atom cluster, which has 1057 basis functions, to compute the spectra in the same amount of time. This choice of 3-Ti atom cluster is consistent with the previously reported restricted open-shell calculations of TiO2 using the B3LYP functional by Maganas and co-workers, who showed that a 3-Ti atom cluster gave nearly converged spectra.85
To elucidate the effects of applied fields on the Ti L-edge XANES spectra and the resulting d-orbitals, the system was converged in the ground state in the presence of static electric fields ranging from 0 V/nm to 0.07 V/nm applied in the x-axis direction. The largest magnitude of field roughly corresponds to half of material’s breakdown voltage (∼0.14 V/nm) for this cluster, i.e., field at which the bandgap disappears. In order to create a quasi-uniform electric field on the Ti3O14H24 cluster, two point charges, with magnitudes ranging from ±43e to 300e, were separated at a distance of ±100 Å in the x-axis of the cluster. Similar calculations with fields applied in y and z-axes were performed but did not show any significant differences. To better resolve the spectral features, the Ti L-edge XANES spectra shown in Fig. 3(a) is uniformly broadened instead of using energy-dependent broadening. Looking at the LIII edge [Fig. 3(b)], first, we notice a subtle red shift of the onset t2g peak. Second, we observe a splitting of the (#) peak at higher fields (0.07 V/nm), and third, an increase in the intensity of the (*) peak with increasing field is observed. The red shift of the t2g peak at higher fields suggests non-degeneracy of the dxy, dyz, and or dzx orbitals123 upon the application of the fields. Analogous to this, the (#) peak was also observed to split at higher fields. This phenomenon of d-orbital Stark splitting in the presence of applied static fields is similar to that of rare earth 4f orbitals, which split in the presence of electromagnetic radiation.124,125 An increase in the intensity of the (*) features with electric field amplitude is indicative of less hybridization with the p-states and, thus, higher oscillator dipole strengths.126 This is attributed to the overlap of eg with the ligand p-states that are along the axes, resulting in a variation in hybridization between the Ti ion and the surrounding ligand. Additionally, the value of crystal field splitting energy (10 Dq), which is calculated as the energy difference between the t2g and (#) peaks in the LIII edge, is observed to slightly increase when applied external fields are above a critical strength, in this case between 0.03 V/nm and 0.07 V/nm (see the supplementary material).
Qualitatively, our calculations are consistent with surface-functionalization experiments on Ni2+-doped TiO2 films using polarized ligands (μ = ±5 D).12 Due to the shallow penetration (<5 nm) of the ligand-induced fields127 and soft x rays,128 as well as the surface segregation of the Ni2+ ions, these experimental observations are essentially surface selective. In these experiments, the eg peak in the Ni LII-edge was observed to increase in intensity with the increasing ligand-induced electric field. In this regard, our calculated Ti LIII-edge XANES spectra show a similar effect, albeit at a ten times stronger field strength compared to the experiments. This is likely due to the partially filled d-orbitals in Ni2+ ions, which are strongly influenced by the field than the unoccupied d-orbitals in Ti4+.129 While we cannot directly compare the values of 10 Dq from our calculations (TiO2) and experiments (Ni2+-doped TiO2), the 10 Dq calculated from the Ti LIII-edge spectra shows a similar trend as the experiments,12,23,130 i.e., an increase in its value for fields above 0.03 V/nm.
In these ligand-bonding experiments, the Ni2+ L-edge spectra were fit to an empirical model using ligand field multiplet theory.39,131 Based on the fitting, the changes in spectra upon ligand-bonding were interpreted as a slight field-induced elongation of the axial bonds around the transition metal ions at the inorganic–organic interface.12 To differentiate between the roles of electronic effects vs geometry on the spectra, we use a single-Ti atom cluster as an extreme case of possible distortions in the presence of fields. While there was an axial bond distortion upon relaxation of this cluster, these geometry relaxation effects counteracted the electronic effects on the spectra. In other words, the shifts in the spectra due to geometry relaxations are opposite to that of the static field-induced changes (see the supplementary material). Since the single-Ti atom cluster grossly overestimates the distortions, it is reasonable to conclude that the geometry would be less perturbed in the bulk, thus having less effect on the spectra. Therefore, the observed changes in the spectra can be attributed to primarily electronic (hybridization) effects.
In summary, we have developed a spin–orbit real-time TDDFT method to compute the XANES spectra of TiO2 under the presence of external fields, which captures both the LIII and LII edges of Ti. Spin–orbit coupling is crucial for the calculation of L-edge spectra in first-row transition metal oxide systems as the coupling is of the order of a few eVs. For this purpose, bulk-mimicking anatase clusters were developed, and the bandgap was evaluated to verify the accuracy of the cluster. This finite cluster method offers some advantages over other ab initio methods as it uses all electron basis sets and allows for the use of hybrid functionals to improve the quality of the XANES spectra. This technique was used to elucidate the field-induced changes in the electronic structure of TiO2 for static electric fields varying from 0 V/nm to 0.07 V/nm. Although it is experimentally challenging to disentangle electronic effects from geometric field-induced effects, our calculations indicate that fields can modify the electronic structure without geometry distortions. Critically, in the limit that geometry relaxation effects are negligible, these changes in the d-orbital hybridization can be probed via XANES. In particular, the onset of the t2g peaks is red shifted and the eg peaks are blue shifted with increasing fields, along with an increase in the intensity of the peak. While these spectral changes are specific to the anatase TiO2 system, which has a distorted octahedral Ti4+ site (D2d), similar effects are likely to be observed for other transition metal oxides for which first-principles calculations may assist in interpretation.
The supplementary material for this article contains the geometries of 107 atom and 41 atom clusters and the effect of geometry optimization on the single-Ti cluster.
A.M.M. and P.D. developed the TiO2 cluster models. A.M.M. performed the RT-TDDFT calculations. P.D. and A.M.M. analyzed the computed XANES spectra. A.B. and M.S. assisted in the Padé analysis and energy-dependent broadening. P.D. and A.M.M. contributed equally to this work.
The data that support the findings of this study are available from the corresponding author upon reasonable request.
K.L. acknowledges support by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Atomic, Molecular and Optical Sciences, under Contract No. DE-SC0017868. J.A.D. acknowledges the National Science Foundation (NSF) under Grant No. CHE-1709902 for financial support. P.D. is thankful to the U.S. Department of Energy (DOE) under EPSCOR Grant No. DE-SC0012432 for financial support and the Graduate School of Louisiana State University for the Dissertation Year Fellowship. This research was conducted with high performance computational resources provided by Louisiana State University (http://www.hpc.lsu.edu). We also acknowledge the support of Dr. Orhan Kizilkaya and the staff of the CAMD synchrotron light source.