Mean-trajectory approximation for electronic and vibrational-electronic nonlinear spectroscopy

Multidimensional vibrational spectroscopy reveals both structural and dynamical properties of molecular systems by probing interactions among nuclear degrees of freedom.^1–5 Numerically exact quantum dynamical evaluations of the relevant nonlinear response functions are generally impractical for large systems, motivating the development of approximate methods. These include treating a significant subsystem quantum mechanically and the rest of the system classically or semiclassically^6–15 or treating all degrees of freedom with a semiclassical approximation to quantum dynamics.^16–21 

Evaluation of nonlinear vibrational response functions with semiclassical propagators requires computing pairs of classical trajectories representing the ket and bra aspects of the density operator.²² Calculating quantum interferences from trajectory pairs requires high dimensional integration over strongly oscillatory integrands, limiting this approach to small systems. The mean-trajectory method^23–26 is based on approximate integration over differences between trajectory pairs, resulting in the computation of a single mean trajectory, subject to the semiclassical quantization of classical action variables. The optimized mean-trajectory (OMT) approximation^27–31 yields improved results relative to the original form of the mean-trajectory approximation and is derived more simply, in terms of a mapping between quantum double-sided Feynman diagrams (2FD)^4,32,33 and semiclassical OMT diagrams.

The essential feature of the OMT approximation for a single nuclear degree of freedom is the representation of the dynamics of the operator $|nl⟩⟨nr|$ with a classical trajectory with action $ℏ(n¯+1/2)$⁠, with $n¯=(nl+nr)/2$ the mean quantum number. Population dynamics are represented by a classical trajectory with action equal to a half-odd-integral multiple of ℏ, while a single quantum coherence maps to a trajectory with action equal to an integral multiple of ℏ. The OMT has been demonstrated to reproduce line broadening from both dissipation and pure dephasing for a chromophore interacting with a bath,²⁹ as well as vibrational coherence transfer between near-resonant chromophore modes.³⁰ 

In addition to the infrared version, multidimensional spectroscopy has proven its importance in electronic spectroscopy^34–38 and has been extended into a mixed vibrational-electronic regime,^39–45 combining infrared and optical excitation. Since the OMT is based on the semiclassical quantization of action variables for nuclear degrees of freedom, an extension to include electronic degrees of freedom may be based on describing electronic levels in terms of classical action-angle variables. The mapping of discrete states onto classical phase-space variables, both exactly and approximately, is well established.^46–53 This extended OMT can provide a consistent semiclassical framework for calculating purely vibrational, purely electronic, and mixed vibrational-electronic multidimensional spectra.

In Sec. II, we specify our model with the quantum Hamiltonian of a two-level electronic system coupled to a nuclear degree of freedom and construct an approximate corresponding classical Hamiltonian, expressible in terms of electronic and nuclear action and angle variables. The third-order optical response function in a measurement with electronic excitation is determined within the OMT in Sec. III. Section IV presents the OMT treatment of a measurement employing an initial pair of infrared pulses that excite vibrations, followed by an optical pulse that excites an electronic transition, followed in turn by the detection of an optical signal. The response functions for the purely electronic measurement and the vibrational-electronic experiment are evaluated in Sec. V for an analytically solvable model of displaced harmonic surfaces. OMT calculations are compared to exact results, derived in  Appendices A and  B. Conclusions are summarized in Sec. VI.

Interaction of matter with electromagnetic fields of wavevectors $k→1$⁠, $k→2$⁠, and $k→3$ generates four-wave mixing signals with wavevectors $±k→1±k→2±k→3$⁠. The third-order response function component^4,32 associated with signal wavevector $αk→1+βk→2+γk→3$ is

The radiation-matter interaction is assumed to be of the electric dipole form with $μ^$ the dipole operator. Superscripts label dipole components that create or destroy excitations, and δ is chosen so that there are two dipole components of each sense. The tetradic propagator³² of the material system in the absence of radiation is denoted $G(t)$⁠, and $ρ^$ is the initial density operator. An absorptive two-dimensional spectrum⁴ is determined from Fourier transforms with respect to t₁ and t₃ of the rephasing or echo amplitude R₊₊₋(t₁, t₂, t₃) and the particular non-rephasing amplitude R₊₋₊(t₁, t₂, t₃), so these are the two signals considered explicitly here.

The material system has the vibronic Hamiltonian

with electronic states labeled 0 and 1 and Ω the electronic transition frequency. For simplicity, a single nuclear degree of freedom with phase-space variables z = (q, p) is treated. The electronically adiabatic nuclear Hamiltonian for state j is $H^j=p^2/2m+Vj(q^)$⁠. For interactions with optical fields resonant with an electronic transition, the dipole in Eq. (1) has the form $μ^el=μ^el++μ^el−$⁠, with $μel−=(μel+)†$ and

with m₀ real-valued. We make the Condon approximation in treating $μ^el$ as a purely electronic operator. For interactions with infrared radiation resonant with a purely vibrational transition, $μ^$ in Eq. (1) is taken to be proportional to the nuclear coordinate $q^$⁠, with the coefficient of proportionality suppressed.

Our semiclassical approximation to nonlinear response functions is based on representing the electronic states with continuously valued classical degrees of freedom.^46–53 This representation is based on the Heisenberg correspondence relation^54–57 connecting a quantum matrix element to a Fourier integral of a classical function of action and angle variables. In particular, for a two-dimensional state space, as in the present model with nuclear degrees of freedom suppressed, a matrix element of a time-evolved operator $Anlnr(t)≡⟨nl|A^(t)|nr⟩$ is written as

Here $ℏn$ and ϕ are the classical action and angle variables representing the electronic degree of freedom. These phase-space variables, propagated for time t with the classical Hamiltonian

are denoted $ℏnt$ and $ϕt$ in Eq. (6). Matrix elements of the quantum Hamiltonian are denoted H_ij in Eq. (8), and exact quantum dynamics are recovered for two states by setting for $λ=3/4$ in Eq. (7).⁵³ $A^$ is propagated in time by evolving the classical variables n and ϕ, using those results to determine $A(t;n,ϕ)$ in Eq. (6), and then performing the Fourier series inversion in Eq. (4). Since the dimensionless action n in Eq. (6) corresponds to the mean quantum number $n¯$ in Eq. (5), it takes on the values 0 or 1 for diagonal matrix elements and 1/2 for off-diagonal elements.

Our semiclassical treatment of the vibronic two-level system in the absence of radiation is based on extending the classical Hamiltonian in Eq. (8) without coupling between states to include the nuclear degrees of freedom in Eq. (2),

According to Eq. (6), the electronic dipole operator component in Eq. (3) and its conjugate are represented by the dynamical variables

The nuclear dipole operator is represented by $μnuc=q$⁠, the classical nuclear coordinate subject to the Hamiltonian in Eq. (9). Equations of motion for electronic $(ℏn,ϕ)$ and nuclear (q, p) phase-space variables are

The electronic action is time-independent, while the time-evolution of the angle is determined by the nuclear energy gap U, the difference between adiabatic potentials. The time-dependence of U arises from that of the nuclear coordinate, which evolves on the potential surface V_n, a superposition of adiabatic surfaces V₀ and V₁. Since the electronic action is constant, the time-dependence of the electronic angle is determined by the nuclear trajectory, which is itself independent of the electronic angle.

Integrating these equations gives the classical electronic and nuclear dipoles,

In Eqs. (15) and (17), q(t; n₀, z₀) denotes the nuclear coordinate propagated from initial nuclear phase point z₀ with potential V_n0 in Eq. (14). For the three possible values of the dimensionless electronic action, n₀ = 0, 1, and 1/2, these dynamics take place on the ground-state, excited-state, and mean potential surfaces, respectively. The optimized mean trajectory approximation to nonlinear response functions described in Sec. III is based on these dynamics together with additional quasiclassical quantization of the nuclear action, as in the purely vibrational version of the OMT.^27–31 

A four-wave mixing measurement with three applied optical pulses resonant with electronic transitions is described by the response function components given in Eq. (1) with the electronic dipole operator of Eq. (3). For a given phase-matching condition, the three nested commutators in Eq. (1) give rise to eight terms that are conventionally represented^4,32,33 by double-sided Feynman diagrams (2FD). In the absence of thermal excitation of electrons, two of these terms are nonzero for the rephasing signal R₊₊₋ and for the particular non-rephasing signal R₊₋₊. These nonzero components are represented by 2FD in Figs. 1 and 2, respectively.

Each 2FD is approximately represented by a semiclassical OMT diagram on the right side of each figure, according to the rules previously developed for nonlinear infrared spectroscopy.^27–30 An OMT diagram may be interpreted as a qualitative plot of electronic energy versus time. A horizontal solid line represents a classical trajectory for electronic and nuclear variables at fixed electronic action. A vertical solid line represents a transition in electronic action by $ℏ/2$⁠, with electronic angle and nuclear phase-space variables unchanged. The electronic action for a two-level system may take on the values 0, representing the ground state, $ℏ/2$⁠, representing a coherence, and $ℏ$⁠, representing the excited state. Each filled circle labeled ± represents a factor of $μel±$ in Eq. (16). We treat the case T = 0; generalization to a thermal distribution of nuclear degrees of freedom at nonzero temperature may be carried out as in the OMT for infrared spectroscopy.³¹ 

The uppermost OMT diagram in Fig. 1, representing $ER↑$⁠, has value expressible as a double integral over nuclear and electronic angles

The four phase-space points at which the dipole is evaluated are numbered 0–3 with increasing absolute time and are represented by the four filled circles of the OMT diagram. Each OMT diagram has an overall algebraic sign determined by the commutators in Eq. (1), which can be identified from the sign of the 2FD from which it is derived.^27–30 The diagram is evaluated by propagating classical trajectories in the combined space of electronic and nuclear variables, linked by transitions in electronic action. The nuclear action-angle variables are denoted $z≡(I,θ)$ and the electronic variables are $Z=(ℏn,ϕ)$⁠. The initial values of the phase-space variables at the leftmost filled circle in Fig. 1 are $z0=(ℏ/2,θ0)$ and $Z0=(ℏ/2,ϕ0)$⁠. The initial nuclear action is set to $I0=ℏ/2$⁠, representing a ground-state population, to describe T = 0. The generalization to finite T would include a thermal average over values of I₀ that are restricted to half-odd-integer multiples of $ℏ$⁠.³¹ The nuclear degrees of freedom are propagated on the mean potential surface for time t₁ to produce z₁ = z(t₁; 1/2, z₀). As a result, the electronic phase space variables evolve to $Z1=(ℏ/2,ϕ1)$⁠, with $ϕ1=ϕ0+Δϕ(t1;1/2,z0)$ and with $Δϕ$ defined in Eq. (17). The system then evolves with dimensionless electronic action n = 1, corresponding to nuclear dynamics on the excited-state potential surface, producing the phase-space point z₂ = z(t₂; 1, z₁), $Z2=(ℏ/2,ϕ2)$⁠, with $ϕ2=ϕ1+Δϕ(t2;1,z1)$⁠. Further propagation for t₃ on the mean surface gives the final phase-space point z₃ = z(t₃; 1/2, z₂), $Z3=(ℏ/2,ϕ3)$⁠, with $ϕ3=ϕ2+Δϕ(t3;1/2,z2)$⁠. The four dipole functions in Eq. (18) are then given by

Applying Eqs. (16) and (17) and performing the integral over electronic angle in Eq. (18) give this contribution to the echo signal as an integral over nuclear angle,

Evaluating this quantity requires propagating classical trajectories for nuclear degrees of freedom on the mean potential surface for time t₁, on the excited-state surface for t₂ and again on the mean surface for t₃.

The value of $ER↓$⁠, the remaining rephasing OMT diagram in Fig. 1, may also be written in the form of Eq. (18), but with the nuclear phase-space point z₂ in Eqs. (22) and (23) replaced by z₂ = z(t₂; 0, z₁), signifying that the propagation during the t₂ time takes place on the ground-state rather than the excited-state surface. The non-rephasing OMT diagrams in Fig. 2 may be evaluated by the same rules. $ENR↑$ is calculated from the same trajectories as $ER↑$⁠, but with $±$ superscripts interchanged in Eqs. (19) and (20). The result differs from Eq. (23) only in the sign of the phase associated with evolution during the t₁ period,

The value of the remaining non-rephasing diagram, $ENR↓$ in Fig. 2, has the form in Eq. (24), but with z₂ = z(t₂, 0, z₁) as for $ER↓$⁠. The accuracy of the approximations in Eqs. (23) and (24) is assessed in Sec. V for displaced harmonic oscillators.

In the notation of Ref. 45, VE denotes the vibrational-electronic four-wave-mixing measurement in which the first and second pulses are resonant with vibrational transitions, the third pulse is resonant with an electronic transition, and an optical signal is detected. For a vibronic two-level system at zero temperature, the two nonzero contributions to the rephasing signal $R+⁣+−(3)$ are shown as 2FD labeled $VER↑$ and $VER↓$ in Fig. 3, and their counterparts for the non-rephasing signal $R+⁣−+(3)$⁠, $VENR↑$⁠, and $VENR↓$ are shown in Fig. 4. The initial pair of interactions that induce vibrational transitions are represented by dashed arrows. As in Figs. 1 and 2, radiation-matter interactions inducing electronic transitions are represented by solid arrows. Next to each 2FD is the corresponding OMT diagram. As in the OMT diagrams in Figs. 1 and 2, horizontal solid lines represent continuous classical trajectories of both electronic and nuclear variables. Shorter dashed vertical lines represent transitions in nuclear action by $ℏ/2$⁠, and longer solid vertical lines indicate transitions in electronic action by $ℏ/2$⁠.

To illustrate the rules for evaluating VE diagrams, we treat $VER↑$ in Fig. 3. The leftmost filled circle represents the phase-space point $z0=(ℏ,θ0)$⁠, $Z0=(0,ϕ0)$⁠. The initial nuclear action $I0=ℏ$ represents a vibrational coherence while the dimensionless electronic action n₀ = 0 labels the electronic ground state. From these initial conditions, the nuclear degrees of freedom evolve on the electronic ground state to reach z₁ = z(t₁; 0, z₀) and the electronic variables evolve to $Z1=(0,ϕ1)$ with the electronic angle $ϕ1=ϕ0+Ωt1+Δϕ(t1;0,z0)$⁠, and with $Δϕ$ given in Eq. (17). The nuclear action is then incremented by $ℏ/2$ resulting in the nuclear phase-space point $z1↑=(3ℏ/2,θ1)$⁠, and the system evolves for t₂ to reach $z2=z(t2;0,z1↑)$ and $Z2=(ℏ/2,ϕ2)$⁠, with $ϕ2=ϕ1+Ωt2+Δϕ(t2,0,z1↑)$⁠. The final time evolution for t₃ on the mean adiabatic potential ends in z₃ = z(t₃; 1/2, z₂), $Z3=(ℏ/2,ϕ3)$⁠, with $ϕ3=ϕ2+Ωt3+Δϕ(t3;1/2,z2)$⁠.

The OMT diagram $VER↑$ has the value

with

Substituting these nuclear and electronic dipole functions into Eq. (25) and integrating over electronic angle give the value of $VER↑$ as an integral over nuclear angle,

Evaluation of this integral requires nuclear trajectories run on both the ground-state and mean potential surfaces.

The value of the second contribution to the rephasing response $VER↓$ in Fig. 3 differs by an overall sign and by propagation during the t₂ interval in the vibrational ground state, which replaces $z1↑≡(3ℏ/2,θ1)$ by $z1↓≡(ℏ/2,θ1)$⁠,

with $z2=z(t2;0,z1↓)$⁠.

OMT diagrams representing non-rephasing processes are shown in Fig. 4. $VENR↑$ is evaluated with the same trajectories as $VER↑$ but with superscripts $±$ interchanged on the two factors of the nuclear dipole moment in Eq. (30). $VENR↓$ is similarly evaluated from the same trajectories as $VER↓$ in Eq. (31) but with nuclear dipoles of opposite sense. The accuracy of these approximations for VE spectroscopy is assessed for a displaced harmonic oscillator model in Sec. V.

We assess the accuracy of the OMT approximations for nonlinear electronic spectroscopy in Sec. III and for vibrational-electronic spectroscopy in Sec. IV by comparison to exact results^14,58 for a two-level displaced-harmonic-oscillator (DHO) model with a single nuclear degree of freedom,

The electronic-nuclear coupling for this model is quantified by the Huang-Rhys factor

Calculating electronic spectra in the OMT approximation requires propagating nuclear trajectories on the potential surface V_n(q) in Eq. (14), a superposition of the adiabatic surfaces depending on the dimensionless electronic action n. For the DHO model, the coordinate subject to these dynamics has the form for a forced harmonic oscillator

with n₀ the initial dimensionless electronic action and $(I0,θ0)$ the initial nuclear action and angle. The nuclear trajectory determines the time-dependence of the electronic angle in Eq. (17),

This quantity is required to evaluate the electronic dipole in Eq. (16), necessary for calculating the nonlinear response for both electronic and vibrational-electronic spectroscopies.

To illustrate the OMT calculation of electronic spectra, we evaluate the rephasing diagram $ER↑$ in Fig. 1, with a value given for general potential surfaces in Eq. (23). Evaluating Eq. (23) requires a nuclear trajectory for time t₁ with electronic action n₀ = 1/2, that is, on the mean adiabatic potential, followed by a trajectory for time t₂ with electronic action n₀ = 1, that is, on the excited-state potential surface, followed by a trajectory for time t₃ with n₀ = 1/2. This sequence of trajectories is determined as a function of initial nuclear angle $θ0$ from Eq. (34) and its counterpart for momentum. At T = 0, the initial nuclear action in Eq. (35) is $I0=ℏ/2$⁠. These trajectories determine the phases in Eq. (23) from Eqs. (17) and (36). The average over nuclear angle in Eq. (23) is performed, giving

The Bessel function of the first kind J₀ in Eq. (37) follows from integrating over $θ0$⁠.

The second contribution to the rephasing response, $ER↓$⁠, in Fig. 1 is computed in a similar fashion, but with the evolution during t₂ with electronic action n₀ = 0, signifying the ground electronic state. $ER↓$ has the form in Eq. (37) with $h23→−h23$⁠. The total rephasing response is the sum of these expressions,

Exact results for $ER↑$ and $ER↓$ within the DHO model are derived in  Appendix A and given in Eqs. (A6) and (A7) for nonzero T. At T = 0, the exact results differ from the OMT expressions only by the replacement of the Bessel function with an exponential

These two factors agree to order S.

The non-rephasing terms in Fig. 2, $ENR↑$ and $ENR↓$⁠, are computed, respectively, from the same sets of trajectories as $ER↑$ and $ER↓$⁠. The results are also related to the exact forms in Eqs. (A8) and (A9) by the replacement of a Bessel function with an exponential as in Eq. (42). Their sum gives the non-rephasing response,

The intensity of the rephasing signal is defined as

For purely electronic spectroscopy, we denote this quantity $I+⁣+−E(t1,t2,t3)$⁠. To make a rigorous quantitative assessment of the accuracy of the OMT approach for the DHO model, we compare curves showing one-dimensional time slices $I+⁣+−E(t,0,t)$ on the same plot, rather than showing 2D spectra on separate plots. Comparisons in other time dimensions or in the frequency domain or of the nonrephasing signals do not alter the conclusions. $I+⁣+−E(t,0,t)$ is plotted in Fig. 5 in units of $m08/ℏ6$⁠. Solid lines show OMT results from Eq. (41) and dashed lines show exact results from Eqs. (A6)–(A7). The electronic-nuclear coupling is increased from S = 0.1 (upper frame) to S = 0.5 (middle frame) to S = 1.0 (lower frame). The OMT reproduces the exact result well at low and moderate coupling, but introduces spurious high frequencies as S increases.

The calculations in Fig. 5 illustrate for the DHO model that the extension of the OMT to electronic transitions provides a reasonable approximation to nonlinear spectra. We next perform the corresponding assessment for the OMT expressions for vibrational-electronic (VE) spectroscopy in Eqs. (30) and (31). To illustrate the application of the OMT to VE spectroscopy for the DHO model, we consider the diagram $VER↑$ in Fig. 3 with a value in Eq. (30). This term is evaluated from a trajectory in which the initial nuclear action is $I0=ℏ$ and the initial dimensionless electronic action is n₀ = 0. After propagation for t₁, the nuclear action is incremented to $I0=3ℏ/2$ and the electronic action remains n₀ = 0. After further propagation for t₂, the dimensionless electronic action is incremented to n₀ = 1/2, and the system evolves for t₃. As the nuclear dipole is taken to be proportional to the vibrational coordinate q, the two factors of $μnuc$ in Eq. (30) are evaluated from Eqs. (34) and (35). The electronic phase factor in Eq. (30) is determined from Eq. (36). Performing the average over nuclear angle in Eq. (30) gives

The exact form is given in Eq. (B9) at finite temperature. Setting T = 0 gives a result related to the OMT expression in Eq. (47) by the replacement

As in the corresponding expression for purely electronic spectroscopy in Eq. (42), these two forms agree to order S.

The second contribution to the rephasing signal in Fig. 3 is evaluated similarly but with the nuclear action decreased by $ℏ/2$ at time t₁, resulting in

Similarly to Eq. (47), this result differs from the exact expression in Eq. (B10) by replacing the Bessel function with an exponential agreeing with it to order S. Summing Eqs. (47) and (50) gives the rephasing response for T = 0 in the DHO model,

For the DHO model, the rephasing amplitude is independent of t₂. The rephasing signal vanishes for S = 0, so that the VE echo in the DHO model may be viewed as induced by electronic-nuclear coupling. Cho⁴⁰ has derived rules constraining this and related measurements. For the VE echo, to first order in linear electronic-nuclear coupling and for nuclear-radiation interaction operators linear in coordinates, potential surfaces must be anharmonic for the signal to be nonzero. That implies that R₊₊₋ in Eq. (51) must vanish to first order in d₀. Indeed, this expression is linear in S and thus quadratic in d₀ to lowest order, so that the OMT approximation preserves this property of the exact solution in Eqs. (B9) and (B10).

The non-rephasing contributions to the VE signal in Fig. 4, $VENR↑$ and $VENR↓$⁠, are computed from the same trajectories just described for $VER↑$ and $VER↓$ and satisfy $VENR↑(t1,t2,t3)=VER↑(−t1,t2,t3)$ and $VENR↓(t1,t2,t3)=VER↓(−t1,t2,t3)$⁠, so that

The intensity of the rephasing signal $I+⁣+−VE(t1,t2,t3)$ is given by these results with Eq. (46). The OMT approximation to $I+⁣+−VE(t,0,t)$ in units of $(m02/mωℏ2)2$ is assessed in Fig. 6. Dashed lines show exact results from Eqs. (B9)–(B10) and solid lines show the OMT approximation from Eq. (51). The electronic-nuclear coupling is increased from S = 0.2 (upper frame) to S = 0.5 (middle frame) to S = 1.0 (lower frame). As for the electronic response shown in Fig. 5, the OMT performs well at low or moderate coupling but becomes less accurate as the coupling is increased.

The OMT method, as originally formulated for multidimensional infrared spectroscopy,^27–30 allows the calculation of vibrational response functions from classical trajectories with quantized action variables. This procedure is broadly applicable to a variety of measurements as it is based on the mapping of exact quantum mechanical density operator diagrams onto semiclassical OMT diagrams. By describing electronic degrees of freedom in terms of classical phase-space variables, we have extended the OMT protocol to describe any nonlinear measurement involving sequences of infrared and optical excitation pulses. The resulting method requires quantizing nuclear action variables, as for vibrational response functions, and running sequences of classical trajectories on ground or excited state adiabatic surfaces or on mean surfaces.

Computing nonlinear electronic spectra from energy gaps that evolve through classical trajectories is an established procedure,^59,60 representing a mechanical realization of stochastic line shape theory.^32,61 Fried and Mukamel^62,63 developed semiclassical treatments of nonlinear electronic spectroscopy requiring propagation of energy gaps on ground or excited state surfaces. Shemetulskis and Loring⁶⁴ expressed the electronic photon echo amplitude for two levels in terms of Wigner transforms and resummed an expansion in ℏ to obtain a procedure in which electronic coherences are represented by propagation on the mean potential surface. Mixed quantum-classical Ehrenfest calculations of multidimensional electronic spectra⁶⁵ similarly involve propagating nuclear degrees of freedom on adiabatic surfaces.

The relationship between the OMT procedure and linearized approximations involving semiclassical propagators in the initial value representation, sometimes denoted LSC-IVR^14,66–68 merits further discussion. Both methods have been applied to vibrational spectra on a single adiabatic surface and to electronic spectra on multiple surfaces, and these distinct applications must be considered separately. For purely vibrational spectra, the central approximation of the LSC-IVR method is a simplification of the total phase associated with semiclassical representations of pairs of quantum propagators of opposite time sense. This phase depends on a pair of classical trajectories with different initial conditions. The phase is linearized in differences in initial phase-space variables, as would be correct for propagation on a harmonic potential. To this order the phase vanishes, resulting in an approximation that neglects quantum recurrences in nonlinear response functions.^23,68 By contrast, the semiclassical quantization condition of the OMT may be derived by continuing the expansion of this phase to higher order,²³ thereby incorporating the anharmonic effects that lead to quantum recurrences in response functions.⁶⁸ In calculating vibrational spectra on a single adiabatic surface, the OMT is distinct from linearized semiclassical methods. The application of the LSC-IVR method to electronic spectra involving multiple adiabatic surfaces invokes the linearization of trajectories with respect to initial conditions and also a linearization with respect to the change in the potential surface. Rather than vanishing, the semiclassical phase is proportional to a time integral of the energy gap. Dynamics of electronic populations are represented by energy gap propagation on the appropriate adiabatic surface, and electronic coherence dynamics are represented by propagation on a mean surface.¹⁴ The application of the OMT to electronic spectroscopy of a pair of adiabatic surfaces in Sec. II generates the same dynamical approximations. However, the applications of the OMT to a purely vibrational measurement on a single surface or to a measurement with both vibrational and electronic excitation as treated in Sec. IV go beyond a linearized theory.

McRobbie and Geva have shown that the LSC-IVR treatment of linear and nonlinear electronic spectra for the two-level displaced harmonic oscillator model of Sec. V is exactly correct, in contrast to the OMT treatment of the same model as seen in Fig. 5. The difference in results lies in the sampling of initial conditions, as we illustrate with reference to the comparison in Eq. (42) between OMT and exact results for the electronic rephasing amplitude at zero temperature. The Bessel function in Eq. (42) results from carrying out the integration over initial nuclear angle in Eq. (23) and setting the initial nuclear action to $I0=ℏ/2$ according to OMT quantization rules. It is useful to express this quantity as an average over nuclear action with generalized weight F(I₀),

The OMT result is obtained by setting $F(I0)→fOMT(I0)$⁠,

The LSC-IVR result, exactly correct for this model, is obtained by setting the weight F(I₀) in Eq. (53) to the zero-temperature Wigner distribution for the harmonic ground state, expressed in terms of action,

It is noted in Sec. V that the right-hand sides of Eqs. (55) and (57) agree to order S and hence that response functions calculated with the OMT are correct to second order in electron-nuclear coupling. Expanding the integrand in Eq. (53) in powers of S shows that this agreement originates from the two distributions f_OMT(I₀) and $ρW(I0)$ having the same first moment.

The thermal weight for OMT response functions³¹ was not derived as an approximation to the Wigner distribution but has the form of an approximate Wigner distribution in action-angle variables discussed by Miller and Cotton.⁵⁷ These authors draw the distinction between the continuous Wigner distribution obtained exactly in Cartesian coordinates and momenta and then transformed to action-angle variables, and a discretized approximation to the Wigner distribution in action-angle variables derived from semiclassical wavefunctions. The latter form embodies semiclassical quantization rules, and its advantages in other contexts are discussed in Ref. 57.

Replacement of the OMT distribution, e.g., Eq. (54), by the exact Wigner function, e.g., Eq. (56), in OMT expressions for electronic spectroscopy of multiple surfaces yields an exact result for the DHO model of Sec. V. The approximate Bessel function in Eq. (55) is replaced by the exactly correct exponential in Eq. (57). The OMT expressions for VE response functions in this model also differ from exact results by the replacement of an exponential function by a Bessel function as in Eq. (49). However, it must be emphasized that replacing the OMT distribution by the Wigner distribution in this VE response function does not produce a correct result for the DHO model. The reason is that the effects of infrared excitation are treated differently in the OMT than in a linearized theory. In summary, the OMT does not exactly reproduce the observable in a purely electronic measurement for the DHO model but goes beyond a linearized treatment for purely vibrational and mixed vibrational-electronic experiments. The extended OMT presented here generalizes previous treatments of nonlinear spectroscopy by providing a systematic way of treating any sequence of infrared or optical pulses and unambiguously prescribing the classical trajectories to be employed for each time interval.

As for the OMT applied to vibrational spectroscopy, the extended OMT requires approximately identifying nuclear action and angle variables, in order that a quantization condition can be applied to the action variables. In principle, this raises concern over the applicability of the method to large anharmonic systems, for which good action and angle variables generally cannot be defined. In practice,^29,30 we have found in computing two-dimensional infrared spectra of anharmonic systems with the OMT, that on the short time scales relevant to ultrafast spectroscopy we require only that the system behaves approximately as though good action and angle variables exist. In the applications we have considered, low order classical mechanical perturbation theory in anharmonicity⁶⁹ suffices to define the canonical transformations between Cartesian coordinates and approximate action-angle variables.

In the model treated here, the only nonadiabatic interaction between electronic states is that driven by radiation, and so is treated perturbatively in a response function calculation. The effective classical Hamiltonian in Eq. (9) describing propagation in the absence of radiation contains no further coupling between states. It remains to be determined how the method treats models with multiple electronic excited states, subject to nonadiabatic interactions.^70–72 Also to be performed is a quantitative assessment of electronic and vibrational-electronic spectroscopy with anharmonic nuclear surfaces, or, indeed harmonic surfaces with different frequencies for which linearized theories are no longer exact.¹⁴ The successes of the OMT in describing anharmonic effects in vibrational spectroscopy suggest that the method will yield qualitatively reasonable results for nonlinear electronic and vibrational-electronic spectra of large anharmonic systems.

This material is based upon the work supported by the National Science Foundation under CHE-1361484.

The rephasing and non-rephasing amplitudes represented by 2FD in Figs. 1 and 2, respectively, are exactly written in terms of quantum equilibrium dipole correlation functions as

Angle brackets indicate the thermal average over both nuclear and electronic degrees of freedom. Tracing over electronic degrees of freedom yields the nuclear operators $⟨0|μ^el−(t)|1⟩=(⟨1|μ^el+(t)|0⟩)†$⁠. For the DHO model,⁵⁸ 

Here $b^†$ is the boson creation operator. The equality in Eq. (A5) follows from the exact truncation⁵⁸ of the Baker-Campbell-Hausdorff expansion for operators linear in $b^$ and $b^†$⁠. Repeated application of this identity together with the exact truncation⁵⁸ of a cumulant expansion of the nuclear average at second order gives closed form expressions for the response function components in Eqs. (A1)–(A4),

The functions g_R, g_NR, h₂₃, f_R, and f_NR are defined in Eqs. (38)–(40), (44), and (45). Summing the rephasing and non-rephasing terms, respectively, yields R₊₊₋ and R₊₋₊ for electronic spectroscopy.

The four 2FD in Figs. 3 and 4 have values

As in Eqs. (A1)–(A4), angle brackets denote averaging over electronic and nuclear degrees of freedom. Performing the electronic trace for the DHO model yields the nuclear operator in Eq. (A5) as well as

Evaluating the thermal averages for $VER↑$ and $VER↓$⁠, respectively, requires the identities