On-the-fly ab initio semiclassical evaluation of time-resolved electronic spectra

We present a methodology for computing vibrationally and time-resolved pump-probe spectra, which takes into account all vibrational degrees of freedom and is based on the combination of the thawed Gaussian approximation with on-the-fly ab initio evaluation of the electronic structure. The method is applied to the phenyl radical and compared with two more approximate approaches based on the global harmonic approximation - the global harmonic method expands both the ground- and excited-state potential energy surfaces to the second order about the corresponding minima, while the combined global harmonic/on-the-fly method retains the on-the-fly scheme for the excited-state wavepacket propagation. We also compare the spectra by considering their means and widths, and show analytically how these measures are related to the properties of the semiclassical wavepacket. We find that the combined approach is better than the global harmonic one in describing the vibrational structure, while the global harmonic approximation estimates better the overall means and widths of the spectra due to a partial cancellation of errors. Although the full-dimensional on-the-fly ab initio result seems to reflect the dynamics of only one mode, we show, by performing exact quantum calculations, that this simple structure cannot be recovered using a one-dimensional model. Yet, the agreement between the quantum and semiclassical spectra in this simple, but anharmonic model lends additional support for the full-dimensional ab initio thawed Gaussian calculation of the phenyl radical spectra. We conclude that the thawed Gaussian approximation provides a viable alternative to the expensive or unfeasible exact quantum calculations in cases, where low-dimensional models are not sufficiently accurate to represent the full system.


I. INTRODUCTION
Vibrationally resolved electronic spectroscopy provides a powerful tool for studying both electronic and nuclear motions in molecular systems.Among other things, this experimental technique gives insight into the photoinduced molecular dynamics and helps elucidate the shape of molecular potential energy surfaces.Its steady-state version, used for many years, has been interpreted mostly through the time-independent Franck-Condon picture, but an alternative, time-dependent approach, [1][2][3] which examines the evolution of the nuclear wavepacket before evaluating the spectra, is equally valid and, indeed, more fundamental.
With the advent of ultrafast spectroscopy, this time-dependent approach became much more popular because it is a natural way to interpret time-resolved spectra.The best-known experimental setup for measuring time-resolved spectra is the pump-probe technique, consisting of a pump laser pulse, which initiates a dynamical process in the system and is followed by a probe pulse, whose response is measured.5][6][7][8] These spectra are commonly analyzed with fitting procedures to recover the kinetics, 9,10 with Fourier transformation to identify the main frequency components, 11 or with simple low-dimensional models to simulate the vibronic features. 7However, computationally feasible full-dimensional methods for simulating vibrationally and time-resolved spectra of large polyatomic molecules are missing.
The evaluation of vibrationally resolved electronic spectra is often limited to approximate potential energy surfaces, e.g., global harmonic models, 12 or exact quantum calculations taking into account only a few degrees of freedom.To include anharmonicity effects without having to evaluate the global potential energy surfaces, a number of trajectory-based on-thefly ab initio (also called "direct dynamics" or "first-principles dynamics") methods have been developed, ranging from Gaussian basis methods, such as ab initio multiple spawning, 13 coupled coherent states, 14 variational multiconfigurational Gaussians, 15 , multiconfigurational Ehrenfest, 16 and Gaussian dephasing representation, 17 to semiclassical methods based on the Herman-Kluk propagator 18,19 , including versions approximating the prefactor [20][21][22] or the time-averaged semiclassical initial-value representation, 23 and its multiple coherent states 24 or divide-and-conquer extensions. 25Despite their success, these methods are rather expensive when it comes to computing the vibrational structure of electronic spectra, since they require a large number of classical trajectories or basis functions to converge.As an alternative, we have recently implemented an on-the-fly ab initio version 26,27 of the singletrajectory thawed Gaussian approximation 28 as a computationally efficient compromise between rather expensive multiple-trajectory semiclassical methods and cheap, but restrictive and often flawed, global harmonic methods.The thawed Gaussian approximation, based on propagating a single Gaussian wavepacket, is commonly regarded as a starting point for developing more accurate approaches, such as the exact wavepacket propagation of Lee and Heller 29 , generalized Gaussian wavepacket dynamics, 30 off-center-guided wavepackets, 31 Hagedorn wavepacket propagation, 32 time-sliced thawed Gaussian propagation, 33 and coupled wavepackets for nonadiabatic molecular dynamics. 34In the hybrid semiclassical approach of Grossmann, 35,36 the original thawed Gaussian approximation is applied to model the "harmonic" degrees of freedom, while treating the rest by more sophisticated semiclassical methods.In this context, the success of the on-the-fly ab initio thawed Gaussian approximation in reproducing steady-state absorption, emission, and photoelectron spectra of rather large systems 26,27 is somewhat surprising.This motivated us to further extend ab initio thawed Gaussian approximation to include the non-Condon effects, e.g., in symmetryforbidden electronic transitions, where the initial wavepacket is not a Gaussian, but a Gaussian multiplied by a polynomial. 29In all cases studied so far, the methods based on the on-the-fly ab initio thawed Gaussian approximation proved to be more accurate than those based on the global harmonic models.
Here, we adapt the on-the-fly ab initio thawed Gaussian approximation in order to evaluate vibrationally and time-resolved electronic spectra, and compare it with the commonly used adiabatic global harmonic approach, which approximates the two potential energy surfaces to the second order about the corresponding minima.Furthermore, we propose and validate a "combined" global harmonic/on-the-fly method, which approximates the ground potential energy surface with a harmonic potential, but propagates the nuclear wavepacket on the excited state with the on-the-fly ab initio thawed Gaussian approximation.The properties of the wavepackets propagated using the on-the-fly, global harmonic, and combined approaches are analyzed through the differences in the spectra they generate.In particular, we discuss the effect of anharmonicity in both electronic states.
As a proof of principle, we apply the methods to the phenyl radical, a species studied for its rich photochemistry induced by visible and ultraviolet light.8][39] Here, we focus our attention on the first excited state Ã2 B 1 , which does not undergo photodissociation and is responsible for the visible spectrum of the phenyl radical.
The n ← π character of the Ã2 B 1 ← X2 A 1 transition from the ground to the first excited state is in agreement with both experimental and theoretical studies. 40,41The high-resolution cavity ring-down spectroscopic study of the phenyl radical predicted a planar geometry of the excited electronic state, involving an expansion of the carbon ring, and excited-state lifetime of 96 ps. 41The vibronic band corresponding to this transition, first characterized by Porter and Ward 42 and later by Radziszewski, 43 exhibits a rich structure arising from the considerable differences between the potential energy surfaces of the ground and excited electronic states.Ab initio calculations 12,40 found strong coupling between the modes, also known as the Duschinsky rotation or mixing, and a large displacement of multiple normal modes 44,45 -the two most displaced modes correspond to the in-plane carbon ring vibrations (shown in Fig. 1 of the Supporting Information).Due to their significant displacement, the anharmonicity of the potential affects the dynamics and the resulting spectrum, as demonstrated using the on-the-fly ab initio semiclassical dynamics. 44,45In addition, we have recently analyzed the Herzberg-Teller contribution, due to the dependence of the electronic transition dipole moment on nuclear coordinates, to the steady-state Ã2 B 1 ← X2 A 1 electronic spectrum of the phenyl radical: 44,45 by comparing the spectra evaluated using the Condon and Herzberg-Teller approximations, we have shown that the Herzberg-Teller correction to the transition dipole moment does not affect the spectrum significantly, justifying the Condon approximation, which will be also employed below for the time-resolved spectra.
Because the time-resolved electronic spectra of phenyl radical have not yet been measured, our calculations provide a prediction of the experimental signal corresponding to the time-resolved stimulated emission (and ground-state bleach).Moreover, our calculations, although not yet confirmed by an experiment, suggest new ways of analyzing time-resolved spectra.Spectra with low vibrational resolution can be described well by their position and width, which are measured by the mean and standard deviation of the normalized spectral lineshape 46,47 -we show how the mean frequency and width of the spectra are directly related to the properties of the thawed Gaussian wavepacket and can be evaluated using a single classical trajectory on the excited electronic surface.
Clearly, the single-Gaussian ansatz imposes a severe constraint on the shape of the nu-clear wavepacket, making the thawed Gaussian approximation susceptible to failure in the presence of nonadiabatic effects or wavepacket splitting.In such situations, one is forced to use exact quantum methods or semiclassical methods based on multiple guiding trajectories, such as those mentioned earlier in this Section, or multiple thawed Gaussian wavepackets.
In addition, the on-the-fly ab initio version of the thawed Gaussian approximation has only been validated on linear spectra, which depend on the wavepacket autocorrelation function.
Such spectra, unless highly-resolved, can validate the time-dependent wavepacket only at times when it returns to the Franck-Condon region because linear spectra depend weakly on the form of the wavepacket outside this region.In general, the approximation is expected to be qualitatively correct in single-well potentials, where the wavepacket does not split during the short-time propagation needed to describe electronic spectra, and in the absence of nonadiabatic couplings.Indeed, the time-resolved spectroscopic signals have been successfully evaluated using the thawed Gaussian approximation on such model potentials. 48,49To justify our prediction of time-resolved electronic spectra of the phenyl radical, we construct a Morse potential corresponding to the single vibrational mode and evaluate the pump-probe spectra using both the thawed Gaussian approximation and exact quantum dynamics.The results imply that the thawed Gaussian approximation is rather accurate at this level of anharmonicity, but also that more-dimensional calculations are required to reproduce the on-the-fly ab initio spectrum.

A. Vibrationally and time-resolved electronic spectroscopy
Let us consider a system described by the total Hamiltonian where Ĥ is the molecular Hamiltonian and Vint (t) represents the interaction potential with the external field E(t) within the electric-dipole approximation: where ˆ µ stands for the molecular electric dipole operator.Above, we have introduced notation, in which the bold font is used for the matrix representation of electronic operators, the nuclear operators are distinguished by the hat ˆsymbol, and the arrows denote threedimensional vectors.
In a pump-probe experiment, the electric field consists of two components, corresponding to the pump and probe pulses: where is the polarization vector and Here, E pu e (t) and E pr e (t) represent the rapidly changing factors of the pump and probe electric fields (i.e., the factors that oscillates on the electronic time scale), whereas δ n is a slowly varying envelope; the pump pulse is centered at time zero and the probe pulse is centered at a delay time τ .For simplicity, we assume that both envelopes are the same and that the fields are linearly polarized.Therefore, we can first keep the molecular orientation fixed and consider only the projections of the dipole operator onto the polarizations of the fields: keeping in mind that in an isotropic sample, in the end one has to average over all possible orientations of the molecule with respect to the external electric fields.

B. Perturbative expressions for the time-resolved electronic spectra
Differential absorption cross-section, which is typically the time-resolved spectrum measured in the pump-probe experiment, is the difference between the pump-on and pump-off spectra of the probe pulse Numerous approaches can be taken to evaluate the differential absorption spectrum, including both perturbative and non-perturbative methods. 2,50The perturbative approach 2,3,46,51,52 requires the third-order time-dependent perturbation theory; we employ this approach here, adopting the notation and derivation of Ref. 53.While higher-order perturbative expressions are simpler for the density operators, the thawed Gaussian approximation relies on the wavefunction representation.The connection between the density-matrix and wavefunction formalisms is described, e.g., in the monographs by Mukamel 2 and Tannor; 3 therefore, in the main text we only mention the relevant approximations, deferring the supporting equations to the Appendix.
Several approximations simplify the analysis of pump-probe experiments substantially: 1. Nonoverlapping pulses-all interactions with the probe pulse follow all interactions with the pump pulse.
2. Ultrashort pulse (or impulsive) approximation-the pulse envelopes δ n are assumed to be delta functions on the nuclear time scale, but long on the electronic time scale; thus, we assume that the wavepacket is created only in one electronic state.
3. Resonance condition-laser fields are almost in resonance with the electronic transitions.
4. Phase matching-overall momentum transfer between the molecule and the pump pulse is approximately zero, i.e., k pu 1 − k pu 2 ≈ 0 because we are interested in the signal in the direction of the probe pulse.Within these approximations, the third-order time-dependent perturbation theory expression for the pump-probe spectrum is 53 where t is the time after the probe pulse, i.e., t := t − τ and C µ,τ is the third-order dipole time autocorrelation with and the molecular dipole operator in the Heisenberg picture In Eq. ( 11), ρ is the initial (stationary) density and μpu RC is the electric dipole operator which incorporates the resonant condition, i.e., its matrix elements corresponding to the transitions between the electronic states j and k are scaled by the Fourier transform of the pump electric field evaluated at the transition frequency ω jk between those electronic states: Equation (10) for C µ,τ (t ) is expanded using Eqs.( 11)- (12) in Appendix A. To derive the final and more specific expressions for the pump-probe spectra used in this paper, we also assume: 1. Three-state system-only the ground and two excited electronic states are relevant for the description of the whole system, whereas other states are either dark or not in resonance; the three states are approximately equally spaced, so that the electric fields are in resonance with both 0-1 and 1-2 transitions, but not with the direct 0-2 electronic transition.2. Zero-temperature limit-only the ground vibrational state of the ground electronic state is initially populated, i.e., ρ00 = |ψ 0,g ψ 0,g |, and all other density matrix elements are zero.
3. Born-Oppenheimer approximation-there is no population transfer between electronic states during the evolution with the molecular Hamiltonian.
4. Condon approximation-the electronic transition dipole moment matrix elements are independent of nuclear coordinates.
The first approximation is an expanded form of the usual two-state model used for the description of linear absorption and emission spectra; here, a third state is needed to account for the excited-state absorption.(Of course, in practice, there may be many more relevant excited states.)Our molecular Hamiltonian is a diagonal 3×3 matrix in electronic states with nuclear Hamiltonians Ĥ0 , Ĥ1 , and Ĥ2 on the diagonal.As for the zero-temperature limit, the assumption that only the ground electronic state is populated is, of course, satisfied to much higher accuracy than the assumption that only the ground vibrational level is occupied.
Nevertheless, the main contributions to the spectrum are often captured even within the zero-temperature limit.As for the validity of the Born-Oppenheimer approximation, the nonadiabatic effects can introduce a decay of the spectrum intensity with increasing probe delay time τ . 53,54For the short times discussed in this work, such decay can be neglected or subsequently added phenomenologically.However, as pointed out in Ref. 54, this does not hold for longer probe delay times.Finally, we note that the Condon approximation is used only at the last step of the derivation; similar expressions are obtained even beyond the Condon approximation.
As shown in Appendix B, the main contributions to the absorption spectrum of Eq. ( 9) are the excited-state absorption, stimulated emission, and ground-state bleach: The ground-state bleach is, within the given approximations, equal to the scaled continuouswave absorption spectrum.If the ultrashort pulse approximation is relaxed, the groundstate bleach is better described as (impulsive) stimulated resonant Raman scattering 46,55 because the short-time dynamics on the excited state, which can proceed between the two interactions with the pump pulse, creates a wavepacket on the ground electronic state.
However, the dynamics of such wavepacket is limited because it does not have time to move far away from the Franck-Condon region during the short pump pulse.In contrast, the stimulated emission and excited-state absorption signals can drift significantly from their original positions because the excited-state dynamics occurs during the delay time τ .
Consequently, the effect of a finite pump pulse can be safely neglected.
Finally, the spectrum of Eq. ( 14) must be averaged over all possible orientations of the molecule with respect to the electric fields.Within the Condon approximation, the averaging results only in a constant scaling factor, which depends on the angle between the polarizations of the pump and probe pulses, ranging from 1/15 for perpendicular to 1/5 for parallel polarizations (see Appendix C).

C. Wavepacket approach to time-resolved stimulated emission
In the following, we concentrate on the time-resolved stimulated emission spectrum and only in the final result we add the ground-state bleach; the excited-state absorption signal can be computed in a similar way as the stimulated emission.Although the excited-state absorption is not negligible, it is commonly found at different frequencies from the other two contributions and can be analyzed separately, unlike the ground-state bleach and stimulated emission, which often overlap. 56e stimulated emission contribution to the correlation function can be rewritten in terms of a correlation function, 2,46 sometimes referred to as fidelity amplitude, 57,58 by defining a non-stationary wavepacket To evaluate Eq. ( 18), one must first propagate the nuclear wavepacket on the excited-state surface for time τ [according to Eq. ( 19)], and then propagate it simultaneously on both ground-and excited-state surfaces for time t .Because the nonstationary state ψ τ depends on τ , every choice of the probe delay time τ requires a new full t propagation on the two surfaces and every t propagation must be long enough to obtain vibrational resolution.In contrast, one or only several delays have to be considered, since the time τ is not related to the vibrational resolution of the spectra.

D. On-the-fly ab initio thawed Gaussian approximation
The thawed Gaussian approximation is one of the simplest approaches to evolving a wavepacket in a general potential.It considers a Gaussian wavepacket, described by its time-dependent position q t , momentum p t , complex symmetric matrix A t , and a complex number γ t , In a harmonic potential, which can even be time-dependent, the Gaussian form of the wavepacket is exactly preserved-the propagation of the wavepacket only requires evaluating the four time-dependent parameters A t , q t , p t , and γ t . 28,59This forms the basis of the thawed Gaussian approximation, where the potential is given by the local harmonic approximation of the exact potential V (q) about the center q t of the wavepacket at time t.The equations of motion for the time-dependent parameters are where m is the mass matrix and L t the Lagrangian.Since the thawed Gaussian approximation requires only a single classical trajectory, along with the corresponding Hessians, it is suitable for an on-the-fly implementation, where the potential information is obtained as-needed from an ab initio electronic structure code.The details of our on-the-fly implementation and, in particular, the conversion from the Cartesian to normal mode coordinates, can be found in Refs.26 and 27.
The local harmonic approximation holds as long as the wavepacket is compact-the spreading of the wavepacket leads to incorrect description of the potential at the tails.Consequently, the single Gaussian wavepacket ansatz cannot describe wavepacket splitting or tunneling through barriers.Nevertheless, the thawed Gaussian approximation has two extremely important properties for molecular quantum dynamics and, in particular, evaluation of vibrationally resolved electronic spectra.First, it is exact in a global harmonic potential, which is often a decent first approximation to the exact potential energy surface, and performs well for harmonic-like potentials, outperforming the global harmonic approaches due to the partial inclusion of anharmonicity.Second, it is valid for short times before the wavepacket splitting takes effect, which is useful for describing low-resolution vibronic spectra and ultrafast time-resolved experiments.
The computational cost of the Hessian evaluation is typically significantly higher than the cost of the corresponding gradients.To address this issue, we (i) evaluate some of the Hessian the potential is not needed for the classical trajectory propagation, which allows us to first propagate the classical trajectory and then evaluate the Hessians at different points along the trajectory in parallel.In addition, the Hessian of the potential changes more slowly than does the gradient, allowing us to evaluate the Hessian only every several steps, and interpolate in between.These two tricks reduce the computational time significantly, with almost no effect on the accuracy.Note that several alternative approaches, based on Hessian updates 60,61 or Gaussian process regression, 62,63 can be used as well.

III. COMPUTATIONAL DETAILS
For the electronic structure calculations, i.e., the evaluation of the energies, gradients, and Hessians, we used density functional theory for the ground state and time-dependent density functional theory for the excited electronic state, as implemented in the Gaussian09 64 electronic structure package, with B3LYP functional and SNSD basis set; 12 this ab initio method was validated in a previous work. 44e time step was 8 a.u.(≈ 0.194 fs).The time-resolved spectra were computed for 65 equidistant probe delays with a time separation of 12 steps (≈ 2.3 fs), resulting in a total time of about 150 fs.The corresponding correlation functions were computed according to Eq. ( 18); the wavepackets were propagated for 1500 steps (≈ 290 fs) on the groundand excited-state surfaces.The ab initio data needed for the wavepacket propagation were acquired in the following way (see Fig. 1): First, a longer "pump" classical trajectory was propagated on the excited electronic state to cover all probe times, i.e., for time 150 + 290 = 440 fs.Then, 65 probe trajectories were run, each with its own initial conditions taken from the snapshots of the pump trajectory at the probe delay times.Finally, the Hessians were evaluated in parallel for all trajectories every 8 steps of the calculation.The described procedure avoids recomputing the same part of the pump trajectory.Once the ab initio data were collected, the wavepacket was propagated in the mass-scaled normal mode coordinates of the ground electronic state.The initial wavepacket was taken as the ground vibrational state of the electronic ground-state potential energy surface approximated by a harmonic potential, i.e., a Gaussian centered at the ground-state minimum geometry with a width determined by the Hessian of the electronic ground-state surface.
The population decay and the interactions with the environment were neglected during the short-time excited-state dynamics, while the intramolecular nuclear dynamics was fully accounted for with our full-dimensional wavepacket propagation method.We applied a phenomenological inhomogeneous Gaussian broadening (half-width at half-maximum of 100 cm −1 ) of the spectra by multiplying the correlation functions by a Gaussian decay function.
Since no particular choice of the pump field is taken, the absolute scaling of the spectrum is not well defined-all spectra were rescaled as indicated in the corresponding figure captions.

IV. RESULTS AND DISCUSSION
The time-resolved electronic spectrum of phenyl radical is evaluated as a sum of the stimulated emission and ground-state bleach contributions; the excited-state absorption term is not included.The stimulated emission term determines the overall dynamics of the spectra, whereas the ground-state bleach is, within our assumptions, independent of time and proportional to the scaled linear absorption spectrum.We compare three methods for evaluating the spectra: the first employs the on-the-fly ab initio thawed Gaussian approximation, which includes the anharmonicities of both ground-and excited-state surfaces, the secondglobal adiabatic harmonic method-expands the two potential energy surfaces to the second order about the corresponding minima, whereas the third method combines the on-the-fly approach for the excited state with the global adiabatic harmonic approach for the ground state.
To analyze a calculated spectrum in more detail, we evaluate the mean (i.e., the first moment) ω τ := ωσ 0 (ω, τ )dω (26)   and width, measured by the standard deviation (i.e., the second central moment) of the normalized spectral lineshapes 46,47 σ 0 (ω, τ ) = σ(ω, τ )/ω (σ(ω, τ )/ω)dω (28)   as functions of the delay time τ .The two spectral moments are closely related to the properties of the nuclear wavepacket and its dynamics on the excited-state surface.The mean is determined by the expectation value of the difference ∆ V between the excited and ground potential energies: whereas the standard deviation of the spectrum is given by the standard deviation of ∆ V : at time τ .
To simplify the expressions for the mean and width of spectra, let us introduce the shorthand notation for the inverse of the real part of the width matrix, and for the differences of energies, forces, and Hessians of the two potential energy surfaces at the current position of the wavepacket.Then, within the local harmonic approximation for the potential and assuming the thawed Gaussian wavepacket, one can analytically evaluate Eq. ( 29) for the spectral mean, to obtain and Eq. ( 30) for the spectral width, to obtain If the two potentials have similar curvatures, i.e. ∆V τ ≈ 0, the mean of the spectrum corresponds to the energy gap evaluated at the center of the Gaussian wavepacket, i.e., while Eq. ( 36) for the width reduces to showing that the width of the spectrum depends both on the width of the nuclear wavepacket and the gradient of the energy gap evaluated at the center of the wavepacket.Equations (35)   and ( 36) are exact for general harmonic potentials, while Eqs.( 37) and ( 38) are exact for displaced harmonic oscillators with the same force constants.In addition, the evaluation of the four expressions for the mean and width does not require running any probe trajectories; all that is needed is the excited-state thawed Gaussian propagation and the additional information on the ground-state potential evaluated along the excited-state trajectory.
Finally, let us point out that the broadening of spectra introduced as a decay of the correlation function can affect the width of the spectra in a nonlinear way.In particular, for the Gaussian broadening of the form exp(−αt 2 ), equations (36) and (38) for the widths must be modified to where k stands for either "LHA" or "approx.",depending on the level of approximation.
Equation ( 39) is used in the calculations presented below.See Appendices D and E for the derivations of the expressions for ω τ and ∆ω τ given here.
A. On-the-fly ab initio time-resolved electronic spectrum The time-resolved stimulated emission spectrum (Fig. 2) shows oscillations over a broad range of frequencies-a 20000 cm −1 window has to be probed to capture all the features.Although the absorption spectrum covers only about 5000 cm −1 , the mean of the time-resolved stimulated emission spectrum moves according to the energy gap between the two electronic states [Eq.( 29)], thus covering different frequency regions.These oscillations reflect directly the dynamics of the wavepacket on the excited-state surface: As the wavepacket leaves the Franck-Condon region, the energy gap decreases, shifting the spectrum towards lower frequencies.Along with the oscillations in the mean frequency of the spectrum, the vibrational resolution also changes periodically as a function of the delay time.In particular, as the wavepacket moves away from the Franck-Condon region, the spectra become broader and less resolved.The period of ≈ 36 fs observed in the time-resolved spectrum corresponds to the frequency of the most-displaced mode in the excited electronic state (924 cm −1 ).Similar observations, namely the high-amplitude oscillations of the time-resolved stimulated emission spectra and nearly constant ground-state bleach signal, have been found in calculations based on the harmonic models 46,55 as well as in experiments. 4,7e means and widths of the spectra can be reproduced directly from the parameters of the wavepacket and the difference between the ground-and excited-state potentials.Figure 3 shows that the means and widths of the spectra are almost perfectly reproduced within the local harmonic approximation, while slightly larger errors are observed with further approximating ∆Hess q V ≈ 0. Since the spectra themselves are computed within the local harmonic approximation, the means and widths evaluated from wavepacket properties using Eqs.( 35) and ( 39) with k ≡ "LHA" should, indeed, match exactly the means and widths

Spectra
Local harmonic approximation  26)] and standard deviation [Eq.( 27)], the green dots are evaluated from the semiclassical trajectory within the local harmonic approximation [Eq.(35)   and Eq. ( 39) with k ≡ "LHA"], and the green empty squares are evaluated from the semiclassical trajectory using Eq. ( 37) for the means and Eq. ( 39) with k ≡ "approx."for the widths.
calculated from spectra [using Eqs. ( 26) and ( 27)]; however, slight numerical errors result in differences of the order of 80 cm −1 for the means and 40 cm −1 for widths.The errors due to neglecting the difference between the Hessians of the two potential energy surfaces are an order of magnitude greater (≈ 800 cm −1 for the means and ≈ 300 cm −1 for the widths).
Nevertheless, both the means and widths are well described even within this additional approximation.Although the two moments do not represent the full, vibrationally resolved spectra completely, they are very useful for describing low-resolution electronic spectra or spectra, where the vibrational resolution is lost (e.g. in solution).Finally, we show that including the difference between the Hessians of the two potentials modifies the results only slightly; the main determinant of the time dependence of the spectral means is the energy gap evaluated along the classical trajectory, while the spectral widths require only the differences between the gradients of the two surfaces.

B. Comparison with the global harmonic methods
The on-the-fly approach is the most costly, yet also the most accurate, 26,27,44 of the three methods discussed here.The least expensive of the three methods and also the most commonly used approach for evaluating spectra is the adiabatic harmonic approximation, in which the ground-and excited-state potentials are expanded to the second order about the corresponding minima.When accurate, this approximation reduces the computational cost both by avoiding the evaluation of multiple Hessians and by employing closed-form solutions.
Finally, since the pump-probe spectrum is determined by the dynamics on both potential energy surfaces, one can also consider combining the on-the-fly and global harmonic methods.
For example, if excited-state dynamics is not feasible due to the size of the system and/or if the excited-state potential energy surface is well approximated by a harmonic potential, one can use the on-the-fly method only for the ground-state dynamics.In contrast, if the number of the ground-state trajectories is the bottleneck for evaluating the full time-frequency map, the on-the-fly ab initio trajectory may be computed only for the excited state, and the ground-state surface approximated by a harmonic potential.This "combined" approach is the third of the methods whose merits are analyzed in Fig. 4.
The differences observed in Fig. 4 between the global harmonic and more accurate onthe-fly spectra can be assigned to the anharmonicity effects, which play a role at all probe delay times.When the wavepacket is far from the Franck-Condon region, the vibrational resolution is lost, so the difference can be described by the means and widths of the spectra (see Fig. 4, left panels, and Fig. 5).In contrast, when the wavepacket returns to the initial position, the differences are in the relative intensities of the individual vibronic bands in the spectrum (see Fig. 4, right panels).This pattern is repeated over time and follows the oscillatory motion of the wavepacket.The errors are, however, enhanced at later times due to accumulation of the anharmonicity effects during the dynamics on the excited-state surface.
Nevertheless, the global harmonic approximation offers a decent description of the means and widths of the spectra, but at a cost of losing the accuracy of the frequency-resolved features.Interestingly, opposite results are obtained with the combined on-the-fly/global harmonic approach (see Fig. 5), which fails to reproduce the on-the-fly means and widths of the spectra at times when the wavepacket is far from the Franck-Condon region, while describing almost perfectly the vibrationally resolved spectra when the wavepacket is in the vicinity of the initial position.This can be understood easily: When the wavepacket returns to the Franck-Condon region, it is near the ground-state minimum, where the ground-state potential is well described by the global harmonic approximation.At the other turning point, far from the Franck-Condon region, the lack of anharmonicity in the global harmonic approximation of the ground-state potential leads to the underestimation of the energy gap between the two surfaces, and, consequently, shifts the combined on-the-fly/global harmonic spectrum towards lower frequencies.This does not happen when global harmonic models are applied to both surfaces because the errors due to omitting anharmonicity in both potential energy surfaces cancel out approximately.

C. Comparison with the exact solution of a one-dimensional model
Due to the Gaussian ansatz, the thawed Gaussian approximation can break down at the turning points, where the wavepackets tend to become very asymmetrical and oscillatory.
Such effects are difficult to investigate using steady-state spectra, which are determined by the autocorrelation function and, therefore, probe the wavepacket's shape only near the Franck-Condon region.Time-resolved spectroscopy poses another challenge to wavepacket propagation methods: when the probe pulse at time τ sends the nuclear wavepacket from the excited state back to the ground state, the wavepacket often reaches regions of higher anharmonicity than in the case of linear, spontaneous or stimulated emission spectroscopy.To investigate the accuracy of the thawed Gaussian approximation throughout the wavepacket propagation, we constructed a simplified one-dimensional model based on the ab initio data for the ground-and excited-state energies of the phenyl radical.The model potential depends only on the most displaced vibrational mode, neglecting the coupling to and dynamics in all other modes.Clearly, such a simplified model cannot describe all spectral features; however, the goal of these calculations is to provide a similar extent of anharmonicitywhich is why we base it on the ab initio data-while allowing for a simple interpretation of the results.Another motivation for the one-dimensional model is the simple structure of the pump-probe spectra in Fig. 4 based on full-dimensional calculations; at first sight, the single dominant progression suggests that only one degree of freedom may be responsible for the major features of the spectra.Exact quantum calculations for the one-dimensional model whose details are available in the Supporting Information will help shed light on this apparent simplicity.the exact quantum dynamics, on-the-fly thawed Gaussian approximation, and adiabatic global harmonic approach.The spectra are presented at the probe delay times which follow the wavepacket during one period of oscillation between the second and third recurrences, i.e., after a rather long anharmonic dynamics on the excited-state surface.The thawed Gaussian approximation performs very well at this level of anharmonicity and for the given time scales, almost perfectly reproducing the exact spectra at all probe delay times, even at times when the wavepacket is far from the Franck-Condon region.Such accuracy requires that the excited-state wavepacket be well described by a Gaussian-indeed, we find only minor difference between the exact and thawed Gaussian wavepackets (see Fig. 3 of the Supporting Information).In contrast to the thawed Gaussian approximation, the global harmonic model fails to recover the exact result.Although the global harmonic model yields only slightly incorrect peak intensities in the steady-state absorption spectrum (see Fig. 4 of the Supporting Information), these differences are enhanced in the time-resolved spectra in Fig. 6, confirming that pump-probe spectroscopy is a more challenging task for approximate methods.
The spectra of Fig. 6 exhibit a similar pattern as the full-dimensional calculations in Fig. 4, namely the periodic spreading and narrowing of the spectra, as well as the timedependent displacement of the spectra following the energy gap between the two surfaces.Nevertheless, the spectral features differ substantially from those observed in the full-dimensional calculation (Fig. 4).First, the one-dimensional model cannot reproduce the loss of vibrational resolution at the delay times when the wavepacket is far from the Franck-Condon region.Second, the broad baseline of the spectra computed using the fulldimensional approach is missing.Third, the double-headed progression, found in the steadystate absorption spectrum of the phenyl radical 43 and observed at some delay times (e.g. at τ = 72 fs, see Fig. 4) does not appear in the one-dimensional results.
Finally, to probe the limits of the thawed Gaussian approximation, we constructed a more challenging example, where the parameters of the Morse potential were modified in order to enhance the anharmonicity.Whereas the resulting thawed-Gaussian spectra (see Fig. 5 of the Supporting Information) are still qualitatively correct, they do differ from the exact spectra due to the increased anharmonicity.The thawed Gaussian approximation provides good estimates of the center and width of the wavepacket evolving in the excited electronic state (see Fig. 6 of the Supporting Information); however, even the minor differences between the exact and thawed Gaussian wavepackets affect the propagation on the ground-state surface.
As expected, the differences in the propagated wavepackets arise due to the inability of the thawed Gaussian approximation to describe the asymmetric deformations of the wavepacket evolved in an anharmonic potential.The time-resolved spectrum evaluated using the adiabatic harmonic model fails completely, even though its linear absorption spectrum agrees with the exact result at least qualitatively (see Fig. 7 of the Supporting Information).In summary, these results confirm that the thawed Gaussian approximation starts to break down much later than does the global harmonic approach; in particular, the thawed Gaussian approximation can remain valid in realistic anharmonic systems, in which the global harmonic approach already fails even qualitatively.

V. CONCLUSION
To conclude, we presented an efficient on-the-fly ab initio method for computing vibrational structure of time-resolved electronic spectra.The method, based on the thawed Gaussian approximation, does not require expensive construction of potential energy surfaces, and is, therefore, applicable to high-dimensional problems.In addition, the method avoids making any assumptions about the global shape of the potential energy surfaces and uses only local dynamical information.Yet, the inherent and severe approximation, assuming a Gaussian form of the wavepacket, is expected to fail in more anharmonic potentials.
Therefore, we complemented our ab initio calculations of the time-resolved spectra of the phenyl radical by comparing the thawed Gaussian approximation with the exact solution of a one-dimensional model based on the most-displaced vibrational mode of the phenyl radical.The one-dimensional results support the ab initio calculations by confirming that (i) the thawed Gaussian approximation remains valid in a one-dimensional system with comparable anharmonicity and (ii) that the on-the-fly ab initio spectra can only be reproduced by including multiple degrees of freedom and couplings between them.
To evaluate the effect of anharmonicity on the pump-probe spectra, we compared the on-the-fly approach with two approximate methods based on the global harmonic approximation.The results confirm that the anharmonicity of the excited-state potential has different effects on the spectrum, depending on the probe delay time-at delays for which the wavepacket is near the Franck-Condon region, the anharmonicity influences the relative intensities of well-resolved vibronic peaks, whereas at delays for which the wavepacket is located far from the initial geometry, the vibrational resolution is lost and the anharmonicity only modifies the mean and width of the spectrum.Interestingly, the two approximate approaches-one assuming quadratic potential for both surfaces, the other using the global harmonic approximation only for the ground electronic state-break down in different regimes, depending on the position of the nuclear wavepacket at the probe delay time.
The combined on-the-fly/global harmonic method presented here serves as a measure of the anharmonicity effects in both potential energy surfaces; while the on-the-fly and global harmonic approaches give qualitatively similar results, indicating only minor anharmonicity of the potential energy surfaces, the combined approach breaks down completely for time delays when the wavepacket is far from the Franck-Condon point, thus confirming that the In the first step, we have applied the three-state system assumption, zero-temperature limit (constraint on ρ), and Born-Oppenheimer approximation (there is no population transfer during the field-free evolution), while in the second step, we have used the separation Eq. (B1).The first term of (B3) gives the time-resolved stimulated emission-its Fourier transform is centered near ω 10 and modulated by the dynamics on the ground and first excited states.In contrast, the second term of Eq. (B3) vanishes within the rotating-wave approximation-its Fourier transform would lead to a spectrum at negative frequencies around −ω 21 .We recall now that the full expression for the spectrum [Eq.( 9)] involves also the complex conjugate of the correlation function, for which the vanishing term is the complex conjugate of the first term of Eq. (B3), whereas the complex conjugate of the second term of Eq. ( B3) is nonzero and gives rise to the time-resolved excited-state absorption spectrum.
Similarly, the second term in Eq. (A2) is The Fourier transforms of both terms of Eq. (B5) are zero within the rotating-wave approximation.However, this does not hold for their complex conjugates-the complex conjugate of the first term of Eq. (B5) gives rise to the ground-state bleach signal-it is centered near ω 10 , but unlike the time-resolved stimulated emission, does not involve wavepacket propagation on the excited electronic state during the delay time τ .Regarding the complex conjugate of the second term in Eq. (B5), we recall that the pulses are short on the nuclear time scale, which allows the use of the ultrashort pulse approximation, but not on the electronic scale-the averaging of the probe spectra over the short probe pulse duration leads to the cancellation of the terms with a highly oscillatory phase exp(±iω 20 τ ).
For the last term of Eq. (A2) we obtain Again, the first term of Eq. (B7) vanishes due to the rotating-wave approximation, while its complex conjugate does not vanish and corresponds to the ground-state bleach.The second term of Eq. (B7) vanishes due to the highly oscillatory term exp(−iω 20 τ ) and its complex conjugate vanishes both due to the rotating-wave approximation and the oscillatory term exp(iω 20 τ ).
We see that the first term of Eq. (A2) gives the time-resolved stimulated emission and excited-state absorption contributions, while the other two terms contribute to the groundstate bleach through two distinct expressions.If one further assumes that the electronic transition dipole moment elements are independent of the nuclear coordinates, i.e., the Condon approximation, the two contributions to the ground-state bleach are equal and correspond to the scaled linear absorption spectrum.Finally, summing all non-zero terms and reintroducing the nuclear Hamiltonians from Eq. (B1) gives the expressions ( 14)- (17).
Appendix C: Orientational averaging of pump-probe spectra Here we derive the constant scaling factors for the time-resolved spectrum of Eq. ( 14) arising from averaging over all possible orientations of molecules in an isotropic sample.We first consider a rank-4 three-dimensional tensor ← → σ (ω) related to the spectrum by where Einstein's summation convention was used (i.e., a sum over I, J, K, L = 1, 2, 3 is implied).For fixed polarization vectors of the pump ( pu ) and probe ( pr ) pulses in the laboratory frame, one can easily evaluate the spectrum from Eq. (C1).However, the spectrum tensor is a molecular property, and as such, is more naturally given in the molecular frame as ← → σ ijkl (ω), where lowercase indices are used to distinguish the molecular from laboratory frame.Tensors in the two frames are related by where the coefficients λ Ii , λ Jj , ... represent the cosines of angles between the laboratory and molecular axes indicated in the subscript.The orientational average of ← → σ IJKL (ω) is obtained by averaging the coefficients λ Ii λ Jj λ Kk λ Ll .Following the procedure of Craig and Thirunamachandran, 66 but using a slightly more general expression from Ref. 67, one obtains the average Now we fix the laboratory frame so that the pump pulse is polarized along the Z-axis, and let the probe pulse be polarized along an axis, denoted α, tilted by an angle α with respect to the Z-axis.Then the sum (C1) reduces to a single term The molecular frame is set so that the z-axis is parallel to the transition dipole moment vector.Within the Condon approximation, only one component of the spectrum tensor in the molecular frame is nonzero and the Eq.(C2) simplifies to In general, the transition dipole moment is a complicated function of the geometry and can change its orientation during the probe delay time, which gives rise to other nonzero spectrum tensor elements.The reorientation of the transition dipole moment can also occur due to the rotation of the molecule, which is, however, neglected here because of the ultrashort time scales involved.Applying Eq. (C3) gives: Thus, the orientationally averaged time-resolved spectrum is 5 to 15 times weaker than the spectrum that would be obtained if the polarizations of both the pump and probe electric fields were aligned with the transition dipole moment of the molecule-the exact ratio depends on the angle α between the pump and probe polarization vectors.Within the Condon approximation, the orientation average is needed only for absolute cross sections; if only relative cross sections are required, the orientational averaging is unnecessary because the effect of averaging is equivalent to a scaling of the spectra by a constant factor.
Appendix D: Derivation of the expressions for the means and widths of the time-resolved stimulated emission spectra Let a, κ, κ 1 , and κ 2 be D × D real symmetric matrices and let a be positive definite; in our setting, D denotes the number of degrees of freedom.In this Appendix, we will use the following analytical expressions for three Gaussian integrals: 17,68,69 Within the local harmonic approximation, the difference between the effective timedependent potentials of the two surfaces is where we used definitions (32)-(34) of ∆V τ , ∆V τ , and ∆V τ .To evaluate the mean of the normalized time-resolved stimulated emission spectrum (28), we start from Eq. ( 29) and insert the expression for the potential difference (D4): Within the thawed Gaussian approximation the probability density |ψ τ | 2 is even, so the terms containing odd powers of (q − q τ ) vanish in the integral over all space and the expression that its zeroth moment is 1.The final expressions for the first and second moments without broadening are equal to those from Ref. 46 up to constant factors which vanish after dividing by the zeroth moment.
The nth moment of the normalized spectral lineshape σ 0 (ω, τ ) at a probe delay time τ is given by the derivative of the corresponding correlation function evaluated at time zero: We now take our correlation function to be the stimulated emission correlation function (18)   multiplied by a decay function f (t ) responsible for the phenomenological broadening of the spectra.Then the expression for the nth moment of the normalized spectral lineshape reads Before continuing, we note that for the commonly used exponential decay, which results in the Lorentzian broadening of the spectral lines, the standard deviation of the spectrum is not defined.We, therefore, used the Gaussian decay function f (t ) = exp(−αt 2 ).Since f (0) = 1 and f (0) = 0, the first moment is not affected by this decay function.However, the second moment changes to To obtain Eq. (E5) from Eq. (E4), we used the result of Pollard et al. for the first term, and the second derivative of the Gaussian function, for the second term.Then, the data was fitted to a Morse potential

I. MOST DISPLACED MODES OF THE PHENYL RADICAL
The fitting parameters for the two states are given in Table I and the potentials are shown in Fig. 2.
The initial state was obtained as the ground vibrational state of the harmonic fit to the ground-state Morse potential, i.e., it was a Gaussian with the initial position corresponding to the ground-state minimum (ground-state q 0 in Table I), zero initial momentum, and standard deviation σ q = 14.42 (in atomic units).We note that the vibrational modes experience significant coupling and that the chosen coordinate is not a vibrational mode of the ground-state potential; a multi-dimensional model would be required to accurately describe the dynamics of the phenyl radical.Therefore, the initial wavepacket discussed here, which is obtained from a harmonic fit to the ground-state model potential, is not physical, i.e., it does not correspond to any of the degrees of freedom of the full-dimensional initial wavepacket.Nevertheless, the goal of this calculation was to present the accuracy of the thawed Gaussian approximation in moderately anharmonic systems by comparing it to the exact dynamics-we did not aim to recover the spectrum of the phenyl radical.Yet, our model was based on the ab initio data, and so has a realistic extent of anharmonicity.
Exact quantum dynamics was performed using the second-order split-operator method; the same time step and total propagation times were used as for the thawed Gaussian propagation (see Section III. of the main text).The spatial grid consisted of 16384 equidistant points between -200 and 160 atomic units.

IV. TIME-RESOLVED SPECTRA, EXCITED-STATE WAVEPACKET, AND ABSORPTION SPECTRUM OF THE MODEL WITH INCREASED ANHARMONICITY
Here, we compare the thawed Gaussian approximation with the exact result for a more challenging Morse potential, which is constructed from the ab initio model described in the Supporting Information Section II by increasing the anharmonicity.This is accomplished by modifying the a and d parameters so that the product 2da 2 , equal to the force constant of the global harmonic approximation at the minimum of the true potential, remains unchanged.
In this example, we multiply the ground-and excited-state a parameters by 2 and divide the d parameters by 4 (the original parameters are given in Table I).

FIG. 1 .
FIG.1.We first evaluate a long "pump trajectory" moving on the excited-state surface and then use its snapshots as initial conditions for "probe trajectories" moving on the ground-state surface and for evaluating the correlation functions (left panel), instead of recomputing the pump trajectory for each probe delay time (right panel).

1 FIG. 2 .
FIG. 2. Time-resolved stimulated emission spectrum (left) and the pump-probe spectrum including both the stimulated emission and ground-state bleach (right) of phenyl radical evaluated with the on-the-fly ab initio thawed Gaussian approximation.Both spectra were rescaled according to the maximum of the time-resolved spectrum in the right panel.

1 ]FIG. 3 .
FIG.3.Means (top panel) and widths (measured by the standard deviation, bottom panel) of the time-resolved stimulated emission spectra: blue lines are computed from the spectra, using the general definitions of the mean [Eq.(26)] and standard deviation [Eq.(27)], the green dots

FIG. 4 . 1 ]FIG. 5 .
FIG.4.Time-resolved stimulated emission spectra of phenyl radical, evaluated using on-the-fly ab initio thawed Gaussian approximation (TGA), global harmonic, and combined on-the-fly/global harmonic methods, are shown at different delay times τ of the probe pulse.All spectra are scaled by the maximum of the on-the-fly spectrum at zero delay.In addition, as done previously in Ref.44, we introduce frequency shifts by matching the absorption spectra with the experiment: −297 cm −1 for the global harmonic spectrum and −437 cm −1 for the on-the-fly and combined on-the-fly/global harmonic spectra.

FIG. 6 .
FIG.6.Time-resolved spectra of the one-dimensional model of phenyl radical evaluated using the exact quantum dynamics, on-the-fly thawed Gaussian approximation (TGA), and global harmonic approach.Probe delay times are taken as to follow the wavepacket motion during one period of ≈ 36 fs between the second and third recurrences.

Figure 6
Figure 6 compares the time-resolved spectra of the one-dimensional model evaluated with

FIG. 1 .
FIG. 1.Most displaced modes of the phenyl radical corresponding to the totally symmetric inplane vibrations; we show the eigenvectors of the excited-state Hessian evaluated at the optimized geometry of the excited electronic state (ab initio method described in the Section III of the main text).Left: Mode of frequency 924 cm −1 .Right: Mode of frequency 590 cm −1 .The full table of the ground-and excited-state frequencies, along with the relative displacements of the corresponding modes can be found in the Supporting Information of Ref. 1.

FIG. 2 .
FIG. 2. Ab initio data (blue points) and the fitted Morse potentials (blue lines) of the ground and excited electronic states of the phenyl radical along the excited-state vibrational mode of frequency 924 cm −1 ; zero position is the minimum of the excited-state potential energy surface (obtained by the ab initio geometry optimization).

ω
FIG.5.Time-resolved spectra of the modified Morse potentials with increased anharmonicity, defined in Section IV of the Supplementary Information.Analogue of Fig.6of the main text.

FIG. 6 .FIG. 7 .
FIG. 6. Probability density of the excited-state wavepacket at two different probe delay times τ .Analogue of Fig. 3 of Supporting Information, but for the modified Morse potentials with increased anharmonicity.

TABLE I .
Parameters (in atomic units) for the Morse potentials [Eq.(1)] of the ground and excited electronic states.