Investigating ultrafast two-pulse experiments on single DNQDI fluorophores: a stochastic quantum approach

Abstract

Ultrafast two-pulse experiments on single molecules are invaluable tools to investigate the microscopic dynamics of a fluorophore. The first pulse generates electronic or vibronic coherence and the second pulse probes the time-evolution of the coherence. A protocol that is able to simulate ultrafast experiments on single molecules is applied in this study. It is based on a coupled quantum-mechanical description of the fluorophore and real-time dynamics of the system vibronic wave packet interacting with an electric field, described by means of the stochastic Schrödinger equation within the Markovian limit. This approach is applied to the DNQDI fluorophore, previously investigated experimentally [D. Brinks et al., Nature, 2010, 465, 905–908]. We find this to be in good agreement with the experimental outcomes and provide microscopic and atomistic interpretation.

---

1 Introduction

In the last decade, ultrafast spectroscopy on single molecules has become a powerful tool to investigate different features of microscopic dynamics as quantum coherence.^1–8 A possible strategy to study processes at the atomic scale is to set up techniques which have appropriate sensitivity. Ultrafast spectroscopy employs ultrashort pulses whose duration has the same time scale (ideally shorter) as the processes one wants to investigate.⁹ Pulses resolved at the femtosecond timescale are essential to detect very fast molecular processes such as electronic dephasing, charge transfer, excitation energy transfer, and vibrational energy relaxation.^10,11 Many of these ultrafast processes play a relevant role in biological processes: indeed one of the applications of ultrafast experiments is the detection of electronic and vibrational coherences^12,13 in light-harvesting complexes involved in photosynthesis.^14–16

In order to detect femtosecond molecular processes, different ultrafast techniques based on ensemble measurements have been developed. In particular in the framework of non-linear optics, multidimensional electronic spectroscopy has allowed the investigation of the vibrational structure of a system and the detection of coherence decay.^17,18 Different from multidimensional spectroscopy, single-molecule techniques^{1–3,19–23} have the advantage of preventing ensemble-average effects, detecting the intrinsic inhomogeneity due to different conformations and local environments experienced by each molecule. The outcomes of two-pulse ultrafast experiments may be strongly affected by quantum coherence whose persistence depends on the interaction between the system (i.e., the molecule under study) and the environment. For this reason the effects of the surrounding environment cannot be neglected, and the system is described as open.¹⁷

Approaches based on the time-dependent Schrödinger equation have been widely employed to perform a time-resolved analysis of the optical processes of interest, also coupled with a quantum chemical description of the molecule,^24–26 usually focusing on a close system. A description of the spectroscopic target in terms of open quantum systems^27,28 is required in order to include the presence of a surrounding environment. Most of the approaches present in the literature are based on density matrix master equations^29–32 involving phenomenological parameters,^33,34 in some cases including a vibronic analysis.^35,36 The density matrix approach has been also coupled with a TDDFT^37,38 or quantum-chemistry-based³⁹ description of the spectroscopic target. The protocol here employed is based on existing methodologies that, when combined, allow the investigation of ultrafast processes preserving a refined description of the fluorophore. The fluorophore time propagation is computed through the stochastic Schrödinger equation (SSE)^40–42 within the Markovian limit.^43,44 In this model, dissipation and fluctuations due to the environment are included by means of stochastic terms in every realization of the wave function evolution. Averaging over a number of realizations leads to results equivalent to the Lindblad master equation solution with the advantage of less computational effort.⁴⁵ The density matrix terms (i.e. populations and coherences), computed starting from the coefficients of wave function expansion, are used to interpret the experimental results.

The SSE is coupled with a quantum mechanical description of the vibronic structure of the molecular target, such as vibronic energies and transition dipole moments. This information can be obtained by fitting experimental results. However, this method usually leads to no univocal sets of input data. On the other hand a quantum chemistry computation of the fluorophore leads to a univocal system description. In the present work the electronic ground and first-excited states of DNQDI have been characterized in terms of energies, dipole moments and normal mode frequencies employing DFT/TDDFT methods. A vibrational structure has been included by means of harmonic vibrational levels and related Franck–Condon factors.

The whole protocol, schematized in Fig. 1, has been applied in the present work to a two-pulse experiment performed on the DNQDI fluorophore,¹ in order to clarify the interpretation of the experimental results. In the experiment,¹ fluorophore molecules embedded in a thin polymeric matrix interact with two identical pulses 15 fs long separated by a controlled delay time Δt of the order of femtoseconds. Moreover a phase difference Δϕ between the two pulses is applied. Ultrashort pulses have a bandwidth large enough in the frequency domain to excite various vibrational levels of the electronic excited state. The fluorescence signal of single molecules (proportional to the population in the excited state after the two pulses) has been detected for fixed phase differences. Their profile as a function of delay time has a dumped oscillating behaviour with frequencies distributed between 670 and 1500 cm⁻¹ and a peak at 1070 cm⁻¹.¹ The signal decay is the result of pure electronic and/or vibronic dephasing of the wave packet:^20,46 typical dephasing times T₂ for this molecule under the experimental conditions reported in ref. 1 are around 60 fs. Oscillations were explained in terms of wave packet interference and interpreted as an excited-state vibrational signature of DNQDI.¹ In detail, the first pulse populates a manifold of vibrational levels of the excited state, generating a vibronic wave packet and it creates a coherence between the electronic ground state and the excited state, which evolves with time. The second pulse interacts with the system generating another wave packet that can interfere constructively or destructively with the first one, depending on Δϕ and Δt values. The excited state population after the second pulse depends on such interference that in turn is sensitive to the decoherence taking place during Δt. Thus the emission signal provides signature of the coherence dynamics.

Fig. 1 Scheme of the protocol employed to simulate ultrafast experiments on single molecules.

This paper is organized as follows: the theoretical framework employed to simulate the ultrafast dynamics induced by two pulses is described in Section 2, computational details are briefly summarized in Section 3, the present results are reported and discussed in Section 4, and our final comments and perspectives are summarized in conclusions.

2 Theory

Our theoretical and computational approach is hierarchically based on three steps (Fig. 1): (i) description of the electronic structure of the spectroscopic target by means of DFT/TDDFT; (ii) electronic properties (excitation energies and transition dipole moments) are then “dressed” by a proper vibrational structure, using the quantities computed at step (i); (iii) a real-time propagation of the vibronic wave packet is performed, in terms of close (“standard” time-dependent Schrödinger equation) and open quantum (SSE) systems, using the quantities from steps (i) and (ii) as input parameters.^43,44 Using DFT/TDDFT to characterize states for many-electron dynamics has been performed before.^24,25,47

In this section, we briefly review the SSE protocol and the methodology to vibrationally dress the electronic states.

2.1 Stochastic Schrödinger equation

The stochastic Schrödinger equation (SSE)^48–50 is an approach equivalent to master equation for the density matrix, as Lindblad equation,⁵¹ that is able to perform the real-time propagation of the system wave function.⁴⁵ The SSE employed in this work is thoroughly described in ref. 43, while we report here only the main equations and the new features included. The SSE in a Markovian regime readswhere c runs over the number M of relaxation and dephasing channels of the system. The first term in the right-hand side of the equation describes the unitary evolution of the time-dependent Schrödinger equation while the second term describes the dissipation due to the presence of the bath. The stochastic behaviour of eqn (1) lies in the third term which depends on the differential noise in a Wiener process l_c(t), simulating fluctuations of the system due to the environment.⁵² The average over a great number N_traj of wave function realizations j reproduces the time evolution of the reduced density matrix (only containing the system degrees of freedom)The system wave function is expanded on a basis of vibronic states |Φ_i〉 aswith N the number of field-free eigenstates |Φ_i〉 included in the expansion. N is equal to , with N_el, N_m and N_vib being the number of electronic states, the number of included normal modes, and the number of vibrational states per normal mode, respectively. Eqn (1) returns the coefficients of the wave function expansion, allowing one to compute all the elements of the density matrix at a certain time t, averaging over the N_traj SSE trajectories.⁴³ Diagonal elements of the density matrix are the population of stateswhile off-diagonal elements are defined as coherences. Time-dependent coefficients of the wave function are propagated by a deterministic second-order Euler algorithm, coupled to a quantum jump algorithm^53–55 simulating the stochastic term in eqn (1).⁴³

The Ŝ_q operators in eqn (1) describe the relaxation and dephasing processes. In the model developed before⁴³ relaxation has been included via radiative (spontaneous) emission and nonradiative decay processes. According to previous studies,^43,44,48,56 relaxation has been included by the operator

where |Φ^g₀〉 and |Φ_q〉 are respectively the vibronic ground state (product between the vibrational and the electronic ground state) and a generic vibronic excited state q, i.e., the product of N_m vibrational states and an electronic excited state. The population ([small rho, Greek, circumflex]_S(t))_qq of the q-th state exponentially decays with the rate constant Γ_q. Such decay rates can be included as phenomenological parameters for non-radiative processes or according to Fermi's Golden rule in fluorescence processes. Here we choose Γ_q = T₁⁻¹, with T₁ phenomenological relaxation time, for all the q states.

Vibrational and electronic coherences decay with different time scales.⁵⁷ Following the treatment in ref. 58, the (pure) dephasing of the quantum state of the molecular target originates from the random time fluctuations of the energy difference between the states of interest, due to the interaction with degrees of freedom not treated explicitly. For electronic levels of a molecule in a matrix, the differential interaction of the ground and excited-state charge densities with the environment fluctuating polarization is an important source of excitation energy fluctuations and thus electronic dephasing. Vibrational dephasing may have various origins,⁵⁹ but the anharmonic coupling among vibrational modes (of the composite molecule-environment system) is certainly relevant in modulating vibrational frequencies and thus inducing dephasing.⁶⁰

Here we include them with different operators. For vibrational dephasing we exploit the generic operators defined before,⁴³ specialized for vibronic states k and q:

where M_kk terms are equal to 1 with the exception of the element M_qq = −1. With this choice, vibronic state populations are not affected by the dephasing operator even at the level of single trajectory. This definition of dephasing operator produces coherences (ρ_S(t))_pq that decay with a rate γ_pq = γ_p + γ_q. Here we choose a single dephasing time for all the vibronic state pairs, i.e., we set γ_pq = γ_vib = T₂(vib)⁻¹ for any p and q, where T₂(vib) is the vibrational dephasing time.

To account for a purely electronic dephasing (i.e., a dephasing that does not affect the coherence between different vibrational states belonging to the same electronic state), in the present work we also introduce the following dephasing operator:

where Φ^e_k is the vibronic state obtained as the product of the k vibrational state (for a given normal mode) and the electronic excited state (e), and Φ^g_k is the vibronic state obtained as the product of the k vibrational state (for a given normal mode) and the electronic ground state (g). The operator in eqn (7) is specific for two electronic states (ground and excited), but it has been implemented for a general number N_el of electronic states. The sum in eqn (7) extends to each normal mode. γ_el in eqn (7) is related to the electronic dephasing time T₂(el): γ_el = T₂(el)⁻¹.

The form of the two operators in eqn (6) and (7) assures that coherences between vibrational levels of the same electronic state exponentially decay with a rate equal to γ_vib (i.e., only vibrational dephasing matters) whereas, for coherences between vibrational levels of different electronic states the decay rate is γ_vib + γ_el (i.e., both vibrational and electronic decoherences are at work).

2.2 Vibronic structure

Electronic-state potential energies are approximated by harmonic oscillators, as shown in the upper panel of Fig. 2.⁶¹ The transition dipole moments are calculated assuming the Franck–Condon (FC) approximation and the integrals between vibrational states are computed using the Sharp and Rosenstock approach.^61–65 It allows the explicit calculation of the 〈0_g|0_e〉 FC integral taking into account the displacement of the equilibrium position due to all the normal modes. The FC integrals between the vibrational ground state of the electronic ground state and any other vibrational levels in an excited electronic state are evaluated recursively from 〈0_g|0_e〉.

Fig. 2 Schematic structure of the DNQDI vibronic system and molecular structure of the reduced model of DNQDI employed in the simulations.

When more than one normal mode is considered to describe the molecule, all the possible combinations between different quanta of different normal modes give a total number of states with N equal to , as already shown. Including the complete set of normal mode frequencies in the vibrational structure of a large system translates into prohibitive computational costs of the wave function propagation. Therefore, a selected number of vibronic states can be included based on the analysis of the equilibrium position displacement vector (within the harmonic approximation), as it will be shown in Section 4. Approaches able to describe the molecular vibronic structure beyond the harmonic approximations and including non-adiabatic coupling can be found in the literature.^66–68 However, since pulses employed in the target experiment are temporarily very close (largest delay time of 120 fs), slower processes driven by non-adiabatic coupling are unlikely. For DNQDI in particular, nonradiative decay induced by non-adiabatic coupling takes place on the much longer ns time scale. For this reason, neglecting non-adiabatic coupling is considered accurate for our goal.

3 Computational details

In order to vibrationally characterize the molecule, ground and excited state energy and vibrational normal modes of DNQDI have been calculated using DFT and TDDFT methods at the B3LYP/6-31G(d) level of theory performed by means of Gaussian16.⁶⁹ Quantum mechanical calculations were performed for a reduced model of the fluorophore (lower panel of Fig. 2) where alkyl lateral chains have been replaced with hydrogen atoms in order to decrease the computational cost (see Fig. S1 in the ESI†).

The vibrationally resolved optical spectrum has been computed with the FCclasses code.⁶¹ FCclasses allows the calculation of the displacements of excited-state normal modes coordinates with respect to the ground state normal modes coordinates and FC integrals. FC integrals have been computed between the vibrational ground states of the electronic states 〈0_g|0_e〉 and also between the ground and the vibrational excited states of the excited singlet state 〈0_g|ν_e〉, up to 19th vibrational quantum number.

SSE simulations were performed by using the homemade WaveT code.⁴³ A time step equal to δt = 1.21 as has been employed while the whole dynamics takes 1 ps. A shorter dynamics, e.g., around 200 fs, would have been large enough to properly reproduce the ultrafast dynamics of DNQDI, since delay times are limited to 120 fs. Yet, in the initial numerical assessment of our approach we tested larger delay times that required 1 ps duration, and we kept this value for all the simulations. When specified, relaxation is included in the simulation viaeqn (5) assuming a fluorescence emission time T₁ = 3 ns (corresponding to Γ_q = 3.3 × 10⁻¹ ns⁻¹ for each q). Unless specified otherwise, dephasing is included viaeqn (6) with γ_q = 2.01 × 10⁻⁴ fs⁻¹ for each q. This decay rate leads to a T₂ of 60 fs, chosen on the basis of a phenomenological analysis of experimental data.¹

100 SSE trajectories were calculated for each delay time and then the averaged populations of the excited states were computed. When compared with the outcomes of Lindblad master equation (Fig. S2 of ESI†), the results with 100 SSE trajectories show a relative error of less than 1%. A numerical test carried out with 200 SSE trajectories (Fig. S2, ESI†) confirms that accurate and reliable results are obtained by using 100 SSE trajectories.

All the dynamics start from the vibrational ground state, since vibrational excited states for the normal modes of interest (reported in Table 1) are not thermally populated.

Table 1 Normal mode frequencies in the ground (GS) and excited (ES) electronic states for normal modes with the largest displacement

GS frequency (cm^−1)	ES frequency (cm^−1)	Displacement (Å)
1368	1370	7.52 × 10^−2
1629	1632	4.53 × 10^−2
1265	1271	3.56 × 10^−2
1647	1638	3.14 × 10^−2
1583	1570	2.65 × 10^−2

---

The field parameters are modeled after the experimental set-up.¹ Pulses have an enveloped Gaussian sinusoidal shape, with the FWHM equal to 15 fs to have a bandwidth of 120 nm while the maximum intensity of the electric field is I_max = 1 × 10⁹ W cm⁻², unless otherwise specified. The pulse frequency ω_exc employed in the experiment is shifted from the maximum absorption energy by 507 cm⁻¹. Assuming that the calculated vertical absorption energy is the maximum absorption energy, a detuning of 1090 cm⁻¹ between ω_exc and the |0_g〉 − |0_e〉 transition can be estimated in the simulations. We have used such detuning to set the central frequency of the pulse (1.76 eV, 703 nm).

4 Results

In this section our results are presented. We first show the simulation of the vibrationally resolved optical absorption spectrum of DNQDI and its comparison with the experimental one. Then we report the results of the simulations of the real-time dynamics of the system interacting with two delayed pulses, compared with the experimental findings. Simulations have been performed by assuming either a pure electronic system or including a vibronic structure. We have systematically investigated several key aspects: the vibronic structure of the fluorophore, the possible interplay of electronic and vibrational dephasing effects, and the roles of pulse shape and pulse intensity.

4.1 Optical spectrum

Fig. 3 reports the comparison between the computed (our work) and the experimental⁷⁰ vibrationally resolved absorption spectrum of DNQDI (the experimental one has been shifted in order to superpose the two spectra). Although the simulated spectrum (in vacuo) is rigidly red-shifted by 0.13 eV with respect to the experimental one⁷⁰ (in toluene), the general features are properly maintained: the distance between the two peaks is around 1370 cm⁻¹ in both spectra, which is a typical value of vibrational modes in aromatic compounds.⁷¹ Generally speaking, the lower-energy peak may be related to the |0_g〉 − |0_e〉 transition (from the vibrational ground state of the electronic ground state to the vibrational ground state of the electronic excited state) and the other one to a |0_g〉 − |1_e〉 transition (from the vibrational ground state of the electronic ground state to the first vibrational state of the electronic excited state). We also calculated the absorption energy of DNQDI in toluene within the PCM model, red-shifted to 0.24 eV with respect to the maximum absorption in the experimental spectrum. The PCM calculation of absorption energy of DNQDI in PMMA (the matrix polymer employed in the ultrafast experiment¹) provided the same result as in toluene. The analysis of the displacements of the excited-state parabolas along each one of the 3N − 6 normal mode coordinates (with N being the number of atoms) allows us to establish which normal modes contribute to the |0_g〉 − |1_e〉 band of the absorption spectrum. Indeed only if a normal mode is characterized by a sufficiently large displacement of the excited-state parabola with respect to the ground state, the vibrational ν = 1 state (|1_e〉) of the electronic excited state has a non-negligible overlap with the ground state |0_g〉. As one can see in Table 1, the normal mode with the largest displacement (7.52 × 10⁻² Å) has a frequency of 1368 cm⁻¹ (1370 cm⁻¹) in the electronic ground (excited) state, fully compatible with the separation of the two peaks in the absorption spectrum of Fig. 3. The two normal modes with the largest displacement (with frequencies 1368 and 1629 cm⁻¹ in the GS) have been considered for the vibronic dynamics, reported below in Section 4.3.

Fig. 3 Experimental⁷⁰ (exp., in toluene) and calculated (calc.) absorption spectrum of DNQDI. The calculated spectrum has been blue-shifted by 0.13 eV in order to superpose the experimental and simulated spectra.

4.2 Pure electron dynamics

We first propagated SSE trajectories, only involving the ground and the first-excited electronic states; the pulse bandwidth, taken from the experimental work,¹ allows us to appreciably populate only the first excited state (as reported in Fig. S3 of the ESI†) since the other electronic excited states are too high in energy (Table S1, ESI†). Two-photon absorption under the conditions of the experiment is also unlikely, as discussed in the ESI.† The delay time has been varied between 0 fs and 120 fs and Δϕ is equal to 0 or π. SSE real-time dynamics were performed including dephasing and relaxation effects. The excited-state population, after the pulse sequence, as a function of delay time has a damped oscillating profile as shown in panel (a) of Fig. 4. The decay of population is due to the dephasing^1,43 with decay time T₂(el) = 60 fs, as a consequence of the interference (constructive or destructive) between the two pulses interacting with the system. On the other hand the dynamics is short enough to prevent population relaxation, and thus the effect of decay time T₁ = 3 ns is negligible. Moreover the excited-state population oscillates with a frequency (1090 cm⁻¹, 0.13 eV) equal to detuning δ = |ω_exc − ω_0–0| (ω_0–0 is the frequency of the 0–0 excitation). The FT of the signal is reported in panel (b) of Fig. 4. The profile of excited state population is analogous to the experimental fluorescence signal reported in panel (c) of Fig. 4. A vibrational structure is not therefore needed to produce fluorescence oscillations: also for a pure electronic two-level system the population of the excited state oscillates with the detuning frequency. This was already noted in ref. 72.

Fig. 4 (a) Population of excited states at the end of simulation as a function of delay time in pure electronic dynamics. (b) FT of (a). (c) Fluorescence emission signal as a function of delay time from ref. 1. ω_0–0 is the excitation energy of the excited state, ω_exc is the pulse frequency.

4.3 Vibronic dynamics

The next step in our investigation was the inclusion of vibrational states in the molecular description. We took into account only the two normal modes responsible for the largest displacement along the nuclear coordinate, as shown in Table 1. These two normal modes, schematically represented in Fig. S4 and S5 of ESI† and labeled as NM1 and NM2, are characterized by the “breathing” of the aromatic structure and by the C–H bending.

For each normal mode, 20 vibrational levels were considered to compute FC integrals. Transition dipole moments between the state levels and ground vibrational states of electronic ground and excited states were calculated within the FC approximation. According to the values of the FC integrals for the two selected normal modes (Table 2), vibrational levels higher in energy were neglected in the SSE propagations. The 〈0_g|0_e〉 contribution is dominant for both normal modes.

Table 2 FC integrals between |0_g〉 and |v_e〉 with ν = 0,…,6 for the two normal modes responsible for the largest displacement

State (νe)	NM1 〈0g|νe〉	NM2 〈0g|νe〉
|0e〉	0.222	0.222
|1e〉	7.43 × 10^−2	4.89 × 10^−2
|2e〉	1.76 × 10^−2	7.78 × 10^−3
|3e〉	3.40 × 10^−3	1.03 × 10^−3
|4e〉	5.69 × 10^−4	1.21 × 10^−4
|5e〉	8.5 × 10^−5	1.3 × 10^−5
|6e〉	1.2 × 10^−5	1 × 10^−5

---

4.3.1 Vibronic wave packet with NM1. The two electronic states were dressed with 4 vibrational levels each for the normal mode NM1, according to the negligible value of Franck–Condon integrals for ν_e ≥ 4 in Table 2. Wave packet dynamics was computed with and without dephasing and relaxation effects. Fig. 5 shows the total excited-state population (i.e. P_0e + P_1e + P_2e + P_3e, where P_ke is the population of the k vibrational level of the electronic excited state) as a function of delay time for Δϕ = 0 and π, without dephasing (panel (a)) and with dephasing (panel (b)) with their corresponding FTs. The damped profile of excited-state populations as a function of delay time, when SSE is employed, is due to the presence of dephasing, as already seen for the pure electronic case in Fig. 4. The influence of higher-energy vibronic states is very small: indeed, the transition dipole moment is significantly lower (one order of magnitude) than (see Table 2). Including the population of state |1_e〉 produces a modulation of the original |0_e〉 population as could be noticed in panel (a) of Fig. 5.

Fig. 5 Sum of excited-state populations at the end of simulation as a function of delay time without (panel (a)) and with (panel (b)) relaxation and dephasing effects. (c) FT of (a). (d) FT of (b). Only the NM1 normal mode was added to the wave packet of eqn (3). ω_0g−0e (ω_0g−1e) is the excitation energy of state |0_e〉 (|1_e〉), and ω_exc is the pulse frequency.

The FT of the sum of populations (panel (c) of Fig. 5) reveals a frequency peak of around 1090 cm⁻¹, namely the already observed detuning δ discussed previously. A smaller frequency peak (about 280 cm⁻¹) is also present, related to detuning between pulse frequency and excitation frequency of state |1_e〉: it has a lower intensity because the |1_e〉 state is much less populated than |0_e〉. The same results are observed with the FT of the SSE signal (panel (d) of Fig. 5), with only a broadening of the frequency peaks. Signals with Δϕ = 0 and π give approximately the same information.

4.3.2 Vibronic wave packet with NM1 and NM2. The same analysis was performed by also including the normal mode NM2, dressing each electronic state with vibrational levels. Fig. 6 shows the sum of excited state populations at the end of simulation as a function of delay time for simulations obtained including dephasing and relaxation effects. The result is very close to panel (b) of Fig. 5 since the contribution of NM2 is smaller in terms of transition dipole moment than the contribution of NM1 (see Table 2). Moreover, the profile of the sum of excited-state populations is not very different with respect to that of the pure electronic system shown in panel (a) of Fig. 4: FC integrals between state |0_g〉 and excited vibrational levels of the electronic excited state are much smaller than the FC integral 〈0_g|0_e〉. Vibronic dynamics does not change the interpretation: the observed oscillations in the fluorescence profile are due to detuning rather than to a normal-mode frequency of DNQDI.

Fig. 6 Sum of excited-state populations at the end of the SSE simulation as a function of delay time. Normal modes NM1 and NM2 were added to the wave packet of eqn (3).

4.3.3 Vibrational and electronic dephasing. In the vibronic dynamics reported in Sections 4.3.1 and 4.3.2, dephasing was included via the operators defined in eqn (6), i.e., coherences between vibrational states belonging to the same or different electronic states decayed with the same rate. To obtain a deeper insight into the ultrafast vibronic processes occurring in the DNQDI fluorophore, when interrogated by a sequence of two pulses, we distinguished between electronic and vibrational dephasing⁵³ by using both the operators defined in eqn (6) and (7). Our aim here is to investigate if the decay of the oscillations in the population of the electronic excited state after the pulses is mainly due to vibrational dephasing^73,74 or electronic dephasing. We investigated this aspect by applying a vibrational dephasing time T₂(vib) = 60 fs and an electronic dephasing time T₂(el) = 10 fs to the molecular system in the SSE simulations. A second simulation has been performed by employing T₂(vib) = 150 fs and T₂(el) = 60 fs.

The sum of populations of excited states as a function of delay time is shown in Fig. 7. In both cases the decay is dominated by the electronic dephasing time. The result may be easily understood in terms of the faster electronic loss of coherence between the states |0_g〉 and |0_e〉 or |0_g〉 and |1_e〉, created by the first pulse. These numerical results show that under the conditions of the experiment, vibrational coherences (between |0_g〉 and |1_g〉 or |0_e〉 and |1_e〉), i.e., the coherent motions of the nuclear wave packet, do not play any appreciable role. Again, for both cases reported in Fig. 7 the fluorescence oscillations are produced by the dominant contribution of the detuning δ between pulse frequency and the excitation frequency for the |0_g〉 →|0_e〉 transition. A different value of detuning was also tested and the results of the corresponding simulations are reported in the ESI† (Fig. S6–S8).

Fig. 7 Sum of excited-state populations at the end of simulation as a function of delay time including dephasing (upper panel) T₂(el) = 10 fs and T₂(vib) = 60 fs, (lower panel) T₂(el) = 60 fs and T₂(vib) = 150 fs.

4.3.4 Role of pulse shape. We also explored the possible effect of the pulse shape on the emission profile. To better mimic the experimental shape of the pulse (Fig. 1a of ref. 1), we considered a more complicated time profile of the pulse (Fig. S9 of ESI†), whose FT (panel (a) of Fig. 8) is closer to the experimental one. The FT of the experimental pulse has a bandwidth of 120 nm centered at 676 nm while another peak at about 620 nm is present.¹ The field shape in these simulations is composed by the superposition of two Gaussian-enveloped sinusoidal functions, with central frequencies displaced by 1340 cm⁻¹ as in the experimental electric field. The first pulse is analogous to the one used before while the second pulse has a wavelength of 649 nm, FWHM = 49.2 fs and maximum intensity I_max = 1.1 × 10⁸ W cm⁻². Panel (a) of Fig. 8 shows the pulse spectrum, where the two pulses (called ω_exc1 and ω_exc2) and DNQDI frequencies of states |0_e〉 and |1_e〉 are highlighted: comparison with the experimental laser field of Fig. 1a of ref. 1 shows good agreement in terms of general shape and relative peak position.

Fig. 8 (a) Applied electric field of the pulse as a function of time, composed by sum of two Gaussian-enveloped sinusoidal functions, see text for details. (b) Sum of excited-state populations at the end of simulation as a function of delay time without including dephasing and relaxation, using the time-resolved pulse of the upper panel. (c) FT of (b). ω_0g−0e (ω_0g−1e) is the excitation energy of state |0_e〉 (|1_e〉), ω_exc1 and ω_exc2 are the pulse frequencies.

Vibronic dynamics in the presence of this pulse, only including NM1, was carried out. Since the goal is to investigate the possible effect of the pulse shape on the fluorescence oscillations, the calculation was performed without dephasing and relaxation effects. Panel (b) of Fig. 8 shows the total excited-state population (i.e. P_0e + P_1e + P_2e + P_3e) at the end of the simulation as a function of delay time: the profile is very close to panel (a) of Fig. 5. The observed oscillations have a frequency equal to the detuning between ω_exc1 and ω_0g−0e (ω_0g−1e), as already observed with the “simple” pulse shape (Fig. 8). Detuning between the added pulse frequency ω_exc2 and the transition frequencies |0_g〉 → |0_e〉 and |0_g〉 → |1_e〉 was not detected probably because of out of the considered frequency domain, and because the excitation to state |1_e〉 is included only in the tail of pulse frequency bandwidth, respectively. As a conclusion, varying the shape of the applied field does not modify our physical interpretation.

4.4 Higher-order responses

All the conditions explored in the previous paragraphs do not allow us to directly detect vibrational signatures of DNQDI encoded in the time profile of the fluorescence. Of course when the detuning is zero, then the vibrational frequency can be detected as a fluorescence signal due to the excitation of the off-resonance vibrational states, but can vibrational frequencies be revealed by two-pulse experiments on a single molecule even in the presence of a detuning? All the simulations reported so far have been carried out in a linear response regime (I_max = 1 × 10⁹ W cm⁻², used in the experimental work¹). We focus on the higher-order responses of DNQDI to the applied field, going beyond the linear regime. Additional components should be added to the present model in order to take into account high-intensity phenomena such as ionization, nonlinear optical responses and multi-photon processes. However, our aim is the interpretation of the experimental work,¹ performed in the linear low-intensity regime, and the study of the possible role of vibrational coherence in the modulation of the fluorescence signal. In all the following simulations a four-level model for DNQDI was considered: two electronic states and two vibrational levels from NM1 normal mode for each electronic state. Dephasing and relaxation effects were neglected.

4.4.1 Role of the field intensity. The maximum field intensity has been varied so that higher order responses appear in the population profile. Pulse frequency and field amplitudes were kept the same as described in computational details. By increasing the field intensity up to I_max = 2.5 × 10¹² W cm⁻² and I_max = 1 × 10¹³ W cm⁻² we obtained the fluorescence profiles (i.e., the sum of populations) reported in panels (a) and (b) of Fig. 9. Rabi oscillations become particularly clear at small delay times: indeed for delay time Δt = 0 fs the Δϕ = 0 peak is much lower. The appearance of Rabi oscillations points out that the nonlinear regime is achieved during excitation by the pulses.

Fig. 9 (a) Sum of excited-state populations as a function of delay time employing a field with I_max = 2.5 × 10¹² W cm⁻² (b) Sum of populations of the excited states as a function of delay time employing a field with I_max = 1 × 10¹³ W cm⁻². (c) FT of (a) Δϕ = π, (d) FT of (b). Δϕ = π.

The excited-state populations after the pulses are larger when the higher-intensity field is applied, as expected. The FT of the sum of excited-state populations is shown in panel (c) (I_max = 2.5 × 10¹² W cm⁻²) and panel (d) (I_max = 1 × 10¹³ W cm⁻²) of Fig. 9. In both cases, two peaks at 1090 and 280 cm⁻¹ are present, which are due to the detuning between pulse frequency and |0_g〉 − |0_e〉 transition frequency, and between pulse frequency and |0_g〉 − |1_e〉 transition frequency, as in the linear response. Only for I_max = 1 × 10¹³ W cm⁻², the FT of the time-resolved signal reveals the presence of a peak at around 1380 cm⁻¹, very close to the normal mode frequency in both electronic ground and excited states. An additional simulation on a model system neglecting vibrational levels in the ground state confirms that the detected frequency in the nonlinear regime is the excited-state vibrational frequency (Fig. S10 in ESI†).

4.4.2 Perturbative analysis of the non-linear results. We want now to rationalize the results of the numerical simulations using an analytical model, in the framework of time-dependent perturbation theory. First-order excited-state coefficients are given bywhere E⃑(τ) is the electric field (which includes ω_exc, labelled as ω in the following equations), ω_0e and ω_1e are respectively the excitation frequencies of states |0_e〉 and |1_e〉, and are the transition dipole moments from |0_g〉 and C_0g(0) is the coefficient of the ground state at the beginning of the simulation (t = 0). In the approximation of an electric field composed of two delayed Dirac functions, the first-order response of excited-state coefficients is a function of delay time and it depends on the detuning between pulse frequency and excitation frequency (see ref. 72) and on the corresponding transition dipole moment, as resulting from the integral solutionC_0g(0) = 1 is assumed as the initial condition (as done for all the simulations of the present work) and E₀ is the amplitude of electric field E⃑ related to field intensity as . Population of states |0_e〉 and |1_e〉 are given byassuming rotating wave approximation (RWA), for which very fast oscillating terms are neglected since they are not detectable. Within the linear-order response, populations of states |0_e〉 and |1_e〉 explicitly depend only on an oscillating term with detuning frequency as a function of delay time, thus confirming our numerical investigation.

Field-induced transitions between states |0_e〉 and |1_e〉 are not possible (the related frequency is not included in the pulse frequency band), therefore only in the third-order correction of the excited-state coefficients, with respect to the electric field, the vibrational frequency of the excited state appears. Indeed the third-order coefficients of excited states depend on the second-order coefficients of the ground state, including oscillating terms with detuning frequency with respect to excitation energy of both states |0_e〉 and |1_e〉:

In particular the terms of |0_e〉 coefficients oscillating with the normal mode frequency depend on the product μ_0g−1e²·μ_0g−0e and on the third power of the electric field amplitude E₀: as a consequence, the population of excited states depends on the sixth power of the electric field (third power of the intensity) and transition dipole moments.

The excited-state populations at the sixth order with respect to the electric field amplitude (within the RWA approximation) include oscillating terms with detuning frequency with respect to excitation energy of both states |0_e〉 and |1_e〉, and also include oscillations at normal mode frequency (ω_1e − ω_0e) in the electronic excited state:

In conclusion, perturbative analytical results agree with the SSE numerical result simulations and confirm that excited-state normal mode frequencies of a single molecule can be only detected in the nonlinear regime, except when, trivially, the pulse and the vibronic transition frequencies coincide.

5 Conclusions

We have set-up a protocol that is able to simulate ultrafast experiments on single molecules by coupling the DFT/TDDFT description of the fluorophore and the real-time dynamics of the system by means of SSE. Based on quantum mechanical calculations, we have introduced a vibrational structure in the system description in terms of harmonic-oscillator levels and FC factors. Analysis on single molecules is a powerful tool to study the effects which are usually hindered by the ensemble average,^35,36 as reported in Fig. S11 of the ESI.†

The protocol has been applied to simulate a two-pulse experiment performed on DNQDI.¹ Our results are in agreement with the experimental findings, providing also microscopic interpretation: the observed oscillations of the fluorescence signal as a function of delay time are due to detuning rather than the vibrational frequency under the conditions of the original work.¹ In this particular case the two frequencies are of the same order of magnitude (around 1090 and 1370 cm⁻¹, respectively). SSE simulations and the perturbative analysis show that the oscillation at the detuning frequencies is given by the linear response of excited-state populations as a function of delay time,⁷² while the excited state vibrational frequency can be directly detected only at the nonlinear regime.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

Funding from the ERC under the Grant ERC-CoG-2015 No. 681285 “TAME-Plasmons” is gratefully acknowledged. Computational support from C3P high-performance computing facility of the Department of Chemistry of the University of Padova and from CINECA (Iscra C project “VIBCOSTO”) is also acknowledged.

References

1. D. Brinks, F. D. Stefani, F. Kulzer, R. Hildner, T. H. Taminiau, Y. Avlasevich, K. Müllen and N. F. Van Hulst, Nature, 2010, 465, 905–908.
2. D. Brinks, R. Hildner, F. D. Stefani and N. F. van Hulst, Faraday Discuss., 2011, 153, 51–60.
3. M. Liebel, C. Toninelli and N. F. van Hulst, Nat. Photonics, 2018, 12, 45–49.
4. D. Brinks, R. Hildner, E. M. H. P. van Dijk, F. D. Stefani, J. B. Nieder, J. Hernando and N. F. van Hulst, Chem. Soc. Rev., 2014, 43, 2476–2491.
5. R. Hildner, D. Brinks, J. B. Nieder, R. J. Cogdell and N. F. van Hulst, Science, 2013, 340, 1448–1451.
6. A. De Sio, F. Troiani, M. Maiuri, J. Réhault, E. Sommer, J. Lim, S. F. Huelga, M. B. Plenio, C. A. Rozzi, G. Cerullo, E. Molinari and C. Lienau, Nat. Commun., 2016, 7, 13742.
7. W. E. Moerner and L. Kador, Phys. Rev. Lett., 1989, 62, 2535–2538.
8. J. J. Macklin, J. K. Trautman, T. D. Harris and L. E. Brus, Science, 1996, 272, 255–258.
9. A. H. Zewail, J. Phys. Chem. A, 2000, 104, 5660–5694.
10. C. V. Shank, Science, 1986, 233, 1276–1280.
11. C. V. Shank, Science, 1983, 219, 1027–1031.
12. M. B. Plenio, J. Almeida and S. F. Huelga, J. Chem. Phys., 2013, 139, 235102.
13. F. Schweighöfer, L. Dworak, M. Braun, M. Zastrow, J. Wahl, I. Burghardt, K. Rück Braun and J. Wachtveitl, Sci. Rep., 2015, 5, 1–6.
14. R. Croce, R. van Grondelle, H. van Amerongen and I. van Stokkum, Light Harvesting in Photosynthesis, CRC Press, Boca Rata, 2018.
15. H. Lee, Y. C. Cheng and G. R. Fleming, Science, 2007, 316, 1462–1465.
16. G. S. Engel, T. R. Calhoun, E. L. Read, T. K. Ahn, T. Mančal, Y. C. Cheng, R. E. Blankenship and G. R. Fleming, Nature, 2007, 446, 782.
17. E. Collini, Chem. Soc. Rev., 2013, 42, 4932–4947.
18. E. Collini, C. Y. Wong, K. E. Wilk, P. M. G. Curmi, P. Brumer and G. D. Scholes, Nature, 2010, 463, 644.
19. X. S. Xie and J. K. Trautman, Annu. Rev. Phys. Chem., 1998, 49, 441–480.
20. R. Hildner, D. Brinks and N. F. van Hulst, Nat. Phys., 2010, 7, 172.
21. R. Hildner, D. Brinks, F. D. Stefani and N. F. van Hulst, Phys. Chem. Chem. Phys., 2011, 13, 1888–1894.
22. W. E. Moerner and M. Orrit, Science, 1999, 283, 1670–1676.
23. M. Orrit and J. Bernard, Phys. Rev. Lett., 1990, 65, 2716–2719.
24. J. A. Sonk, M. Caricato and H. B. Schlegel, J. Phys. Chem. A, 2011, 115, 4678–4690.
25. J. A. Sonk and H. B. Schlegel, J. Phys. Chem. A, 2011, 115, 11832–11840.
26. J. A. Sonk and H. B. Schlegel, J. Phys. Chem. A, 2012, 116, 7161–7168.
27. H. P. Breuer and F. Petruccione, The Theory of Open Quantum Systems, Oxford University Press, Oxford, 2006.
28. R. Kosloff, M. A. Ratner and W. B. Davis, J. Chem. Phys., 1997, 106, 7036–7043.
29. J. C. Tremblay, T. Klamroth and P. Saalfrank, J. Chem. Phys., 2008, 129, 84302.
30. G. Hermann and J. C. Tremblay, J. Phys. Chem. C, 2015, 119, 25606–25614.
31. J. C. Tremblay, P. Krause, T. Klamroth and P. Saalfrank, Phys. Rev. A: At., Mol., Opt. Phys., 2010, 81, 63420.
32. L. Uranga-Piña and J. C. Tremblay, J. Chem. Phys., 2014, 141, 74703.
33. W. Domcke and G. Stock, Theory of Ultrafast Nonadiabatic Excited State Processes and their Spectroscopic Detection in Real Time, 1997.
34. M. J. Tao, Q. Ai, F. G. Deng and Y. C. Cheng, Sci. Rep., 2016, 6, 1–10.
35. E. Palacino González, M. F. Gelin and W. Domcke, Phys. Chem. Chem. Phys., 2017, 19, 32296–32306.
36. E. Palacino González, M. F. Gelin and W. Domcke, Phys. Chem. Chem. Phys., 2017, 19, 32307–32319.
37. Y. Shinohara, K. Yabana, Y. Kawashita, J. I. Iwata, T. Otobe and G. F. Bertsch, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 82, 155110.
38. D. G. Tempel, M. A. Watson, R. Olivares Amaya and A. Aspuru Guzik, J. Chem. Phys., 2011, 134, 74116.
39. C. A. Rozzi, F. Troiani and I. Tavernelli, J. Phys.: Condens. Matter, 2017, 30, 13002.
40. K. Mølmer, Y. Castin and J. Dalibard, J. Opt. Soc. Am. B, 1993, 10, 524–538.
41. J. Dalibard, Y. Castin and K. Mølmer, Phys. Rev. Lett., 1992, 68, 580–583.
42. N. Gisin and I. C. Percival, J. Phys. A: Math. Gen., 1992, 25, 5677.
43. E. Coccia, F. Troiani and S. Corni, J. Chem. Phys., 2018, 148, 204112.
44. E. Coccia and S. Corni, J. Chem. Phys., 2019, 151, 044703.
45. R. Biele and R. D'Agosta, J. Phys.: Condens. Matter, 2012, 24, 273201.
46. G. Stock and W. Domcke, Phys. Rev. A: At., Mol., Opt. Phys., 1992, 45, 3032–3040.
47. G. Hermann, V. Pohl and J. C. Tremblay, J. Comput. Chem., 2017, 38, 2378–2387.
48. N. G. Van Kampen, Stochastic Processes in Physics and Chemistry, North Holland, 2007.
49. H. J. Carmichael, An Open Systems Approach to Quantum Optics, Springer, Berlin, 1993.
50. H. J. Carmichael, Statistical Methods in Quantum Optics 1: Master Equations and Fokker Planck Equations, Physics and Astronomy Online Library, Springer, Berlin, 1999.
51. G. Lindblad, Commun. Math. Phys., 1976, 48, 119–130.
52. C. W. Gardiner and P. Zoller, Quantum Noise, Springer, Berlin, 2005.
53. M. B. Plenio, J. Almeida and S. F. Huelga, J. Chem. Phys., 2013, 139, 235102.
54. M. B. Plenio and P. L. Knight, Rev. Mod. Phys., 1998, 70, 101–144.
55. B. Wolfseder and W. Domcke, Chem. Phys. Lett., 1995, 235, 370–376.
56. J. C. Tremblay, T. Klamroth and P. Saalfrank, J. Chem. Phys., 2008, 129, 84302.
57. G. D. Scholes, G. R. Fleming, L. X. Chen, A. Aspuru-Guzik, A. Buchleitner, D. F. Coker, G. S. Engel, R. van Grondelle, A. Ishizaki, D. M. Jonas, J. S. Lundeen, J. K. McCusker, S. Mukamel, J. P. Ogilvie, A. Olaya-Castro, M. A. Ratner, F. C. Spano, K. B. Whaley and X. Zhu, Nature, 2017, 543, 647–656.
58. S. Mukamel, Principles of nonlinear optical spectroscopy, Oxford University Press, New York, 1995, vol. 6.
59. A. H. Zewail, Advances in Laser Chemistry, 1978.
60. R. Shelby, C. Harris and P. Cornelius, J. Chem. Phys., 1979, 70, 34–41.
61. F. Santoro, R. Improta, A. Lami, J. Bloino and V. Barone, J. Chem. Phys., 2007, 126, 84509.
62. T. E. Sharp and H. M. Rosenstock, J. Chem. Phys., 1964, 41, 3453–3463.
63. F. Santoro, A. Lami, R. Improta, J. Bloino and V. Barone, J. Chem. Phys., 2008, 128, 224311.
64. R. Improta, V. Barone and F. Santoro, J. Phys. Chem. B, 2007, 111, 14080–14082.
65. R. Improta, V. Barone and F. Santoro, Angew. Chem., Int. Ed., 2007, 46, 405–408.
66. H. Köppel, W. Domcke, L. S. Cederbaum and W. V. Niessen, J. Chem. Phys., 1978, 69, 4252–4263.
67. D. Egorova, M. F. Gelin and W. Domcke, J. Chem. Phys., 2005, 122, 134504.
68. M. F. Gelin, D. Egorova and W. Domcke, J. Chem. Phys., 2009, 131, 124505.
69. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery Jr., J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman and D. J. Fox, Gaussian 16 Revision B.01, 2016.
70. Y. Avlasevich, S. Müller, P. Erk and K. Müllen, Chem. – Eur. J., 2007, 13, 6555–6561.
71. D. Padula, M. H. Lee, K. Claridge and A. Troisi, J. Phys. Chem. B, 2017, 121, 10026–10035.
72. N. F. Scherer, R. J. Carlson, A. Matro, M. Du, A. J. Ruggiero, V. Romero Rochin, J. A. Cina, G. R. Fleming and S. A. Rice, J. Chem. Phys., 1991, 95, 1487–1511.
73. S. Mukamel, A. Piryatinski and V. Chernyak, Acc. Chem. Res., 1999, 32, 145–154.
74. S. Mukamel, Annu. Rev. Phys. Chem., 2000, 51, 691–729.

---

Footnote

† Electronic supplementary information (ESI) available: Comparison between full and reduced models of DNQDI (Fig. S1). Comparison between SSE and Lindblad master equation (Fig. S2). Results of electronic dynamics of DNQDI with 5 excited states (Fig. S3). Excitation energies at TDDFT/B3LYP level of theory (Table S1). Analysis of high-energy processes (Table S2). NM1 and NM2 normal modes (Fig. S4 and S5). Results of simulations performed on DNQDI with a different value of detuning with respect to the one employed in the main text (Fig. S6–S8). Time profile of the pulse employed in Section 4.3.4 (Fig. S9). Results for the three level system in the nonlinear regime (Fig. S10). Simulation of ensemble effects (Fig. S11). See DOI: 10.1039/d0cp02557g

---