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

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–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–17 

Understanding optical signals from current-carrying single molecular systems^18–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.¹⁷ 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,4-dithiol 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.

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₀ = 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,

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_MF(t) is written in the rotating wave approximation (RWA) where $E s ( t )$ represents the positive frequency component of the field and $V ˆ ( 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 ± ϕ₁ ± ϕ₂ ± ϕ₃ ± ϕ₄. The various quantum pathways of the signal can be separated out by the phase cycling protocol which selects a given phase.¹⁴ 

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 τ₁, τ₂, and τ₃ 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³ = 8 ladder diagrams (quantum pathways) as shown in Fig. 2. In our previous work,²¹ 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_LX ≡ AX, and A_RX ≡ XA where X is an arbitrary Hilbert space operator. Furthermore we define linear combinations of L and R superoperators as A₊X ≡ (AX + XA), A₋X ≡ (AX − XA). 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 $G ( t ) ≡θ ( t ) exp [ − i ħ ( H 0 − + H M L − ) t ]$ where θ(t) 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 ρ_M(t₀) = Tr_A,B[ρ_T(t₀)]. The signal is now expressed in terms of the reduced density matrix ρ_M(t₀) 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 L t$ 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₀) and the propagator $G ( t )$⁠.²¹ 

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 Ω₁ and Ω₃ are Fourier frequency conjugates corresponding to the delay times T₁ and T₃ (often called coherence times), respectively. This gives

where $G ( ω ) = ∫ − ∞ ∞ dt e i ω t G ( t ) =i ( ω I − i L ) − 1$⁠. The superoperator correlation function appearing in the signal can be expanded in terms of the many body states of the isolated molecule.²¹ For a general steady state density matrix $ρ= ∑ a b ρ a b |ab〉〉$⁠, the signal is given by

Note that, unlike equilibrium spectroscopy, Eq. (10) depends on population and coherence. Here, $ρ a b =〈a| ρ ˆ |b〉$ is the steady state density matrix element between the many-body states |a〉 and |b〉, $V a b =〈a| V ˆ |b〉$ is the transition dipole matrix element, and $G a b ; c d ( ω ) =i〈〈ab| ( ω I − i L ) − 1 |cd〉〉$ is the non-unitary 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.

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.²⁴ 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,28 Geometry is optimized at the density functional theory (DFT) level with the B3LYP functional^29–31 by using the Gaussian 09 program.³² The Au-Au distance is held at 2.88 Å, taken from the distance in the Au(111) surface.²⁸ The LANL2DZ basis set and pseudopotential³³ 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.³⁴ The time-dependent DFT (TDDFT) method with the Tamm-Dancoff approximation (TDA)³⁵ is used. This method has been successfully employed in simulations of current carrying molecular junctions.^36–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²〉 − 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⁴¹ 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_BT) 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 | ρ s s 〉〉=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₀. 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 ω₁, ω₂, ω₃ 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 Ω₁ and Ω₃ for different delay times T₂. 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 $ρ a a s s$ is small compared to that of $ρ a ′ a ′ s s$ (see Fig. 4) and 〈a′|V|d′〉 > 〈a|V|d〉. Also since for zero delay time (T₂ = 0), the second and third interactions occur simultaneously and therefore the evolution of states during T₁ and T₃ 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₂ = 0 can give rise to off-diagonal peaks. For T₂≠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₁ and T₃. This gives rise to the off-diagonal peaks in the signal when displayed in the Fourier domain as a function of Ω₁ and Ω₃.

Fig. 6 presents different pathways involving only the dominant dipole transitions. Diagrams (a) and (b) generate diagonal peaks for all delay times T₂, whereas the two pathways (c) and (d) are responsible for the cross peaks for T₂≠0. These two pathways are not allowed for T₂ = 0. With increasing T₂, 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₂ 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 stim ( a 1 , coh ) ( Ω 1 , T 2 , Ω 3 )$⁠. The signal due to initial population $S stim ( a 1 , pop ) ( Ω 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₂ = 0, diagram (a) gives rise to peak for |a′〉 − |d′〉 transition. For finite T₂, the states evolving during T₁ are different from those of  T₃ (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₂ = 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 eV). As before, cross peaks appear as the delay time T₂ (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.

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–50 and the quantum master equation^51–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 Hamiltonian $H ̄ = 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 ϕ = ϕ₁ − ϕ₂ + ϕ₃ − ϕ₄, we obtain

We first trace over the lead degrees of freedom and then take the time derivative with respect to t to obtain $Q ̇ ( 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$⁠, and $t ̄ =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₁ and T₃ 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⁴ = 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 Ω₁ and Ω₃ for fixed T₂, 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₂ = 0 and plot the signal for various $t ̄$⁠. For $t ̄ =0fs$ (see panel (a)), the signal shows only one strong diagonal peak (Ω₁ = Ω₃) corresponding to the dominant dipole transition for anion state |a′〉 − |d′〉. This superposition state evolving during T₃ 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 finite $t ̄$ (see panel (b), (c), and (d)), the signal shows significant change and display three strong diagonal peaks 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₁ and T₃ 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 increasing $t ̄$⁠, 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 A pop ( t ̄ ; Ω 1 , T 2 , Ω 3 )$ and coherences $Q A coh ( t ̄ ; Ω 1 , T 2 , Ω 3 )$⁠.

In Fig. 11, we show the charge flow signal Q_A for fixed $t ̄ =10fs$ with different values of T₂. For finite T₂, the crucial observation is the appearance of cross peaks, in addition to diagonal peaks, which was significantly weak for T₂ = 0 fs, as shown in Fig. 10. Moreover similar to the photon detected signal, we also observe here that with increasing T₂ these cross peaks become stronger implying population transfer.

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 consider $t ̄ =0$ limit for the charge current which imply that the observation time t coincides with the arrival time of the fourth pulse $E 4$⁠. The stimulated optical signal shows resonances for both anion and neutral charged states. The integrated charge-current signal for lead A (Q_A), on the other hand, shows resonance peak at Ω₃ which corresponds to the transition only between anion states |a′〉 − |d′〉 for all T₂. As mentioned earlier, the diagonal peaks (Ω₁ = Ω₃) in this case imply that, during T₁ and T₃ 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₂, because of the intermediate propagator $G ( T 2 )$⁠, the coherence state |a′〉 − |d′〉 during T₃ can be reached via various additional quantum pathways starting initially from both the anion and neutral states. Therefore, for T₂≠0, the additional cross peaks show up in the signal for different Ω₁ 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.

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 time-scale 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.

U.H. acknowledges the support from Indian Institute of Science, Bangalore, India. S.M. gratefully acknowledges the support of the National Science Foundation (NSF) through Grant No. CHE-1361516 and from Chemical Sciences, Geosciences, and Biosciences Division, Office of Basic Energy Sciences, Office of Science, (U.S.) Department of Energy (DOE). B.K.A. was supported by DOE. He acknowledges the hospitality at Inorganic and Physical Chemistry department at Indian Institute of Science, Bangalore, where most of the numerical work was performed.

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 Hamiltonian $H ̄ = 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.,

ρ_T(t₀) is the density matrix, at time t₀, which has evolved with the Hamiltonian $H ̄$ from some distant past up to time t₀. We further assume that at time t₀ the reduced molecular part of the full density matrix ρ_M(t₀) = Tr_A,B[ρ_T(t₀)] 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₀. The interaction picture superoperators in Eq. (A2) evolve as

where $G ( t , t 0 ) = e − i ħ ( H 0 − + H M L − ) ( t − t 0 )$⁠.

To derive the stimulated signal, we first expand the time-ordered exponential in Eq. (A3) to third order in light-matter interaction and substitute in Eq. (A2),