On-the-fly ab initio semiclassical evaluation of third-order response functions for two-dimensional electronic spectroscopy

Ab initio computation of two-dimensional electronic spectra is an expanding field, whose goal is improving upon simple, few-dimensional models often employed to explain experiments. Here, we propose an accurate and computationally affordable approach, based on the single-trajectory semiclassical thawed Gaussian approximation, to evaluate two-dimensional electronic spectra. Importantly, the method is exact for arbitrary harmonic potentials with mode displacement, changes in the mode frequencies, and inter-mode coupling (Duschinsky effect), but can also account partially for the anharmonicity of the involved potential energy surfaces. We test its accuracy on a set of model Morse potentials and use it to study anharmonicity and Duschinsky effects on the linear and two-dimensional electronic spectra of phenol. We find that in this molecule, the anharmonicity effects are weak, whereas the Duschinsky rotation and the changes in the mode frequencies must be included in accurate simulations. In contrast, the widely used displaced harmonic oscillator model captures only the basic physics of the problem but fails to reproduce the correct vibronic lineshape.


I. INTRODUCTION
Electronic spectroscopy allows us to study excited electronic states and light-induced nuclear dynamics.To track this ultrafast dynamics on femtosecond time scales, a range of time-resolved and two-dimensional spectroscopic techniques were developed.The complex signals obtained in these experiments are, however, difficult to interpret without the help of theoretical modeling. 1st models for describing two-dimensional electronic spectra treat the electronic states as the system and include nuclear dynamics only approximately, as bath effects.The simplest approach assumes that nuclear degrees of freedom induce a Gaussian-like or Lorentzian-like broadening, neglecting completely the coherent nuclear dynamics.Such models are appropriate only if the coherent dynamics is strongly suppressed by the surrounding solvent dynamics.
6][17] Such methods are accurate when the curvatures of different potential energy surfaces are similar and in the limit of strong dephasing, i.e., for short times.Even in this case, however, if the on-the-fly dynamics is performed with ab initio electronic structure methods, evolving the full ensemble of classical trajectories can quickly become prohibitively expensive.To account for both intramolecular and intermolecular nuclear dynamics in two-dimensional spectroscopy, one can use the multimode Brownian oscillator model, 11,[18][19][20] which considers a few primary harmonic modes coupled to a large number of low-frequency bath oscillators.For this model, the spectra can be computed analytically; moreover, the parameters of the model can be computed efficiently from a single ab initio classical trajectory, as demonstrated in Refs.21   and 22.This allows one to perform electronic structure computations at a high level of theory, using, for example, post-Hartree-Fock multiconfigurational wavefunction methods.
Yet, the approach is limited to modeling the molecule as a set of few uncoupled displaced harmonic oscillators.Such a simple description is inadequate in systems that exhibit strong tronic states, or anharmonicity effects.A generalization to a set of uncoupled harmonic oscillators with both displacement and a possible change in the force constant was proposed by Fidler and Engel, who used the approximate third-order cumulant expansion. 23Very recently, the third-order cumulant expansion was also studied as an approximate way to account for the mode-mode coupling (Duschinsky effect) in both linear 24 and nonlinear 25 spectra.
There has been little development in the ab initio simulation of two-dimensional electronic spectra beyond the standard semiclassical methods or displaced harmonic models.To account for anharmonicity effects 26,27 or more general coupled oscillators, one is forced to employ computationally expensive exact quantum dynamics methods, [28][29][30][31][32][33] such as different flavors of the multiconfigurational time-dependent Hartree (MCTDH) method 34,35 , or the hierarchical equations of motion. 36,37These methods require the pre-computation of the full potential energy surfaces and are not suitable for a first-principles on-the-fly implementation.
Here, we propose an efficient semiclassical method to evaluate vibrationally resolved twodimensional electronic spectra.The approach, based on Heller's single-trajectory thawed Gaussian approximation, 57 accounts for inter-mode coupling, changes in the force constants, and, at least partially, for the anharmonicities of the ground-and excited-state potential energy surfaces.First, we study how the accuracy of the method depends on the degree of anharmonicity in the one-dimensional Morse system.The results are compared with the exact benchmark and with the harmonic approximation, which neglects the anharmonicity completely.Second, we analyze the effects of Duschinsky coupling and anharmonicity on the linear absorption and two-dimensional spectra of phenol.

A. Third-order response function
The central object in all types of third-order electronic spectroscopy is the third-order polarization 11,58 (1)   where E(t) is the electric field of light (without the polarization vector) and is the third-order response function, expressed in terms of correlation functions and In Eq. ( 7), Ĥi are the vibrational Hamiltonians of the ground ("1") and excited ("2") electronic states, ρ = exp(−β Ĥ1 )/Tr[exp(−β Ĥ1 )] is the vibrational density operator in the ground electronic state, and μ = μ21 = ˆ µ 21 • is the electronic transition dipole moment projected on the polarization unit vector of the external electric field.Equations ( 3)- (6)   rely on the following assumptions: (i) due to the large gap between the electronic states, only the ground electronic state is initially populated; (ii) Born-Oppenheimer approximation, i.e., there is no population transfer under field-free evolution; (iii) light pulses are linearly polarized in the same direction ; (iv) only two electronic states are involved.In the following, we discuss how to evaluate the components R α of the response function (2).

B. Zero-temperature limit: Wavepacket picture
In the zero-temperature limit, we assume that only the ground ("g") vibrational state |1, g of the ground electronic state is populated initially, i.e., ρ = |1, g 1, g|.Then, we may rewrite Eq. ( 7) in terms of nuclear wavepackets: = 1, g|μe i Ĥ2 τa/ μe i Ĥ1 τ b / μe −i Ĥ2 τc/ μ|1, g e −iω 1,g (τa+τ b −τc) = φ τ b ,τa |φ 0,τc , where ω 1,g = 1, g| Ĥ1 |1, g , Ĥ i = Ĥi − ω 1,g , and The result (11) has an appealing interpretation in terms of bra and ket wavepackets, which we represent pictorially for R 3 [Eq.( 5)] in Fig. 1.The bra wavepacket is first evolved in the excited electronic state for a time τ a = t 1 and then for a time τ b = t 2 +t 3 in the ground state, where it is a non-stationary wavepacket due to the initial t 1 dynamics on the excited-state potential energy surface.The ket wavepacket "waits" during the t 1 and t 2 times and is only evolved for a time τ c = t 3 in the excited-state.This simple picture has been discussed in the literature in the context of pump-probe 59 and two-dimensional 28 spectroscopy.In general, during the t 1 (coherence) and t 3 (detection) times, the bra and ket wavepackets evolve on different potential energy surfaces, i.e., the system is in a state of electronic coherence; during t 2 , also called population or waiting time, both nuclear wavepackets are in the same electronic state, i.e., the system is in an electronic population state. 19e evaluation of R α functions requires only one excited-state wavepacket evolution up to time t 1 + t 2 + t 3 and, in addition, wavepackets propagated in the ground electronic state starting from the snapshots along the excited-state trajectory.Since such calculations would be difficult to perform with multiple-trajectory direct dynamics methods, we employ the efficient, single-trajectory thawed Gaussian approximation.
where q t and p t are D-dimensional position and momentum vectors, A t is a complex and symmetric D × D matrix with positive-definite imaginary part, and γ t is a complex number whose imaginary part ensures the normalization.D is the number of degrees of freedom.
The wavepacket (12) solves exactly the time-dependent Schrödinger equation where is the local harmonic approximation of the true potential V (q) about q t , if the Gaussian's parameters satisfy the system 57 In Eq. ( 18), L t = T (p t ) − V (q t ) is the Lagrangian of the classical trajectory (q t , p t ).The above equations are interpreted as follows: the phase-space center (q t , p t ) of the Gaussian (12) evolves classically with the exact classical Hamiltonian, the complex matrix A t evolves according to the Hessian computed at the current q t , and the complex number γ t is updated according to the Lagrangian of the classical trajectory (q t , p t ) and the matrix A t .Since the only source of error is the local harmonic approximation ( 14), the thawed Gaussian propagation is exact for arbitrary, multi-dimensional harmonic potentials.
The method was originally proposed for problems involving short-time dynamics, such as photodissociation spectra. 60However, its accuracy appears to be surprisingly satisfactory in molecular systems even at longer times, because many molecules are only weakly to moderately anharmonic. 61,62Using a single thawed Gaussian wavepacket, which is the essence of Heller's thawed Gaussian approximation, is rather restrictive but also very efficient for on-the-fly dynamics coupled to ab initio electronic structure.The on-the-fly ab initio thawed Gaussian approximation 63 proved useful in treating efficiently anharmonicity effects on linear absorption, 61,62,64,65 emission, 65,66 and photoelectron spectra, 64 as well as in understanding nuclei-induced electronic decoherence in attosecond experiments. 67Recently, we extended our on-the-fly implementation of the single-Gaussian approach to simulate frequency-and time-resolved pump-probe spectra, 68 similar to the earlier work by Rohrdanz and Cina 69 and Braun et al. 70 on model potentials.Although ensembles of thawed Gaussians were largely discarded and replaced by frozen Gaussians due to the numerical instabilities that often appear in nonadiabatic dynamics simulations, thawed Gaussians are being reintroduced, e.g., in the semiclassical hybrid dynamics [71][72][73] or Gaussian-based MCTDH, [74][75][76] and especially for spectroscopic applications, [77][78][79][80] due to their ability to describe couplings between different degrees of freedom. 81

D. Two-dimensional electronic spectroscopy
A variety of different third-order experiments can be simulated through the computation of the response function (2). 11For example, the full response function is needed for evaluating transient absorption spectra with finite-duration pulses. 82Here, we focus on the two-dimensional electronic spectra obtained from the individual correlation functions R α with t 2 = 0 (i.e., at zero delay time).
Spectra at nonzero delay could be obtained by using t 2 > 0 in Eq. (19).Spectra S α represent the ideal signals obtained in the limit of ultrashort pulses.In a more general setting with finite pulses, the two-dimensional spectra are computed from the time-dependent polarization (1), which involves explicitly the electric fields. 19,83To ensure that all spectra appear at positive frequencies ω 1 , nonrephasing spectra (α = 1, 4) are computed with the positive sign in the exponent of Eq. ( 19), while the negative sign is used for the rephasing spectra (α = 2, 3). 19Furthermore, it is easy to see from Eqs. ( 3)-( 6) that R 1 (t 3 , 0, and R 2 (t 3 , 0, t 1 ) = R 3 (t 3 , 0, t 1 ); hence, we will show only two sets of spectra: S 1 ≡ S 4 and S 2 ≡ S 3 .In general (i.e., for arbitrary t 2 ), correlation functions R 1 and R 2 are associated with the stimulated emission process because the system evolves in the excited state during the population time t 2 ; functions R 3 and R 4 correspond to the ground-state bleaching because the system is in the ground electronic state during the t 2 delay time.For t 2 = 0, one cannot distinguish between these two processes.
To analyze the accuracy of different approximate approaches, we introduce the spectral contrast angle between the reference [S (ref) ] and approximate (S) spectra, where is the inner product of two two-dimensional spectra and S = √ S • S is the associated norm.

III. COMPUTATIONAL DETAILS
A. One-dimensional models: Harmonic and Morse potentials An arbitrary one-dimensional harmonic potential, is described by the equilibrium position q eq , energy minimum V eq , and frequency ω.We set the mass m = 1 in all of our model calculations.Let us also define a one-dimensional Morse potential, in terms of the anharmonicity parameter χ and the parameters V eq , q eq , and ω, which relate to the harmonic potential ( 22) fit to the Morse potential at q eq .We construct a set of one-dimensional systems composed of the ground-state harmonic potential, and the excited-state Morse potentials, V 2 (q) = V Morse (q; V 2,eq = 0, q 2 = 1.5, ω 2 = 0.9, χ), of variable anharmonicity χ ranging from 0.006 to 0.02.The initial vibrational state, i.e., the ground vibrational state of the ground electronic state, is a Gaussian due to the ground-state harmonic potential.The exact two-dimensional electronic spectra are compared to approximate spectra evaluated either with the harmonic approximation or with the thawed Gaussian approximation.Within the harmonic approximation, the excited-state Morse potential is replaced by the harmonic potential Note that the harmonic result does not depend on the anharmonicity parameter χ of the Morse potential.
Next, we compare the harmonic and thawed Gaussian approximations for a one-dimensional system composed of two Morse potentials V 2 (q) = V Morse (q; V 2,eq , q 2 , ω 2 , χ = 0.01), (28)   with the same degree of anharmonicity.The exact initial state is no more a Gaussian.
However, in the thawed Gaussian simulations, we approximate it by the vibrational ground state of the harmonic potential (24) fitted to the ground-state Morse potential ( 27) at its minimum.The harmonic approximation replaces both ground-state and excited-state potential energy surfaces by the harmonic potentials; the result is the same as for the harmonic-Morse system described above.
Wavepacket propagation was performed for 150 steps in both t 1 and t 3 times and with a time step of 0.2.The transition dipole moment was set to 1 (Condon approximation).
Correlation functions R 1 , R 2 , R 3 , and R 4 were multiplied by a Gaussian damping function exp[−a(t 2 1 + t 2 3 )] with a = 0.014427, resulting in the Gaussian broadening (half-width at half-maximum of 0.2) of the spectra along both frequency axes.The exact spectra were computed in the eigenstate representation, which is feasible for these one-dimensional systems since both harmonic and Morse eigenfunctions are known; 84 the associated Franck-Condon overlaps were computed numerically.

B. On-the-fly ab initio calculations
The electronic structure of phenol was modeled using the density functional theory with the PBE0 functional and 6-311G(d, p) basis set, as implemented in the Gaussian 16 quantum chemistry package. 85Excited-state calculations were performed with the time-dependent density functional theory.This choice of electronic structure theory provides ground-state frequencies similar to those computed at the MP2/aug-cc-pVDZ level (see Table II of the supplementary material and Ref. 86) and transition energies along the excited-state trajectory that agree, up to an approximately constant shift (which results only in a shift of the computed spectrum but does not affect its shape), to those evaluated at the EOM-CCSD/6-311G(d, p) level (Fig. 1 of the supplementary material).A single ab initio excited-state trajectory was run for 1000 steps starting from the ground-state optimized geometry; subsequent ground-state classical trajectories were propagated for 500 steps.Overall, the calculations allow the evaluation of the correlation functions with 500 steps in both t 1 and t 3 delay times; t 2 delay was set to zero.All dynamics simulations used a time step of 0.25 fs and a standard second-order Verlet integrator.The ab initio calculations evaluated not only the energies and gradients at each step but also the Hessians of the electronic energy.These potential energy data were transformed to ground-state normal mode coordinates and used to propagate the 33-dimensional wavepacket according to Eqs. ( 15)- (18).After evolving the ground-and excited-state Gaussian wavepackets, the correlation functions were computed using Eqs.( 3)-( 6), (11), and the expression for the overlap of two thawed Gaussian wavepackets with parameters q i , p i , A i , and γ i (i = 1, 2).In Eq. ( 29), we defined vectors and scalars as well as the notation δΛ := Λ 2 − Λ * 1 for Λ = A, ξ, η.To construct the harmonic model, also known as the generalized Brownian oscillator model, 24 of phenol, an additional Hessian was computed at the optimized excited-state geometry.This corresponds to the so-called adiabatic Hessian or adiabatic harmonic model. 63,87o more approximate models were also studied: The uncoupled harmonic model was obtained by neglecting the off-diagonal terms of the excited-state Hessian expressed in terms of the ground-state normal modes.The displaced harmonic oscillator model, also called the Brownian oscillator model, was constructed by replacing the excited-state Hessian in the adiabatic harmonic model by the ground-state Hessian; this specific way of constructing the displaced harmonic oscillator parameters is called the adiabatic shift approach. 63,87When applied to any of these different harmonic potentials, the thawed Gaussian propagation is exact and enables an efficient evaluation of linear and two-dimensional spectra.Although explicit expressions are available for the evaluation of linear absorption and emission spectra of harmonic systems, 88 no such analytical approaches have been presented for the two-dimensional spectra of arbitrarily shifted, distorted, and rotated harmonic potentials.Spectra simulations assumed Condon approximation for the transition dipole moment.
Linear absorption spectra were computed from the first 500 steps of the excited-state wavepacket autocorrelation function (see Fig. 2 of the supplementary material, where the convergence is confirmed, and Ref. 63 for more details) and were broadened by a Gaussian with half-width at half-maximum of 120 cm −1 ; same broadening was used for the twodimensional spectra along both ω 1 and ω 3 frequency axes.This corresponds to a phenomenological inhomogeneous broadening; homogeneous broadening due to direct system-bath interactions is neglected.The system-bath coupling would be needed for spectra at later delay times t 2 > 0, as the system would have time to relax and dissipate energy to the environment; we assume that the response functions with t 2 = 0 and t 1 , t 3 < 125 fs are only weakly affected by the system-bath coupling.

IV. RESULTS AND DISCUSSION
A. Model potentials

Harmonic-Morse system
Two-dimensional spectra for the harmonic ground-state potential and Morse excitedstate potential are shown in Fig. 2. The exact nonrephasing spectrum appears only along the diagonal, whereas the rephasing spectrum exhibits a characteristic checkerboard pattern due to vibronic transitions that involve various ground-and excited-state vibrational states.proximation reproduce the exact spectra well, which is not the case for the harmonic results.
The nonrephasing harmonic spectrum flattens out at higher frequencies (spectral region A indicated in the top right panel of Fig. 2), unlike the exact and thawed Gaussian spectra, which exhibit clear vibronic peaks at these frequencies.Similar effects are seen in the rephasing spectra, mostly in the spectral region labeled B (see Fig. 2, bottom right).Again, the exact spectrum is composed of a long vibronic progression up to ω 3 = 6, which, in the harmonic approximation, is truncated around ω 3 = 4.In the region C, the harmonic spectrum is missing negative vibronic peaks, which are reproduced well by the thawed Gaussian approximation.The thawed Gaussian approximation, however, suffers from another form of error: as in linear spectroscopy (see, e.g., Ref. 64), artificial negative peaks may also appear in the two-dimensional spectra, which is most obvious in the nonrephasing spectrum of Fig. 2 around (ω 1 , ω 3 ) ≈ (−1, −1).
We now compare the exact and approximate spectra at different levels of anharmonicity by measuring the error (see Fig. 3) through the spectral contrast angle [Eq.(20)].The thawed Gaussian approximation exhibits smaller errors in the computed spectra than the harmonic approximation at all levels of anharmonicity and for both rephasing and nonrephasing spectra.As expected, the accuracy of both approximate approaches deteriorates as the anharmonicity of the system increases.

Morse-Morse system
Two-dimensional rephasing spectra of the system composed of two Morse potentials, both with the anharmonicity parameter χ = 0.01, are shown in Fig. 4. As in the harmonic-Morse system, the errors of the harmonic spectrum are observed in the spectral regions B and C; the accuracy of the thawed Gaussian spectrum is not much affected by the additional anharmonicity in the ground-state potential surface.To analyze further the two approximate methods, we inspect one-dimensional cuts of the two-dimensional spectra along two different values of ω 1 frequency (Fig. 5).We see clearly that the thawed Gaussian approximation recovers the positions and intensities of the vibronic peaks both at low ω 1 ≈ 1 and high ω 1 ≈ 4 frequencies.0][91] If the simulation, for example, based on a model harmonic potential, cannot reproduce the vibronic peaks found in the experimental spectra, these peaks might end up incorrectly assigned to another electronic state or another excitation process.

B. Two-dimensional electronic spectrum of phenol
Phenol is an ultraviolet chromophore present in proteins as the residue of the naturally occurring amino acid tyrosine.Recently, accurate electronic structure methods were employed to simulate its two-dimensional electronic spectrum 94,95 in an attempt to explore theoretically the capabilities of this spectroscopic technique to resolve the features specific to chromophore-chromophore interactions in oligopeptides and, more generally, in proteins. 3,96ese recent calculations included multiple electronic states but neglected the vibronic structure of the individual electronic transitions.Here, we present a complementary result: we focus only on the ground and first excited electronic states, i.e., we neglect the excited-state absorption process, but study in detail the vibronic lineshape of the ground-state bleaching/stimulated emission signal.The methods we use neglect the nonadiabatic effects; this is ) and the spectra computed with the on-the-fly ab initio thawed Gaussian approximation, harmonic approximation, uncoupled harmonic model ("Uncoupled"), and displaced harmonic oscillator (DHO) model.For ease of comparison, the computed spectra are shifted in frequency and scaled in intensity so that they all match at the maximum of the experimental spectrum (see Sec. IV of the supplementary material for details).
an acceptable approximation for the dynamics in the first excited state of phenol, as demonstrated by the MCTDH simulations performed on a vibronic-coupling Hamiltonian model of phenol. 86e linear absorption spectrum of phenol was computed with four different approximate methods: the on-the-fly ab initio thawed Gaussian approximation, harmonic approximation, uncoupled harmonic model, and displaced harmonic oscillator model (Fig. 6).Harmonic and on-the-fly thawed Gaussian spectra (Fig. 6, top) are similar in accuracy for this specific system: while the thawed Gaussian propagation results in more accurate intensities of the low-frequency peaks, namely, the 0-0 transition at ≈ 36350 cm −1 and the shoulder at ≈ 36800 cm −1 , the harmonic approximation gives a better estimate of the high-frequency region and the tail of the spectrum.One of the main disadvantages of the thawed Gaussian approximation, the appearance of artificial negative spectral intensities, shows up clearly in the absorption spectrum of phenol.Although the simulated harmonic and thawed Gaussian spectra resemble the experiment, there are remaining differences, most notably in the intensities of the spectral peaks.These errors could be either due to anharmonicity effects not captured by the approximate thawed Gaussian wavepacket propagation or due to the errors in the potential energy data evaluated with an approximate electronic structure method.
Fairly small difference between the harmonic and thawed Gaussian spectra suggests that the anharmonicity effects are weak and that the remaining errors in our simulation are due to the inaccuracies of the electronic structure theory.In the bottom panel of Fig. 6, we show the results of two more approximate approaches-the uncoupled harmonic and displaced harmonic oscillator models.These spectra clearly deviate from the experiment, indicating the importance of both mode distortion (changes in mode frequencies) and intermode couplings (Duschinsky effect).When going from the displaced harmonic model, which neglects mode distortion, to the uncoupled harmonic model, which includes mode distortion, the peaks broaden but still exhibit inaccurate intensities.Additional inclusion of the Duschinsky effect, which is achieved by moving to the (coupled) harmonic model, improves the intensities.
The two-dimensional spectra simulated with different approximate methods are shown in Fig. 7. Again, the uncoupled harmonic and displaced harmonic oscillator models predict spectra that differ substantially from the harmonic and thawed Gaussian results, which are, in turn, similar to each other.Therefore, based on both linear and two-dimensional spectra simulations, we may conclude that the anharmonicity effects are truly weak in the ground and first excited states of phenol, at least in the region explored by the nuclear wavepacket for short time after the photoexcitation.More precisely, the anharmonicity effects are much weaker than the effects of the Duschinsky rotation and frequency changes, which are, in contrast, significant, as demonstrated by the simulations based on the uncoupled or displaced harmonic models.The anharmonicity could, however, play a role at longer simulation times, needed, for example, to simulate high-resolution spectra.We note that most simulations supporting experimental results are nowadays performed with the simplified displaced harmonic oscillator model, which captures the basic physics of the problem, but is inadequate in certain cases, such as the presented example of phenol.
Interestingly, the spectrum spans a broad range of frequencies in both ω 1 and ω 3 , which

V. CONCLUSIONS AND OUTLOOK
We have presented a new method for simulating vibrationally resolved two-dimensional electronic spectra that is exact for any shifted, distorted, and coupled harmonic model and, in addition, can approximately account for anharmonicity effects.The method, based on the thawed Gaussian approximation, is shown to be superior to the harmonic approximation for a series of Morse models of varying anharmonicity.On the example of phenol, we show that inter-mode couplings and changes in the mode frequencies, both of which are frequently neglected in simulations, can be crucial for recovering the correct vibronic shape of the twodimensional electronic spectra.In this specific case, the anharmonicity is shown to be weak, which could allow further studies on the nonlinear spectra of phenol based on the harmonic approximation.For example, our results could be augmented by constructing harmonic models with more accurate electronic structure methods, in order to simulate excited-state absorption signals.For systems that do exhibit anharmonicity effects, we propose the onthe-fly ab initio thawed Gaussian approximation as a computationally affordable approach beyond harmonic approximation.
Finally, let us also give a short outlook on how to include features that are missing in the current method.First, as a wavefunction method, the thawed Gaussian approximation is not suitable for treating systems at non-zero temperature.We have shown recently that this limitation can be overcome efficiently with the so-called thermo-field dynamics theory. 97rrently, we are exploring the application of this idea to the computation of nonlinear spectra.Second, the method is originally constructed for isolated systems.An obvious, "ab initio way" to augment the system with an environment would be to include a number of solvent molecules directly into the system.To account for inhomogeneous broadening, the dynamics would have to be repeated for different conformations of the solute-solvent system.
Alternatively, the bath effects could be treated through a number of low-frequency harmonic oscillators coupled to the system; the procedures for computing the parameters of the bath oscillators are well-studied in the literature.The extensions that include temperature and environment effects would enable accurate and efficient first-principles simulation of timeresolved (t 2 > 0) two-dimensional electronic spectra in the condensed phase.FIG. 1.The energy gaps between the excited and ground electronic states of phenol evaluated along the first 500 steps of the excited-state ab initio trajectory at the PBE0/6-311G(d, p) level of theory, compared to the transition energies computed using the EOM-CCSD/6-311G(d, p) method.
The nearly constant shift in the transition energies induces only a constant frequency shift in the computed spectra, but does not affect their shape.TABLE IV.Overall frequency shifts (in cm −1 ) introduced into the calculated linear and twodimensional spectra of phenol.Same shifts were used along both frequency axes in two-dimensional electronic spectra.

FIG. 2 .
FIG. 2. Exact, thawed Gaussian, and harmonic two-dimensional electronic spectra for the harmonic-Morse system described in Sec.III A with the anharmonicity of the excited-state Morse potential χ = 0.01.Spectral regions A, B, and C, discussed in the main text, are indicated on the harmonic spectra.

FIG. 7 . 3 of
FIG.7.Rephasing two-dimensional electronic spectra of phenol at t 2 = 0 computed with the onthe-fly ab initio thawed Gaussian approximation, harmonic approximation, uncoupled harmonic model, and displaced harmonic oscillator (DHO) model.The spectra correspond to the stimulated emission (S 2 ) and ground-state bleaching (S 3 ) processes (S 2 = S 3 for t 2 = 0).Computed spectra were shifted along both frequency axes and scaled in intensity as in Fig.6.

FIG. 2 .
FIG. 2. Top: Undamped (red, dashed) and damped (blue, solid) wavepacket autocorrelation function 1, g|e −i Ĥ 2 t/ |1, g , evaluated with the on-the-fly ab initio thawed Gaussian approximation, and the damping (black, dotted), which is a Gaussian function exp(−t 2 /2σ 2 t ) with σ t = √ 2 ln 2/HWHM, where HWHM is the half-width of half-maximum of the broadening Gaussian function in the frequency domain (HWHM= 120 cm −1 , as reported in the main text).Bottom: Absorption spectra computed as the Fourier transforms of the first 500 (blue, solid) or 1000 (red, dashed) steps of the damped autocorrelation function (blue, solid line in the top panel).For the given damping, 500 steps (125 fs) of the autocorrelation function are sufficient to obtain a converged absorption spectrum, as demonstrated by the comparison of simulated spectra in the bottom panel.

TABLE I .
Optimized geometries (in Å) of the ground and first excited electronic states of phenol

TABLE II .
Ground-and excited-state frequencies and dimensionless relative displacements between the two states.Frequencies are given in cm −1 .Dimensionless displacements are defined as | ω i / q 2,i |, where ω i is the i-th ground-state frequency and q 2,i is the excited-state position in the mass-scaled normal mode coordinate i of the ground state.Both ground-and excited-state frequencies are listed in the descending order.