Simulated field-modulated x-ray absorption in titania

Transition metal oxides possess a wide range of optical,¹ electrical,² and magnetic properties^3–5 that stem from the character and occupations of the d-orbitals.^6–9 These properties can be tailored using external stimuli including electric,¹⁰ optical,¹¹ and chemical fields¹² 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,¹⁷ and non-volatile memory devices^18,19 that demand charge mobility and reversibility⁴ 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 oxides^26–28 due to its atomic specificity and ability to capture subtle changes in the unoccupied electronic states (conduction band)²⁹ that result from the changes in the lattice geometry,³⁰ oxidation state,³¹ band spacing,³² and band populations.³³ XANES has been applied to measure changes due to weak fields, such as the effect of surface ligands on the d-orbitals of TiO₂,¹² as well as strong field processes such as band-tunneling and transient metallization.³⁴ 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)³⁸ and charge-transfer multiplet (CTM)^39,40 can give transition metal spectra that match with experiments quite well^12,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),⁵³ and multireference coupled-cluster (MRCC).⁵⁴ 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.⁵⁵ Additionally, they may require modifications to explicitly enforce orthogonality.⁵⁶ 

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.⁵⁹ 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),⁷⁹ capture relaxation well but typically require transitions to be calculated separately, which can make them inconvenient for broadband spectroscopy.⁸⁰ 

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,⁴⁸ 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 TiO₂ using DFT⁸⁵ and multiple-scattering,⁸⁶ CaF₂ using ROCIS,⁸⁵ and iron oxides using MRCI.⁵³ 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.⁸⁷ In TiO₂, for example, the SO splitting in the 2p orbitals is roughly 6 eV,^12,13 meaning one cannot simply do a separate L_III and L_II simulation as for larger Z elements.⁵⁹ 

TDDFT-based methods are convenient for valence properties of condensed-phase systems^88–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,⁶² alkaline-earth oxides,⁹⁴ and titania.⁹⁵ 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.⁹⁸ LR-TDDFT can be challenging due to the convergence issues of the large number of roots required for the calculation of the XANES spectra.⁹⁹ 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-linear^63,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 TiO₂.^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,¹⁵ (bio)sensing,^106,107 tunable displays,¹⁰⁸ and capacitors,¹⁰⁹ etc., where application of external fields can be used to modify the optoelectronic properties.

All calculations were performed using a customized development version of NWChem¹¹⁰ with Gaussian-type orbital (GTO) basis sets and relativistic effects described using zeroth order regular approximation (ZORA).¹¹¹ 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 TiO₂.^96,112 All the geometry optimizations and similar calculations, along with the XANES calculations, were performed with the following basis sets:¹¹³ 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.¹¹⁴ Herein, Ti_xO_y clusters were developed using a covalent embedding procedure.^84,89 The experimental bulk anatase TiO₂ 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 TiO₂, 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.⁸⁹ 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 module⁶⁵ 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.¹¹⁵ 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.⁶⁵ 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⁻⁴ 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 t_o = 20 a.u = 0.48 fs, and the width was w = 0.024 fs. $d^=x,y,z$ 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 analysis¹¹⁶ 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.¹¹⁷ 

To compute the spectra, the time-dependent dipole moment is first written as a sum of occupied-virtual pair dipoles,

where i = 1, …, m_occ are the virtual orbitals and a = m_occ+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 (⁠$αddia$⁠) is obtained by taking Padé approximant to the Fourier transform for each dipole contribution signal separately,

Here, $d^=x,y,z$ 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 L_III and L_II 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, E_o stands for edge energy and is equal to 458.1 eV for the L_III edge, and τ varies from 70 s⁻¹ to 60 s⁻¹. For the L_II edge, E_o = 463.4 eV and τ varies from 47 s⁻¹ to 25 s⁻¹. 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.¹¹⁸ 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 TiO₂ without an applied field. To confirm convergence of results with the cluster size, we checked the optical gaps and orbital character for a Ti₉O₃₈H₆₀ (107 atom) and a smaller Ti₃O₁₄H₂₄ (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,¹¹⁹ 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 Ti₃O₁₄H₂₄ 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 L_III (Ti 2p_3/2 → Ti 3d) and L_II (Ti 2p_1/2 → Ti 3d) edges, as well as the energy splitting between them (∼6 eV). The peak splitting of the e_g peak in the L_III edge is attributed to $dz2$ and $dx2−y2$ orbitals due to the deviation of the Ti from O_h symmetry in anatase.¹²² However, in the simulated spectra, the peak intensities of these e_g 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 TiO₂ using the B3LYP functional by Maganas and co-workers, who showed that a 3-Ti atom cluster gave nearly converged spectra.⁸⁵ 

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 Ti₃O₁₄H₂₄ 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 L_III edge [Fig. 3(b)], first, we notice a subtle red shift of the onset t_2g peak. Second, we observe a splitting of the $dz2$(#) peak at higher fields (0.07 V/nm), and third, an increase in the intensity of the $dx2−y2$ (^*) peak with increasing field is observed. The red shift of the t_2g peak at higher fields suggests non-degeneracy of the d_xy, d_yz, and or d_zx orbitals¹²³ upon the application of the fields. Analogous to this, the $dz2$(#) 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 $dx2−y2$ (^*) features with electric field amplitude is indicative of less hybridization with the p-states and, thus, higher oscillator dipole strengths.¹²⁶ This is attributed to the overlap of e_g 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 t_2g and $dz2$ (#) peaks in the L_III 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 Ni²⁺-doped TiO₂ films using polarized ligands (μ = ±5 D).¹² Due to the shallow penetration (<5 nm) of the ligand-induced fields¹²⁷ and soft x rays,¹²⁸ as well as the surface segregation of the Ni²⁺ ions, these experimental observations are essentially surface selective. In these experiments, the e_g peak in the Ni L_II-edge was observed to increase in intensity with the increasing ligand-induced electric field. In this regard, our calculated Ti L_III-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 Ni²⁺ ions, which are strongly influenced by the field than the unoccupied d-orbitals in Ti⁴⁺.¹²⁹ While we cannot directly compare the values of 10 Dq from our calculations (TiO₂) and experiments (Ni²⁺-doped TiO₂), the 10 Dq calculated from the Ti L_III-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 Ni²⁺ 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.¹² 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 TiO₂ under the presence of external fields, which captures both the L_III and L_II 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 TiO₂ 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 t_2g peaks is red shifted and the e_g peaks are blue shifted with increasing fields, along with an increase in the intensity of the $dx2−y2$ peak. While these spectral changes are specific to the anatase TiO₂ system, which has a distorted octahedral Ti⁴⁺ site (D_2d), 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 TiO₂ 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.