Transient-Absorption Pump-Probe Spectra as Information-Rich Observables: Case Study of Fulvene

1. Introduction

Nonadiabatic processes taking place near crossing areas of potential energy surfaces (PESs) are ubiquitous in photochemistry and photobiology [1,2,3,4,5,6]. When two electronic states have the same multiplicity, conical intersections (CIs), i.e. degeneracy areas between adjacent adiabatic PESs, facilitate radiationless transitions [7] and population transfer [8]. CI-driven nonadiabatic processes occur, usually, on the femtosecond timescale [3,8]. It is thus instrumental but challenging to determine microscopic mechanisms and molecular motions governing the CI-driven wavepacket dynamics.

Nowadays, femtosecond spectroscopy is the main experimental source of information on ultrafast photoinduced processes [9,10,11,12,13,14,15,16,17], and transient-absorption (TA) pump-probe (PP) spectroscopy is one of the most widespread techniques [18,19,20,21,22,23,24]. TA PP spectroscopy delivers information on the wavepacket dynamics in the (coupled) ground and lower-lying excited electronic states. The information is encoded into, e.g., integral TA-PP spectra $I_{int}(\tau ,{\omega }_{pr})$[25] which are monitored as a function of the time time delay between the pump and probe pulses (so-called population time $\tau$) and probe frequency ${\omega }_{pr}$, but depend implicitly on the excitation frequency ${\omega }_{pu}$, durations of the PP pulses, and parameters describing molecular system under study. However, it is impossible to comprehensively characterize photoinduced dynamics at CIs by using TA PP spectra alone: theoretical support is essential for their interpretation. Thus, simulation of time-resolved spectra is important because it helps us to extract the essential information on nonadiabatic dynamics near CIs[26,27,28,29,30]. Consequently, considerable efforts were devoted to simulate nonlinear time-resolved signals within the framework of quantum dynamics[31,32], hierarchical equations of motion (HEOM)[33], mixed quantum-classical Liouville equations[34,35], surface hopping[36,37,38,39], and Ehrenfest dynamics[40]. These efforts were aimed to build a more straightforward connection between nonadiabatic dynamics driven by CIs and time-resolved signals[25,41,42,43,44,45,46,47,48].

With the advancement of computational facilities and theoretical algorithms, efforts of many groups have been shifted to the construction of on-the-fly protocols for simulations of spectroscopic signals for realistic molecular systems [45,48,49,50,51,52,53,54,55,56,57]. These efforts open the possibility to microscopically understand ultrafast time-resolved spectra of polyatomic systems and to get insight into nonadiabatic dynamics involving large-amplitude molecular motions.

Usually, the on-the-fly simulations are conducted within the trajectory-based dynamics formulism. Among them, the quasiclassical or semiclassical dynamics methods based on the mapping Hamiltonian [58,59,60] received considerable attentions. Different versions of the mapping approach were implemented to study nonadiabatic dynamics at CIs in recent decades[61,62,63,64,65,66,67,68,69,70,71,72]. Within the mapping framework, different theoretical approaches were proposed to simulate different types of nonlinear spectroscopic signals [73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89]. It is to be noted that the direct simulation of spectroscopic third-order response functions by the mapping approaches is quite tedious even for model systems [46]. On the other hand, the direct inclusion of laser fields at different time delays in the propagation is also computationally costly [90,91,92,93]. In simple words, more rigorous methods generally have to face larger computational costs. It is thus worthwhile to invoke some alternative, practical ideas for the simulation of time-resolved spectra[48].

The doorway-window (DW) approximation[94,95,96,97,98] accounting for finite duration and arbitrary spectral shapes of the pump and probe pulses provides a practical approach for the evaluation of spectroscopic signals. This approximation can conveniently be interfaced with any ab initio trajectory method: As the input, the DW protocol requires exclusively electronic energies and transition dipole moments (TDMs) along trajectories. So far, the ab initio DW methodology has been combined with two commonly used simulation techniques: surface-hopping (Tully’s fewest switches [99,100,101], Landau-Zenner [102,103], and machine-learning-enhanced Landau-Zenner [104]) and Ehrenfest [105].

Recently, a combination of symmetrical quasi-classical/Meyer-Miller-Stock-Thoss (SQC/MMST) dynamics with the DW approximation has been proposed for the simulation of TA PP signals [106]. Clearly, the SQC/MMST approach permits one to calculate nonadiabatic wavepacket dynamics with (usually) sufficient accuracy and at affordable computational costs [84,107,108,109,110,111,112]. Furthermore, the mapping approach may have advantages over the surface hopping and Ehrenfest counterparts in reproducing electronic coherences [64,113], which can become essential for adequate simulation of TA PP and other nonlinear spectroscopic signals. In addition, the mapping approach can be reformulated in the adiabatic representation [114] and conveniently used in the on-the-fly simulations to describe realistic molecular systems [114,115,116,117].

In this work, we took a practical idea to simulate TA PP signals by ab initio methods by interfacing the DW approximation and the SQC/MMST dynamics (DW-SQC/MMST). It has both theoretical and methodological dimensions. On the methodological side, we further extend the application base of the on-the-fly DW simulation protocol [102].

On the theoretical side, we explore the possibility of abstracting the key information about the nuclear motions involved in the CI-driven nonadiabatic dynamics by combined on-the-fly simulation of the traditional dynamical observables (electronic populations, bond lengths, and bond angles) and TA PP spectra.

The molecule of fulvene was chosen as a typical CI-driven molecular system, which is frequently used as a real-world representative of "Tully’s third model" [118]. The photoinduced dynamics in fulvene is fairly well understood at the ab initio level [119,120,121,122,123,124,125,126,127,128,128]. This allows us to concentrate on squeezing additional information from the TA PP signals. It essential that surface hopping methods do not give very precise description of fulvene’s population dynamics [119,120,121,122], while the SQC/MMST method with triangular window function provides results which are consistent with ab initio multiple spawning (AIMS) data [122]. This lends support to the choosing of the SQC/MMST method for simulations of TA PP signals of fulvene. Our results show that the TA PP signals directly reflect the key molecular motions in the nonadiabatic dynamics of fulvene. This work thus provides useful ideas on how to connect ultrafast spectroscopic signals with the essential features in nonadiabatic dynamics.

2. Theoretical Methods and Computational Details

2.1. SQC/MMST Approach

Mapping models, in which an F-level quantum system is mapped to a system with F coupled continuous degrees of freedom (DoFs), provide practical ideas for the construction of classical-like dynamics approaches for the simulation of quantum evolutions of realistic large system[64,67]. In the MMST model[58,59,60], a set of discrete quantum states is mapped to a group of coupled classical oscillators. With the SQC protocol, the initial sampling and the final assignment of the electronic DoFs is performed through the symmetrical triangle window function [129,130]. More rigorous derivations can be found in the previous works [59,110,131,132]. By combining the MMST mapping model with the SQC protocols, the SQC/MMST approach allows us to simulate the nonadiabatic dynamics [107,108] with acceptable accuracy and affordable computational costs[109,111,112]. Since the SQC/MMST method is easy to reformulate in the adiabatic representation,[114] it can be conveniently combined with the ab into quantum chemistry software to describe the nonadiabatic dynamics in realistic molecular systems[115,116,117] by using the trajectory-adjusted zero-point energy (ZPE) correction $\gamma$ [114,133].

The classical MMST Hamiltonian in the adiabatic representation reads

$$\begin{array}{cc}\hfill H& ={\displaystyle \frac{1}{2\mathbf{M}}}{\mathbf{P}}_{kin}^2+\sum _i^F\left[{\displaystyle \frac{1}{2}}{\left(x_i\right)}^2+{\displaystyle \frac{1}{2}}{\left(p_i\right)}^2-\gamma \right]V_i\left(\mathbf{R}\right).\hfill \end{array}$$

(1)

Here the kinematic momentum ${\mathbf{P}}_{kin}$ is defined as

$$\begin{array}{cc}\hfill {\mathbf{P}}_{kin}& =\mathbf{P}+\sum _{ij}^F\left(x_i{p}_j\right){\mathbf{d}}_{ij}\left(\mathbf{R}\right)\hfill \end{array}$$

(2)

where $\mathbf{R}$ and $\mathbf{P}$ are nuclear coordinates and momenta; $x_i$ and $p_i$ are coordinates and momenta of the mapped electronic states $i=1,...,F$; F is the total number of these states; $V_i\left(\mathbf{R}\right)$ is the potential energy of the electronic state i; ${\mathbf{d}}_{ij}$ is the first-order non-adiabatic coupling vector.

The effective potential

$$\begin{array}{cc}\hfill V_{\mathrm{eff}}\left(\mathbf{R}\right))& =\sum _i^F{\displaystyle \frac{1}{2}}(x_i^2+p_i^2-2\gamma )V_i\left(\mathbf{R}\right)\hfill \end{array}$$

(3)

is customarily rearranged into the symmetrized form

$$\begin{array}{cc}\hfill V_{\mathrm{eff}}\left(\mathbf{R}\right))& ={\displaystyle \frac{1}{F}}\sum _i^F{V}_i\left(\mathbf{R}\right)+{\displaystyle \frac{1}{F}}\sum _{ij}^F{\displaystyle \frac{1}{4}}(x_i^2+p_i^2-x_j^2-p_j^2)(V_i\left(\mathbf{R}\right)-V_j\left(\mathbf{R}\right)),\hfill \end{array}$$

(4)

which enhances numerical stability.

The classical equations of motion generated by the Hamiltonian of Eq. (1) assume the form [114]:

$$\begin{array}{cc}\hfill {\dot{x}}_i& =p_i{\displaystyle \frac{1}{F}}\sum _j^F(V_i\left(\mathbf{R}\right)-V_j\left(\mathbf{R}\right))+\sum _j^F\left(x_j{\mathbf{d}}_{ji}\left(\mathbf{R}\right)\right)·{\displaystyle \frac{{\mathbf{P}}_{kin}}{\mathbf{M}}}\hfill \\ \hfill {\dot{p}}_i& =-x_i{\displaystyle \frac{1}{F}}\sum _j^F(V_i\left(\mathbf{R}\right)-V_j\left(\mathbf{R}\right))+\sum _j^F\left(p_j{\mathbf{d}}_{ji}\left(\mathbf{R}\right)\right)·{\displaystyle \frac{{\mathbf{P}}_{kin}}{\mathbf{M}}}\hfill \\ \hfill \dot{\mathbf{R}}& ={\displaystyle \frac{{\mathbf{P}}_{kin}}{\mathbf{M}}}\hfill \\ \hfill {\dot{\mathbf{P}}}_{kin}& =-{\displaystyle \frac{\partial V_{\mathrm{eff}}\left(\mathbf{R}\right)}{\partial \mathbf{R}}}+{\displaystyle \frac{1}{2}}\sum _{ij}^F(p_i{p}_j+x_i{x}_j)(V_j\left(\mathbf{R}\right)-V_i\left(\mathbf{R}\right)){\mathbf{d}}_{ji}\left(\mathbf{R}\right).\hfill \end{array}$$

(5)

In the SQC/MMST simulations, the window function (which should not be confused with the window function of the DW representation) is employed for the initial sampling and final assignment of discrete electronic states [108,129]. The so-called triangle window function is defined as [130]

$$\begin{array}{c}\hfill W_i(\mathbf{n}=n_1,...n_i,...n_f)=w_1\left(n_i\right)\prod _{i\ne j}^F{w}_0(n_i,n_j),\end{array}$$

(6)

where

$$\begin{array}{c}\hfill w_1\left(n_i\right)=\left\{\begin{array}{cc}{(2-\gamma -n_i)}^{2-F}\hfill & \mathrm{for}1-\gamma <n_i<2-\gamma \hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array}\right.\end{array}$$

(7)

and

$$\begin{array}{c}\hfill w_0(n_i,n_j)=\left\{\begin{array}{cc}1\hfill & \mathrm{for}{n}_j<2-2\gamma -n_i\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array}\right.\end{array}$$

(8)

In other words, $w_i=1$ if the system is in the state i, and $w_i=0$ otherwise.

The action-angle sampling method produces initial electronic positions and momenta

$$\begin{array}{c}\hfill x_i^0=\sqrt{2(n_i^0+\gamma )}cos\theta ,\\ \hfill p_i^0=\sqrt{2(n_i^0+\gamma )}sin\theta ,\end{array}$$

(9)

where $n_i^0$ is obtained by using the triangle window function and $\theta \in [-\pi ,\pi ]$. For the final assignment, the triangle window function recasts the action variables $n_i$ into integers 0 or 1,

$$\begin{array}{c}\hfill n_i={\displaystyle \frac{1}{2}}{\left(x_i\right)}^2+{\displaystyle \frac{1}{2}}{\left(p_i\right)}^2-\gamma .\end{array}$$

(10)

Recently, Miller and Cotton proposed a trajectory-adjusted ZPE scheme, wherein each trajectory incorporates its specific ZPE correction $\gamma$[134]. Once the electronic coordinates and momenta are initially sampled, the ZPE correction is defined as

$${\gamma }_i={\displaystyle \frac{1}{2}}{\left(x_i^0\right)}^2+{\displaystyle \frac{1}{2}}{\left(p_i^0\right)}^2-{\delta }_{ij},$$

(11)

where j is the adiabatic electronic state in the initial sampling. Naturally, the value for ${\gamma }_i$ varies for each trajectory. Using this trajectory-specific ${\gamma }_i$, the action variable of the electronic DoF is defined as

$$n_i={\displaystyle \frac{1}{2}}{\left(x_i\right)}^2+{\displaystyle \frac{1}{2}}{\left(p_i\right)}^2-{\gamma }_i$$

(12)

during the trajectory propagation [134].

Meanwhile, the MMST Hamiltonian in the adiabatic representation is expressed in terms of the trajectory-adjusted ${\gamma }_i$ as follows:

$$H={\displaystyle \frac{1}{2\mathbf{M}}}{\mathbf{P}}_{kin}^2+\sum _i^F\left[{\displaystyle \frac{1}{2}}{\left(x_i\right)}^2+{\displaystyle \frac{1}{2}}{\left(p_i\right)}^2-{\gamma }_i\right]E_i\left(\mathbf{R}\right).$$

(13)

More details on the trajectory-adjusted ZPE approach can be found in Ref. [134]. Recently, further interpretations and physical insights into the adjustable ZPE correction $\gamma$ were presented[66,110,132,135,136].

2.2. DW Representation of TA PP Signals

The integral TA PP signal $I_{int}(\tau ,{\omega }_{pr})$ can be obtained as a sum of non-rephasing ($\alpha =\mathrm{NR}$, ${\xi }_{\alpha }=1$) and rephasing ($\alpha =\mathrm{R}$, ${\xi }_{\alpha }=-1$) terms, each of which consists of the ground-state bleach (GSB), stimulated emission (SE), and excited-state absorption (ESA) contributions. $I_{int}(\tau ,{\omega }_{pr})$ can be evaluated in terms of the third-order response functions $R_{\alpha k}(t_3,t_2,t_1)$ as follows [25,42,47,48]:

$$I_{int}(\tau ,{\omega }_{pr})\sim \mathrm{Re}\sum _{\alpha =\mathrm{R},\mathrm{NR}}\sum _{k=0,\mathrm{I},\mathrm{II}}{\int }_{-\infty }^{\infty }dt{\int }_0^{\infty }d{t}_3{\int }_0^{\infty }d{t}_2{\int }_0^{\infty }d{t}_1\times$$

(14)

$$E_{pu}(t+\tau -t_3-t_2-t_1)E_{pu}(t+\tau -t_3-t_2)E_{pr}(t-t_3)E_{pr}\left(t\right)e^{i{\xi }_{\alpha }{\omega }_{pu}{t}_1}{e}^{i{\omega }_{pr}{t}_3}{R}_{\alpha k}(t_3,t_2,t_1).$$

Here $k=0,\mathrm{I},\mathrm{II}$ correspond to the GSB, SE, and ESA contributions, $E_j\left(t\right)$ and ${\omega }_j$ are the dimensionless envelopes of the pump ($j=pu$) and probe ($j=pr$) pulses, and $\tau$ is the time delay between them.

To evaluate the response functions $R_{\alpha k}(t_3,t_2,t_1)$, it is convenient to subdivide all molecular electronic states into three manifolds 0, I and II which correspond to the ground, lower-lying, and upper-lying states. Then the molecular Hamiltonian of Eq. (1) assumes the block-diagonal form,

$$\begin{array}{cc}\hfill H=& \left(\begin{array}{ccc}{H}_0& H_0\mathrm{I}& 0\\ H_{\mathrm{I}}0& H_{\mathrm{I}}& 0\\ 0& 0& H_{\mathrm{II}}\end{array}\right).\hfill \end{array}$$

(15)

Here the Hamiltonians $H_0$, $H_{\mathrm{I}}$ and $H_{\mathrm{II}}$ – each of which can be represented as a sum of kinetic and potential parts – describe the (nonadiabatic) dynamic in the corresponding manifolds. $H_0\mathrm{I}=H_{\mathrm{I}}^{†}0$ are the inter-manifold coupling terms, which are relevant to the nonadiabatic couplings in the adiabatic representation. The rising and lowering TDM operators can also be recast in the block-diagonal form

$$\begin{array}{c}\hfill {\mu }^{\uparrow }=\left(\begin{array}{ccc}0& {\mu }_0,\mathrm{I}& 0\\ 0& 0& {\mu }_{\mathrm{I}},\mathrm{II}\\ 0& 0& 0\end{array}\right),{\mu }^{\downarrow }=\left(\begin{array}{ccc}0& 0& 0\\ {\mu }_{\mathrm{I}},0& 0& 0\\ 0& {\mu }_{\mathrm{II}},\mathrm{I}& 0\end{array}\right).\end{array}$$

(16)

As explained in Refs. [100,102], the TA PP signal of Eq. (14) can be represented in the form suitable for on-the-fly simulations by performing a series of approximations. First, the quantum DW approximation is applied, according to which

(i) The pump and probe pulses are well temporally separated. Hence the time delay $\tau$ between the pulses is (much) longer than the pump (${\tau }_{pu}$) and probe (${\tau }_{pr}$) pulse durations.

(ii) The pump and probe pulses are short on the nuclear dynamics timescale.

With these assumptions, the integral TA PP signal takes the DW form [100,102]

$$\begin{array}{c}\hfill I_{int}(\tau ,{\omega }_{pr})={\omega }_{pr}\mathrm{Tr}\left[D\left({\omega }_{pu}\right)\left(e^{i{H}_0\tau }{W}_0\left({\omega }_{pr}\right)e^{-i{H}_0\tau }+e^{i{H}_{\mathrm{I}}\tau }(W_{\mathrm{I}}\left({\omega }_{pr}\right)-W_{\mathrm{II}}\left({\omega }_{pr}\right))e^{-i{H}_{\mathrm{I}}\tau }\right)\right],\end{array}$$

(17)

where

$$\begin{array}{cc}\hfill D\left({\omega }_{pu}\right)& ={\int }_{-\infty }^{\infty }d{t}_2^{\prime }{\int }_0^{\infty }d{t}_1{E}_{pu}\left(t_2^{\prime }\right)E_{pu}(t_2^{\prime }-t_1)e^{i{\omega }_{pu}{t}_1}{e}^{-i{H}_{\mathrm{I}}{t}_1}{\mu }_{\mathrm{I}},0{\rho }_0{e}^{i{H}_0{t}_1}{\mu }_0,\mathrm{I}+h.c.\hfill \end{array}$$

(18)

is the doorway operator, ${\rho }_0$ is the initial nuclear distribution in manifold 0, which is determined by either the thermal equilibrium or a particular vibrational level, and

$$\begin{array}{cc}\hfill W_0\left({\omega }_{pr}\right)& ={\int }_{-\infty }^{\infty }d{t}^{\prime }{\int }_0^{\infty }d{t}_3{E}_{pr}\left(t^{\prime }\right)E_{pr}(t^{\prime }+t_3)e^{i{\omega }_{pr}{t}_3}{e}^{i{H}_0{t}_3}{\mu }_0,\mathrm{I}{e}^{-i{H}_{\mathrm{I}}{t}_3}{\mu }_{\mathrm{I}},0+h.c.\hfill \\ \hfill W_{\mathrm{I}}\left({\omega }_{pr}\right)& ={\int }_{-\infty }^{\infty }d{t}^{\prime }{\int }_0^{\infty }d{t}_3{E}_{pr}\left(t^{\prime }\right)E_{pr}(t^{\prime }+t_3)e^{i{\omega }_{pr}{t}_3}{\mu }_{\mathrm{I}},0e^{i{H}_0{t}_3}{\mu }_0,\mathrm{I}{e}^{-i{H}_{\mathrm{I}}{t}_3}+h.c.\hfill \\ \hfill W_{\mathrm{II}}\left({\omega }_{pr}\right)& ={\int }_{-\infty }^{\infty }d{t}^{\prime }{\int }_0^{\infty }d{t}_3{E}_{pr}\left(t^{\prime }\right)E_{pr}(t^{\prime }+t_3)e^{i{\omega }_{pr}{t}_3}{\mu }_{\mathrm{I}},\mathrm{II}{e}^{-i{H}_{\mathrm{II}}{t}_3}{\mu }_{\mathrm{II}},\mathrm{I}{e}^{i{H}_{\mathrm{I}}{t}_3}+h.c.\hfill \end{array}$$

(19)

are the window operators.

Then, a series of classical DW approximations is applied.

(iii) The DW operators $D\left({\omega }_{pu}\right)$ and $W_k\left({\omega }_{pr}\right)$ become functions of the nuclear coordinates $\mathit{R}$ and momenta $\mathit{P}$, $D({\omega }_{pu},\mathit{R},\mathit{P})$ and $W_k({\omega }_{pr},\mathit{R},\mathit{P})$.

(iv) Radiative transitions are treated within the classical Condon approximation [74,137].

(v) The initial vibrational distribution ${\rho }_0$ in the doorway function (18) is replaced by the Wigner distribution ${\rho }_g^{Wig}(\mathit{R},\mathit{P})$ [138].

(vi) The trace $\mathrm{Tr}[...]$ in Eq. (17) is evaluated by Monte Carlo averaging $〈...〉$ over classical initial conditions.

(vii) The Heisenberg propagators in Eq. (17) are replaced by the evolution along (quasi) classical trajectories $\mathit{R}\left(T\right),\mathit{P}\left(T\right)$.

For converting these approximations into operational expressions and simulation protocols, we denote: electronic states in manifolds 0, I, II as $|g\rangle$, $|e\rangle$, $|f\rangle$; potential energy functions in these states as $V_g\left(\mathit{R}\right)$, $V_e\left(\mathit{R}\right)$, $V_f\left(\mathit{R}\right)$; transition frequencies as $U_{eg}\left(\mathit{R}\right)=V_e\left(\mathit{R}\right)-V_g\left(\mathit{R}\right)$ and $U_{fe}\left(\mathit{R}\right)=V_f\left(\mathit{R}\right)-V_e\left(\mathit{R}\right)$; matrix elements of TDMs as ${\mu }_{ge}\left(\mathit{R}\right)=\langle g|{\mu }_{0,\mathrm{I}}|e\rangle$ and ${\mu }_{fe}\left(\mathit{R}\right)=\langle f|{\mu }_{\mathrm{II},\mathrm{I}}|e\rangle$. Then approximations (iii) and (iv) yield [74,137]

$$\langle g|e^{i{H}_0{t}_3}{\mu }_{0,\mathrm{I}}{e}^{-i{H}_{\mathrm{I}}{t}_3}|e\rangle \approx e^{i{U}_{ge}\left(\mathit{R}\right)t_3}{\mu }_{ge}\left(\mathit{R}\right),$$

(20)

$$\langle f|e^{i{H}_{\mathrm{II}}{t}_3}{\mu }_{\mathrm{II},\mathrm{I}}{e}^{-i{H}_{\mathrm{I}}{t}_3}|e\rangle \approx e^{i{U}_{fe}\left(\mathit{R}\right)t_3}{\mu }_{fe}\left(\mathit{R}\right)$$

(21)

and we arrive at the final quasiclassical DW expression for the integral TA PP signal:

$$\begin{array}{cc}\hfill I_{int}(\tau ,{\omega }_{pr})& =I_{int}^{GSB}(\tau ,{\omega }_{pr})+I_{int}^{SE}(\tau ,{\omega }_{pr})+I_{int}^{ESA}(\tau ,{\omega }_{pr}),\hfill \\ \hfill I_{int}^{GSB}(\tau ,{\omega }_{pr})& ={\omega }_{pr}\sum _e\int d{\mathbf{R}}_{g}d{\mathbf{P}}_{g}D({\omega }_{pu},{\mathbf{R}}_g,{\mathbf{P}}_g)W_0({\omega }_{pu},{\mathbf{R}}_g\left(\tau \right),{\mathbf{P}}_g\left(\tau \right)),\hfill \\ \hfill I_{int}^{SE}(\tau ,{\omega }_{pr})& ={\omega }_{pr}\sum _e\int d{\mathbf{R}}_{g}d{\mathbf{P}}_{g}D({\omega }_{pu},{\mathbf{R}}_g,{\mathbf{P}}_g)W_{\mathrm{I}}({\omega }_{pu},{\mathbf{R}}_e\left(\tau \right),{\mathbf{P}}_e\left(\tau \right)),\hfill \\ \hfill I_{int}^{ESA}(\tau ,{\omega }_{pr})& =-{\omega }_{pr}\sum _e\int d{\mathbf{R}}_{g}d{\mathbf{P}}_{g}D({\omega }_{pu},{\mathbf{R}}_g,{\mathbf{P}}_g)W_{\mathrm{II}}({\omega }_{pu},{\mathbf{R}}_e\left(\tau \right),{\mathbf{P}}_e\left(\tau \right)).\hfill \end{array}$$

(22)

Here ${\mathbf{R}}_g$ and ${\mathbf{P}}_g$ represent the initial nuclear coordinates and momenta in the electronic ground state sampled according to the Wigner distribution ${\rho }_g^{Wig}({\mathbf{R}}_g,{\mathbf{P}}_g)$, ${\mathbf{R}}_g\left(\tau \right)$ and ${\mathbf{P}}_g\left(\tau \right)$ denote the nuclear coordinates and momenta propagated in the electronic ground state up to $t=\tau$, and ${\mathbf{R}}_e\left(\tau \right)$ and ${\mathbf{P}}_e\left(\tau \right)$ denote the nuclear coordinates and momenta propagated up to $t=\tau$ in manifold I of the lower-lying excited electronic states. The quasiclassical doorway function reads

$$\begin{array}{cc}\hfill D({\omega }_{pu},{\mathbf{R}}_g,{\mathbf{P}}_g)& =|{\mu }_{ge}\left({\mathbf{R}}_g\right){|}^2{E}_{pu}^2({\omega }_{pu}-U_{eg}\left({\mathbf{R}}_g\right)){\rho }_g^{Wig}({\mathbf{R}}_g,{\mathbf{P}}_g),\hfill \end{array}$$

(23)

and the quasiclassical window functions are defined as

$$\begin{array}{cc}\hfill W_0({\omega }_{pu},{\mathbf{R}}_g\left(\tau \right),{\mathbf{P}}_g\left(\tau \right))& =\sum _e{\left|{\mu }_{ge}\left({\mathbf{R}}_g\left(\tau \right)\right)\right|}^2{E}_{pr}^2({\omega }_{pr}-U_{eg}\left({\mathbf{R}}_g\left(\tau \right)\right)),\hfill \\ \hfill W_{\mathrm{I}}({\omega }_{pu},{\mathbf{R}}_e\left(\tau \right),{\mathbf{P}}_e\left(\tau \right))& =|{\mu }_{ge\left(\tau \right)}\left({\mathbf{R}}_e\left(\tau \right)\right){|}^2{E}_{pr}^2({\omega }_{pr}-U_{e\left(\tau \right)g}\left({\mathbf{R}}_e\left(\tau \right)\right)),\hfill \\ \hfill W_{\mathrm{II}}({\omega }_{pu},{\mathbf{R}}_e\left(\tau \right),{\mathbf{P}}_e\left(\tau \right))& =\sum _f{\left|{\mu }_{e\left(\tau \right)f}\left({\mathbf{R}}_e\left(\tau \right)\right)\right|}^2{E}_{pr}^2({\omega }_{pr}-U_{fe\left(\tau \right)}\left({\mathbf{R}}_e\left(\tau \right)\right)).\hfill \end{array}$$

(24)

Here $E_{pu}\left(\omega \right)$ and $E_{pr}\left(\omega \right)$ are the Fourier transforms of $E_{pu}\left(t\right)$ and $E_{pr}\left(t\right)$, and the notion $e\left(\tau \right)$ means that a trajectory initiated in an excited state e may jump into another electronic state.

The above derivations have been made under the assumption that nonadiabatic coupling operators $H_0\mathrm{I}=H_{\mathrm{I}}^{†}0$ can be neglected on the timescale of interest. Once the internal conversion $e\left(\tau \right)\to g$ is allowed, Eq. (17) remains valid, but the window functions $W_k$ defined per Eq. (24) have to be replaced by the new window functions ${\overline{W}}_k$ which are defined in terms of the original window functions as follows [100]:

$${\overline{W}}_0=W_0,$$

(25)

$${\overline{W}}_{\mathrm{I}}=\left[\begin{array}{c}{W}_{\mathrm{I}},\mathrm{if}\mathrm{trajectory}\mathrm{stays}\mathrm{within}\mathrm{I}\\ -W_0,\mathrm{if}\mathrm{trajectory}\mathrm{jumps}\mathrm{from}\mathrm{I}\mathrm{to}0,\end{array}\right.$$

(26)

$${\overline{W}}_{\mathrm{II}}=\left[\begin{array}{c}{W}_{\mathrm{II}},\mathrm{if}\mathrm{trajectory}\mathrm{stays}\mathrm{within}\mathrm{I}\\ 0,\mathrm{if}\mathrm{trajectory}\mathrm{jumps}\mathrm{from}\mathrm{I}\mathrm{to}0.\end{array}\right.$$

(27)

If a trajectory jumps from manifold I back to the ground state, then $e\left(\tau \right)$ changes to g in Eq. (26). Owing to the $e\left(\tau \right)\to g$ internal conversion, two GSB contributions arise: the cold and the hot. The cold contribution is the conventional GSB signal which reveals the nuclear wavepacket on the electronic ground state. It is described by the window function $W_0$ of Eq. (25). The hot contribution is the internal-conversion-induced GSB signal, that reveals the manifold-I trajectory which jumps to the electronic ground state after the $e\left(\tau \right)\to g$ internal conversion. This trajectory contributes to the GSB signal with a minus sign, $-W_0$ in Eq. (26) [100]. In the SQC/MMST simulations, the DW functions are straightforwardly calculated by "binning" trajectories to the windows corresponding to the currently occupied electronic states [130].

2.3. Computational Details

The optimized structures of fulvene are shown in Figure 1. The $S_0$ minimum was optimized at the B3LYP/6-31G(d) level by using the Gaussian 16 package [139]. In the SQC/MMST simulations, two state-averaged CASSCF(6,6) method with the basis set 6-31G(d) was employed in the MOLPRO2022 package [140]. Initial nuclear coordinates and momenta were sampled according to the Wigner distribution of the lowest vibrational level on the ground electronic state[141]. Initial sampling and final assignments of electronic states were performed by using the symmetrical triangle window function of Eq. (6), while coordinates and momenta of mapping electronic degrees of freedom were generated according to the action-angle distribution of Eq. (9). The time steps were fixed at 0.2 fs and 0.002 fs for nuclear and mapping electronic motions, respectively. All dynamics calculations were carried out with the JADE package [116,141].

For fulvene, manifold I contains a single state $S_1$, while manifold II consists of four states $S_2$-$S_5$. 100 trajectories were initiated on the $S_1$ state and propagated for 600 fs in manifold I to simulate the hot GSB, SE, and ESA signals. Another 100 trajectories were initiated on the $S_0$ state and propagated for 600 fs to simulate the cold GSB signal. To evaluate ESA signals, snapshots from every trajectory were taken to obtain vertical excitation energies and TDMs for the states of manifold II.

Gaussian pump and probe pulses were chosen, $E_a\left(t\right)=exp\{-{(t/{\tau }_a)}^2\}$ and $E_a\left(\omega \right)=exp\{-{\left(\omega {\tau }_a\right)}^2/4\}$, with ${\tau }_a=2$ fs. This gives the pulse bandwidth of 1.04 eV (full width at half maximum).

3. Results and Discussion

To setup the stage, we begin with the populations of the $S_0$ and $S_1$ states of fulvene which are shown in Figure 2. These populations correlate with those obtained by the numerically accurate AIMS method [119,120,121,122] as well as with those computed by the SQC/MMST method [122]. No decay of the $S_1$ population occurs within the first 5 fs, during which the wavepacket travels from the Franck-Condon to the CI region. Then, the $S_1$ population decreases rapidly, reaching half of its initial value at 9 fs and less than $20\%$ within 10-15 fs. After that, a weak-amplitude population recurrence occurs within 15-25 fs. This flat recurrence is caused by the torsional motions which bring the wavepacket back to the PES crossing each 20-25 fs. After the recurrence, the $S_1$ population remains small, but its quenching slows down considerably.

Panels (a, b, c) of Figure 3 show the SE, ESA, GSB contributions to the integral TA PP spectrum, while panel (d) displays the total TA PP spectrum. The spectra are initiated with the pump pulse with ${\omega }_{pu}=4.39$ eV, which is in resonance with the $S_0-S_1$ transition (see Table A1). Note that all observables of the present section are evaluated at the SA6-CASSSCF level. The SE signals obtained at the SA2-CASSCF level are given in Appendix A. These two levels of the electronic structure theory give very similar results. The long-time TA PP spectra are also shown in the Appendix A.

At $\tau =0$, the SE spectrum is concentrated around ${\omega }_{pu}$, as expected. Within the first 10 fs, it rapidly moves to the red. This redshift is accompanied by the significant loss of intensity. Then $I_{int}^{SE}(\tau ,{\omega }_{pr})$ exhibits a pronounced recurrence around $\tau \approx 25$ fs, but the overall intensity of the spectrum after 10 fs is roughly an order of magnitude lower than the initial intensity. This behavior correlates with the $S_1$ population evolution in Figure 2. However, the $S_1$ population recurrence tells us nothing about the details of the wavepacket motion on $S_1$. The SE signal, on the contrary, demonstrates that a fraction of the wavepacket returns almost to the initial position in the Franck-Condon region. This indicates that the SE signal can be used as an informative observable containing more information than the traditional excited-state population. Note that the low-intensity high-amplitude oscillations in $I_{int}^{SE}(\tau ,{\omega }_{pr})$ last for several hundreds of femtoseconds (see Appendix A), indicating a partially reversible wavepacket motion at the CI. This kind of behavior is highly unusual for CI systems, which normally exhibit irreversible ultrafast population transfer [1,2].

The ESA contribution is shown in Figure 3 (b). At short times, $I_{int}^{ESA}(\tau ,{\omega }_{pr})$ has two components: the upper high-intensity component starting at ${\omega }_{pr}\approx 4.2$ eV and the lower low-intensity component starting at ${\omega }_{pr}\approx 2$ eV. These two branches of the ESA spectrum reveal the $S_1-S_4/S_5$ and the $S_1-S_2$ transitions, respectively (see Appendix A). Initially, both ESA components move to the higher energy domain and reach their maxima around 10 fs. Then the intensity of the upper component decreases rapidly, exhibits a minimum around 20 fs, and continues oscillatory evolution at longer times. The lower component vanishes after 10 fs on the intensity scale of the figure. Both SE and ESA signals reflect the same wavepacket motion on the $S_1$ state. They differ in the "spectator states" (manifold 0 for SE and manifold II for ESA) which are used for the detection of the signals. It is not surprising therefore that the SE and ESA signals behave similarly. More precisely, the spectra in panels (a) and (b) look like mirror images of each other: while the SE spectrum moves to the red, reaches a minimum, exhibits a maximum and keeps in oscillating at much lower intensity, the ESA spectrum moves to the blue, reaches a maximum, exhibits a minimum and keeps in oscillating at much lower intensity. This mirror behavior indicates that the energy gap between the states contributing to the SE (ESA) spectrum initially decreases (increases) with $\tau$.

Figure 3 (c) shows the GSB signal which is composed of the cold (positive) and hot (negative) components. The cold spectrum exhibits undamped oscillations with a $\sim 2$ eV amplitude, $\sim 20$ fs period, and $\sim 4$ eV baseline. The spectrum is dominated by the $S_0-S_1$ transitions in the Franck-Condon region (see Appendix A). The hot spectrum emerges at $\sim 10$ fs (that correlates with the initial decay of the $S_1$ population in Figure 2) and shows high-amplitude ($\sim 7$ eV) oscillations with a period of $\sim 20$ fs around the baseline of $\sim 4$ eV. Such a large oscillation amplitude of the hot spectrum is caused by the excess kinetic energy acquired by the wavepacket in the course of its CI-driven internal conversion. The same baseline of 4 eV is a signature of the dominant role of the $S_0-S_1$ transitions in the formation of both cold and hot signals. Similar oscillation periods of the cold and hot components seem to be a remarkable feature of fulvene, caused by the fact that various nuclear motions in this molecule have similar periods around 20 fs.

The total TA PP signal is shown in Figure 3 (d). It is dominated by the ESA contribution. This is a consequence of the large values of TDMs to higher-lying excited states of manifold II.

The above discussions clearly show the close correspondence between the TA PP signals in Figure 3 and the population dynamics in Figure 2. However, it is even more important to ask whether it is possible to establish correlation between the spectral patterns in the TA PP signals and the key molecular motions in the nonadiabatic dynamics. For this purpose, several critical structures in the nonadiabatic dynamics were optimized at the SA2-CASSCF(6,6)/6-31G(d) level, as shown in Figure 1.

The $S_0$ equilibrium in Figure 1(a) shows an exocyclic bond with a single bond double bond alternation. For the $S_1$ equilibrium (b), however, the alternation is reversed and reveals a biradical character[127]. The $S_0-S_1$ CI (c) is similar to the $S_1$ minimum, which implies that the quick nonadiabatic decay may easily take place after the initial vibrational relaxation to the $S_1$ minimum [127]. Importantly, the C1-C6 bond distance is the leading coordinate affecting the CI-driven internal conversion and TA PP spectra. Figure 1 shows clearly that C1-C6, C2-C4, and C3-C5 bond distances are longer near the CI region than those at the $S_0$ and $S_1$ minima. On the other hand, C1-C2 and C1-C3 bond distances are shorter in the vicinity of the CI region, but longer near the Frank-Condon region. In addition, the C2-C1-C3 angle is smaller at the $S_1$ and $S_0$ equilibrium and larger near the $S_0-S_1$ CI. Hence, extension of the bonds C2-C4 and C3-C5, contraction of the bonds C1-C2 and C1-C3, and decrease of the angle C2-C1-C3 accompany the wavepacket motion from the Frank-Condon region toward the CI region (see Figure 1). These skeleton motions finally comprise the symmetric tuning modes governing the CI-driven internal conversion.

The linear-interpolated reaction pathways between the optimized $S_0$ minimum and $S_0/S_1$ CI and the relevant TDMs obtained by the single point calculations at both SA2-CASSCF(6,6)/6-31G(d) and SA6-CASSCF(6,6)/6-31G(d) levels are presented in Figure 4. Both levels of electronic structure theory give very similar results for potential energy profiles and TDMs. It is thus safe to use any of them for the simulation of spectroscopic signals. This is corroborated by the comparison of Figure 3 and its counterpart in the Appendix A.

Along the interpolated pathways that are driven by the above reactive coordinates and shown in Figure 4 (a) and (c), no barrier exists between the Franck-Condon region and the CI seam. In the Franck-Condon region, a large gradient is observed. Clearly, there is a $S_1$ minimum in the pathway toward the CI, and this minimum lies lower than the CI. As a result, the so-called sloped CI is created [127].

The $S_0-S_1$ TDM decreases along the reaction pathways (Figure 4 (b)). Among electronic transitions from the $S_1$ state to higher-lying states (Figure 4 (d)), $S_1-S_4$ and $S_1-S_5$ TDMs are larger than the $S_0-S_1$ TDM at the Franck-Condon region (Figure 4 (b)). This explains why the ESA signals are stronger than the SE signals. Interestingly, an abrupt rise of $S_1-S_3$ TDM is related to the $S_3/S_4$ surface crossing, while the exchange of the TDM values between $S_1-S_4$ and $S_1-S_5$ transitions is due to the $S_4/S_5$ PES crossing (shown in Figure 4 (d)).

Ab initio simulations permit one to explore how the time evolution of TA PP signals is correlated with the dynamics of specific molecular bond lengths and angles. The results are presented in Figure 5. Panel (a) shows how the C1-C6 bond distance evolves in time on the $S_1$ (a) and $S_0$ (e) states. In the first 10 fs, elongation of the C1-C6 bond is observed. When the C1-C6 distance reaches 1.5-1.7 Å, some trajectories jump back to the ground state, as is clear from Figure 5 (a) and (e). The C1-C6 bond displays pronounced oscillations on the $S_1$ state. Starting from 7-10 fs, a large-amplitude stretching motion of the C1-C6 bond emerges also on the ground state. Evolutions of the bond distances in panels (b), (c), (f) and (g) on the $S_1$ state also exhibit qualitatively similar behaviors. Panels (d) and (h) display evolutions of the bond angle C2-C1-C3 on the excited and ground electronic states, respectively. In panel (d), the C2-C1-C3 angle increases during the first 10 fs. When the angle reaches 110$°$-120$°$, only a few trajectories remain on $S_1$. Afterwards, the remaining trajectories produce high-amplitude oscillations of the C2-C1-C3 angle. There are no trajectories on the $S_0$ state before 10 fs (panel (h)). At $\approx 10$ fs, abundant ground-state trajectories emerge for the C2-C1-C3 angles around 110$°$-120$°$. At longer times, a large amplitude of the C2-C1-C3 bending was observed. Note that the evolution of the C2-C1-C3 angle (panels (d) and (h)) is synchronous with the C1-C6 stretching motion (panels (a) and (e)). Such unusually fast bending motion can be explained by the fact that the C1-C6 stretching motion and the C2-C1-C3 bending motion are strongly coupled with each other since both involve the C1 atom. In contrast, other motions (e.g., the C4-C5 stretching, see Appendix A) are only weakly involved in the dynamic propagation, as they either show small oscillation amplitudes or different oscillation periods.

Interestingly, the nuclear motions on $S_0$ are more erratic than those on $S_1$. This is manifested in a larger spread of trajectories in the lower panels in comparison with that in the upper panels. The reason is that internal conversion produces large-amplitude wavepackets which experience vibrational energy redistribution induced by the mode-mode couplings.

By the combination of the above discussions, the relationship between the population dynamics, the key molecular motions, and the TA PP signals is well established. As the carrier frequency of the pump pulse is resonant with the $S_0-S_1$ transition, the nonadiabatic dynamics of fulvene start from the first excited state. In the Franck-Condon region, the large gradient of the $S_1$ state drives ultrafast nuclear motions, leading to pronounced stretching (C1-C6, C1-C2/C1-C3, C2-C4/C3-C5) and bending (C2-C1-C3) motions. Following these motions, the system first arrives in the $S_1$ minimum and then moves close to the CI on the excited state. This motion also leads to the TDM decrease (Figure 4). As a consequence, the SE signal moves to the low-frequency domain quickly, which is consistent with the time scale of the aforementioned nuclear motions along the internal coordinates. In addition, all trajectories follow similar propagation patterns and no dispersion is observed. Thus the SE signal does not show bifurcation in the early stage of the dynamics. Similar patterns were observed in the ESA signal. The only difference is that the ESA signal moves to the higher energy domain because energy gaps between the $S_1$ state and the higher-lying excited electronic states responsible for the relevant transitions increase with time.

When the CI is reached, many trajectories jump back to the ground state, resulting in significant population transfer. Hence, both the SE and ESA signals suddenly become weaker at $\sim 7-10$ fs, which directly reflects the ultrafast nonadiabatic transitions. After that, a few trajectories stay on the excited state for a rather long time. As the $S_1$ minimum is lower than the CI and the shape of the CI is sloped, these trajectories oscillate near the $S_1$ minimum. Due to the lack of dissipation, these oscillations endure for quite a long time and induce the oscillatory patterns in the weak SE and ESA signals (see Appendix A). Whenever the trajectories access the CI, there is a minor fluctuation of the electronic population due to the nonadiabatic transition, which is consistent with the population dynamics in Figure 2. As some trajectories jump back to the ground state, they create the GSB hot signal. Here, the molecular vibrations are highly excited by the excessive energy emerged after internal conversion. This leads to diversity in the trajectory propagation patterns and results in quite erratic hot GSB signals. Due to the small system size and lack of dissipation, these signals show large-amplitude oscillations, and their periods are consistent with the vibrational periods of the few key coordinates, such as the C1-C6 stretching motion. At the same time, the mode-mode coupling is also responsible for spreading of trajectories, producing a broad hot GSB signal. Overall, the TA PP signals provide clear fingerprints which allow us to monitor the population decay and identify the key tuning modes in the nonadiabatic dynamics.

In addition, evolutions of the asymmetric stretching modes are given in the Appendix A. The modes exhibit oscillations with a $\sim 20$ fs period and different initial phases. As a result, these modes produce characteristic pulsating wavepacket, which is symmetric with respect to positive and negative elongations. There are several demonstrations on how to extract frequencies of the coupling modes from the femtosecond spectra [142,143,144]. In the present case of fulvene, however, both coupling and tuning modes exhibit oscillations with a period of $\sim 20$ fs. Hence, we cannot unambiguously obtain information on the coupling modes from the TA PP spectra.

Unfortunately, experimental femtosecond TA PP signals of fulvene are not available. However, the TA PP spectra can be used as important theoretical observables containing useful information on the CI-driven wavepacket dynamics. This is further illustrated by Figure 6 which shows Fourier transforms of the SE (a, c) and ESA (b, d) spectra with respect to the population time $\tau$,

$$I_{int}^k(\Omega ,{\omega }_{pr})=\int d\tau e^{i\Omega \tau }{I}_{int}^k(\tau ,{\omega }_{pr})$$

(28)

($k=SE,ESA$). These analyses are inspired by the successful use of beating map for disentangling electronic 2D spectra [145,146]. Let us consider the SE spectrum $I_{int}^{SE}(\Omega ,{\omega }_{pr})$ (a) first. Its value around ${\omega }_{pr}\sim 4$ eV gives a broad featureless spectrum along $\Omega$ which is a mere indication of the short lifetime of the system in $S_1$. If we move to ${\omega }_{pr}\sim 2.5$ eV, we see a series of relatively broad peaks which are dominated by the high-intensity peak around $\Omega \approx 1500-1700$ ${\mathrm{cm}}^{-1}$, that yields a 20 fs period which dominates all temporal evolutions of the spectra. The peak widths indicate significance of anharmonic effects produced by nonadiabatic couplings. If we inspect the spectrum around ${\omega }_{pr}\sim 2.5$ eV, we observe splitting of the majority of the broad peaks into pairs of narrower subpeaks. This is a clear signature of the mode-mode coupling. The ESA spectrum $I_{int}^{ESA}(\Omega ,{\omega }_{pr})$ (b) also reveals multiple broad peaks dominated by the one at $\Omega \approx 1500-1700$ ${\mathrm{cm}}^{-1}$. The peaks do not exhibit the doublet splitting which is characteristic of the SE peaks, owing to the larger number of transitions contributing to the ESA spectra. Exclusion of the first 10 fs from the Fourier transform (panels b, d) removes the most intense part of the TA PP, and the corresponding $I_{int}^{SE}(\Omega ,{\omega }_{pr})$ and $I_{int}^{ESA}(\Omega ,{\omega }_{pr})$ exhibit much simpler peak structures featuring a well separated $\Omega \approx 1500-1700$ ${\mathrm{cm}}^{-1}$ peak.

4. Conclusions

By combining the SQC/MMST dynamics and the DW representation, we build an efficient protocol for the on-the-fly simulation of femtosecond TA PP signals. The protocol can readily be extended to electronic 2D and other third-order spectroscopic signals.

In the present work, the protocol has been employed to simulate integral TA PP spectrum of fulvene. On the one hand, the SE and ESA contributions to the TA PP spectrum reveal the characteristic times of the CI-driven $S_1-S_2$ internal conversion which can be extracted from the $S_1$ population dynamics. On the other hand, analyses of the GSB-cold, GSB-hot, SE, and ESA components of the TA PP spectrum give richer and much more comprehensive information on the nuclear motions contributing to the nonadiabatic dynamics. This indicates that ultrafast spectroscopic signals can also be considered as information-rich purely theoretical observables, which are complementary to the standard set of dynamic observables which includes evolutions of populations, bond lengths, and bond angles. Taken together, the observables contain information on the peculiarities of the wavepacket motion on the ground and lower-lying excited electronic states, mode-coupling effects, and the leading active modes responsible for the CI-driven $S_1-S_0$ internal conversion in fulvene. For instance, oscillations of the TA PP signals in the time domain can be attributed to the symmetric stretching and bending modes on the $S_1$ and $S_0$ states, and amplitudes of these oscillations provide information on the relative shifts of the corresponding PESs.

Fulvene is a molecular representative of Tully’s third model [118]. Recently, Worth and coworkers compared electronic populations of three Tully’s molecules predicted by the fewest switches surface hopping method, variational multi-configuration Gaussian (vMCG) method, and (numerically accurate) multi-configuration time-dependent Hartree (MCTDH) method [147]. All calculations employed the linear vibronic coupling (LVC) Hamiltonians developed for the three molecules. As was demonstrated in our work, electronic populations are just the most conventional observables, while much deeper and stringent test would require comparison of nonlinear spectroscopic signals. Obviously, construction of LVC Hamiltonians for higher-lying excited states is a difficult and not fully microscopically justified task, but the available LVC Hamiltonians can be used for the evaluation and comparison of SE and GSB contributions to all third-order spectroscopic signals. The next step may involve testing of on-the-fly and LVC-based signals evaluated by different quantum and quasiclassical methods. Our work can be regarded as a first step in this direction, and the DW methodology can become an efficient tool for the evaluation of all spectroscopic signals.