Finite-Temperature, Anharmonicity, and Duschinsky Effects on the Two-Dimensional Electronic Spectra from Ab Initio Thermo-Field Gaussian Wavepacket Dynamics.

Accurate description of finite-temperature vibrational dynamics is indispensable in the computation of two-dimensional electronic spectra. Such simulations are often based 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 of thermo-field dynamics to derive an exact finite-temperature expression that lends itself to an intuitive wavepacket-based interpretation. Furthermore, an efficient method for computing finite-temperature two-dimensional spectra is obtained by combining the exact thermo-field dynamics approach with the thawed Gaussian approximation for the wavepacket dynamics, which is exact for any displaced, distorted, and Duschinsky-rotated harmonic potential but also accounts partially for anharmonicity effects in general potentials. Using this new method, we directly relate a symmetry breaking of the two-dimensional signal to the deviation from the conventional Brownian oscillator picture.

Multidimensional optical spectroscopy is an emerging experimental method for studying molecular photochemistry and photophysics, but its further development and the interpretation of new experiments rely heavily on theoretical modeling.[1−7] To this end, a number of theoretical methods[8−15] were developed to account for typical vibrational–electronic effects occurring in molecular systems, such as anharmonicity, different curvatures of the ground- and excited-state potential energy surfaces, or mode–mode mixing (Duschinsky rotation).[16−20] In its original formulation, the second-order cumulant expansion[21−24] is exact only for the Brownian oscillator (i.e., displaced harmonic) model and cannot treat the intermode coupling in the excited state. Although this basic molecular model shaped our understanding of steady-state, ultrafast, and multidimensional electronic spectroscopy in the past decades, it is inadequate for many molecules that exhibit Duschinsky and anharmonicity effects.[12,14] Similar limitations are met when using the semiclassical phase averaging,[22,25] also known as the Wigner-averaged classical limit[26−30] or dephasing representation.[31−34] The recently developed third-order cumulant approach seems to overcome these limitations,[10,13,14] yet it is accurate only in systems with weakly coupled or distorted modes.[14]

In contrast, quantum dynamics methods[35−40] are well suited for describing the evolution of nuclear wavepackets but often neglect temperature effects. To avoid the impractical Boltzmann averaging over the initial states, a number of alternative strategies for including temperature in wavepacket-based methods have been proposed.[41−47] We turn to the so-called thermo-field dynamics,[48,49] which transforms the von Neumann evolution of a density matrix to a Schrödinger equation with a doubled number of degrees of freedom. This approach has only recently been introduced in chemistry for solving the electronic structure,[50−52] vibronic,[53−56] and spectroscopic[57] problems at finite temperature. Here, we show how it could be used to compute two-dimensional vibronic spectra. The finite-temperature treatment is combined with the thawed Gaussian approximation,[58] an efficient first-principles[59,60] method for wavepacket propagation, and applied to the stimulated emission and ground-state bleach signals of azulene.

In two-dimensional spectroscopy, a nonlinear time-dependent polarization[22,61]is induced in the sample through interaction with the electric field E(t) comprised of three light pulses centered at times –t2 – t1, −t2, and 0, where t1 is the delay between the first two pulses, t2 is the delay between the second and third pulses, and R(3)(t‴, t″, t′) is the third-order response function.[22] In a heterodyne detection scheme, the measured signal is[61]where ELO(t) is the fourth, local oscillator pulse centered at time t3 after the third pulse. The two-dimensional spectrum is obtained by scanning S(t3, t2, t1) as a function of the three time delays and Fourier transforming over t1 and t3. We focus on the absorptive two-dimensional spectrum[62,63]and assume ultrashort and nonoverlapping pulse approximations, where the rephasing and nonrephasing spectra[61]are defined throughandIn eq , Ĥ terms are the vibrational Hamiltonians corresponding to the ground (j = 1) and excited (j = 2 or 3) electronic states, ρ̂ = exp(−βĤ1)/Tr[exp(−βĤ1)] is the vibrational density operator at temperature T = 1/kBβ, and is the electronic transition dipole moment between electronic states i and j projected on the polarization unit vector ε⃗ of the external electric field. Correlation function C1 corresponds to the stimulated emission and ground-state bleach processes, while C3, which involves a higher excited electronic state, corresponds to the excited-state absorption (see section 1 of the Supporting Information). Although the excited-state absorption term involves, in general, a sum over several higher excited states (i ≥ 3), here, for the sake of brevity, we consider only one such state.

An intuitive physical interpretation of eq is available in the zero-temperature limit, where the density operator ρ̂ = |1,g⟩⟨1,g| is given in terms of the ground vibrational state |1,g⟩ of the ground electronic state. Then[40]whereĤ′ = Ĥ – ℏω1,, and ℏω1, = ⟨1,g|Ĥ1|1,g⟩. In Figure , we illustrate how eq is evaluated for stimulated emission contribution C1(t1 + t2, t3, t2 + t3) (eq ) to the rephasing signal (eq ). The bra nuclear wavepacket is first evolved for a time τ = t1 + t2 in the excited electronic state and then for a time τ = t3 in the ground state; the ket wavepacket is in the ground electronic state during t1 and evolves in the excited state for a time τ = t2 + t3. In general, during time delays t1 and t3, also known as coherence and detection times, the bra and ket wavepackets evolve on different potential energy surfaces; during the so-called population time t2, the two wavepackets are in the same electronic state: in the ground state for the ground-state bleach contribution and in the excited electronic state for the stimulated emission and excited-state absorption components.

Figure 1

Evolution of the bra (a, dotted line) and ket (b, solid line) wavepackets of eq for τ = t1 + t2, τ = t3, and τ = t2 + t3. Their overlap (c) is stimulated emission term C1(t1 + t2, t3, t2 + t3) (eq ) of rephasing signal SR(t3, t2, t1) (eq ).

We now address the question of whether it is possible to retain the simple wavepacket picture without 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 finite temperature to the evolution of a wave function with a doubled number of coordinates. In the thermo-field dynamics theory,[48] the thermal vacuum is defined aswhere |kk̃⟩ = |k⟩|k̃⟩ is the basis vector of the tensor-product space obtained from the physical (with basis {|k⟩}) and “fictitious” (with basis {|k̃⟩}) Hilbert spaces. We note that physical operators (denoted only by a hat ^, such as ρ̂ or μ̂) act only on the physical subspace. With these definitions, eq can be rewritten aswhereis the analogue of |ϕ(⟩ from eq is the Hamiltonian acting in the full, tensor-product space, and is the ground-state vibrational Hamiltonian acting in the fictitious space only. The proof of eq goes as follows:Equation is obtained from eq by inserting the definition (eq ) of , while eq results upon substituting eq for |0̅(β)⟩; in going from eq to eq , we used the fact that operators acting in different subspaces commute. In going from eq to eq , we used the conjugation rules relating the physical and fictitious spaces (see section 2 of the Supporting Information). The resolution of identity and commutation of ρ̂1/2 with Ĥ1 were used to obtain eq , and the definition and cyclic property of the trace to obtain eq .

Remarkably, the result (eq ) has exactly the same form as the zero-temperature expression (eq ) and can be interpreted as in Figure . It also allows finite-temperature effects to be included in regular wave function-based codes, by modifying only the definition of the initial state and the Hamiltonians under which this state is evolved. In section 3 of the Supporting Information, we prove that the same wavepacket picture can be justified even beyond the Born–Oppenheimer approximation, which was invoked implicitly in eqs –8. To avoid exponentially scaling exact quantum methods on precomputed potential energy surfaces[64−67] or computationally demanding multiple-trajectory[68−82] approaches, we propose using the simple, yet efficient, single-trajectory thawed Gaussian approximation, which can be interfaced with on-the-fly ab initio evaluation of potential energy information.[83]

Let us consider a Gaussian wavepacketwhere q and p are the real, D-dimensional expectation values of the position and momentum, respectively, A is a D × D complex symmetric matrix with a positive-definite imaginary part, γ is a complex scalar whose imaginary part ensures normalization of the wavepacket, and D is the number of coordinates. Within the thawed Gaussian approximation,[58] one replaces true potential energy V(q) by its local harmonic approximationabout the center q of the wavepacket, which leads to the following equations of motion for the Gaussian’s parameters:[58,84]where L = p·(2m)−1·p – V(q) is the Lagrangian along the trajectory (q, p) and m is the symmetric mass matrix. According to eqs –26, the position and momentum of the Gaussian wavepacket evolve classically, while matrix A depends on the Hessians along the classical trajectory. The described evolution of the Gaussian wavepacket is exact for a harmonic potential because the local Taylor expansion of eq becomes exact in this case. For more general, anharmonic potentials, the method is only approximate, but typically accurate for moderate anharmonicity and short times, which makes it practical in spectroscopic applications.[58−60,85,86] Although the thawed Gaussian propagation is not suited for nonadiabatic dynamics, it can treat accurately the effects that arise due to different force constants of the ground- and excited-state potential surfaces: mode distortion, i.e., the change in the frequency of a normal mode, and intermode coupling or Duschinsky rotation. The on-the-fly ab initio thawed Gaussian approximation, which uses electronic structure calculations to compute potential energies, gradients, and Hessians only when needed, was recently validated for the simulation of finite-temperature linear[57] and zero-temperature two-dimensional spectra.[40]

To construct the initial state, we approximate the ground-state potential energy surface by a harmonic potential and use the corresponding mass-scaled normal mode coordinates. Then, in the zero-temperature limit, the initial state ψ0(q) = ⟨q|1,g⟩ is a Gaussian (eq ) and D = F, where F is the number of vibrational degrees of freedom. In the thermo-field dynamics formulation, D = 2F, the initial state ψ̅0(q̅) = ⟨q̅|0̅(β)⟩ is also a Gaussian, and q̅ = (q, q̃) is the 2F-dimensional coordinate vector.[57] To solve the equations of motion in the finite-temperature picture, we need the potential energies, gradients, and Hessians in the extended 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 (57). Remarkably, the thermo-field dynamics under Hamiltonian (eq ) requires exactly the same classical trajectory, in electronic state j, as the conventional, zero-temperature thawed Gaussian propagation with Hamiltonian Ĥ.[57] No further ab initio evaluations are needed for the finite-temperature implementation, meaning that, within the thawed Gaussian approximation, the temperature effects can be included almost for free. The only difference in the computational cost is in solving the equations of motion with 2F rather than F coordinates, which is approximately 23 = 8 times more expensive due to the roughly cubic scaling of the involved matrix operations, including matrix–matrix multiplication and matrix inverse. This cost is, however, negligible compared to the cost of electronic structure calculations.

Formally, the propagation of the wavepacket according to eqs –26 requires not only the potential energies and gradients but also the Hessians at each step of the dynamics. In this work, we employed the single-Hessian method,[87] which further approximates V″(q) ≈ V″(qref) in eq , where qref is a reference geometry at which the Hessian of the excited-state potential surface is evaluated once and reused during the excited-state dynamics. Because the 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 (87), the accuracy of this method was shown to be similar to that of the thawed Gaussian approximation. Here, we chose qref as the excited-state minimum. The ground-state potential surface was assumed to be harmonic in all simulations.

To analyze the effects of the excited-state anharmonicity, we compare the anharmonic calculations, based on the on-the-fly single-Hessian thawed Gaussian approximation for the excited-state propagation, with the harmonic model (also called the generalized Brownian oscillator model), where the excited-state potential surface is approximated by a harmonic potential fitted to the surface at its minimum (so-called adiabatic harmonic or adiabatic Hessian scheme). In the mass-scaled normal mode coordinates of the ground state, the excited-state force constant is a symmetric, nondiagonal matrix, whose off-diagonal terms reflect intermode couplings, also known as Duschinsky mixing. To study the effects of the difference between the excited- and ground-state force constants on linear and two-dimensional spectra, we construct the displaced harmonic model (also called the Brownian oscillator model), where the excited-state force constant is approximated by the force constant in the ground electronic state. This model neglects mode distortion and Duschinsky effects. The two-dimensional spectra can be computed exactly with the thawed Gaussian propagation, as described above, for both harmonic and displaced harmonic oscillator models. Whereas the exact solution to the displaced harmonic oscillator model was known before in the form of the second-order cumulant expansion,[22] to the best of our knowledge, no method has been published for computing exactly the two-dimensional spectra of the global harmonic (or generalized Brownian oscillator) model.[14]

Azulene is a well-known example of a Kasha-violating molecule,[88] as it emits light from the second, rather than first, excited electronic state. This is due to the interplay of two factors:[89] (i) weak nonadiabatic coupling between states S1 and S2 and (ii) fast (≈1 ps) nonradiative decay from S1 to S0. These properties make azulene one of the key building blocks in the synthesis of novel optoelectronic materials.[90] Although nonadiabatic couplings between the ground and first excited states play an important role in the photoinduced dynamics of azulene,[89] they do not affect its vibrationally resolved S1 ← S0 absorption spectrum. Indeed, the linear absorption spectrum can be reproduced well using adiabatic, Born–Oppenheimer approaches that neglect nonadiabatic effects.[89,91,92] Here, we also ignore the nonadiabatic effects on the two-dimensional spectra, which we compute only at short t2 delay times. In the results, we focus on the ground-state bleach and stimulated emission contributions to the two-dimensional spectrum (the first two terms on the right-hand sides of eqs and 7); according to the oscillator strengths of the S1–S0 (0.009)[93−95] and S2–S1 (≈10–5)[94,95] transitions, the excited-state absorption is expected to be ∼3 orders of magnitude weaker.

In Figure (top), we compare linear absorption spectra simulated at 300 and 0 K with the experimental spectrum recorded at room temperature. One of the main effects of temperature is the broadening of the spectral features, which also affects the relative intensities of vibronic peaks, namely, those at 14300 and 15800 cm–1. These intensities are overestimated in the zero-temperature spectrum but corrected by the finite-temperature treatment.

Figure 2

S1 ← S0 absorption spectra of azulene (top) computed with the on-the-fly ab initio single-Hessian thawed Gaussian approximation at zero temperature (red, dashed) and at 300 K (blue, dotted), compared with the experimental spectrum (black, solid) recorded at room temperature in cyclohexane.[89] Absorptive two-dimensional electronic spectra (bottom) (eq ) at zero delay time (t2 = 0), computed at zero temperature (left) and 300 K (right). Each two-dimensional spectrum shows the sum of the ground-state bleach and stimulated emission terms (first two terms on the right-hand sides of eqs and 7) corresponding to the S1–S0 electronic transition in azulene. See Figure S1 for the rephasing and nonrephasing contributions to these spectra and Figures S3 and S4 for the spectra at delays t2 > 0.

A non-zero temperature has an even stronger effect on the two-dimensional spectrum (Figure , 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 linear spectrum, the temperature effects modify not only the resolution of the spectrum but also the relative intensities of the peaks. However, in contrast to the linear absorption spectrum, where these differences affect only a few peaks and could still be considered acceptable, the two-dimensional spectrum is strongly affected due to the increased complexity of spectral features.

To investigate the effects of anharmonicity, mode distortion, and mode–mode coupling, we first compare the linear absorption spectra computed using three models with different accuracies (see Figure ). The spectrum computed with the displaced harmonic oscillator model displays a highly regular intensity pattern, as opposed to the irregular intensities found in the experiment, and overestimates the frequency spacing between the peaks. The results are largely improved by including Duschinsky coupling and changes in the mode frequencies through the global harmonic model. However, the harmonic approximation suffers from an overly broad tail in the high-frequency region. This is further corrected by accounting for the anharmonicity effects with the on-the-fly thawed Gaussian approximation.

Figure 3

S1 ← S0 absorption spectra of azulene computed with the on-the-fly ab initio single-Hessian thawed Gaussian approximation (“Anharmonic”, blue, dotted), harmonic approximation (red, dashed), and displaced harmonic oscillator (DHO) model (green, dashed–dotted) at 300 K, compared with the experimental spectrum (black, solid) recorded at room temperature in cyclohexane.[89]

The corresponding two-dimensional spectra (Figure , top) exhibit similar differences, which we can conveniently analyze in the time domain [see Figure , bottom, for |SR(t3, 0, t1)| and Figure S6 for |SNR(t3, 0, t1)|]. The displaced harmonic oscillator model results in stronger recurrences after 45 fs (in t1, t3, or both t1 and t3) than the harmonic or anharmonic approaches. This translates into sharper peaks in the two-dimensional spectrum. The anharmonic spectrum extends less into the high-frequency region, compared to the harmonic and displaced harmonic oscillator models, because the thawed Gaussian propagation gives a slower initial decay (for t1 and t3 < 6 fs) in the time domain than the models that neglect anharmonicity (see Figure S5). Subtle differences between the harmonic and anharmonic excited-state dynamics affect the peak intensities in the region between 15000 and 18000 cm–1.

Figure 4

Absorptive two-dimensional electronic spectra (top) (eq ) at zero delay time (t2 = 0), computed with the on-the-fly ab initio single-Hessian thawed Gaussian approximation (“Anharmonic”, left), harmonic approximation (middle), and displaced harmonic oscillator (DHO) model (right) at 300 K. Each spectrum shows the sum of the ground-state bleach and stimulated emission terms (first two terms on the right-hand sides of eqs and 7) corresponding to the S1–S0 electronic transition in azulene. See Figure S2 for the rephasing and nonrephasing contributions to these spectra and Figures S3 and S4 for the spectra at delays t2 > 0. First 60 fs of |SR(t3, 0, t1)| (bottom) (eq ). See Figure S6 for |SNR(t3, 0, t1)|.

Interestingly, for the displaced harmonic oscillator model, |SR(t3, 0, t1)| is symmetric with respect to the diagonal (Figure , bottom right), which does not hold when mode distortion, rotation, and anharmonicity are included (Figure , bottom left and middle). We prove this analytically in sections 7 and 8 of the Supporting Information, where we also demonstrate that the asymmetry can appear only in rephasing signal |SR(t3, 0, t1)|. Moreover, we show that the (incorrect) symmetry of |SRDHO(t3, 0, t1)| with respect to the diagonal t1 = t3 is, more generally, imposed by the second-order cumulant approximation,[22] which is exact for the displaced harmonic oscillator model and is employed regularly to model two-dimensional spectra.[23,24,96,97] Hence, the second-order cumulant method cannot account for the asymmetry induced by the deviation from the displaced harmonic oscillator model. This erroneous qualitative behavior was difficult to study in the past, partly due to the absence of practical methods that could easily go beyond the second-order cumulants or Brownian oscillators.

In conclusion, we derived a general and exact expression for computing finite-temperature vibrationally resolved two-dimensional electronic spectra with wave function-based methods. The inclusion of temperature is the key to simulating spectra of larger systems or solvated molecules, due to the multitude of low-frequency modes that are thermally excited at room temperature. By combining the exact expression with the thawed Gaussian approximation, we developed a practical and efficient method for computing two-dimensional spectra beyond zero temperature and beyond the displaced harmonic oscillator model. With the help of the newly developed method, we identified an asymmetry in the time domain signal that could serve as evidence for the changes in mode frequencies, mode–mode coupling, or anharmonicity. This asymmetry cannot be described with the conventional and widely used second-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-order Laplace-transformed density-fitted local algebraic diagrammatic construction [LT-DF-LADC(2)] scheme,[98−101] as implemented in the Molpro 2015 package.[102] The cc-pVDZ basis set was used throughout (see ref (89)). We first evaluated the Hessians in the ground and excited states at the respective optimized geometries. Then, starting from the minimum of the ground state, an on-the-fly ab initio classical trajectory was evolved in the excited electronic state for 1130 steps with a time step of 8 au ≈ 0.19 fs (total time of ≈219 fs).

Linear spectra were computed by Fourier transforming the first 500 steps of the wavepacket autocorrelation function (see ref (83)). With regard to the simulation of two-dimensional spectra, t1 and t3 times were propagated up to ≈106 fs (500 steps); t2 delays ranged from 0 (results shown in the main text) to 25 fs (130 steps), in intervals of 5 fs or 26 steps. Condon approximation, which was justified for the S1 ← S0 absorption of azulene in ref (89), was employed. Gaussian broadening with a half-width at half-maximum of 90 cm–1 was used in both linear and two-dimensional spectra. Linear spectra were shifted in frequency and scaled in intensity to match at the maximum intensity peak of the experiment; two-dimensional spectra were shifted by the same frequency shifts as the linear absorption spectra and scaled according to the maximum of the fully absorptive two-dimensional spectrum (eq ).

Data supporting this publication can be found at 10.5281/zenodo.4552858.