Simulation of Multi-Dimensional Signals in the Optical Domain

Article pubs.acs.org/JCTC

Simulation of Multi-Dimensional Signals in the Optical Domain: Quantum-Classical Feedback in Nonlinear Exciton Propagation Martin Richter and Benjamin P. Fingerhut* Max-Born-Institut für Nichtlineare Optik und Kurzzeitspektroskopie, D-12489 Berlin, Germany S Supporting Information *

ABSTRACT: We present an algorithm for the simulation of nonlinear 2D spectra of molecular systems in the UV−vis spectral region from atomistic molecular dynamics trajectories subject to nonadiabatic relaxation. We combine the nonlinear exciton propagation (NEP) protocol, that relies on a quasiparticle approach with the surface hopping methodology to account for quantum-classical feedback during the dynamics. Phenomena, such as dynamic Stokes shift due to nuclear relaxation, spectral diffusion, and population transfer among electronic states, are thus naturally included and benchmarked on a model of two electronic states coupled to a harmonic coordinate and a classical heatbath. The capabilities of the algorithm are further demonstrated for the bichromophore diphenylmethane that is described in a fully microscopic fashion including all 69 classical nuclear degrees of freedom. We demonstrate that simulated 2D signals are especially sensitive to the applied theoretical approximations (i.e., choice of active space in the CASSCF method) where population dynamics appears comparable.

1. INTRODUCTION The advent of ultrashort laser technology recently allowed the extension of multidimensional spectroscopic techniques to the visible1,2 and ultraviolet3−8 spectral region where organic molecules, such as, DNA bases and aromatic compounds, possess prominent electronic absorption bands. Especially upon UV excitation the short time dynamics of these systems is governed by ultrafast excited state deactivation channels that are mediated by the strong coupling of intramolecular motion with electronic degrees of freedom, that is, nonadiabatic relaxation via conical intersections (CI). As such diphenylmethane (DPM) can be regarded as a prototypical bichromophore in the UV spectral region where two distinct aromatic phenyl moieties are coupled by a bridging CH2 group with the potential for ultrafast energy transfer dynamics among excited states. Early experimental investigations reported a single broad peak in the UV−vis spectrum in gas phase and solution9 which was assigned to the vibrationally coupled and nearly degenerate S1 and S2 electronic singlet excited states.10−13 Coherent 2D spectroscopy allows to resolve couplings between different chromophores by spreading the spectral information content in 2D14 with the ability to resolve vibronic structure15 on the femtosecond time scale. The technique thus has the potential to selectively map energy flow between chromophores in realtime and to resolve congested spectra, making DPM a promising candidate to obtain microscopic insight into ultrafast relaxation dynamics of electronically excited molecules. Ab initio molecular dynamics (MD) allows to simulate the aforementioned relaxation phenomena in full coordinate space by employing high-level quantum mechanical methods that © XXXX American Chemical Society

describe excited electronic state properties on-the-fly during the dynamics simulations while treating nuclear degrees of freedom by means of classical mechanics. These simulations are commonly interpreted on the basis of state populations to compare to time-resolved experimental data.16−19 Nevertheless the direct simulation of the nonlinear spectroscopic signals allows for a thorough interpretation of spectral features20−27 and could provide a rigorous test of the employed force fields and approximations.28,29 A common approach for the simulation of 2D spectroscopic signals relies on the separation of the total system in a quantum and classical subsystem where the classical subsystem acts as bath that defines time-dependent phenomena of a preparametrized Hamiltonian of the quantum subsystem.30 In the case of Gaussian fluctuations the cumulant expansion provides an efficient route for the description of transport dynamics31,32 beyond the weak system-bath coupling limit employed in Redfield theory. Numerically exact methods like the hierarchical equations of motion (HEOM)33,34 allow to recover the coupling of populations and coherences but the extension to complex spectral densities including underdamped vibrations, as common for molecular systems is challenging.35 In the case of large structural rearrangements apparent in light-induced chemical reactions36 and nonadiabatic relaxation via CI the reduced density matrix approach focusing on the electronic degrees of freedom only reaches the limits of its applicability. An alternative approach that accounts for the microscopic origin of (non-Gaussian) fluctuations is commonly employed Received: April 12, 2016

A

DOI: 10.1021/acs.jctc.6b00371 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

Figure 1. Left: Effective NEP ladder diagram arising from the contraction of ESE, GSB, and ESA contributions to the kI signal. Middle: Contributions of populations and coherences to the effective NEP diagram. Right: Schematic of the model consisting of two electronic states coupled to classical harmonic oscillator and a stochastic heatbath.

methodology we show that back-coupling between the classicaland quantum subsystem can be efficiently accounted for allowing for the investigation of time-dependent dynamics in excited electronic states (e.g., population branching and dynamic Stokes shift). The capabilities of the method are demonstrated on a system of two coupled electronic states in turn coupled to a harmonic nuclear coordinate and phenomena like dynamic Stokes shift, spectral diffusion and population transfer are compared to accurate NISE and HEOM results.41,42 Further, we demonstrate that the molecular bichromophore DPM can be described in a fully microscopic fashion by abinitio MD to simulate 2D spectra upon electronic excitation. The microscopic description via trajectory dynamics allows for the investigation of energy fluctuations due to nuclear dynamics and population transfer without beforehand knowledge of the molecular PES. Linear absorption and nonlinear 2D spectra are presented revealing ultrafast population dynamics. Special care is given to the convergence of the active space in the employed CASSCF method. Our simulations demonstrate and reiterate that the choice of the active space can have crucial effects on the simulated spectra while population dynamics of both approximations is similar. The remainder of this work is organized as follows: Section 2 recaptures the main aspects of the NEP protocol and TSH molecular dynamics and describes the combination of both methods in the QF-NEP protocol. Section 3 presents benchmark results on a model system of two coupled electronic states, each coupled to a harmonic classical bath coordinate. Section 4 presents results of the QF-NEP protocol applied to DPM where the nonadiabatic population dynamics is evaluated via TSH. Section 5 gives a discussion and conclusions.

for the simulation of 2D spectroscopic signals in the IR spectral region and relies on the direct propagation of a time-dependent Hamiltonian along the configurations imposed by a classical trajectory.37−39 This exciton propagation approach allows to treat arbitrary lineshapes in 2D spectroscopic signals by alleviating the Gaussian fluctuation criterium while its applicability to electronic excited states is limited. The limitations arise from the fact that backcoupling between the quantum system and the classical system is absent and thus nuclear relaxation effects like the dynamic Stokes shift are not accounted for. Moreover the exciton propagation approach does not guarantee thermal equilibrium for the branching of populations, that is, the detailed balance limit between reaction channels. Recently Jansen and co-workers reported the combination of Tully’s fewest switches surface hopping (TSH) algorithm40 in the framework of a sum-over-states (SOS) Ansatz for the simulation of third-order nonlinear 2D spectra by numerical integration of the Schrödinger equation (NISE).41,42 Their methodology relies on a pathway decomposition scheme that separates between population and coherence terms and allowed to recover the detailed balance limit and dynamic Stokes shift for a generic coupled two-site model. Applications employ preparametrized transition densities for 2D spectra simulations of molecular aggregates43 and the calculation of time-resolved fluorescence spectra44 in J- and H-aggregates. Within this work, we report a quasi-particle Ansatz that accounts for the feedback between the quantum and classical subsystems in the calculation of nonlinear 2D spectra thereby alleviating some limitations of a SOS Ansatz.30,45,46 SOS methods describe the system in terms of quantum mechanical eigenstates requiring diagonalization of the biexciton Hamiltonian that makes the SOS Ansatz prohibitively expensive for large systems where the number of biexcitons increases with N2. In a quasi-particle description, one switches to a basis of coupled electronic oscillators where two-exciton resonances are introduced by exciton−exciton scattering47 thus avoiding the diagonalization of the biexciton Hamiltonian. Arbitrary fluctuations can be treated by a direct nonlinear exciton propagation (NEP) protocol39 where all nonlinearities of the 2D signal are accounted for within a single diagram. The quasiparticle description offers thus a computationally efficient algorithm for the simulation of nonlinear 2D spectra that exploits the cancellations of interfering terms in the SOS Ansatz.47 By combining the NEP protocol with the TSH

2. THEORY We start by briefly recapturing the essence of the NEP algorithm (for a detailed derivation see ref 39), followed by introducing the quantum-classical backcoupling that allows to describe nonadiabatic excited state relaxation. 2.1. Nonlinear Exciton Propagation (NEP) Algorithm. The NEP algorithm contracts contributions of distinct diagrams of the SOS Ansatz (i.e., excited state emission (ESE), ground state bleach (GSB), excited state absorption (ESA)) into one effective diagram. This single diagram contains the observable nonlinearities which arise in a specific four-wave mixing setup experiment (cf., kI = −k1 + k2 + k3 signal, Figure 1). B


Article


has been created at time τ1 and propagated until τ2 where the biexciton has been created as symmetrized product

The calculation of the kI signal can then be expressed in the compact form Sνk4Iν3ν2ν1(τ4 , τ3 , τ2 , τ1) ⎛i⎞ = 2⎜ ⎟ ⎝ℏ⎠ ×

τ4

∫τ3

ν2 (1) ν2 (1) Ψ(2) mn(τ2 ; τ2 ; τ1) = μn (τ2)Ψ m (τ2 ; τ1) + μm (τ2)Ψ n (τ2 ; τ1)

4

∑ ∑

(7)

μnν1(τ1)μnν2 (τ2)μnν3(τ3)μnν4 (τ4) 1

n1 − 4 m1 − 5

2

3

4

In eqs 5−7 polarization dependence of the wave functions (Ψ(2);ν1,ν2, Ψ(1);ν1) has been dropped for clarity. Taking into account anharmonicities due to exciton−exciton scattering, the response function can be expressed with help of the auxiliary vectors XkmI4 and Rkn4I

dsGn4m5(τ4 , s)Um5m1m4m3(s)Gm4m3n3m2(s , τ3)Gm2n2(τ3 , τ2)Gm*1n1(s , τ1)

(1)

where τi denote the interaction times of the system with light pulses, μnνii(τi) are the time dependent transition dipole moments of exciton ni generated by the pulse with polarization νi. The matter response can be directly read of the diagram (Figure 1, kI-NEP) where solid arrows represent the forward (backward) propagation of single-exciton Green’s functions and double arrows describe the free evolution of biexciton Green’s function: the interaction with laser pulse 2 creates a singleexciton at time τ2 that propagates forward from τ2 to τ3 described by the single-exciton Green’s function Gm2n2(τ3, τ2). At τ3 the interaction with the laser pulse generates a biexciton Gm4m3n3m2(s, τ3). Nonlinearities are accounted for by the time integral over s between τ3 and τ4 due to the scattering interaction Um5m1m4m3(s). Exciton scattering splits the biexciton into two single-excitons, one propagating forward from s until detection time τ4, described by the single-exciton Green’s function Gn4m5(τ4,s) and the other propagating backward in time from s to interaction time of the first laser pulse τ1 (singleexciton retarded Green’s function G*m1n1(s, τ1)). Single- and biexciton Green’s function are defined by Gmn(τ2 , τ1) = ⟨g |BmU (τ2 , τ1)Bn†|g ⟩

X mkI4(s ; τ3 , τ2 , τ1) =

(8)

R nk4I(τ4 , τ3 , τ2 , τ1) =

1

τ2

⎛ i ⎞4 = 2⎜ ⎟ ∑ μnν4 (τ4) × R nk4I(τ4 , τ3 , τ2 , τ1) 4 ⎝ℏ⎠ n4

(10)

The integral in eq 9 is numerically integrated using trapezoid rule and the calculation of the Green’s function Gn4,m4(τ4,s) (eq 5) is again avoided by propagating RkmI(τ4 + Δτ, τ3, τ2, τ1) in a wave function like fashion. 2.2. Quantum-Classical Feedback in the Nonlinear Exciton Propagation Algorithm. We consider nonadiabatic relaxation in the framework of mixed quantum classical dynamics40 where the complete evolving system is partitioned into the quantum system of electronic degrees of freedom and classical nuclear degrees of freedom following a classical trajectory R̅ (t). The time-dependence of the Hamiltonian H(t) arises from the instantaneous solution of the stationary electronic Schrödinger equation along R̅ (t) thus accounting for fluctuations of the instantaneous energy levels. The resulting response function include the dynamic characteristics of the underlying molecular system at the employed ab-initio level of theory defining lineshapes during t2 and t3. Energy fluctuations during t2 allow to access time-dependent spectral features such as spectral diffusion and dynamic Stokes shift. The time evolution of the single-exciton wave function (eq 5) is governed by the time-dependent Schrödinger equation

(4)

and exp+ is the time ordered exponential. and Bm denote boson creation and annihilation operators of the m-th excitation ([Bm, B†n] = δmn). The respective expression for the kII = k1 − k2 + k3 technique is given in the Supporting Information (SI). Instead of propagating the single- and biexciton Green’s function directly, single- and biexciton wave functions are employed for numerical simulation of the nonlinear response. This representation allows to incorporate nonadiabatic transitions between multiple electronic states via CIs (see below). The single- and biexciton wave functions are defined as

∑ Gmn(τ2 , τ1)μnν (τ1)

dsGn4 , m4(τ4 , s) × X mk4I (s ; τ3 , τ2 , τ1)

Sνk4νI 3ν2ν1(τ4 , τ3 , τ2 , τ1)

B†m

Ψ(1) m (τ2 ; τ1) =

τ3

With these definitions, the final response function of the kItechnique reads

(3)

⎞ H(τ )dτ ⎟ ⎠

τ4

(9)

respectively, where U(τ2, τ1) denotes the time evolution operator

∫τ

∑∫ m4

and

⎛ i U (τ2 , τ1) = exp+⎜ − ⎝ ℏ

(1) * Um4m1m3m2(s) × Ψ(2) m3m2(s ; τ3; τ2)Ψ m1 (s ; τ1)

m1m2m3

(2)

Gmn , kl(τ2 , τ1) = ⟨g |BmBnU (τ2 , τ1)Bk†Bl†|g ⟩

∑

iℏ

∂ (1); ν |Ψ (t ; t0)⟩ = H(t )|Ψ(1); ν(t ; t0)⟩ ∂t

(11)

By expanding the single-exciton wave function in a set of instantaneous (adiabatic) eigenfunctions ϕ̃ k(t) that depend only parametrically on time

1

|Ψ(1); ν(t ; t0)⟩ =

(5)

n

k

and Ψ(2) mn(τ3; τ2 ; τ1) =

ν (t ; t0)B†|g ⟩ = ∑ ck(t )|ϕk̃ (t )⟩ ∑ Ψ(1); k k

(12)

∑ Gmn,kl(τ3 , τ2) × μkν (τ2)Ψ(1) l (τ2 ; τ1)

the time evolution of expansion coefficients can be propagated in infinitesimal small steps along a time dependent Hamiltonian:

2

kl

(6) (1);ν

Ψm 1(τ2; τ1) describes the free evolution of a single exciton (2);ν ,ν generated at time τ1 until τ2 and Ψmn 1 2(τ3; τ2; τ1) describes the time evolution of biexcitons during τ3 whose single exciton

∂ck(t ) i = − ∑ Hkl(t )cl(t ) ∂t ℏ l C

(13) DOI: 10.1021/acs.jctc.6b00371 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article


4. Perform subpathway decomposition according to active states (m = n are interpreted as populations and m ≠ n as interstate coherences) (cf., Figure 1). Classical bathmodes of population states (m = n) experience the quantum feedback force FQ, resulting in a dynamic Stokes shift of the exciton energies during t2. Interstate coherences (m ≠ n) are propagated during t2 with the ground state reference Hamiltonian, neglecting quantum feedback. 5. At τ1 + t1 + t2, force the one-exciton wave functions to collapse to current active state and generate a biexciton wave function (eq 7) that is propagated until τ1 + t1 + t2 + t3. 6. Propagate auxiliary vector RkmI during t3 used for the evaluation of the response function S(t3, t2, t1) (eq 10) and weight subpathway contributions to response function by corresponding coefficient Pmn. The procedure is repeated for a swarm of trajectories with different initial conditions where for each initial condition individual trajectories of n excited states and the ground states are propagated (population and interstate coherence contributions). The stochastic surface hopping algorithm accounts for ensemble averaging and correct population branching of an approximate wave packet and enters the 2D spectra simulation during the propagation of the two one-exction wave functions during the t2 delay time. Accordingly populations can be transformed into interstate coherences and vice versa. The reference Hamiltonian in the outlined procedure can be arbitrarily chosen during t1 and t3 as well as for interstate coherence pathways during t2.30 We follow the suggestion of Tempelaar et al.42 where the ground state Hamiltonian is used as reference for all appearing coherence terms. We note that different approaches exist49 (see section 5).

with Hkl(t ) = ϵk (t ) + Kkl(t )

(14)

Here, ϵk denotes the time dependent eigenvalue matrix that is diagonal in the adiabatic representation and Kkl denotes timederivative couplings responsible for nonadiabatic relaxation between electronic states: Kkl(t ) = ⟨ϕk̃ (t )|∂/∂t |ϕl̃ (t )⟩ ∂R ̃ ⟨ϕ (R(t ))|∂/∂R |ϕl̃ (R(t ))⟩ ∂t k = Ṙ(t )Dkl (R(t ))

=

(15)

with the matrix of nonadiabatic coupling vectors Dkl. During the exciton propagation, the quantum system is affected by the bathmodes (nuclear degrees of freedom), inducing spectral diffusion and energy transport due to nonadiabatic couplings between electronic states of the time dependent Hamiltonian. Nevertheless self-consistent back coupling of quantum and classical subsystems is neglected and fluctuations are determined by the evolution on a single potential energy surface.30,46 The motion of classical bath modes is then dictated by the quantum mechanically evaluated force of a single quantum state a using the Hellmann− Feynman theorem F Q = −∇⟨ϕa|H |ϕa⟩ = −⟨ϕa|∇H |ϕa⟩

(16) 40

Introducing the idea of surface hopping allows to lift this limitation by intuitively following the motion of classical bath modes on a single potential energy surface but allowing for hops between different electronic states induced by nonadiabatic coupling. The probability for a hop from the currently active state a to another state i within the time interval Δt is40,48 pa → i =

( (

2ℜ ca*ci

i H ℏ ai

ca*ca

+ Kai

))

3. MODEL DIMER WITH A SINGLE HARMONIC NUCLEAR COORDINATE To demonstrate the capabilities of the QF-NEP algorithm, we investigate a model system composed of two single-exciton states each coupled to a classical nuclear coordinate. To allow for energy dissipation and equilibration, the classical coordinates are coupled to a stochastic heat bath. Model parameters are adopted from refs 41 and 42 to allow comparison to reported NISE and HEOM results. 2D signals are presented as absorptive 2D spectra SA calculated from the respective kI and kII signals14,46 (simulation details are given in the SI). We start by discussing the linear absorption spectrum of the model (Figure 2a). Two resonances can be identified corresponding to the eigenstates of the single-exciton Hamiltonian (adiabatic state energies of ±269 cm−1, cf. eq SI.2). Broadening of resonances arises due to the trajectory inhomogeneity. As we assumed parallel dipole moments in the local site representation the intensity is redistributed toward the high-energy exciton state (corresponding to a H-aggregate). Additionally, we investigated the convergence of the linear absorption spectrum with the number of employed trajectories (N = 256, 4096, 60000). We find that for τl = 50 fs the global line shape of the linear absorption spectrum is reproduced qualitatively already with 256 trajectories. Figure 2b presents the branching of population between reaction channels due to the surface hopping algorithm. We observe the transfer from higher to lower exciton state where thermal equilibrium is approached after ∼5 ps (0.92 for an energy splitting of 538 cm−1, cf. gray line). For ensembles of N = 100−200 trajectories the short time equilibration dynamics can be captured only

Δt (17)

The state amplitudes of the quantum subsystem are propagated in small steps along a Hamiltonian (eq 13), that parametrically depends on the configuration of the classical bath modes. During the propagation, a stochastic algorithm is used to decide whether instantaneous switches between the active state a and inactive states i occur. The signal simulation within the nonlinear exciton propagation algorithm with quantum-classical feedback (QFNEP) is evaluated with the response function SkI(t3, t2, t1) in time intervals t1 = τ2 − τ1, t2 = τ3 − τ2 and t3 = τ4 − τ3 and involves the following steps: 1. Initialize the single exciton wave function |Ψ(1);ν1(t1 + τ1; τ1)⟩ at time τ1 and propagate until τ1 + t1 (eq 12). 2. At τ1 + t1, initialize a second single exciton wave function |Ψ(1) ;ν2(τ1 + t1 + t2; τ1 + t1)⟩. 3. At t2 = 0, evaluate single exciton correlation coefficient matrix Pmn Pmn = c1*m(τ1)c1m(τ2)c 2*n(τ2)c 2n(τ2)

(18)

with c1m being the quantum amplitude of exciton m of Ψ generated at τ1 and c2n of exciton n of Ψ(1) ;ν2 generated at time τ2, respectively. Afterward, reinitialize the one-exciton wave functions to their respective active state.

(1) ;ν1

D


Article


shoulder above the diagonal. Contribution from biexciton states to the absorptive 2D spectra appear anharmonically shifted by 2000 cm−1 along the Ω3 axis and represent the mirror image of single-exciton contributions to the 2D signal. At t2 = 300 fs (Figure 2d), spectral diffusion induces a change in peak shapes from ellipsoid to nearly round where spectral diffusion dynamics reflects the employed 220 fs correlation time of the model. Additionally, a shift of diagonal peaks toward lower detection frequencies Ω3 can be idenitifed, a clear signature of the dynamic Stokes shift induced by the nuclear relaxation that is recaptured by the QF-NEP algorithm. The onset of population transfer from the higher to the lower exciton state induces a pronounced crosspeak below the diagonal and to a lesser extend above the diagonal (t2 = 1.5 ps, Figure 2e). For t2 = 15 ps, we observe thermally equilibrated spectra, where, because of population transfer, a cross-peak of higher intensity than diagonal peaks is developed below the diagonal, in agreement with the findings of ref 42. This crosspeak is further shifted from the Ω1 excitation frequency due to nuclear relaxation (dynamic Stoke’s shift). The spectral signatures arising due to biexciton contributions similarly indicate the population transfer from the higher to the lower adiabatic state. Population transfer between excitons is quantified by integrating the population transfer induced peak intensities of single-excitons (0 cm−1 < Ω1 < 500 cm−1 and −750 cm−1 < Ω3 < 0 cm−1) and biexcitons (0 cm−1 < Ω1 < 500 and 1250 cm−1 < Ω3 < 2000 cm−1) (symbols in Figure 2b). The integrated peak intensities allow one to follow nicely the underlying population transfer of the trajectory ensemble from the higher adiabatic state to the low energy state; the thermal equilibrium is approximately recovered. For comparison, we present the calculated spectrum neglecting back-coupling between the quantum and classical subsystem (inlay Figure 2f, t2 = 15 ps). In the resulting spectrum, only a weak off-diagonal peak is developed below the diagonal demonstrating that population transfer does not fulfill thermal equilibrium for long delay times t2. Moreover the diagonal peaks appear located at Ω1 = Ω3 due to the missing Stokes shift, while the peak shapes are modulated by spectral diffusion.

4. AB-INITIO-BASED SIGNAL SIMULATION ON DIPHENYLMETHANE We further demonstrate the simulation of 2D signals of nonadiabatic relaxation of the bichromophore DPM directly from ab-initio on-the-fly surface hopping MD simulations on the basis of the QF-NEP protocol. DPM (see inlay of Figure 3) is composed of two aromatic phenyl moieties coupled by a bridging CH2 group and is therefore an excellent candidate molecule for investigating intramolecular population relaxation dynamics. All quantum mechanical calculations were performed with the MOLPRO 2012 program package,50 surface hopping MD was carried out using the SHARC software.48,51,52 (Computational details are given in the SI.) 4.1. Excited States and Adiabatic State Populations. The vertical excitation energies of DPM at the ground state minimum (C2 symmetry) are listed in Table 1. In resonant twophoton ionization and single vibronic level fluorescence experiments, the first two excited states of DPM are reported at 4.63−4.64 eV. The S1/S2 splitting at the CASSCF(8,8) level of theory that includes the π system of both phenyl rings but excludes the complete bonding and antibonding orbitals (for details see SI), is in good agreement with the experimentally

Figure 2. (a) Linear absorption of a model system of two coupled quantum states for 256, 4096, and 60000 trajectories. (b) Population transfer for 10−4096 trajectories. Thin gray line: Boltzmann limit. Green circles: Integrated peak intensities of ESE signal (0 < Ω1 < 500 and −750 < Ω3 < 0 cm−1). Blue squares: Integrated peak intensities of ESA signal (0 < Ω1 < 500 and 1250 < Ω3 < 2000 cm−1). ESE and ESA (c−f) calculated 2D spectra for different delay times t2 = 0−15 ps (256 trajectories); the inset in panel f depicts the spectrum without quantum feedback.

qualitatively (2−4 ps), while the long time equilibrium is reproduced within a few percent error. Accordingly, the surface hopping scheme in the adiabatic representation allows to reproduce the Boltzmann limit for population branching, a key requirement for accurate spectra predictions. Absorptive 2D spectra calculated for ensembles of N = 256 trajectories at different delay times t2 are presented in Figure 2c−f. At t2 = 0 ps, the two peaks arising from exciton eigenstates can be identified along the diagonal of the spectrum. We note that the large inhomogeneity of the model (cf., Ftherm, eq SI.4) leads to ellipsoid linehapes. The asymmetry is introduced by excitonic coupling J that shows up as weak E


Article


al.55 demonstrated that the Born−Oppenheimer approximation employed in static ab-initio calculations systematically overestimates the S1/S2 energy gap of nearly degenerate, weakly coupled dimers, whereas a vibronic treatment reduces the S1/S2 energy splittings. In that sense, DPM can be considered as a phenyl dimer, weakly coupled by the bridging CH2 group. Additionally, we find that the S1/S2 splitting is extremely sensitive to the employed molecular structure. Slight distortions out of symmetry increase the S1/S2 energy gap substantially. Therefore, the lowest energy gap was obtained when calculating the energies at the level of theory used during geometry optimiziation, that is, CASSCF(8,8). The averaged adiabatic state populations of three lowest electronic states of DPM are shown in Figure 3 where we assume that electronic excitation is initiated by a δ-pulse populating the S1 and S2 state (for fitted lifetimes and fitting procedure see SI, Table SI.2). At the CASSCF(8,8) level of theory (Figure 3, solid line), population decay from S2 to S1 proceeds immediately, characterized by a 51 fs time scale. Population of S0 starts delayed (208 fs) which can be interpreted as the required time to reach nuclear configurations that mediate S1 → S0 decay and proceeds on a picosecond time scale. Population transfer happens in the vicinity of CoIns that are characterized in the SI (section SI3.4). For S2 → S1 population transfer, the observed CoIn structure coincides with the S2 minimum geometry showing a delocalized excitation of both phenyl moieties. S1 → S0 relaxation occurs localized on a single phenyl moiety via structures similar to the reported CoIn seam of benzene.56,57 At the CASSCF(12,12) level of theory (Figure 3, dashed line), population decay between adiabatic states is slightly retarded: a monoexponential fit yields S2 → S1 decay characterized by a 79 fs timeconstant while a biexponential fit reveals two components, an ultrafast 17 fs and a slower 199 fs component where 53% of the population decays via the ultrafast channel. Within the 300 fs propagation time, population of the electronic ground state was not observed. The increased S2 lifetime is a consequence of the superior description at the CASSCF(12,12) level. The accurate description of excitation energies (cf. Table 1) results in a lower kinetic energy during the trajectory dynamics. Accordingly, a smaller phase space is sampled and displaced S1 → S0 hopping geometries are not reached. 4.2. Spectral Signatures of Nonradiative Relaxation in DPM. In the following we present the simulated absorptive 2D spectra of DPM employing the QF-NEP algorithm. According to the algorithm introduced in section 2.2 the signal of a single initial condition incorporates the dynamic effects of three MD trajectories run in either participating electronic state thus accounting for population and coherence contributions to the signal. To keep the computational cost feasible we neglect the effects of system dynamics on lineshapes during t1 and t3 in spectra simulation, while MD results affect t2 dynamics (for computational details see SI). 4.2.1. Single Initial Condition 2D Spectra. Even though not an observable, we start by discussing the 2D signal of an exemplary initial condition to demonstrate the microscopic information contained in the 2D spectra (Figure 4). At t2 = 0 fs, the diagonal peaks arise from population contributions, whereas offdiagonal peaks originate from coherence contributions to the signal (Figure 4a).58 Peaks at Ω1 = Ω3 = 37500 and 40500 cm−1 mark excitation energies of S1 and S2, respectively, that have a particular large energy gap for this initial condition. This again

Figure 3. Time evolution of adiabatic state populations for trajectories launched in S1 and S2. Solid lines: CASSCF(8,8) level of theory (120 trajectories). Dotted lines: CASSCF(12,12) level of theory (68 trajectories). Light colored lines: Fitted population. Inset: Ground state equilibrium structure of DPM.

Table 1. DPM Excitation Energies (in eV) and S1/S2 Splittings (ΔS1/2 in cm−1) at the CASSCF(8,8) Equilibrium Geometry (C2 Symmetry, 6-31G*Basis Set) S1 CASSCF(8,8) CASSCF(12,12) SS-CASPT2(8,8) MS-CASPT2(8,8) SS-CASPT2(12,12) LT-DF-LCC2† DPMb DPM-d12a,c CASSCF(8,8)c CISc TD-DFT/B3LYPc EOM-CCSDd a

this work 6.27 5.03 4.81 5.07 4.56 5.14 experiment 4.63 4.65 theory 6.87 5.86 5.04

S2

ΔS1/2

6.28 5.05 4.85 5.10 4.64 5.18

113 159 270 265 591 311

4.64 4.66

123 116

6.88 5.95 5.22

136 697 1429 430

DPM-d12: Deuterated derivative. † Basis: cc-pVDZ. bRef 54. cRef 10. Ref 11.

d

reported splitting (113 vs 123 cm−1, Table 1). However, because of the restricted active space and lack of dynamic electron correlation the absolute excitation energies are about 1.6 eV to high. Increasing the active space by considering all π orbitals of both phenyl (CASSCF(12,12)) accounts for 75% of electron correlation with respect to experiments (S1/S2 splitting: 159 cm−1). CASPT2(8,8) calculations further improve the excitation energy by 0.2 eV, while MS-CASPT2(8,8) and LT-DF-LCC2 overestimate the excitation energies by about 0.45−0.50 eV and are not superior to the CASSCF(12,12) method. It is interesting to note that the formally superior MS-CASPT2(8,8) excitation energies are outperformed by the SS-CASPT2(8,8) method. Given the moderate quality of the employed basis set (6-31G*) we assign this behavior to error cancelation effects also observed in other molecules.53 Finally, SS-CASPT2(12,12) calculations are in almost quantitative agreement with experiment. The S1/S2 energy splittings at CASPT2 and LCC2 level of theory (265−591 cm−1) are generally larger than the observed experimental value. This is in agreement with refs 10 and 11, which find a strong dependence of S1/S2 splittings on the level of theory up to 1429 cm−1 (cf., Table 1). Recently, Ottiger et F


Article


Figure 4. Exemplary initial condition absorptive 2D spectra for delay times t2 = 0 fs (a) and (b) t2 = 10 fs. (CASSCF(12,12) level of theory for population contributions; CASSCF(8,8) level of theory for coherence contributions shifted by 8680 cm−1 to match CASSCF(12,12) excitation energies; exponential lifetime broadening τl = 20 fs; nuclear decoherence is neglected; for details see SI).

demonstrates the high sensitivity of the S2/S1 splitting on the employed geometry. The peak intensities are a result of the utilized z-polarization of pulses that favors S2 contributions over S1 for configurations near the Smin 0 . Accordingly, the diagonal peak at Ω1 = Ω3 = 40500 cm−1 carries the highest intensity, whereas coherence and S1-population contributions appear less pronounced. Already at t2 = 10 fs the evolution of molecular geometry affects the spectrum (Figure 4b). Because of nuclear relaxation the S2-population contribution appears shifted along the detection axis at (Ω1, Ω3) = (40500, 39000) cm−1, while the S1-population contribution appears as distinct resonance at (Ω1, Ω3) = (37500, 37000) cm−1. As in our model, the coherence contributions ((Ω1, Ω3) = (37500, 43500) and (40500, 40000) cm−1) follow energy fluctuations of S0, the respective contributions to the 2D signal do not show the relaxation behavior of population contributions. 4.2.2. Ensemble Averaged Linear Absorption and Nonlinear 2D Spectra. Figure 5a presents the computed ensemble averaged linear absorption spectra calculated on various levels of theory (color lines) together with the experimental spectrum9 that exhibits a broad asymmetric peak around 38600 cm−1 and a small peak around 37400 cm−1 in gas phase and cyclohexane solution (black dotted and solid lines, respectively). Incoherent averaging over 500 initial conditions (blue line, CASSCF(8,8)) yields a symmetric spectrum with a maximum around 38500 cm−1 (note that because of the limited active space simulated CASSCF(8,8) spectra are shifted by 10000 cm−1, cf. Table 1). Taking into account ground state energy gap fluctuations introduces asymmetry for higher energy components of the linear absorption spectrum now peaking at 39100 cm −1 (light green, 60 MD trajectories). The corresponding spectrum obtained from propagation in excited states S1 and S2, respectively (dark green line, 120 trajectories) better recaptures the asymmetric peak shape observed in the experiment with the peak maximum shifted toward lower energies. The inclusion of excited state relaxation dynamics leads to dispersive features on the single trajectory level that average out upon ensemble averaging, resulting in the observed asymmetry of the absorption peak, in agreement with results reported in refs.59,60 The different levels of signal simulation demonstrate the effect of ensemble averaging and intramolecular system dynamics already to the linear absorption signal. While the incoherent averaging over initial conditions only allows to recover the width of the experimental spectrum,

Figure 5. Linear absorption and nonlinear 2D spectra of DPM. (a) Blue: Lorenzian broadening (500 Wigner distributed initial conditions). Light green: NEP spectra of ground state trajectories. Dark green: QF-NEP spectrum of excited state trajectories (simulated spectra are shifted by 10000 cm−1). Black solid and dotted lines: Experimental spectra in solution and gas phase.9 (b−d) Absorptive 2D spectra for t2 = 0 fs: (b) CASSCF(8,8) (τl = 20 fs, 60 initial conditions), (c) CASSCF(8,8) (τl = 5 fs, 60 initial conditions), and (d) CASSCF(12,12) (τl = 5 fs, 34 initial conditions).

the inclusion of dynamic effects systematically improves the envelope peak shape. Ensemble averaged absorptive 2D spectra of DPM for t2 = 0 fs are presented in Figure 5b−d, where we compare signals calculated on the CASSCF(8,8) level with the signal calculated on the more accurate CASSCF(12,12) level of theory. Using a dephasing lifetime τl = 20 fs (Figure 5b) still allows to reveal contributions from individual trajectories. The elliptic shape along the diagonal is caused by the inhomogeneous excitation energy distribution of different initial conditions. Offdiagonal peaks are determined by S1/S2 splitting at the respective initial condition leading to an antidiagonal broadening of 2D spectra. Here we want to point out that the S1/S2 gap itself is strongly affected by nuclear configuration and is thus not a direct measure of excitonic coupling between phenyl moieties. A reduction in lifetime (Figure 5c, τl = 5 fs) results in a single, broad resonance in the 2D signal showing an elliptic elongation along the diagonal. Antidiagonal peak broadening due to S1/S2 energy gap is now obscured by lifetime broadening but still noticeable (note the asymmetry of signal). The 2D spectrum obtained from CASSCF(12,12) trajectories (Figure 5d) shows a peak maximum located around Ω1 = Ω3 = 40000 cm−1 (unshifted value) in good agreement with experimental absorption spectrum that appears more compact along the diagonal reflecting the reduced inhomogeneous excitation energy distribution. The laser polarization (zzzz) is chosen such that S2 contributions to the spectrum are highlighted, reflected by an asymmetry of the peak with the maximum slightly shifted to the blue end of the spectrum. 4.2.3. Population Dynamics For Varying t2 Delays. Increasing t2 affects the spectra via a multitude of timedependent effects induced by nonadiabatic relaxation and G


Article


Figure 6. Absorptive 2D spectra of DPM for t2 delays 3−300 fs. Top: CASSCF(8,8) (τl = 5 fs, 60 initial conditions). Bottom: CASSCF(12,12) (τl = 5 fs, 34 initial conditions).

Figure 7. Excitation energy gap distributions and average of S1 (red) and S2 (blue) excitation energies at the CASSCF(8,8) (a) and CASSCF(12,12) level of theory (b). (c) Peak position along Ω3 of absorptive 2D spectra (CASSCF(12,12), τl = 2 fs). (d) Fourier transform (FT) of oscillatory component of (a). (e) FT of panel b. (f) FT of panel c. Black sticks mark frequencies of vibrational modes of ground state DPM.

onset of spectral diffusion that alters the initially diagonal tilted peak to a roundish shape can already be recognized. Energy gap fluctuations originating from fast C−H stretching vibrations around 3300−3400 cm−1 are responsible for the ultrafast modulation of the peak shape. This ultrafast spectral diffusion

nuclear dynamics. Figure 6 presents the absorptive 2D signal for t2-delay times 3, 10, 100, and 300 fs. From the time evolution of adiabatic state populations (cf., Figure 3) we expect the onset of population transfer at ∼10 fs which after 100 fs should substantially affect the spectra. For t2 = 3 fs the H


Article


DPM.63 The corresponding Smin 0 harmonic vibrational frequencies were calculated to be 1056 and 1061 cm−1 (CASSCF(8,8)). By a Fourier transform of the oscillatory component of energy gap fluctuations of the trajectory ensemble these frequencies can be approximately recovered with dominant peaks around 1000−1100 cm−1 (Figure 7d and e). Further important modes appear as side bands at 800−850 and 1250− 1300 cm−1 and are assigned as ring distortion and CH2−modes affecting the bending angle between the phenyl moieties (cf., Figure SI.4). By analyzing the time-dependent peak shifts along Ω3 observed in the 2D spectra the same oscillatory behavior can be recovered (Figure 7c, CASSCF(12,12)). The respective Fourier transform yields a clear peak around 1000 cm−1 and weaker peaks at 800 and 1300 cm−1 (Figure 7f), in close agreement with modes extracted from the underlying energy gap oscillations. Activated ring breathing/distortion modes necessary to reach the relaxed nuclear configurations of S1 and S2 are thus predicted to have clear signatures in the 2D signals.

can be identified similarly on both, the CASSCF(8,8) and CASSCF(12,12) level of theory. At t2 = 10 fs, relaxation of excited states by nuclear motion affects the signal by shifting the main resonance below the diagonal. For signals evaluated at the CASSCF(8,8) level the corresponding 2D spectra spreads over a range of probe frequencies (Ω3), spanning about 10000 cm−1. Accordingly two distinct peaks can be identified around (Ω1, Ω3) = (48000, 40000) cm−1 and (Ω1, Ω3) = (50000, 45000) cm−1 where the lower energy resonance corresponds to hot S1 and the higher energy resonance corresponds to hot S2. In contrast, signals evaluated at the CASSCF(12,12) level show a narrower, compact signal with a single resonance that covers

Simulation of Multi-Dimensional Signals in the Optical Domain

Recommend Documents