Coherent ( photon ) vs incoherent ( current ) detection of multidimensional optical signals from single molecules in open junctions

The nonlinear optical response of a current-carrying single molecule coupled to two metal leads and driven by a sequence of impulsive optical pulses with controllable phases and time delays is calculated. Coherent (stimulated, heterodyne) detection of photons and incoherent detection of the optically induced current are compared. Using a diagrammatic Liouville space superoperator formalism, the signals are recast in terms of molecular correlation functions which are then expanded in the many-body molecular states. Two dimensional signals in benzene-1,4-dithiol molecule show cross peaks involving charged states. The correlation between optical and charge current signal is also observed. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4919955]


I. INTRODUCTION
Ultrafast picosecond to femtosecond processes determine the elementary photophysical and photochemical properties of molecules, for example, the (200 fs) isomerization of rhodopsins 1,2 which is the primary step in visual signalling, vibrational relaxation, charge transfer in molecules, and coherent energy transfer in chromophore complexes. 3,4 It is of interest to understand the ultrafast dynamics at the single molecule level since ensemble averaging often obscures the microscopic dynamics.
Time-domain multidimensional spectroscopy 5-8 uses a sequence of well separated impulsive laser pulses and records the signal versus the time delays between these pulses. In separate developments, it is now possible to prepare a single molecule in contact with metal leads (single molecular junctions) and study its conduction properties. Molecular junctions are characterized by their response to electric as well as optical radiation fields. [9][10][11] Single molecule spectroscopy is markedly different from bulk samples. Due to small size of the molecule (comparable to the wavelength of the electromagnetic field), the wavevector matching condition cannot be used to separate different contributions to coherent signals. This can be done however by phase cycling protocols. [12][13][14][15][16][17] Understanding optical signals from current-carrying single molecular systems [18][19][20][21] have drawn much fundamental interest with potential applications to sensors and optoelectronic devices. Time-resolved optical measurements in open junctions could provide information about neutral and various charged states of the molecule and their relaxation due to interactions with the surrounding. 17 Here, we use a diagrammatic Liouville space approach to calculate time-resolved signals from a molecular junction. Both photon and electronic (current) detections are compared. After the junction reaches a steady state, characterized by a constant electron flux, it is driven by a sequence of temporally well-separated optical pulses. The photon and electron signals are then expressed in terms of molecular correlation functions and expanded in molecular eigenstates, obtained from standard quantum chemistry calculations. We apply this formulation to compute signals from a junction made of Benzene-1,4-dithiol molecule.
We first introduce the model in Sec. II and calculate coherent (heterodyne) detected nonlinear stimulated signal in response to three temporally well separated pulses. In Sec. III, we give numerical results for coherent signal for benzene-1,4dithiol molecular junction setup. In Sec. IV, we simulate the photo-induced transient incoherent charge current signal and compare with the optical stimulated signal. We conclude in Sec. V.

II. COHERENT HETERODYNE DETECTED SIGNALS IN MOLECULAR JUNCTION
We consider a single molecule (M) coupled to two metallic electrodes (A and B) held at different chemical potentials. The system is described by the Hamiltonian where H 0 = H M + H A + H B + H F . H M represents the Hamiltonian for isolated molecule with electronic and vibrational degrees of freedom. The non-interacting metallic leads A and B and the optical field Hamiltonian H F are given as Here, ω s is the frequency for radiation field mode, ϵ ν represents the ν-th electron energy state of the leads. c † (c) and a † (a) denote single particle Fermi and Bose creation (annihilation) operators, respectively. The last two terms in (1) represent the molecule-lead tunnelling, FIG. 1. Laser pulse sequence for time domain experiment. The system first interacts with the pulse E 1 atτ 1 , then E 2 atτ 2 , and finally with E 3 atτ 3 . The fourth pulse E 4 atτ 4 detects the stimulated T 3 . Each field is distinguished by its frequency ω s and phase φ s , s = 1, 2, 3, 4. and the molecule-field coupling Hamiltonian given as (in the interaction picture with respect to the field Hamiltonian H F ) where T iν is the tunnelling matrix element between the i-th orbitals of the molecule and ν-th state of the electrode. H M F (t) is written in the rotating wave approximation (RWA) where E s (t) represents the positive frequency component of the field andV (V † ) is the dipole lowering (raising) operator which brings the molecule down (up) from an excited (lower) state to a lower (excited) state. The molecule is initially prepared in a steady state by keeping the leads at different chemical potential. The steady state is characterized by a constant electron flux between the molecule and the leads. The molecule is subsequently driven by a sequence of impulsive optical pulses with controllable phases and time delays. The nonlinear response is then recorded as a function of the time delays between these pulses.
We first consider heterodyne detected stimulated signal from a junction driven by four temporally well separated pulses. The positive frequency component of the total field is written as where E s (t −τ s ) = E s (t −τ s ) a s exp(−iω s t) is the complex envelope of the s-th pulse, centred aroundτ s . Time-ordering is fully maintained: pulse E 1 comes first atτ 1 , followed E 2 at τ 2 , E 3 atτ 3 , and finally the detected pulse E 4 atτ 4 . We denote the time delays between the pulses T s =τ s+1 −τ s , s = 1, 2, 3 ( Fig. 1). The emission pattern from a single molecule is significantly different from the directional beam obtained in four-wave mixing in bulk samples. Since the molecule is small compared to the optical wavelength, the emission is isotropic and depends on all possible combination of phases ± φ 1 ± φ 2 ± φ 3 ± φ 4 . The various quantum pathways of the signal can be separated out by the phase cycling protocol which selects a given phase. 14 Here, we derive the stimulated signal for a particular phase Other phases can be calculated in a similar way. Since the molecule at steady state is exchanging electrons with the leads, several neutral, and excited charged states may have appreciable populations. The interaction with the three time ordered incoming pulses at times τ 1 , τ 2 , and τ 3 may correspond to either emission or absorption of a photon. Note that this is different from the equilibrium case where the initial interaction always corresponds to an absorption from the ground state. Each molecular-field interaction can act on the molecular density matrix either from the left (ket) or from the right (bra), which generates 2 3 = 8 ladder diagrams (quantum pathways) as shown in Fig. 2. In our previous work, 21 we considered stimulated signals generated by continuous waves, which depend on many additional Liouville space pathways as the fields are not time-ordered. However, time-ordered pulse sequences are more selective.
Expressions for the signal can be obtained in terms of superoperators using the diagrams in Fig. 2. We first briefly introduce the basic elements of Liouville space superoperator notation. With each Hilbert space operator A, we associate a left (L) and a right (R) superoperator, defined by, A L X ≡ AX, and A R X ≡ X A where X is an arbitrary Hilbert space operator. Furthermore we define linear combinations of L and FIG. 2. Ladder diagrams for stimulated signal from a molecular junction in the fourth order of optical field with phase signature φ 1 − φ 2 + φ 3 − φ 4 . The molecule is driven with three impulsive laser pulses with incoming frequencies ω 1 , ω 2 , and ω 3 and heterodyne detected with the fourth impulsive pulse with frequency ω 4 . |a⟩, |b⟩, |c⟩. . . denote the neutral as well as singly charged many body states of the molecule. Time is increasing from bottom to top. The propagation in each branch is determined by the Hamiltonian H M + H A + H B + H M L which lead to the evolution between different charged states.

R superoperators as
The stimulated signal is then written in a compact form as where . The signal in Eq. (6) contains contributions from all 8 diagrams shown in Fig. 2. We assume that the fields are in coherent states and are impulsive, i.e., E j (t) = E j δ(t), j = 1, 2, 3, 4. The derivation is given in Appendix A. The delays between the pulses T i , i = 1, 2, 3 enter through the retarded propagators in the molecule + lead space is the Heaviside step function. The expectation value in Eq. (6) is with respect to the molecule + lead density operator.
We perform partial trace over the leads' degrees of freedom and obtain the correlation function in terms of reduced dynamics of the molecule, where in the second line we perform the trace over leads A and B, and The signal is now expressed in terms of the reduced density matrix ρ M (t 0 ) for the molecule at steady state. Another effect of the partial trace is to replace the unitary evolution G(t) by the non-unitary Liouville space evolution G(t). There are several techniques to compute the reduced evolution. We adopt a quantum master equation (QME) approach (Appendix B) and write G(t) = θ(t)e Lt where L is the effective Liouvillian. The dependence on the molecule-lead interactions and the bias in the signal enters through the steady state density matrix ρ M (t 0 ) and the propagator G(t). 21 Alternatively, the multidimensional signal in Eq. (6) may be recast in the frequency domain by performing Fourier transformation with respect to any two selected delays, e.g., where Ω 1 and Ω 3 are Fourier frequency conjugates corresponding to the delay times T 1 and T 3 (often called coherence times), respectively. This gives where The superoperator correlation function appearing in the signal can be expanded in terms of the many body states of the isolated molecule. 21 For a general steady state density matrix ρ =  ab ρ ab |ab⟩⟩, the signal is given by Note that, unlike equilibrium spectroscopy, Eq. (10) depends on population and coherence. Here, ρ ab = ⟨a|ρ|b⟩ is the steady state density matrix element between the many-body states |a⟩ and |b⟩, V ab = ⟨a|V |b⟩ is the transition dipole matrix element, and G ab;cd (ω) = i⟨⟨ab|(ωI − iL) −1 |cd⟩⟩ is the nonunitary time-evolution propagator in frequency domain. Since the many body states a, b, c . . . can represent neutral or charged (cation and anion) states of the molecule, the signal can show resonances involving charge states. The matrix elements of the evolution operators ⟨⟨ab|G(ω)|cd⟩⟩ and ⟨⟨ab|G(t)|cd⟩⟩ are evaluated by solving the QME given in Appendix B. The QME approach ignores the effect of bias on the many-body energy states. Such effects are important for strongly coupled systems and require a self-consistent calculation, as was done in Refs. 22 and 23.

III. COHERENT 2D SIGNAL OF BENZENE-1,4-DITHIOL
We consider a benzene-1,4-dithiol molecule connected to two gold leads as our model system (Fig. 3). This is a widely studied system similar to the mechanical break junction of Reed et al. 24 To describe the strong coupling between the molecule and the electrodes. 25,26 We used an extended molecule that includes 3 gold atoms bonded to each sulfur. The three gold atoms represents the Au(111) surface and the sulfur atom resides in its hollow site. 27 optimized at the density functional theory (DFT) level with the B3LYP functional 29-31 by using the Gaussian 09 program. 32 The Au-Au distance is held at 2.88 Å, taken from the distance in the Au(111) surface. 28 The LANL2DZ basis set and pseudopotential 33 are set for gold, and the 6-31G* basis set for the rest atoms. In the optimized geometry, the distance between the sulfur and the gold surface is around 2.3 Å. At the optimized geometry, we calculated 10 low-lying excited states of the neutral system as well as two charged systems with −1 and +1 charges by using the NWChem package. 34 The timedependent DFT (TDDFT) method with the Tamm-Dancoff approximation (TDA) 35 is used. This method has been successfully employed in simulations of current carrying molecular junctions. [36][37][38] A larger 6-311+ G * basis set is used for the dithiol part, for better description of the charged states. For the doublet anion and cation, calculated states with large spin contaminations were filtered out (a threshold of |⟨S 2 ⟩ − 3/4| > 1.0 was used where 3/4 is the reference value for a doublet state). 39,40 This results in only 5 anion and 7 cation physical valence states. Calculated ground and valence excited state energies for neutral and anion molecule and the dipole moment values are given in Appendix C. Details for calculating generalized overlap amplitudes 41 between different many-electronic states of the molecule are given in Appendix D.
All simulations were made at room temperature T = 300 K keeping the bias difference between the two leads as 4 eV (µ A = −5 eV and µ B = −1 eV). For this bias only the neutral and anion states of the molecule are populated. The cation states are much higher in energy (> 9 eV ≫ k B T) and are ignored in numerical calculation. In this applied bias, we find six anion and eleven neutral populated states. Moreover, we consider the dipole transitions only between the ground and other excited states for the same charged species. To compute the signal, we first obtain steady state population and coherence for the neutral and charge states by numerically solving the master equation (Eq. (B3)) for ∂ ∂t | ρ ss ⟩⟩ = 0. The populations and absolute values of dominant coherences are shown in Fig. 4. We observed that the coherences between same charged states are finite, however Fock space coherence between different charged states in the steady state vanishes. The optical fields are considered to be polarized in the Y-Z plane with E x = E y = E 0 . We therefore consider only the Y and Z component of the dipole vector (V y , V z ). They are given in Table II in Appendix C.
We first examine the stimulated signal for pathway (a1) of Fig. 2. This pathway can be selected under the resonant condition by tuning the external laser frequencies ω 1 ,ω 2 ,ω 3 corresponding to a particular set of transitions |a⟩ → |r⟩, |q⟩ → |h⟩, and |g⟩ → |e⟩. For this pathway, the first interaction with the field is an absorption of a photon from the left (ket) side. Therefore, the anion and neutral ground states give rise to distinct pathways as shown in Fig. 6. Fig. 5 shows the real, imaginary, and the absolute value of the two-dimensional (2D) signal as a function of Ω 1 and Ω 3 for different delay times T 2 . The real (φ = 0) and imaginary parts (φ = π/2) of the signal are obtained by controlling the phases of the optical fields. The signal shows two strong peaks whose positions are marked along the X and Y-axis corresponding to the dipole transitions between anion states |a ′ ⟩ and |d ′ ⟩ and neutral states |a⟩ and |d⟩ (state labelling is given in Appendix C: Table I). These correspond to the dominant transition dipole matrix elements (Appendix C: Table II). The signal coming from the |d⟩ → |a⟩ transition is weak compared to the transition from |d ′ ⟩ → |a ′ ⟩ as the steady state population ρ ss aa is small compared to that of ρ ss a ′ a ′ (see Fig. 4) and ⟨a ′ |V |d ′ ⟩ > ⟨a|V |d⟩. Also since for zero delay time (T 2 = 0), the second and third interactions occur simultaneously and therefore the evolution of states during T 1 and T 3 periods are the same which gives rise to diagonal peaks. Note that this is because we are considering only most dominant transitions for neutral and anion states. In principle, even T 2 = 0 can give rise to off-diagonal peaks. For T 2 0, population is transferred between the charge and neutral states due to the interaction with the leads, which results in optical transitions of two different charge states during T 1 and T 3 . This gives rise to the off-diagonal peaks in the signal when displayed in the Fourier domain as a function of Ω 1 and Ω 3 . Fig. 6 presents different pathways involving only the dominant dipole transitions. Diagrams (a) and (b) generate diagonal peaks for all delay times T 2 , whereas the two pathways (c) and (d) are responsible for the cross peaks for T 2 0. These two pathways are not allowed for T 2 = 0. With increasing T 2 , the weak peak (panel (a) in Fig. 5) is significantly enhanced indicating population transfer from state |a ′ ⟩ to state |a⟩ which occurs in the time scale of ≈30 fs. Therefore, the change in intensity of peaks with delay time T 2 carry information about the relative magnitude of population transfer due to the interaction with the leads as well as the time-scale for population relaxation which also provides a measure for system-lead coupling strength for particular molecular states. In Fig. 7, we display the contribution to pathway (a1) starting with only coherences S (a1,coh) stim (Ω 1 ,T 2 , Ω 3 ). The signal due to initial population S (a1,pop) stim (Ω 1 ,T 2 , Ω 3 ) is shown in Appendix E (Fig. 14). The total signal is dominated by the population contribution. Therefore, the signal in Fig. 5 is nearly indistinguishable from Fig. 14. However, the signal coming due to quantum coherence in Fig. 7 shows different contributions. The appearance of a single diagonal peak for pathway (a1) indicates that the coherence between ground state and the higher excited state is much stronger for a particular charge state. In Fig. 4, we see that the absolute value of coherence ρ a ′ b ′ is much higher compared to ρ ab and therefore the state |a ′ ⟩ gets excited by absorbing a photon and generates the signal. In Fig. 8, we show the different pathways starting with the dominant coherence element ρ a ′ b ′. For T 2 = 0, diagram (a) gives rise to peak for |a ′ ⟩ − |d ′ ⟩ transition. For finite T 2 , the states evolving during T 1 are different from those of T 3 (diagrams (b) and (c)) which gives rise to cross peaks. Other coherence elements do not contribute to the pathway (a1) as only the ground states can be optically excited.
In Fig. 9, we display the total stimulated signal by incorporating all possible quantum pathways (a1-a8) shown in Fig. 2. In this case, for T 2 = 0, the signal shows three strong peaks along the diagonal (see panel (a)). The extra peak appears due to |i⟩ → |a⟩ optical transition. This peak appears to be stronger as the transition frequency is near resonance with the second and third incoming field frequencies (ω 2 = ω 3 = 2.3 eV). As before, cross peaks appear as the delay time T 2 (panel (d)) is increased. In addition, the weak diagonal peaks (|a⟩ − |d⟩ and |a⟩ − |i⟩) are enhanced which indicates that the population from anion state |a ′ ⟩ is transferred to neutral states |a⟩ and |i⟩. Note that the broadening of the resonance peaks is due to the coupling with the leads.

IV. INCOHERENT ELECTRONIC CURRENT DETECTION OF MULTIDIMENSIONAL SIGNALS
An alternative detection of optical response in a molecular junction setup is provided by the charge current. The steady state charge current and its higher order fluctuations (noise) in the presence of a fixed chemical bias have been studied extensively. These were mostly formulated using nonequilibrium Green's function [42][43][44][45][46][47][48][49][50] and the quantum master equation [51][52][53] approaches. However, much less attention has been given to the photoinduced transient charge current. This incoherent detection of coherent signals can reveal information about nonequilibrium dynamics and many body interaction effects. 54,55 In this section, we apply the diagrammatic formalism to calculate the response to charge current to optical pulses. As in Sec. II, the molecular junction at steady state is subjected to four temporally well separated impulsive optical pulses. The total charge of the system at any given time t is where N =  s c † s c s is the charge operator. | ρ T (t)⟩⟩ is the total density matrix in the molecule + lead + field space. We compute the signal perturbatively in the incoming pulses by transforming to the interaction picture with respect to the molecule + lead HamiltonianH = H 0 + H M L . We then write where Expanding the time-ordered exponential to fourth order in the field-matter interaction, assuming the field involves as impulsive and selecting a particular phase component φ = φ 1 − φ 2 + φ 3 − φ 4 , we obtain We first trace over the lead degrees of freedom and then take the time derivative with respect to t to obtainQ(t) = I A (t) + I B (t) where I A and I B are instantaneous charge currents from the left and right leads, T 1 =τ 2 −τ 1 , T 2 =τ 3 −τ 2 , T 3 =τ 4 −τ 3 , andt = t −τ 4 are the delay times between the pulses. L X is the contribution to the full Liouvillian L from lead X such that L = L A + L B . L X evolves the molecular density matrix due to the coupling with lead X. | ρ M ⟩⟩ is the steady state density matrix for the molecule. We perform Fourier transformation with respect to the two time delays T 1 and T 3 to get For better comparison with the stimulated optical signal (Eq. (9)), we integrate the t variable and obtain the total charge flow from the molecule to lead X, There are 2 4 = 16 (2 V − and 2 V † − ) pathways that contribute to the signal and can be written in terms of the many-body states, as was done in Sec. II. Note that Eq. (17) involves four Green's functions as opposed to three in Eq. (9). However, the density matrix is the same up to the third evolution in both cases.
In Fig. 10, we show 2D plots of the total charge flow between the molecule and the lead A, in presence of optical perturbations. We plot the 2D signal as a function of Ω 1 and Ω 3 for fixed T 2 , and the delay between final pulse and observation time,t. We compute the total signal (Eq. (17)) by summing over all possible pathways as we did for stimulated signal. We first fix the delay T 2 = 0 and plot the signal for varioust. Fort = 0 fs (see panel (a)), the signal shows only one strong diagonal peak (Ω 1 = Ω 3 ) corresponding to the dominant dipole transition for anion state |a ′ ⟩ − |d ′ ⟩. This superposition state evolving during T 3 indicates that after fourth optical interaction the molecule resides in an anionic population state and finally the operator L A transfer an electron from the molecule to the lead A and the molecule jumps from anion to neutral state and contribute to the charge current. for the optical transitions between anionic states |a ′ ⟩ − |d ′ ⟩ and neutral states |a⟩ − |d⟩ and |a⟩ − |i⟩, similar to what we observed for the photon detection. In this case, because of the presence of extra time evolution propagator G(t), the molecule after the fourth optical interaction could reside in neutral or anionic population state, and G(t) further evolves this state to the anionic state and finally L A convert the molecule from anion to a neutral state. Therefore, the molecule during T 1 and T 3 evolutions could be in coherence states involving both neutral or charged states and thus we see three dominant transitions along the diagonal, one transition due to anionic and two transitions due to neutral states. Also note that with increasingt, the intensity of the peak for the transitions |a ′ ⟩ − |d ′ ⟩ becomes weak whereas peak due to the transition |a⟩ − |i⟩ is enhanced. This clearly indicates that the extra evolution G(t) transfer population from neutral states to anionic states and contribute more to the charge current signal. As done for coherent stimulated signal, in Appendix E (Figs. 15 and  16), we display the contribution to the charge current signal due to initial populations Q pop A (t; Ω 1 ,T 2 , Ω 3 ) and coherences Q coh A (t; Ω 1 ,T 2 , Ω 3 ).
In Fig. 11, we show the charge flow signal Q A for fixed t = 10 fs with different values of T 2 . For finite T 2 , the crucial observation is the appearance of cross peaks, in addition to diagonal peaks, which was significantly weak for T 2 = 0 fs, as shown in Fig. 10. Moreover similar to the photon detected signal, we also observe here that with increasing T 2 these cross peaks become stronger implying population transfer.

A. Comparison between coherent stimulated signal and incoherent charge current signal
In Fig. 12, we compare the absolute values for the total stimulated signal and the integrated charge current signal for both the leads. For better comparison between the two signals, we considert = 0 limit for the charge current which imply that the observation time t coincides with the arrival time of the fourth pulse T 2 . As mentioned earlier, the diagonal peaks (Ω 1 = Ω 3 ) in this case imply that, during T 1 and T 3 evolutions, the molecule is present in a superposition of anionic electronic states |a ′ ⟩ − |d ′ ⟩ which reduces to a population state (|a ′ ⟩ − |a ′ ⟩, |d ′ ⟩ − |d ′ ⟩) after the fourth optical interaction. Finally, the operator L A (Eq. (17)) transfers of an electron from the molecule to the lead A and the molecule jumps to a neutral state. For finite T 2 , because of the intermediate propagator G(T 2 ), the coherence state |a ′ ⟩ − |d ′ ⟩ during T 3 can be reached via various additional quantum pathways starting initially from both the anion and neutral states. Therefore, for T 2 0, the additional cross peaks show up in the signal for different Ω 1 corresponding to the transitions within the neutral states |a⟩ − |d⟩ and |a⟩ − |i⟩. Similar conclusion can be drawn for the integrated charge current signal from lead B (Q B ) which shows two diagonal peaks for resonances within neutral states |a⟩ − |d⟩ and |a⟩ − |i⟩ and one off-diagonal peak for resonance within anion states |a ′ ⟩ − |d ′ ⟩. We note that the sum of two signals, Q A + Q B is identical to the stimulated optical signal. This holds true for large range of parameters. It therefore shows an interesting correlation between the optical and charge current signal in molecular junction.

V. CONCLUSIONS
We have employed a Liouville space diagrammatic approach to calculate nonlinear coherent stimulated signal and incoherent charge current response from a molecular junction driven by impulsive optical pulses. The signals are calculated to fourth order in the molecule-field interaction. They are expressed in terms of superoperator correlation functions which are expanded in terms of many body states of the molecule. The nonunitary evolution of the system due to coupling to the leads is accounted for by the QME approach. We apply this formulation to study a molecular junction made of benzene-1,4-dithiol. By controlling the delay times between the pulses both photon and charge current signals reveal information about the timescale for population transfer, dissipation, (relaxation) induced by the charge transfer due to leads which can also be used to estimate the lead-molecule coupling strength. We also find that the measurement of time integrated left and right lead charge currents contain information about different charged states of the molecule and the sum of this two incoherent signal is identical to the coherent optical signal which further indicates a strong correlation between optical and charge response of the molecular junction.

APPENDIX A: DERIVATION OF COHERENT HETERODYNE DETECTED STIMULATED SIGNAL
Here, we give the derivation of Eq. (6). In the RWA, the coherent heterodyne detected stimulated signal 56,57 is given as where E 4 refers to the detected field mode which is assumed to be in a coherent state and is centered aroundτ 4 . ρ T (t) is the total density matrix in the molecule + leads space and is evolving with respect to the full Hamiltonian H(t) given in Eq. (1). The signal can be computed perturbatively in the incoming electric field by transforming to the interaction picture with respect to the molecule + lead HamiltonianH = H 0 + H M L . Equation (A1) can then be written in the Liouville space as where Here,T is now the time ordered superoperator which orders the product of superoperators such that their time arguments increase from right to left, i.e., is the density matrix, at time t 0 , which has evolved with the HamiltonianH from some distant past up to time t 0 . We further assume that at time t 0 the reduced molecular part of the full density matrix ρ M (t 0 ) = Tr A, B [ρ T (t 0 )] reaches a steady state which we consider as the initial preparation for the molecule. The interaction with the optical field is switched on after time t 0 . The interaction picture superoperators in Eq. (A2) evolve as To derive the stimulated signal, we first expand the timeordered exponential in Eq. (A3) to third order in light-matter interaction and substitute in Eq. (A2), Now writing down H M F− explicitly and selecting a particular phase component φ = φ 1 − φ 2 + φ 3 − φ 4 with keeping in mind that the pulses are well separated in time, we obtain Time ordering is now maintained by the integration limits and therefore we no longer need the time ordering operator. Writing down the propagators G(t,t ′ ) explicitly, we get Now setting the limit t 0 → −∞ and making the change of time variables to a new set of variables which represents the intervals between the successive light-matter interactions, i.e., t − τ 3 = t 3 , τ 3 − τ 2 = t 2 , τ 2 − τ 1 = t 1 , we get Now assuming that the pulse envelopes are impulsive, i.e., E j (t) = E j δ(t) we can get rid of the integrals and obtain t 1 = T 1 , t 2 = T 2 and t 3 = T 3 . Finally, the integrated stimulated signal S stim =  dtS stim (t) can be solely expressed in terms of these delays and is given in Eq. (6).
The stimulated signal expressions can be further simplified if we assume that the molecular evolution between the optical interaction does not change the state of the molecule and the interaction with the leads contributes only to the life-time of the many-body states. This approximation can be justified in the limit when the energy gaps between different electronic many-body states are much larger than the applied bias. In such a situation, the propagator G will be diagonal, i.e., G ab;c d = G ab δ ac δ bd , where G ab = ⟨⟨ab|G|ab⟩⟩. It is then clear that the resonances due to the charged states will be absent. In terms of the diagrams only (a1), (a3), and (a7) in Fig. 2 contributes to the signal. For the level scheme shown in Fig. 13, we write the following compact expression for the signal: Here, G ab (ω) = i/(ω − ω ab + iΓ), ω ab = ω a − ω b is the energy difference between the state |a⟩ and |b⟩ and Γ = Γ A + Γ B contributes to the life-time of the states.

APPENDIX B: THE QUANTUM MASTER EQUATION FOR THE MANY-BODY MOLECULAR STATES
We write down the quantum master equation in terms of many body states of the isolated molecule |m⟩ where m includes the set of quantum numbers to characterize a particular state of the molecule. The molecular Hamiltonian is given by where E m is the energy eigenvalues for the isolated molecule. The state |m⟩ form an orthonormal basis set in the Fock space. For transport process to be possible, the state |m⟩ could represent different charged states. In the weak molecule-lead coupling limit (second order), the number of electrons in the molecule changes by ±1. Therefore, it is sufficient to consider the molecular states for N and N ± 1 electrons. We write the lead-molecule coupling Hamiltonian given in Eq. (3) in terms of many-body states as where A, B denote two non-interacting metallic leads and T ν mm ′ =  i T iν ⟨m|c i |m ′ ⟩. The overlap amplitudes ⟨m|c i |m ′ ⟩ between different many body states are computed using DFT scheme and the details are given in Appendix D. For simplicity, in our simulation, we assume T iν to be a constant.
The reduced molecular density-matrix evolves according to the Liouville equation or in terms of many-body stateṡ The Liouvillian is obtained assuming weak molecule-lead coupling (second-order), wide-band limit of the leads, and the coupling elements to be real. It can be written as where with f A and f B representing the Fermi functions for lead A and lead B, respectively, and The formal solution of the master equation is given as The matrix element of this propagator is calculated by diagonalising the Liouvillian, where we have used the resolution of identity operator in terms of the left and right eigenvectors of the nonunitary operator L. The index ν runs over all eigenvalues (λ ν ) of L; L|R ν ⟩⟩ = λ ν |R ν ⟩⟩.

APPENDIX C: MOLECULAR JUNCTION PARAMETERS USED FOR NUMERICAL CALCULATION
Here we report the energy eigenvalues and dipole moments for anion and neutral molecule.

The coherent stimulated signal due to initial population
In Fig. 14 we display the stimulated signal for pathway (a1) starting only with population states. In Figs. 15 and 16, we show the contribution due to the population Q pop A (t; Ω 1 ,T 2 , Ω 3 ) and coherence Q coh A (t; Ω 1 ,T 2 , Ω 3 ), i.e., initial state is either a population or in a superposition, to the charge current signal. Similar to the optical signal in this case also we find that the total charge signal is dominated by the contribution from population. The signal due to initial coherences shows one dominant diagonal peak coming from the initial state ρ a ′ b ′ which is much larger compared to the other coherences (see Fig. 4).