The simulation of X-ray absorption spectra including vibronic coupling: application of the QD-DFT/MRCI(2) method

Abstract

We here demonstrate the accuracy and computational efficiency of the recently introduced QD-DFT/MRCI(2) formalism [Neville et al., J. Chem. Phys., 2024, 160, 234109.] for the direct calculation of quasi-diabatic core-excited electronic states via the simulation of the X-ray absorption spectra for a trio of unsaturated hydrocarbons: ethylene, allene, and butadiene. This is achieved via the construction of vibronic coupling Hamiltonians, including anharmonicity through 6th-order in the one-mode terms, as well as bilinear two-mode coupling terms. The models include all vibronically-coupled core-excited states below the core-ionization potential: 24, 27, and 26 states of the C K-edge for ethylene, allene, and butadiene, respectively. These results, which furnish excellent agreement with recently measured X-ray absorption spectra, demonstrate the effectiveness of the QD-DFT/MRCI(2) method. Further, they highlight the importance of accounting for vibronic coupling in first principles simulation of X-ray absorption spectra. Specifically, we demonstrate that simulations constructed from best estimate vertical excitation energies and oscillator strengths alone generally fail to predict the positions of experimental peak maxima, particularly for transition energies above the band origin.

---

1 Introduction

X-ray absorption spectroscopy (XAS) is a powerful technique that is sensitive to the nuclear and electronic structure of molecules.^1–3 Since the core electrons probed by this approach are highly localized in space, XAS can offer not only element-specific insights, but also site-specific information regarding local chemical environments. This makes core-level spectroscopies such as XAS ideally suited for tracking the movement of electronic charge that results from photo-excitation.^4–7 Recent developments in X-ray laser sources^8–10 have introduced new possibilities for ultrafast time-resolved pump–probe experiments^11–14 able to probe vibronic wavepacket dynamics in molecules.^15–18

In a manner completely analogous to long-established valence-level pump–probe experiments, an accurate accounting of both the manifold of vibronic states into which the initial state is pumped as well as those onto which the wave packet is projected by the probe is essential for the correct interpretation of the time-evolving signal. Unlike valence-level spectroscopy, however, which generally involves transitions from a few high-lying frontier orbitals to a small number of unoccupied valence orbitals, the pre-edge structure of X-ray spectra are characterized by transitions from the core orbitals of a specific element (C 1s, S 2p, etc.) on each atomic center in a molecule, resulting in dense manifolds of electronic states. The presence of these dense manifolds all but ensures that strong vibronic coupling will manifest. This means that the spectral envelope will, in general, be comprised of strongly vibronically mixed states, the analysis of which typically requires advanced theoretical and computational approaches to interpret and assign. Furthermore, given that each core-level transition is a resonance embedded in a continuum far above the lowest ionization potential, the lifetime of these core-excited states will be on the order of femtoseconds. The spectral envelope will thus be a result of pronounced lifetime broadening of the vibronic energy levels, in addition to non-Born–Oppenheimer coupling effects.^19–21

Typically, these spectra are simulated using quantum chemistry methods via the computation of vertical excitation energies and the corresponding oscillator strengths at the ground-state minimum geometry. The spectral envelope is then constructed via the convolution of the resultant purely electronic stick spectrum with a Gaussian or Lorentzian function with an empirically chosen full-width at half maximum (FWHM).^22,23 More sophisticated simulations may include vibronic effects within the Born–Oppenheimer approximation (BOA) via the consideration of vibrational structure in the core-excited states.^24–29 However, the high computational cost of calculating the highly dense manifolds of states required for X-ray spectroscopy simulations, coupled with the inclusion of vibrational degrees of freedom, presents a significant challenge and limits these methods to the treatment of small molecules. Furthermore, the approximation that the core-excited states are well-described within the Born–Oppenheimer approximation will, in general, not hold.

The first study to consider non-Born–Oppenheimer vibronic coupling in the simulation of XAS using techniques related to those discussed here was reported by Gadea et al. in their 1991 work on the C K-edge spectrum of ethylene.^20,30–32 The authors showed that the C K-edge origin band was a dense manifold of vibronic states arising from coupling between the two 1sπ* electronic states. Those simulations involved the construction of a two-state, six-mode model and included the nearly degenerate ¹B_1u and ¹B_2g states in conjunction with each of the three a_g modes, both b_3u-symmetry coupling modes, and the a_u-symmetry torsional mode. The vibronic wave functions and state energies were determined using a time-independent approach in which a linear vibronic Hamiltonian (LVC)^30–32 augmented with selected anharmonic terms involving the torsional mode was used to simulate the XAS using a time-independent methodology. The authors identified the importance of the coupling modes of b_3u symmetry on the overall spectral envelope of the origin band of the C K-edge XAS of ethylene. Moreover, a complementary investigation extended the theoretical description of vibronic coupling in ethylene by incorporating additional electronic states and vibrational modes and applying the model to the C K-edge XAS of a number of deuterated isotopomers.^33,34

Subsequent studies have employed LVC model Hamiltonians to explore the role of vibronic coupling effects in the XAS of acetylene and related isotopomers.³⁵ To the best of our knowledge, simulations involving vibronic coupling in K-shell absorption spectra have so far been limited to small molecules (involving a very small number of chemically equivalent cores), typically involving a reduced-dimensionality vibrational subspace, and a focus on the band origins,^20,33–39 highlighting the computational challenges of extending such approaches to larger systems.

The determination of expansion coefficients for vibronic coupling models, essential for exploring vibronic effects in XAS, is a nontrivial task. Typically, these coefficients are obtained from quantum chemistry calculations, which are performed in the adiabatic electronic basis, and then transformed to the required diabatic basis via an adiabatic-to-diabatic transformation. A more recent approach is QD-DFT/MRCI(2) method,^40,41 which constructs diabatic potentials and couplings directly using a perturbative approximation to block-diagonalization diabatization^42,43via an effective Hamiltonian formalism.^44–48 This diabatization method utilizes the machinery of DFT/MRCI(2)^49–51 and can thus describe a large range of electronic states, including those of multi-reference, Rydberg, doubly excited and charge transfer character. Given that it utilizes DFT/MRCI(2) wave functions, it is largely black-box, computationally efficient and generally furnishes predictive accuracy. Furthermore, and germane to the present work, the direct construction of diabatic potentials and couplings means that higher-order potential terms are readily determined. The direct computations of diabatic states, coupled with the DFT/MRCI(2) formalism also enables the construction of vibronic models with a very large number of electronic states. This will be of particular import in the simulation of core-excitation spectra which are generally comprised of dense manifolds of electronic states for which standard diabatization techniques may become impractical.

In this work, we employ the QD-DFT/MRCI(2) method to directly compute quasi-diabatic potential matrices for core-excited states, including the vibronic coupling between them. The direct computation of these quantities enables the construction of vibronic coupling Hamiltonians involving large numbers of electronic states. The corresponding absorption spectra are then determined within a time-dependent framework via wave packet dynamics simulations employing the multi-layer multiconfigurational time-dependent Hartree (ML-MCTDH) approach.^31,32 To demonstrate the efficacy of the approach, we here simulate the pre-edge structure of the X-ray absorption spectra of ethylene, allene, and butadiene, for which C K-edge X-ray absorption spectra have recently been reported.⁵²

Furthermore, a secondary aim of the present work will involve the comparison of spectral simulations obtained from a hierarchy of approximations, moving from purely electronic spectra (employing only vertical electronic excitation energies and oscillator strengths) to Born–Oppenheimer harmonic vibrational models, to full vibronic coupling models. In addition to demonstrating the computational efficacy of our approach, the ability to generate benchmark-level spectral simulations will enable a detailed analysis of the relative importance of the various factors that give rise to the overall spectral envelope. In particular, the present results will enable comment on the extent to which the line shapes of the main absorption features in XAS spectra are influenced by the significant vibronic coupling of the underlying core-excited states versus the lifetime broadening of these transitions, all of which are significantly above the ionization potential.

The rest of the paper is arranged as follows. In Section 2, we will provide a brief overview of the computational methods employed to simulate XAS, including the electronic structure calculations, the construction of the vibronic coupling Hamiltonians, and the ultimate simulation of the absorption spectra using quantum dynamics calculations. In Section 3, we demonstrate the accuracy and applicability of the QD-DFT/MRCI(2) method by simulating the absorption spectrum of ethylene, allene, and butadiene. In Section 4, we will provide an analysis and assignment of the XAS, as well as compare the results obtained using different levels of approximation. Finally, in Section 5, we will provide our concluding remarks on the advantages of the QD-DFT/MRCI(2) method for calculating diabatic potentials and couplings, where we emphasize the significance of vibronic coupling in shaping the XAS spectral envelope and the limitations of simulations employing only vertical excitation energies.

2 Methods

2.1 Calculation of vibronic absorption spectra

2.1.1 Vertical excitation approximation. In this work, XAS are simulated employing a hierarchy of approximations. The simplest, and most widely employed, is to approximate the absorption spectrum using electronic vertical excitation energies and oscillator strengths only. In practice, in order to phenomenologically account for lifetime and vibronic broadening effects, the resulting stick spectrum is convoluted with a Lorentzian function with a FWHM chosen empirically so as to furnish the best agreement with experiment. We will denote the resulting XAS as I_VE(E), where VE denotes “vertical excitation”. This quantity may be written as follows:where N is the number of core-excited states included, ε_i are the vertical excitation energies, f_i are the oscillator strengths, L(E; Γ) is a Lorentzian function of FWHM of Γ, and * denotes the convolution operation. Here, the parametric dependence of I_VE(E) on Γ is to be implicitly understood. A significant broadening mechanism will generally be lifetime broadening due to Auger decay of the core-excited states. For the molecules in this study, the C K-edge transition lifetimes are around 10 fs or less, which corresponds to a transform-limited peak width of around 0.1 eV.⁵² However, the values of Γ employed in the following simulations are significantly larger than this value, as this convolution accounts for additional effects (i.e., vibronic coupling) that are not included in the vertical excitation approximation. In the vertical excitation simulations, Γ was set to 0.3 eV for each molecule.

2.1.2 Harmonic vibrational approximation. The next level in the hierarchy of approximation incorporates vibrational structure into the XAS within the Born–Oppenheimer and harmonic approximations. The core-excited state manifold is assumed to be comprised of a set of displaced and energetically shifted ground state harmonic potentials, where the magnitude of the energetic shift is determined by the vertical excitation energy and the geometric displacement is wholly determined by the gradient on the core-excited state computed at the ground state minimum geometry. The energy levels of this Hamiltonian have an analytic solution and generates the so-called Poisson distribution spectrum.³¹ We denote the resulting absorption spectrum as I_HV(E), with the HV standing for “harmonic vibrational”. At this level of approximation, the absorption spectrum assumes an analytical form, and, after convolution by the line shape L(E; Γ), is given bywhere α indexes the ground state mass- and frequency-scaled vibrational normal modes, with corresponding frequencies ω_α, and M_α denotes the maximum excitation degree of each normal mode. is the Franck–Condon factor corresponding to the overlap of the ground vibrational state in the ground electronic state with the vibrational state with n_α quanta in each of the modes α in the core-excited electronic state i, and d denotes the total number of vibrational normal modes included in the model. Within the harmonic vibrational approximation, these are given bywhere A_iα denotes the so-called Poisson parameter,with g⁽ⁱ⁾_α being the gradient of the i'th core-excited state with respect to the normal mode α evaluated at the Franck–Condon point.

In these simulations, the Lorentzian line shape L(E; Γ) was chosen to have a FWHM Γ that accounts primarily for lifetime broadening. The reason for this is so that, when compared to the results of the simulations performed using the vibronic coupling approximation models (see Section 2.1.3), the effects of vibronic coupling may be more readily discerned. Thus, the values of Γ chosen reflect the expected core–hole life times of the states involved, with some small adjustments being made to bring the results into agreement with the experimental results. Specifically, we assumed core–hole lifetimes of 5.0 fs for the origin band, and 7.5 fs for all other bands in ethylene, 7.5 fs for all transitions in allene, and 15 fs for all transitions in butadiene. These values correspond to FWHMS of 0.26 eV, 0.18 eV, and 0.09 eV, respectively.

2.1.3 Non-adiabatic vibronic coupling approximation. The final, and most sophisticated, level of approximation used corresponds to the use of a vibronic coupling Hamiltonian, going beyond the BOA and accounting for the non-adiabatic coupling between the core-excited states. The quasi-diabatic states employed as an electronic basis are denoted {|ϕ_i(Q)〉}, where the parametric dependence on the nuclear coordinates Q is explicitly given. In a quasi-diabatic representation, neglecting the small residual derivative couplings, the non-adiabatic coupling is contained within the now non-diagonal potential matrix, with the total, vibronic Hamiltonian given bywhere Ĥ_el denotes the electronic Hamiltonian, and T̂_n the nuclear kinetic energy operator.

The quasi-diabatic potential matrix elements W_ij(Q) = 〈ϕ_i(Q)|Ĥ_el|ϕ_j(Q)〉 are approximated using the vibronic coupling Hamiltonian model of Köppel et al.^31,32,53 This corresponds to the truncated Taylor expansion of the quasi-diabatic potential matrix W(Q) in terms of the ground state mass- and frequency-scaled normal modes Q_α about the Franck–Condon point Q₀. In this work, we consider an expansion in which only one- and two-mode terms are retained, with the one-mode terms appearing through 6th-order, and the two-mode terms through 2nd-order:

The global gauge of the adiabatic-to-diabatic transformation was fixed by taking the adiabatic and quasi-diabatic states to be equal at the point of expansion Q₀, giving rise to the relation τ^(i,j)₀ = ε_iδ_ij for the zeroth-order expansion coefficients, where ε_i is the vertical excitation energy of the ith adiabatic core-excited state. The remaining expansion coefficients were determined using a least-square fit to quasi-diabatic potential matrix values computed directly using the QD-DFT/MRCI(2) method, using the fitting procedure described in ref. 54.

The absorption spectrum obtained using the vibronic coupling Hamiltonian will be denoted by I_NA(E), where NA stands for “non-adiabatic” (see Table S1). This could be computed in a manner analogous to the electronic vertical excitation spectrum (eqn (1)), with the substitution of the vibronic eigenvalues and transition moments for their electronic counterparts. However, as the complexity of this problem scaled exponentially with system size, and for vibronic coupling systems the density of states is usually so high that only the overall band envelope is particularly meaningful, we choose to follow a time-dependent approach. Specifically, the vibronic absorption spectrum was computed via the Fourier transforms of wave packet autocorrelation functions corresponding to initial wave packets formed from excitation of the vibronic ground state to each of the optically bright core-excited electronic states:

where the wave packet autocorrelation function a_i(t) is given by

a_i(t) = 〈Ψ₀|Ô_ie^−iĤtÔ_i|Ψ₀〉, (8)

with |Ψ₀〉 denoting the ground state of the vibronic Hamiltonian Ĥ, and

Ô_i = |ϕ_i〉〈ϕ₀| + h.c. (9)

being the vertical excitation operator to core-excited state |ϕ_i〉 from the ground state |ϕ₀〉. We note that in this approach, we have assumed the validity of the Condon approximation, which is taken to be justified due to the slowly varying nature of the quasi-diabatic states with the nuclear coordinates.

The required wave packet propagations were performed using the powerful ML-MCTDH method^55–61 as implemented within the Heidelberg MCTDH package.⁶² For each molecule, all internal nuclear degrees of freedom were included, corresponding to 12, 15, and 24 vibrational modes for ethylene, allene, and butadiene, respectively. A propagation time of 100 fs was used in each case, allowing access to the autocorrelation function up to 200 fs.⁶³

In practice, a finite wave packet propagation time T has to be used. This introduces artifactual high frequency components to the Fourier transforms of the wave packet autocorrelation functions (the so-called Gibbs phenomenon).⁶⁴ To minimize these nonphysical contributions to the spectrum, the autocorrelation function is pre-multiplied by a window function

prior to being Fourier transformed. Here, Θ(t − T) denotes the Heaviside step function centered at time t = T.

Lastly, in order to account for phenomenological broadening present in the experimental spectra, the autocorrelation functions are further scaled by an exponential damping function,

where τ is a chosen damping time. The effect of multiplication of the autocorrelation function by f(t) is to convolute the spectrum with a Loretzian line shape with a FWHM of Γ = 2/τ.^65,66 Here, due to the accounting of vibronic coupling effects, this convolution will predominantly be accounting for lifetime broadening of the core-excited states.^19–21 For ethylenic molecules, the core–hole lifetime is typically around 10 fs.⁵²

2.2 Computational details

The ground state geometries Q₀, normal modes Q_α and frequencies ω_α were computed at the aug-cc-pVTZ/B3LYP^67–70 level of theory using the Turbomole set of programs.⁷¹ The self-consistent field (SCF) procedure was converged with an energy threshold of 10⁻⁶ a.u. Geometries were optimized with a Cartesian gradient threshold of 10⁻⁶ and a quadrature grid value of “m4” in the Turbomole package.⁷¹ Numerical quadrature for the DFT integration was performed with level = 5 in PySCF. The results of these calculations are provided in the SI.

To compute the core-excited state quasi-diabatic potential matrices W(Q) used in the fitting of the vibronic coupling Hamiltonians, the QD-DFT/MRCI(2) method was used within the core-valence separation (CVS) approximation.^72–75 In the QD-DFT/MRCI(2) calculations, the CVS-QE12 Hamiltonian⁷⁶ was used for the core-excited states, while the QE8⁷⁷ Hamiltonian was used for the initial (ground) state. In all calculations, the aug-cc-pCVTZ basis set was used.^68,69 The density fitting approximation^78,79 was used in the evaluation of the two-electron integrals, along with the J&K aug-cc-pVTZ auxiliary basis set.⁸⁰ The number of core-excited states included for ethylene, allene, and butadiene was 24, 27, and 26, respectively, so chosen such that the highest energy state approaches the vertical core-ionization potential. All QD-DFT/MRCI(2) calculations were performed using the General Reference Configuration Interaction (GRaCI) program,⁸¹ while the prerequisite DFT calculations were performed using an interface to the PySCF package.^82–84

In the calculation of the Poisson distribution spectrum, eqn (2), the energy gradients g⁽ⁱ⁾_α at the Franck–Condon point Q₀ are required. For technical reasons, the calculation of analytical gradients is difficult using DFT/MRCI methodologies.⁵¹ Instead, these values are computed from the fitted vibronic coupling Hamiltonian potentials (eqn (6)). As we fix the global gauge of the quasi-diabatic representation by taking the adiabatic and quasi-diabatic representations to be equal Q₀, we simply have that g⁽ⁱ⁾_α = τ^(i,i)_1α. For the harmonic vibrational approximation calculations, a maximum excitation degree M_α = 5 was used for all totally symmetric modes, with the non-totally symmetric modes not contributing given elementary symmetry considerations.

3 Results

Fig. 1 shows the nuclear geometries of the ground state minima for ethylene, allene, and butadiene. Note that the atomic numbering employed in Fig. 1 will be assumed in the subsequent discussion.

Fig. 1 Nuclear structures of the ground state minima for (a) ethylene, (b) allene, and (c) butadiene.

To aid comparison with the experiment, the computed spectra were shifted such that the spectrum computed at the highest level of approximation, I_NA(E), was brought into maximum agreement with its experimental counterpart. Further, this allows us to provide theoretical best estimates of the vertical core-excitation energies involved. The applied shifts were in the range of 0.15 to 0.5 eV, depending on the molecule and individual state, and each shift is given in Tables S2, S4, and S6 in the SI. The maximum shift of 0.5 eV demonstrates the accuracy of the vertical excitation energies yielded by the combination of the DFT/MRCI(2) method with the CVS-QE12 Hamiltonian.⁷⁷

3.1 Ethylene

Fig. 2 shows the X-ray absorption spectrum of ethylene calculated using the hierarchy of approximations described above: vertical excitation energies and oscillator strengths only, I_VE(E) (green line), inclusion of Born–Oppenheimer harmonic vibrational states, I_HV(E) (blue line), and the non-adiabatic full vibronic coupling model, I_NA(E) (red line). The I_VE(E) spectrum in Fig. 2 is annotated with a series of electronic transition assignments (A–F), which are summarized in Table 1, along with the dominant hole-particle natural transition orbital (NTO) pairs for the optically allowed transitions.⁸⁵

Fig. 2 The experimental X-ray absorption of ethylene compared to simulations performed employing a hierarchy of approximations. From bottom, the I_VE is a purely electronic spectrum employing vertical excitation energies and oscillator strengths only (green line), I_HV is Poisson distribution spectrum corresponding to vibrational structure arising from treating the excited states as a series of displaced harmonic oscillators (blue line), and I_NA, the full non-adiabatic vibronic coupling simulation (red line). The experimental spectrum is at top (black line).⁵² The bottom panels evince a detailed comparison of the various approximations with the experimental band origin.

Table 1 Symmetry (Sym), theoretical best estimate vertical excitation energies (ε_i), oscillator strengths f_i, and electronic characters and dominant natural transition orbital (NTO) pair for core-excited states of ethylene. All energies calculated at the ground-state minimum (FC point) are given in eV

Sym	Transition	ε i (eV)	f i	Assignment	Hole	Particle
B1u	A	284.86	0.098	1sπ*
B3u	B	287.29	0.009	1s3s
B2u	C	287.92	0.028	1s3pz
B1u	D	288.58	0.010	1s3px
B2u	E	288.75	0.005	1s3dxy
B2u	F	289.74	0.050	1s3dxz

---

At the D_2h Franck–Condon point, ethylene possesses two symmetry equivalent C 1s atomic orbitals (AOs), the symmetric and anti-symmetric linear combinations of which give rise to two near-degenerate C 1s molecular orbitals (MOs) of a_g and b_3u symmetry, respectively. Excitation from these core-level MOs into the b_2g π* orbital gives rise to the B_1u(1sπ*) and B_2g(1sπ*) states, which lie under the first band of ethylene's X-ray absorption spectrum, centered at around 285 eV. Only the B_1u(1sπ*) state is optically bright. However, strong vibronic coupling to the B_2g(1sπ*) state is known to give rise to vibronic states of mixed electronic character.^20,34 The broad peak at higher energies, in the 287–289 eV region, originates from transitions into the n = 3 Rydberg manifold. As far as we are aware, the role of vibronic coupling in this region of the spectrum has not been studied before.

As mentioned above, the bright band origin is identified as arising due to a transition from the b_3u 1s orbital into the b_2g π* orbital, while the corresponding transition from the a_1g 1s orbital is optically forbidden. While the overall features of the experimental XAS are well reproduced by the vertical excitation energies, the origin band is poorly described by a purely electronic treatment, given the prominent vibrational structure visible in the experimental spectrum.

Describing each core-excited state as a displaced harmonic oscillator yields the I_HV (blue) spectrum in Fig. 2. This level of approximation furnishes much improved agreement with the experimental result. Firstly, it reproduces the shoulder of the origin band, which corresponds to a vibrational progression in the symmetric C–H stretching mode, indicating that this spectral feature arises primarily to vibrational structure in the optically bright state. Previous experimental X-ray absorption studies on ethylene by Ma et al.⁸⁶ and Rabus et al.⁸⁷ interpreted the vibrational structure observed in the C 1s → π* photo-absorption origin band within the BOA. In ethylene, they identified contributions primarily from the C=C and C–H stretching modes, as well as the out-of-plane bending mode. Specifically, they attributed the small shoulder near the band origin to excitation of the C–H stretching mode, which is confirmed by the present simulation. Lastly, the inclusion of the vibrational state contributions to the spectrum also improves the intensity corresponding to the electronic transition F, bringing it into better agreement with the experimental result.

Finally, the full non-adiabatic vibronic coupling model spectral simulation shown in Fig. 2, denoted I_NA and given by the red curve, yields a spectrum very similar to the I_HV simulation and in near quantitative agreement with the experimental result. In contrast to previous work, the current model incorporates 24 electronic states, sufficient to reach the core-ionization potential, along with all 12 vibrational degrees of freedom. The effects of adding non-adiabatic effects and inclusion of anharmonicity via the inclusion of higher-order interstate coupling coefficients reveals a slight broadening and subtle changes in the intensities of the origin band, centered at 285.1 eV, as shown in Fig. S3 in the SI, consistent with the findings of Gadea et al.²⁰ and Köppel et al.³⁴ As shown in Fig. S4 in the SI, the comparison of a model including only anharmonic terms, with no inter-state coupling (orange line), with a model that includes both inter-state coupling and anharmonicity (cyan line), indicates that the spectral broadening arises primarily from the addition of non-adiabatic effects rather than anharmonicity. Furthermore, comparing the cyan line with the full non-adiabatic spectrum, which also includes bilinear coefficients, indicates that the inclusion of bilinear terms has only a minor impact on the spectrum, suggesting that this effect is minimal for ethylene.

As discussed in Section 2.1, the simulation of absorption spectra using a time-dependent approach will generally involve the application of a damping function to the time auto-correlation function to account for peak broadening mechanisms not explicitly included in the Hamiltonian (eqn (11)). In order to achieve optimal agreement with the experimental spectrum, our simulations employ a damping time of 5.0 fs for the first band and 7.5 fs for all other transitions in ethylene, which correlates well with the estimated core–hole lifetime of 10 fs for band origin.⁵² That the empirical damping time that results in the best agreement with experiment aligns well with the core–hole lifetime not only suggests that the broadening due to vibronic coupling is well described by our model, but also unsurprisingly confirms that lifetime-broadening due to Auger decay is the primary broadening mechanism.

Finally, we note that, based on the shifts applied to bring the simulated and experimental spectra into maximum agreement, our theoretical best estimate for the vertical excitation energy of the B_1u(1sπ*) is 284.86 eV. We note that this value is in good agreement with recent benchmark-level equation-of-motion coupled cluster singles, doubles, and triples (EOM-CCSDT) results,⁷⁷ which furnish a vertical excitation energy of 285.06 eV.

3.2 Allene

Similar to ethylene, the experimental X-ray absorption spectrum of allene is characterized by a high intensity band origin, around 285.2 eV, with a prominent shoulder at higher excitation energies on the first band. The simulated spectra are shown in Fig. 3 and the dominant NTO hole-particle pairs for the optically bright transitions are shown in Table 2. The ground state minimum structure of allene has D_2d point group symmetry, and the three C 1s orbitals transform as 1a₁, 1b₂ and 2a₁. The 1a₁ orbital is localized on the central C atom, the 1b₂ and 2a₁ orbitals correspond to the anti-symmetric and symmetric combinations, respectively, of the terminal C 1s orbitals. This latter pair of orbitals are near-degenerate and exhibit a splitting of 0.002 eV as computed with the QTP17 exchange–correlation functional.

Fig. 3 The experimental X-ray absorption of ethylene compared to simulations performed employing a hierarchy of approximations. From bottom, the I_VE is a purely electronic spectrum employing vertical excitation energies and oscillator strengths only (green line), I_HV is Poisson distribution spectrum corresponding to vibrational structure arising from treating the excited states as a series of displaced harmonic oscillators (blue line), and I_NA, the full non-adiabatic vibronic coupling simulation (red line). The experimental spectrum is at top (black line).⁵² The bottom panels evince a detailed comparison of the various approximations with the experimental band origin.

Table 2 Symmetry (Sym), theoretical best estimate vertical excitation energies (ε_i), oscillator strengths f_i, and electronic characters and dominant natural transition orbital (NTO) pair for core-excited states of allene. All energies calculated at the ground-state minimum (FC point) are given in eV

Sym	Transition	ε i (eV)	f i	Assignment	Hole	Particle
E	A	285.58	0.0564	1sπ*
		285.58	0.0564	1sπ*
E	B	285.79	0.0417	1sπ*
		285.79	0.0417	1sπ*
E	C	286.56	0.0055	1s3py
		286.56	0.0055	1s3px
E	D	288.51	0.0146	1s3px
		288.51	0.0146	1s3py

---

The origin band in Fig. 3 is comprised of multiple electronic transitions, as shown by the I_VE simulation (green line). The lowest energy transition corresponds to an excitation out of the 1a₁ 1s orbital on the central C atom into the degenerate e-symmetry π* (LUMO) orbital. The degenerate NTO pairs are labeled A and are shown in Table 2. The second transition is less than 0.2 eV higher in energy and corresponds to coupled excitations out of the near-degenerate 1b₂ and 2a₁ C 1s orbitals into the LUMO e-symmetry π*. This transition is labeled B in Table 2. There is also a minor contribution to this band, labeled C in Table 2, also an E-symmetry state. This state corresponds to excitations from the terminal C 1s orbitals into a degenerate e(3p) Rydberg orbital.

Incorporating vibrational effects into the absorption spectrum simulation (Fig. 3, blue line) results in a number of quantitative changes to the XAS simulation, beyond the addition of vibrational structure. Specifically, the I_HV simulation (Fig. 3, blue line) yields a pronounced red shift of the leading edge of the origin band and accurately reproduces the peak positions of the main feature at 285.3 eV and the shoulder at 285.6 eV, thereby improving agreement with the experimental spectrum in terms of both peak positions and overall spectral shape. Although the I_HV spectrum evinces improved peak positions and overall spectral shape, it fails to accurately reproduce the experimentally observed bandwidths, particularly for the band origin as evident in the middle plot of the bottom panel in Fig. 3. These results are in qualitative agreement with previous simulations in which vertical excitation energies were computed for an ensemble of nuclear configurations generated from ground state nuclear density, thereby describing Born–Oppenheimer vibronic effects.⁸⁸ Although increasing the Lorentzian broadening to FWHM = 0.26 eV (5 fs lifetime) or even 0.53 eV (2.5 fs lifetime) yields somewhat broader bands, the experimentally observed widths are still not well reproduced, and essential spectral features such as the shoulder at 285.6 eV disappear; thus, increasing Lorentzian broadening does not provide a meaningful improvement as shown in Fig. S11 of SI.

The inclusion of non-Born–Oppenheimer vibronic coupling yields an XAS in near quantitative agreement with experiment. This simulation, given by the red curve in Fig. 3 furnishes an origin band that is significantly broader than the I_HV simulation and is in excellent agreement with the experimental result. The calculated excitation energy spectra was adjusted by values ranging from +0.05 to 0.35 eV to align absorption band with experiment. The un-shifted energies and spectra are shown in Table S4 and Fig. S7 as SI, with the theoretical best estimate vertical excitation energies, ε_i, given in Table 2. A damping time of 7.5 fs was independently determined for allene to furnish the best agreement with the experimental spectrum.⁵² Similar to ethylene, the comparison in Fig. S8 (SI) of a model which includes anharmonic terms, but not inter-state coupling (orange line), with one that also includes inter-state coupling (cyan line), shows that the additional broadening originates predominantly from non-adiabatic effects rather than anharmonicity. Moreover, including bilinear coupling terms in the full non-adiabatic model (red line) produces only negligible changes, indicating that their influence on the spectrum of allene is minor.

The vibrational modes that most strongly couple the electronic states of the origin band are the degenerate C–C–C bending (Q_10x and Q_10y) and CH₂ twisting (Q₄) vibrational modes (see Table S5 and Fig. S6 in SI). A detailed comparison between spectra obtained with and without vibronic coupling provides valuable insights into the dominant mode influencing the absorption spectrum. When the large inter-state couplings between excited states along these modes are removed the resulting spectrum (see Fig. S9 in SI) closely resembles HV simulation in terms of the width of the origin band. This indicates that these three modes are the primary contributors to significant changes in the spectral shape when including vibronic coupling. Importantly, this also indicates that the additional broadening observed when transitioning from the HV to NA model arises primarily from non-adiabatic effects, rather than from excited-state mode coupling or the inclusion of anharmonicity in the diabatic potentials.

3.3 Butadiene

Fig. 4 presents a comparison of experimental and simulated static X-ray absorption spectra for butadiene. To evaluate the present simulations, we compared the simulated XAS of trans-1,3-butadiene with two experimental spectra: a recently measured high-resolution XAS spectrum (obtained from energy-dependent ion-yield measurements),⁵² and an older result obtained from energy loss spectroscopy (EELS).⁸⁹ The high-resolution XAS spectrum exhibits what appears to be vibrational structure overlaying the electronic transitions that comprise the origin band. This structure is absent in the lower-resolution EELS spectrum. The more fundamental difference, however, can be observed in a significantly broader origin band in the XAS compared to the EELS spectrum, with the former also exhibiting a shoulder near 285.5 eV that is absent in the former. The relative intensities between the origin band and the feature around 287.5 eV are also visibly different. Further comparison of these experimental spectra with the simulated result will be discussed below.

Fig. 4 The experimental X-ray absorption of ethylene compared to simulations performed employing a hierarchy of approximations. From bottom, the I_VE is a purely electronic spectrum employing vertical excitation energies and oscillator strengths only (green line), I_HV is Poisson distribution spectrum corresponding to vibrational structure arising from treating the excited states as a series of displaced harmonic oscillators (blue line), and I_NA, the full non-adiabatic vibronic coupling simulation (red line). The experimental XAS spectrum (black line)⁵² and EELS spectrum (magenta line)⁸⁹ are at top. The bottom panels evince a detailed comparison of the various approximations with the experimental band origin.

The XAS spectrum of butadiene exhibits significantly more structure than that of ethylene and allene. The experimental spectrum is comprised of three prominent bands, with the two lower-energy bands, between 284–286 eV, strongly overlapped. Turning first to the electronic spectrum, it is observed that the vertical excitation energies and oscillator strengths alone roughly reproduce the dominant features in the spectrum in terms of position and relative intensities of peaks. Table 3 lists the lowest optically bright core-excited states of butadiene, along with their corresponding oscillator strengths and dominant NTO pairs. In contrast to allene and ethylene, bringing the experimental and simulated spectra into maximal agreement required no shift of the ab initio vertical excitation energies for the states comprising the band origin. Rather, an energy shift of +0.5 eV was applied exclusively to all Rydberg transitions (unshifted spectra are shown in Fig. S13 in SI). The most intense spectral features are the first two peaks, labeled A and B. These transitions originate from excitations involving in-phase combinations of C 1s AOs located on the terminal and central carbon atoms, respectively, to the π* orbital. As shown in Table 3, the energy difference between peaks A and B is approximately 0.6 eV. The third band in the overall simulated XAS spectrum is a broad, featureless peak arising from the overlapping contributions of at least three distinct electronic transitions (Peaks C, D, and E). These transitions involve excitations into the n = 3 Rydberg manifold.

Table 3 Symmetry (Sym), theoretical best estimate vertical excitation energies (ε_i), oscillator strengths f_i, and electronic characters and dominant natural transition orbital (NTO) pair for core-excited states of butadiene. All energies calculated at the ground-state minimum (FC point) are given in eV

Sym	Transition	ε i (eV)	f i	Assignment	Hole	Particle
B2	A	284.43	0.099	1sπ*
B1	B	285.06	0.084	1sπ*
B2	C	287.16	0.010	1s3s
B2	D	287.53	0.030	1s3px
B1	E	287.64	0.011	1s3s

---

As shown in Fig. 4, the experimental XAS origin band maxima are clearly red-shifted relative to the TBE vertical excitation energies (Fig. 4, green line). Furthermore, the I_VE simulation predicts the lower energy transition to be higher intensity, in contrast to what is observed in the experimental XAS. When the vibrational structure is considered, as shown by the I_HV spectrum (Fig. 4, blue line) agreement with XAS experiment is improved, particularly the splitting between the first two bands. However, the relative intensities of the vibrational peaks for these two bands do not match the XAS experimental observations. Furthermore, while the relative positions of the first two peaks in the I_HV spectrum are accurate, band widths are significantly narrower than those in the experimental spectrum. Additionally, the harmonic vibrational structure model produces three distinct peaks in the 287–289 eV region. In the XAS spectrum, a three-peak structure is indeed discernible, although the features are broader and less pronounced than in the HV model. In contrast, the EELS spectrum displays mainly a single broad peak with a comparatively weaker third feature. Thus, while the HV model captures aspects of the band structure seen in XAS, it also yields excessive vibrational structure within each peak.

In contrast, the full vibronic coupling model reproduces the main features of the pre-edge region more accurately. As shown in Fig. 4 (red line), this approach attenuates the overly pronounced vibronic structure of the harmonic model, leading to broader and smoother spectral features. The resulting spectrum shows a mixed agreement with experiment: for the first band around 286 eV, the simulated profile is narrower and lacks the shoulder feature, in closer agreement with the EELS spectrum. In contrast, for the second band around 287.5 eV, the calculation captures the multi-peak structure more consistently with the XAS measurement. Comparison of an anharmonic model with the HV model (Fig. S14 of SI) shows that anharmonicity moderately alters the spectral shape and broadening, but does not significantly improve agreement with experiment. By contrast, the inclusion of non-adiabatic effects (cyan line) enhances the relative intensity of the first band, leading to a better match with the experimental data. Furthermore, unlike in ethylene and allene where bilinear contributions had only a minor impact, in butadiene the addition of bilinear coupling terms not only improves the relative intensity of the third peak, yielding closer agreement with the EELS spectrum, but also broadens the origin band, resulting in better consistency with both EELS and XAS measurements.

A damping time of 15 fs was employed in the simulation of the absorption spectrum of butadiene to account for the effects of finite experimental resolution as well as lifetime broadening. Knowing that lifetime broadening is a primary factor in XAS, our simulations suggest broadening due to the vibronic coupling and excited-state mode-mixing, unlike in allene, where vibronic coupling is dominant, both contribute equally to spectral broadening in the XAS of butadiene, as shown in Fig. S14 in the SI.

The present simulations well reproduce the main spectral features of the XAS and EELS spectra, but interestingly, are in better agreement with the EELS results. Specifically, the origin band in the XAS spectrum is significantly broader than that evinced in the EELS spectrum and also exhibits a shoulder around 286 eV that is absent in both the vibronic simulation and EELS result. Additionally, the second absorption band, centered around 287.5 eV, is of approximately equal intensity to the origin band in the XAS spectrum, but is of lesser intensity in the I_NA(E) simulation and EELS spectrum, significantly so for the latter.

For butadiene, the most significant vibronic coupling between states is observed along modes Q₁₈ (CH stretch), Q₁₉ (CH₂ stretching), Q₂₀ (C=C stretch), Q₂₂ (CH bending), and Q₂₄ (C–C–C deformation) within b_u symmetry. A schematic representation of all vibrational normal modes of butadiene is shown in Fig. S12 in the SI. Turning off the inter-state coupling along these modes results in an absorption spectrum that closely matches the I_HV model, particularly in the relative intensity of the vibrational progression of the first band. This indicates their high impact on spectral broadening and reflect their strong activation due to the vibronic coupling between electronic excited states with suitable symmetry within photoexcitation (Fig. S15, maroon line, in the SI).

4 Discussion

The XAS simulations in the preceding section evince some general trends with respect to the level of approximation employed in the computation. Beginning with the XAS simulation of ethylene, it is apparent that while the electronic spectrum predicts the position of the first absorption peaks quite accurately, it significantly underestimates the band widths, the relative intensities, and the energetic positions of the peaks at higher energy. While increasing the broadening parameter could improve the bandwidth of the band origin, it would lead to over-broadening and loss of spectral structure at higher energies. Most notably, and unsurprisingly, the I_VE spectral fail to reproduce the prominent shoulders in the origin bands that arise to the vibrational structure. Furthermore, as the number of equivalent cores increases, concomitantly with the density of core-excited states, the ability of the I_VE to even qualitatively reproduce the structure of the absorption bands degrades significantly. In the I_VE simulation of butadiene, the separation between the two prominent peaks in the band origin and their relative intensity are incorrectly predicted. Lastly, we note that for the examples studied here, even when the ab initio vertical excitation energies are quantitatively accurate, the band maxima may differ by up to 0.5 eV (typically red-shifted) from the empirically corrected theoretical best estimate transition energies obtained from the spectral simulations employing the full vibronic coupling Hamiltonian. As the number of electronic structure methods for computing core-excitation energies continues to increase, this has implications regarding the use of experimental spectral data for benchmarking quantum chemistry approaches.^90–92 Specifically, this practice should be avoided if high accuracy is required, given that vibronic structure may result in significant differences between the computed vertical excitation energies and the band maximum of the corresponding transition band in the XAS.

Incorporation of vibrational structure within the Born–Oppenheimer approximation furnishes a significant improvement over the purely electronic spectra, particularly in the ability to predict peak positions. In ethylene, the I_HV spectrum reproduces the shoulder in the origin band at 285.5 eV as well as providing close agreement with experimental peak widths. Similarly, the inclusion of vibrational structure in the I_HV spectrum of allene induces a red shift, improving its alignment with the experimental results, and predicts a shoulder at 285.6 eV in the origin band, thereby enhancing the overall spectral shape. However, the I_HV spectrum still underestimates the broadening of the origin band, centered at 285.3 eV, and fails to fully reproduce the intensity ratio between the origin band and the second peak at 288.5 eV, which remains more pronounced in the experiment. For butadiene, the I_HV spectrum predicts better separation between peaks A and B and captures the vibrational progression observed for both peaks. However, the relative intensities of the vibrational peaks at peaks A and B are overestimated compared to the experimental spectrum. Moreover, the I_HV spectrum fails to replicate the broad third band observed in the experimental spectrum in the 287–289 eV range, instead predicting three distinct peaks within this region. That said, given that Poisson distribution simulations only require the equivalent of a single gradient evaluation for each excited state (to determine the shift of the harmonic potential), the present results indicate that this is a computationally cost-effective approach to dramatically improve XAS spectral simulations.

The XAS simulations of ethylene based on the full vibronic coupling models (I_NA) employed here uniformly furnish near quantitative agreement with the experimental spectra. For the allene simulation, incorporating non-BO vibronic coupling effects corrects the underestimation of the experimentally observed first-band width at 286.0 eV and the underestimated relative intensity of the second band at 288.5 eV, which were observed when using the harmonic vibrational model. Notably, the simulated bandwidth now closely matches the experimental observations, and the total absorption spectrum, shown by the red line in Fig. 3, exhibits excellent agreement with the experimental spectrum.

Similarly, for butadiene, the inclusion of non-Born–Oppenheimer vibronic coupling yields an XAS spectrum in very good agreement with experimental results, particularly with the EELS spectrum. As shown in Fig. 4, the full vibronic coupling model (red line) accurately reproduces the broad third band observed in the experimental spectra, and corrects the pronounced vibrational intensities in the first two peaks (centered at 284.7 and 285.4, respectively) present in the harmonic model (Fig. 4, blue line), and enhances spectral broadening, leading to a simulated spectrum that more closely resembles experimental observations. Both of the experimental spectra shown in Fig. 4 are ‘pseudo-XAS’ and were measured via changed particle detection rather than transmittance: ion-yields in the case of ref. 52 and in-elastically scattered electrons for ref. 89. A convincing explanation for the difference between these measurements and the simulation (and why, for example, the I_NA seems to be in better agreement with the EELS result) cannot be given by the authors in this work. However, the ability to perform highly accurate first-principles simulations of X-ray spectroscopy will ideally encourage more detailed comparisons between not only theory and experiment, but also between different experimental techniques.

Analysis of the XAS spectral simulations reveals that vibronic coupling occurs predominantly within sub-manifolds of the full set of core-excited states. Specifically, vibronic coupling is operative primarily between those states involving transitions from sets of core orbitals localized on symmetry equivalent C atoms, while coupling between transitions in which the dominant hole NTOs are localized on symmetry inequivalent core-orbitals is observed to be, for the examples considered here, negligible. To visualize this empirical observation, Fig. 5 shows simulations for allene and butadiene in which the vibronic coupling between manifolds of core-excited states involving transitions from symmetry in-equivalent carbon atoms has been removed. For these simulations, the total spectrum can be constructed as a sum of the two non-interacting manifolds of core-excited states. As Fig. 5 evinces, when these approximate spectra are compared to the full I_NA spectra, the agreement is near quantitative.

Fig. 5 Comparison of absorption spectra for allene and butadiene, illustrating the contributions of excited states based on their origins. Excited states are categorized by transitions originating from the central carbon atom (green line) versus those originating from terminal carbon atoms (orange line). The solid black area represents the sum of these contributions, while the dashed blue line shows the exact spectra obtained by employing the full vibronic coupling Hamiltonian. Top panel: allene; Bottom panel: butadiene.

Recent work has shown that vibronic coupling between core-excited states can lead to extremely rapid wave packet dynamics that gives rise to charge transfer between different core–holes.⁹³ Indeed, these charge-transfer dynamics can occur on sub-10 fs timescales, thereby potentially occurring prior to Auger decay. However, a prerequisite for these dynamics is the presence of strong vibronic coupling between core-excited states involving transitions from different core-orbitals. The present calculations suggest that these dynamics can only occur between states corresponding to excitation from core orbitals localized on symmetry-equivalent atoms, as only transitions from these core-orbitals give rise to strongly vibronically coupled states. While this result is evocative, further simulation is required to confirm the generality of this observation.

5 Conclusion

We have here evaluated the accuracy and reliability of our recently proposed method for the determination of quasi-diabatic potentials and couplings within the DFT/MRCI(2) framework, termed QD-DFT/MRCI(2), by simulating the XAS of a set of unsaturated small-chain hydrocarbons. The QD-DFT/MRCI(2) approach provides a route to the direct determination of quasi-diabatic potential surfaces and potential couplings, greatly facilitating the evaluation of the Taylor expansion coefficients required to parameterize the vibronic coupling Hamiltonians employed here. This is particularly true for expansions involving many vibrational degrees of freedom, and/or very many electronic states. The dense electronic manifolds that result in the X-ray absorption spectra of polyatomic molecules thus serve as particularly compelling use-case for the present method.

Furthermore, an additional motivation for this work is to present a comparison of frequently employed approximation schemes for the simulation of these spectroscopic quantities. This study aims to quantify to the different effects of vibrational structure, anharmonicity, Dushinsky mixing, and non-adiabatic vibronic coupling, in addition to lifetime broadening, of the spectral envelope of molecular X-ray absorption spectra. To that end, the present work shows that while purely electronic spectra, determined solely from vertical excitation energies and oscillator strengths, may provide reasonable zeroth order descriptions of the XAS, they in general fail to capture differential band broadening effects (due to vibrational/vibronic structure). More importantly, vibrational and vibronic structure was observed to result in significant shifts of band maxima, up to 0.5 eV for the examples considered here. Thus, one recommendation that emerges from the present study is to urge caution in the use of experimental data, and peak maxima in particular, to benchmark vertical excitation energies from electronic structure methods.

A second notable result of the present simulations is the extent to which inclusion of vibrational structure, even at the highly approximate level of a Poisson distribution function, dramatically improves the agreement with experimental spectra. Furthermore, the computational cost of such simulations is not high, as they require a single additional gradient evaluation for each core-excited state. These computations are readily performed for most of the electronic structure methods currently employed to simulate XAS.

The vibronic Hamiltonians employed here, which describe the anharmonicity of the diabatic potentials, as well as mode-coupling in both the diagonal (Dushinsky rotation) and off-diagonal blocks (quadratic inter-state coupling), furnished near quantitative agreement with experimental results for each of the molecules considered here. The vibronic spectra were particularly successful in reproducing the energetic width of the (broad) bands in the XAS. However, given that these spectra were determined from time-dependent wave packet simulations, the spectral envelope will in practice only be determined by the dynamics simulations up to an empirically chosen damping time. That said, it was observed that the damping time that furnished the best agreement with experiment corresponded to the expected core–hole lifetime (7.5–15 fs) for these molecules. This suggests that the vibronic structure was well described by the current Hamiltonians and the remaining band broadening could be attributed to the (short) lifetime of the core-excited states.

The present results show that quantitative agreement with contemporary experimental measurements of the pre-edge structure of molecular X-ray spectra is achievable with high-level simulation that accounts for all electronic and vibrational degrees of freedom. The ability to determine these spectra from first principles enables not only enables more refined analysis of existing experimental results, but also facilitates spectroscopic simulation with useful predictive accuracy.

Conflicts of interest

There are no conflicts to declare.

Data availability

In addition to Cartesian coordinates of all optimized structures, tables summarizing the quantum chemistry computations, plots of vibrational normal modes and additional spectra, all MCTDH operator and input files are provided as supplementary information (SI). Supplementary information: cartesian coordinates of all optimized structures. Additional Tables on the energies and properties of the excited states for optimized geometries. Plot of vibrational normal modes. Vibronic absorption spectra of molecules including individual state component computed. See DOI: https://doi.org/10.1039/d5cp02832a.

Acknowledgements

All authors thank R. Ingle for providing experimental spectra and helpful discussion, and A. Stolow (AS) for helpful discussions. MSS thanks the Natural Sciences and Engineering Research Council of Canada (NSERC) Discovery Grant Program and the National Research Council Canada (NRC) QSP Challenge Program, project QSP-075-01 for financial support.

References

1. M. Chergui, Struct. Dyn., 2016, 3, 031001.
2. C. Milne, T. Penfold and M. Chergui, Coord. Chem. Rev., 2014, 277, 44–68.
3. M. Chergui and E. Collet, Chem. Rev., 2017, 117, 11025–11065.
4. I. Gabalski, F. Allum, I. Seidu, M. Britton, G. Brenner, H. Bromberger, M. Brouard, P. H. Bucksbaum, M. Burt and J. P. Cryan, et al. , J. Phys. Chem. Lett., 2023, 14, 7126–7133.
5. D. Mayer, F. Lever and M. Gühr, Photochem. Photobiol., 2024, 100, 275–290.
6. M. E. Reinhard, B. K. Sidhu, I. B. Lozada, N. Powers-Riggs, R. J. Ortiz, H. Lim, R. Nickel, J. V. Lierop, R. Alonso-Mori and M. Chollet, et al. , J. Am. Chem. Soc., 2024, 146, 17908–17916.
7. S. Oghbaiee, M. Gisselbrecht, N. Walsh, B. Oostenrijk, J. Laksman, E. P. Månsson, A. Sankari, J. H. Eland and S. L. Sorensen, Phys. Chem. Chem. Phys., 2022, 24, 28825–28830.
8. A. Malyzhenkov, Y. P. Arbelo, P. Craievich, P. Dijkstal, E. Ferrari, S. Reiche, T. Schietinger, P. Juranić and E. Prat, Phys. Rev. Res., 2020, 2, 042018.
9. J. Duris, S. Li, T. Driver, E. G. Champenois, J. P. MacArthur, A. A. Lutman, Z. Zhang, P. Rosenberger, J. W. Aldrich and R. Coffee, et al. , Nat. Photonics, 2020, 14, 30–36.
10. C. Callegari, A. N. Grum-Grzhimailo, K. L. Ishikawa, K. C. Prince, G. Sansone and K. Ueda, Phys. Rep., 2021, 904, 1–59.
11. S.-C. Lin, C.-C. Chang, S.-Y. Chiu, H.-T. Pai, T.-Y. Liao, C.-S. Hsu, W.-H. Chiang, M.-K. Tsai and H. M. Chen, Nat. Commun., 2020, 11, 3525.
12. M. D. Waters, W. Du, A. M. Carrascosa, B. Stankus, M. Cacciarini, P. M. Weber and T. I. Sølling, J. Phys. Chem. Lett., 2021, 12, 9220–9225.
13. F. Ota, K. Yamazaki, D. Sébilleau, K. Ueda and K. Hatada, J. Phys. B: At., Mol. Opt. Phys., 2022, 54, 244002.
14. F. Ota, S. Abe, K. Hatada, K. Ueda, S. Daz-Tendero and F. Martn, Phys. Chem. Chem. Phys., 2021, 23, 20174–20182.
15. C. Wang, M. Gong, X. Zhao, Q. W. Nan, X. Y. Yu, Y. Cheng, V. Kimberg, X.-J. Liu, O. Vendrell and K. Ueda, et al. , Commun. Phys., 2024, 7, 1.
16. S. M. Cavaletto, Y. Nam, J. R. Rouxel, D. Keefer, H. Yong and S. Mukamel, J. Chem. Theory Comput., 2023, 19, 2327–2339.
17. D. Garratt, L. Misiekis, D. Wood, E. W. Larsen, M. Matthews, O. Alexander, P. Ye, S. Jarosch, C. Ferchaud and C. Strüber, et al. , Nat. Commun., 2022, 13, 3414.
18. K. S. Zinchenko, F. Ardana-Lamas, I. Seidu, S. P. Neville, J. van der Veen, V. U. Lanfaloni, M. S. Schuurman and H. J. Wörner, Science, 2021, 371, 489–494.
19. H. Hayashi, Y. Udagawa, W. A. Caliebe and C.-C. Kao, Chem. Phys. Lett., 2003, 371, 125–130.
20. F. Gadea, H. Köppel, J. Schirmer, L. Cederbaum, K. Randall, A. M. Bradshaw, Y. Ma, F. Sette and C. Chen, Phys. Rev. Lett., 1991, 66, 883.
21. K. Hämäläinen, J. Hastings, D. P. Siddons and L. Berman, X-ray absorption spectroscopy beyond the core–hole lifetime, Brookhaven national lab., upton, ny (united states) technical report, 1992.
22. T. Fransson, I. Zhovtobriukh, S. Coriani, K. T. Wikfeldt, P. Norman and L. G. Pettersson, Phys. Chem. Chem. Phys., 2016, 18, 566–583.
23. D. Hait, K. J. Oosterbaan, K. Carter-Fenk and M. Head-Gordon, J. Chem. Phys., 2022, 156, 201104.
24. M. Ehara and K. Ueda, Journal of Physics: Conference Series, 2010, p. 012020.
25. M. Ehara, T. Horikawa, R. Fukuda, H. Nakatsuji, T. Tanaka, M. Hoshino, H. Tanaka and K. Ueda, Journal of Physics: Conference Series, 2011, p. 012024.
26. A. Trofimov, T. Moskovskaya, E. Gromov, N. Vitkovskaya and J. Schirmer, J. Struct. Chem., 2000, 41, 483–494.
27. F. Frati, F. De Groot, J. Cerezo, F. Santoro, L. Cheng, R. Faber and S. Coriani, J. Chem. Phys., 2019, 151, 064107.
28. T. Moitra, D. Madsen, O. Christiansen and S. Coriani, J. Chem. Phys., 2020, 153, 234111.
29. D. M. Holland, J. Suchan, J. Janoš, C. Bacellar, L. Leroy, T. R. Barillot, L. Longetti, M. Coreno, M. de Simone and C. Grazioli, et al. , Phys. Chem. Chem. Phys., 2024, 26, 15130–15142.
30. H. Köppel, W. Domcke and L. Cederbaum, J. Chem. Phys., 1981, 74, 2945–2968.
31. H. Köppel, W. Domcke and L. S. Cederbaum, Adv. Chem. Phys., 1984, 59–246.
32. L. Cederbaum, H. Köppel and W. Domcke, Int. J. Quantum Chem., 1981, 20, 251–267.
33. B. Kempgens, B. Itchkawitz, K. Randall, J. Feldhaus, A. Bradshaw, H. Köppel, F. Gadea, D. Nordfors, J. Schirmer and L. Cederbaum, Chem. Phys. Lett., 1995, 246, 347–352.
34. H. Köppel, F. Gadea, G. Klatt, J. Schirmer and L. Cederbaum, J. Chem. Phys., 1997, 106, 4415–4429.
35. B. Kempgens, B. Itchkawitz, J. Feldhaus, A. Bradshaw, H. Köppel, M. Döscher, F. Gadea and L. Cederbaum, Chem. Phys. Lett., 1997, 277, 436–442.
36. J. Schirmer, A. Barth and F. Tarantelli, Chem. Phys., 1988, 122, 9–15.
37. A. Trofimov, E. Gromov, H. Köppel, J. Schirmer, K. Prince, R. Richter, M. De Simone and M. Coreno, J. Phys. B: At., Mol. Opt. Phys., 2003, 36, 3805.
38. F. Segatta, D. Aranda, F. Aleotti, F. Montorsi, S. Mukamel, M. Garavelli, F. Santoro and A. Nenov, J. Chem. Theory Comput., 2023, 20, 307–322.
39. A. Freibert, D. Mendive-Tapia, N. Huse and O. Vendrell, J. Phys. B: At., Mol. Opt. Phys., 2022, 54, 244003.
40. S. P. Neville and M. S. Schuurman, J. Chem. Phys., 2024, 160, 234109.
41. S. P. Neville and M. S. Schuurman, J. Chem. Phys., 2022, 157, 164103.
42. T. Pacher, L. Cederbaum and H. Köppel, J. Chem. Phys., 1988, 89, 7367–7381.
43. T. Pacher, L. Cederbaum and H. Köppel, Adv. Chem. Phys., 1993, 84, 293–391.
44. D. Klein, J. Chem. Phys., 1974, 61, 786–798.
45. I. Shavitt and L. T. Redmon, J. Chem. Phys., 1980, 73, 5711–5717.
46. C. E. Soliverez, Phys. Rev. A: At., Mol., Opt. Phys., 1981, 24, 4.
47. P. Durand, Phys. Rev. A: At., Mol., Opt. Phys., 1983, 28, 3184.
48. V. Hurtubise and K. F. Freed, Adv. Chem. Phys., 1992, 465–541.
49. S. Grimme and M. Waletzke, J. Chem. Phys., 1999, 111, 5645–5655.
50. I. Lyskov, M. Kleinschmidt and C. M. Marian, J. Chem. Phys., 2016, 144, 034104.
51. C. M. Marian, A. Heil and M. Kleinschmidt, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2019, 9, e1394.
52. R. Ingle, A. Banerjee, C. Bacellar, T. Barillot, L. Longetti, M. Coreno, M. De Simone, F. Zuccaro, L. Poletto and P. Miotti, et al. , J. Phys. B: At., Mol. Opt. Phys., 2022, 55, 044001.
53. H. Köppel, Conical intersections: electronic structure, dynamics and spectroscopy, World Scientific, 2004, pp. 175–204.
54. S. P. Neville, I. Seidu and M. S. Schuurman, J. Chem. Phys., 2020, 152, 114110.
55. H.-D. Meyer, U. Manthe and L. S. Cederbaum, Chem. Phys. Lett., 1990, 165, 73–78.
56. U. Manthe, H.-D. Meyer and L. S. Cederbaum, J. Chem. Phys., 1992, 97, 3199–3213.
57. M. H. Beck, A. Jäckle, G. A. Worth and H.-D. Meyer, Phys. Rep., 2000, 324, 1–105.
58. H.-D. Meyer and G. A. Worth, Theor. Chem. Acc., 2003, 109, 251–267.
59. H. Wang and M. Thoss, J. Chem. Phys., 2003, 119, 1289–1299.
60. U. Manthe, J. Chem. Phys., 2009, 130, 054109.
61. O. Vendrell and H.-D. Meyer, J. Chem. Phys., 2011, 134, 044135.
62. G. A. Worth, M. H. Beck, A. Jäckle, H.-D. Meyer and O. Vendrell, The MCTDH Package, Version 8.5, 2019, See https://mctdh.uni-heidelberg.de.
63. V. Engel, Chem. Phys. Lett., 1992, 189, 76–78.
64. H. S. Carslaw, Introduction to the Theory of Fourier's Series and Integrals, Macmillan and Company, limited, 1921, vol. 1.
65. A. Raab, G. A. Worth, H.-D. Meyer and L. Cederbaum, J. Chem. Phys., 1999, 110, 936–946.
66. G. A. Worth, H.-D. Meyer, H. Köppel, L. Cederbaum and I. Burghardt, Int. Rev. Phys. Chem., 2008, 27, 569–606.
67. A. D. Becke, J. Chem. Phys., 1993, 98, 5648–5652.
68. R. A. Kendall, T. H. Dunning and R. J. Harrison, J. Chem. Phys., 1992, 96, 6796–6806.
69. D. E. Woon and T. H. Dunning Jr, J. Chem. Phys., 1994, 100, 2975–2988.
70. P. J. Stephens, F. J. Devlin, C. F. Chabalowski and M. J. Frisch, J. Phys. Chem., 1994, 98, 11623–11627.
71. TURBOMOLE, TURBOMOLE V6.1, 2009.
72. L. S. Cederbaum, W. Domcke and J. Schirmer, Phys. Rev. A: At., Mol., Opt. Phys., 1980, 22, 206–222.
73. A. Barth and L. S. Cederbaum, Phys. Rev. A: At., Mol., Opt. Phys., 1981, 23, 1038–1061.
74. A. Barth and J. Schirmer, J. Phys. B: At. Mol. Phys., 1985, 18, 867.
75. I. Seidu, S. P. Neville, M. Kleinschmidt, A. Heil, C. M. Marian and M. S. Schuurman, J. Chem. Phys., 2019, 151, 144104.
76. T. S. Costain, J. B. Rolston, S. P. Neville and M. S. Schuurman, J. Chem. Phys., 2024, 161, 224106.
77. T. S. Costain, V. Ogden, S. P. Neville and M. S. Schuurman, J. Chem. Phys., 2024, 160, 224106.
78. B. I. Dunlap, Phys. Chem. Chem. Phys., 2000, 2, 2113–2116.
79. M. Feyereisen, G. Fitzgerald and A. Komornicki, Chem. Phys. Lett., 1993, 208, 359–363.
80. F. Weigend, Phys. Chem. Chem. Phys., 2002, 4, 4285–4291.
81. S. Neville and M. Schuurman, GRaCI: General Reference Configuration Interaction, 2021, https://github.com/schuurman-group/graci.git.
82. Q. Sun, T. C. Berkelbach, N. S. Blunt, G. H. Booth, S. Guo, Z. Li, J. Liu, J. D. McClain, E. R. Sayfutyarova, S. Sharma, S. Wouters and G. K.-L. Chan, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2018, 8, e1340.
83. Q. Sun, X. Zhang, S. Banerjee, P. Bao, M. Barbry, N. S. Blunt, N. A. Bogdanov, G. H. Booth, J. Chen, Z.-H. Cui, J. J. Eriksen, Y. Gao, S. Guo, J. Hermann, M. R. Hermes, K. Koh, P. Koval, S. Lehtola, Z. Li, J. Liu, N. Mardirossian, J. D. McClain, M. Motta, B. Mussard, H. Q. Pham, A. Pulkin, W. Purwanto, P. J. Robinson, E. Ronca, E. R. Sayfutyarova, M. Scheurer, H. F. Schurkus, J. E. T. Smith, C. Sun, S.-N. Sun, S. Upadhyay, L. K. Wagner, X. Wang, A. White, J. D. Whitfield, M. J. Williamson, S. Wouters, J. Yang, J. M. Yu, T. Zhu, T. C. Berkelbach, S. Sharma, A. Y. Sokolov and G. K.-L. Chan, J. Chem. Phys., 2020, 153, 024109.
84. Q. Sun, J. Comput. Chem., 2015, 36, 1664–1671.
85. F. Plasser, M. Wormit and A. Dreuw, J. Chem. Phys., 2014, 141, 024106.
86. Y. Ma, F. Sette, G. Meigs, S. Modesti and C. Chen, Phys. Rev. Lett., 1989, 63, 2044.
87. H. Rabus, D. Arvanitis, M. Domke, A. Puschmann, L. Wenzel, G. Comelli, G. Kaindl and K. Baberschke, Phys. Scr., 1990, 1990, 131.
88. E. Muchova, D. Hollas, D. M. Holland, C. Bacellar, L. Leroy, T. R. Barillot, L. Longetti, M. Coreno, M. de Simone and C. Grazioli, et al. , Phys. Chem. Chem. Phys., 2023, 25, 6733–6745.
89. R. N. Sodhi and C. Brion, J. Electron Spectrosc. Relat. Phenom., 1985, 37, 1–21.
90. F. J. Avila Ferrer, J. Cerezo, E. Stendardo, R. Improta and F. Santoro, J. Chem. Theory Comput., 2013, 9, 2072–2082.
91. S. Bai, R. Mansour, L. Stojanović, J. M. Toldo and M. Barbatti, J. Mol. Model., 2020, 26, 1–9.
92. B. Lasorne, J. Jornet-Somoza, H.-D. Meyer, D. Lauvergnat, M. A. Robb and F. Gatti, Spectrochim. Acta, Part A, 2014, 119, 52–58.
93. S. P. Neville, M. Yaghoubi Jouybari and M. S. Schuurman, J. Phys. Chem. Lett., 2025, 16, 8885–8890.

---