Time-resolved x-ray/optical pump-probe simulations on N₂ molecules

The dynamics associated with the Auger decay of nitrogen molecules has been the subject of recent studies based on ultrafast x-ray sources.^1–5 In fact, N₂ was the subject of the first time-resolved x-ray pump–IR probe experiments at the Linac Coherent Light Source (LCLS).¹ After removal of a core electron via electron impact ionization or photoionization, N₂ decays primarily into dicationic nitrogen molecules $N22+$⁠.⁶ The molecular dications may end up as metastable $N22+$ in quasi-bound electronic states featuring a barrier towards molecular fragmentation or undergo further symmetric (⁠$N22+→N++N+$⁠) or asymmetric fragmentation (⁠$N22+→N+N2+$⁠).^5,6 Many studies on the Auger process in nitrogen molecules have been devoted to uncover these fragmentation dynamics focusing on Auger spectra,^7–12 the core-hole states from which the Auger process begins,^2,11,13–15 and the characterization of the quasi-bound and dissociative dicationic states populated by the Auger decay.^1,5,16 Most of the available works are based on measurements of properties of the fragmentation products, e.g., the kinetic energy release of the fragments or Auger electron energies, while the real-time nuclear dynamics triggered by x-ray absorption followed by Auger decay remains comparatively much less explored.

A time-resolved experimental study on the dynamics after the Auger process was conducted by Glownia et al.¹ using an x-ray pump–IR probe setup at LCLS. They found an enhancement of $N22+$ dissociation as a function of the pump-probe delay when the IR pulse was set up to operate after the x-ray pulse. Unfortunately, the timing uncertainty of a few hundred femtoseconds of the delay between both pulses¹ rendered more precise analysis of the data extremely challenging. Later on, those experimental data were re-analyzed based on a newly developed singular-value-decomposition scheme where the timing uncertainty of the original data could be largely compensated.⁴ A recent experimental and theoretical study was carried out to investigate the x-ray-induced wave-packet dynamics of $N22+$ after Auger decay.⁵ In this study, the kinetic-energy release of the ion fragments was measured. It was found to decrease as a function of the pump-probe delay time signaling the progress of dissociation in dissociative channels, whereas it was observed to remain constant for the quasi-bound channels of $N22+$⁠. A separate experimental study by Iwayama et al.¹⁷ employed coincidence measurements of the Auger electron and the ions following the Auger decay and indicated that the metastable states featuring quasi-bound potential energy curves and leading to $N22+$ formation are found for binding energies up to about 50 eV, in agreement with the window of final states considered in our previous theoretical study¹⁶ and in this work.

The quasi-bound states of $N22+$ and the vibrational dynamics on their potential energy curves have been, on their own, the subject of numerous theoretical and experimental investigations.^18–25 Early theoretical studies of $N22+$ by Wetmore and Boyd provided potential energy curves of the dissociative and quasi-bound states and discussed the tunneling dissociation properties of the latter.¹⁸ The theoretical study of Taylor and Patridge concluded definitively that the ground state of the $N22+$ system is $X1Σg+$⁠, ending a previously unresolved dispute.¹⁹ A recent theoretical investigation by Fantuzzi et al. focused on the electronic-structure origin of quasi-bound potential energy curves in homonuclear dicationic molecules including $N22+$ and found that the local minima of the dication states arise from a strong stabilization by polarization-aided interference effects.²⁵ The spectroscopic properties of vibrational states of quasi-bound potential energy curves in $N22+$ were the subject of theoretical and experimental investigations,^20–22 and a theoretical study by Jiang et al. used a Floquet picture to study the wave-packet dynamics in the quasi-bound potentials and their dissociation in the presence of a strong laser field.²³ 

In this study, the possibility to obtain dynamical information on the vibrational wave packets of $N22+$ after Auger decay through modulations in the yield of decay products in an x-ray pump–IR probe setup is investigated theoretically. We introduce a model that includes the relevant initial, intermediate, and final electronic states and the pump and probe pulses. The nuclear and electronic dynamics in this model are described quantum mechanically in a time-dependent picture. Based on this model, we calculate the variation of the $N22+$ yield as a function of the pump-probe delay time and show that the Fourier transform of this time-signal yields information on the vibrational wave packets created after the Auger decay process in the quasi-bound potential energy curves. Since the real-time x-ray-pump–IR-probe simulations are computationally expensive (cf. Sec. II A), a perturbative treatment of the probe pulse is introduced and described as well. Our results can be directly compared to the study by Fung et al.,⁴ where the delay-dependent $N22+$ yield obtained experimentally¹ was analyzed in detail. The potential energy curves and transition dipole moments used in this study, as well as the quantum-mechanical model of the pump and probe steps, have been described in detail in Ref. 16.

This paper is organized as follows. The theoretical and computational methods for full pump-probe simulations and the perturbative analysis are described in Secs. II A and II B, respectively. The results from the perturbative model are then discussed in Sec. III A, and the results from the pump-probe simulations are discussed in Sec. III B. Finally, the work is summarized in Sec. IV.

The theoretical model used to describe the Auger decay of $N2+$ and the subsequent wave-packet dynamics of the final states was presented elsewhere,¹⁶ including a detailed description of the potential energy surfaces, transition dipole moments, and decay width parameters. The relevant set of differential equations in atomic units is given here for the sake of completeness

In Eqs. (1)–(3), $|Φi(t)⟩, |Φd(t)⟩$⁠, and $|Φf(E,t;τ)⟩$ are time-dependent nuclear wave functions of the initial, intermediate, and final states, respectively, and E is the Auger electron energy. Note that the wave functions of the final states are Auger-energy-dependent,¹⁶ which is a consequence of total energy conservation during the Auger decay process, and they parametrically depend on the delay time of the probe pulse, τ. $Ĥ(i/d/f)=T̂nuc+V̂(i/d/f)$ is the nuclear Hamiltonian of the initial, intermediate, and final states, respectively, where $T̂nuc$ is the kinetic energy of the nuclei and $V̂(i/d/f)$ is the potential energy curve of the corresponding electronic states. $Âx(t), Γ̂, Ŵf$⁠, and $μ̂fg$ are the x-ray pulse envelope, the energy decay width of the intermediate state, the transition operators between the intermediate and the final states, and the transition dipole moment operators between final states f and g, respectively. ε(t − τ) is the electric field amplitude of the IR probe pulse centered at delay time τ and e_IR indicates its polarization direction. The potential energy curves of the final dicationic states considered in this study consist of the quasi-bound electronic states $X1Σg+, 21Σg+, 11Πu$⁠, and $11Σu+$ and the dissociative states $11Δg, 11Πg, 21Πg$⁠, and 3¹Π_g.¹⁶ One core-ionized intermediate state is considered in the model.¹⁶ 

The nuclear wave functions are described using two degrees of freedom, namely, the interatomic distance R and the angle θ between the molecular axis and the IR polarization direction.¹⁶ The wave-packet propagations were performed with the multi-configuration time-dependent Hartree method using the Heidelberg package.²⁶ The constant mean-field integrator^27–29 was used with an initial stepsize of 0.05 fs and an error threshold of 10⁻⁹. This gave an accurate propagation significantly faster in comparison to the more accurate variable mean-field propagator used in our previous work.¹⁶ The wave-packet propagations were carried out using shared-memory parallelization on 12-core Intel^® Xeon^® X5650 CPUs with 2.67 GHz and 96 GB of memory.

In this work, we are interested in the dynamical information that can be extracted from metastable $N22+$ after core-shell ionization and subsequent Auger decay, as a function of the delay time of an IR probe pulse. The IR pulse interacts with the system after the Auger decay and modifies the amount of unfragmented $N22+$ that reaches the detector. The yield of $N22+$ reaching the detector as a function of the pump-probe delay thus contains information on the metastable vibrational wave packets created after the Auger decay and is defined as

where $Θ̂(Rb−R)$ is 1 for R < R_b and 0; otherwise, R_b is the interatomic distance after which all potential energy curves become dissociative, and t_D is the time of arrival at the detector.

In propagating $|Φf(E,t;τ)⟩$⁠, which is later on used to evaluate Ω(τ), we use a Gaussian IR pulse of 800 nm central wavelength and 42.4 fs full-width at half-maximum of the intensity with a peak optical intensity I₀ of 5.62 × 10¹³ W cm⁻². The delay time between the pump and the probe pulse is varied from −100 fs up to 1500 fs in intervals of 0.5 fs each. This dense delay-time sampling is necessary to avoid aliasing effects when performing the Fourier transform due to the highly oscillatory $N22+$ yield with respect to the delay time (cf. Fig. 3). Each wave-packet propagation for a particular delay time and Auger energy is carried out up to a final simulation time t_F of 1700.0 fs. As for the pump pulse, we assume a Gaussian pulse of 5.0 fs duration.

The initial wave packets and the time-dependent basis function description to solve Eqs. (1)–(3) are taken to be the same as in Ref. 16. The propagation of $|Φf(E,t;τ)⟩$ and the evaluation of Ω(τ) using Eq. (4) involve 3201 delay times, and for each delay time, the Auger energy is discretized into 24 energy steps.¹⁶ This results in 76 824 separate wave-packet propagations that need to be carried out.

In the actual pump-probe experiments, the atomic ions N⁺ and the remaining molecular dication $N22+$ after the probe pulse are detected using an ion time-of-flight device with a time-of-flight on the order of μs.¹ Hence, the actual t_D is orders of magnitude larger than the final wave-packet propagation time t_F. This means that the possibility that the vibrational wave packets tunnel out of the quasi-bound potential energy curves after t_F has to be considered. Tunneling corrections were calculated using the method developed by Le Roy and Bernstein,^30,31 which is based on the WKB approximation. For a given quasi-bound vibrational energy E_ν in electronic state f, the tunneling probability per collision against the dissociation energy barrier is, in atomic units

where μ is the reduced mass of the system and V_f(R) is the potential energy curve of the corresponding electronic state. The integral in Eq. (5) is carried out between the second classical turning point and the exit point along the tunneling region below the potential energy barrier [R₂(E_ν), R₃(E_ν)]. On the other hand, the vibrational angular frequency is

where the integral is performed in the classical region [R₁(E_ν), R₂(E_ν)] between the first two classical turning points. The tunneling rate k_f(E_ν) is then given by

and the tunneling coefficient for vibrational state ν in electronic state f follows from first-order kinetics

where t_D was defined earlier as the time needed to reach the detector (the time-of-flight of the ion in question). For $N22+ tD=3.5 μs$ is considered after Ref. 1.

The tunneling projector for each electronic state f is now defined as

where ν runs over all quasi-bound vibrational states in the final electronic state f. The total $N22+$ yield as a function of delay time τ is now

where Ω(τ) in Eq. (4) has been redefined to include the tunneling correction. The vibrational states for the tunneling correction are computed by artificially continuing the potential energy curves horizontally after the corresponding dissociation barrier. $|ϕfν⟩$ in Eq. (9) corresponds to the ν-th vibrational eigenstate in such continued potentials. The corresponding eigenenergies are presented in Table I and their associated tunneling coefficients at the designated t_D are shown in Fig. 1. The final propagation time t_F is set 200 fs longer than the largest pump-probe delay time considered, such that all electronic populations have become stable after the IR pulse, and the newly promoted dissociative parts of the nuclear wave packets have been absorbed by a complex absorbing potential.¹⁶ 

As described in the previous section, the computational cost of the pump-probe signal is high due to the many propagations involved. Thus, we introduce a perturbative approximation of the effect of the IR probe on the $N22+$ yield to extract the dynamical information. After time t₀, which corresponds to a sufficiently long time after Auger decay such that all direct fragmentation has been completed (neglecting the tunneling effects) leaving only the quasi-bound $N22+$⁠, Eq. (3) is reduced to

In the interaction picture

Equation (11) reads

after some algebraic manipulations to simplify the expression. After going back to the Schrödinger picture, the first-order perturbation wave-function approximation to Eq. (13) is

where $|Φf(0)(E,t)⟩$ is the unperturbed wave-function solution, i.e., the field-free wave function. We note that $|Φ̃f(E,t;τ)⟩$ may contain contributions from all electronic states g after the perturbative interaction with the probe pulse.

To simplify Eq. (14), we introduce an IR pulse with a broad bandwidth that covers all possible transitions and a duration much shorter than the vibrational periods on the quasi-bound potential energy curves, namely, a δ-pulse ε(t) = ε₀ δ(t). For t > τ, this results in

Using the wave-function approximation in Eq. (15) and neglecting the tunneling correction, the approximated observable can be calculated as

where the double-numbered superscripts refer to the respective perturbation order of the bras and kets plugged into Eq. (4). Specifically, these terms are

As can be appreciated in Eq. (17), the (0, 0) term is delay-independent. At t > t₀, the only remaining wave-function amplitudes within the unfragmented region are from quasi-bound final states. Since $Ω̃(τ)$ concerns only the unfragmented wave packets at t_F, we assume that the $N22+$ yield from a pulse-induced transition remains unfragmented only when the transition is between quasi-bound states, i.e., $Θ̂(Rb−R)=1̂$ when f, g, and $g′$ are quasi-bound electronic states, and 0 otherwise. By applying this assumption to Eq. (18), the time-evolution operator can act to the left resulting in a transition dipole-coupled term between $|Φf(0)(E,τ)⟩$ and $|Φg(0)(E,τ)⟩$⁠, and its imaginary part cancels that of its complex conjugate, which also appears in the double sum in Eq. (18). Thus, the cross-term $Ω(0,1)(τ)+Ω(1,0)(τ)$ vanishes. By applying this assumption to Eq. (19), the time-evolution operators cancel each other and $Ω̃(τ)$ in Eq. (16) is simplified to

To calculate $Ω̃(τ)$⁠, field-free nuclear wave packets were propagated for τ up to 4000.0 fs using 0.5 fs intervals. The convolution integral with the x-ray pulse envelope and the Fourier transform of the resulting signal were then performed for τ > t₀, where we set t₀ = 200 fs.

We first present the results of the perturbative analysis. Figure 2(a) and Table II show the Fourier spectrum and peaks of $Ω̃(τ)$ from the perturbative model in Eq. (20). Moreover, $Ω̃(τ)$ can also be decomposed into the contributions from each quasi-bound state f

The Fourier peaks of $Ω̃f(τ)$ for each quasi-bound final electronic state are shown in Figs. 2(b)–2(e). Upon inspection of terms in Eq. (21), one can show that the extracted frequencies are, in fact, energy differences between vibrational states of $N22+$ (cf. Table I). To indicate the two vibrational levels of a quasi-bound electronic state involved in a given vibrational energy difference, we use the notation (ν₁, ν₂). For $g′=g$⁠, each term in Eq. (21) is a combination of expansion coefficients of the involved vibrational states, transition dipole assisted coupling, and oscillating functions corresponding to vibrational energy differences within a quasi-bound electronic state. Meanwhile, for $g′≠g$⁠, the vibrational states involved are from two different quasi-bound electronic states associated with the same final state f when interacting with the probe pulse.

As one can see from Eq. (21), for $f=X1Σg+$ and $21Σg+$⁠, the Fourier peaks shown in Figs. 2(b) and 2(c) originate from wave-packet oscillations in 1¹Π_u and $11Σu+$ states, which are coupled to the former states by transition dipole moments. Thus, these two spectra are similar in frequency, whereas the discrepancies in amplitude are from the transition dipole moments involved in those two states (cf. Fig. 3 of Ref. 16). The trend in peaks A1–A7 of Fig. 2(b) indicates that the population of middle-lying vibrational states in 1¹Π_u is higher than for the corresponding few lowest and highest vibrational states. Apart from these two electronic states, peaks A3, A5, and A7 also appear in the spectra of 1¹Π_u and $11Σu+$ with slightly different frequencies [cf. Figs. 2(d), 2(e), and Table II]. These multiple peak assignments explain why the corresponding peaks in the spectrum of the total yield in Fig. 2(a) are higher than neighboring ones. Peaks A11 and A12 shown in Fig. 2(c) indicate that $11Σu+$ has predominant population in the corresponding few lowest energy vibrational eigenstates. The $g′≠g$ terms between 1¹Π_u and $11Σu+$ states, which contribute to higher frequencies, are washed out by the finite duration of the x-ray pulse.

Meanwhile, the sources of the Fourier peaks for 1¹Π_u and $11Σu+$ states shown in Figs. 2(d) and 2(e) are the wave-packet oscillations in $X1Σg+$ and $21Σg+$ states, and thus, their frequencies are similar. The trend in peaks A3, A5, A7–A10, and A18 in Figs. 2(d) and 2(e) indicates that the higher vibrational states of $21Σg+$ are more populated than the corresponding lower states. On the other hand, peaks A13 and A14 show that the few lowest vibrational states of $X1Σg+$ are predominantly populated during the Auger decay. This combination of small population in high vibrational states of $X1Σg+$ and in low vibrational states of $21Σg+$ has the effect that the contributions of $g′≠g$ terms between these two states are unobserved.

The predominant population of the low-lying vibrational states of $X1Σg+$ and $11Σu+$ is to be expected given that the Franck-Condon regions of those two electronic states overlap with those of neutral and core-ionized N₂, whereas those of $21Σg+$ and 1¹Π_u are shifted to a longer interatomic distance.¹⁶ Consequently, $21Σg+$ and 1¹Π_u states contribute to the dissociation as reported by Iwayama et al.¹⁷ 

Since the contribution of the $g′≠g$ terms from pairs of distinct quasi-bound final electronic states is unobserved, Eq. (20) implies that the Fourier peaks in Fig. 2 come from the varying transition dipole moment functions within the unfragmented interatomic distance R_b, that is, within the well of the quasi-bound final electronic states. If the transition dipole moment functions were constants, then $Ω̃(τ)$ would be delay-independent and the Fourier peaks in Fig. 2 would be unobserved. Later, we will also find an effect of this kind in the pump-probe simulations in Sec. III B.

Despite this simple approach, the frequencies corresponding to the vibrational energy differences agree, within the experimental uncertainty, with the recent data analysis⁴ of the pump-probe experiments described in Ref. 1. Table II compares the experimental frequencies shown in Fig. 7 of the supplementary material of Ref. 4 to the frequencies corresponding to the vibrational energy differences in Fig. 2 for the frequency range 10–70 THz. Our results permit an assignment of each peak to the quasi-bound dicationic electronic states. However, since the frequencies are related to vibrational energy differences in each electronic state, which can be close to each other, there are multiple possible assignments to some of the experimental peaks. We note, however, that we consider only singlet dication electronic states,¹⁶ while it is possible that other electronic states (e.g., triplet states) may also contribute to the vibrational energy difference spectra.

Since we assumed in Eq. (15) the probe pulse to be represented by a delta function, more peaks appear in the Fourier spectrum of the $N22+$ yield than would be expected for a more spectrally selective probe pulse (as used in experiment). The pump-probe simulations in the following subsection confirm this expectation.

Figure 3 shows the Auger-energy integrated $N22+$ yield normalized by the total Auger yield

The $N22+$ yield with tunneling correction is calculated using Eq. (10), while the one without tunneling correction is calculated using Eq. (4) evaluated at t_D = t_F. The summations in those two equations can also be decomposed into the contributions from each quasi-bound state f to obtain electronically resolved yield. At negative delay times, the IR pulse does not interact with the final states. For these negative delays, the blue curve indicates that, of all Auger decayed molecules, about 23% remain bound as $N22+ 1700 fs$ after the x-ray pulse. At the detector (after 3.5 μs), the black curve, which includes the tunneling corrections, indicates a further decrease in the $N22+$ yield to about 17% of the Auger yield Ω_A. At positive delay times, the $N22+$ yield is further reduced by the interaction with the IR probe pulse which promotes amplitude from quasi-bound to dissociative potential energy curves. This accounts for the further reduction in the $N22+$ yield with and without tunneling correction by about 11% and 14% of the Auger yield Ω_A, respectively. A large portion of this IR-induced dissociation at positive delay times arises from population transfer from $11Σu+$ to 1¹Π_g. This is because the equilibrium distance of the former state matches the equilibrium distance of the N₂ ground state and the $N2+$ core-ionized intermediate state. Furthermore, the energy spacing between the $11Σu+$ and 1¹Π_g states is close to resonant.¹⁶ The slope of the yield reduction in Fig. 3 around the pulse overlap region is steeper than that reported experimentally [cf. Fig. 5(c) of Ref. 1]. The difference is caused by the relatively large pump-probe timing jitter in the experiment.¹ 

The signal of interest consists of the small modulations at delay times beyond 200 fs, for which the probe pulse interacts only with the metastable vibrational wave packets in the quasi-bound potential energy curves. The Fourier transform of the pump-probe delay signal, the $N22+$ yield, was obtained similar to that in the perturbative model. The raw signal is first convolved with the x-ray pulse envelope and the Fourier transform is taken for delay times τ ≥ 200.0 fs. Note that now the IR probe pulse is already included in each single time-propagation. The Fourier spectra of the $N22+$ yields with tunneling correction for each final quasi-bound electronic state as well as for the total yield are shown in Fig. 4. The peaks of the spectra are summarized in Table III.

Peaks of the total Fourier spectrum shown in Fig. 4 indicate that the unfragmented $N22+$ population reaching the detector is oscillating at those frequencies as the delay time advances. This oscillation in population is dominated by the frequency of peak B3 with the second major contribution from peak B5. There are also small contributions coming from peaks B2, B4, and B6, while peak B1 does not appear in the total spectrum. By resolving the detected $N22+$ yield electronically, i.e., decomposing the summation in Eq. (10) into the contribution from each quasi-bound state f

and obtaining spectra from each Ω_f(τ), we can see that peaks B2, B3, and B4 are from the $X1Σg+$ state, while peaks B5 and B6 originate from the 1¹Π_u and $11Σu+$ states, respectively. The electronically resolved spectra also show that the contributions from the $X1Σg+$ and 1¹Π_u states cancel each other out in peak B1 because the population transferred by the IR pulse among the two electronic states remains bound in this case, leading to no modulation of the total $N22+$ yield.

Comparison of the A peaks in Table II and B peaks in Table III shows that only beating frequencies originating from 1¹Π_u and $11Σu+$ states are observed in the spectra of the numerical pump-probe simulations, whereas $X1Σg+$ and $21Σg+$ states do not contribute to the B peaks. Since one of the major differences between the perturbative model giving rise to the A peaks and the numerical pump-probe simulations is the δ-pulse assumption [cf. Eq. (15)], these discrepancies between A and B peaks must also be related to this assumption. As has been noted earlier, since the use of this broad bandwidth limit of the IR pulse allows any possible transition, peaks from any vibrational energy differences with significant population and coupling function can be extracted. On the other hand, the use of an IR pulse of a certain duration, like the one used in the numerical pump-probe simulations, allows transitions only within its bandwidth, thus limiting the peaks that can be extracted from the $N22+$ yield, and only a few frequencies in the A peaks in Table II appear in the B peaks in Table III. Moreover, the purely dissociative potential energy curves (cf. Fig. 2 of Ref. 16) are of g symmetry. Thus, within one-photon transitions, the quasi-bound electronic states that contribute to the B peaks are mostly of u symmetry. Nonetheless, $X1Σg+$⁠, which is quasi-bound, contributes as well to modulations to the total $N22+$ yield through coupling to u symmetry potential energy surfaces, as seen for example for peak B3 below. To further understand these phenomena, in the following, we discuss the origins of B peaks from the $N22+$ yield before the tunneling dissociation takes place in an electronically and vibrationally resolved manner.

The total and electronically resolved spectra without tunneling correction shown in Fig. 5(a) can be obtained from the $N22+$ yields using the same equations as those with tunneling correction, i.e., Eqs. (10) and (23), respectively, by setting the tunneling coefficients ζ_f(E_ν) = 1. Furthermore, the vibrationally resolved spectra without tunneling correction shown in Figs. 5(b) and 5(c) can be obtained from the vibrationally resolved yields

evaluated at the final simulation time t_F, that is, prior to the tunneling dissociation.

As shown in Fig. 5, the most prominent peak B3 at 29.0 THz corresponds to the vibrational energy difference of 1¹Π_u (5, 6) (cf. Fig. 2 and Table II). Unlike the broad bandwidth limit in the perturbative model, the 800 nm Gaussian IR pulse with limited bandwidth in these pump-probe simulations couples these vibrational states of 1¹Π_u (ν = 5) and 1¹Π_u (ν = 6), which oscillate at a beating frequency of 29.0 THz, to the dissociative 1¹Π_g state more effectively at R close to the outer turning point of the quasi-bound potential well. There, the energy gap is close to resonant. Moreover, these two vibrational states are also resonantly coupled to the $X1Σg+ (ν=4)$ state by the probe pulse. From the previous discussion in Sec. III A, we know that the vibrational states 1¹Π_u (ν = 5) and 1¹Π_u (ν = 6) are initially populated. The incoming probe pulse at varying delay times modulates these two initially populated vibrational states and couples them to the dissociative 1¹Π_g state and to the $X1Σg+ (ν=4)$ state. Thus, these vibrational states involved exhibit the same peak B3 at 29.0 THz as shown in Figs. 5(b) and 5(c). Moreover, since the population transfers between the states involved differ in phase, the total modulation in the $N22+$ yield is smaller than the corresponding peak evaluated for each electronic state separately.

A fairly similar mechanism occurs as well for peak B2 shown in Fig. 5. In peak B2, the initially populated 1¹Π_u (ν = 6) state and a much less populated 1¹Π_u (ν = 7) state (as discussed in Sec. III A and shown in Fig. 2), which together oscillate at a frequency of 24.9 THz, are coupled by the probe pulse to the dissociative 1¹Π_g state and to $X1Σg+ (ν=5)$ state giving rise to the peak B2 in the spectra of the corresponding electronic and vibrational states. Due to low initial population of the 1¹Π_u (ν = 7) state, however, the amplitude of modulation of peak B2 is consequently lower than that of peak B3 (cf. Fig. 5).

The beating frequency indicated by peak B4 originates from the vibrational energy difference of 1¹Π_u (4, 5) [cf. Fig. 5(c)]. These two states, which are initially populated after the Auger decay, are coupled to $X1Σg+ (ν=4)$ by the probe pulse, as discussed previously, thus giving rise to the similar peak in the spectra of the corresponding electronic and vibrational states. Meanwhile, the vibrational states 1¹Π_u (ν = 3) and 1¹Π_u (ν = 4), which oscillate with the beating frequency indicated by peak B5, are coupled by the IR probe pulse to $X1Σg+ (ν=3)$ state [see Figs. 5(a) and 5(b)]. Apart from the IR-induced coupling to the quasi-bound vibrational states, these two pairs of vibrational states are also weakly coupled to the dissociative channels giving rise to the total spectrum shown in Fig. 5(a).

So far, we have seen how the oscillation mechanism originates from the shape of the potential energy curves of the states involved. For the $11Σu+$ state, however, the main source of the oscillation in the delay time signal is its transition dipole moment to the dissociative 1¹Π_g state. This is similar to that of the perturbation model in Sec. III A. Within the quasi-bound potential well of $11Σu+$ (roughly between 1.0 Å and 1.3 Å), the transition dipole moment function that couples $11Σu+$ to 1¹Π_g exhibits a large slope and, hence, drastically changes as a function of R (see Figs. 2(b) and 3(b) of Ref. 16). When the IR pulse arrives, the wave packet on this electronic state experiences a greater coupling function when $⟨R⟩$ of the wave packet is close to the inner turning point of the quasi-bound potential well than when its $⟨R⟩$ is close to the outer turning point of the well. Therefore, its coupling function is also oscillating as the wave packet oscillates in the quasi-bound well when the varying IR probe interacts with and couples it to the dissociative 1¹Π_g state. Since essentially only the three lowest vibrational states of $11Σu+$ are populated (cf. Table II), its population is oscillating with the frequency of the vibrational energy difference of $11Σu+ (0,1)$ manifested in peak B6 with a small side peak at 53.1 THz [cf. Fig. 5(b)]. This is the origin of the peak B6.

If one compares the Fourier spectra of the full simulations with and without tunneling correction in Figs. 4 and 5(a), respectively, one finds that peaks B2, B3, and B4 disappear in the 1¹Π_u final state and peak B5 changes in the tunneling-corrected spectra. During the flight towards the detector, the vibrational 1¹Π_u states for ν = 4 and higher undergo dissociation by tunneling (cf. Fig. 1). This has the effect that the corresponding modulations of the $N22+$ yield and, thus, the resulting Fourier peaks B2, B3, and B4 are absent in the 1¹Π_u (ν ≥ 4) final state in the spectra with tunneling correction. Peak B5 in the total spectrum, however, is more visible in the tunneling-corrected case. This is because one component in peak B5, namely, 1¹Π_u (ν = 4), undergoes tunneling dissociation during the flight to the detector and does not contribute to the signal, while 1¹Π_u (ν = 3) and $X1Σg+ (ν=3)$ remain unfragmented.

As expected, not all peaks among those in the perturbative model appear in the spectra of the actual pump-probe simulations. The discrepancies between the former and the latter spectra can be associated with the δ-pulse assumption made in the perturbative model, as well as the inclusion of tunneling predissociation in the pump-probe simulations. In the perturbative model, any dipole allowed transition contributes to and appears in the spectra, whereas in the full pump-probe simulations, the IR-induced transitions only contribute partially to the spectra. The discrepancies can also be associated with the nonperturbative nature of the pump-probe simulations. Although the A peaks originate from the leading-order of the perturbative model [cf. Eq. (20)], the nonperturbative nature of the full pump-probe simulations which generate the B peaks may very well exhibit features incompatible with the perturbative approximation, including the existence of peak B1 in the partial spectra.

A comparison of our result with the analysis by Fung et al. (cf. Fig. 4 of the main article and Fig. 7 of the supplementary material of Ref. 4), which is presented in Table II, shows that two peaks can be assigned from the simulations reported here, namely, peaks B3 (29.0 THz) and B4 (32.1 THz). We note, however, that in our simulations, we underestimated the pulse energy and did not use the same pulse parameters as the experiment in Ref. 1. We also note that the actual IR pulse intensity in the experiment of Ref. 1 may have been higher than the nominal value by a factor of two, possibly even by a factor of four.³² The appearance of peaks originating mostly from quasi-bound electronic states of u-type symmetry indicates that at peak optical intensity I₀, the IR-induced transitions in the pump-probe simulations arise from one-photon transitions and two-photon processes are mostly excluded. It cannot be ruled, however, that an increased peak optical intensity of the IR pulse beating frequencies involving g-type quasi-bound electronic states, which are partly observed in the data analysis, can be observed in the $N22+$ yield via two-photon transitions to dissociative g-symmetry states from the ground state electronic state of the dication.

We report real-time wave-packet simulations of x-ray pump-IR probe experiments on molecular nitrogen. As it is well known, several electronic states of $N22+$ are quasi-bound, namely, $X1Σg+, 21Σg+, 11Πu$⁠, and $11Σu+$⁠.¹⁶ The N₂ molecules are first core-ionized by the x-ray pulse and subsequently undergo Auger decay. The IR probe pulse then interrogates the nuclear wave-packet dynamics of $N22+$ in the quasi-bound potential energy curves of the dicationic states mentioned above. This work is motivated by the first time-resolved pump-probe experiments at LCLS¹ and the sophisticated data-analytical work that followed.⁴ We assume a Gaussian x-ray pump pulse of duration 5.0 fs and an IR probe pulse of 800 nm central wavelength and 42.4 fs duration. We also vary the peak optical intensity of the IR probe pulse to qualitatively determine its effect on the frequencies of the pump-probe spectrum.

On the basis of a perturbative treatment of the probe step, we first demonstrate that it is possible to extract dynamical information from the time-resolved $N22+$ yield modulated by the delay time of the IR probe pulse. This perturbative treatment yields all frequencies participating in the vibrational wave packets that can potentially be revealed by the probe at a fraction of the computational cost of the complete simulations. The frequencies extracted from this treatment are in good agreement with the frequencies extracted from the experimental data.⁴ 

The complete simulations include the probe pulse explicitly as well as the tunneling corrections of the populations of metastable vibrational levels. The simulated spectrum shows three main peaks at 24.9, 29.0, and 34.1 THz with several lower peaks around them. These peaks are mainly due to population transfer among 1¹Π_u and 1¹Π_g states induced by the IR probe pulse at R close to the inner barrier turning points. The assignment of these peaks to vibrational energy differences between specific pairs of vibrational states in a particular electronic state is reported as well. Two common frequencies (29.0 and 32.1 THz) are found by comparison of the data-analytical and theoretical results. From the results of the peak optical intensity variation, we have an indication that there may have been some inaccuracy between the measured peak optical intensity and the one experienced by the $N22+$ molecules in the experiment.¹ 

This work paves the way for a more sound understanding of the nuclear dynamics in molecules after Auger decay. It theoretically demonstrates that the delay-dependent yield of metastable species formed after the decay of core-ionized molecules and probed with a resonant laser can be used to reveal dynamical information on the time evolution of those decayed species. This approach is, in principle, not limited to small systems such as N₂ and may potentially be applied to molecules beyond diatomics after their initial decay and possible partial fragmentation. Further studies in this direction will be required.