We present a time-dependent theory for non-resonant x-ray emission spectrum (XES) and normal Auger spectrum (NAS) calculation, based on a fully quantum description of nuclear dynamics using the vibrational wave packet concept. We compare two formulations of the time-dependent theory, either employing a two-time propagation scheme or using spectral integration over the electron energy continuum. We find that the latter formulation is more efficient for numerical simulations, providing a reasonable accuracy when the integration step is shorter than the lifetime broadening of the core-ionized state. We demonstrate our approach using the example of non-resonant x-ray emission from a water molecule, considering the lowest core-ionized K−1 and first core-ionized shake-up K−1V−1V1 intermediate states. These channels exemplify the developed theory on bound–bound, bound–continuum, continuum–bound, and continuum–continuum transitions. Our results suggest that the time-dependent approach is efficient for simulating XES involving dissociative states, whereas the time-independent approach, based on Franck–Condon factors, is more efficient for bound–bound transitions expressed as discrete frequency dependence in the energy domain. The methods and discussion have general applicability, including both NAS and more complex systems, such as liquid water.
I. INTRODUCTION
In modern x-ray spectroscopy, much attention has been paid to two akin spectroscopies, namely resonant inelastic x-ray scattering (RIXS) and resonant Auger spectroscopy (RAS), allowing for high spectral resolutions and, thus, studies of nuclear dynamics in various systems.1 RIXS and RAS processes can be modeled theoretically in a very similar way using two qualitatively different approaches for the description of the vibrational structure. The first, time-independent, approach is based on the implementation of the Kramers–Heisenberg equation, which allows us to obtain spectra directly in the energy domain. This method is well suited for systems with a limited number of degrees of freedom and for transitions between bound electronic states. However, quite often, intermediate core-excited and final valence-excited states of RIXS/RAS processes manifest a dissociative character. In those cases, the required integration over the vibrational continuum in a dissociative potential makes the energy-domain calculations numerically expensive. It was recognized, however, that the intrinsic difficulty with a high density of nuclear states can be overcome using a time-dependent wave packet technique,2–5 which constitutes an alternative approach for the modeling of vibrationally resolved RIXS and RAS spectra. Besides its numerical benefits, this method also brings a clear physical picture of the time-dependent nuclear dynamics induced by x rays.6,7
Besides the resonant x-ray emission processes, non-resonant processes, involving x-ray induced core-ionization with a photon energy tuned above the ionization threshold, play an important role in modern x-ray spectroscopy. As in the case of resonant processes, the created core-vacancy is filled by a valence electron with the simultaneous emission of a secondary x-ray photon or an Auger electron but here giving rise to non-resonant x-ray emission spectroscopy (XES)8 or normal Auger spectra (NAS).9 In addition, measurements of XES and NAS intrinsically include fingerprints of nuclear dynamics in the intermediate (core-ionized) and final (valence-ionized) states. To the best of our knowledge, the theory of these non-resonant x-ray and Auger emission processes, including the nuclear dynamics, has so far been described only using the time-independent method based on the Kramers–Heisenberg equation and the calculation of the Franck–Condon amplitudes.10 However, quite often, the core-ionized state also shows a dissociative character, making the energy-domain calculations expensive and inappropriate. For example, this is the case of liquid water11–13 and the HCl molecule.5 The goal of this study is to present a numerical implementation of a time-dependent theory developed for the simulation of XES and NAS. The energy-domain calculations for non-resonant x-ray emission processes are simpler than the corresponding calculations of RIXS and RAS. The implementation of the time-dependent formulation for XES, however, is non-trivial and becomes more complex than the wave packet theory of resonant x-ray scattering. For the sake of generality, we discuss here two qualitatively different time-domain approaches to how the wave packet technique can be implemented for XES and NAS, either based on a two-time propagation scheme or based on integration over the energy continuum of the emitted photoelectron. From a simple estimation of the numerical costs for the two schemes, we note that the energy integration scheme shows a reasonable convergence already when the integration step is equal to or smaller than the lifetime broadening of the core-ionized state. For the full convergence, however, the integration step needs to be smaller than the typical experimental resolution. Due to this, the latter scheme was implemented in our numerical protocols and used in the case study of XES of the gas-phase water molecule. The approach that we developed is general and can also be applied for NAS simulations. Moreover, it can be very useful in the case of more complex multi-dimensional systems with dissociative core-ionized potential energy curves (PECs), e.g., ice,14 liquid water,8,11,13 and biomolecules in aqueous solutions.15 It is worth noting that in the case of core-ionization in the vicinity of shape resonances, the XES shows some features of RIXS spectroscopy, as it was shown for N2 in theory16 and experiment.17 In the present study, we consider ionization far from the threshold, where the electron continuum has no specific features.
This paper is organized as follows: Our theoretical approach is outlined in Sec. II, where we start with a general discussion of approaches to resonant x-ray scattering (Sec. II A), which are then adapted for the simulation of non-resonant x-ray emission in energy (Sec. II B) and time (Sec. II C) domains. The numerical costs for different time-dependent approaches are estimated and discussed in Sec. II D. The results and discussions for the showcase of a water molecule are summarized in Sec. III, in which the results of ab initio simulations are presented in Sec. III A, the vibrationally resolved XES via two core-ionized states is discussed in Sec. III B, and the numerical convergence is considered in Sec. III C. Our findings are summarized in Sec. IV.
II. GENERAL THEORY OF NON-RESONANT X-RAY EMISSION AND NORMAL AUGER ELECTRON EMISSION
Before formulating the expression for the XES cross section, let us remind the reader about the basic equations for the RIXS cross section. In what follows, we assume the validity of the Born–Oppenheimer approximation and consider channels via a core-excited (c) electronic state, which are energetically well separated from other core-excited states. We label the ground electronic state with 0 and the final electronic states with f. Atomic units are used throughout this paper unless otherwise stated.
A. Time-independent vs time-dependent formalism for resonant x-ray scattering
B. Non-resonant XES cross sections in the energy domain
C. Non-resonant XES cross sections in the time domain
Let us now introduce a representation of XES cross section fully in the time domain. Here, we discuss two possible approaches: (i) the strict derivation of the cross section and simulations using two-time propagation and (ii) the calculation of the XES cross section in RIXS mode using integration over electron kinetic energy to treat the non-resonant behavior in XES (lack of energy conservation in the fluorescence decay) due to electron ejection.
- Two-time forward–backward propagation scheme. To solve the problem mentioned above, we follow the idea outlined in Refs. 5 and 19 and introduce a two-time wave packet asHere, the forward propagation of the packet |Ψc(t)⟩ with the Hamiltonian hc of the core-ionized state along the time t and initial condition |Ψc(t)⟩ = |0⟩ is followed by the backward propagation of the wave packet Φ(τ, t) in the final state with the Hamiltonian hf along the time τ from τ = 0 to τ = −t (t ≥|τ|≥ 0) and initial condition |Φc(0, t)⟩ = dfc|Ψc(t)⟩; see Fig. 1(a). Then, the one-time wave packet |Ψfc(t)⟩ (10) can be obtained as(11)and the XES cross section (9) can now be computed as half-Fourier transform of the auto-correlation function(12)with the wave packet defined in (12); the ways to implement the two-time propagation scheme are outlined in Appendix A.(13)
- XES profile as integral RIXS cross section. The XES cross section, however, can also be directly computed by the numerical integration of the RIXS cross section over incident photon energy, as it was discussed above,where the partial RIXS cross section σ(ω′, ω) can be computed using the time-dependent representation (4). Using this protocol, one has to compute RIXS cross section for a large number of incoming photon energies in order to get smooth numerical integration in (14).(14)
(a) Schematic illustration of trajectories related to the two-time propagation technique (11) with forward propagation in the core-ionized state and backward propagation in the final state with Hamiltonians hc and hf, respectively. (b) Illustration of the number of time loops N = T/Δt needed to compute the wave packet (12); T is the final propagation time.
(a) Schematic illustration of trajectories related to the two-time propagation technique (11) with forward propagation in the core-ionized state and backward propagation in the final state with Hamiltonians hc and hf, respectively. (b) Illustration of the number of time loops N = T/Δt needed to compute the wave packet (12); T is the final propagation time.
D. Comparisons of numerical costs for the two-time propagation and integral RIXS cross section schemes
Let us now analyze and compare the computational costs for the two time-dependent XES implementations outlined above: (1) two-time propagation (11)–(13) and (2) XES in RIXS mode (14). The two-time propagation framework represents the strict transfer of the energy-domain Eq. (6) to the time domain. We clearly see that despite the simple energy-domain representation of the XES cross section, the time-dependent representation becomes considerably more complex. Moreover, the two-time propagation also requires a large computation time, as compared to RIXS simulations. Indeed, in the case of the time-dependent RIXS scheme (4), a single propagation of the wave packet in the core-excited state and a single propagation in the final states are required, which provides the numerical costs for RIXS simulations to scale linearly with the number N of time steps Δt. In the case of XES simulation within the two-time propagation scheme (12), the numerical costs scale as N2/2 [see the illustration in Fig. 1(b)], where N = T/Δt and T is the final propagation time.
Let us estimate the costs for a typical XES simulation. As one can see from the definition of the XES cross section (9), the upper limits for numerical time integration can be estimated by the lifetime of the core-ionized state 1/Γc. In order to fully account for the nuclear dynamics in the core-hole and final states, one has to use a full propagation time T, which is larger (an order of magnitude is sufficient) than the core-ionized lifetime T ≳ 10/Γc. Assuming T ≈ 100 fs and a time step Δt ≈ 0.01 fs [required for a smooth convergence of the time-dependent solution of the Schrödinger Eq. (11)], we get N ≈ 104. This simple estimation shows that a XES simulation in the two-time propagation scheme typically requires calculations of about 108 time points in order to obtain the auto-correlation function (13). Moreover, due to the numerical coupling between the forward and backward propagation [Fig. 1(b)], an efficient parallelization procedure cannot be applied for this simulation scheme. Similar scaling of the numerical costs is estimated for other alternative two-time propagation schemes (A3) and (A4).
Let us now evaluate the numerical costs for XES simulations in RIXS mode (14). Contrary to the narrow vibrational resonances in the RIXS profile, defined by the lifetime broadening of the final state Γf, the broadening of a single vibrational resonance of XES profile is significantly larger and is given by Γc + Γf (see Appendix B). Γf is usually negligibly small (nanoseconds for valence-excited states) in comparison with Γc (femtoseconds). In numerical simulations of XES, however, one can safely apply Γf ≪Γc so that Γc + Γf ≈ Γc, providing a linear numerical cost scale for single RIXS spectra at each incident photon energy ω in (14) , N ≈ 10/(ΓcΔt). In order to perform numerical integration in (14), one has to compute the RIXS profile for M energy points, giving the total scale for XES in RIXS mode calculations as N × M. However, in all practical applications, M ≪ N. Indeed, the energy range of integration in (14) is defined by the width of the x-ray photoemission band, which is usually not larger than 1.0 eV,20 and an integration step Δω ≈ Γf is sufficient in this case, which for the choice of Γf = 0.01 eV provides M ≈ 103.
These estimations show that the numerical simulations of XES in RIXS mode are N/M times cheaper than the two-time propagation scheme (N/M ≈ 10 in the typical example considered here). Moreover, the RIXS spectra for each excitation energy ω are fully independent, which allows for using an efficient parallel computational protocol and, thus, gives an enormous improvement on high-performance computer clusters. These estimations clearly show the benefits of XES in RIXS mode scheme (14), which was used for all time-dependent XES simulations presented below.
III. RESULTS AND DISCUSSIONS: XES OF H2O MOLECULE
Potential energy curves and transitions for Channel 1 (blue) and Channel 2 (red). The ionic ground state involved in both the channels is shown in black. The equilibrium bond length of the ground state of neutral water is shown by the vertical dashed line. The energy of the potentials is related to the minimum in the ground state potential of the neutral water molecule. The assignment of the valence states is consistent with the photoelectron spectra.21–24
Potential energy curves and transitions for Channel 1 (blue) and Channel 2 (red). The ionic ground state involved in both the channels is shown in black. The equilibrium bond length of the ground state of neutral water is shown by the vertical dashed line. The energy of the potentials is related to the minimum in the ground state potential of the neutral water molecule. The assignment of the valence states is consistent with the photoelectron spectra.21–24
A. Ab initio calculations of the PECs and R-dependent transition dipole moments
One-dimensional cuts in the potential energy surfaces of the electronic states involved in the transitions in (15) were calculated along the bond distance of one of the O–H groups associated with R (keeping the other O–H bond fixed) of the isolated water molecule, based on the same geometry, scanning grid, and computational framework as in our previous study of RIXS spectra.25 Important differences are that (i) OpenMolcas, version 22.02,26 was now used, (ii) no additional diffuse Rydberg basis set was employed only atomic ANO-RCC basis sets,27 and (iii) only cationic intermediate and cationic final states were considered for the XES simulations. Since the calculations are performed in Cs symmetry, we can distinctly separate A′ and A″ states, and an active space RAS (1, 9, 0; 1, 13, 0) with the O1s core-orbital in RAS1 and 13 valence orbitals in RAS2 allows for an accurate description of the bond dissociation of ground and ionized states. Separate state-averaging of eight valence-ionized states in each symmetry was performed, and the core-ionized states were calculated state-specifically. One hole is allowed in RAS1, and a total number of nine electrons are used for the cationic states. The ground state potential was taken from a previous study.25 The final states of the XES process (e.g., cationic valence-excited states) can also be studied within the direct photoemission process (photoelectron spectra),21,22 which has recently been applied to dynamical studies of H2O molecule,23,24 providing general agreement with PECs obtained in the present article.
Using this approach, we compute R-dependent transition dipole moments between the core-ionized and final states and PECs of all considered electronic states. The electronic XES (neglecting the vibrational structure) are shown in Figs. 3(a) and 3(b) for Channel 1 and Channel 2, respectively. The assignments of the final ionic states of the water molecule are shown23,28 in the figure; the molecular frame geometry is shown in Fig. 4(a). The spectra were obtained using the convolution with a Lorentzian function of width 0.4 eV to mimic the resolving power of the x-ray spectrometer29 used in Ref. 28.
Electronic (vibrational contribution is omitted) non-resonant XES via (a) K−1 core-ionized state (Channel 1) and (b) K−1V−1V1 shake-up state (Channel 2). The electronic XES is convoluted with a Lorentzian function of width 0.4 eV. The binding energy scale BE = ω′ − I1s is used, where I1s is the ionization potential of the O1s orbital. The vertical axes are scaled as the squared total decay transition dipole moment (a.u.) at the ground state equilibrium; see Fig. 4 for the components.
Electronic (vibrational contribution is omitted) non-resonant XES via (a) K−1 core-ionized state (Channel 1) and (b) K−1V−1V1 shake-up state (Channel 2). The electronic XES is convoluted with a Lorentzian function of width 0.4 eV. The binding energy scale BE = ω′ − I1s is used, where I1s is the ionization potential of the O1s orbital. The vertical axes are scaled as the squared total decay transition dipole moment (a.u.) at the ground state equilibrium; see Fig. 4 for the components.
Molecular frame geometry (a). R-dependent components of the transition dipole moment on Channel 1 (b)–(d) and Channel 2 (e)–(f) to those considered ionic final states (see the panel titles). Only components with a significant contribution are shown.
Molecular frame geometry (a). R-dependent components of the transition dipole moment on Channel 1 (b)–(d) and Channel 2 (e)–(f) to those considered ionic final states (see the panel titles). Only components with a significant contribution are shown.
PECs of the electronic states for the two channels (15) discussed here are presented in Fig. 2, where the transitions corresponding to Channel 1 and Channel 2 are shown with the blue and red arrows, respectively. One can see that Channel 1 involves only bound–bound transitions for the ionization and emission, while Channel 2 involves bound–continuum (core-excitation), continuum–continuum, and continuum–bound transitions. The R-dependent transition dipole moment components for all considered final states are shown in Fig. 4. Taking into account the R-dependence of the dipole moment plays an important role in accurately calculating some spectral features, as it will be discussed below.
B. Vibrational XES simulations
Let us now present the vibrationally resolved XES modeled using the time-dependent wave packet approach (14). In Fig. 5(a), the individual contributions of the three final states of Channel 1 are shown: the ground cationic state (red line), the first valence excited state (blue line), and the second valence excited state (green line); the total spectrum is shown as the black curve. For comparison, we also performed simulations in the energy domain employing Eq. (6), which is presented by the black dashed curve in Fig. 5(a). Time-dependent wave packet simulations are in very good agreement with the simulations based on the computation of the Franck–Condon amplitudes. It is worth noting that time-independent simulations for the case of bound–bound transitions are much cheaper than the wave packet approach. The R-dependence of the transition dipole moment [Figs. 4(b)–4(d)] was taken into account in both the types of simulations presented in Fig. 5(a); however, it plays a minor role in the XES profile formation. Moreover, single dipole moment components dx, dy, and dz majorly contribute to the ground state equilibrium geometry for [Fig. 4(b)], [Fig. 4(c)], and [Fig. 4(d)], respectively. Note that the molecular orientation is shown in Fig. 3.
(a) Channel 1. Experimental XES of gas-phase water (light-blue dotted line)28 compared with simulated vibrationally resolved XES obtained using the energy-domain (black dashed line) and time-domain (black line) approaches. The contribution decay channels to the three final states are outlined (see the legends). (b) Channel 2. Total theoretical XES computed in the time-domain approach (black line). The molecular band (MB) and fragment peaks (FPs) are marked in the spectra. The contributions from dy and dz components to each final state are outlined (see the legends).
(a) Channel 1. Experimental XES of gas-phase water (light-blue dotted line)28 compared with simulated vibrationally resolved XES obtained using the energy-domain (black dashed line) and time-domain (black line) approaches. The contribution decay channels to the three final states are outlined (see the legends). (b) Channel 2. Total theoretical XES computed in the time-domain approach (black line). The molecular band (MB) and fragment peaks (FPs) are marked in the spectra. The contributions from dy and dz components to each final state are outlined (see the legends).
The upper profile in Fig. 5(a) presents the experimental XES of a gas-phase water molecule reproduced from Ref. 28. One can see good qualitative agreement between experiment and theory, while some discrepancies deserve a comment. The experimental XES profile shows a broader structure for each final electronic state, without a vibrational resolution. There are two possible reasons for this: (i) the limited spectral resolution in the experiment ( eV28,29) and (ii) the contributions from other vibrational degrees of freedom are not taken into account in our simulations. Indeed, in our illustrative simulations, we consider only stretching nuclear dynamics along a single OH bond, while all other degrees of freedom are kept frozen. However, the bending motion can give a strong contribution to the XES profile. The bending mode is not excited in the core-ionized state as it follows from our ab initio simulations: the equilibrium H–O–H angle of the core-ionized state is the same as in the molecular ground state (see Appendix B, Fig. 8). In the ground cationic states, the equilibrium H–O–H angle is preserved unchanged compared to the ground state of the neutral molecule ( Appendix B, Fig. 8), which results in a narrow XES profile for this state in the experiment [Fig. 5(a), thin black line]. The additional broadening observed in the experiment can be explained by the limited spectral resolution. For the cationic excited states and , the equilibrium bond angle is 180° and 65°, respectively, according to our calculations ( Appendix B, Fig. 8) and the literature.30 Such a strong variation of the bond angle, relative to the neutral ground state equilibrium angle of 104°, results in inducing a strong vibration excitation of the bending mode, which is manifested as broadening of the experimental XES profile for these excited states. The vibrational structure is unresolved due to the limited experimental spectral resolution.
Let us now discuss the theoretical XES profile arising via the dissociative core-ionized shake-up state, corresponding to Channel 2 in Eq. (15) presented in Fig. 5(b). The overall structure of the XES profile consists of a broad feature near 541 eV and three strong and narrow peaks in the region 525–530 eV. This profile can be fully understood from the consideration of the nuclear dynamics involved in the continuum–continuum and continuum–bound decay transitions.
The strongly unbound core-ionized shake-up state results in an ultra-fast dissociation along the OH bond, leading to the appearance of the sharp fragment peaks in the XES profile. The profile structure for the cationic ground and shake-up states are qualitatively different, reflecting the difference between the continuum–bound and continuum–continuum transitions, respectively. Moreover, the contribution from two R-dependent components of the transition dipole moment plays a crucial role here. In the case of the continuum–bound transition to the final state, one can see a weak and broad “molecular band” near 541 eV emission energy and a strong fragment peak near 529.3 eV. Due to a large gradient of the core-ionized state, the vibrational wave packet propagates very fast toward the dissociation region ( Å), resulting in dominant contributions from the dissociation fragment OH. Let us note that, in general, the fragment peak may have a vibrational structure due to the excitation of the OH stretching mode,25 which is, however, neglected in our one-dimensional model. Near the neutral ground state equilibrium position, the dy transition dipole moment component is zero, resulting in no contribution to the molecular band [see the green line for the spectra in Fig. 5(b)]. At the elongated bond distance, the contribution to the fragment line from dy and dz components becomes comparable [see the green and blue lines for the spectra in Fig. 5(b)]. It is worth noting, however, that for both channels, there is a small very asymmetric contribution to the molecular band peak from the dy component.
The structure of the XES profile to the dissociative final state is qualitatively different, showing two narrow structures near 527.4 and 528.4 eV, related to the fragment and pseudo-fragment contributions, similar to the recently observed analog in RIXS on the continuum–continuum transition.25 Due to the similarity of the PECs of the core-ionized and final states, the molecular band shrinks to a narrow peak at 527.4 eV, which can be assigned to the transitions in the bound system R ≤ 2.0 Å. The separate contributions of the dy and dz transition dipole moments clearly confirm this fact [Fig. 5(b)]. Indeed, the dy value is close to zero near the ground state equilibrium, which results in no contribution of this component to the pseudo-fragment peak [see the green line for the state, Fig. 5(b)]. The dipole moment dy increases with an increase in the OH bond length, giving a significant contribution to the fragment peak near 528.4 eV from that component.
C. Remarks on energy integration convergence
The time-dependent wave packet simulation of XES using (14) requires computation of a series of RIXS spectra with changing incoming photon energy and further integration over the energy, which is the main source of numerical expenses. This motivates us to study the general rules for the convergence of the integral over incoming energy. In our simulations, we use core-ionized state lifetime broadening [half-width half-maximum (HWHM)], Γc = 0.08 eV.25,28 For the lifetime broadening (HWHM) of the final states, we choose Γf = 0.01 eV (Γf ≪ Γc), which allows us to properly mimic the XES broadening (see Appendix B) and to reduce the computational costs significantly. We integrate over incoming photon energy in a 2 eV energy range of the core-ionized state varying the integration step. The summary of the convergence study is presented in Fig. 6 for all four types of the decay transitions: (a) bound–bound, the final state of Channel 1; (b) bound–continuum, the final state of Channel 1; (c) continuum–bound, the final state of Channel 2; and (d) continuum–continuum, the final state of Channel 2. For generality, we consider here a dipole forbidden bound–continuum decay transition (b) , which we omitted in total XES simulations [Fig. 5(a)] as its intensity is negligibly small compared to other decay transitions of ionization Channel 1.
Convergence of the XES in RIXS mode integration (14) for energy steps 0.02 eV (orange curves) and 0.005 eV (blue curves); the “resonant” spectra are shown as black lines. (a) Bound–bound , (b) bound–continuum , (c) continuum–bound , and (d) continuum–continuum decay transitions (see the text for details).
Convergence of the XES in RIXS mode integration (14) for energy steps 0.02 eV (orange curves) and 0.005 eV (blue curves); the “resonant” spectra are shown as black lines. (a) Bound–bound , (b) bound–continuum , (c) continuum–bound , and (d) continuum–continuum decay transitions (see the text for details).
In Fig. 6, the black curves show single point RIXS mode calculations on top of the core-ionization resonance (the resonant energy is shown in the figure). The orange curves show the results of integration with an energy step Δω = 0.02 eV, which is smaller than the core-ionized lifetime broadening Γc = 0.08 eV, while the blue curves correspond to Δω = 0.005 eV, which is smaller than the final state broadening Γf = 0.01 eV. In the case of the bound–bound transition [panel (a)], the spectrum of a single incoming energy shows well resolved vibrational peaks (of width Γf), which converges to the vibrational structure broadened with Γc after integration over incoming energy. One can see that already at Δω = Γc/4, the main features of the XES profile are reproduced, while for the full convergence, one has to use an integration step smaller than or about Γf. In the case of limited experimental spectral resolution, an integration step smaller than the spectral resolution is sufficient. A similar behavior is observed for the case of the continuum–bound transition illustrated in Fig. 6(c). Indeed, the general spectral shape in the spectrum is reproduced already at Δω < Γc, while the full convergence is reached when Δω ≤ Γf.
The situation changes qualitatively for the bound–continuum [panel (b)] and continuum–continuum [panel (c)] transitions. Indeed, in these cases, the spectrum computed on top of core-ionization resonance fully coincides with the spectra integrated over incoming photon energy with the integration steps Δω = 0.02 and 0.005 eV. The reason for this is the high energy broadening in the final state, according to the reflection principle,31,32 related to a high gradient of the potential energy curve ΔE/ΔR ≈ 10 eV/Å. Taking into account that the size of the ground state wave function for water in harmonic approximation is Å, with a harmonic frequency of stretching mode ω0 ≈ 0.45 eV, we get the estimation of energy broadening for the bound–continuum transition as 1 eV, which is much higher than the natural XES broadening of Γc and corresponds to the one observed in Fig. 6(b). In the case of continuum–continuum transition, the value of energy broadening for molecular band is similar to the one estimated above [see the “MB” peak in Fig. 6(d)]. The broadening is obviously smaller for the decay far from the equilibrium, where the fragment peak is formed [see the “FP” peak in Fig. 6(d)] due to the smaller potential gradient in that region. Nevertheless, it is still much larger than Γc and, thus, the spectral integration in (14) does not affect the XES profile and can, in principle, be omitted. This is a very interesting observation, which allows us to optimize the computation costs for that kind of transitions with a dissociative final state.
IV. CONCLUSIONS
We developed a time-dependent theory for the vibrational non-resonant x-ray emission spectra. Two ways of computation were proposed: (i) using a two-time propagation scheme and (ii) using XES in RIXS mode, where integration over incoming photon energy was employed. The latter scheme was shown to have the advantage of lower computational costs. Using the water molecule as a showcase, we apply a time-dependent scheme for XES simulations on bound–bound, bound–continuum, continuum–bound, and continuum–continuum transitions. For this, we consider two XES channels: (1) via the lowest core-ionized bound state K−1 and (2) via a core-ionized dissociative shake-up state K−1V−1V1. The final states are the cationic ground state and three lowest cationic valence excited states.
Using high-level ab initio simulations, the PECs and R-dependent transition dipole moments were computed for all considered transitions. The theoretical XES of Channel 1 was compared with the available experimental data, providing a reasonable agreement. For the case of bound–bound transition, our time-dependent simulation was also compared with the energy-domain XES calculations using the generalized Franck–Condon factors (R-dependence of the transition dipole moment was taken into account), showing full agreement between those two methods.
We show that using XES simulations in RIXS mode, the general spectral profile is already obtained when the integration step over incoming photon energy is smaller than the core-ionized state lifetime broadening, while for the full convergence, the integration step smaller than the spectral resolution or the final state lifetime broadening has to be used. We show that on bound–continuum and continuum–continuum transitions, integration over incoming photon energy can be omitted (or use a large integration step) due to the smoothness of the continuum wave functions, which allows for a significant reduction in computational costs.
We conclude that for bound–bound transitions, the energy-domain simulations employing generalized Franck–Condon factors are most efficient, while XES simulations in RIXS mode are preferable for continuum–bound, continuum–continuum, and bound–continuum transitions. The proposed technique developed for XES can be fully applied to the simulation of the normal Auger spectra and extended to a larger number of nuclear degrees of freedom. In particular, the proposed protocol can be applied to simulations of XES from liquid water with a high accuracy and vibrational resolution, using the hybrid quantum–classical approach for treating highly dimensional systems developed earlier.7
ACKNOWLEDGMENTS
V.K. acknowledges the financial support from the Swedish Research Council (Grant No. 2019-03470). M.O. acknowledges the funding from the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie Grant Agreement No. 860553 and the Swedish Research Council (VR Contract No. 2021-04521). A.B. and M.O. acknowledge the funding from the Carl Tryggers Foundation (Contract No. CTS18:285). NI acknowledges Russian Science Foundation (Grant No. 21-12-00193). The computations were enabled by resources provided by the National Academic Infrastructure for Supercomputing in Sweden (NAISS) and the Swedish National Infrastructure for Computing (SNIC) at PDC and NSC partially funded by the Swedish Research Council through Grant Agreement Nos. 2022-06725 and 2018-05973.
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
Author Contributions
Viktoriia Savchenko: Conceptualization (equal); Data curation (equal); Investigation (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal). Michael Odelius: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Resources (equal); Supervision (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal). Ambar Banerjee: Data curation (equal); Methodology (equal); Software (equal); Validation (equal). Nina Ignatova: Formal analysis (equal); Investigation (equal); Methodology (equal). Alexander Föhlisch: Formal analysis (equal); Funding acquisition (equal); Methodology (equal); Resources (equal); Supervision (equal). Faris Gelmukhanov: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal). Victor Kimberg: Conceptualization (equal); Funding acquisition (equal); Project administration (equal); Resources (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal).
DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.
APPENDIX A: ALTERNATIVE WAYS TO COMPUTE XES USING TWO-TIME PROPAGATION METHOD
- One can apply a two-time forward propagator scheme [Fig. 7(a)] by permutation of and in the cross section (9),and computing the cross section as half-Fourier transform of the autocorrelation function σ(t, t1),To find the autocorrelation function, we need to solve the following time-dependent Schrödinger equations:(A2)(A3)
- Yet another approach to compute two-time integrals is to use two-dimensional (2D) time integration. Indeed, let us rewrite the cross section (9) asThe 2D integral does not depend on the order of integration, and we can permute the order with the corresponding modification of the integration limits, which results inwhere f(t, t1) = f*(t1, t) andIntroducing a new variable τ = t1 − t ≥ 0, we getwhere the wave packet Φfc(t, τ) is the solution of the Schrödinger equation in the final state with the initial condition of the wave packet in the core-ionized state |Ψc(t)⟩. This equation shows that the non-resonant XES cross section is defined by the interference of the two paths [see Fig. 7(b)].(A4)
Schematic illustration of various techniques using the two-time propagation method. (a) Two-time forward (hc and hf) propagators; (b) interference of two-time forward propagators.
Schematic illustration of various techniques using the two-time propagation method. (a) Two-time forward (hc and hf) propagators; (b) interference of two-time forward propagators.
Bending mode PECs for the states involved in Channel 1 (15) computed at the O–H distance taken at the equilibrium of the ground state of the neutral H2O. The vertical dashed line shows the equilibrium bending angle of the neutral molecule at the ground state.
Bending mode PECs for the states involved in Channel 1 (15) computed at the O–H distance taken at the equilibrium of the ground state of the neutral H2O. The vertical dashed line shows the equilibrium bending angle of the neutral molecule at the ground state.