Probing quantum coherence in ultrafast molecular processes: an ab initio approach to open quantum systems
Abstract
Revealing possible longliving coherence in ultrafast processes allows detecting genuine quantum mechanical effects in molecules. To investigate such effects from a quantum chemistry perspective, we have developed a method for simulating the time evolution of molecular systems, based on ab initio calculations that includes relaxation and environmentinduced dephasing of the molecular wave function, whose rates are external parameters. The proposed approach combines a quantum chemistry description of the molecular target with a realtime propagation scheme within the timedependent stochastic Schrödinger equation. Moreover, it allows a quantitative characterization of the state and dynamics coherence, through the norm of coherence and the linear entropy, respectively. To test the approach, we have simulated femtosecond pulseshaping ultrafast spectroscopy of terrylenediimide, a well studied fluorophore in singlemolecule spectroscopy. Our approach is able to reproduce the experimental findings [R. Hildner et al.,Nature Phys., 7, 172 (2011)], confirming the usefulness of the approach and the correctness of the implementation.
University of Padova] Department of Chemical Sciences, University of Padova, via Marzolo 1, Padova, Italy \alsoaffiliationS3 Center, CNR Institute of Nanoscience, via Campi 213/A, Modena, Italy CNR]S3 Center, CNR Institute of Nanoscience, via Campi 213/A, Modena, Italy University of Padova] Department of Chemical Sciences, University of Padova, via Marzolo 1, Padova, Italy \alsoaffiliationS3 Center, CNR Institute of Nanoscience, via Campi 213/A, Modena, Italy \SectionNumbersOnQuantum coherence in molecules] Probing quantum coherence in ultrafast molecular processes: an ab initio approach to open quantum systems
1 Introduction
Ultrafast spectroscopy is a powerful tool to investigate, control and manipulate quantum coherence in molecules and complex systems ^{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20}. Detection of electronic and vibrational coherence in biological systems, as lightharvesting complexes involved in photosynthesis, and interpretation of the experimental evidences, is a matter of stimulating and open debate ^{21, 22, 23, 19, 24}. Recently, specific ultrafast spectroscopy techniques could probe a single molecule ^{17, 18, 19, 20}. To understand the outcomes of such experiments, theoretical and computational approaches are required, able to include all the important features of the simulated system. One such feature is certainly the ubiquitous coupling with the environment, that qualifies the probed system as open. In principle, every system must be considered as open, namely interacting with a surrounding, larger environment ^{25, 26, 27, 28, 29}.
Abstraction process leading to closed and isolated systems can indeed represent a crude approximation of the intrinsic features of the microscopic world.
For this reason, coupling between quantum systems and an external environment is essential to give a complete description of physical phenomena ^{30, 31}, also in the ionization regime ^{32}.
Timedependent Schrödinger equation (or, equivalently, von Neumann equation for the time evolution of the density matrix) describes a coherent dynamics, i.e. a dynamics with a well defined phase relation between the eigenstates of the system.
In principle, one could extend the boundary of the system and include the environment in a larger system, but this is typically impractical computationally
due to the enormous number of degrees of freedom involved in describing the environment ^{33}. Moreover, one usually is not interested in the microscopic description of the environment, which in most cases can be regarded as an external bath.
A feasible way to treat open systems is to reduce the number of degrees of freedom, by tracing out those of the bath, and defining the socalled reduced density matrix ^{25, 26, 27, 28, 29, 34, 35, 36, 37}.
Assuming that bath relaxation timescales are much faster than those of the system ^{33}, i.e. the bath is seen unchanged from the system standpoint, the time evolution of the molecular wave function is not affected by memory effects. This corresponds to the Markovian limit, an approximation used in this work^{33, 38, 39}.
Supposing a weak coupling between system and bath, and that bath degrees of freedom are in thermal equilibrium at any time, one obtains the Lindblad master equation for (in the Markovian limit)^{34, 40}. An alternative and less explored approach to the same problem is given by the Markovian stochastic Schrödinger equation (SSE). Within SSE, one directly follows the time evolution of the system wave function, , in presence of dissipation and fluctuation effects induced by the environment. SSE approach overcomes lack of microscopic knowledge about the environment by including stochastic terms in any single realization of the wave function evolution. In short, SSE time propagation of the system wave function has been seen to be fully equivalent to solve Lindblad equation for , in the limit of infinite number of quantum trajectories ^{33}. The main computational advantage of using SSE lies in the fact that the system wave function only depends linearly on the number of states of the system , while shows a quadratic dependence. On the other hand, although the single SSE realization is characterized by a linear dependence on , a large number of trajectories has to be produced: this issue can be tackled computationally by exploiting the inherent parallel nature of the procedure.
Several theoretical and computational protocols for the numerical propagation of SSE have been defined over the years ^{41, 42, 43, 44, 45, 46, 47, 48}. Details on our choice will be given below.
It is worth to mention that wellestablished approaches are present in literature, that include nonMarkovian effects, and make less restrictive assumptions about the quantum nature of the dynamics ^{49, 50, 51, 52, 53, 54}. The hierarchical equation of motion method^{49, 52, 53, 54}, for instance, is based on a hierarchy of auxiliary density matrices to account for nonMarkovian dynamics.
In the context of pigmentprotein complexes approaches based on theory of nonMarkovian open quantum systems has been successfully combined with QM/MM techniques for the study of quantum effects in photosynthesis ^{55, 56, 57, 58, 59, 60}.
In this work we present a computational approach based on an ab initio description of fluorophores, that includes the effect of dephasing and relaxation via the Markovian SSE. In particular, we show how this method can be applied to simulate and interpret ultrafast spectroscopy measurements. Coupling SSE to an ab initio representation of the electronic Hamiltonian of the molecular target is the key ingredient of the approach proposed here, allowing us to define dephasing and relaxation effects in the Hilbert space defined by the specific quantumchemistry method adopted.
A step forward in the definition of a computational protocol based on SSE is also the application of global and rigorous quantifiers to analyse coherence in (bio)chemical and physical systems. A growing interest on the development of a systematic theory of quantum coherence as a physical resource has recently arisen ^{61}. By reconstructing the density matrix from the ensemble of SSE trajectories, we can investigate the time evolution of quantum coherence by means of well established quantifiers, as norm ^{62} and linear entropy ^{63}. A quantitative analysis of quantum coherence is therefore possible, starting from an ab initio description of the molecular target. Enabling such kind of analysis, common for model Hamiltonians but rather original in an ab initio framework ^{64}, is an important result of the present work.
We have included in the realtime model developed recently ^{65}
the dephasing due to the environment surrounding the molecule and the relaxation, i.e., the spontaneous decay from the excited states to the ground state. The latter may effectively simulate nonradiative decay due, e.g., to internal conversion or radiative decay. Timedependent is expanded into the set of timeindependent eigenstates of the fluorophore, obtained in this work at Configuration Interaction with singly excited configurations (CIS) ^{66, 67, 68}, with a perturbative correction for energies involving doubly excited configurations (CIS(D)).
Only few alternative examples of coupling an ab initio description with a treatment of dephasing and relaxation are present in literature. Most of them use ab initio results as input for masterequation approaches ^{69, 70}. Stochastic approaches have been used in the context of TDDFT, at various level of complexity ^{71, 72, 73, 45}. The proposed approach is alternative to both, avoiding the possible artifacts of TDDFT beyond the linear regime ^{74, 75, 76, 77} and defining a seamless integration of ab initio and openquantum systems approaches.
As a realistic test case, we have considered pulseshaping spectroscopy on the terrylenediimide (TDI) fluorophore ^{17} (Figure 1). The authors of the original work reported how to control and manipulate electronic coherence in single molecules, by studying the interplay between variations of coherence and emission of TDI. The choice of TDI to validate our approach gives us the opportunity to compare our results with well established experimental findings.
The paper is organized as follows. In Section 2, first SSE is briefly reviewed, then we discuss how relaxation and pure dephasing channels ^{78} are introduced in the chosen quantum chemistry framework (CIS) using the version of quantum jump algorithm proposed in Refs. ^{41, 42}, and finally we provide the definition of the quantum coherence quantifiers used in this work. The numerical results are shown and discussed in Section 3, while conclusions and perspectives are collected in Section 4.
2 Theory
2.1 Stochastic Schrödinger equation
As reported in the Introduction, the fluorophore interacts with an environment causing dephasing in the wave function. In the theory of open quantum systems ^{25}, the total Hamiltonian is defined as the sum of the Hamiltonian of the system (S), of the bath (or environment, B) and their mutual interaction:
(1)  
(2) 
with and being the identity in the system and bath Hilbert space, respectively. In the present study, the system is given by the fluorophore. The strength factor modulates the interaction between the system and the bath. The interaction term has a bilinear form, characterized by the two sets of operators and , operating on the Hilbert space of the system and of the bath, respectively. Operator describes the effect of the bath on the system, while the operator describes how the bath is affected by the presence of the system. The sum in runs over the number of interaction channels between the system and the bath. Since our final goal is to simulate optical processes, in the timedependent interaction term with the external (classical) electromagnetic field is also included
(3) 
where is the Hamiltonian of the isolated system and is the system dipole interacting with the external electric field .
In this work, only the Markovian limit will be explicitly taken into account, corresponding to a delta approximation of the time autocorrelation function of the bath, on the timescale of the system response ^{33}. In other words, the loss of information from the system is irreversible ^{79}.
As pointed out in the Introduction, since the resolution of the full problem defined in Eq. 1 is impractical for realistic cases, the total number of degrees of freedom is drastically reduced when those of the bath are traced out from the total (i.e. ) density matrix : the effect of the bath is treated in an effective way by introducing the reduced density matrix of the system (time dependence is made explicit):
(4) 
Under the assumption of Markovian behaviour and of weak coupling between system and (up to the second order in the interaction term ), the time evolution of the reduced density matrix follows the Lindblad master equation ^{34}:
(5)  
The first term in the rhs of Equation 5 corresponds to the von Neumann timeevolution of the density matrix of a closed system. The interaction with the bath is described by the Lindblad superoperator .
The analytical solution of Equation 5 is known only for few model systems, otherwise a numerical realtime propagation is mandatory. The main drawback of solving the Lindblad equation lies in the dimension of reduced density matrix, i.e. , which could make numerical simulations computationally demanding.
An alternative but equivalent approach is given by the stochastic Schrödinger equation (SSE) ^{33}, written in the Markovian limit as:
(6) 
Starting from the fluctuationdissipation theorem ^{80}, one can interpret the nonHermitian term as the dissipation due to the environment, whereas is the fluctuation term, modeled by a Wiener process , i.e. a white noise associated to the Markov approximation.
The main advantage in using SSE is that one directly treats the system wave function, that only depends linearly on , thus saving computational time with respect to the propagation of .
Diagonal and offdiagonal elements of the reduced density matrix , respectively populations and coherences of the states of the system at time , are obtained by averaging on the number of independent realizations of propagating SSE. Given the definition of the reduced density matrix
(7) 
where is the system wave function corresponding to th realization, and expanding into stationary eigenstates of the system
(8) 
one defines population and coherences as the following:
(9)  
(10) 
In the limit of large number of quantum trajectories, SSE reproduces the same coming out from the Lindblad master equation ^{33}.
For sake of clarity, we will use in the following the definition: , with .
2.2 CIS expansion of the wave function
SSE, which has been used in the past for model (typically twostate) systems^{81, 82}, is here coupled to a quantumchemistry description of the molecular target. In this work, in the expansion of Eq. 8 are timedependent expansion coefficients, and represents the th timeindependent CIS eigenstate of the isolated system, with eigenvalue . The CIS eigenstates are symmetryadapted linear combinations of singlyexcited Slater determinants ^{83}:
(11) 
where is the reference HF state, is the configuration obtained by the single excitation from the occupied HF orbital to the virtual HF occupied , and the coefficients and are obtained by diagonalizing the Hamiltonian of the isolated molecule in this space. Each Slater determinant is defined as the antisymmetrized product of singleelectron molecular orbitals
(12) 
expanded on Gaussian basis functions, where r is the collective electronic coordinate.
However, our model is general and not limited to a CIS expansion of the wave function.
Our computational protocol is therefore articulated in two main steps: first, a quantumchemistry calculation; second, a realtime SSE propagation with these ab initio quantities.
The matrix form of the SSE is formally given by
(13) 
where is the vector of the timedependent expansion coefficients and is the matrix representation at time of in the basis of the CIS eigenstates (.
2.3 Including relaxation decay
Relaxation refers to the decay from an electronic excited state of the fluorophore to its ground state ^{66}. It can be due to photon emission (radiative decay) or to nonradiative decay (e.g. through internal conversion). The former can be seen as an effect of the electromagnetic field seen as an environment, the latter is more molecular based, although the ladder of vibrational levels may also be seen as an environment for the electronic level.
Nonadiabatic vibronic coupling provides a nonradiative decay channel, which plays an essential role in the molecular relaxation process^{84, 85}. In the SSE framework, such process is accounted by the operator
(14) 
which induces an exponential decay of the population and quantum jumps corresponding to the collapse of the system wavefunction into the ground state . The nonradiative relaxation rate can be regarded as a phenomenological parameter or obtained by ab initio nonadiabatic simulations ^{86}. The same operator can be used to account for the radiative relaxation. Here, given the matrix element of the optical transition , the decay rate is derived through the Fermi’s Golden Rule. In the presence of both decay channels, the overall decay rate is given by the sum of the radiative and nonradiative ones.
2.4 Including pure dephasing
Dephasing acts on the decay of the offdiagonal elements of the reduced density matrix, i.e. the coherences of the system. The choice of the form of the dephasing operator is not univocal. Dephasing operators are usually given for twostate (2s) systems as proportional to the Pauling matrix
(15) 
with the associated dephasing rate.
Here, we extend this form to a generic multistate system. In fact, we define an operator for the pure dephasing that changes the sign of the element ; more specifically, it is
(16) 
where is equal to 1 if or equal to 1 otherwise.
The operator in Eq. 16 guarantees that the population of the various states remains unchanged during the propagation (and not only on average), i. e. for each , with being the initial time. Moreover, the specific definition of the dephasing operator in Eq. 16 keeps the population unchanged also in presence of a quantum jump. Offdiagonal elements of the reduced density matrix exponentially decay with a rate equal to ; this results directly comes out from the analysis of and in the Lindblad superoperator. and are introduced as phenomenological parameters, since treating dephasing at ab initio level is challenging, due to the interaction between the system and the many degrees of freedom of the environment.
Given the definitions of the operators in Eqs. 14 and 16, one can verify that the number of interaction channels coincides with .
2.5 Quantum jump algorithm
SSE propagation can be performed by means of a quantum jump algorithm ^{41, 42, 43, 44, 45}: a deterministic nonHermitian dynamics is coupled to a number of random jumps (simulating fluctuation induced by the bath), obtained with Monte Carlo techniques ^{47}. Alternatively, one can use the quantum state diffusion model ^{87, 48}, in which the propagation is performed in terms of continuous stochastic differential equations, in linear or nonlinear form ^{71}.
Continuous propagators ^{46}, as the EulerMaruyama ^{46} or the LeihmulkerMatthews ^{88, 89} methods, can be used.
In the present investigation we have used the quantum jump algorithm, as proposed in Ref. ^{42}. The Hamiltonian of the deterministic nonHermitian part of the realtime propagation is given by
(17) 
A secondorder version of the Euler algorithm is used to propagate the coefficients ^{65, 90} of the wave function expansion, leading to
(18) 
where is the finite time step used for the numerical propagation of the deterministic part of SSE.
The time evolution based on does not conserves the norm of the wave function.
At first order in , the norm of the timedependent wave function at the time can be written as
(19) 
with
(20)  
(21) 
The quantity represents the probability that a generic quantum jump occurs, while defines the probability that the quantum jump involves the specific interaction channel , of relaxation or dephasing type. Both probabilities are imposed using Monte Carlo techniques. In detail, , which therefore corresponds to the amount of lost norm in the time step from time to , is compared at each step with a random number uniformly distributed between 0 and 1:

if no quantum jump occurs, and the wave function is then normalized;

if , a quantum jump occurs, and the new function is defined as the following
(22) with probability , determined using again the same Monte Carlo technique.
Within the CIS expansion of and using the relaxation and dephasing operators defined in Eqs 14 and 16, one finds for the relaxation channel:
(23)  
(24) 
and for the dephasing one:
(25)  
(26) 
Clearly Eq. 24 corresponds to an exponential decay of populations, as anticipated before. Consequences of Eq. 26 are less obvious, but still an exponential decay of the coherences is obtained. The dephasing operator defined in Eq. 16 maintains the population unchanged for each trajectory because is equal to the identity.
2.6 Quantifying coherence
As reported in the Introduction, different quantifiers of quantum coherence have been introduced in the last years ^{61}, fulfilling a number of fundamental requirements. In the following, we refer to the socalled norm of coherence, defined as ^{62}
(27) 
which is timedependent in our simulations due to the time evolution of the reduced density matrix. One can show that varies from zero for a fully incoherent state to 1 (where is the dimension of the Hilbert space) for a maximally coherent state.
This quantity corresponds to the smallest distance, as quantified by the norm (i.e., the least absolute deviation), between the density matrix operator at a given time and that of any incoherent state, and describes the wavelike character of the state of the system.
The norm coherence has an intuitive interpretation, since it is given by the sum of the moduli of the offdiagonal elements of the reduced density matrix , which is built from the ensemble of SSE trajectories.
The second quantifier used here is the linear entropy ^{63}
(28) 
which refers to another interpretation of coherence, related to the system dynamics. Dynamics is indeed defined coherent if the time evolution at any time 0 can be expressed as a unitary transformation of the initial state. If the system is initialized in a pure state, typically the ground state, it remains in a pure state in the presence of a coherent evolution. The degree of coherence of the dynamics can thus be quantified in terms of the purity of the density operator throughout the time evolution. The trace of varies from 1, for pure states, to 1/, for a maximally mixed state. Correspondingly, the linear entropy varies from 0 to 1  1/.
3 Numerical tests: application to fs pulseshaping spectroscopy on TDI
In this section we consider fs pulseshaping spectroscopy on TDI, simulated by means of the theoretical model described above. The goal is to validate the proposed approach by reproducing the experimentally observed emission properties of TDI ^{17}. The latter is interrogated with a sequence of two pulses (Figure 1): the first pulse generates a coherent superposition of ground and first excited state, the second one (switched on with a given time delay ) probes the phase memory in TDI. In the following, we report both the values of the excitedstate population, which is proportional to the observed fluorescence signal ^{17}, and the time evolution of the offdiagonal matrix elements. The interest in the coherences is twofold: on one hand, they are responsible for the interference between the excitations produced by the two pulses; on the other hand, the coherences are related to the wavelike character of the system state, and thus its ”quantumness” as quantified, e.g., by the norm of coherence and by the linear entropy. These quantities are computed as a function of the time delay and of the phase shift between the two pulses, for several dephasing times . We have also explored the effect of the detuning , corresponding to the difference between the central frequency of the pulsed field and the energy of the molecular transition between the ground and the first excited state of TDI, which turns out to be the most relevant one in the present case.
TDI has been extensively studied from the experimental side ^{17}. Our aim is to test the approach reported above on the detection of quantum coherence of this fluorophore by reproducing the twopulse spectroscopy results.
3.1 Computational details
The geometry of the TDI molecule has been optimized at the DFT level, with the B3LYP functional and the 631G(d) basis set.
Timeindependent CIS and CIS(D) calculations on the optimized TDI structure have been carried out using a locally modified version of Gamess ^{91, 92}. A 631G(d) basis set has been employed;
10 excited states have been kept in the expansion of the timedependent wave function, corresponding to excitation energies up to 5 eV. CIS and CIS(D) excitation energies are collected in Table S1 of Supporting Information. A nonradiative decay time of 3.5 ns ^{17} and a pure dephasing time of 30, 60 and 120 fs have been chosen. The nonradiative decay time is taken from the experimental work ^{17}, and we have selected the three values of dephasing time from the experimental distribution ^{17}: 60 fs represents the maximum of the distribution, while 30 and 120 fs represent the extreme values detected in the experiment. All these values are input parameters in our model. Detuning values of 80 and 160 cm have been employed in some of the simulations.
For all the cases studied here, dynamics of 1 ps have been considered. A time step of 1.21 as has been employed in all the simulations.
In general, any pulse shape could be chosen, since pulses are encoded numerically. In this case, the two pulses are shaped with a Gaussian envelope function
(29) 
where is the maximum field amplitude, is the center of the first pulse, is the width of the Gaussian and the carrier frequency.
FWHM has been chosen to be equal to 49 fs, corresponding to an energy bandwidth of 48 cm. The most part of the calculations has been carried out with an intensity I==5x10W/cm, while for a direct comparison with the experimental results we have also used I=6.6x10 W/cm.
The wavelength is equal to 501 nm, coinciding with the CIS(D) transition from the ground to the first excited state .
Time delays of 0, 10, 30, 50, 100, 200, 400 and 600 fs have been used, together with a phase shift of 0 and .
We have used 512 quantum trajectories for each of the simulations reported in this paper: this number assured acceptable statistical errors. In Figure S1 of the Supporting Information the convergence of the SSE results with respect to the number of trajectories is reported. Quantum jumps associated to the relaxation do not occur along the 1 ps dynamics. SSE propagation with only dissipation due to the relaxation (while quantum jumps associated to the dephasing have been observed in our simulations) produces the right time evolution of the system wave function, which refers to the first excited state as an example. Since the absolute value of is much smaller than 1, the expected exponential decay is indeed recovered without intervention of quantum jumps. Indeed, at first order in and only considering dissipation from relaxation (no jumps) one obtains
(30) 
The realtime propagation of the wave function, with the addition of relaxation and dephasing through SSE, has been performed using the homemade WaveT code ^{65}.
3.2 Results
In the following, we show the simulated time evolution of the system state, and analyze in some detail both the excitedstate populations and the coherence between ground and the first excited state.
In order to verify the reliability of the proposed approach and to test its ability to provide the right physical insight,
we focus on the following aspects: i) the overall effect of dephasing: ii) the quantitative changes in the emission of TDI, as a function of the dephasing time : iii) the effect of the detuning and of the intensity; iv) quantifying coherence. Unless otherwise specified, all the SSE calculations have been performed with an intensity of 5x10 W/cm.
We start by considering the qualitative effect of dephasing on the emission signal of TDI. In the upper panel of Figure 2, the population of the first excited state (at the end of the simulation, 1 ps) is reported as a function of the time delay (in fs) between the two pulses, for and fs. The nonradiative decay time is estimated to be around 3.5 ns, which makes the effects of relaxation negligible ^{17}. In this specific case, indeed, loss of phase memory is largely dominated by the pure dephasing time (much shorter than the relaxation one), even though in principle also relaxation plays a role in the decoherence process, as one can argue from the application of to the reduced density matrix.
According to the analysis reported in Ref. ^{17}, electronic coherence in TDI is interrogated by the second pulse, and gives rise to interference features, provided that the time delay is smaller than the dephasing time. In fact, the phase coherence is gradually erased by the interaction with the environment. As the values of become larger than , the population of the excited state tends to become independent of the phase shift applied to the second pulse. In Figure 2, we report the population of the first excited state, for and . At zero delay, these two values of the phase shift give rise to constructive and destructive interference, respectively, resulting in larger excited state occupation for .
In Figure 2 we have also reported results obtained by solving the master equation for the reduced density matrix (DM), using the same ab initio quantities (energies, transition dipole moments) of the SSE calculations. A quantitative agreement between SSE and DM data is found within the statistical error, thus validating our approach and its implementation.
Even though an inversion in the values of the populations is seen for a delay time of 400 fs, they are identical within the error bars. To clarify this aspect, we repeated the calculation by doubling the number of trajectories (from 512 to 1024): the gap between the population for the two cases ( and ) converges to zero, which is also the DM finding. (see Figure 1 in the Supporting Information).
We note that, if the molecule is considered isolated, i.e. if we switch off dephasing in the calculations, the excitedstate population is independent of the time delay: its value corresponds to twice the value obtained for two excitations summed incoherently for , while excitedstate population vanishes due to destructive interference in the case .
These interference effects are related to the coherence between and , whose time evolution is reported in the lower panel of Figure 2. The shaded region corresponds to the uncertainty produced by the average over a finite number of SSE quantum trajectories. The build up of the coherence is clearly seen at short times, namely for the first 100 fs, and is due to the pump pulse, generating a superposition of ground and excited states. At longer times, the dephasing determines a rather rapid suppression of the coherence, which eventually goes to zero (within the error bar). As explained in Section 2, the decay is exponential with a rate equal to , since the contribution to decoherence due to the relaxation is negligible here.
In order to get a more quantitative insight into the emission properties of the fluorophore, we have repeated the calculations with different values of the dephasing time. Figure 3 collects time evolutions of the excitedstate populations for fs (panel A), fs (panel B, same data of Figure 2) and fs (panel C). Using three values of allowed us to test our approach more effectively in different dephasing regimes. The value of corresponding to the merge of the curves with opposite phases increases with the dephasing time. In fact, if is much larger than , the excitations induced by the two pulses sum up incoherently, and the final population of the excited state approximately coincides with twice that induced by each pulse. The excited state population for fs and is seen to increase with the dephasing time, approaching the value obtained in the absence of dephasing (upper panel of Figure 2) in the limit where is much larger than the duration of the laser pulse. This reflects the effect of dephasing already during the excitation of the system by means of a single laser pulse.
In panel D of Figure 3 we report the experimental emission for a single TDI molecule ^{17} as a function of the time delay and for and , together with SSE results obtained with a low (I=5x10 W/cm) and a high intensity (I=6.6x10 W/cm), within the range of values used in the experiments. SSE data have been scaled to match the value at 600 fs in the two cases and to superimpose the experimental profile for phaseindependent values, since the proportionality factor between the (experimental) fluorescence count and the (computed) value of the excitedstate population at 1 ps is unknown. Regarding the curves, SSE points at high intensity are smaller for short delay times because of the occurrence of Rabi oscillations ^{17}. For both intensities, a dephasing time of 60 fs has been used.
Comparison with DM results for = 30 and 120 fs is shown in Figure S2 and S3 of the Supporting Information. Furthermore, Figures S4, S5 and S6 in the Supporting Information collect the SSE and DM time evolution of the population for and 120 fs, = 100 and 600 fs, and = 0, indicating that SSE and DM profiles coincide within the error bar. Results in panels A, B and C of Figure 3 have been obtained in resonance conditions.
Comparison between the SSE reported reported in panel B of Figure 3 (where we set fs) and the results in panel D of Figure 3 shows a nice qualitative agreement. However, we also note a large difference in the value of where the two curves meet: indeed, they meet at around 100130 fs in the experiment, while at around 300 fs in the simulations. SSE results in Figure 3D have been obtained with a detuning cm.
Indeed, in order to reproduce quantitatively the measured fluorescence, we have to take into account that the molecule in the experiment is not excited at perfect resonance, and thus we have to include a finite detuning in the simulations. Based on the information in Ref. ^{17}, reasonable values of range from few cm to around 500 cm. The effect of detuning is reported in Figure 4, with and 160 cm; populations obtained resonantly are also reported for comparison. The presence of detuning induces oscillations of the final excited state populations. Besides, the larger the detuning, the earlier a crossing between the populations corresponding to and is observed.
Using cm, we can reproduce the observed crossing position between the excitedstate populations at around 100150 fs (Figure 3D). The value of 80 cm has also been extracted from the fitting procedure in Ref. ^{17}. Effect of the detuning has been also investigated at DM level (Figures S7 and S8 of the Supporting Information): comparison between SSE and DM profiles shows a good agreement in reproducing oscillations and the detuningdependent position of the first crossing between = 0 and populations.
For TDI, in the excitation regime studied in the present work, multistate effects do not play a significant role, as verified by comparing the present results with a twostate calculation for = 60 fs in resonant conditions (Figure S9 and S10 in the Supporting Information): for delay times of 100, 200 and 400 fs populations coincide within the error.
Including relaxation between excited states, i.e. going beyond the relaxation operator defined in Eq. 14, does not significantly change the emission pattern in this case, as reported in Figure S11 of the Supporting Information.
The dependence of the excitedstate population on the waiting time and on the phase reflects the time evolution of the coherences, especially that occurring between the two laser pulses. However, the simulation of the system dynamics allows us a more general and direct characterization of the overall quantum coherence of the system state. A quantitative analysis of such coherence is reported in Figure 5, where the norm of coherence (Eq. 27) is shown as a function of time, for different values of the dephasing time and of the detuning (the time delay is fs, respectively in the upper and lower panels, and ). In the absence of dephasing (green curves) is increased by an equal amount by each laser pulse, and its final value is independent on the time delay . Coherence decreases with decreasing dephasing times . In particular, values of that are comparable to, or smaller than the pulse duration result in a reduction of the coherence generated by each pulse. Besides, if the phase memory is longer than the time delay (upper panel), the second laser pulse is followed by a maximum in that is higher than that produced by the first pulse. This corresponds to the occurrence of constructive interference in the system excitation, and thus in the observed fluorescence. If instead (lower panel), the coherence generated by the first pulse is completely suppressed by dephasing before the arrival of the second pulse. As a result, the amount of coherence generated by the two (identical) pulses coincides, up to fluctuations related to the averaging over the different trajectories, and no interference shows up in the final occupation of the excited state.
We finally note that the close correspondence that we have established between the dependences on and of the fluorescence on the one hand and of the state coherence on the other results from the fact that both are essentially related to the same, single excited state.
In general, interference features in the molecule emission result specifically from the coherence between the ground state and the excited states that are involved in the radiative recombination, while accounts for the coherence between any two eigenstates.
Time evolution of the linear entropy for the same cases as for the norm of coherence is reported in Figure 6a for a delay time of 100 fs, and in Figure 6b for a delay time of 400 fs. Phase shift is set to zero in both cases.
As mentioned above, linear entropy is a quantitative measure of the purity of the state. If the system is intialized in the ground state (or in any pure state) and in the absence of dephasing, remains a pure state throughout its time evolution, and is always zero. In the presence of dephasing, the linear entropy evolves as follows. Each of the two laser pulses tends to populate the excited states and to create a linear superposition between these and the ground state. Dephasing tends to turn such a linear superposition into a statistical mixture. The asymptotic value of the linear entropy is an increasing function of the excited state population. This explains why, rather counterintuitively, the highest values of the entropy are obtained for the larger values of the dephasing time .
In fact, with a slow dephasing (large ) populating the excited states is more efficient: this contributes to the ”disorder” of the probability distribution in the final occupations, eventually increasing . At fixed dephasing time, fs, detuning makes excitation less efficient, resulting in a smaller value of the linear entropy.
Error bars for both and (not shown) are, along the dynamics, at least one order of magnitude smaller than the corresponding value.
4 Conclusions
We have presented a new approach, which combines SSE and quantum chemical description of a fluorophore to include dephasing and relaxation, thus enabling the study of ultrafast processess on an ab initio footing. The proposed model combines standard quantum chemistry treatment of the molecular target with the study of ultrafast processes occurring at the femtosecond scale, with the goal to study possible longliving coherence effects in timedependent molecular properties. Relaxation and dephasing operators have been defined in the space of CIS timeindependent eigenstates of the fluorophore. A key point of the proposed approach is the use of an ab initio formulation of the electronic problem related to the molecular target.
As a test case, we have chosen to reproduce the experimentally detected emission signal of TDI, as a function of the delay time and phase shift between two pulses ^{17}.
The interplay between electronic quantum coherence and emission signal of the TDI fluorophore has been theoretically investigated by means of a realtime propagation of the stochastic Schrödinger equation, which includes effects from a surrounding environment, such as dephasing, and relaxation. Quantitative analysis of the coherence has been also performed, using the norm of coherence and linear entropy.
We have analyzed effects of the dephasing time, the detuning and of considering more than two states. Coherence has been quantitatively investigated.
The approach was able to correctly reproduce the experimental behavior reported in Ref. ^{17}.
Perspectives for the future work will move along three main research lines: i) the investigation of vibrational coherence by including the vibronic structure in our approach; ii) inclusion of the effects of solvent and/or nanostructures; iii) extending the method to a nonMarkovian formulation of SSE within the quantum chemistry framework.
i)
Only by including vibrational structure in our model, i.e. by going beyond the pure electronic transition, we could investigate both electronic and vibrational (i.e, vibronic) coherence. For example, recently a densitymatrix approach with a single vibrational mode has been applied to theoretically reproduce the vibrational modulations of emission intensity of the DNQDI fluorophore ^{93, 94}.
ii) The proposed model already includes the option to treat the time evolution of molecular properties embedded in a solvent and/or in presence of a nanostructure, using a polarizable continuum model ^{95, 96, 97, 98, 99, 100, 101}.
iii)
Our aim is to extend our approach to the nonMarkovian SSE ^{102, 103, 104, 54, 105, 106, 107, 108, 81, 109, 110}, in the form given in Ref. ^{102}, using the polarizable continuum model, to set the proper environment response functions and fluctuations.
In conclusion, the presented approach represents the first promising step of a long term research line that aims at integrating all the aspects of the time evolution of molecular systems probed by ultrafast spectroscopy in a complex environment, into a ab initio based simulation framework.
Funding from the ERC under the grant ERCCoG2015 No. 681285 ”TAMEPlasmons” is gratefully acknowledged. Computational support from HPC Lab of University of Modena and Reggio Emilia and from CINECA (Iscra C project ”QNANO”) is also acknowledged. Marta Rosa is gratefully acknowledged for some tests on timepropagation calculations.
Supporting Information
CIS and CIS(D) excitation energies in Table S1; comparison between results with 512 and 1024 SSE trajectories and DM data for the case reported in Figure 2 (Figure S1); comparison between SSE and DM results (Figures S2S8); comparison between multi and twostate calculations for = 60 fs (Figures S9 and S10); comparison between SSE results using only the relaxation to the ground state (Eq. 14) or also relaxation between excited states (Figure S11).
References
 Engel, G. S.; Calhoun, T. R.; Read, E. L.; Ahn, T.K.; Mancal, T.; Cheng, Y.C.; Blankenship, R. E.; Fleming, G. R. Nature 2007, 446, 782.
 Lee, H.; Cheng, Y. C.; Fleming, G. R. Science 2007, 316, 1462.
 Ishizaki, A.; Fleming, G. R. Proc. Natl. Acad. Sci. USA 2009, 106, 17255.
 Sarovar, M.; Ishizaki, A.; Fleming, G. R.; Whaley, K. B. Nature Phys. 2010, 6, 462.
 Ishizaki, A.; Fleming, G. R. Annu. Rev. Condens. Matter Phys. 2012, 3, 333.
 Cerullo, G.; De Silvestri, S. Rev. Sci. Instrum. 2003, 74, 1.
 Polli, D.; Altoé, P.; Weingart, O.; Spillane, K. M.; Manzoni, C.; Brida, D.; Tomasello, G.; Orlandi, G.; Kukura, P.; Mathies, R. A.; Garavelli, M.; Cerullo, G. Nature 2010, 467, 440.
 Rozzi, C. A.; Falke, S. M.; Spallanzani, N.; Rubio, A.; Molinari, E.; Brida, D.; Maiuri, M.; Cerullo, G.; Schramm, H.; Christoffers, J.; Lienau, C. Nat. Comm. 2013, 4, 1062.
 Falkeand, S. M.; Rozzi, C. A.; Brida, D.; Maiuri, M.; Amato, M.; Sommer, E.; De Sio, A.; Rubio, A.; Cerullo, G.; Molinari, E.; Lienau, C. Science 2014, 344, 1001.
 Cerullo, G.; Garavelli, M. Nat. Chem. 2017, 9, 506.
 Stoll, T.; Branchi, F.; Rehault, J.; Scotognella, F.; Tassone, F.; Kriegel, I.; Cerullo, G. J. Phys. Chem. Lett. 2017, 8, 2285.
 Mukamel, S. Annu Rev Phys Chem. 2000, 51, 691.
 Abramavicius, D.; Palmieri, B.; Voronine, D. V.; Sanda, F.; Mukamel, S. Chem. Rev. 2009, 109, 2350.
 Collini, E.; Wong, C. Y.; Wilk, K. E.; Curmi, P. M. G.; Brumer, P.; Scholes, G. D. Nature 2010, 463, 644.
 Kim, J.; Mukamel, S.; Scholes, G. D. Acc. Chem. Res. 2009, 42, 1375.
 Scholes, G. D. et al. Nature 2017, 543, 647.
 Hildner, R.; Brinks, D.; van Hulst, N. Nature Phys. 2011, 7, 172.
 Brinks, D.; Stefani, F. D.; Kulzer, F.; Hildner, R.; Taminiau, T. H.; Avlasevich, Y.; Müllen, K.; van Hulst, N. Nature 2010, 465, 905.
 Hildner, R.; Brinks, D.; Nieder, J. B.; Cogdell, R. J.; van Hulst, N. F. Science 2013, 340, 1448.
 Brinks, D.; Hildner, R.; van Dijk, E. M. H. P.; Stefani, F. D.; Nieder, J. B.; Hernando, J.; van Hulst, N. Chem. Soc. Rev. 2014, 43, 2476.
 van Grondelle, R.; Novoderezhkin, V. I. Phys. Chem. Chem. Phys. 2006, 8, 793.
 Scholes, G. D. J. Phys. Chem. Lett. 2010, 1, 2.
 Collini, E. Chem. Soc. Rev. 2013, 42, 4932.
 Chen, L.; Gelin, M. F.; Domcke, W.; Zhao, Y. J. Chem. Phys. 2015, 142, 164106.
 Breuer, H.P.; Petruccione, F. The Theory of Open Quantum Systems; Oxford University Press, Oxford, 2006.
 Gardiner, C.; Zoller, P. Quantum Noise; Springer, Berlin, 2005.
 Carmichael, H. An Open Systems Approach to Quantum Optics; Springer, Berlin, 1993.
 Carmichael, H. Statistical Methods in Quantum Optics 1: Master Equations and FokkerâPlanck Equations; Physics and Astronomy Online Library, Springer, Berlin, 1999.
 Weiss, U. Quantum Dissipative Systems; World Scientific, Singapore, 2012.
 Schlosshauer, M. Rev. Mod. Phys. 2005, 76, 1267.
 Daley, A. J. Adv. Phys. 2014, 63, 77.
 Tremblay, J.; Klinkusch, S.; Klamroth, T.; Saalfrank, P. J. Chem. Phys. 2011, 134, 044311.
 Biele, R.; D’Agosta, R. J. Phys.: Condens. Matter 2012, 24, 273201.
 Lindblad, G. Commun. Math. Phys. 1976, 48, 119.
 Redfield, A. G. Adv. Magn. Res. 1965, 1, 1.
 Egorova, D.; Gelin, M. F.; Thoss, M.; Wang, H.; Domcke, W. J. Chem. Phys. 2008, 129, 214303.
 Chen, L.; Gelin, M. F.; Chernyak, V. Y.; Domcke, W.; Zhao, Y. Faraday Discuss. 2016, 194, 61.
 Haken, H.; Reineker, P. Z. Phys. A 1972, 249, 253 .
 Haken, H.; Strobl, G. Z. Phys. A 1973, 262, 135 .
 Gorini, V.; Kossakowski, A.; Sudarshan, E. C. G. J. Math. Phys. 1976, 17, 821.
 Dalibard, J.; Castin, Y.; Mølmer, K. Phys. Rev. Lett. 1992, 68, 580.
 Mølmer, K.; Castin, Y.; Dalibard, J. J. Opt. Soc. Am. B 1993, 10, 524.
 Dum, R.; Zoller, P.; Ritsch, H. Phys. Rev. A 1992, 45, 4879.
 Plenio, M. B.; Night, P. L. Rev. Mod. Phys. 1998, 70, 101–144.
 HofmannMees, D.; Appel, H.; Ventra, M. D.; Kümmel, S. J. Phys. Chem. B 2013, 117, 14408.
 Higham, D. J. SIAM Rev. 2001, 643, 525.
 Makarov, D. E.; Metiu, H. J. Chem. Phys. 1999, 11, 10126.
 Gisin, N.; Percival, I. C. J. Phys. A: Math. Gen. 1993, 26, 2245.
 Tanimura, Y.; Kubo, R. J. Phys. Soc. Jpn. 1989, 58, 1199 .
 Makarov, D.; Makri, N. Chem. Phys. Lett. 1994, 221, 482 .
 Makri, N. J. Math. Phys. 1995, 36, 2430 .
 Yan, Y.; Xu, R. Annu. Rev. Phys. Chem. 2005, 56, 187 .
 Schröder, M.; Kleinekathöfer, U.; Schreiber, M. J. Chem. Phys. 2006, 124, 084903 .
 Roden, J.; Eisfeld, A.; Wolff, W.; Strunz, W. T. Phys. Rev. Lett. 2009, 103, 058301 .
 Damjanovic, A.; Kosztin, I.; Kleinekathöfer, U.; Schulten, K. Phys. Rev. E 2002, 65, 031919 .
 Olbrich, C.; Jansen, T. L. C.; Liebers, J.; Aghtar, M.; Strümpfer, J.; Schulten, K.; Knoester, J.; Kleinekathöfer, U. J. Phys. Chem. B 2011, 115, 8609 .
 Shim, S.; Rebentrost, P.; Valleau, S.; AspuruGuzik, A. Biophys. J. 2012, 102, 649 .
 Valleau, S.; Eisfeld, A.; ; AspuruGuzik, A. J. Chem. Phys. 2012, 137, 224103 .
 Rivera, E.; Montemayor, D.; Masia, M.; Coker, D. F. J. Phys. Chem. B 2013, 117, 5510 .
 Wang, X.; Ritschel, G.; Wüster, S.; Eisfeld, A. Phys. Chem. Chem. Phys. 2015, 17, 25629 .
 Streltsov, A.; Adesso, G.; Plenio, M. B. Rev. Mod. Phys. 2017, 89, 041003.
 Baumgratz, T.; Cramer, M.; Plenio, M. B. Phys. Rev. Lett. 2014, 113, 140401.
 Zurek, W. H.; Habib, S.; Paz, J. P. Phys. Rev. Lett. 1993, 70, 1187.
 Rozzi, C. A.; Troiani, F.; Tavernelli, I. J. Phys.: Condens. Matter 2017, 30, 013002.
 Pipolo, S.; Corni, S. J. Phys. Chem. C 2016, 120, 28774–28781.
 Tremblay, J. C.; Klamroth, T.; Saalfrank, P. J. Chem. Phys. 2008, 129, 084302.
 Hermann, G.; Tremblay, J. C. J. Phys. Chem C 2015, 119, 25606.
 Klinkusch, S.; Tremblay, J. C. J. Chem. Phys. 2016, 144, 184108.
 Zenichowski, K.; Dokic, J.; Klamroth, T.; Saalfrank, P. J. Chem. Phys. 2012, 136, 094705.
 Rivalta, I.; Nenov, A.; Cerullo, G.; Mukamel, S.; Garavelli, M. Int. J. Quant. Chem. 2014, 114, 85.
 D’Agosta, R.; Di Ventra, M. Phys. Rev. B 2008, 78, 165105.
 Appel, H.; Di Ventra, M. Chem. Phys. 2011, 391, 27.
 Tempel, D. G.; AspuruGuzik, A. Chem. Phys. 2011, 391, 130.
 Raghunathan, S.; Nest, M. J. Chem. Theory Comput. 2011, 7, 2492.
 Raghunathan, S.; Nest, M. J. Chem. Theory Comput. 2012, 8, 806.
 Raghunathan, S.; Nest, M. J. Phys. Chem. A 2012, 116, 8490.
 Habenicht, B. F.; Tani, N. P.; Provorse, M. R.; Isborn, C. M. J. Chem. Phys. 2014, 141, 184112.
 Tempel, D. G.; Watson, M. A.; OlivaresAmaya, R.; AspuruGuzik, A. J. Chem. Phys. 2011, 134, 074116 .
 Piilo, J.; Maniscalco, S.; Härkönen, K.; Suominen, K.A. Phys. Rev. Lett. 2008, 100, 180402.
 Van Kampen, N. G. Stochastic Processes in Physics and Chemistry; North Holland, 2007.
 Ritschel, G.; Suess, D.; Möbius, S.; Strunz, W. T.; Eisfeld, A. J. Chem. Phys. 2015, 142, 034115.
 Zhang, P.P.; Eisfeld, A. J. Phys. Chem. Lett. 2016, 7, 4488 .
 Dreuw, A.; HeadGordon, M. Chem. Rev. 2005, 105, 4009.
 Toniolo, A.; Persico, M. J. Comput. Chem. 2001, 22, 968.
 Santoro, F.; Lami, A.; Improta, R.; Barone, V. J. Chem. Phys. 2007, 126, 184102.
 Curchod, B. F. E.; Martinez, T. J. Chem. Rev. 2018, 118, 3305.
 Gisin, N.; Percival, I. C. J. Phys. A: Math. Gen. 1992, 25, 5677.
 Leimkuhler, B.; Matthews, C. J. Chem. Phys. 2013, 138, 174102.
 Leimkuhler, B.; Matthews, C. Appl. Math. Res. Express 2013, 1, 34–56.
 Pipolo, S.; Corni, S.; Cammi, R. J. Chem. Phys. 2017, 146, 064116.
 Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; et al., J. Comput. Chem. 1993, 14, 13347.
 Gordon, M. S.; Schmidt, M. In In Theory and Applications of Computational Chemistry: the first forty years; Dykstra, C., Frenking, G., Kim, K., Scuseria, G. E., Eds.; Elsevier: Amsterdam, 2005; p 1167.
 PalacinoGonzalez, E.; Gelin, M. F.; Domcke, W. Phys. Chem. Chem. Phys. 2017, 19, 32296.
 PalacinoGonzalez, E.; Gelin, M. F.; Domcke, W. Phys. Chem. Chem. Phys. 2017, 19, 32307.
 Tomasi, J.; Mennucci, B.; Cammi, R. Chem. Rev. 2005, 105, 2999.
 Tomasi, J.; Cammi, R.; Mennucci, B.; Cappelli, C.; Corni, S. Phys. Chem. Chem. Phys. 2002, 4, 5697.
 Caricato, M.; Mennucci, B.; Tomasi, J.; Ingrosso, F.; Cammi, R.; Corni, S.; Scalmani, G. J. Chem. Phys. 2006, 24, 124520.
 Corni, S.; Tomasi, J. J. Chem. Phys. 2001, 114, 3739.
 Andreussi, O.; Corni, S.; Mennucci, B.; Tomasi, J. J. Chem. Phys. 2004, 121, 10190.
 Caricato, M.; Andreussi, O.; Corni, S. J. Phys. Chem. B 2006, 110, 16652.
 Vukovic, S.; Corni, S.; Mennucci, B. J. Phys. Chem. C 2009, 113, 121.
 Gaspard, P.; Nagaoka, M. J. Chem. Phys. 1999, 111, 5676.
 Breuer, H.P.; Kappler, B.; Petruccione, F. Phys. Rev. A 1999, 59, 1633.
 Stockburger, J.; Grabert, H. Phys. Rev. Lett. 2002, 88, 170407.
 Roden, J.; Strunz, W. T.; Eisfeld, A. J. Chem. Phys. 2011, 134, 034902.
 Roden, J.; Strunz, W. T.; Whaley, K. B.; Eisfeld, A. J. Chem. Phys. 2012, 137, 204110.
 Suess, D.; Eisfeld, A.; Strunz, W. T. Phys. Rev. Lett. 2014, 113, 150403.
 Biele, R.; Timm, C.; D’Agosta, R. J. Phys.: Condens. Matter 2014, 26, 395303.
 de Vega, I.; Alonso, D. Rev. Mod. Phys. 2017, 89, 015001.
 CohenTannoudji, C.; DupontRoc, J.; Grynberg, G. AtomPhoton Interactions: Basic Process and Applications; WILEYVCH, 2004.