Finite-temperature, anharmonicity, and Duschinsky effects on the two-dimensional electronic spectra from ab initio thermo-field Gaussian wavepacket dynamics
FFinite-temperature, Anharmonicity, andDuschinsky Effects on the Two-dimensionalElectronic Spectra from Ab Initio Thermo-fieldGaussian Wavepacket Dynamics
Tomislav Beguˇsi´c ∗ and Jiˇr´ı Van´ıˇcek ∗ Laboratory of Theoretical Physical Chemistry, Institut des Sciences et Ing´enierie Chimiques,Ecole Polytechnique F´ed´erale de Lausanne (EPFL), CH-1015, Lausanne, Switzerland
E-mail: tomislav.begusic@epfl.ch; jiri.vanicek@epfl.ch a r X i v : . [ phy s i c s . c h e m - ph ] F e b bstract Accurate description of finite-temperature vibrational dynamics is indispensablein the computation of two-dimensional electronic spectra. Such simulations are oftenbased on the density matrix evolution, statistical averaging of initial vibrational states,or approximate classical or semiclassical limits. While many practical approaches exist,they are often of limited accuracy and difficult to interpret. Here, we use the concept ofthermo-field dynamics to derive an exact finite-temperature expression that lends itselfto an intuitive wavepacket-based interpretation. Furthermore, an efficient method forcomputing finite-temperature two-dimensional spectra is obtained by combining theexact thermo-field dynamics approach with the thawed Gaussian approximation forthe wavepacket dynamics, which is exact for any displaced, distorted, and Duschinsky-rotated harmonic potential but also accounts partially for anharmonicity effects ingeneral potentials. Using this new method, we directly relate a symmetry breaking ofthe two-dimensional signal to the deviation from the conventional Brownian oscillatorpicture.
Graphical TOC Entry
To this end, a numberof theoretical methods were developed to account for typical vibrational-electronic ef-fects occurring in molecular systems, such as anharmonicity, different curvatures of theground- and excited-state potential energy surfaces, or mode-mode mixing (Duschinsky rota-tion).
In its original formulation, the second-order cumulant expansion is exact onlyfor the Brownian oscillator (i.e., displaced harmonic) model and cannot treat the inter-modecoupling in the excited state. Although this basic molecular model shaped our understandingof steady-state, ultrafast, and multidimensional electronic spectroscopy in the past decades,it is inadequate for many molecules that exhibit Duschinsky and anharmonicity effects.
Similar limitations are met when using the semiclassical phase averaging, also knownas the Wigner-averaged classical limit or dephasing representation.
The recentlydeveloped third-order cumulant approach seems to overcome these limitations.
Yet, itis only accurate in systems with weakly coupled or distorted modes. In contrast, quantum dynamics methods are well suited for describing the evolu-tion of nuclear wavepackets but often neglect temperature effects. To avoid the impracticalBoltzmann averaging over the initial states, a number of alternative strategies for includingtemperature in wavepacket-based methods have been proposed.
We turn to the so-calledthermo-field dynamics, which transforms the von Neumann evolution of a density ma-trix to a Schr¨odinger equation with a doubled number of degrees of freedom. This approachhas only recently been introduced in chemistry for solving the electronic structure, vi-bronic, and spectroscopic problems at finite temperature. Here, we show how it couldbe used to compute two-dimensional vibronic spectra. The finite-temperature treatment iscombined with the thawed Gaussian approximation, an efficient first-principles methodfor wavepacket propagation, and applied to the stimulated emission and ground-state bleachsignals of azulene. 3n two-dimensional spectroscopy, a nonlinear time-dependent polarization P (3) ( t, t , t ) = (cid:18) i ~ (cid:19) Z ∞ dt Z ∞ dt Z ∞ dt R (3) ( t , t , t ) E t ,t ( t − t ) E t ,t ( t − t − t ) E t ,t ( t − t − t − t )(1)is induced in the sample through interaction with the electric field E t ,t ( t ) comprised ofthree light pulses centered at times − t − t , − t , and 0, where t is the delay between thefirst two pulses, t is the delay between the second and third pulses, and R (3) ( t , t , t ) is thethird-order response function. In a heterodyne detection scheme, the measured signal is S ( t , t , t ) ∝ i Z ∞−∞ E LO ( t ) P (3) ( t, t , t ) dt, (2)where E LO ( t ) is the fourth, local oscillator pulse centered at time t after the third pulse. Thetwo-dimensional spectrum is obtained by scanning S ( t , t , t ) as a function of the three timedelays and Fourier transforming over t and t . We focus on the absorptive two-dimensionalspectrum ˜ S ( ω , t , ω ) = ˜ S R ( ω , t , ω ) + ˜ S NR ( ω , t , ω ) , (3)and assume ultrashort and nonoverlapping pulse approximations, where the rephasing andnon-rephasing spectra, ˜ S R ( ω , t , ω ) = Re Z ∞ dt e iω t Z ∞ dt e − iω t S R ( t , t , t ) , (4)˜ S NR ( ω , t , ω ) = Re Z ∞ dt e iω t Z ∞ dt e iω t S NR ( t , t , t ) , (5)are defined through S R ( t , t , t ) = C ( t + t , t , t + t ) + C ( t , t + t , t ) − C ( t , t , t + t + t ) ∗ , (6) S NR ( t , t , t ) = C ( t , t , t + t + t ) + C ( − t , − t , t ) − C ( t + t , t , t + t ) ∗ , (7)4nd C i ( τ a , τ b , τ c ) = Tr[ ˆ ρ ˆ µ e i ˆ H τ a / ~ ˆ µ i e i ˆ H i τ b / ~ ˆ µ i e − i ˆ H τ c / ~ ˆ µ e − i ˆ H ( τ a + τ b − τ c ) / ~ ] . (8)In Eq. (8), ˆ H i are the vibrational Hamiltonians corresponding to the ground ( i = 1)and excited ( i = 2 ,
3) electronic states, ˆ ρ = exp( − β ˆ H ) / Tr[exp( − β ˆ H )] is the vibrationaldensity operator at temperature T = 1 /k B β , and ˆ µ ij = ˆ ~µ ij · ~(cid:15) is the electronic transitiondipole moment between electronic states i and j projected on the polarization unit vector ~(cid:15) of the external electric field. The correlation function C corresponds to the stimulated emis-sion and ground-state bleach processes, while C , which involves a higher excited electronicstate, corresponds to the excited-state absorption (see Sec. 1 of the Supporting Information).Although the excited-state absorption term involves, in general, a sum over several higherexcited states ( i > = 3), here, for brevity, we consider only one such state.An intuitive physical interpretation of Eq. (8) is available in the zero-temperature limit,where the density operator ˆ ρ = | , g ih , g | is given in terms of the ground vibrational state | , g i of the ground electronic state. Then, C i ( τ a , τ b , τ c ) = h φ ( i ) τ b ,τ a | φ ( i )0 ,τ c i , (9)where | φ ( i ) τ,t i = e − i ˆ H i τ/ ~ ˆ µ i e − i ˆ H t/ ~ ˆ µ | , g i , (10)ˆ H i = ˆ H i − ~ ω ,g , and ~ ω ,g = h , g | ˆ H | , g i . In Fig. 1, we illustrate how Eq. (9) is evaluatedfor the stimulated emission contribution C ( t + t , t , t + t ) [Eq. (8)] to the rephasing signal[Eq. (6)]. The bra nuclear wavepacket is first evolved for a time τ a = t + t in the excitedelectronic state and then for a time τ b = t in the ground state; the ket wavepacket is in theground electronic state during t and evolves in the excited state for a time τ c = t + t . Ingeneral, during the time delays t and t , also known as coherence and detection times, thebra and ket wavepackets evolve on different potential energy surfaces; during the so-calledpopulation time t , the two wavepackets are in the same electronic state—in the ground state5or the ground-state bleach contribution and in the excited electronic state for the stimulatedemission and excited-state absorption components. V V t + t t 〈ϕ t , t + t |( a ) t + t V V | ϕ t + t 〉 ( b ) V V | ϕ t + t 〉〈ϕ t , t + t |( c ) Figure 1: Evolution of the bra (a, dotted line) and ket (b, solid line) wavepackets of Eq. (9)for τ a = t + t , τ b = t , and τ c = t + t . Their overlap (c) is the stimulated emission term C ( t + t , t , t + t ) [Eq. (8)] of the rephasing signal S R ( t , t , t ) [Eq. (6)].We now address the question “Is it possible to retain the simple wavepacket picturewithout neglecting finite-temperature effects?” To answer this question in the affirmative,we employ thermo-field dynamics, which maps the evolution of a density operator at finitetemperature to the evolution of a wavefunction with a doubled number of coordinates. Inthe thermo-field dynamics theory, the thermal vacuum is defined as | ¯0( β ) i = ˆ ρ / X k | k ˜ k i , (11)where | k ˜ k i = | k i| ˜ k i is the basis vector of the tensor-product space obtained from the physical(with basis {| k i} ) and “fictitious” (with basis {| ˜ k i} ) Hilbert spaces. We note that physicaloperators (denoted only by hat ˆ, such as ˆ ρ or ˆ µ ) act only on the physical subspace. Withthese definitions, Eq. (8) can be rewritten as C i ( τ a , τ b , τ c ) = h ¯ φ ( i ) τ b ,τ a | ¯ φ ( i )0 ,τ c i , (12)where | ¯ φ ( i ) τ,t i = e − i ˆ¯ H i τ/ ~ ˆ µ i e − i ˆ¯ H t/ ~ ˆ µ | ¯0( β ) i (13)6s the analogue of | φ ( i ) τ,t i from Eq. (10),ˆ¯ H i = ˆ H i − ˆ˜ H ( i = 1 , ,
3) (14)is the Hamiltonian acting in the full, tensor-product space, and ˆ˜ H is the ground-statevibrational Hamiltonian acting in the fictitious space only. The proof of Eq. (12) goes asfollows: C i ( τ a , τ b , τ c ) = h ¯0( β ) | ˆ µ e i ˆ¯ H τ a / ~ ˆ µ i e i ˆ¯ H i τ b / ~ ˆ µ i e − i ˆ¯ H τ c / ~ ˆ µ | ¯0( β ) i (15)= X k ,k h k ˜ k | ˆ ρ / ˆ µ e i ( ˆ H − ˆ˜ H ) τ a / ~ ˆ µ i e i ( ˆ H i − ˆ˜ H ) τ b / ~ ˆ µ i e − i ( ˆ H − ˆ˜ H ) τ c / ~ ˆ µ ˆ ρ / | k ˜ k i (16)= X k ,k h k | ˆ ρ / ˆ µ e i ˆ H τ a / ~ ˆ µ i e i ˆ H i τ b / ~ ˆ µ i e − i ˆ H τ c / ~ ˆ µ ˆ ρ / | k ih ˜ k | e − i ˆ˜ H ( τ a + τ b − τ c ) / ~ | ˜ k i (17)= X k ,k h k | ˆ ρ / ˆ µ e i ˆ H τ a / ~ ˆ µ i e i ˆ H i τ b / ~ ˆ µ i e − i ˆ H τ c / ~ ˆ µ ˆ ρ / | k ih k | e − i ˆ H ( τ a + τ b − τ c ) / ~ | k i (18)= X k h k | ˆ ρ / ˆ µ e i ˆ H τ a / ~ ˆ µ i e i ˆ H i τ b / ~ ˆ µ i e − i ˆ H τ c / ~ ˆ µ e − i ˆ H ( τ a + τ b − τ c ) / ~ ˆ ρ / | k i (19)= Tr[ ˆ ρ ˆ µ e i ˆ H τ a / ~ ˆ µ i e i ˆ H i τ b / ~ ˆ µ i e − i ˆ H τ c / ~ ˆ µ e − i ˆ H ( τ a + τ b − τ c ) / ~ ] . (20)Equation (15) is obtained from Eq. (12) by inserting the definition (13) of ¯ φ ( i ) τ,t , while Eq. (16)results by substituting relation (11) for | ¯0( β ) i ; in going from Eq. (16) to (17) we used thefact that operators acting in different subspaces commute; in going from (17) to (18) we usedthe conjugation rules relating the physical and fictitious spaces (see Sec. 2 of the SupportingInformation). Resolution of identity and commutation of ˆ ρ / with ˆ H were used to obtainEq. (19), and the definition and cyclic property of the trace to get Eq. (20).Remarkably, the result (12) has exactly the same form as the zero-temperature expres-sion (9) and can be interpreted as in Fig. 1. It also allows finite-temperature effects to beincluded in regular wavefunction-based codes, by modifying only the definition of the initialstate and the Hamiltonians under which this state is evolved. In Sec. 3 of the Support-7ng Information, we prove that the same wavepacket picture can be justified even beyondthe Born–Oppenheimer approximation, which was invoked implicitly in Eqs. (6)–(8). Toavoid exponentially scaling exact quantum methods on precomputed potential energy sur-faces or computationally demanding multiple-trajectory approaches, we proposeusing the simple, yet efficient, single-trajectory thawed Gaussian approximation, which canbe interfaced with on-the-fly ab initio evaluation of potential energy information. Let us consider a Gaussian wavepacket ψ t ( q ) = e i ~ [( q − q t ) T · A t · ( q − q t )+ p Tt · ( q − q t )+ γ t ] , (21)where q t and p t are the real, D -dimensional expectation values of position and momentum, A t is a D × D complex symmetric matrix with positive-definite imaginary part, γ t is acomplex scalar whose imaginary part ensures normalization of the wavepacket, and D is thenumber of coordinates. Within the thawed Gaussian approximation, one replaces the truepotential energy V ( q ) by its local harmonic approximation V LHA ( q ) = V ( q t ) + V ( q t ) T · ( q − q t ) + 12 ( q − q t ) T · V ( q t ) · ( q − q t ) (22)about the center q t of the wavepacket, which leads to the following equations of motion forthe Gaussian’s parameters: ˙ q t = m − · p t , (23)˙ p t = − V ( q t ) , (24)˙ A t = − A t · m − · A t − V ( q t ) , (25)˙ γ t = L t + i ~ m − · A t ) , (26)where L t = p Tt · (2 m ) − · p t − V ( q t ) is the Lagrangian along the trajectory ( q t , p t ) and m is the symmetric mass matrix. According to Eqs. (23)–(26), the position and momentum8f the Gaussian wavepacket evolve classically, while the matrix A t depends on the Hessiansalong the classical trajectory. The described evolution of the Gaussian wavepacket is exactfor a harmonic potential because the local Taylor expansion of Eq. (22) becomes exact inthis case. For more general, anharmonic potentials, the method is only approximate, buttypically accurate for moderate anharmonicity and short times, which makes it practical inspectroscopic applications. Although the thawed Gaussian propagation is not suitedfor nonadiabatic dynamics, it can treat accurately the effects that arise due to different forceconstants of the ground- and excited-state potential surfaces: mode distortion, i.e., thechange in the frequency of a normal mode, and inter-mode coupling or Duschinsky rotation.The on-the-fly ab initio thawed Gaussian approximation, which uses electronic structurecalculations to compute potential energies, gradients, and Hessians only when needed, wasrecently validated for the simulation of finite-temperature linear and zero-temperaturetwo-dimensional spectra. To construct the initial state, we approximate the ground-state potential energy surfaceby a harmonic potential and use the corresponding mass-scaled normal mode coordinates.Then, in the zero-temperature limit, the initial state ψ ( q ) = h q | , g i is a Gaussian (21)and D = F , where F is the number of vibrational degrees of freedom. In the thermo-fielddynamics formulation, D = 2 F , the initial state ¯ ψ (¯ q ) = h ¯ q | ¯0( β ) i is also a Gaussian, and¯ q = ( q, ˜ q ) is the 2 F -dimensional coordinate vector. To solve the equations of motion in thefinite-temperature picture, we need the potential energies, gradients, and Hessians in the ex-tended coordinate space, which can be easily formulated in terms of the energies, gradients,and Hessians of the two potential energy surfaces, as shown in Ref. 56. Remarkably, thethermo-field dynamics under Hamiltonian ˆ¯ H i [Eq. (14)] requires exactly the same classicaltrajectory—in the electronic state i —as the conventional, zero-temperature thawed Gaus-sian propagation with the Hamiltonian ˆ H i . No further ab initio evaluations are needed forthe finite-temperature implementation, meaning that, within the thawed Gaussian approxi-mation, the temperature effects can be included almost for free. The only difference in the9omputational cost is in solving the equations of motion with 2 F rather than F coordinates,which is approximately 2 = 8 times more expensive due to the roughly cubic scaling of theinvolved matrix operations, including matrix-matrix multiplication and matrix inverse. Thiscost is, however, negligible compared to the cost of electronic structure calculations.Formally, the propagation of the wavepacket according to Eqs. (23)–(26) requires not onlythe potential energies and gradients but also the Hessians at each step of the dynamics. Inthis work, we employed the single-Hessian method, which further approximates V ( q t ) ≈ V ( q ref ) in Eq. (25), where q ref is a reference geometry at which the Hessian of the excited-state potential surface is evaluated once and reused during the excited-state dynamics. Sincethe center of the wavepacket still follows the fully anharmonic classical trajectory, the single-Hessian version partially includes anharmonicity effects; in several examples studied in Ref.86, this method was shown to be of similar accuracy as the thawed Gaussian approximation.Here, we chose q ref as the excited-state minimum. The ground-state potential surface wasassumed to be harmonic in all simulations.To analyze the effects of the excited-state anharmonicity, we compare the anharmoniccalculations, based on the on-the-fly single-Hessian thawed Gaussian approximation for theexcited-state propagation, with the harmonic model (also called the generalized Brownianoscillator model), where the excited-state potential surface is approximated by a harmonicpotential fitted to the surface at its minimum (so-called adiabatic harmonic or adiabatic Hes-sian scheme). In the mass-scaled normal mode coordinates of the ground state, the excited-state force constant is a symmetric, non-diagonal matrix, whose off-diagonal terms reflectinter-mode couplings, also known as Duschinsky mixing. To study the effects of the differ-ence between the excited- and ground-state force constants on linear and two-dimensionalspectra, we construct the displaced harmonic model (also called the Brownian oscillatormodel), where the excited-state force constant is approximated by the force constant in theground electronic state. This model neglects mode distortion and Duschinsky effects. Thetwo-dimensional spectra can be computed exactly with the thawed Gaussian propagation,10s described above, for both harmonic and displaced harmonic oscillator models. Whereasthe exact solution to the displaced harmonic oscillator model was known before in the formof the second-order cumulant expansion, to the best of our knowledge, no method hasbeen published for computing exactly the two-dimensional spectra of the global harmonic(or generalized Brownian oscillator) model. Azulene is a well-known example of a Kasha-violating molecule, as it emits light fromthe second, rather than first, excited electronic state. This is due to the interplay of twofactors: (i) weak nonadiabatic coupling between S and S states and (ii) fast ( ≈ to S . These properties make azulene one of the key buildingblocks in the synthesis of novel optoelectronic materials. Although nonadiabatic couplingsbetween the ground and first excited states play an important role in the photoinduceddynamics of azulene, they do not affect its vibrationally resolved S ← S absorptionspectrum. Indeed, the linear absorption spectrum can be well reproduced using adiabatic,Born-Oppenheimer approaches that neglect nonadiabatic effects. Here, we also ignorethe nonadiabatic effects on the two-dimensional spectra, which we compute only at short t delay times. In the results, we focus on the ground-state bleach and stimulated emissioncontributions to the two-dimensional spectrum [the first two terms on the right-hand sidesof Eqs. (6) and (7)].In Fig. 2 (top), we compare linear absorption spectra simulated at 300 K and 0 K withthe experimental spectrum recorded at room temperature. One of the main effects of tem-perature is the broadening of the spectral features, which also affects the relative intensitiesof vibronic peaks, namely, those at 14 300 cm − and 15 800 cm − . These intensities are over-estimated in the zero-temperature spectrum, but are corrected by the finite-temperaturetreatment.Non-zero temperature has an even stronger effect on the two-dimensional spectrum(Fig. 2, bottom). The zero-temperature spectrum is composed of sharp vibronic peaks,which are broadened and less resolved in the spectrum computed at 300 K. As in the linear11 xp T = T = [ cm - ] I n t en s i t y M a x . i n t en s i t y Figure 2: Top: S ← S absorption spectra of azulene computed with the on-the-fly ab initiosingle-Hessian thawed Gaussian approximation at zero temperature (red, dashed) and at300 K (blue, dotted), compared with the experimental spectrum (black, solid) recorded atroom temperature in cyclohexane. Bottom: Two-dimensional electronic spectra [Eq. (3)]at zero delay time ( t = 0), computed at zero temperature (left) and at T = 300 K (right).Each spectrum shows the sum of the ground-state bleach and stimulated emission terms [firsttwo terms on the right-hand sides of Eqs. (6) and (7)] corresponding to the S –S electronictransition in azulene. See Fig. S1 for the rephasing and non-rephasing contributions to thesespectra and Figs. S3 and S4 of the Supporting Information for the spectra at delays t > Exp Anharmonic Harmonic DHO14 000 16 000 18 000 20 000 22 000 24 0000.00.20.40.60.81.0 Wavenumber [ cm - ] I n t en s i t y M a x . i n t en s i t y Figure 3: S ← S absorption spectra of azulene computed with the on-the-fly ab initiosingle-Hessian thawed Gaussian approximation (“Anharmonic”, blue, dotted), harmonic ap-proximation (red, dashed), and the displaced harmonic oscillator (DHO) model (green, dash-dotted) at 300 K, compared with the experimental spectrum (black, solid) recorded at roomtemperature in cyclohexane. The corresponding two-dimensional spectra (Fig. 4, top) exhibit similar differences, whichwe can conveniently analyze in the time domain (see Fig. 4, bottom, for | S R ( t , , t ) | andFig. S6 for | S NR ( t , , t ) | ). The displaced harmonic oscillator model results in stronger13ecurrences after 45 fs (in t , t , or in both t and t ) than the harmonic or anharmonicapproaches. This translates into sharper peaks in the two-dimensional spectrum. The an-harmonic spectrum extends less into the high frequency region, compared to the harmonicand displaced harmonic oscillator models, because the thawed Gaussian propagation givesa slower initial decay (for t , t < − .Figure 4: Top: Two-dimensional electronic spectra [Eq. (3)] at zero delay time ( t = 0),computed with the on-the-fly ab initio single-Hessian thawed Gaussian approximation (“An-harmonic”, left), harmonic approximation (middle), and the displaced harmonic oscillator(DHO) model (right) at 300 K. Each spectrum shows the sum of the ground-state bleachand stimulated emission terms [first two terms on the right-hand sides of Eqs. (6) and (7)]corresponding to the S –S electronic transition in azulene. See Fig. S2 for the rephasingand non-rephasing contributions to these spectra and Figs. S3 and S4 of the SupportingInformation for the spectra at delays t >
0. Bottom: First 60 fs of | S R ( t , , t ) | [Eq. (6)].See Fig. S6 for | S NR ( t , , t ) | .Interestingly, for the displaced harmonic oscillator model, | S R ( t , , t ) | is symmetric with14espect to the diagonal (Fig. 4, bottom right), which does not hold when mode distortion,rotation, and anharmonicity are included (Fig. 4, bottom left and middle). We prove thisanalytically in the Supporting Information (Secs. 7 and 8), where we also demonstrate thatthe asymmetry can appear only in the rephasing signal | S R ( t , , t ) | . Moreover, we showthat the (incorrect) symmetry of | S DHOR ( t , , t ) | with respect to the diagonal t = t is,more generally, imposed by the second-order cumulant approximation, which is exact forthe displaced harmonic oscillator model and is employed regularly to model two-dimensionalspectra. Hence, the second-order cumulant method cannot account for the asymme-try induced by the deviation from the displaced harmonic oscillator model. This erroneousqualitative behavior was difficult to study in the past, partly due to the absence of practicalmethods that could easily go beyond the second-order cumulants or Brownian oscillators.To conclude, we derived a general and exact expression for computing finite-temperaturevibrationally resolved two-dimensional electronic spectra with wavefunction-based methods.The inclusion of temperature is the key to simulating spectra of larger systems or solvatedmolecules, due to the multitude of low-frequency modes that are thermally excited at roomtemperature. By combining the exact expression with the thawed Gaussian approximation,we developed a practical and efficient method for computing two-dimensional spectra be-yond zero temperature and beyond displaced harmonic oscillator model. With the help ofthe newly developed method, we identified an asymmetry in the time-domain signal thatcould serve as evidence for the changes in mode frequencies, mode-mode coupling, or an-harmonicity. This asymmetry cannot be described with the conventional and widely usedsecond-order cumulant approach.
Computational Methods
The ground electronic state of azulene was modeled at the second-order Møller–Plesset(MP2) perturbation theory level; the first excited state was modeled using the second-15rder Laplace-transformed density-fitted local algebraic diagrammatic construction [LT-DF-LADC(2)] scheme, as implemented in the Molpro 2015 package. cc-pVDZ basis setwas used throughout (see Ref. 88). We first evaluated the Hessians in the ground and ex-cited states at the respective optimized geometries. Then, starting from the minimum of theground state, an on-the-fly ab initio classical trajectory was evolved in the excited electronicstate for 1130 steps with a time step of 8 a.u. ≈ .
19 fs (total time ≈
219 fs).Linear spectra were computed by Fourier transforming the first 500 steps of the wavepacketautocorrelation function (see Ref. 82). Regarding the simulation of two-dimensional spec-tra, t and t times were propagated up to ≈
106 fs (500 steps); t delays ranged from zero(results shown in the main text) to 25 fs (130 steps), in intervals of 5 fs or 26 steps. Condonapproximation, which was justified for the S ← S absorption of azulene in Ref. 88, wasemployed. Gaussian broadening with a half-width at half-maximum of 90 cm − was used inboth linear and two-dimensional spectra. Linear spectra were shifted in frequency and scaledin intensity to match at the maximum intensity peak of the experiment; two-dimensionalspectra were shifted by the same frequency shifts as the linear absorption spectra and scaledaccording to the maximum of the fully absorptive two-dimensional spectrum [Eq. (3)].Data supporting this publication can be found at http://doi.org/10.5281/zenodo.4552858. Acknowledgement
The authors acknowledge the financial support from the European Research Council (ERC)under the European Union’s Horizon 2020 research and innovation programme (grant agree-ment No. 683069 – MOLEQULE) and from the Swiss National Science Foundation throughthe NCCR MUST (Molecular Ultrafast Science and Technology) Network.
Supporting Information Available
Derivation of Eqs. (6)–(8), conjugation rules in thermo-field dynamics, thermo-field dynam-16cs expression beyond the Born–Oppenheimer approximation, rephasing and non-rephasingcontributions to the spectra of Figs. 2 (bottom) and 4 (top), two-dimensional spectra at t > | S R ( t , , t ) | at short times, (a)symmetry of the rephasing and non-rephasing spec-tra in the time domain, and proof of the symmetry of the two-dimensional signal within thesecond-order cumulant approximation. References (1) Yuen-Zhou, J.; Aspuru-Guzik, A. Quantum process tomography of excitonic dimersfrom two-dimensional electronic spectroscopy. I. General theory and application to ho-modimers.
J. Chem. Phys. , , 134505.(2) Kreisbeck, C.; Kramer, T.; Aspuru-Guzik, A. Disentangling electronic and vibroniccoherences in two-dimensional echo spectra. J. Phys. Chem. B , , 9380–9385.(3) Yuen-Zhou, J.; Krich, J. J.; Kassal, I.; Johnson, A. S.; Aspuru-Guzik, A. UltrafastSpectroscopy ; IOP Publishing, 2014.(4) Zhou, N.; Chen, L.; Huang, Z.; Sun, K.; Tanimura, Y.; Zhao, Y. Fast, Accurate Simu-lation of Polaron Dynamics and Multidimensional Spectroscopy by Multiple DavydovTrial States.
J. Phys. Chem. A , , 1562–1576.(5) Ikeda, T.; Tanimura, Y. Probing photoisomerization processes by means of multi-dimensional electronic spectroscopy: The multi-state quantum hierarchical Fokker-Planck equation approach. J. Chem. Phys. , , 014102.(6) Xiang, B.; Ribeiro, R. F.; Dunkelberger, A. D.; Wang, J.; Li, Y.; Simpkins, B. S.;Owrutsky, J. C.; Yuen-Zhou, J.; Xiong, W. Two-dimensional infrared spectroscopy ofvibrational polaritons. Proc. Nat. Acad. Sci. USA , , 4845–4850.(7) Conti, I.; Cerullo, G.; Nenov, A.; Garavelli, M. Ultrafast Spectroscopy of Photoactive17olecular Systems from First Principles : Where We Stand Today and Where We AreGoing. J. Am. Chem. Soc. , , 16117.(8) Kano, H.; Saito, T.; Kobayashi, T. Observation of Herzberg-Teller-type wave packetmotion in porphyrin J-aggregates studied by sub-5-fs spectroscopy. J. Phys. Chem. A , , 3445–3453.(9) Kobayashi, T.; Wang, Z.; Otsubo, T. Classification of dynamic vibronic couplings invibrational real-time spectra of a thiophene derivative by few-cycle pulses. J. Phys.Chem. A , , 12985–12994.(10) Fidler, A. F.; Engel, G. S. Nonlinear spectroscopic theory of displaced harmonic os-cillators with differing curvatures: A correlation function approach. J. Phys. Chem. A , , 9444–9453.(11) Bizimana, L. A.; Carbery, W. P.; Gellen, T. A.; Turner, D. B. Signatures of Herzberg-Teller coupling in three-dimensional electronic spectroscopy. J. Chem. Phys. , ,084311.(12) Anda, A.; Abramaviˇcius, D.; Hansen, T. Two-dimensional electronic spectroscopy ofanharmonic molecular potentials. Phys. Chem. Chem. Phys. , , 1642–1652.(13) Zuehlsdorff, T. J.; Montoya-Castillo, A.; Napoli, J. A.; Markland, T. E.; Isborn, C. M.Optical spectra in the condensed phase: Capturing anharmonic and vibronic featuresusing dynamic and static approaches. J. Chem. Phys. , , 074111.(14) Zuehlsdorff, T. J.; Hong, H.; Shi, L.; Isborn, C. M. Nonlinear spectroscopy in the con-densed phase: The role of Duschinsky rotations and third order cumulant contributions. J. Chem. Phys. , , 044127.(15) Fuji, T.; Saito, T.; Kobayashi, T. Dynamical observation of Duschinsky rotation bysub-5-fs real-time spectroscopy. Chem. Phys. Lett. , , 324–330.1816) Bizimana, L. A.; Brazard, J.; Carbery, W. P.; Gellen, T.; Turner, D. B. Resolving molec-ular vibronic structure using high-sensitivity two-dimensional electronic spectroscopy. J. Chem. Phys. , , 164203.(17) Galestian Pour, A.; Lincoln, C. N.; Perl´ık, V.; ˇSanda, F.; Hauer, J. Anharmonic vi-brational effects in linear and two-dimensional electronic spectra. Phys. Chem. Chem.Phys. , , 24752–24760.(18) Zhu, R.; Zou, J.; Wang, Z.; Chen, H.; Weng, Y. Electronic state-resolved multimode-coupled vibrational wavepackets in oxazine 720 by two-dimensional electronic spec-troscopy. J. Phys. Chem. A , , 9333–9342.(19) Fumero, G.; Schnedermann, C.; Batignani, G.; Wende, T.; Liebel, M.; Bassolino, G.;Ferrante, C.; Mukamel, S.; Kukura, P.; Scopigno, T. Two-dimensional impulsively stim-ulated resonant Raman spectroscopy of molecular excited states. Phys. Rev. X , , 11051.(20) Mukamel, S. Dressed-cluster and hydrodynamic expansions for line broadening in simplefluids. Phys. Rev. A , , 617.(21) Mukamel, S. Principles of nonlinear optical spectroscopy , 1st ed.; Oxford UniversityPress: New York, 1999.(22) Segarra-Mart´ı, J.; Mukamel, S.; Garavelli, M.; Nenov, A.; Rivalta, I. Towards AccurateSimulation of Two-Dimensional Electronic Spectroscopy.
Top. Curr. Chem. , ,24.(23) Picchiotti, A.; Nenov, A.; Giussani, A.; Prokhorenko, V. I.; Miller, R. J. D.;Mukamel, S.; Garavelli, M. Pyrene, a Test Case for Deep-Ultraviolet Molecular Photo-physics. J. Phys. Chem. Lett. , , 3481–3487.1924) Mukamel, S. On the semiclassical calculation of molecular absorption and fluorescencespectra. J. Chem. Phys. , , 173–181.(25) Egorov, S. A.; Rabani, E.; Berne, B. J. Vibronic spectra in condensed matter: A com-parison of exact quantum mechanical and various semiclassical treatments for harmonicbaths. J. Chem. Phys. , , 1407–1422.(26) Egorov, S. A.; Rabani, E.; Berne, B. J. Nonradiative relaxation processes in condensedphases: Quantum versus classical baths. J. Chem. Phys. , , 5238–5248.(27) Shi, Q.; Geva, E. A derivation of the mixed quantum-classical Liouville equation fromthe influence functional formalism. J. Chem. Phys. , , 3393–3404.(28) McRobbie, P. L.; Hanna, G.; Shi, Q.; Geva, E. Signatures of nonequilibrium solvationdynamics on multidimensional spectra. Acc. Chem. Res. , , 1299–1309.(29) McRobbie, P. L.; Geva, E. A benchmark study of different methods for calculating one-and two-dimensional optical spectra. J. Phys. Chem. A , , 10425–10434.(30) Mollica, C.; Van´ıˇcek, J. Beating the Efficiency of Both Quantum and Classical Simula-tions with a Semiclassical Method. Phys. Rev. Lett. , , 214101.(31) Wehrle, M.; ˇSulc, M.; Van´ıˇcek, J. Time-Resolved Electronic Spectra with EfficientQuantum Dynamics Methods. Chimia , , 334–338.(32) ˇSulc, M.; Van´ıˇcek, J. Accelerating the calculation of time-resolved electronic spectrawith the cellular dephasing representation. Mol. Phys. , , 945–955.(33) Zambrano, E.; ˇSulc, M.; Van´ıˇcek, J. Improving the accuracy and efficiency of time-resolved electronic spectra calculations: Cellular dephasing representation with a pref-actor. J. Chem. Phys. , , 054109.(34) Schubert, A.; Engel, V. Two-dimensional vibronic spectroscopy of coherent wave-packetmotion. J. Chem. Phys. , , 104304.2035) Krˇcm´aˇr, J.; Gelin, M. F.; Domcke, W. Calculation of third-order signals via drivenSchr¨odinger equations: General results and application to electronic 2D photon echospectroscopy. Chem. Phys. , , 53–62.(36) Krˇcm´aˇr, J.; Gelin, M. F.; Domcke, W. Simulation of femtosecond two-dimensionalelectronic spectra of conical intersections. J. Chem. Phys. , , 074308.(37) Picconi, D.; Cina, J. A.; Burghardt, I. Quantum dynamics and spectroscopy of dihalo-gens in solid matrices. I. Efficient simulation of the photodynamics of the embeddedI Kr cluster using the G-MCTDH method. J. Chem. Phys. , , 064111.(38) Picconi, D.; Cina, J. A.; Burghardt, I. Quantum dynamics and spectroscopy of di-halogens in solid matrices. II. Theoretical aspects and G-MCTDH simulations of time-resolved coherent Raman spectra of Schr¨odinger cat states of the embedded I Kr cluster. J. Chem. Phys. , , 064112.(39) Beguˇsi´c, T.; Van´ıˇcek, J. On-the-fly ab initio semiclassical evaluation of third-orderresponse functions for two-dimensional electronic spectroscopy. J. Chem. Phys. , , 184110.(40) Matzkies, F.; Manthe, U. Accurate reaction rate calculations including internal and ro-tational motion: A statistical multi-configurational time-dependent Hartree approach. J. Chem. Phys. , , 88–96.(41) Manthe, U.; Huarte-Larra˜naga, F. Partition functions for reaction rate calculations:Statistical sampling and MCTDH propagation. Chem. Phys. Lett. , , 321–328.(42) Gelman, D.; Kosloff, R. Simulating dissipative phenomena with a random phase thermalwavefunctions, high temperature application of the Surrogate Hamiltonian approach. Chem. Phys. Lett. , , 129–138. 2143) Nest, M.; Kosloff, R. Quantum dynamical treatment of inelastic scattering of atoms ata surface at finite temperature: The random phase thermal wave function approach. J. Chem. Phys. , , 134711.(44) Lorenz, U.; Saalfrank, P. Comparing thermal wave function methods for multi-configuration time-dependent Hartree simulations. J. Chem. Phys. , , 044106.(45) Wang, L.; Fujihashi, Y.; Chen, L.; Zhao, Y. Finite-temperature time-dependent varia-tion with multiple Davydov states. J. Chem. Phys. , , 124127.(46) Werther, M.; Grossmann, F. Including temperature in a wavefunction description ofthe dynamics of the quantum Rabi model. J. Phys. A: Math. Theor. , 014001.(47) Suzuki, M. Thermo Field Dynamics in Equilibrium and Non-Equilibrium InteractingQuantum Systems.
J. Phys. Soc. Jap. , , 4483–4485.(48) Takahashi, Y.; Umezawa, H. Thermo Field Dynamics. Int. J. Mod. Phys. B , ,1755–1805.(49) Harsha, G.; Henderson, T. M.; Scuseria, G. E. Thermofield Theory for Finite-Temperature Coupled Cluster. J. Chem. Theory Comput. , , 6127–6136.(50) Harsha, G.; Henderson, T. M.; Scuseria, G. E. Thermofield theory for finite-temperaturequantum chemistry. J. Chem. Phys. , , 154109.(51) Shushkov, P.; Miller, T. F. Real-time density-matrix coupled-cluster approach for closedand open systems at finite temperature. J. Chem. Phys. , , 134107.(52) Borrelli, R.; Gelin, M. F. Quantum electron-vibrational dynamics at finite temperature:Thermo field dynamics approach. J. Chem. Phys. , , 224101.(53) Borrelli, R.; Gelin, M. F. Simulation of Quantum Dynamics of Excitonic Systems atFinite Temperature: An efficient method based on Thermo Field Dynamics. Sci. Rep. , , 9127. 2254) Gelin, M. F.; Borrelli, R. Thermal Schr¨odinger Equation: Efficient Tool for Simulationof Many-Body Quantum Dynamics at Finite Temperature. Annalen der Physik , , 1700200.(55) Chen, L.; Zhao, Y. Finite temperature dynamics of a Holstein polaron: The thermo-field dynamics approach. J. Chem. Phys. , , 214102.(56) Beguˇsi´c, T.; Van´ıˇcek, J. On-the-fly ab initio semiclassical evaluation of vibronic spectraat finite temperature. J. Chem. Phys. , , 024105.(57) Heller, E. J. Time-dependent approach to semiclassical dynamics. J. Chem. Phys. , , 1544–1555.(58) Wehrle, M.; ˇSulc, M.; Van´ıˇcek, J. On-the-fly Ab Initio Semiclassical Dynamics: Identi-fying Degrees of Freedom Essential for Emission Spectra of Oligothiophenes. J. Chem.Phys. , , 244114.(59) Wehrle, M.; Oberli, S.; Van´ıˇcek, J. On-the-fly ab initio semiclassical dynamics of floppymolecules: Absorption and photoelectron spectra of ammonia. J. Phys. Chem. A , , 5685.(60) Schlau-Cohen, G. S.; Ishizaki, A.; Fleming, G. R. Two-dimensional electronic spec-troscopy and photosynthesis: Fundamentals and applications to photosynthetic light-harvesting. Chem. Phys. , , 1–22.(61) Khalil, M.; Demird¨oven, N.; Tokmakoff, A. Obtaining absorptive line shapes in two-dimensional infrared vibrational correlation spectra. Phys. Rev. Lett. , , 4.(62) Johnson, P. J. M.; Koziol, K. L.; Hamm, P. Intrinsic phasing of heterodyne-detectedmultidimensional infrared spectra. Opt. Express , , 2928.(63) Meyer, H.-D., Gatti, F., Worth, G. A., Eds. Multidimensional Quantum Dynamics:MCTDH Theory and Applications , 1st ed.; Wiley-VCH: Weinheim, 2009.2364) Roulet, J.; Choi, S.; Van´ıˇcek, J. Efficient geometric integrators for nonadiabatic quan-tum dynamics. II. The diabatic representation.
J. Chem. Phys. , , 204113.(65) Choi, S.; Van´ıˇcek, J. Efficient geometric integrators for nonadiabatic quantum dynam-ics. I. The adiabatic representation. J. Chem. Phys. , , 204112.(66) Choi, S.; Van´ıˇcek, J. A time-reversible integrator for the time-dependent Schr¨odingerequation on an adaptive grid. J. Chem. Phys. , , 234102.(67) Tatchen, J.; Pollak, E. Semiclassical on-the-fly computation of the S → S absorptionspectrum of formaldehyde. J. Chem. Phys. , , 041103.(68) Conte, R.; Ceotto, M. In Quantum chemistry and dynamics of excited states ;Gonz´alez, L., Lindh, R., Eds.; John Wiley & Sons, Ltd, 2020; Chapter 19, pp 595–628.(69) Buchholz, M.; Grossmann, F.; Ceotto, M. Mixed semiclassical initial value represen-tation time-averaging propagator for spectroscopic calculations.
J. Chem. Phys. , , 094102.(70) Buchholz, M.; Grossmann, F.; Ceotto, M. Application of the mixed time-averagingsemiclassical initial value representation method to complex molecular spectra. J. Chem. Phys. , , 164110.(71) Curchod, B. F. E.; Mart´ınez, T. J. Ab Initio Nonadiabatic Quantum Molecular Dy-namics. Chem. Rev. , , 3305–3336.(72) Makhov, D. V.; Glover, W. J.; Martinez, T. J.; Shalashilin, D. V. Ab initio multiplecloning algorithm for quantum nonadiabatic molecular dynamics. J. Chem. Phys. , , 054110. 2473) Makhov, D. V.; Symonds, C.; Fernandez-Alberti, S.; Shalashilin, D. V. Ab initio quan-tum direct dynamics simulations of ultrafast photochemistry with MulticonfigurationalEhrenfest approach. Chem. Phys. , , 200–218.(74) Thompson, A. L.; Mart´ınez, T. J. Time-resolved photoelectron spectroscopy from firstprinciples: Excited state dynamics of benzene. Faraday Discuss. , , 293–311.(75) ˇSulc, M.; Hern´andez, H.; Mart´ınez, T. J.; Van´ıˇcek, J. Relation of exact Gaussian basismethods to the dephasing representation: Theory and application to time-resolvedelectronic spectra. J. Chem. Phys. , , 034112.(76) Zimmermann, T.; Van´ıˇcek, J. Efficient on-the-fly ab initio semiclassical method for com-puting time-resolved nonadiabatic electronic spectra with surface hopping or Ehrenfestdynamics. J. Chem. Phys. , , 134102.(77) Worth, G. A.; Robb, M. A.; Burghardt, I. A novel algorithm for non-adiabatic directdynamics using variational Gaussian wavepackets. Faraday Discuss. , , 307–323.(78) Richings, G. W.; Polyak, I.; Spinlove, K. E.; Worth, G. A.; Burghardt, I.; Lasorne, B.Quantum Dynamics Simulations Using Gaussian Wavepackets: the vMCG Method. Int. Rev. Phys. Chem. , , 269–308.(79) Bonfanti, M.; Petersen, J.; Eisenbrandt, P.; Burghardt, I.; Pollak, E. Computationof the S1 S0 vibronic absorption spectrum of formaldehyde by variational Gaussianwavepacket and semiclassical IVR methods. J. Chem. Theory Comput. , , 5310–4323.(80) Polyak, I.; Richings, G. W.; Habershon, S.; Knowles, P. J. Direct quantum dynam-ics using variational Gaussian wavepackets and Gaussian process regression. J. Chem.Phys. , , 041101. 2581) Chen, L.; Sun, K.; Shalashilin, D. V.; Gelin, M. F.; Zhao, Y. Efficient simulation of time-and frequency-resolved four-wave-mixing signals with a multiconfigurational Ehrenfestapproach. J. Chem. Phys. , , 054105.(82) Van´ıˇcek, J.; Beguˇsi´c, T. In Molecular Spectroscopy and Quantum Dynamics ; Mar-quardt, R., Quack, M., Eds.; Elsevier, 2021; pp 199–229.(83) Lasser, C.; Lubich, C. Computing quantum dynamics in the semiclassical regime.
ActaNumerica , , 229–401.(84) Rohrdanz, M. A.; Cina, J. A. Probing intermolecular communication via lattice phononswith time-resolved coherent anti-Stokes Raman scattering. Mol. Phys. , , 1161–1178.(85) Beguˇsi´c, T.; Roulet, J.; Van´ıˇcek, J. On-the-fly ab initio semiclassical evaluation of time-resolved electronic spectra. J. Chem. Phys. , , 244115.(86) Beguˇsi´c, T.; Cordova, M.; Van´ıˇcek, J. Single-Hessian thawed Gaussian approximation. J. Chem. Phys. , , 154117.(87) Beer, M.; Longuet-Higgins, H. C. Anomalous light emission of azulene. J. Chem. Phys. , , 1390–1391.(88) Prlj, A.; Beguˇsi´c, T.; Zhang, Z. T.; Fish, G. C.; Wehrle, M.; Zimmermann, T.; Choi, S.;Roulet, J.; Moser, J.-E.; Van´ıˇcek, J. Semiclassical Approach to Photophysics BeyondKasha’s Rule and Vibronic Spectroscopy Beyond the Condon Approximation. The Caseof Azulene. J. Chem. Theory Comput. , , 2617–2626.(89) Xin, H.; Gao, X. Application of azulene in constructing organic optoelectronic materials:New tricks for an old dog. Chempluschem , , 945–956.(90) Dierksen, M.; Grimme, S. Density functional calculations of the vibronic structure ofelectronic absorption spectra. J. Chem. Phys. , , 3544–3554.2691) Niu, Y.; Peng, Q.; Deng, C.; Gao, X.; Shuai, Z. Theory of excited state decays andoptical spectra: Application to polyatomic molecules. J. Phys. Chem. A , ,7817–7831.(92) Nenov, A.; Giussani, A.; Fingerhut, B. P.; Rivalta, I.; Dumont, E.; Mukamel, S.; Gar-avelli, M. Spectral lineshapes in nonlinear electronic spectroscopy. Phys. Chem. Chem.Phys. , , 30925–30936.(93) Farfan, C. A.; Turner, D. B. Interference among multiple vibronic modes in two-dimensional electronic spectroscopy. Mathematics , , 157.(94) Kats, D.; Sch¨utz, M. A multistate local coupled cluster CC2 response method based onthe Laplace transform. J. Chem. Phys. , , 124117.(95) Lederm¨uller, K.; Kats, D.; Sch¨utz, M. Local CC2 response method based on the Laplacetransform: Orbital-relaxed first-order properties for excited states. J. Chem. Phys. , , 084111.(96) Lederm¨uller, K.; Sch¨utz, M. Local CC2 response method based on the Laplace trans-form: Analytic energy gradients for ground and excited states. J. Chem. Phys. , , 164113.(97) Sch¨utz, M. Oscillator strengths, first-order properties, and nuclear gradients for localADC(2). J. Chem. Phys. , upporting Information:Finite-temperature, Anharmonicity, andDuschinsky Effects on the Two-dimensionalElectronic Spectra from Ab Initio Thermo-fieldGaussian Wavepacket Dynamics Tomislav Beguˇsi´c ∗ and Jiˇr´ı Van´ıˇcek ∗ Laboratory of Theoretical Physical Chemistry, Institut des Sciences et Ing´enierie Chimiques,Ecole Polytechnique F´ed´erale de Lausanne (EPFL), CH-1015, Lausanne, Switzerland
E-mail: tomislav.begusic@epfl.ch; jiri.vanicek@epfl.ch
S-1 a r X i v : . [ phy s i c s . c h e m - ph ] F e b Rephasing and non-rephasing signals for a three-level system within the Born-Oppenheimer approx-imation
Since the Eqs. (6)–(8) of the main text are rarely found in that form, let us summarizehow they can be obtained from the expressions that are easily found in the literature. Thethird-order response function
S1,S2 R (3) ( t , t , t ) = X α =1 [ R α ( t , t , t ) − R α ( t , t , t ) ∗ ] (1)can be expressed in terms of correlation functions R ( t , t , t ) = C ( t , t , t + t + t ) , (2) R ( t , t , t ) = C ( t + t , t , t + t ) , (3) R ( t , t , t ) = C ( t , t + t , t ) , (4) R ( t , t , t ) = C ( − t , − t , t ) , (5)where C ( τ a , τ b , τ c ) = Tr[ ˆ ρ ˆ µ e i ˆH τ a / ~ ˆ µ e i ˆH τ b / ~ ˆ µ e − i ˆH τ c / ~ ˆ µ e − i ˆH ( τ a + τ b − τ c ) / ~ ] . (6)In Eq. (6), we use the bold font to denote S × S matrices representing operators act-ing on the electronic subspace and hat ˆ to denote operators acting on the nuclear sub-space; the trace is taken over both nuclear and electronic degrees of freedom. Hence, ˆ ρ = exp( − β ˆH ) / Tr[exp( − β ˆH )] is the full molecular density matrix, ˆH is the molecularHamiltonian, and ˆ µ is the dipole moment matrix. Equations (2)–(6) can be easily comparedto those found in the literature, for example, with Eqs. (3.5a)–(3.6) of Ref. S2.We now consider only three electronic states and invoke the following approximations:S-2. Born-Oppenheimer approximation: ˆH = diag( ˆ H , ˆ H , ˆ H ).2. Resonance condition: Only transitions between states 1 and 2 or 2 and 3 are allowed: ˆ µ = µ µ µ µ . (7)3. The system is initially in the ground electronic state, i.e., the temperature effects on theelectronic subspace are neglected: ˆ ρ = diag( ˆ ρ, , ρ = exp( − β ˆ H ) / Tr[exp( − β ˆ H )].Then, we take the trace over the electronic subspace in Eq. (6) to arrive at C ( τ a , τ b , τ c ) ≈ C ( τ a , τ b , τ c ) + C ( τ a , τ b , τ c ) , (8)where C and C are defined according to Eq. (8) of the main text.Finally, as discussed in Ref. S2, under the rotating-wave approximation, the rephasingand non-rephasing signals are given by S R ( t , t , t ) = R ( t , t , t ) + R ( t , t , t ) − R ( t , t , t ) ∗ , (9) S NR ( t , t , t ) = R ( t , t , t ) + R ( t , t , t ) − R ( t , t , t ) ∗ . (10)Although each R α splits into two terms of the right-hand side of Eq. (8), only one of thoseterms remains under the rotating-wave approximation. Therefore, Eqs. (9) and (10) lead toEqs. (6) and (7) of the main text. In the main text, we have introduced the fictitious Hilbert space with an orthonormal basis {| ˜ k i} , which is related to the physical system (with the orthonormal basis {| k i} ) throughS-3he conjugation rules S3 ( ˆ A ˆ B )˜= ˆ˜ A ˆ˜ B, (11)( ˆ A | α i )˜= ˆ˜ A | ˜ α i , (12)( a ˆ A + b ˆ B )˜= a ∗ ˆ˜ A + b ∗ ˆ˜ B, (13)( a | α i + b | β i )˜= a ∗ | ˜ α i + b ∗ | ˜ β i , (14)where the operators ( ˆ A , ˆ B ) act on and states ( | α i , | β i ) belong to the original (without tilde)or fictitious (with tilde) Hilbert spaces; a and b are complex numbers.Let | α i = P k a k | k i . Then, | ˜ α i = P k a ∗ k | ˜ k i , h ˜ α | ˜ β i = X k a k b ∗ k = h β | α i , (15)and h ˜ α | ˆ˜ A | ˜ α i = h ˜ α | ˆ˜ A ˜ α i = h ˆ Aα | α i = h α | ˆ A † | α i , (16)where in the second step of Eq. (16) we used Eqs. (12) and (15).We now complete the proof of the main text, namely the step connecting Eqs. (17) and(18) of the main text, as follows: h ˜ k | e − i ˆ˜ H ( τ a + τ b − τ c ) / ~ | ˜ k i = h ˜ k | h e i ˆ H ( τ a + τ b − τ c ) / ~ i ˜ | ˜ k i (17)= h k | h e i ˆ H ( τ a + τ b − τ c ) / ~ i † | k i (18)= h k | e − i ˆ H ( τ a + τ b − τ c ) / ~ | k i , (19)where in (17) we used the conjugation rules (11) and (13), in going from (17) to (18) weused Eq. (16), and in the last step we used the fact that the Hamiltonian is Hermitian.S-4 Thermo-field dynamics expression for C ( τ a , τ b , τ c ) be-yond the Born-Oppenheimer approximation In the main text, we invoked the Born–Oppenheimer approximation so that the thermo-field dynamics expression could be easily combined with the thawed Gaussian wavepacketdynamics. Of course, in many systems of interest, such as excitonic or multichromophoricsystems, or photochemically active molecules, the coupling between electronic states mustbe included in the calculations. Here, we present the result beyond the Born-Oppenheimerapproximation.As in Eq. (11) of the main text, we define the thermal vacuum, now in the extendedmolecular (electronic and vibrational) Hilbert space: | ¯ ( β ) i = ˆ ρ / X k | k ˜ k i . (20)In other words, the state | ¯ ( β ) i is described not only by a doubled number of vibrationalcoordinates, but also a doubled number of electronic degrees of freedom. Following similarsteps as in the derivation of the main text, we can rewrite Eq. (6) (see Sec. 1 of the SupportingInformation) as C ( τ a , τ b , τ c ) = h ¯ φ τ b ,τ a | ¯ φ ,τ c i , (21)where | ¯ φ τ,t i = e − i ˆ¯ H τ/ ~ ˆ µ e − i ˆ¯ H t/ ~ ˆ µ | ¯ ( β ) i (22)and ˆ¯ H = ˆ H − ˆ˜ H . The result (21) accounts for coupled electronic states and thermal popula-tion of excited electronic states. Importantly, Eq. (21) justifies the thermo-field wavepacketpicture even beyond the Born–Oppenheimer approximation.S-5 Rephasing and non-rephasing two-dimensional spec-tra
Figure S1: Rephasing (top) and non-rephasing (bottom) contributions to the two-dimensional spectra from Fig. 2 of the main text at zero temperature (left) and at T = 300 K(right). S-6igure S2: Rephasing (top) and non-rephasing (bottom) contributions to the two-dimensional spectra from Fig. 4 of the main text, computed with the on-the-fly ab initiosingle-Hessian thawed Gaussian approximation (“Anharmonic”, left), harmonic approxima-tion (middle), and the displaced harmonic oscillator (DHO) model (right).S-7 Two-dimensional spectra at t > Figure S3: Two-dimensional electronic spectra [Eq. (3) of the main text] at delay times t = 0 , ,
10 fs, computed with the on-the-fly ab initio single-Hessian thawed Gaussian ap-proximation (“Anharmonic”) at 0 K (first column) and 300 K (second column), harmonicapproximation at 300 K (third column), and the displaced harmonic oscillator (DHO) at300 K model (fourth column). Each spectrum shows the sum of the ground-state bleach andstimulated emission terms [first two terms on the right-hand sides of Eqs. (6) and (7) ofthe main text] corresponding to the S –S electronic transition in azulene. At zero delay,the ground-state bleach and stimulated emission contributions are indistinguishable. How-ever, already after 5 or 10 fs, the stimulated emission signal moves away from the diagonal,whereas the ground-state bleach spectrum remains close to the diagonal. At later t de-lays (see Fig. S4), the stimulated emission starts returning towards the diagonal, reflectingcoherent wavepacket dynamics in the excited electronic state. Indeed, the period of about20–25 fs corresponds to the faster displaced modes (see Table S1). However, the recurrenceis incomplete due to the slower dynamics in the other highly displaced modes.S-8igure S4: Same as Fig. S3 but for delay times t = 15 , ,
25 fs.Table S1: The most displaced modes of azulene. Ground-state frequencies ω i of the normalmodes are reported in terms of the wavenumber. The dimensionless displacements are definedas ∆ i = | p ω i / ~ q ,i | , where q is the excited-state minimum geometry expressed in the mass-scaled normal mode coordinates of the ground electronic state and q ,i is its component alongmode i . Only the modes with ∆ i > .
25 are shown.Wavenumber / cm − Displacement1660 0 . . . . . . . | S R ( t , , t ) | at short t and t times Figure S5: First 6 fs of | S R ( t , , t ) | (see Fig. 4 of the main text) computed with the on-the-fly ab initio thawed Gaussian approximation (“Anharmonic”), harmonic approximation,and the displaced harmonic oscillator (DHO) model.S-10 (A)symmetry of the rephasing and non-rephasingsignals in the time domain Figure S6: Absolute value of the non-rephasing time-domain signal ( | S NR ( t , , t ) | ) at zerotime delay ( t = 0), including only the ground-state bleach and stimulated emission terms[first two terms on the right-hand side of Eq. (7) of the main text], computed with the on-the-fly ab initio single-Hessian thawed Gaussian approximation (“Anharmonic”, left), har-monic approximation (middle), and the displaced harmonic oscillator (DHO) model (right)at 300 K.In Fig. S6, we show that the absolute value of the non-rephasing signal S NR ( t , , t ) at t = 0is symmetric with respect to the diagonal t = t for all three—anharmonic, harmonic, andthe displaced harmonic oscillator—models, whereas the rephasing signal | S R ( t , , t ) | (Fig. 4of the main text) exhibits this symmetry only for the displaced harmonic oscillator model.Indeed, if we assume the Condon approximation (ˆ µ ≈ µ = const.) and include only theground-state bleach and stimulated emission contributions, the non-rephasing signal at t = 0S-11an be rewritten as a function of t + t : S NR ( t , , t ) = C (0 , t , t + t ) + C ( − t , , t ) (23)= 2 | µ | Tr[ ˆ ρe i ˆ H ( t + t ) / ~ e − i ˆ H ( t + t ) / ~ ] (24)= S NR ( t , , t ) (25)and is, therefore, symmetric with respect to the exchange of t and t . In going from (23)to (24), we used the definition of C [Eq. (8) of the main text], the Condon approximation,the fact that ˆ ρ commutes with the ground-state vibrational Hamiltonian ˆ H , and the cyclicproperty of the trace. The same, however, does not hold for the rephasing contribution: S R ( t , , t ) = 2 C ( t , t , t ) (26)= 2 | µ | Tr[ ˆ ρe i ˆ H t / ~ e i ˆ H t / ~ e − i ˆ H t / ~ e − i ˆ H t / ~ ] (27) = 2 | µ | Tr[ ˆ ρe i ˆ H t / ~ e i ˆ H t / ~ e − i ˆ H t / ~ e − i ˆ H t / ~ ] (28)= S R ( t , , t ) . (29)Yet, within the displaced harmonic oscillator model, the symmetry is recovered for theabsolute value of S R ( t , , t ). This is discussed in the following Section. Within the second-order cumulant method, S1 assuming only two electronic states (i.e., ne-glecting the excited-state absorption) and the Condon approximation (ˆ µ ≈ µ = const.), theS-12omponents R α ( t , t , t ) of the third-order response function (1) are approximated as R ( t , t , t ) = | µ | e − iω t − iω t e − g ( t ) ∗ − g ( t ) − f + ( t ,t ,t ) , (30) R ( t , t , t ) = | µ | e iω t − iω t e − g ( t ) ∗ − g ( t ) ∗ + f + ( t ,t ,t ) ∗ , (31) R ( t , t , t ) = | µ | e iω t − iω t e − g ( t ) − g ( t ) ∗ + f − ( t ,t ,t ) ∗ , (32) R ( t , t , t ) = | µ | e − iω t − iω t e − g ( t ) − g ( t ) − f − ( t ,t ,t ) , (33)where ω = 1 ~ Tr[( ˆ V − ˆ V ) ˆ ρ ] (34)is the thermally averaged electronic energy gap, ˆ V i are the potential energy operators in theground ( i = 1) or excited ( i = 2) electronic states, g ( t ) = 1 ~ Z t dτ Z τ dτ Tr[ e i ˆ H τ / ~ ∆ ˆ V e − i ˆ H τ / ~ ∆ ˆ V ˆ ρ ] , (35)∆ ˆ V = ˆ V − ˆ V − ~ ω , and f − ( t , t , t ) = g ( t ) − g ( t + t ) − g ( t + t ) + g ( t + t + t ) , (36) f + ( t , t , t ) = g ( t ) ∗ − g ( t + t ) ∗ − g ( t + t ) + g ( t + t + t ) . (37)For simple systems, such as a pair of harmonic potentials, exact quantum-mechanical g ( t ) can be found analytically. S1,S4
The second-order cumulant expansion is exact only inthe special case of the displaced harmonic (Brownian) oscillator model. For general poten-tials, efficient classical approximations to g ( t ) can be employed, at the price of introducingadditional approximations in the evaluation of the response function (1). Remarkably, it iseasy to show that, no matter how g ( t ) is computed, the symmetry relationship | R α ( t , t , t ) | = | R α ( t , t , t ) | , α = 1 , . . . , , (38)S-13olds for the four second-order cumulant expressions (30)–(33). For example, for R we have | R ( t , t , t ) | = R ( t , t , t ) ∗ R ( t , t , t ) (39)= e − g ( t ) ∗ − g ( t )+ f − ( t ,t ,t ) e − g ( t ) − g ( t ) ∗ + f − ( t ,t ,t ) ∗ (40)= e − g ( t )+ g ( t ) − f − ( t ,t ,t )] (41)= e − g ( t )+ g ( t ) − f − ( t ,t ,t )] (42)= | R ( t , t , t ) | , (43)where (42) follows from (41) because f − ( t , t , t ) = f − ( t , t , t ) . (44)Similar procedure can be followed for other R α functions, where we also need to use therelation: Re f + ( t , t , t ) = Re f + ( t , t , t ) . (45)Since the exact R α for the displaced harmonic oscillator model can be written in the formof Eqs. (30)–(33), the proof holds in this special case as well.Equation (38) implies that the symmetry with respect to exchanging t and t is presentfor arbitrary delay t . However, in experiments, one cannot measure individual componentsof the response function but only their linear combination. For example, the signal measuredin the rephasing phase-matching direction is S2 S R ( t , t , t ) = R ( t , t , t ) + R ( t , t , t ) . (46)S-14hen, | S R ( t , t , t ) | = | R ( t , t , t ) | + | R ( t , t , t ) | + 2Re[ R ( t , t , t ) ∗ R ( t , t , t )] (47) = | S R ( t , t , t ) | (48)because R ( t , t , t ) ∗ R ( t , t , t ) = e − g ( t ) ∗ − g ( t )+ f − ( t ,t ,t ) e − g ( t ) ∗ − g ( t ) ∗ + f + ( t ,t ,t ) ∗ (49)= e − g ( t ) ∗ − g ( t )]+ f − ( t ,t ,t )+ f + ( t ,t ,t ) ∗ (50)= e − g ( t ) ∗ − g ( t )]+2 g ( t ) − g ( t + t ) − g ( t + t )]+2Re[ g ( t + t + t )] (51) = e − g ( t ) ∗ − g ( t )]+2 g ( t ) − g ( t + t ) − g ( t + t )]+2Re[ g ( t + t + t )] (52)= e − g ( t ) ∗ − g ( t )]+ f − ( t ,t ,t )+ f + ( t ,t ,t ) ∗ (53)= R ( t , t , t ) ∗ R ( t , t , t ) . (54)Nevertheless, the equality is restored at t = 0 because R ( t , , t ) = R ( t , , t ): | S R ( t , , t ) | = 4 | R ( t , , t ) | = 4 | R ( t , , t ) | = | S R ( t , , t ) | . (55) References (S1) Mukamel, S.
Principles of nonlinear optical spectroscopy , 1st ed.; Oxford UniversityPress: New York, 1999.(S2) Schlau-Cohen, G. S.; Ishizaki, A.; Fleming, G. R. Two-dimensional electronic spec-troscopy and photosynthesis: Fundamentals and applications to photosynthetic light-harvesting.
Chem. Phys. , , 1–22.(S3) Suzuki, M. Thermo Field Dynamics in Equilibrium and Non-Equilibrium InteractingQuantum Systems. J. Phys. Soc. Jap. , , 4483–4485.S-15S4) Zuehlsdorff, T. J.; Hong, H.; Shi, L.; Isborn, C. M. Nonlinear spectroscopy in the con-densed phase: The role of Duschinsky rotations and third order cumulant contributions. J. Chem. Phys. ,153