Multidimensional spectroscopy with entangled light: loop vs ladder delay scanning protocols

In coherent nonlinear optical spectroscozpy the applied optical pulses induce a polarization in the matter system which is then measured. There are two types of bookkeeping representations for computing an observable (such as the polarization) in a quantum system subjected to time dependent perturbations. Both are exact and should yield the same final results provided no approximations are made. However they offer a very different physical picture and suggest different types of approximations that lead to different predictions.

In the first representation we follow the evolving density matrix in real time. This representation is most suitable for impulsive experiments involving sequences of short, temporally well-separated pulses ranging from NMR to the x-ray regimes [1]. The time variables used to represent the delays between successive pulses [2] ${{t}_{1}}$, ${{t}_{2}}$, ${{t}_{3}}$, ... serve as the primary control parameters. Spectra are displayed vs the Fourier conjugates ${{\overset{\sim}{\mathop{\Omega }}\,}_{1}}$, ${{\overset{\sim}{\mathop{\Omega }}\,}_{2}}$, ${{\overset{\sim}{\mathop{\Omega }}\,}_{3}}$, ... to these variables. Such signals can be represented by ladder diagrams (see figures 1(b) and 2). We shall denote this way of displaying the multidimensional signals as the ladder delay scanning protocol (LAP). The signals with different phase matching directions are distinct when displayed vs ladder delays. The density matrix further allows for reduced descriptions where bath degrees of freedom which cause pure dephasing and relaxation processes are eliminated.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Alternatively we can follow the evolving wave function. Rather than keeping track of both the bra and the ket we can place the entire burden of the time evolution on the ket. In that case we must use artificial time variables where the ket first evolves forward and then backward in time, eventually returning to the initial time. This is represented by loop diagrams [3] as is commonly done in many body theory [4]. This gives a more compact description (fewer terms). It is harder to visualize impulsive experiments in this language. However it proves most useful for frequency domain techniques involving long pulses where the time evolution is not monitored directly [3]. In this picture we give up the full control over time ordering between pulses. We will denote the delays along the loop as ${{\tau }_{1}}$, ${{\tau }_{2}}$, ${{\tau }_{3}}$, ... (see figures 1(a), (c)). By displaying the spectra vs the Fourier conjugates to the loop times ${{\Omega }_{1}}$, ${{\Omega }_{2}}$, ${{\Omega }_{3}}$ we obtain the loop delay scanning protocol (LOP).

In this paper we compare the two display protocols for multidimensional spectroscopy in molecular aggregates with fluorescence detection. Since the two protocols use different time variables the resulting multidimensional signals obtained by Fourier transforms conjugate to these variables appear very different and highlight different resonances. This can be exploited for highlighting desired features in optical signals. We further show some advantages of the loop representation for describing measurements with quantum light, i.e. entangled broadband photons which have intermediate time/frequency character. We should emphasize that these protocols offer two languages for describing the same physics. However the translation is somewhat tricky making them suitable for different applications. We show how such LOP signals can be realized experimentally and compare it to the LAP. Narrow exciton resonances have been predicted recently in two dimensional spectra of aggregates and were attributed to photon entanglement [39]. We show in contrast that these resonances are unrelated to entanglement. Instead they are signatures of the LOP and can be observed with classical light.

The utility of each protocol depends on experimental details including e.g. the system dynamics, bath effects and the specific light field configuration. For instance when the system is in a pure state and the fields are classical, the loop delays ${{\tau }_{j}}$, $j=1,2,3$ which represent forward and backward time propagation periods of the wave function are the natural independent variables and it makes sense to adopt their conjugate frequencies for display, thus using the LOP. If pure dephasing processes due to a bath are added the signal may no longer factorize into a product of terms each depending on a single delay ${{\tau }_{j}}$ when calculating the optical response. In this case the ladder variables ${{t}_{j}}$ which represent the LAP delays in real time and correspond to propagation of a density matrix become more natural since the signal can be recast as a product of individual terms each depending on a single ${{t}_{j}}$ variable. Stochastic or entangled light fields cause additional coupling between the interaction times imposing that the signal may not generally be factorized in either protocol since the field correlation functions depend on products of factors that depend on pairs of times. In that case neither protocol allows the observed signals to be factorized in a simple way discussed above. The two protocols highlight different resonances and processes. In the following we demonstrate what type of information can be extracted from each protocol for Frenkel excitons in a model molecular aggregate.

We further compare signals obtained with classical vs quantum light (entangled photons). The LAP and LOP denote the protocols for displaying multidimensional signals. Calculations performed with either the wavefunction or the density matrix can be displayed using either protocol. In earlier studies ladder diagrams were denoted as double-sided Feynman diagrams, and loop diagrams were denoted as close-time-path-loops (CTPL) [3].

We investigate the multidimensional signals in a molecular aggregate obtained by incoherent two-photon absorption (TPA) detection. Incoherent detection is often more sensitive than heterodyne as the latter is limited by the pulse duration so there are fewer constraints on the laser system. In addition the low intensity requirements for biological samples limit the range of heterodyne detection setups. This have been demonstrated [5–7] even in single molecule spectroscopy [8]. Historically Ramsey fringes constitute the first example of incoherent detection [9–11]. Information similar to coherent spectroscopy can be extracted from the parametric dependence on various pulse sequences applied prior to the incoherent detection [12, 13]. Possible incoherent detection modes include fluorescence [14–16], photoaccoustic [17–19], AFM [20–23] or photocurrent detection [24, 25].

Quantum spectroscopy which utilizes the quantum nature of light to reveal matter properties is an emerging field. Entangled photons is one notable example and offer several advantages. First, the signals scale to lower order in the incoming intensity [26]. The pump-probe signal for example scales linearly rather than quadratically. This allows to to perform nonlinear spectroscopy with much lower intensity limiting damage in e.g. imaging applications [26–34]. Second, time-and-frequency entanglement often allows to obtain higher temporal and spectral resolutions since the two are not Fourier conjugates. Namely, the temporal resolution $\Delta t$ depends on the length of the nonlinear crystal, that is, the entanglement time T, while spectral resolution $\Delta \omega$ is determined by the pump envelope. These are independent control variables, not Fourier conjugates and not bound by the uncertainty $\Delta \omega \Delta t1$. We show that entangled photons allow to observe narrow spectral features even in the limit of fast dephasing where the classical line shapes are broad. Elaborate pulse shaping techniques that involve standard prisms compressors and spatial light modulators [35–38] can be used to control the amplitude and phase modulation of entangled photon pairs necessary for creating the desired pulse sequence. This can be done using e.g. the Franson interferometer with variable phases and delays in both arms of the interferometer as proposed in [39]. The beam splitters in two arms allow to create four pulses using a single entangled photon pair. In the following we do not specify the experimental details of shaping the pulses, rather we assume a generic sequence of shaped entangled photons.

We consider a model system of an aggregate described by the Frenkel exciton Hamiltonian

$$\begin{eqnarray}&&H={{H}_{0}}+{{H}^{\prime }},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{{H}_{0}}=\hbar \underset{m}{\mathop{\mathop{\sum }}}\,{{\epsilon }_{m}}B_{m}^{\dagger }{{B}_{m}}+\hbar \underset{m\ne n}{\mathop{\mathop{\sum }}}\,{{J}_{mn}}B_{m}^{\dagger }{{B}_{n}}+\hbar \underset{m}{\mathop{\mathop{\sum }}}\,\frac{{{\Delta }_{m}}}{2}B_{m}^{\dagger }B_{m}^{\dagger }{{B}_{m}}{{B}_{m}},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{{H}^{\prime }}=E(t){{V}^{\dagger }}+{{E}^{\dagger }}(t)V,\quad {{V}^{\dagger }}=\underset{m}{\mathop{\mathop{\sum }}}\,V_{m}^{\dagger }B_{m}^{\dagger },\end{eqnarray} \tag{ 3 }$$

where ${{H}_{0}}$ is the excitonic part, ${{\epsilon }_{m}}$ are site energies, ${{J}_{mn}}$ are hopping and ${{\Delta }_{m}}$ is an onsite repulsion (Hubbard type), and ${{B}_{m}}$ is an exciton Pauli annihilation operator at site m (e.g. pigment or quantum dot). ${{H}^{\prime }}$ is the dipole interaction with the optical field E in the rotating wave approximation. E is the electric field operator. The eigenstate of equation (2) form distinct exciton bands (see figure 1(d)). In the diagonal eigenstate representation the Hamiltonian for the lowest three manifold of states which are relevant for the present study - ground g, single excited e and double excited f manifolds (see figure 1(d)) reads

$$\begin{eqnarray}&&{{H}_{0}}=\hbar {{\omega }_{g}}|g\rangle \langle g|+\hbar \underset{e}{\mathop{\mathop{\sum }}}\,{{\omega }_{e}}|e\rangle \langle e|+\hbar \underset{f}{\mathop{\mathop{\sum }}}\,{{\omega }_{f}}|f\rangle \langle f|,\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{{V}^{\dagger }}=\underset{e}{\mathop{\mathop{\sum }}}\,V_{ge}^{*}|e\rangle \langle g|+\underset{e,f}{\mathop{\mathop{\sum }}}\,V_{ef}^{*}|f\rangle \langle e|.\end{eqnarray} \tag{ 5 }$$

We consider the following experiment: a sequence of four pulses centered at times ${{T}_{a}}$, ${{T}_{b}}$, ${{T}_{c}}$, and ${{T}_{d}}$ with phases ${{\phi }_{a}}$, ${{\phi }_{b}}$, ${{\phi }_{c}}$, and ${{\phi }_{d}}$ [40], given by the electric field $E={{\mathop{\sum }}_{\alpha =a,b,c,d}}{{E}_{\alpha }}{{e}^{i{{\phi }_{\alpha }}}}$ brings the molecule into its doubly-excite state [41] (see figures 1(a), (b)) and the population of f states is detected. This can be done by fluorescence $f\to e$ or after a rapid internal conversion process the molecule can be deexcited from f to e and fluorescence from e to g is then detected. We assume that the $e\to g$ and $f\to e$ channels can be distinguished in time or frequency and therefore we can isolate the TPA contributions. Thus, we define the signal as the sum of populations of states f.

$$\begin{eqnarray}&&S(\Gamma )=\underset{f}{\mathop{\mathop{\sum }}}\,{{\rho }_{ff}}(\Gamma ),\end{eqnarray} \tag{ 6 }$$

where $\Gamma$ represents collectively the set of parameters of the incoming pulses. These depend on the protocol and will be specified later.

The signal (6) for our model is given by the single unrestricted loop diagram in figure 1(c) (for diagram rules see [3]). $a,b,c,d$ denote the pulse sequence ordered along the loop (not in real time); a represents 'first' on the loop etc. Pulses chronologically-ordered in real time will be denoted $1,2,3,4$ which are permutations of $a,b,c,d$ determined by the time arguments, as will be shown below. One can scan various delays ${{T}_{\alpha }}-{{T}_{\beta }}$, $\alpha ,\beta =a,b,c,d$ and control the phases $\pm {{\phi }_{a}}\pm {{\phi }_{b}}\pm {{\phi }_{c}}\pm {{\phi }_{d}}$. Phase cycling techniques have been successfully demonstrated as a control tool for the selection of fixed-phase components of optical signals generated by multiwave mixing [42–46]. Phase cycling can be easily implemented using a pulse shaper by varying the relative inter-pulse phases, which are cycled over 2 radians in a number of equally spaced steps [42, 43]. To realize the LOP experimentally the indices a, b, c, d are assigned as follows: first by phase cycling we select a signal with phase ${{\phi }_{a}}+{{\phi }_{b}}-{{\phi }_{c}}-{{\phi }_{d}}$. The two pulses with positive phase detection are thus denoted a, b and with negative phase—c, d. In the a, b pair pulse a comes first. In the c, d pair pulse d comes first. The time variables in figure 1(c) are ${{\tau }_{1}}={{T}_{b}}-{{T}_{a}}$, ${{\tau }_{2}}={{T}_{c}}-{{T}_{b}}$, ${{\tau }_{3}}={{T}_{c}}-{{T}_{d}}$. With this choice ${{\tau }_{1}}$ and ${{\tau }_{3}}$ are positive whereas ${{\tau }_{2}}$ can be either positive or negative. This completely defines the LOP experimentally.

2.1. Pure states and the loop representation

In figure 1(c) two interactions with bra- and two - with ket- promote the system to the state described by a population density matrix element ${{\rho }_{ff}}$. In the following we omit the phase factor ${{e}^{i({{\phi }_{a}}+{{\phi }_{b}}-{{\phi }_{c}}-{{\phi }_{d}})}}$, keeping in mind that all the signals contain it. The corresponding signal (6) can be read-off the diagrams

$$\begin{eqnarray}\begin{matrix} S(\Gamma ) & = & \frac{1}{{{\hbar }^{4}}}\mathop{\int }_{-\infty }^{\infty }d{{r}_{a}}\mathop{\int }_{-\infty }^{\infty }d{{r}_{b}}\mathop{\int }_{-\infty }^{\infty }d{{r}_{c}}\mathop{\int }_{-\infty }^{\infty }d{{r}_{d}}\langle E_{d}^{\dagger }({{r}_{d}})E_{c}^{\dagger }({{r}_{c}}){{E}_{b}}({{r}_{b}}){{E}_{a}}({{r}_{a}})\rangle \\ & & \times \langle \mathcal{T}V({{r}_{d}})V{{r}_{c}}){{V}^{\dagger }}({{r}_{b}}){{V}^{\dagger }}({{r}_{a}})\rangle . \\ &&\end{matrix}\end{eqnarray} \tag{ 7 }$$

Here ${{r}_{\alpha }}$, $\alpha =a,b,c,d$ are the interaction times of our four pulses with the aggregate, $\mathcal{T}$ denotes the time ordering operator along the loop [47]. Equation (7) can be recast using the loop intervals figure 1(c) ${{s}_{j}}$, $j=1,2,3$

$$\begin{eqnarray}\begin{matrix} S(\Gamma ) & = & \mathcal{R}\frac{1}{{{\hbar }^{4}}}\mathop{\int }_{-\infty }^{\infty }d{{r}_{b}}\mathop{\int }_{0}^{\infty }d{{s}_{1}}\mathop{\int }_{0}^{\infty }d{{s}_{2}}\mathop{\int }_{0}^{\infty }d{{s}_{3}}\langle E_{d}^{\dagger }({{r}_{b}}-{{s}_{2}}-{{s}_{3}})E_{c}^{\dagger }({{r}_{b}}-{{s}_{2}}){{E}_{b}}({{r}_{b}}){{E}_{a}}({{r}_{b}}-{{s}_{1}})\rangle \\ & & \times \langle V({{r}_{b}}-{{s}_{2}}-{{s}_{3}})V({{r}_{b}}-{{s}_{2}}){{V}^{\dagger }}({{r}_{b}}){{V}^{\dagger }}({{r}_{b}}-{{s}_{1}}). \\ &&\end{matrix}\end{eqnarray} \tag{ 8 }$$

Time ordering is now explicitly specified by the integration limits and we no longer need the time ordering operator. In this expression ${{s}_{2}}$ is positive (interaction with pulse b is chronologically the last). The contribution where the field c is the last is included by taking the real part $\mathcal{R}$.

One can alternatively recast equation (8) in frequency-domain using the electric field operators ${{\mathcal{E}}_{\alpha }}(t)=\mathop{\int }_{-\infty }^{\infty }\frac{d\omega }{2\pi }E(\omega ){{e}^{-i\omega (t-{{T}_{\alpha }})}}$, $\alpha =a,b,c,d$

$$\begin{eqnarray}\begin{matrix} {{S}_{LOP}}({{\tau }_{1}},{{\tau }_{2}},{{\tau }_{3}}) & = & \mathcal{R}\frac{i}{{{\hbar }^{4}}}\mathop{\int }_{-\infty }^{\infty }\frac{d{{\omega }_{a}}}{2\pi }\mathop{\int }_{-\infty }^{\infty }\frac{d{{\omega }_{b}}}{2\pi }\mathop{\int }_{-\infty }^{\infty }\frac{d{{\omega }_{d}}}{2\pi }\langle {{E}^{\dagger }}({{\omega }_{d}}){{E}^{\dagger }}({{\omega }_{a}}+{{\omega }_{b}}-{{\omega }_{d}})E({{\omega }_{b}})E({{\omega }_{a}})\rangle \\ & & \times \underset{e,{{e}^{\prime }},f}{\mathop{\mathop{\sum }}}\,\mu _{g{{e}^{\prime }}}^{a}\mu _{{{e}^{\prime }}f}^{b}\mu _{fe}^{c*}\mu _{eg}^{d*}G_{e\prime }^{\dagger }({{\omega }_{d}}){{G}_{e}}({{\omega }_{a}})[{{G}_{f}}({{\omega }_{a}}+{{\omega }_{b}})\theta ({{\tau }_{2}})-G_{f}^{\dagger }({{\omega }_{a}}+{{\omega }_{b}})\theta (-{{\tau }_{2}})] \\ & & \times {{e}^{-i{{\omega }_{a}}{{\tau }_{1}}+i{{\omega }_{d}}{{\tau }_{3}}-i({{\omega }_{a}}+{{\omega }_{b}}){{\tau }_{2}}}}\theta ({{\tau }_{1}})\theta ({{\tau }_{3}}), \\ &&\end{matrix}\end{eqnarray} \tag{ 9 }$$

where the LOP control variables ${{\tau }_{1}}={{T}_{b}}-{{T}_{a}}$, ${{\tau }_{2}}={{T}_{c}}-{{T}_{b}}$, ${{\tau }_{3}}={{T}_{c}}-{{T}_{d}}$ are the delays between pulse centers and $G(\omega )=1/[\omega +{{\omega }_{g}}-H/\hbar +i\epsilon ]$ is a frequency domain Green's function. $\mu _{j{{j}^{\prime }}}^{\alpha }={{V}_{j{{j}^{\prime }}}}\cdot {{\sigma }_{\alpha }}$ is the projection of the transition dipole moment ${{V}_{j{{j}^{\prime }}}}$, $j,{{j}^{\prime }}=g,e,{{e}^{\prime }},f$ onto the polarization vector ${{\sigma }_{\alpha }}$ of the corresponding field $\alpha =a,b,c,d$.

In the frequency-domain the field correlation function is defined as a Fourier transform of the time-domain field correlation function

$$\begin{eqnarray} &&\langle {{E}^{\dagger }}({{\omega }_{d}}){{E}^{\dagger }}({{\omega }_{c}})E({{\omega }_{b}})E({{\omega }_{a}})\rangle \\ &&=\mathop{\int }_{-\infty }^{\infty }dt_{1}^{\prime }\mathop{\int }_{-\infty }^{\infty }dt_{2}^{\prime }\mathop{\int }_{-\infty }^{\infty }dt_{3}^{\prime }\mathop{\int }_{-\infty }^{\infty }dt_{4}^{\prime }{{e}^{i{{\omega }_{a}}t_{1}^{\prime }+i{{\omega }_{b}}t_{2}^{\prime }-i{{\omega }_{c}}t_{3}^{\prime }-i{{\omega }_{d}}t_{4}^{\prime }}}\langle E_{d}^{\dagger }(t_{3}^{\prime })E_{c}^{\dagger }(t_{4}^{\prime }){{E}_{b}}(t_{2}^{\prime }){{E}_{a}}(t_{1}^{\prime })\rangle . \\ \end{eqnarray} \tag{ 10 }$$

In equation (9) we used equation (10) and the time translation invariance symmetry which implies ${{\omega }_{a}}+{{\omega }_{b}}-{{\omega }_{c}}-{{\omega }_{d}}=0$.

In the absence of a bath, the matter correlation function is given by

$$\begin{eqnarray} &&\underset{e,{{e}^{\prime }},f}{\mathop{\mathop{\sum }}}\,\mu _{g{{e}^{\prime }}}^{a}\mu _{{{e}^{\prime }}f}^{b}\mu _{fe}^{c*}\mu _{eg}^{d*}G_{e\prime }^{\dagger }({{\omega }_{d}}){{G}_{e}}({{\omega }_{a}})[{{G}_{f}}({{\omega }_{a}}+{{\omega }_{b}})\theta ({{\tau }_{2}})-G_{f}^{\dagger }({{\omega }_{a}}+{{\omega }_{b}})\theta (-{{\tau }_{2}})] \\ &&=\mathop{\int }_{-\infty }^{\infty }dt_{1}^{\prime }\mathop{\int }_{-\infty }^{\infty }dt_{2}^{\prime }\mathop{\int }_{-\infty }^{\infty }dt_{3}^{\prime }\mathop{\int }_{-\infty }^{\infty }dt_{4}^{\prime }{{e}^{-i{{\omega }_{a}}t_{1}^{\prime }-i{{\omega }_{b}}t_{2}^{\prime }+i{{\omega }_{c}}t_{3}^{\prime }+i{{\omega }_{d}}t_{4}^{\prime }}}\langle \mathcal{T}{{V}_{g{{e}^{\prime }}}}(t_{3}^{\prime }){{V}_{{{e}^{\prime }}f}}(t_{4}^{\prime })V_{fe}^{\dagger }(t_{2}^{\prime })V_{eg}^{\dagger }(t_{1}^{\prime })\rangle . \\ \end{eqnarray} \tag{ 11 }$$

A Fourier transform of (9) with respect to loop delays then gives a 3D signal

$$\begin{eqnarray}&&{{S}_{LOP}}({{\Omega }_{1}},{{\Omega }_{2}},{{\Omega }_{3}})=\mathop{\int }_{-\infty }^{\infty }d{{\tau }_{1}}\mathop{\int }_{-\infty }^{\infty }d{{\tau }_{2}}\mathop{\int }_{-\infty }^{\infty }d{{\tau }_{3}}{{e}^{i{{\Omega }_{1}}{{\tau }_{1}}+i{{\Omega }_{2}}{{\tau }_{2}}+i{{\Omega }_{3}}{{\tau }_{3}}}}{{S}_{LOP}}({{\tau }_{1}},{{\tau }_{2}},{{\tau }_{3}}).\end{eqnarray} \tag{ 12 }$$

Combining equations (7)–(12) gives

$$\begin{eqnarray}&&{{S}_{LOP}}({{\Omega }_{1}},{{\Omega }_{2}},{{\Omega }_{3}})=S_{LOP}^{\prime }({{\Omega }_{1}},{{\Omega }_{2}},{{\Omega }_{3}})+S_{LOP}^{\prime *}(-{{\Omega }_{1}},-{{\Omega }_{2}},-{{\Omega }_{3}}),\end{eqnarray} \tag{ 13 }$$

where

$$\begin{eqnarray} S_{LOP}^{\prime }({{\Omega }_{1}},{{\Omega }_{2}},{{\Omega }_{3}}) & = & \frac{2}{{{\hbar }^{4}}}\mathop{\int }_{-\infty }^{\infty }\frac{d{{\omega }_{a}}}{2\pi }\mathop{\int }_{-\infty }^{\infty }\frac{d{{\omega }_{d}}}{2\pi }\langle E_{d}^{\dagger }({{\omega }_{d}})E_{c}^{\dagger }({{\omega }_{fg}}-{{\omega }_{d}}){{E}_{b}}({{\omega }_{fg}}-{{\omega }_{a}}){{E}_{a}}({{\omega }_{a}})\rangle \\ & & \times \underset{e,{{e}^{\prime }},f}{\mathop{\mathop{\sum }}}\,\frac{\mu _{g{{e}^{\prime }}}^{a}\mu _{{{e}^{\prime }}f}^{b}\mu _{fe}^{c*}\mu _{eg}^{d*}{{G}_{e}}({{\omega }_{a}})G_{e\prime }^{\dagger }({{\omega }_{d}})}{[{{\omega }_{a}}-{{\Omega }_{1}}-i\epsilon ][{{\omega }_{d}}+{{\Omega }_{3}}+i\epsilon ]}\left[ \frac{{{G}_{f}}({{\omega }_{a}}+{{\omega }_{b}})}{{{\omega }_{a}}+{{\omega }_{b}}-{{\Omega }_{2}}-i\epsilon } \right.\left. +\frac{G_{f}^{\dagger }({{\omega }_{a}}+{{\omega }_{b}})}{{{\omega }_{a}}+{{\omega }_{b}}-{{\Omega }_{2}}+i\epsilon } \right]. \\ & & \\ \end{eqnarray} \tag{ 14 }$$

where the limit $\epsilon \to 0$ is understood. One can then evaluate the remaining frequency integrals in equation (14) for a given light field correlation function using residue calculus.

So far we did not specify the nature of the field, and equations (9)–(14) hold for arbitrary type of field, be it classical, stochastic or entangled. All relevant field information is contained in its four point field correlation function which must be evaluated separately. For classical coherent fields this function factorizes (in time or frequency) into a product of four amplitudes. Otherwise for entangled or stochastic fields the correlation function causes a coupling between two interaction times, which affects the signals.

2.2. Pure dephasing, bath effects and the ladder representation

When the exciton system is coupled to a bath, it can no longer be described by a wavefunction once the bath is eliminated. To evaluate the loop diagram it must be broken into several ladder diagrams (for notation see [3]) which represent the density matrix. The unrestricted loop diagram in figure 1(b) is split into the six ladder diagrams shown in figure 1(c) and the signal (9) is given by sum of all six terms ${{S}_{LOP}}({{\tau }_{1}},{{\tau }_{2}},{{\tau }_{3}})=\mathop{\sum }_{j=1}^{6}S_{LOP}^{(j)}({{\tau }_{1}},{{\tau }_{2}},{{\tau }_{3}})$ where

$$\begin{eqnarray}&&S_{LOP}^{(j)}({{\tau }_{1}},{{\tau }_{2}},{{\tau }_{3}})=\mathop{\int }_{-\infty }^{\infty }\frac{d{{\omega }_{a}}}{2\pi }\frac{d{{\omega }_{b}}}{2\pi }\frac{d{{\omega }_{d}}}{2\pi }D_{LOP}^{(j)}({{\tau }_{1}},{{\tau }_{2}},{{\tau }_{3}};{{\omega }_{a}},{{\omega }_{b}},{{\omega }_{d}}){{\overset{\sim}{\mathop{S}}\,}^{(j)}}({{\omega }_{a}},{{\omega }_{b}},{{\omega }_{d}})-c.c.,\end{eqnarray} \tag{ 15 }$$

where

$$\begin{eqnarray}&&D_{LOP}^{(j)}({{\tau }_{1}},{{\tau }_{2}},{{\tau }_{3}};{{\omega }_{a}},{{\omega }_{b}},{{\omega }_{d}})=\theta ({{\tau }_{1}}){{\theta }_{j}}(\pm {{\tau }_{2}})\theta ({{\tau }_{3}}){{e}^{-i{{\omega }_{a}}{{\tau }_{1}}+i{{\omega }_{d}}{{\tau }_{3}}-i({{\omega }_{a}}+{{\omega }_{b}}){{\tau }_{2}}}}\end{eqnarray} \tag{ 16 }$$

is a display function which depends on the control parameters specific to the chosen protocol. In ${{\theta }_{j}}(\pm {{\tau }_{2}})$ the 'minus' sign applies for diagrams $j=1,2,3$ and the 'plus' sign for $j=4,5,6$,

$$\begin{eqnarray}&&{{\overset{\sim}{\mathop{S}}\,}^{(j)}}({{\omega }_{a}},{{\omega }_{b}},{{\omega }_{d}})=\langle {{E}^{\dagger }}({{\omega }_{d}}){{E}^{\dagger }}({{\omega }_{a}}+{{\omega }_{b}}-{{\omega }_{d}})E({{\omega }_{b}})E({{\omega }_{a}})\rangle {{R}^{(j)}}({{\omega }_{a}},{{\omega }_{b}},{{\omega }_{d}}),\end{eqnarray} \tag{ 17 }$$

and

$$\begin{eqnarray}\begin{matrix} R^{\left(1\right)}\left(\omega_a,\omega_b,\omega_d\right)&=&\hbar^{-4}\sum_{e,e^{\prime},f}\mu_{ge^{\prime}}^a\mu_{e^{\prime}f}^b\mu_{fe}^{c^{\ast}}\mu_{eg}^{d^{\ast}}\mathcal{G}_{ge^{\prime}}\left(-\omega_d\right)\mathcal{G}_{ef}\left(-\omega_b\right)\mathcal{G}_{gf}\left(-\omega_a-\omega_b\right),\\ R^{\left(2\right)}\left(\omega_a,\omega_b,\omega_d\right)&=& \hbar^{-4}\sum_{e,e^{\prime},f}\mu_{ge^{\prime}}^a\mu_{e^{\prime}f}^b\mu_{fe}^{c^{\ast}}\mu_{eg}^{d^{\ast}}\mathcal{G}_{eg}\left(\omega_a\right)\mathcal{G}_{ef}\left(-\omega_b\right)\mathcal{G}_{ee^{\prime}}\left(\omega_a-\omega_d\right),\\ R^{\left(3\right)}\left(\omega_a,\omega_b,\omega_d\right)&=& \hbar^{-4}\sum_{e,e^{\prime},f}\mu_{ge^{\prime}}^a\mu_{e^{\prime}f}^b\mu_{fe}^{c^{\ast}}\mu_{eg}^{d^{\ast}}\mathcal{G}_{ge^{\prime}}\left(-\omega_d\right)\mathcal{G}_{ef}\left(-\omega_b\right)\mathcal{G}_{ee^{\prime}}\left(\omega_a-\omega_d\right),\\ R^{\left(4\right)}\left(\omega_a,\omega_b,\omega_d\right)&=& \hbar^{-4}\sum_{e,e^{\prime},f}\mu_{ge^{\prime}}^a\mu_{e^{\prime}f}^b\mu_{fe}^{c^{\ast}}\mu_{eg}^{d^{\ast}}\mathcal{G}_{eg}\left(\omega_a\right)\mathcal{G}_{fe^{\prime}}\left(\omega_a+\omega_b-\omega_d\right)\mathcal{G}_{fg}\left(\omega_a+\omega_b\right),\\ R^{\left(5\right)}\left(\omega_a,\omega_b,\omega_d\right)&=& \hbar^{-4}\sum_{e,e^{\prime},f}\mu_{ge^{\prime}}^a\mu_{e^{\prime}f}^b\mu_{fe}^{c^{\ast}}\mu_{eg}^{d^{\ast}}\mathcal{G}_{ge^{\prime}}\left(-\omega_d\right)\mathcal{G}_{fe^{\prime}}\left(\omega_a+\omega_b-\omega_d\right)\mathcal{G}_{ee^{\prime}}\left(\omega_a-\omega_d\right),\\ R^{\left(6\right)}\left(\omega_a,\omega_b,\omega_d\right)&=&\hbar^{-4}\sum_{e,e^{\prime},f}\mu_{ge^{\prime}}^a\mu_{e^{\prime}f}^b\mu_{fe}^{c^{\ast}}\mu_{eg}^{d^{\ast}}\mathcal{G}_{eg}\left(\omega_a\right)\mathcal{G}_{fe^{\prime}}\left(\omega_a+\omega_b-\omega_d\right)\mathcal{G}_{ee^{\prime}}\left(\omega_a-\omega_d\right). \end{matrix}\end{eqnarray} \tag{ 18 }$$

Here we had introduced the Liouville space Green's function

$$$ {{\mathcal{G}}_{\alpha \beta }}(\omega )=-i\mathop{\int }_{0}^{\infty }dt{{e}^{i\omega t}}{{\langle {{G}_{\alpha }}(t)G_{\beta }^{\dagger }(t)\rangle }_{B}}, \tag{ 19 }$$

where ${{\langle \ldots \rangle }_{B}}=\text{Tr}[\ldots {{\rho }_{B}}]$ represents the trace over the bath degrees of freedom. $G(t)=-i\theta (t){{e}^{-iHt/\hbar }}$ is the Hilbert space Green's function, $\theta (t)$ is the Heaviside step function. The bra and the ket evolutions (and the corresponding time variables) are now coupled by the bath. The effect of couplings between interaction times due to nonclassical field can be demonstrated by evaluating the frequency integrals in equation (15) using time-domain display function in equation (16). The result for entangled photons is given in appendix A. To see how this effect translates on the mixing of the frequency variables we take a Fourier transform of equation (15) with respect to loop delay variable ${{\tau }_{j}}$, $j=1,2,3$ and obtain the signal

$$\begin{eqnarray}\begin{matrix} {{S}_{LOP}}({{\Omega }_{1}},{{\Omega }_{2}},{{\Omega }_{3}}) & = & \underset{6}{\overset{j=1}{\mathop{\mathop{\sum }}}}\,\mathop{\int }_{-\infty }^{\infty }\frac{d{{\omega }_{a}}}{2\pi }\frac{d{{\omega }_{b}}}{2\pi }\frac{d{{\omega }_{d}}}{2\pi }[D_{LOP}^{(j)}({{\Omega }_{1}},{{\Omega }_{2}},{{\Omega }_{3}};{{\omega }_{a}},{{\omega }_{b}},{{\omega }_{d}}){{\overset{\sim}{\mathop{S}}\,}^{(j)}}({{\omega }_{a}},{{\omega }_{b}},{{\omega }_{d}}) \\ & & -D_{LOP}^{(j)*}(-{{\Omega }_{1}},-{{\Omega }_{2}},-{{\Omega }_{3}};{{\omega }_{a}},{{\omega }_{b}},{{\omega }_{d}}){{\overset{\sim}{\mathop{S}}\,}^{(j)*}}({{\omega }_{a}},{{\omega }_{b}},{{\omega }_{d}})], \\ &&\end{matrix}\end{eqnarray} \tag{ 20 }$$

where

$$\begin{eqnarray}&&D_{LOP}^{(j)}({{\Omega }_{1}},{{\Omega }_{2}},{{\Omega }_{3}};{{\omega }_{a}},{{\omega }_{b}},{{\omega }_{d}})=\frac{\mp i}{[{{\omega }_{a}}-{{\Omega }_{1}}-i\epsilon ][{{\Omega }_{3}}+{{\omega }_{d}}+i\epsilon ][{{\Omega }_{2}}-{{\omega }_{a}}-{{\omega }_{b}}\mp i\epsilon ]},\end{eqnarray} \tag{ 21 }$$

and minus (plus) sign corresponds to contributions of diagrams 1–3 (4–6). The coupling between interaction times now translates into a mixing of their conjugate frequency variables ${{\Omega }_{j}}$, $j=1,2,3$. The 3D signals (20) are given by a 3D spectral overlap between Green's functions of the matter and field, where the latter are governed by ${{[{{\Omega }_{j}}-{{\omega }_{\alpha }}\pm i\epsilon ]}^{-1}}$ dressed by a four point field correlation function which selects the field-matter pathways. The response of the system to classical light fields is given by nonlinear response functions which can be expressed by sums over various quantum pathways of matter. In the case of quantum field the response is typically treated in the joint field-matter space to account for back-reaction and other nonclassical effects of the field. In this case the response is summed over various quantum pathways in the joint field-matter space. Depending on the field parameters some quantum pathways can be suppressed or enhanced. The field correlation function controls the relevant spectral range of the pathways that contribute to the signal. Different integrations may couple various frequencies ${{\omega }_{\alpha }}$, $\alpha =a,b,d$ into a single field-matter Green's function. Upon evaluating the relevant frequency integrations different ${{\Omega }_{j}}$, $j=1,2,3$ will be coupled. This will result in various cross-peaks between ${{\Omega }_{j}}$ variables, as becomes apparent by comparing a field contribution in equation (21) with various responses in equation (18). Together with the bath dephasing effects, the relevant spectral width of these cross-peaks can vary significantly compared to that of the system without bath interacting with classical fields. Below we will investigate the signatures of the bath and the state of field in the signals.

In standard multidimensional techniques, the time variables represent the pulses as they interact with sample in chronological order [1]. These are conveniently given by the ladder delays. In the LOP the time ordering between pulses is maintained only on each branch of the loop but not between branches. The LAP in contrast involves full time-ordering of all four pulses. The arrival time of the various pulses in chronological order is ${{T}_{1}}<{{T}_{2}}<{{T}_{3}}<{{T}_{4}}$. The indices $1,2,3,4$ are some permutation of $a,b,c,d$ depending on the diagram. The ladder delays are defined as ${{t}_{1}}={{T}_{2}}-{{T}_{1}}$, ${{t}_{2}}={{T}_{3}}-{{T}_{2}}$, ${{t}_{3}}={{T}_{4}}-{{T}_{3}}$. Ladder diagrams keep track of chronological delays. Each ladder diagram will have its own set of relations between ${{t}_{j}}$, $j=1,2,3$ and pulse delays ${{T}_{\alpha }}-{{T}_{\beta }}$, $\alpha ,\beta =a,b,c,d$. One can then use the phase cycling to select the diagrams shown in figure 2 e.g. ${{\mathbf{k}}_{I}}=-{{\mathbf{k}}_{1}}+{{\mathbf{k}}_{2}}+{{\mathbf{k}}_{3}}$, ${{\mathbf{k}}_{II}}={{\mathbf{k}}_{1}}-{{\mathbf{k}}_{2}}+{{\mathbf{k}}_{3}}$ and ${{\mathbf{k}}_{III}}={{\mathbf{k}}_{1}}+{{\mathbf{k}}_{2}}-{{\mathbf{k}}_{3}}$. This gives

$$\begin{eqnarray}&&{{S}_{{{\mathbf{k}}_{I}}}}({{t}_{1}},{{t}_{2}},{{t}_{3}})=S_{LAP}^{(2)}({{t}_{1}},{{t}_{2}},{{t}_{3}})+S_{LAP}^{(5)}({{t}_{1}},{{t}_{2}},{{t}_{3}}),\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&{{S}_{{{\mathbf{k}}_{II}}}}({{t}_{1}},{{t}_{2}},{{t}_{3}})=S_{LAP}^{(3)}({{t}_{1}},{{t}_{2}},{{t}_{3}})+S_{LAP}^{(6)}({{t}_{1}},{{t}_{2}},{{t}_{3}}),\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&{{S}_{{{\mathbf{k}}_{III}}}}({{t}_{1}},{{t}_{2}},{{t}_{3}})=S_{LAP}^{(1)}({{t}_{1}},{{t}_{2}},{{t}_{3}})+S_{LAP}^{(4)}({{t}_{1}},{{t}_{2}},{{t}_{3}}).\end{eqnarray} \tag{ 24 }$$

where

$$\begin{eqnarray}&&S_{LAP}^{(j)}({{t}_{1}},{{t}_{2}},{{t}_{3}})=\mathop{\int }_{-\infty }^{\infty }\frac{d{{\omega }_{a}}}{2\pi }\frac{d{{\omega }_{b}}}{2\pi }\frac{d{{\omega }_{d}}}{2\pi }D_{LAP}^{(j)}({{t}_{1}},{{t}_{2}},{{t}_{3}};{{\omega }_{a}},{{\omega }_{b}},{{\omega }_{d}}){{\overset{\sim}{\mathop{S}}\,}^{(j)}}({{\omega }_{a}},{{\omega }_{b}},{{\omega }_{d}})-c.c..\end{eqnarray} \tag{ 25 }$$

Here the LAP display functions are given by

$$\begin{eqnarray}&&D_{LAP}^{(1)}({{t}_{1}},{{t}_{2}},{{t}_{3}};{{\omega }_{a}},{{\omega }_{b}},{{\omega }_{d}})=\theta ({{t}_{1}})\theta ({{t}_{2}})\theta ({{t}_{3}}){{e}^{i{{\omega }_{b}}{{t}_{3}}+i{{\omega }_{d}}{{t}_{1}}+i({{\omega }_{a}}+{{\omega }_{b}}){{t}_{2}}}},\end{eqnarray} \tag{ 26 }$$

where ${{t}_{3}}={{T}_{b}}-{{T}_{a}}$, ${{t}_{2}}={{T}_{a}}-{{T}_{c}}$, ${{t}_{1}}={{T}_{c}}-{{T}_{d}}$.

$$\begin{eqnarray}&&D_{LAP}^{(2)}({{t}_{1}},{{t}_{2}},{{t}_{3}};{{\omega }_{a}},{{\omega }_{b}},{{\omega }_{d}})=\theta ({{t}_{1}})\theta ({{t}_{2}})\theta ({{t}_{3}}){{e}^{i{{\omega }_{b}}{{t}_{3}}-i{{\omega }_{a}}{{t}_{1}}+i({{\omega }_{d}}-{{\omega }_{a}}){{t}_{2}}}},\end{eqnarray} \tag{ 27 }$$

where ${{t}_{3}}={{T}_{b}}-{{T}_{c}}$, ${{t}_{2}}={{T}_{c}}-{{T}_{d}}$, ${{t}_{1}}={{T}_{d}}-{{T}_{a}}$.

$$\begin{eqnarray}&&D_{LAP}^{(3)}({{t}_{1}},{{t}_{2}},{{t}_{3}};{{\omega }_{a}},{{\omega }_{b}},{{\omega }_{d}})=\theta ({{t}_{1}})\theta ({{t}_{2}})\theta ({{t}_{3}}){{e}^{i{{\omega }_{b}}{{t}_{3}}+i{{\omega }_{d}}{{t}_{1}}+i({{\omega }_{d}}-{{\omega }_{a}}){{t}_{2}}}},\end{eqnarray} \tag{ 28 }$$

where ${{t}_{3}}={{T}_{b}}-{{T}_{c}}$, ${{t}_{2}}={{T}_{c}}-{{T}_{a}}$, ${{t}_{1}}={{T}_{a}}-{{T}_{d}}$.

$$\begin{eqnarray}&&D_{LAP}^{(4)}({{t}_{1}},{{t}_{2}},{{t}_{3}};{{\omega }_{a}},{{\omega }_{b}},{{\omega }_{d}})=\theta ({{t}_{1}})\theta ({{t}_{2}})\theta ({{t}_{3}}){{e}^{-i({{\omega }_{a}}+{{\omega }_{b}}-{{\omega }_{d}}){{t}_{3}}-i{{\omega }_{a}}{{t}_{1}}-i({{\omega }_{a}}+{{\omega }_{b}}){{t}_{2}}}},\end{eqnarray} \tag{ 29 }$$

where ${{t}_{3}}={{T}_{c}}-{{T}_{d}}$, ${{t}_{2}}={{T}_{d}}-{{T}_{b}}$, ${{t}_{1}}={{T}_{b}}-{{T}_{a}}$.

$$\begin{eqnarray}&&D_{LAP}^{(5)}({{t}_{1}},{{t}_{2}},{{t}_{3}};{{\omega }_{a}},{{\omega }_{b}},{{\omega }_{d}})=\theta ({{t}_{1}})\theta ({{t}_{2}})\theta ({{t}_{3}}){{e}^{-i({{\omega }_{a}}+{{\omega }_{b}}-{{\omega }_{d}}){{t}_{3}}+i{{\omega }_{d}}{{t}_{1}}+i({{\omega }_{d}}-{{\omega }_{a}}){{t}_{2}}}},\end{eqnarray} \tag{ 30 }$$

where ${{t}_{3}}={{T}_{c}}-{{T}_{b}}$, ${{t}_{2}}={{T}_{b}}-{{T}_{a}}$, ${{t}_{1}}={{T}_{a}}-{{T}_{d}}$.

$$\begin{eqnarray}&&D_{LAP}^{(6)}({{t}_{1}},{{t}_{2}},{{t}_{3}};{{\omega }_{a}},{{\omega }_{b}},{{\omega }_{d}})=\theta ({{t}_{1}})\theta ({{t}_{2}})\theta ({{t}_{3}}){{e}^{-i({{\omega }_{a}}+{{\omega }_{b}}-{{\omega }_{d}}){{t}_{3}}-i{{\omega }_{a}}{{t}_{1}}+i({{\omega }_{d}}-{{\omega }_{a}}){{t}_{2}}}},\end{eqnarray} \tag{ 31 }$$

where ${{t}_{3}}={{T}_{c}}-{{T}_{b}}$, ${{t}_{2}}={{T}_{b}}-{{T}_{d}}$, ${{t}_{1}}={{T}_{d}}-{{T}_{a}}$. The corresponding time-domain signals (25) for entangled photons are given in appendix B.

We now take the Fourier transform with respect to ladder delay variable ${{t}_{j}}$, $j=1,2,3$

$$\begin{eqnarray}&&S_{LAP}^{(j)}({{\overset{\sim}{\mathop{\Omega }}\,}_{1}},{{\overset{\sim}{\mathop{\Omega }}\,}_{2}},{{\overset{\sim}{\mathop{\Omega }}\,}_{3}})=\mathop{\int }_{-\infty }^{\infty }d{{t}_{1}}\mathop{\int }_{-\infty }^{\infty }d{{t}_{2}}\mathop{\int }_{-\infty }^{\infty }d{{t}_{3}}{{e}^{i{{\overset{\sim}{\mathop{\Omega }}\,}_{1}}{{t}_{1}}+i{{\overset{\sim}{\mathop{\Omega }}\,}_{2}}{{t}_{2}}+i{{\overset{\sim}{\mathop{\Omega }}\,}_{3}}{{t}_{3}}}}S_{LAP}^{(j)}({{t}_{1}},{{t}_{2}},{{t}_{3}}).\end{eqnarray} \tag{ 32 }$$

This gives

$$\begin{eqnarray}\begin{matrix} S_{LAP}^{(j)}({{\overset{\sim}{\mathop{\Omega }}\,}_{1}},{{\overset{\sim}{\mathop{\Omega }}\,}_{2}},{{\overset{\sim}{\mathop{\Omega }}\,}_{3}}) & = & \mathop{\int }_{-\infty }^{\infty }\frac{d{{\omega }_{a}}}{2\pi }\frac{d{{\omega }_{b}}}{2\pi }\frac{d{{\omega }_{d}}}{2\pi }\left[ D_{LAP}^{(j)}({{\overset{\sim}{\mathop{\Omega }}\,}_{1}},{{\overset{\sim}{\mathop{\Omega }}\,}_{2}},{{\overset{\sim}{\mathop{\Omega }}\,}_{3}};{{\omega }_{a}},{{\omega }_{b}},{{\omega }_{d}}){{\overset{\sim}{\mathop{S}}\,}^{(j)}}({{\omega }_{a}},{{\omega }_{b}},{{\omega }_{d}}) \right. \\ & & \left. -D_{LAP}^{(j)*}(-{{\overset{\sim}{\mathop{\Omega }}\,}_{1}},-{{\overset{\sim}{\mathop{\Omega }}\,}_{2}},-{{\overset{\sim}{\mathop{\Omega }}\,}_{3}};{{\omega }_{a}},{{\omega }_{b}},{{\omega }_{d}}){{\overset{\sim}{\mathop{S}}\,}^{(j)*}}({{\omega }_{a}},{{\omega }_{b}},{{\omega }_{d}}) \right], \\ &&\end{matrix}\end{eqnarray} \tag{ 33 }$$

where

$$\begin{eqnarray}\begin{matrix} D_{LAP}^{\left(1\right)}\left(\Omega_1,\Omega_2,\Omega_3,\omega_a,\omega_b,\omega_d\right)&=&\frac{1}{[\Omega_3+\omega_b+i\epsilon][\Omega_1+\omega_d+i\epsilon][\Omega_2+\omega_a+\omega_b+i\epsilon]},\notag\\ D_{LAP}^{\left(2\right)}\left(\Omega_1,\Omega_2,\Omega_3,\omega_a,\omega_b,\omega_d\right)&=&\frac{1}{[\Omega_3+\omega_b+i\epsilon][\Omega_1-\omega_a+i\epsilon][\Omega_2+\omega_d-\omega_a+i\epsilon]},\notag\\ D_{LAP}^{\left(3\right)}\left(\Omega_1,\Omega_2,\Omega_3,\omega_a,\omega_b,\omega_d\right)&=&\frac{1}{[\Omega_3+\omega_b+i\epsilon][\Omega_1+\omega_d+i\epsilon][\Omega_2+\omega_d-\omega_a+i\epsilon]},\notag\\ D_{LAP}^{\left(4\right)}\left(\Omega_1,\Omega_2,\Omega_3,\omega_a,\omega_b,\omega_d\right)&=&\frac{1}{[\Omega_3-\omega_a-\omega_b+\omega_d+i\epsilon][\Omega_1-\omega_a+i\epsilon][\Omega_2-\omega_a-\omega_b+i\epsilon]},\notag\\ D_{LAP}^{\left(5\right)}\left(\Omega_1,\Omega_2,\Omega_3,\omega_a,\omega_b,\omega_d\right)&=&\frac{1}{[\Omega_3-\omega_a-\omega_b+\omega_d+i\epsilon][\Omega_1+\omega_d+i\epsilon][\Omega_2+\omega_d-\omega_a+i\epsilon]},\notag\\ D_{LAP}^{\left(6\right)}\left(\Omega_1,\Omega_2,\Omega_3,\omega_a,\omega_b,\omega_d\right)&=&\frac{1}{[\Omega_3-\omega_a-\omega_b+\omega_d+i\epsilon][\Omega_1-\omega_a+i\epsilon][\Omega_2+\omega_d-\omega_a+i\epsilon]}. \end{matrix}\end{eqnarray} \tag{ 34 }$$

We note that the frequency variables ${{\omega }_{\alpha }},$, $\alpha =a,b,d$ in equation (34) are the same combinations that appear in the matter responses in equation (18). This means that for a classical light the signal will factorize into a product of several Green's functions with uncoupled frequency arguments ${{\overset{\sim}{\mathop{\Omega }}\,}_{j}}$, $j=1,2,3$. Of course this holds only in the absence of additional sources of correlating the variables caused by e.g. dephasing (bath) or the state of light. In the LOP, in contrast, the time correlations that result in the frequency mixing are apparent. Different frequency components ${{\omega }_{\alpha }}$, $\alpha =a,b,d$ that enter the Green's function in equation (18) interfere when convoluted with the same display function equation (21).

The state of light that enters the signal via the four-point frequency domain correlation function of the electric field in equation (17) can mix various frequency variables which arise from the coupling between the interaction times. In the following we consider the twin photon entangled state of light and compare it to the classical (coherent) state. Ideal multidimensional techniques use impulsive fields well separated in time with infinite bandwidth. However as shown in the following it is crucial to keep the finite bandwidth.

In the case of classical light, the four-point correlation function simply factorizes into a product of four electric field amplitudes

$$\begin{eqnarray}&&\langle {{E}^{\dagger }}({{\omega }_{d}}){{E}^{\dagger }}({{\omega }_{a}}+{{\omega }_{b}}-{{\omega }_{d}})E({{\omega }_{b}})E({{\omega }_{a}})\rangle ={{\mathcal{E}}^{*}}({{\omega }_{d}}){{\mathcal{E}}^{*}}({{\omega }_{a}}+{{\omega }_{b}}-{{\omega }_{d}})\mathcal{E}({{\omega }_{b}})\mathcal{E}({{\omega }_{a}}).\end{eqnarray} \tag{ 35 }$$

Note, that because classical fields do not impose correlations between various interaction times, either LOP or LAP can be used. In the following simulations we assume Lorentzian pulses and set $\mathcal{E}(\omega )=\langle E(\omega )\rangle ={{A}_{1}}/[\omega -{{\omega }_{p}}+i{{\sigma }_{p}}]$.

Twin photons are created via type-I spontaneous parametric down conversion of a classical pump pulse with frequency $2{{\omega }_{p}}$ into a pair of photons with central frequencies $\omega _{1}^{(0)}$ and $\omega _{2}^{(0)}$. For the degenerate process $\omega _{1}^{(0)}=\omega _{2}^{(0)}={{\omega }_{p}}$ the quantum state of light is given by the wave function

$$\begin{eqnarray}&&|\psi \rangle =\mathop{\int }_{-\infty }^{\infty }\frac{d{{\omega }_{1}}}{2\pi }\frac{d{{\omega }_{2}}}{2\pi }\Phi({{\omega }_{1}},{{\omega }_{2}})a_{{{\omega }_{1}}}^{\dagger }a_{{{\omega }_{2}}}^{\dagger }|0\rangle ,\end{eqnarray} \tag{ 36 }$$

where $a_{\omega }^{\dagger }$ is the photon creation operator in the frequency mode ω and $\Phi({{\omega }_{1}},{{\omega }_{2}})$ is the two-photon amplitude

$$\begin{eqnarray}\begin{matrix} \Phi({{\omega }_{1}},{{\omega }_{2}}) & = & {{A}_{p}}({{\omega }_{1}}+{{\omega }_{2}})sinc\left[ \frac{{{\omega }_{1}}-{{\omega }_{p}}}{2}{{T}_{1e}}+\frac{{{\omega }_{2}}-{{\omega }_{p}}}{2}{{T}_{2e}} \right]{{e}^{i\frac{{{\omega }_{1}}-{{\omega }_{p}}}{2}{{T}_{1e}}+i\frac{{{\omega }_{2}}-{{\omega }_{p}}}{2}{{T}_{2e}}}} & +({{T}_{1e}}\leftrightarrow {{T}_{2e}}), \\ &&\end{matrix}\end{eqnarray} \tag{ 37 }$$

where ${{A}_{p}}(\omega )={{A}_{0}}/[\omega -2{{\omega }_{p}}+i{{\sigma }_{p}}]$ is the Lorentzian classical pump pulse envelope. The variables ${{T}_{1e}}=L/{{v}_{p}}-L/{{v}_{1}}$ and ${{T}_{2e}}=L/{{v}_{p}}-L/{{v}_{2}}$ represent the time delays between the various beams acquired in the course of the propagation through the crystal with length L. Here, ${{v}_{p}}$, ${{v}_{1,2}}$ denote the group velocity of the pump pulse, or beams 1 and 2, respectively. The entanglement time ${{T}_{e}}={{T}_{2e}}-{{T}_{1e}}$ along with the pump bandwidth ${{\sigma }_{p}}$ are the two key parameters that define the degree of correlation between twin photons. The four-point field correlation function in equation (17) is now given by

$$\begin{eqnarray}&&\langle {{E}^{\dagger }}({{\omega }_{d}}){{E}^{\dagger }}({{\omega }_{a}}+{{\omega }_{b}}-{{\omega }_{d}})E({{\omega }_{b}})E({{\omega }_{a}})\rangle =\Phi{{\text{;}}^{*}}({{\omega }_{a}}+{{\omega }_{b}}-{{\omega }_{d}},{{\omega }_{d}})\Phi({{\omega }_{a}},{{\omega }_{b}}).\quad \;\end{eqnarray} \tag{ 38 }$$

It is important to note that since the four-point correlation function of the entangled twin state factorizes into a product of two two-point correlation functions of the form $\langle E({{\omega }_{b}})E({{\omega }_{a}})\rangle$ it only couples different interaction times within the bra- (${{T}_{a}}$, ${{T}_{b}}$) and within the ket- (${{T}_{c}}$, ${{T}_{d}}$). This means that the coupling between the interaction times in this case occurs on one branch of the loop and interaction times on different branches are not coupled. LOP thus offers a natural scanning protocol for quantum spectroscopy with entangled twin-state of light.

We have simulated the signal (20) using the LOP protocol and compared it with the standard fully time ordered LAP protocol given by equation (32) for a model trimer described by the Frenkel exciton Hamiltonian (2). We first present LOP results for classical and entangled light. We then provide the reference by demonstrating the LAP results. The parameters used are ${{\epsilon }_{1}}=1.518\;\text{eV}$, ${{\epsilon }_{2}}=1.530\;\text{eV}$, ${{\epsilon }_{3}}=1.526\;\text{eV}$, ${{J}_{12}}=10\;\text{meV}$, ${{J}_{13}}=2\;\text{meV}$, ${{J}_{13}}=3\;\text{meV}$. All three chromophores have the same transition dipole ${{V}_{1}}={{V}_{2}}={{V}_{3}}$. In the eigenstate basis the Hamiltonian (4) has parameters $\hbar {{\omega }_{{{e}_{1}}}}=1.512\;\text{eV}$, $\hbar {{\omega }_{{{e}_{2}}}}=1.525\;\text{eV}$, $\hbar {{\omega }_{{{e}_{3}}}}=1.537\;\text{eV}$, $\hbar {{\omega }_{{{f}_{1}}}}=3.044\;\text{eV}$, $\hbar {{\omega }_{{{f}_{2}}}}=3.048\;\text{eV}$, and $\hbar {{\omega }_{{{f}_{3}}}}=3.056\;\text{eV}$. We focus on the three exciton bands: g, e, and f. In our model we have two sources of dephasing. First intraband dephasing which is associated with transition within excited state band, e.g. $e-e$. Second interband dephasing that governs the transitions between e.g single and double excited states $e-f$.

5.1. LOP signals

Below we present two-dimensional signals obtained by setting one time interval to zero. Figure 3 shows the simulated ${{S}_{LOP}}({{\Omega }_{1}},{{\tau }_{2}}=0,{{\Omega }_{3}})$ for a trimer using classical light (top row) and entangled light (mid and bottom row). This signal reveals the intraband dephasing rate ${{\gamma }_{e{{e}^{\prime }}}}$ that enters through the resonance ${{\Omega }_{1}}-{{\Omega }_{3}}={{\omega }_{e{{e}^{\prime }}}}-i{{\gamma }_{e{{e}^{\prime }}}}$. We indicate the corresponding states rather than transitions, since in the loop all the transitions are calculated using the ground state as reference, thus ${{e}_{j}}\to {{e}_{j}}g$. This follows from the bookkeeping of the wavefunction. We first discuss the left column for which we set the dephasing rate ${{\gamma }_{e{{e}^{\prime }}}}=1\;\text{meV}$. Figure 3(a) shows the result for a classical light with narrow intraband dephasing ${{\gamma }_{e{{e}^{\prime }}}}=1\;\text{meV}$. It gives a diagonal cross peak $e={{e}^{\prime }}$ and one pair of weak side peaks parallel to the main diagonal at $(e,{{e}^{\prime }})=({{e}_{2}},{{e}_{3}})$. The remaining two pairs of side peaks at $(e,{{e}^{\prime }})=({{e}_{1}},{{e}_{2}})$ and $(e,{{e}^{\prime }})=({{e}_{1}},{{e}_{3}})$ are too weak to be seen. Figure 3(d) shows the signal obtained using entangled photons with short entanglement time ${{T}_{e}}=10$ fs. Panels a and d are very similar. The situation changes as the entanglement time is increased further. For ${{T}_{e}}=100$ fs in figure 3(g) we observe two additional strong side cross peak pairs with $(e,{{e}^{\prime }})=({{e}_{1}},{{e}_{3}})$ and $(e,{{e}^{\prime }})=({{e}_{1}},{{e}_{2}})$. The weak peak at $(e,{{e}^{\prime }})=({{e}_{2}},{{e}_{3}})$ is significantly enhanced as well. Note that the visibility and intensity of the side peaks is enhanced for longer entanglement time. This can be explained as follows: the long entanglement time together with the broad pump bandwidth ${{\sigma }_{p}}$ defines a parameter regime where the entanglement manifests with positive frequency correlation, i.e. the difference between frequencies of entangled photons has a narrow distribution [48]. In this case the narrow resonance occurs for ${{\Omega }_{1}}-{{\Omega }_{3}}$ and the inverse of the entanglement time is an effective bandwidth of the pulse envelope which oscillates as a function of frequency (sinc-function). The oscillating envelope enhances or suppresses certain peaks and the longer entanglement time provides the narrow bandwidth which implies a higher frequency resolution. The other two columns in figure 3 repeat these calculations for larger dephasing rates ${{\gamma }_{e{{e}^{\prime }}}}$. If the intraband dephasing ${{\gamma }_{e{{e}^{\prime }}}}$ is broader then the classical result depicted in figure 3(b) shows broadening of the main $e={{e}^{\prime }}$ peak and the side peaks are significantly suppressed compared to those shown in figure 3(a). Further increase of ${{\gamma }_{e{{e}^{\prime }}}}$ and use of classical fields leads to further broadening of the main diagonal peak whereas the side peaks completely disappear (see figure 3(c)). The same argument applies to the entangled fields with short dephasing time shown in figures 3(e)–(f). Broader dephasing rate covers the side peaks and only the main diagonal peak $e={{e}^{\prime }}$ remains strong and broad. For long entanglement time, intraband dephasing leads to broadening and enhancement of the side peaks. For instance in figure 3(g) the side peaks at $(e,{{e}^{\prime }})=({{e}_{1}},{{e}_{3}})$ are quite weak. Same peaks are broadened and enhanced in figure 3(h) and even more so in figure 3(i). Thus, the display $({{\Omega }_{1}},{{\Omega }_{3}})$ in LOP allows for effective determining of the intraband dephasing for distinct pair of e and $e\prime$ states even if intraband dephasing is broad. The advantage of having cross peaks compared to diagonal resonances is that they allow to distinguish individual states even if ${{\omega }_{e{{e}^{\prime }}}}$ is degenerate for several pairs of states e and $e\prime$. If interband dephasing ${{\gamma }_{eg}}$ which determines the longitudinal dimension of the cross peak is broad, the cross-peaks will remain distinct if properly engineered entangled light is used for probing these states. Note that the above parameter regime is different from the one studied in [34] where a narrow pump bandwidth and short entanglement time give rise to negative frequency correlations and a narrow sum frequency resonance [48]. That regime will be discussed in section 6.

Zoom In Zoom Out Reset image size

We now turn to interband dephasing. The LOP allows to extract the detailed information about ${{\gamma }_{fe}}$. Figure 4 depicts ${{S}_{LOP}}({{\Omega }_{1}},{{\Omega }_{2}},{{\tau }_{3}}=0)$. Figure 4(a) shows the signal using classical light at narrow dephasing rate ${{\gamma }_{fe}}=1\;\text{meV}$. The spectra are dominated by the resonance ${{\Omega }_{2}}-{{\Omega }_{1}}={{\omega }_{fe}}-i{{\gamma }_{fe}}$. There are total nine possible transitions between three states ${{e}_{j}}\to {{f}_{k}}$, $j,k=1,2,3$. For the small dephasing rate as in figure 4(a) one can resolve individual cross peaks and extract the information about the interband dephasing. As the dephasing rate is increased, excitation by classical light does not allow to resolve individual transitions but one can rather see only well resolved group of peaks as per figure 4(b). Further increase the dephasing rate makes the spectra broad and poorly resolved (see figure 4(c)). The short entanglement time used here provides extra selectivity over the distribution of double-excited states via ${{\Omega }_{2}}$ as follows from figure 4(d). Unlike the classical case where selectivity over ${{\Omega }_{2}}$ and ${{\Omega }_{1}}$ is the same and is determined by the interband dephasing ${{\gamma }_{eg}}\sim {{\gamma }_{fg}}$, in the entangled case, the time constraint due to ${{T}_{e}}$ provides better selectivity over ${{\Omega }_{2}}$. As the dephasing rate is increased (figure 4(e)) the ${{\Omega }_{1}}$ resolution decreases similarly to the classical case whereas the selectivity over ${{\Omega }_{2}}$ remains fixed. The same tendency holds if the dephasing is further increased as per figure 4(f). This allows to resolve individual quantum pathway that contain a single f and single e state and the dephasing ${{\gamma }_{fe}}$. Note, that the resolution of ${{\Omega }_{2}}$ is eroded for the longer entanglement time. Therefore the selectivity in both ${{\Omega }_{1}}$ and ${{\Omega }_{2}}$ is eroded quite rapidly with increase of ${{\gamma }_{fe}}$ as illustrated in figures 4(g)–(i).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

5.2. LAP signals

As we did for the LOP we show 2D signals obtained by setting one time interval to zero. Figure 5 depicts the LAP signal ${{S}_{LAP}}({{\overset{\sim}{\mathop{\Omega }}\,}_{1}},{{t}_{2}}=0,2{{\omega }_{p}}-{{\overset{\sim}{\mathop{\Omega }}\,}_{3}})$ (32) (we plotted it versus $2{{\omega }_{p}}-{{\overset{\sim}{\mathop{\Omega }}\,}_{3}}$ for a better comparison with figure 3). As we did for the LOP we investigate the effect of intraband dephasing ${{\gamma }_{e{{e}^{\prime }}}}$. We set ${{\gamma }_{e{{e}^{\prime }}}}=1\;\text{meV}$. Unlike the LOP which contains contributions from all six diagrams in figure 2, LAP allows to distinguish between the ${{\mathbf{k}}_{I}}$, ${{\mathbf{k}}_{II}}$ and ${{\mathbf{k}}_{III}}$ contributions. Figure 5(a) shows the ${{\mathbf{k}}_{I}}$ signal for the narrow intraband dephasing ${{\gamma }_{e{{e}^{\prime }}}}=1\;\text{meV}$. It is dominated by two resonances for ${{\overset{\sim}{\mathop{\Omega }}\,}_{1}}\simeq {{\omega }_{eg}}$, $2{{\omega }_{p}}-{{\omega }_{fe}}$, and ${{\overset{\sim}{\mathop{\Omega }}\,}_{3}}\simeq {{\omega }_{fe}}$, $2{{\omega }_{p}}-{{\omega }_{eg}}$. The ${{\mathbf{k}}_{II}}$ signal shown in figure 5(b) is dominated by a cross-peak at ${{\overset{\sim}{\mathop{\Omega }}\,}_{1}}+{{\overset{\sim}{\mathop{\Omega }}\,}_{3}}=2{{\omega }_{p}}-{{\omega }_{e{{e}^{\prime }}}}-i({{\sigma }_{p}}+{{\gamma }_{e{{e}^{\prime }}}})$. Note, that unlike the LOP, in the case of LAP the width of the $ee\prime$ resonance is affected by the pump pulse bandwidth ${{\sigma }_{p}}$ and the resonance is broadened as can be seen from figure 5(b). The same applies to ${{\mathbf{k}}_{III}}$. For comparison with the LOP we plot the sum of all three techniques in figure 5(d). It resembles ${{\mathbf{k}}_{II}}$ and ${{\mathbf{k}}_{III}}$ and shows that for the same parameters compare to LOP, we get significantly broader resonances and thus, information about intraband dephasing cannot be effectively extracted from this display mode. As shown below it can be done from the $({{t}_{1}}=0,{{\overset{\sim}{\mathop{\Omega }}\,}_{2}},{{\overset{\sim}{\mathop{\Omega }}\,}_{3}})$ display. Unlike the LOP where entanglement at long times ${{T}_{e}}$ plays a crucial role, LAP does not carry extra information about intraband dephasing and essentially gives similar spectra to classical light. Slight changes in peaks intensities can be observed at long entanglement times in ${{\mathbf{k}}_{II}}$ and ${{\mathbf{k}}_{III}}$ signals (see figures 5(j), (k)) compared to short entanglement time in figures 5(f), (g) and classical light in figures 5(b), (c).

As demonstrated above, displaying the LAP signal vs $({{\overset{\sim}{\mathop{\Omega }}\,}_{1}},{{t}_{2}}=0,{{\overset{\sim}{\mathop{\Omega }}\,}_{3}})$ does not allow to extract the intraband dephasing ${{\gamma }_{e{{e}^{\prime }}}}$ since the spectra are broadened by the pulse bandwidth (see figure 5). However we can extract the intraband dephasing by plotting ${{\mathbf{k}}_{I}}$ signal (equation (22)) and ${{\mathbf{k}}_{II}}$ (equation (23)) if displayed versus $({{t}_{1}}=0,{{\overset{\sim}{\mathop{\Omega }}\,}_{2}},{{\overset{\sim}{\mathop{\Omega }}\,}_{3}})$—see figure 6. Note, that here we depicted the ticks along the axes corresponding to the relevant transitions keeping track of the density matrix. The spectra are dominated by ${{\overset{\sim}{\mathop{\Omega }}\,}_{2}}={{\omega }_{e{{e}^{\prime }}}}-i{{\gamma }_{e{{e}^{\prime }}}}$ resonance. The ${{\mathbf{k}}_{III}}$ signal does not show any features in the vicinity of ${{\overset{\sim}{\mathop{\Omega }}\,}_{2}}={{\omega }_{e{{e}^{\prime }}}}$ so it is not shown. For the narrow dephasing ${{\gamma }_{e{{e}^{\prime }}}}$ the spectra of ${{\mathbf{k}}_{I}}$ signal shown in figure 6(a) shows strong diagonal $e={{e}^{\prime }}$ resonance and weak cross peaks at $(e,{{e}^{\prime }})=({{e}_{1}},{{e}_{2}})$ and $(e,{{e}^{\prime }})=({{e}_{2}},{{e}_{3}})$. The peak at $(e,{{e}^{\prime }})=({{e}_{1}},{{e}_{3}})$ is significantly weaker than the other two. Similar spectra is obtained for ${{\mathbf{k}}_{II}}$ signal (figure 8(b)) and the total ${{\mathbf{k}}_{I}}+{{\mathbf{k}}_{II}}$ signal—figure 6(c). Using entangled light with short entanglement time ${{T}_{e}}=10$ fs, the spectra are virtually identical to the classical light as shown figures 6(d)–(f) compared to figures 6(a)–(c). The interesting effect occurs for the long entanglement time as in the case of LOP. Figure 6(g) shows the side peaks $e\ne {{e}^{\prime }}$ in rephasing signal ${{\mathbf{k}}_{I}}$ are suppressed, whereas the nonrephasing contribution ${{\mathbf{k}}_{II}}$ in figure 6(h) has enhanced side peaks including $(e,{{e}^{\prime }})=({{e}_{1}},{{e}_{3}})$ resonance that becomes well pronounced. The total ${{\mathbf{k}}_{I}}+{{\mathbf{k}}_{II}}$ signal depicted in figure 6(i) shows the suppressed side resonances.

For comparison with the LOP and with figure 4 that reveals interband dephasing ${{\gamma }_{fe}}$ we plot the LAP signal versus $({{\overset{\sim}{\mathop{\Omega }}\,}_{2}},{{\overset{\sim}{\mathop{\Omega }}\,}_{3}})$ in figure 7. For a narrow dephasing ${{\gamma }_{fe}}=1\;\text{meV}$ the spectra reveals nine ${{\overset{\sim}{\mathop{\Omega }}\,}_{3}}={{\omega }_{fe}}-i{{\gamma }_{fe}}$ peaks as shown in figure 7(a). For broader dephasing ${{\gamma }_{fe}}=2\;\text{meV}$ (figure 7(b)) and ${{\gamma }_{fe}}=4\;\text{meV}$ (figure 7(c)) the spectra are broadened and various peaks overlap. When using entangled light the short entanglement case with ${{T}_{e}}=10$ fs is very similar to the classical light as can be seen by comparing figures 7(d)–(f) with figures 7(a)–(c). Unlike spectra in figures 3–6, longer entanglement time does not provide any benefit. Rather it makes various peaks suppressed compared to the classical case as can be seen in figures 7(g)–(i).

Zoom In Zoom Out Reset image size

It follows from figures 3 and 5 that entanglement is not necessary to reveal the narrow intraband dephasing ${{\gamma }_{e{{e}^{\prime }}}}$. The narrow resonances can be observed with classical light if displayed using the LOP as shown in figure 3(d) with the proper choice of field parameters. The LAP on the other hand cannot reveal the narrow dephasing irrespective of entanglement, as shown in figure 5.

2D spectra of a model dimer with classical and entangled light were calculated recently in [39]. Narrow resonances found in that work were attributed to photon entanglement. However our analysis shows that these resonances are solely connected with the scanning protocol and are unrelated to entanglement. Figure 8 displays the signals calculated using our approach for the same model dimer parameters of [39]. The LOP spectra for classical and entangled light are compared in the left column. The corresponding LAP spectra are shown in the right column. We see that entanglement makes no difference in this parameter regime (the two rows are virtually identical). However the scanning protocol does as seen by the two columns. The LOP signals are narrow and clearly resolve the ${{e}_{1}}$ and ${{e}_{2}}$ states whereas the corresponding LAP signals are broad and featureless.

Zoom In Zoom Out Reset image size

For more in depth comparison we now describe the signals calculated in [39] using our terminology. In that work the entangled LOP spectrum (bottom row of their figure 8, our figure 8(c)) was compared with the classical LAP spectrum (bottom row of their figure 6 corresponding to our figure 8(b)). In [39] the difference was attributed to entanglement effects. Our results show that the difference is solely due to the different scanning protocol (LAP/LOP) and is unrelated to entanglement. Note that the LAP yields three different signals that can be distinguished by the choice of phase, whereas the LOP combines all six contributions into one signal. Furthermore, in order to recover the expressions in equations (21)–(23) and (24)–(29) of [39] using our model of entangled light we had to take continuous limit for entangled case ${{\sigma }_{p}}\to 0$ with ${{T}_{e1}}={{T}_{e2}}=0$ and the impulsive limit ${{\sigma }_{p}}\to \infty$ in the case of classical light, which corresponds to two completely different parameter regimes. For a consistent comparison of the classical vs entangled light we used in figure 8 the impulsive limit for all four signals.

Zoom In Zoom Out Reset image size

The LOP protocol can be generally realized using a pulse shaper as explained earlier. Using our analysis we conclude that the Franson interferometer proposed in [39] shall provide a convenient method for realizing this protocol experimentally.

So far we investigated multidimensional signals obtained by scanning various time delays between pulses. This time-domain protocol makes sense if the pulses that interact with the system are relatively short. For entangled light this implies that frequencies of the modes corresponding to the twin photons are positively correlated [49]. We demonstrated that this is crucial especially for the long entanglement time where the narrow difference-frequency-resonances can be observed in the spectra of ${{S}_{LOP}}({{\Omega }_{1}},{{\tau }_{2}}=0,{{\Omega }_{3}})$, ${{S}_{LAP}}({{\overset{\sim}{\mathop{\Omega }}\,}_{1}},{{t}_{2}}=0,{{\overset{\sim}{\mathop{\Omega }}\,}_{3}})$, ${{S}_{LOP}}({{\Omega }_{1}},{{\Omega }_{2}},{{\tau }_{3}}=0)$, and ${{S}_{LAP}}({{t}_{1}}=0,{{\overset{\sim}{\mathop{\Omega }}\,}_{2}},{{\overset{\sim}{\mathop{\Omega }}\,}_{3}})$ signals.

In our recent work [34] we have investigated the effects of entanglement on the control of the transport properties in molecular aggregates. Narrow fg resonances were observed when the entangled pair has been generated by narrowband pump and entanglement time is short. In this case narrow pump along with short entanglement time implies negative frequency correlation (the sum of two frequencies is narrowly distributed). This is a different parameter regime than used in section 5. In the following we consider narrowband pump pulse and fg resonances with entangled photons. In this case we can adopt mixed time-and-frequency domain scanning, where we scan one time delay between pulses and the pump frequency ${{\omega }_{p}}$. Again we compare the LOP and LAP protocols.

Figure 9 depicts the corresponding time-and-frequency domain signal. Figure 9(a) shows the signal ${{S}_{LOP}}({{\Omega }_{1}},{{\Omega }_{2}}=3.11\;eV,{{\tau }_{3}}=0,{{\omega }_{p}})$ LOP signal which contains three distinct peaks corresponding to $2{{\omega }_{p}}={{\omega }_{fg}}-i{{\gamma }_{fg}}$ resonances for short entanglement time ${{T}_{e}}=10$ fs. As the entanglement time is increased, the peaks become weaker as shown in figure 9(f). For comparison we depict the corresponding series of LAP signals ${{S}_{{{\mathbf{k}}_{j}}}}({{t}_{1}}=0,{{\overset{\sim}{\mathop{\Omega }}\,}_{2}}=3.11\;eV,{{\overset{\sim}{\mathop{\Omega }}\,}_{3}};{{\omega }_{p}})$ where $j=I,II,III$. For short entanglement time the ${{\mathbf{k}}_{II}}$ signal contains three well pronounced narrow peaks similar to LOP as seen from figure 9(c), whereas resonances in ${{\mathbf{k}}_{I}}$ and ${{\mathbf{k}}_{III}}$ signals are not as clearly seen. For longer entanglement time, all sharp features of LAP spectra become fuzzy (see figures 9(g)–(j)) and even in the case of ${{\mathbf{k}}_{II}}$, ${{\omega }_{fg}}$ resonances become suppressed. This is consistent with earlier results for narrowband pump pulse [34].

Zoom In Zoom Out Reset image size

Multidimensional optical signals are obtained by subjecting the system to sequences of short pulses and generating and analyzing correlation plots between different resonances generated during controlled delay periods. These allow to visualize such an event as a e.g. cross-peak in the space of two frequency variables that are related to Fourier transform of two different delay intervals. Most commonly, the delays are between consecutive chronologically-ordered pulses that can differ by their frequencies, polarizations and wavevectors. Such signals can be naturally described by the density matrix and represented by ladder diagrams. We had presented a new protocol based on the wavefunction description that involves both forward and backward time evolution. This protocol uses different types of delays represented by loop diagrams and can be realized experimentally by phase cycling. This new type of bookkeeping of field-matter interactions that is not based on chronologically time ordered pulses suggests a new way of monitoring and displaying various resonances. We demonstrated it for two photon absorption experiments with incoherent fluorescence detection in a molecular aggregate with classical and entangled light.

Broadband entangled light with long entanglement time allows to selectively reduce the background and reveal certain resonances because of intrinsic frequency correlations due to entanglement. The resonances remain well resolved even for the short dephasing which typically is a source of strong background for the signals measured with classical fields. In particular, entangled light and the loop-based protocol can reveal intra and interband dephasing in the single and double exciton manifold not possible by classical light. We demonstrated better-resolved signals compared to those obtained with standard ladder scanning protocol. Entangled light causes correlations of the various time delay variables thus providing new spectroscopic windows and physical picture of the system dynamics. The current formalism can be readily applied for an arbitrary state of light including stochastic, squeezed or other quantum and classical states. The signals are given by sums of products of four-point correlation functions of the electric field and matter which can be calculated for arbitrary pulse shapes and bandwidths including temporally overlapping pulses. The necessary Liouville space Green's functions can be evaluated by taking bath effects into account, e.g. pure dephasing, inhomogeneous broadening, transport and other dynamical bath effects.

We wish to thank Frank Schlawin for stimulating discussions and gratefully acknowledge the support of the National Science Foundation through grant no. CHE-1058791 and computations are supported by CHE-0840513, the Chemical Sciences, Geosciences and Biosciences Division, Office of Basic Energy Sciences, Office of Science, US Department of Energy and the National Institute of Health grant no. GM-59230.

Here we evaluate the frequency integrals in equation (15). The time-domain LOP signal then reads

$$\begin{eqnarray} S_{LOP}^{\left(1\right)}\left(\tau_1,\tau_2,\tau_3\right)&=&\mathcal{I}\frac{1}{2\sigma_p\hbar^4}\sum_{e,e^{\prime}}\mu_{ge^{\prime}}^a\mu_{e^{\prime}f}^b\mu_{fe}^{c^{\ast}}\mu_{eg}^{d^{\ast}}\theta\left(\tau_2\right)\mathcal{G}_{ef}\left(\tau_1\right)\mathcal{G}_{ge^{\prime}}\left(\tau_3\right)\\ &&\times\left[\theta\left(\tau_2-\tau_1\right)e^{[2i\omega_p-\sigma_p]\left(\tau_2-\tau_1\right)}\mathcal{G}_{gf}\left(-2\omega_p-i\sigma_p\right)\mathcal{F}\left(2\omega_p-\omega_{fe}-i\left(\gamma_{fe}-\sigma_p\right),\omega_{fe}+i\gamma_{fe},\omega_{e^{\prime}g}+i\gamma_{e^{\prime}g}\right)\right.\\ &&\left.+\theta\left(\tau_1-\tau_2\right)e^{-[2i\omega_p+\sigma_p]\left(\tau_1-\tau_2\right)}\mathcal{G}_{gf}\left(-2\omega_p+i\sigma_p\right)\mathcal{F}\left(2\omega_p-\omega_{fe}-i\left(\gamma_{fe}+\sigma_p\right),\omega_{fe}+i\gamma_{fe},\omega_{e^{\prime}g}+i\gamma_{e^{\prime}g}\right)\right.\\ &&\left.+2\sigma_p\mathcal{G}_{gf}\left(\tau_2-\tau_1\right)\mathcal{G}_{gf}\left(-2\omega_p-i\sigma_p\right)\mathcal{G}_{gf}\left(-2\omega_p+i\sigma_p\right)\mathcal{F}\left(\omega_{eg}-i\gamma_{eg},\omega_{fe}+i\gamma_{fe},\omega_{e^{\prime}g}+i\gamma_{e^{\prime}g}\right)\right], \end{eqnarray} \tag{ A.1 }$$

$$\begin{eqnarray}\begin{matrix} S_{LOP}^{\left(2\right)}\left(\tau_1,\tau_2,\tau_3\right)&=&-\mathcal{I}\frac{1}{2\sigma_p\hbar^4}\sum_{e,e^{\prime}}\mu_{ge^{\prime}}^a\mu_{e^{\prime}f}^b\mu_{fe}^{c^{\ast}}\mu_{eg}^{d^{\ast}}\theta\left(\tau_1\right)\theta\left(\tau_2\right)\mathcal{G}_{ee^{\prime}}\left(\tau_3\right)\\ &&\times\left[\theta\left(\tau_2-\tau_1+\tau_3\right)e^{[2i\omega_p-\sigma_p]\left(\tau_2-\tau_1+\tau_3\right)}\mathcal{G}_{ef}\left(\tau_1-\tau_3\right)\mathcal{G}_{gf}\left(-2\omega_p-i\sigma_p\right)\right.\\ &&\left.\times\mathcal{F}\left(2\omega_p-\omega_{fe}-i\left(\gamma_{fe}-\sigma_p\right),\omega_{fe}+i\gamma_{fe},2\omega_p-\omega_{fe^{\prime}}-i\left(\gamma_{fg}-\gamma_{e^{\prime}g}-\sigma_p\right)\right)\right.\\ &&\left.+\theta\left(\tau_1-\tau_3-\tau_2\right)e^{-[2i\omega_p+\sigma_p]\left(\tau_1-\tau_3+\tau_2\right)}\mathcal{G}_{ef}\left(\tau_1-\tau_3\right)\mathcal{G}_{gf}\left(-2\omega_p+i\sigma_p\right)\right.\\ &&\left.\times\mathcal{F}\left(2\omega_p-\omega_{fe}-i\left(\gamma_{fe}+\sigma_p\right),\omega_{fe}+i\gamma_{fe},2\omega_p-\omega_{fe^{\prime}}-i\left(\gamma_{fg}-\gamma_{e^{\prime}g}+\sigma_p\right)\right)\right.\\ &&\left.-\theta\left(\tau_2\right)e^{[2i\omega_p-\sigma_p]\tau_2}\mathcal{G}_{eg}\left(\tau_1-\tau_3\right)\mathcal{G}_{gf}\left(-2\omega_p-i\sigma_p\right)\right.\\ &&\left.\times\mathcal{F}\left(\omega_{eg}-i\gamma_{eg},2\omega_p-\omega_{eg}+i\left(\gamma_{eg}+\sigma_p\right),\omega_{e^{\prime}g}+i\gamma_{e^{\prime}g}\right)\right.\\ &&\left.+2\sigma_p[\mathcal{G}_{ef}\left(\tau_1-\tau_3\right)\mathcal{G}_{gf}\left(\tau_2-\tau_1+\tau_3\right)-\mathcal{G}_{eg}\left(\tau_1-\tau_3\right)\mathcal{G}_{gf}\left(\tau_2\right)]\right.\\ &&\left.\times\mathcal{G}_{gf}\left(-2\omega_p-i\sigma_p\right)\mathcal{G}_{gf}\left(-2\omega_p+i\sigma_p\right)\mathcal{F}\left(\omega_{eg}-i\gamma_{eg},\omega_{fe}+i\gamma_{fe},\omega_{e^{\prime}g}+i\gamma_{e^{\prime}g}\right)\right], \end{matrix}\end{eqnarray} \tag{ A.2 }$$

$$\begin{eqnarray}\begin{matrix} S_{LOP}^{\left(3\right)}\left(\tau_1,\tau_2,\tau_3\right)&=&\mathcal{I}\frac{1}{2\sigma_p\hbar^4}\sum_{e,e^{\prime}}\mu_{ge^{\prime}}^a\mu_{e^{\prime}f}^b\mu_{fe}^{c^{\ast}}\mu_{eg}^{d^{\ast}}\theta\left(\tau_1\right)\theta\left(\tau_2\right)\theta\left(\tau_3\right)\\ &&\times\left[-\theta\left(\tau_1-\tau_2\right)e^{-[2i\omega_p+\sigma_p]\left(\tau_1-\tau_2\right)}\mathcal{G}_{ef}\left(\tau_1\right)\mathcal{G}_{ge^{\prime}}\left(\tau_3\right)\mathcal{G}_{gf}\left(-2\omega_p+i\sigma_p\right)\right.\\ &&\left.\times\mathcal{F}\left(2\omega_p-\omega_{fe}-i\left(\gamma_{fe}+\sigma_p\right),\omega_{fe}+i\gamma_{fe},\omega_{e^{\prime}g}+i\gamma_{e^{\prime}g}\right)\right.\\ &&\left.-\theta\left(\tau_2-\tau_1\right)e^{[2i\omega_p-\sigma_p]\left(\tau_2-\tau_1\right)}\mathcal{G}_{ef}\left(\tau_1\right)\mathcal{G}_{ge^{\prime}}\left(\tau_3\right)\mathcal{G}_{gf}\left(-2\omega_p-i\sigma_p\right)\right.\\ &&\left.\times\mathcal{F}\left(2\omega_p-\omega_{fe}-i\left(\gamma_{fe}-\sigma_p\right),\omega_{fe}+i\gamma_{fe},\omega_{e^{\prime}g}+i\gamma_{e^{\prime}g}\right)\right.\\ &&\left.-2\sigma_p\mathcal{G}_{gf}\left(\tau_2-\tau_1\right)\mathcal{G}_{ef}\left(\tau_1\right)\mathcal{G}_{ge^{\prime}}\left(\tau_3\right)\mathcal{G}_{gf}\left(-2\omega_p-i\sigma_p\right)\mathcal{G}_{gf}\left(-2\omega_p+i\sigma_p\right)\right.\\ &&\left.\times\mathcal{F}\left(\omega_{eg}-i\gamma_{eg},\omega_{fe}+i\gamma_{fe},\omega_{e^{\prime}g}+i\gamma_{e^{\prime}g}\right)\right.\\ &&\left.+\theta\left(\tau_2\right)e^{[2i\omega_p-\sigma_p]\tau_2}\mathcal{G}_{eg}\left(\tau_1\right)\mathcal{G}_{ge^{\prime}}\left(\tau_3\right)\mathcal{G}_{gf}\left(-2\omega_p-i\sigma_p\right)\right.\\ &&\left.\times\mathcal{F}\left(\omega_{eg}-i\gamma_{eg},2\omega_p-\omega_{eg}+i\left(\sigma_p+\gamma_{eg}\right),\omega_{e^{\prime}g}+i\gamma_{e^{\prime}g}\right)\right.\\ &&\left.+2\sigma_p\mathcal{G}_{gf}\left(\tau_2\right)\mathcal{G}_{eg}\left(\tau_1\right)\mathcal{G}_{ge^{\prime}}\left(\tau_3\right)\mathcal{G}_{gf}\left(-2\omega_p-i\sigma_p\right)\mathcal{G}_{gf}\left(-2\omega_p+i\sigma_p\right)\right.\\ &&\left.\times\mathcal{F}\left(\omega_{eg}-i\gamma_{eg},\omega_{fe}+i\gamma_{fe},\omega_{e^{\prime}g}+i\gamma_{e^{\prime}g}\right)\right.\\ &&\left.+\theta\left(\tau_1-\tau_3-\tau_2\right)e^{-[2i\omega_p+\sigma_p]\left(\tau_1-\tau_3-\tau_2\right)}\mathcal{G}_{ef}\left(\tau_1-\tau_3\right)\mathcal{G}_{ee^{\prime}}\left(\tau_3\right)\mathcal{G}_{gf}\left(-2\omega_p+i\sigma_p\right)\right.\\ &&\left.\times\mathcal{F}\left(2\omega_p-\omega_{fe}-i\left(\gamma_{fe}+\sigma_p\right),\omega_{fe}+i\gamma_{fe},2\omega_p-\omega_{fe^{\prime}}-i\left(\gamma_{fg}-\gamma_{e^{\prime}g}+\sigma_p\right)\right)\right.\\ &&\left.+\theta\left(\tau_2+\tau_3-\tau_1\right)e^{[2i\omega_p-\sigma_p]\left(\tau_2+\tau_3-\tau_1\right)}\mathcal{G}_{ef}\left(\tau_1-\tau_3\right)\mathcal{G}_{ee^{\prime}}\left(\tau_3\right)\mathcal{G}_{gf}\left(-2\omega_p-i\sigma_p\right)\right.\\ &&\left.\times\mathcal{F}\left(2\omega_p-\omega_{fe}-i\left(\gamma_{fe}-\sigma_p\right),\omega_{fe}+i\gamma_{fe},2\omega_p-\omega_{fe^{\prime}}-i\left(\gamma_{fg}-\gamma_{e^{\prime}g}-\sigma_p\right)\right)\right.\\ &&\left.+2\sigma_p\mathcal{G}_{gf}\left(\tau_2+\tau_3-\tau_1\right)\mathcal{G}_{ef}\left(\tau_1-\tau_3\right)\mathcal{G}_{ee^{\prime}}\left(\tau_3\right)\mathcal{G}_{gf}\left(-2\omega_p-i\sigma_p\right)\mathcal{G}_{gf}\left(-2\omega_p+i\sigma_p\right)\right.\\ &&\left.\times\mathcal{F}\left(\omega_{eg}-i\gamma_{eg},\omega_{fe}+i\gamma_{fe},\omega_{e^{\prime}g}+i\gamma_{e^{\prime}g}\right)\right.\\ &&\left.-\theta\left(\tau_2\right)e^{[2i\omega_p-\sigma_p]\tau_2}\mathcal{G}_{eg}\left(\tau_1-\tau_3\right)\mathcal{G}_{ee^{\prime}}\left(\tau_3\right)\mathcal{G}_{gf}\left(-2\omega_p-i\sigma_p\right)\right.\\ &&\left.\times\mathcal{F}\left(\omega_{eg}-i\gamma_{eg},2\omega_p-\omega_{eg}+i\left(\gamma_{eg}+\sigma_p\right),\omega_{e^{\prime}g}+i\gamma_{e^{\prime}g}\right)\right.\\ &&\left.-2\sigma_p\mathcal{G}_{gf}\left(\tau_2\right)\mathcal{G}_{eg}\left(\tau_1-\tau_3\right)\mathcal{G}_{ee^{\prime}}\left(\tau_3\right)\mathcal{G}_{gf}\left(-2\omega_p-i\sigma_p\right)\mathcal{G}_{gf}\left(-2\omega_p+i\sigma_p\right)\right.\\ &&\left.\times\mathcal{F}\left(\omega_{eg}-i\gamma_{eg},\omega_{fe}+i\gamma_{fe},\omega_{e^{\prime}g}+i\gamma_{e^{\prime}g}\right)\right], \end{matrix}\end{eqnarray} \tag{ A.3 }$$

$$\begin{eqnarray}&&S_{LOP}^{(4)}({{\tau }_{1}},{{\tau }_{2}},{{\tau }_{3}})={{S}^{(1)*}}({{\tau }_{3}},{{\tau }_{2}},{{\tau }_{1}}),\end{eqnarray} \tag{ A.4 }$$

$$\begin{eqnarray}&&S_{LOP}^{(5)}({{\tau }_{1}},{{\tau }_{2}},{{\tau }_{3}})={{S}^{(2)*}}({{\tau }_{3}},{{\tau }_{2}},{{\tau }_{1}}),\end{eqnarray} \tag{ A.5 }$$

$$\begin{eqnarray}&&S_{LOP}^{(6)}({{\tau }_{1}},{{\tau }_{2}},{{\tau }_{3}})={{S}^{(3)*}}({{\tau }_{3}},{{\tau }_{2}},{{\tau }_{1}}),\end{eqnarray} \tag{ A.6 }$$

where

$$\begin{eqnarray}&&\mathcal{F}({{\omega }_{a}},{{\omega }_{b}},{{\omega }_{d}})=\overset{\sim}{\mathop{\Phi}}\,({{\omega }_{a}},{{\omega }_{b}}){{\overset{\sim}{\mathop{\Phi}}\,}^{*}}({{\omega }_{a}}+{{\omega }_{b}}-{{\omega }_{d}},{{\omega }_{d}}).\end{eqnarray} \tag{ A.7 }$$

and

$$\begin{eqnarray}&&\Phi({{\omega }_{a}},{{\omega }_{b}})={{A}_{p}}({{\omega }_{a}}+{{\omega }_{b}})\overset{\sim}{\mathop{\Phi}}\,({{\omega }_{a}},{{\omega }_{b}}).\end{eqnarray} \tag{ A.8 }$$

Evaluating the frequency integrals in (25) we obtain for the time-domain LAP signal

$$\begin{eqnarray}\begin{matrix} S_{LAP}^{(1)}({{t}_{1}},{{t}_{2}},{{t}_{3}}) & = & \mathcal{I}\frac{1}{2{{\sigma }_{p}}{{\hbar }^{4}}}\underset{e,{{e}^{\prime }}}{\mathop{\mathop{\sum }}}\,\mu _{g{{e}^{\prime }}}^{a}\mu _{{{e}^{\prime }}f}^{b}\mu _{fe}^{c*}\mu _{eg}^{d*}{{\mathcal{G}}_{ef}}({{t}_{3}}){{\mathcal{G}}_{g{{e}^{\prime }}}}({{t}_{1}}) \\ & & \times \left[ \theta ({{t}_{2}}){{e}^{[2i{{\omega }_{p}}-{{\sigma }_{p}}]{{t}_{2}}}}{{\mathcal{G}}_{gf}}(-2{{\omega }_{p}}-i{{\sigma }_{p}})\mathcal{F}(2{{\omega }_{p}}-{{\omega }_{fe}}-i({{\gamma }_{fe}}-{{\sigma }_{p}}), \right. \\ & & \times {{\omega }_{fe}}+i{{\gamma }_{fe}},{{\omega }_{{{e}^{\prime }}g}}+i{{\gamma }_{{{e}^{\prime }}g}}) \\ & & +2{{\sigma }_{p}}{{\mathcal{G}}_{gf}}({{t}_{2}}){{\mathcal{G}}_{gf}}(-2{{\omega }_{p}}-i{{\sigma }_{p}}){{\mathcal{G}}_{gf}}(-2{{\omega }_{p}}+i{{\sigma }_{p}}) \\ & & \times \left. \mathcal{F}({{\omega }_{eg}}-i{{\gamma }_{eg}},{{\omega }_{fe}}+i{{\gamma }_{fe}},{{\omega }_{{{e}^{\prime }}g}}+i{{\gamma }_{{{e}^{\prime }}g}}) \right], \\ &&\end{matrix}\end{eqnarray} \tag{ B.1 }$$

$$\begin{eqnarray}\begin{matrix} S_{LAP}^{\left(2\right)}\left(t_1,t_2,t_3\right)&=&-\mathcal{R}\frac{1}{2\sigma_p\hbar^4}\sum_{e,e^{\prime}}\mu_{ge^{\prime}}^a\mu_{e^{\prime}f}^b\mu_{fe}^{c^{\ast}}\mu_{eg}^{d^{\ast}}\mathcal{G}_{ee^{\prime}}\left(t_2\right)\\ &&\times\left[e^{[2i\omega_p-\sigma_p]t_3}\mathcal{G}_{eg}\left(t_3\right)\mathcal{G}_{eg}\left(t_1\right)\mathcal{G}_{gf}\left(-2\omega_p-i\sigma_p\right)\right.\\ &&\left.\times\mathcal{F}\left(\omega_{eg}-i\gamma_{eg},2\omega_p-\omega_{eg}+i\left(\sigma_p+\gamma_{eg}\right),\omega_{e^{\prime}g}+i\gamma_{e^{\prime}g}\right)\right.\\ &&\left.+e^{-[2i\omega_p+\sigma_p]t_1}\mathcal{G}_{ef}\left(t_3\right)\mathcal{G}_{ef}\left(t_1\right)\mathcal{G}_{gf}\left(-2\omega_p+i\sigma_p\right)\right.\\ &&\left.\times\mathcal{F}\left(2\omega_p-\omega_{fe}-i\left(\gamma_{fe}+\sigma_p\right),\omega_{fe}+i\gamma_{fe},2\omega_p-\omega_{fe^{\prime}}-i\left(\gamma_{fg}-\gamma_{e^{\prime}g}+\sigma_p\right)\right)\right.\\ &&\left.+2i\sigma_p\mathcal{G}_{ef}\left(t_3\right)\mathcal{G}_{eg}\left(t_1\right)\mathcal{G}_{gf}\left(-2\omega_p-i\sigma_p\right)\mathcal{G}_{gf}\left(-2\omega_p+i\sigma_p\right)\right.\\ &&\left.\times\mathcal{F}\left(\omega_{eg}-i\gamma_{eg},\omega_{fe}+i\gamma_{fe},\omega_{e^{\prime}g}+i\gamma_{e^{\prime}g}\right)\right], \end{matrix}\end{eqnarray} \tag{ B.2 }$$

$$\begin{eqnarray}\begin{matrix} S_{LAP}^{(3)}(t_1,t_2,t_3)&=&\mathcal{R}\frac{1}{2\sigma_p\hbar^4}\sum_{e,e^{\prime}}\mu_{ge^{\prime}}^a\mu_{e^{\prime}f}^b\mu_{fe}^{c^{\ast}}\mu_{eg}^{d^{\ast}}\theta(t_1)\theta(t_2)\theta(t_3)\\ &&\times\left[-e^{[2i\omega_p-\sigma_p]t_3}\mathcal{G}_{eg}(t_3)\mathcal{G}_{eg}(t_2)\mathcal{G}_{ge^{\prime}}(t_2+t_1)\mathcal{G}_{gf}(-2\omega_p-i\sigma_p)\right.\\ &&\left.\times\mathcal{F}(\omega_{eg}-i\gamma_{eg},2\omega_p-\omega_{eg}+i(\sigma_p+\gamma_{eg}),\omega_{e^{\prime}g}+i\gamma_{e^{\prime}g})\right.\\ &&\left.-e^{-[2i\omega_p+\Gamma]t_2}\mathcal{G}_{ef}(t_3)\mathcal{G}_{ef}(t_2)\mathcal{G}_{ge^{\prime}}(t_2+t_1))\mathcal{G}_{gf}(-2\omega_p+i\sigma_p)\right.\\ &&\left.\times\mathcal{F}(2\omega_p-\omega_{fe}-i(\gamma_{fe}+\sigma_p),\omega_{fe}+i\gamma_{fe},\omega_{e^{\prime}g}+i\gamma_{e^{\prime}g})\right.\\ &&\left.-2i\sigma_p\mathcal{G}_{ef}(t_3)\mathcal{G}_{eg}(t_2)\mathcal{G}_{ge^{\prime}}(t_2+t_1)\mathcal{G}_{gf}(-2\omega_p-i\sigma_p)\mathcal{G}_{gf}(-2\omega_p+i\sigma_p)\right.\\ &&\left.\times \mathcal{F}(\omega_{eg}-i\gamma_{eg},\omega_{fe}+i\gamma_{fe},\omega_{e^{\prime}g}+i\gamma_{e^{\prime}g})\right.\\ &&\left.+ie^{(2i\omega_p-\sigma_p]t_3}\mathcal{G}_{eg}(t_3-t_1)\mathcal{G}_{ee^{\prime}}(t_2+t_1)\mathcal{G}_{gf}(-2\omega_p-i\sigma_p)\right.\\ &&\left.\times\mathcal{F}(\omega_{eg}-i\gamma_{eg},2\omega_p-\omega_{eg}+i(\sigma_p+\gamma_{eg}),\omega_{e^{\prime}g}+i\gamma_{e^{\prime}g})\right.\\ &&\left. +\theta(t_1-t_3)e^{[2i\omega_p-\sigma_p]t_1-[i\omega_{fe}-\gamma_{fe}](t_1-t_3)}\mathcal{G}_{ee^{\prime}}(t_2+t_1)\mathcal{G}_{gf}(-2\omega_p-i\sigma_p)\right.\\ &&\left.\times\mathcal{F}(2\omega_p-\omega_{fe}-i(\gamma_{fe}-\sigma_p),\omega_{fe}+i\gamma_{fe},2\omega_p-\omega_{fe^{\prime}}-i(\gamma_{fg}-\gamma_{e^{\prime}g}-\sigma_p)\right.\\ &&\left.-ie^{[2i\omega_p-\sigma_p]t_1-[i\omega_{fe}-\gamma_{fe}]t_1}\mathcal{G}_{ef}(t_3)\mathcal{G}_{ee^{\prime}}(t_2+t_1)\mathcal{G}_{gf}(-2\omega_p+i\sigma_p)\right.\\ &&\left.\times\mathcal{F}(2\omega_p-\omega_{fe}-i(\gamma_{fe}-\sigma_p),\omega_{fe}+i\gamma_{fe},2\omega_p-\omega_{fe^{\prime}}-i(\gamma_{fg}-\gamma_{e^{\prime}g}-\sigma_p)\right], \end{matrix}\end{eqnarray} \tag{ B.3 }$$

$$\begin{eqnarray}&&S_{LAP}^{(4)}({{t}_{1}},{{t}_{2}},{{t}_{3}})={{S}^{(1)*}}({{t}_{1}},{{t}_{2}},{{t}_{3}}),\end{eqnarray} \tag{ B.4 }$$

$$\begin{eqnarray}&&S_{LAP}^{(5)}({{t}_{1}},{{t}_{2}},{{t}_{3}})={{S}^{(2)*}}({{t}_{1}},{{t}_{2}},{{t}_{3}}),\end{eqnarray} \tag{ B.5 }$$

$$\begin{eqnarray}&&S_{LAP}^{(6)}({{t}_{1}},{{t}_{2}},{{t}_{3}})={{S}^{(3)*}}({{t}_{1}},{{t}_{2}},{{t}_{3}}).\end{eqnarray} \tag{ B.6 }$$