Time-resolved near-edge X-ray absorption fine structure of pyrazine from electronic structure and nuclear wave packet dynamics simulations

X-ray absorption spectroscopy (XAS) is one of the oldest methods for structural analysis, continuously employed from the dawn of quantum mechanics. Since Sayers and co-workers have shown that the fine structures in XAS can be inverted, based on Fourier analysis,¹ to yield geometric structural information, two important structural analysis techniques emerged out of XAS: near-edge X-ray absorption fine structure (NEXAFS), also known as X-ray absorption near-edge structure (XANES), and the extended X-ray absorption fine structure (EXAFS). EXAFS, which uses the high-energy part where intramolecular electron scattering by the nuclei takes place, is a method for the analysis of local geometry.² NEXAFS, which uses the low-energy part and is therefore dominated by core-valence transitions and core-quasi-bound continuum transitions, provides information on the local electronic structure near the X-ray absorbing atoms. EXAFS has matured in both theory and experiment.^2–4 Theoretical development of NEXAFS is ongoing, while it is widely applied for structural analysis of condensed matter.^5–9 

The idea of using time-resolved (TR-)XAS for tracking dynamical processes, such as intramolecular charge transfer, structural dynamics, solvation dynamics, and catalytic processes, has been around for a long time. However, it has only become feasible with the advent of third generation synchrotrons. Now, TR-XAS with microsecond to picosecond resolutions is routinely used to study crystals, bulk amorphous materials, and molecules in liquids.^7–11 Using TR-XAS with subpicosecond resolution to study the light-induced spin switching dynamics in an iron(ii)-based complex in aqueous solution has been reported.¹² In parallel to the advances at synchrotrons, X-ray free electron lasers (XFELs) were developed.^8,13–19 XFELs provide pulses with high intensity and femtosecond time resolution, enabling TR-XAS not only with shorter time resolution but also for noncrystalline biological systems and isolated gas-phase molecules. Exploiting these advantages of the XFELs, photochemical processes of many transition metal complexes have been investigated with TR-XAS, especially TR-NEXAFS.^20–23 Similar investigations for organic molecules have also begun to emerge.²⁴ The development of table-top laser techniques for X-ray pulses with high intensity and femtosecond time resolution based on high harmonic generation (HHG)²⁵ is further expanding the capabilities of TR-NEXAFS. For example, TR-NEXAFS employing this technique for some organic molecules has been reported.^26–28 

Extracting the mechanistic information from the TR-NEXAFS data relies on computational modeling. The first theoretical TR-NEXAFS spectra were reported by Penfold and co-authors,^29,30 who employed time-dependent density functional theory (TDDFT) and the multiconfiguration time-dependent Hartree (MCTDH) method.^31–34 Later, Neville and co-authors computed the spectra^35–37 with the ADC(2) (second-order algebraic diagrammatic construction) method^38,39 and the ab initio multiple spawning (AIMS) method.⁴⁰ Wolf et al.²⁴ presented experimental results for the TR-NEXAFS of thymine at the oxygen K-edge, interpreted by means of accurate coupled cluster spectra calculations at selected excited-state geometries.²⁴ Attar and co-authors calculated the time-dependent differential absorption of NEXAFS²⁷ with TDDFT and the fewest-switches surface-hopping method⁴¹ to characterize experimental spectra obtained for an electrocyclic ring-opening reaction. Recently, Northey and co-authors^42,43 applied the variational multiconfigurational Gaussian (vMCG) method^44,45 for the nuclear dynamics part of TR-NEXAFS simulation. These calculated spectra show the sensitivity of the pre-edge regions, which reflect core-to-valence electronic transitions, to both the molecular electronic state and the molecular geometry at the instant of X-ray absorption.^35–37 Neville and co-authors also computed the time-resolved photoelectron spectra (TR-PES) to demonstrate that the contribution in TR-NEXAFS of each initial electronic state at the probe process is more clearly resolved, as compared to the TR-PES spectra where many ionized final states contribute. Thus, TR-NEXAFS is a promising tool, which offers capabilities beyond those of ultraviolet absorption, fluorescence, and phosphorescence spectroscopies, the only methods that were available so far for investing photophysical processes.

In this contribution, we present further investigation of the spectroscopic signatures of electronic and nuclear dynamics involved in a typical photophysical process. As a benchmark system, we use pyrazine, a prototypical heterocyclic molecule whose photophysics had been rigorously investigated.^46–51 We report the computed TR-NEXAFS spectra at both carbon and nitrogen K-edges of pyrazine excited from the S₀(¹A_g) state to S₂(¹B_2u(ππ^*)) state with 4.79 eV (259 nm) excitation. The system rapidly decays to the lowest-lying S₁(¹B_3u(nπ^*)) state via internal conversion (IC) driven by the ν_10a vibrational mode. Previous simulations of nuclear dynamics indicate that the lowest-lying dark state, ¹A_u(nπ^*), is involved in the relaxation.^48–51 We first characterize the population dynamics between these three electronic states and investigate the mechanism of the electronic transitions by visualizing the time-resolved reduced nuclear densities obtained with the MCTDH method. We then present the NEXAFS spectra (at the carbon and nitrogen K-edges) for the ground state and the electronic excited states, both computed within the frozen-core/core-valence-separated equation-of-motion coupled cluster singles and doubles (fc-CVS-EOM-CCSD) method.⁵² Finally, we simulate the two-dimensional TR-NEXAFS spectra by summing the spectra at selected time delays of the excited electronic states weighted by their populations.

In the (equation-of-motion) coupled-cluster ansatz, the ground-state wave function is parameterized as

where |Φ₀⟩ is the Hartree-Fock (HF) reference state and T is the cluster operator, with τ_μ being the excitation operators and t_μ the corresponding cluster amplitudes. Here, the indices {i, j, k, …} and {a, b, c, …} refer to occupied and virtual orbitals in the HF reference, respectively. In the following, the indices {p, q, r, …} denote generic orbitals that may be either occupied or unoccupied.

The CC ground-state energy can be obtained from the projected Schrödinger equation,

and the cluster amplitudes satisfy the CC equations

where ⟨Φ_μ|’s represent μ-tuply excited determinants and $H ¯$ is the similarity transformed Hamiltonian,

At the coupled-cluster singles and doubles (CCSD) level, which is employed in this study, the cluster operator T is truncated to the single and double terms indicated in Eq. (3).

The (right) excited states are generated by application of excitation operators R to the ground-state wave function,^53–56 

Since the excitation operators T and R commute, the Schrödinger equation for the excited states can be recast as

By projecting the Schrödinger equation (9) onto ⟨Φ_μ|, it reduces to an eigenvalue problem of the similarity transformed Hamiltonian matrix. Due to the non-Hermiticity of the CC operator T, each root of $H ¯$ is associated with two eigenvectors corresponding to distinct bra and ket states

where L is the de-excitation operator,

The EOM amplitudes r_μ and l_μ are found as stationary points of the EOM functional,

By applying the bivariational principle,^57,58 one arrives at the nonsymmetric eigenvalue problem,

where the eigenvectors of the Hamiltonian are chosen to form a biorthogonal set,⁵⁹ 

At the EOM-CCSD level, the EOM amplitudes are found by diagonalization of effective Hamiltonian $H ¯$ in the basis of the reference and of singly and doubly excited determinants, i.e., solving the following right and left eigenvalue equations:

where ω_i is the energy difference of state i with respect to the reference state (i.e., excitation energy).

To compute the TR-NEXAFS spectra, one needs the energy differences and the dipole moment matrix elements between the valence-excited state involved in the dynamics and the set of final core-excited states generated by the probe pulse. Toward this goal, we use the fc-CVS-EOM-CCSD formalism⁵² in which the CCSD ground-state and the valence-excited EOM target states are obtained using the frozen-core approximation and the core-excited states are computed using the core-valence separation (CVS). CVS entails imposing a restriction that at least one core orbital is involved in R and L amplitudes.⁶⁰ The oscillator strengths between the valence- and core-excited states are then computed according to

where the transition dipole moments $μ α i → j$ and $μ α j → i$ are contractions of the electric dipole property integral matrices with the transition density matrices between state i and j,

In Eq. (19), μ_α refers to a specific Cartesian component of the electric dipole operator, the state i is the valence-excited state and j is the core-excited state. The transition density matrix elements $γ p q i → j$ and $γ p q j → i$ are given by

Applying a singular value decomposition (SVD) to the transition matrices yields natural transition orbitals (NTOs),^61–64 which can be used to describe the electronic excitations in terms of hole-particle pairs,

where $ϕ ̃ K e ( i → j )$ and $ϕ ̃ K h ( i → j )$ are particle and hole orbitals obtained by SVD of the transition density matrix γ^i→j, and $σ K i → j$ are the corresponding singular values. Note that, since the SVD is unitary transform of γ^i→j, the sum of the square of $σ K i → j$ is equal to the squared norm of γ^i→j, which provides a simple metric quantifying the one-electron character of the transition, i.e., for pure single excitations, ∥γ^i→j∥ = 1.

MCTDH is an efficient method to solve the molecular time-dependent Schrödinger equation in multiple dimensions. To carry out the MCTDH calculations,^31–34 we utilize the vibronic-coupling (VC) Hamiltonian in the diabatic electronic basis. In this work, we adopt the Hamiltonian of Ref. 49 that was constructed with knowledge of the excited-state potential energy surfaces along the normal vibrational modes. It is composed of the reference Hamiltonian $H 0 ( v c ) ( Q )$⁠, and the matrix W^(vc)(Q), which expresses the changes in the excited state potential energies with respect to the ground state as a Taylor expansion around the ground state equilibrium geometry [the Franck-Condon (FC) geometry],

where n ≠ n′, with n being the electronic state index. The vector Q represents dimensionless mass-frequency weighted normal coordinates. The reference Hamiltonian $H 0 ( v c ) ( Q )$ is the harmonic vibrational Hamiltonian for the ground state in dimensionless coordinates, where Ω_i are the harmonic vibrational frequencies and I is the identity matrix of the matrix space spanned by the electronic basis set. The diagonal parts of the matrix W^(vc)(Q) are expressed in terms of the vertical excitation energies E_n and the linear (⁠$κ i ( n )$⁠) and quadratic (⁠$Γ j ( n )$⁠) intrastate coupling constants for the nth electronic state. The constant $λ i ( n n ′ )$ in the off-diagonal elements of W^(vc)(Q) couples the nth and n′th electronic states. We utilize the 9-mode 3(4)-state vibronic Hamiltonian of Ref. 49, derived from ab initio data obtained at the XMCQDPT2⁶⁵/aug-cc-pVDZ level.

From a total of 24 normal vibrational modes, we selected 4 totally symmetric modes with large intrastate couplings (Fig. 1) and 5 other modes that strongly couple the 3 electronic states as active modes. The harmonic vibrational frequencies, Ω_i, for the 9 active modes are collected in Table I.

The wave function for the total system in the multiset variant⁶⁶ of the MCTDH method^31–34 is given by

where {|α⟩} denotes the set of diabatic electronic states. The nuclear wave function Ψ^(α)(Q, t) in the electronic state |α⟩ is expanded on the basis of time-dependent functions φ(Q_κ, t) called single-particle functions (SPFs),

The SPFs are further expanded in a primitive basis, built of one-dimensional basis functions for each degree of freedom,

For the calculation of TR-NEXAFS spectra, the populations of the diabatic electronic states, $P α d ( t ) =∫| Ψ ( α ) ( Q , t ) | 2 dQ$⁠, and the expectation values $Q ¯ ( α ) ( t )$ of the position operator [that is, the centroid of the nuclear density |Ψ^(α)(Q, t)|² of the α state], given by $Q ¯ ( α ) ( t ) =∫Q| Ψ ( α ) ( Q , t ) | 2 dQ/ P α d ( t )$⁠, are used.

Following excitation to ¹B_2u(ππ^*), the quantum state of pyrazine evolves into a coherent superposition of the three diabatic states, ¹B_3u(nπ^*), ¹A_u(nπ^*), and ¹B_2u(ππ^*), as per Eq. (26). The ensuing dynamics is probed by X-ray pulses at different time delays. The measured probe intensities are modeled as an incoherent superposition of the core excitations from the ¹B_3u(nπ^*), ¹A_u(nπ^*), and ¹B_2u(ππ^*) states, weighted by the respective oscillator strengths multiplied by the nuclear densities at the instant of the X-ray probe, averaged over the nuclei coordinates,

where ω is the photon energy, α and α^* refer to the initial and final electronic states, and ω_α(Q) and $ω α * ( Q )$ are their respective electronic energy eigenvalues at the nuclear coordinates Q. In the present work, we replace the nuclear densities |Ψ^(α)(Q, t)|² by the δ-functions centered at $Q ¯ ( α ) ( t )$ for simplicity. Then, the expression for the TR-NEXAFS intensities is recast as

Thus, to compute TR-NEXAFS spectra, we first calculate the populations and the centroids of the nuclear wave packets (i.e., the expectation values of the geometries) for the ¹B_3u(nπ^*), ¹A_u(nπ^*), and ¹B_2u(ππ^*) states. We then calculate the energies and oscillator strengths for core-level excitations of these three electronic states at the centroid geometries, weight them by the populations, and sum them together. We note that Eqs. (29) and (30) are commonly used for calculating TR-XAS and TR-PES (see, for example, Refs. 67 and 68). Implications of the approximation from Eqs. (29) to (30) are explained in  Appendixes A and  B.

Because quantum chemistry yields adiabatic states, the state characters can change at different geometries. Consequently, energetic ordering of diabatic electronic configurations can switch depending on the molecular geometry. Yet, the quantum nuclear dynamics simulations use a diabatic basis, i.e., the one which preserves the electronic character, in order to avoid numerical divergence of the VC potential expressed in the adiabatic basis near avoided crossings.⁶⁹ In the present study, since only the centroids of the nuclear wave packets (always D_2h symmetry) are sampled, adiabatic and diabatic states are the same at all sampled points. We ensured the consistency of the electronic states between the electronic structure spectral calculations and quantum dynamics simulations by inspection of the corresponding electronic configurations characterized by NTOs. Hereafter, we frame our discussion in terms of the diabatic electronic basis. The range of validity of this treatment is discussed in  Appendix B.

We note that in the present work, the excitation energies of the initial and final excited states and the oscillator strengths between them were calculated by fc(-CVS)-EOM-CCSD, whereas the vibronic Hamiltonian of Ref. 49 was derived from the XMCQDPT2 data. The energetic ordering of the ¹A_u(nπ^*) and ¹B_2u(ππ^*) states at the FC region is inconsistent between the two levels of theory. However, the geometry optimized for the ground state at the CCSD level agrees well with the one optimized at the MP2 level reported in Ref. 49. The geometries optimized for the ¹B_3u(nπ^*), ¹A_u(nπ^*), and ¹B_2u(ππ^*) states at the EOM-CCSD level are also in reasonable agreement with the potential energy minima at the XMCQDPT2 level, which were derived from the vibronic Hamiltonian and geometrical parameters for the ground state in Ref. 49 (see Table II). Therefore, if the character of the electronic configurations is consistent throughout the photophysical process, it is reasonable to assume that the difference between the theoretical levels used in the core-excitation and construction of the vibronic Hamiltonian should not qualitatively affect the resulting TR-NEXAFS.

We end this section with a few notes specific to pyrazine. Pyrazine belongs to the D_2h point group, and the principal axis is conventionally set along the N–N direction, with the y-axis set in the molecular plane.⁴⁷ The vertical transitions to the adiabatic S₁ and S₂ states⁴⁶ correspond to the transitions to the diabatic ¹B_3u(nπ^*) and ¹B_2u(ππ^*) states,⁴⁷ respectively.

Nuclear quantum dynamics simulations were performed using the Heidelberg MCTDH package⁷⁰ with the model Hamiltonian from Ref. 49, described in Sec. II B. The vibrational modes were combined, a number of SPFs were set, and the grid points were fixed, exactly as given in Table 6 of Ref. 49. The populations $P α d ( t )$ and the centroids of the nuclear densities $Q ¯ ( α ) ( t )$ in Eq. (30) were extracted from the MCTDH outputs. MCTDH outputs give the geometries $Q ¯ ( α ) ( t )$ as displacement from the FC geometry in the dimensionless mass-frequency weighted normal coordinates. The transformation to Cartesian coordinates is given by the eigenvectors of the mass-weighted Hessian matrix.

To obtain the energy eigenvalues of the initial and final core-excited states in Eq. (30), $ω α ( Q ¯ ( α ) ( t ) )$ and $ω α * ( Q ¯ ( α ) ( t ) )$⁠, and the oscillator strength $f ( α → α * ; Q ¯ ( α ) ( t ) )$⁠, the fc-EOM-CCSD and fc-CVS-EOM-CCSD⁵² calculations were carried out using the ccman2 module of the Q-Chem electronic structure package,⁷¹ exploiting the libtensor library.⁷² All spectra were calculated using Pople’s 6-311++G^** basis set. For the calculation of the ground-state NEXAFS spectrum, the basis set was further augmented with uncontracted Rydberg-type functions whose exponents were generated according to the prescription of Kaufmann et al.⁷³ Each spectral stick was broadened with a Lorentzian function (FWHM = 0.4 eV), as was done in Ref. 52 to match the experimental peak widths. Finally, the spectra for the relevant electronic excited states were weighted by the populations and summed together.

To model excited-state dynamics following the photoexcitation, we simply project the nuclear wave packet on the ¹B_2u(ππ^*) potential energy surface at the FC geometry and initiate the nuclear wave packet propagation. The resulting population dynamics, which is illustrated in Fig. 2, agrees well with that of Ref. 49. The latter was obtained by explicitly taking into account the interaction with the pump pulse in the way written in Refs. 74 and 75. The populations in Fig. 2 indicate that the nuclear wave packet in the ¹B_2u(ππ^*) state immediately starts to decay into both the ¹B_3u(nπ^*) and ¹A_u(nπ^*) states. Then, the nuclear wave packet undergoes oscillatory transitions between the ¹B_3u(nπ^*) and ¹A_u(nπ^*) states in a period of ∼60 fs. To interpret this population dynamics, we show in Figs. 3 and 4, the nuclear density reduced to one-dimensional function like |Ψ^(κ,α)(Q_κ, t)|² ≡ ∫|Ψ^(α)(Q, t)|²dQ₁⋯dQ_κ−1dQ_κ+1⋯, along the normal coordinates of the active totally symmetric vibrational modes, Q_6a, Q₁, Q_9a, and Q_8a.

Figure 3 shows the slice cuts of the diabatic potential energy surfaces at the FC geometry and the evolution of the reduced nuclear densities up to 30 fs. The initial dynamics of the nuclear wave packet created on the ¹B_2u(ππ^*) surface proceeds along the Q_6a, Q₁, Q_9a, and Q_8a normal mode coordinates. In all displayed slice cuts of the potential energy surfaces, the ¹B_2u(ππ^*) state curve intersects with that of the ¹A_u(nπ^*) state and the reduced nuclear density always surmounts or reaches the conical intersections. Therefore, the nuclear wave packet in the ¹B_2u(ππ^*) state is always in the vicinity of the conical intersection with the ¹A_u(nπ^*) state, and it continually decays to the ¹A_u(nπ^*) state, despite relatively small interstate couplings [relative to the coupling between the ¹B_3u(nπ^*) and ¹B_2u(ππ^*) states, see Table 5 of Ref. 49]. In the course of the decay, the nuclear wave packet created on the ¹A_u(nπ^*) state evolves, while new amplitudes continue to pour in the region of −2 ≤ Q_6a ≤ 0 and −1 ≤ Q_8a ≤ 2 from ¹B_2u(ππ^*). This locality of the internal conversion results in a split shape of the reduced nuclear density on the ¹A_u(nπ^*) state. In contrast, the transition from ¹B_2u(ππ^*) to ¹B_3u(nπ^*) occurs in a broad range of nuclear coordinates consistently with the large ¹B_2u–¹B_3u coupling. Thus, the nuclear wave packet on the ¹B_3u(nπ^*) surface is created below that on the ¹B_2u(ππ^*) surface. At 10 fs, another transition path from ¹A_u(nπ^*) to ¹B_3u(nπ^*) via the conical intersection along Q_8a opens up. Therefore, our calculations visualize an earlier finding^48–51 that two paths contribute to the decay to the ¹B_3u(nπ^*) state: the direct one and the one via ¹A_u(nπ^*).

The decay process of the ¹B_2u(ππ^*) population is nearly complete at 40 fs. The evolution of the reduced nuclear densities after the decay process is shown in Fig. 4. The nuclear wave packets on the ¹B_3u(nπ^*) and ¹A_u(nπ^*) surfaces continue to propagate along the normal coordinates of the totally symmetric modes, undergoing transitions via the conical intersection along Q_8a. The populations oscillate with a period of about 60 fs, which approximately corresponds to 593 cm⁻¹, the frequency of the ν_6a mode. The transition of ¹B_3u(nπ^*) →¹A_u(nπ^*) occurs at Q_6a > 0, and the inverse occurs at Q_6a < 0. These transitions are driven by the nuclear wave packet motion along Q_6a, which can be rationalized by considering slice cuts of the potential energy surfaces along Q_8a. As shown in Fig. 5, in slice cuts of Q_6a > 0, the relative energy difference between the minima of the ¹B_3u(nπ^*) and ¹A_u(nπ^*) states is larger than in Q_6a < 0. Consequently, the transitions from ¹B_3u(nπ^*) to ¹A_u(nπ^*) occur at Q_6a > 0 and vice versa. In this mechanism, the oscillatory nuclear wave packet evolution along Q_6a results in oscillation of the ¹B_3u(nπ^*) and ¹A_u(nπ^*) populations. In a similar fashion, the small fluctuations of the populations with a period of about 30 fs are due to the transitions induced by the motion along Q_9a. If the potential energy surface of the ¹A_u(nπ^*) state is higher than in the present model, the population of the ¹B_3u(nπ^*) state would be larger, which is exactly the case in Ref. 49: In that work, the vertical excitation energy for the ¹A_u(nπ^*) state was empirically adjusted from 4.45 to 4.69 eV to reproduce the UV absorption spectrum.

The electronic configuration of the ground state of pyrazine, which has the D_2h structure at the FC geometry optimized at the MP2/aug-cc-pVDZ level, is

where

Pyrazine has three occupied π orbitals: 1b_3u, 1b_2g, and 1b_1g, which are similar to the 1a_2u and the degenerate 1e_1g orbitals of benzene. The three virtual π^* orbitals, 2b_3u, 1a_u, and 2b_2g, are similar to the degenerate 1e_2u and 2b_2g orbitals of benzene, respectively. The lowest excited states are dominated by the excitations involving these π, π^*, and 6a_g (n) orbitals.

Figure 6 illustrates the characters of the ¹A_g → ¹B_3u(nπ^*), ¹A_g → ¹A_u(nπ^*), and ¹A_g → ¹B_2u(ππ^*) transitions. The symmetries of NTOs are the same as those of the respective canonical HF orbitals used above to describe the electronic configuration of pyrazine, and their shapes are very similar. The transition between the ground state and the ¹B_3u(nπ^*) state is dominated by 6a_g → 2b_3u. The optically forbidden transition to the ¹A_u(nπ^*) state from the ground state is dominated by 6a_g → 1a_u. The transition from the ground state to the ¹B_2u(ππ^*) state has contributions from two pairs of NTOs: 1b_1g → 2b_3u (dominant) and 1b_2g → 1a_u (minor).

Figure 7 shows the computed carbon K-edge NEXAFS of the ground electronic state (solid line) with the one obtained from the experiment (dashed line, redigitized from Ref. 76). The ionization energy (IE) corresponding to the computed spectrum is 292.9 eV. The NTOs for the underlying transitions are summarized in Fig. 8. Figure 9 shows a simplified MO scheme of the excitation processes responsible for the first two peaks in the spectrum.

The two most intense peaks (at 286.6 eV and 288.1 eV) are due to almost equal contributions from the transitions to the two π^* orbitals: 2a_g → 2b_3u and 1b_3g → 1a_u. The remaining two peaks, one at 289.6 eV dominated by 1b_2u → 10a_g transition, with a minor contribution from 2a_g → 7b_2u, and the other at 290.7 eV dominated by 2a_g → 6b_3u transition, with a minor contribution from 1b_2u → 2b_1g.

The experimental carbon K-edge NEXAFS spectrum obtained by synchrotron radiation has been reported by Vall-llosera et al.,⁷⁶ along with the spectrum calculated with transition potential (TP) DFT. While our computed spectrum agrees well with the experimental one at higher energies, there is a discrepancy in both intensity and position of the first two peaks. Our computed intensity ratio between the first and second peaks is less than 0.1, in contrast to a 0.5 ratio in the experimental spectrum. In addition, our computed energy separation between the first two peaks is 1.5 eV, which is three times larger than the experimental value of 0.5 eV. We confirmed that this discrepancy is not due to vibrational effects, by calculating the spectra of the ground state at the XMCQDPT2 potential minima of the lowest-lying valence-excited states, which may be close to those of the final core-excited states.

To find the source of this discrepancy, we carried out additional calculation at slightly distorted geometry such that the character of each transition can be described by a single NTO pair rather than by two NTO pairs, which are delocalized due to symmetry restrictions. Specifically, we computed the spectrum and the NTOs at an artificially distorted geometry, in which the nitrogen atoms were displaced in the opposite y-direction by 0.01 Å, and one of the C–C bond distances was shortened by 0.02 Å. The resulting spectrum was not much different from the one shown in Fig. 7. The NTOs, on the other hand, transform into two sets of four equivalent hole-particle pairs, which are shown in Fig. 10. From inspection of these NTOs, one can clearly see that peak 1 corresponds to transitions from a localized core orbital on one carbon atom to a π^* orbital with some localized density on the same carbon atom. Peak 2, on the other hand, contains transitions from a localized carbon core orbital to a π^* orbital with electron density delocalized on a C–N bond. This results in lower intensity of the second peak. Thus, the difference in the intensities and the large separation of the peaks arise due to the localization of the excited electron near the core hole, which optimizes the Coulomb interaction between the hole and the electron. We verified [via CCSDR(3) calculations⁷⁷] that the large separation between the two peaks is reduced if triple excitation effects are included in the calculation of the excitation energies. Interestingly, both the peak separation and the intensities of the first two peaks obtained with TP-DFT⁷⁶ agree better with the experimental spectrum. The TP-DFT method describes the final core-excited states by means of single electronic configurations. In addition, we checked that the first and second peaks obtained with TDDFT/B3LYP are also in better agreement with those in the experimental spectrum, and the corresponding core-excited states are dominated by the 2a_g → 2b_3u, and 1b_3g → 1a_u configurations, respectively. Therefore, localization of the excited electron near the core hole at the TP-DFT/TDDFT levels is not as strong as obtained here with the fc-CVS-EOM-CCSD method, which (accidentally) results in better agreement with the experiment. The possible reasons for this are (i) strongly delocalized 1b_3u(π) electrons screen the Coulomb potential from the core hole and (ii) the 2b_3u and 1a_u orbitals of π^* character are energetically close, but not completely degenerate. As a comparison, in the case of C_1s excitations from benzene, the π^* orbitals, corresponding to the 2b_3u and 1a_u orbitals of pyrazine, are completely degenerate, so the main peak in the spectrum does not split, and localization effects have no influence on it.⁷⁸ Splitting of the main peak in the carbon K-edge spectra may be general for aromatic compounds, such as other azabenzenes or substituted benzenes, and reflects the level of symmetry of the compounds.

Before analyzing the time-resolved spectra, let us examine the character of the spectra corresponding to the core-level transitions from electronically excited valence states below the IE.

In the fc-CVS-EOM-CCSD approach⁵² used here for calculating the valence-excited state’s NEXAFS, it is assumed that the final core-excited states of each valence-excited state belong to the manifold of the core-excited states generated by excitation from the ground state. However, since one of the valence electrons is already excited at the instant of X-ray absorption, many of such final excited states cannot be connected by one-electron excitation from the ground state. In other words, these target EOM states are multiply excited with respect to the ground-state reference. Consequently, the quality of description of these states depends on the level of theory. Importantly, the lack of triple excitations in the EOM-CCSD anstaz used here results in deterioration in quality of doubly excited target states.

Alternatively, calculations could be performed using a reference state that has the appropriate electronic configuration (the one of the valence-excited state). This is commonly achieved^{27,28,42,43,79} by using an open-shell reference determinant, which can be computed using the maximum overlap method (MOM) algorithm.⁸⁰ For example, in Ref. 43, the TDDFT N K-edge spectrum for pyrazine in the ¹B_3u(nπ^*) state was obtained using the M_s = 0 nπ^* open-shell MOM Kohn-Sham determinant, which corresponds to the half of a properly spin-adapted excited-state wave function. However, using this strategy with the CCSD ansatz often leads to convergence issues.

To overcome this convergence problem, we utilized instead the lowest-lying high-spin open-shell reference with the same orbital character as one of the target excited states, specifically ¹B_3u(nπ^*). This is illustrated in Fig. 11.

We used this high-spin $( 6 a g ) 1 3 ( 2 b 3 u ) 1$ reference to compute the XAS spectrum at the optimized geometry of the ³B_3u(nπ^*) state. Figure 12 shows the resulting carbon K-edge NEXAFS of the lowest-lying high-spin open-shell reference of nπ^* character. The spectrum is an approximation of the NEXAFS spectrum of the valence-excited ¹B_3u(nπ^*) state over the entire frequency region below the IE.

Partial occupation of the 2b_3u orbital in the initial electronic state hinders the mixing between the 2b_3u and 1a_u orbitals. As a result, the gap between the shoulder and the most intense peak drops to 0.9 eV. The excited-state spectrum almost completely overlaps with the ground-state spectrum in the region between 289 and 291 eV. The decrease of intensity near the IE is likely due to the contributions of the doubly excited configurations, $( 2 a g ) 1 ( 6 a g ) 1 ( 2 b 3 u ) 1 ( 3 b 3 u ) 1$⁠, to the target core-excited states. The peak at 283.8 eV corresponds to the C_1s(β) → 6a_g(β) transition, i.e., the excitation of the core electron into the singly occupied molecular orbital (SOMO).

The core-to-SOMO transition appears as a single peak in this spectrum, whereas it splits into two in the spectrum of the ¹B_3u(nπ^*) state, computed from the reference state that has the ground-state electronic configuration [see Fig. 13(a)]. Since the dominant configuration of the B_3u(nπ^*) state is $( 6 a g ) 1 ( 2 b 3 u ) 1$⁠, irrespective of the spin multiplicity, the dominant configuration of the final core-excited state is supposed to be $( C 1 s ) 7 ( 2 b 3 u ) 1$ (see also Fig. 14).

However, the core excitation changes the effective Coulomb potential that the valence electrons experience so that the canonical orbitals of the initial and final states are not exactly the same. As a result, concomitant valence-electronic excitation of “shake-up”^5,81 type may occur and a contribution to the spectrum of the energetically close configuration $( C 1 s ) 7 ( 1 a u ) 1$ is also expected. Therefore, the true core-to-SOMO transition peaks are supposed to split. The shake-up is a double-electronic excitation from the open-shell reference so that it is not well described within our method using this reference, whereas both the main and shake-up peaks are well described within our method using the conventional ground-state reference because the shake-up configuration is generated by one-electron excitation from the ground state configuration (see schematic Figs. 9 and 14). We note again that, except for the core-to-SOMO transitions, the final states of the core-excitations from valence-excited states in the bleaching region are dominated by doubly excited configurations from the ground state, which are not well described within standard single-reference formalisms, including the fc-CVS-EOM-CCSD method used here. At the same time, from an experimental point of view, the core-to-SOMO peaks are the preferred choice for tracking excited-state dynamics, since they unambiguously occur only in the valence-excited states.

For these reasons, our discussion will hereafter concentrate on the core-to-SOMO transition peaks. We note that our spectrum computed using the high-spin open-shell reference and the TDDFT spectrum using the low-spin open-shell reference of the same orbital occupation are very similar, apart from a total energy shift.

The left column of Fig. 13 shows the carbon K-edge NEXAFS of the ¹B_3u(nπ^*), ¹A_u(nπ^*), and ¹B_2u(ππ^*) excited states at the minima of the XMCQDPT2 potential energy surfaces [Eq. (24)]; the corresponding NTOs are summarized in Fig. 15. The excitation processes from the ¹B_3u(nπ^*) state are schematically illustrated in Fig. 14. The schemes for the corresponding ones from the ¹B_2u(ππ^*) and ¹A_u(nπ^*) states are shown in Fig. 16.

Note that the spectra of the excited states at the EOM-CCSD-optimized geometries and those at the minima of the XMCQDPT2 potential energy surfaces are almost identical. Two peaks appear in the spectra of all three states. For the ¹B_3u(nπ^*) state, the two peaks appear at 282.5 and 284.0 eV with nearly equal intensities. The peaks for the ¹A_u(nπ^*) state appear at 282.9 and 284.7 eV. The first peak is three times more intense than the second one. The peaks for the ¹B_2u(ππ^*) state appear at 282.0 and 283.5 eV. The first peak is of nearly equal intensity as that of the ¹A_u(nπ^*) state. The second peak is three times more intense than the first one, and prominent among the core-to-SOMO excitation peaks. The sum of all $σ K 2$ for each peak given in Fig. 15 is less than 0.30. As mentioned in Subsection II A, this sum equals the square norm of one-particle transition density matrix; therefore, its maximal value is 1, which is achieved when the transition is of purely one-electron excitation character. Thus, these C_1s → SOMO excitations feature significant many-electron excitation character.

The final core-excited states used in the fc-CVS-EOM-CCSD approach to compute the core-to-SOMO transitions are the dark counterparts of the core excitations yielding the NEXAFS spectrum of the ground state. As illustrated in Fig. 8, the core excitations from the ground state have strong one-electron character. Therefore, character of the final core-excited states of the core-to-SOMO transitions may be investigated by the C K-edge NEXAFS spectra for the ground state at the geometries optimized for the ¹A_u(nπ^*) and ¹B_2u(ππ^*) states, which are shown in Fig. 17. [We skipped the calculation of the spectrum at the geometry optimized for the ¹B_3u(nπ^*) state, which is close to the FC geometry, see Table II.] The first two peaks of the spectrum at the energy minimum of ¹B_2u(ππ^*) appear to have and preserve characters similar to those in the spectrum at the FC geometry. In the spectrum obtained at the energy minimum of ¹A_u(nπ^*), however, the oscillator strength corresponding to the second peak in the spectrum at the FC geometry is 500 times less than that corresponding to the first peak; thus, the peak disappears. The value of $σ K 2$ for the peak that disappears is 0.58 for 2a_g → 2b_3u, and 0.16 for 1b_3g → 1a_u, whereas that corresponding to the first peak is 0.55 and 0.15 for 1b_3g → 1a_u and 2a_g → 2b_3u, respectively. These single-configuration characters of the two core-excited states and weak oscillator strength for the second peak may be explained from consideration of the nodal structure of the orbitals (see Fig. 6): As shown in Table II, the C–C bond elongates by 0.1 Å and the C–N bond shortens by 0.03 Å at the energy minimum of the ¹A_u(nπ^*) state. Elongation of the C–C bond leads to a weaker stabilization effect along this bond and lower electron density on the C atoms of the 2b_3u orbital. In addition, the shorter C–N bond leads to energetic destabilization of the 2b_3u orbital. As a result, the 2b_3u orbital becomes separated from the 1a_u orbital toward energetically higher position, and the transition dipole moment between the ground state configuration and $( C 1 s ) 7 ( 2 b 3 u ) 1$ becomes small. Recalling that the dominant configurations of the ¹B_3u(nπ^*) and ¹A_u(nπ^*) states are $( 6 a g ) 1 ( 2 b 3 u ) 1$ and $( 6 a g ) 1 ( 1 a u ) 1$⁠, respectively, this is consistent with the fact that the energetic ordering of these states switches between the potential minima.

As discussed before, the final core-excited states corresponding to the first two peaks of the carbon K-edge NEXAFS for the ground state at the FC geometry maintain strong single-configuration character, possibly owing to a strong screening effect on the Coulomb potential from the core hole. This screening effect by the strongly delocalized 1b_3u(π) electrons is not supposed to be very sensitive to the molecular geometry. Therefore, the single-configuration character of the final core-excited states is also expected to be the case for the final states of core-excitations from the valence-excited electronic states. For this reason, the NEXAFS spectrum for ¹A_u(nπ^*) may accurately show the dominant peak, and the peak corresponding to the shake-up of 1a_u → 2b_3u. On the other hand, for the ¹B_3u(nπ^*) and ¹B_2u(ππ^*) states, mixing of the configurations $( 2 b 1 u ) 1 ( 2 b 3 u ) 1$ and $( 1 b 2 u ) 1 ( 1 a u ) 1$⁠, which are the counterparts of $( 2 a g ) 1 ( 2 b 3 u ) 1$ and $( 1 b 3 g ) 1 ( 1 a u ) 1$⁠, respectively, may be overestimated in the final core-excited states. Thus, the relative intensities and the separation of the two peaks in the NEXAFS spectra for the ¹B_3u(nπ^*) and ¹B_2u(ππ^*) states are not necessarily accurate enough. The protruding character of the second peak in the spectrum for the ¹B_2u(ππ^*) state is especially uncertain.

The right column of Fig. 13 shows the carbon K-edge TR-NEXAFS spectra of the individual electronic states. The position of peaks oscillates in the range of about 0.5 eV, as the nuclear wave packet oscillates on the potential energy surface. The spectra of the ¹A_u(nπ^*) state, which are the most reliable, indicate that the peak intensities are stable against the position of the nuclear wave packet, probably due to the stability of the configurations throughout the nuclear dynamics.

Figure 18 shows the total two-dimensional carbon K-edge TR-NEXAFS spectrum, and its slice cuts at the delay time of 30, 60, 90, and 120 fs. In the slice cuts, the intensities of the core-to-SOMO transitions for the individual electronic states are weighted by the populations, broadened by the Lorentzian functions, and summed together. Then, it is clear that the relative peak positions for the individual electronic states are not fixed throughout the photophysical dynamics. In particular, the peaks from the ¹B_2u(ππ^*) state are separated from those from the ¹A_u(nπ^*) state at the delay time of 30 and 90 fs but are submerged under the dominant peak of the ¹A_u(nπ^*) state at 60 and 120 fs. Due to these different oscillations of the peaks, the contribution of the dark ¹A_u(nπ^*) state is revealed by the TR-NEXAFS spectra at the carbon K-edge. Unfortunately, however, the separation of peaks for the ¹B_3u(nπ^*) and the ¹A_u(nπ^*) states is not sufficient for the extraction of the population dynamics between these two states. It should be noted again that the present qualitative character of the C K-edge TR-NEXAFS spectra is possibly affected by the inaccurate description of the final core-excited states. Therefore, tracking the population dynamics in photoexcited organic molecules with TR-NEXAFS might be hindered by difficulties in computing the final core-excited states and the uncertainties in theoretical protocols. From an experimental point of view, energy resolution of the electron detectors can also be hindrance. Nonetheless, C K-edge TR-NEXAFS spectra should be sensitive to the molecular electronic structure, which depends on the symmetry lowering due to distortions and substitutions. Thus, TR-NEXAFS may be rather convenient for tracking photoinduced isomerization reactions, such as pyrazine → pyrimidine. In this context, both accurate theoretical modeling methods and electron detectors with higher energy resolution are desired.

Figure 19 shows the nitrogen K-edge NEXAFS of the ground electronic state in the solid line (the calculated IE is 407.1 eV), with the one obtained from the experiment as the dashed line, redigitized from Ref. 76. The respective NTOs are given in Fig. 20. The main peak at 400.0 eV corresponds to the 1a_g → 2b_3u transition, with minor contributions from 1b_1u → 2b_2g. The small peaks at 405.1 and 405.8 eV correspond to the transitions to diffuse Rydberg orbitals. In contrast to the carbon K-edge spectrum, the 1a_u orbital is not involved in the N_1s excitations from the ground state because it does not span any nitrogens. Consequently, the electron-hole localization does not occur. As a result, the present calculated spectrum agrees well, apart from some energy shifts, with the experimental one over the entire pre-edge region.⁷⁶ 

The spectrum for the lowest-lying high-spin open-shell reference, which has the same orbital character as the ¹B_3u(nπ^*) state, is also shown in Fig. 21 to estimate the fine structures of the NEXAFS spectrum of the ¹B_3u(nπ^*) state over the entire pre-edge region. Single occupation of the 2b_3u orbital reduces the intensity of the main peak by half. Excitation to the near-ionization region is enhanced by the contribution of the two-electron excited configurations like $( 1 a g ) 1 ( 6 a g ) 1 ( 2 b 3 u ) 1 ( 6 b 1 u ) 1$⁠. The characteristic core-to-SOMO peak appears at 396.0 eV.

The left column of Fig. 22 shows the core-to-SOMO peaks of the NEXAFS spectra of the ¹B_3u(nπ^*), ¹A_u(nπ^*), and ¹B_2u(ππ^*) states at the respective XMCQDPT2 potential energy minima; the corresponding NTOs are summarized in Fig. 23. Dominance of single configurations in the final core-excited states results in relatively strong single-excitation characters compared to the C_1s core excitations. The core excitations from both the ¹B_3u(nπ^*) and ¹A_u(nπ^*) states are dominated by the 1b_1u → 6a_g transition. For the ¹B_2u(ππ^*) state, only the minor $( 1 b 2 g ) 1 ( 1 a u ) 1$ configuration allows the core-to-SOMO excitation of 1b_1u → 1b_2g and the intensity of the peak is low compared to those for the ¹B_3u(nπ^*) and ¹A_u(nπ^*) states. The dominant configurations for the final core-excited states are $( 1 b 1 u ) 1 ( 2 b 3 u ) 1$ for excitation from the ¹B_3u(nπ^*) state and $( 1 b 1 u ) 1 ( 1 a u )$ for excitation from both the ¹A_u(nπ^*) and ¹B_2u(ππ^*) states. In the $( 1 b 1 u ) 1 ( 1 a u ) 1$ configuration, the electron is always separated from the N_1s core hole, whereas the hole and the particle can be on a same nitrogen atom in the $( 1 b 1 u ) 1 ( 2 b 3 u ) 1$ configuration, maximizing the Coulomb interaction between them. This energetic disadvantage of the $( 1 b 2 g ) 1 ( 1 a u ) 1$ configuration relative to $( 1 b 1 u ) 1 ( 2 b 3 u ) 1$ leads to a blue shift of 2.4 and 2.6 eV, respectively, in the positions of the peaks for the ¹A_u(nπ^*) and ¹B_2u(ππ^*) states relative to the ¹B_3u(nπ^*) state (i.e., 398.2 and 398.4 eV vs 395.8 eV). This energetic separation between $( 1 b 1 u ) 1 ( 2 b 3 u ) 1$ and $( 1 b 1 u ) 1 ( 1 a u ) 1$ hinders not only the mixture of the configuration in the final core-excited states but also the shake-up effect. This is the reason for almost complete agreement of the core-to-SOMO transition peak between the spectra for the triplet and singlet B_3u(nπ^*) states [compare Figs. 21 and 22(a)].

The right column of Fig. 22 shows the nitrogen K-edge TR-NEXAFS spectra for the individual electronic states. Due to the relatively strong single-excitation characters, both the intensities and positions are stable throughout the nuclear dynamics. For the ¹B_2u(ππ^*) state, the amplitude of the minor configuration $( 1 b 2 g ) 1 ( 1 a u ) 1$ depends on the molecular geometry so that the peak intensity varies in the course of the nuclear dynamics.

Figure 24 shows the total nitrogen K-edge TR-NEXAFS spectra. Due to insensitivity to the molecular distortion around the FC geometry, the spectra purely reflect the population dynamics. The peaks for the ¹B_2u(ππ^*) and ¹A_u(nπ^*) states overlap almost completely so that the populations of these two states cannot be individually resolved. If the developing process of the peak near 398.0 eV could be observed, one might recognize that a nuclear wave packet is created in another electronic state after pyrazine has been promoted to the ¹B_2u(ππ^*) state. However, in an experiment, the X-ray pulse shots for the probe process will be initiated after the pump pulse has been switched off. Since the ¹B_2u(ππ^*) state has almost completely decayed at around 40 fs (see Fig. 2), tracking excited-state dynamics begins only after the decay process. If the experimental analysis is implemented taking only the optically bright states into account, then it would lead to the incorrect conclusion that the nuclear wave packet oscillates between the ¹B_2u(ππ^*) and ¹B_3u(nπ^*) states. In this context, even simple TR-NEXAFS spectra can be deceptive. Careful consideration of all the electronic states energetically close to or below the state where the molecule is initially promoted, including optically dark states, and a rigorous analysis of the correspondence between the electronic structures and the fine structures of the spectra are thus required.

On the basis of the above results, heteroatom K-edge TR-NEXAFS appears to be particularly effective for tracking the excited-state dynamics involving dark states in the FC region. The present theoretical methods afford the extraction of the population dynamics so that visualization of excited-state dynamics in the FC region is already feasible.

We computed and analyzed the TR-NEXAFS spectra of photoexcited pyrazine to assess the effectiveness of these spectroscopies for tracking excited-state dynamics. We simulated nuclear dynamics using the MCTDH method. The nuclear dynamics simulation upon excitation to the ¹B_2u(ππ^*) state shows contributions of the direct decay to the ¹B_3u(nπ^*) state and from another pathway via the optically dark ¹A_u(nπ^*) state. The simulations also show the oscillations between the potential minima of the ¹B_3u(nπ^*) and ¹A_u(nπ^*) states driven by the totally symmetric vibrations after the decay. The analysis of the electronic structure and the NEXAFS spectra revealed different characters of C and N (heteroatom) K-edge TR-NEXAFS. On one hand, the final core-excited states for the nitrogen K-edge TR-NEXAFS are dominated by single electronic configurations close to the initial valence-excited states. As a result, the spectra are insensitive to molecular distortions near the FC region and reflect the population dynamics of the involved electronic states. On the other hand, for the carbon K-edge TR-NEXAFS, there are not only more symmetry-allowed final core-excited states, but several electronic configurations may also dominantly contribute to one of the final core-excited states. These electronic states are sensitive to even small changes of molecular symmetry and geometry. As a result, these factors complicate the TR-NEXAFS spectra. Therefore, TR-NEXAFS employing the absorption edge of the heteroatom is a better choice for tracking the photophysical process near the FC region, whereas employing the K-edge of the carbon atoms would be better to track photoinduced isomerization. The other significant difference between the C and N K-edge TR-NEXAFS is that separation of the core-to-SOMO excitation peaks for the ¹A_u(nπ^*) and ¹B_2u(ππ^*) states is impossible in the N K-edge spectra, but possible in the C K-edge spectra, and vice versa for the separation of those for the ¹B_3u(nπ^*) and ¹A_u(nπ^*) states. Thus, in the case of pyrazine, it is necessary to employ both the C and N K-edge TR-NEXAFS for complete extraction of the population dynamics. More generally, combination of TR-NEXAFS of different absorption edges is effective for a comprehensive understanding of the photoinduced processes.

When it comes to extracting photophysical dynamics from the experimental TR-NEXAFS spectra, this is already feasible in the FC region: The population dynamics is obtained by selection of the involved electronic states and assignment of the computational NEXAFS spectra at a single geometry to the TR-NEXAFS spectra, only employing the core-to-SOMO excitation peaks. Then, from the population dynamics as extracted from the experimental spectra, the (correct) vibrational modes and potential energy surfaces are obtained by trial and error. Finally, the nuclear wave packet evolution processes are visualized. For photoinduced isomerization, however, the limitations of the theoretical protocols may hinder the extraction of the dynamics. Therefore, further development of theoretical methods is desired for aiding future experiments.

The current theoretical protocols, which only sample the centroids of the nuclear wave packets, may face a limitation when the nuclear wave packets are broad. In such cases, accuracy would be improved if multicentroids are sampled. When an interstate coupling is strong, it may result in a bifurcation on a nuclear wave packet along the coupling normal coordinate. In this case, molecular geometries corresponding to the expectation values of the absolute values of the normal coordinates, |Q|, should be sampled instead. Note again that adiabatic and diabatic electronic states are the same if and only if the molecule maintains the highest symmetry when sampling the centroids of the nuclear wave packets. Therefore, in the case of strong interstate coupling, one should either adiabatize the nuclear wave packets or compute the oscillator strengths for transitions between diabatic states. These limitations are discussed in  Appendixes A and  B.

In contrast to other theoretical studies of TR-NEXAFS, the present work focuses on the analysis of how the electronic structure and the nuclear dynamics manifest themselves in the spectra. We hope that this study can be instructive on how to employ TR-NEXAFS especially to investigate nuclear dynamics, and ultimately chemical reactions.

We thank Dr. Michael Epshtein and Professor Stephen R. Leone at the University of California for discussions from an experimental perspective. We also thank Dr. Thomas Northey and Dr. Thomas J. Penfold at Newcastle University for valuable discussions about electronic structure aspects and Dr. Rasmus Faber at DTU Chemistry for his help with some computational tasks. This research was funded by the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie Grant Agreement No. 713683 (COFUNDfellowsDTU), by DTU Chemistry, and by the Danish Council for Independent Research (now Independent Research Fund Denmark), Grant Nos. 4002-00272, 014-00258B, and 8021-00347B. A.I.K. acknowledges research support from the U.S. National Science Foundation (Grant No. CHE-1856342) and the Simons Foundation, which funded her sabbatical stay in Germany.

Conflicts of interest. The authors declare the following competing financial interest: A.I.K. is a part owner and a board member of Q-Chem, Inc.

We calculated the NEXAFS spectrum of ¹A_u(nπ^*) at 60 fs by replacing the nuclear wave packet of Eq. (29) with the reduced nuclear density along Q_6a (the normal coordinate of the dominant totally symmetric mode),

where $A u *$ means the final core excited states from ¹A_u(nπ^*) and $Q ¯ ′ ( A u ) ( t )$ is the set of expectation values of the other normal coordinates. The result is shown as an impulse plot in Fig. 25, superposed by the spectrum calculated with Eq. (30) and broadened by Lorentzian (FWHM = 0.4 eV).

The reduced nuclear density along Q_6a at 60 fs is asymmetric with respect to the centroid, as shown in Fig. 4. Reflecting this shape, the two peaks of the NEXAFS spectrum calculated with Eq. (A1) are also asymmetric. In addition, one can also see lump features corresponding to the tail of the reduced nuclear density. Importantly, the peak positions agree with those given by Eq. (30) within 0.2 eV.

On the basis of the above, we confidently conclude that the positions and relative intensities of the dominant absorption peaks are reproduced with reasonable accuracy by single-point calculations at the centroids of the nuclear wave packets, as per Eq. (30), and the vibrational broadening is reproduced by a Lorentzian broadening of 0.4 eV. However, the shape of the absorption peaks reflects the asymmetric features of the nuclear wave packets. Furthermore, if the nuclear wave packet bifurcates, the absorption peaks would also develop split features. In these cases, the absorption spectra calculated with Eq. (30) and a Lorentzian broadening may not be quantitatively accurate.

Although Eq. (29) can be easily used in combination with lower-level electronic structure methods, such as TDDFT, the peak positions and intensities would be affected by poor accuracy of the underlying methods. Unfortunately, evaluating the integral of the oscillator strength calculated with more advanced many-body methods over the entire potential energy surface can quickly become unpractical. One possible way to overcome this dilemma is to represent the nuclear wave packets as superpositions of symmetric functions such as Gaussians and to calculate the oscillator strength at the centroids of each symmetric function. In this context, AIMS⁴⁰ or vMCG^44,45 is rather suitable for TR-NEXAFS, provided that the nuclear dynamics simulation is sufficiently accurate.

In the present model vibronic Hamiltonian (see Sec. II B), the interstate coupling is zero at Q_κ = 0, i.e., at the FC geometry, see Eq. (25). Thus, the coupling is linear with Q_κ for all the nontotally symmetric modes. Consequently, the transitions between the relevant electronic states do not occur at Q_κ = 0, but rather at a finite value of Q_κ. Therefore, even though the expectation values of these normal coordinates remain zero, the nuclear wave packet may split along these normal coordinates, as shown in Fig. 26, especially when the coupling constant in Eq. (25) is large, such as λ_10a between ¹B_3u(nπ^*) and ¹B_2u(ππ^*). In addition, molecular distortions along the positive and negative directions of normal coordinates of the nontotally symmetric modes are symmetric with respect to inversion. To retain a one-point geometry-sampling picture at geometries of lower symmetry than D_2h, one could instead sample the absorption spectra at nuclear configurations corresponding to the expectation value of |Q|. These are shown in Fig. 26 as vertical dashed lines.

Figure 27 illustrates the characters of the transitions between the adiabatic states, S₀ → S₁, S₂, at Q_10a = 1 and 4 (all other normal coordinates are fixed at zero). Comparing with Fig. 6, one finds that the nπ^* and ππ^* characters begin to mix in the region where the nuclear wave packets spread. The NEXAFS spectra of S₁ and S₂ of at Q_10a = 4, shown in Fig. 28, reflect such character mixing.

On the basis of these observations, we conclude that the mixing of electronic characters cannot be entirely dismissed throughout the simulations of photoinduced processes. Developing methods for separate extraction of the electron configurations from the TR-NEXAFS spectra remain a task for future work. As illustrated in Sec. IV B 3, the spectra at the absorption edges of heteroatoms are less sensitive to molecular geometries, more clearly reflecting the electronic characters. Therefore, the direct extraction of population dynamics of the diabatic electronic states appears to be possible from the spectra recorded at the absorption edges of heteroatoms.