Chinese Journal of Chemical Physics  2018, Vol. 31 Issue (2): 135-151

#### The article information

Yu-xiang Weng

Detection of Electronic Coherence via Two-Dimensional Electronic Spectroscopy in Condensed Phase

Chinese Journal of Chemical Physics, 2018, 31(2): 135-151

http://dx.doi.org/10.1063/1674-0068/31/cjcp1803055

### Article history

Accepted on: April 21, 2018
Detection of Electronic Coherence via Two-Dimensional Electronic Spectroscopy in Condensed Phase
Yu-xiang Weng
Dated: Received on March 30, 2018; Accepted on April 21, 2018
Laboratory of Soft Matter Physics, Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China
*Author to whom correspondence should be addressed. Weng Yu-xiang, E-mail:yxweng@iphy.ac.cn
Abstract: Two dimensional Fourier transform electronic spectroscopy (2DES) in the visible region enables direct observation of complex dynamics of molecules including quantum coherence in the condensed phase.This review aims to provide a bridge between the principles and intuitive physical description of 2DES for tutorial purpose.Special emphasis is laid upon how 2DES circumvents the restrictions from both uncertainty principle and the wave-packet collapse during the coherent detection, leading to the successful detection of the coherence in terms of energy difference between the eigenstates showing as the quantum beats; then upon the possible mixing among the pure electronic transition, single-mode and multi-mode coupled vibronic transition leading to the observed beating phenomena.Finally, recent advances in experimentally distinguishing between the electronic coherence and the vibrational coherence are briefly discussed.
Key words: Quantum coherence    Uncertainty principle    Wave-packet collapse    Multi-mode vibronic coupling    Two dimensional electronic spectroscopy
Ⅰ. INTRODUCTION

Electronic coherence arising from the superposition of the excited-state wave functions of the atoms or molecules has been gradually shaped from a theoretical concept of quantum to the experimentally detectable reality, thanks to the advent of the ultrafast pulsed lasers. Among the various time-resolved spectroscopic methods [1-4], two-dimensional electronic nonlinear spectroscopy is one of the powerful tools in detection of electronic coherence, which has been widely used to study the coherent electronic energy transfer in a variety of photosynthetic light-harvesting systems [5-10], organic coherent intrachain energy migration in a conjugated polymer [11], exciton valley coherence in monolayer WSe$_2$ [12]. Interestingly, many quantum coherent effects were observed at room temperature in condensed phase, for the time scales involved in the relaxation and dephasing processes of excited electronic states are in the femtosecond regime in the condensed phase systems [13]. In solution, a spectroscopically labeled molecule constantly interacts with its surroundings, leading to fluctuations of its eigenstate energies on a time scale that is comparable to or longer than the "homogeneous" dephasing time, and the spectroscopically "labeled" subpopulation would be randomized after excitation due to the system-bath coupling. The time scale of frequency randomization varies, but for typical condensed phase systems it is generally below a few picoseconds, and often it is subpicosecond [14]. Thus, in principle the electronic coherence should be detected with a high resolution in both frequency and time domains. However, due to the problem arising from the Fourier transformation limit imposed by the laser pulse, i.e., a short pulse has a broad spectrum [15], both pump and probe frequency resolved ultrafast relaxation process can hardly be realized by means of conventional spectrometer, hence multi (two)-dimensional electronic spectroscopy has been developed. The basic idea of multi-dimensional spectroscopy is to perform multi-dimensional optical Fourier transform spectroscopy to reduce the spectral congestion that results from strong system-bath interactions by spreading frequency information on two or more axes. Furthermore, owing to the coherent coupling of the excited and ground state of the molecules, the spectral information regarding excitation and the detection is always entangled. The obvious advantage of 2DES over traditional one dimensional spectroscopy is that it can disentangle congested spectrum by spreading the spectrum out over two dimensions, i.e., excitation and probe frequencies, and the coupling between resonances can be identified by the presence of cross-peaks in the 2D spectrum, which provides an ideal tool for unraveling coupling among resonances and dynamics in complex systems [16]. 2D spectroscopy has been developed experimentally by Jonas's group in the near IR region [17, 18]. In the visible region, 2D spectroscopy was implemented by the groups of Fleming [19], Miller [20], Scholes [6], Ogilvie [21], et al., and 2DES has become a powerful tool suitable for study of electronic couplings in multichromophore molecular aggregates [7], coupled quantum dots [22], and other systems. 2D spectroscopy can trace the spectral evolution caused by energy fluctuation, energy transfer and delocalization with a high temporal resolution by detecting the correlation between the excitation and probe light, it becomes a powerful tool for detection of the electronic, vibronic and vibrational coherences within the molecules generated by coherent excitation. By recording 2D spectra as functions of the "waiting" or "population" time, the pathways by which a complex dynamical system evolves may be observed directly. Finally, because the 2D method is a form of photon echo spectroscopy, the inhomogeneous broadening is removed in the anti-diagonal direction [23]. In this review, we confine our scope to the 2D electronic spectroscopy based on the three-photon echo technique.

Ⅱ. PHYSICAL LIMITATIONS IMPOSED ON COHERENCE DETECTION

In measuring quantum coherence with ultrafast laser pulses, the information expected from the measurement at least should include (ⅰ) coherence time which can be as short as tens of femtosecond, and (ⅱ) the multiple energy levels comprising the coherence states. However, it seems that, in principles, there are two physical barriers that have to be overcome. One is the uncertainty principle; the other is the collapse of the wave-packet when subjected to any external perturbation like measurement. For the uncertainty principle, roughly speaking, it states that one cannot simultaneously determine the exact values of a pair of conjugate observables such as momentum and position of a physical system. Conjugate variables are pairs of variables mathematically defined in such a way that they become Fourier transform duals. The duality relations lead naturally to an uncertainty relation between them, in physics called the Heisenberge uncertainty principle [24]. While in matrix mechanics, this corresponds to that pair of observables whose operators do not commute [25]. The fact that short pulses probe dynamics and long pulses probe energy levels is indicative of the uncertainty principle, that time and energy resolution are related to each other through the Fourier transform [26].

A. Uncertainty Principle for $\Delta E$ and $\Delta t$ [27]

The Heisenberg uncertainty principle is a relationship between certain types of physical variables like position and momentum. The most well-known expression takes the position and momentum to be the conjugate variables. Another uncertainty relation which is often referenced in discussion of quantum mechanics is the energy-time uncertainty principle, which is directly related to the coherence detection. It is tempting to interpret this energy-time relation as the statement that a system may fluctuate in energy by an arbitrarily large amount over a sufficiently short time scale. This explanation is often given as a description for particle-antiparticle production and annihilation, where a particle and its antiparticle appear spontaneously from the vacuum briefly via "borrowed" energy before colliding and returning to vacuum. However, this explanation is not very precise and the given inequality is not so well-defined in quantum mechanics despite the nice physical interpretation. The reason that is not well-defined is because there is no operator in quantum mechanics corresponding to the measurement of time, although the Hamiltonian is the operator corresponding to energy. Nevertheless, there are some ways to make sense of an energy-time uncertainty principle by considering how the measurement of an arbitrary operator changes in time.

Since time is not an operator, it is unclear how time enters quantum mechanics at all. The answer is that time is incorporated into the Schrödinger equation, where it describes the time rate of change of a wave function. Physically, the passage of time is recorded by noting that certain physical observables are changing over time: for instance, perhaps the position of a particle is changing, which one interprets as motion over time, or the momentum of a particle is changing, which one interprets as accelerating or decelerating over time.

To quantify this statement, consider the Ehrenfest theorem governing the dynamics of the expectation value of an operator in terms of the commutator with the Hamiltonian [27]:

 $\begin{eqnarray} \frac{\textrm{d}}{{\textrm{d}t}}\left\langle A \right\rangle = \frac{i}{\hbar }\left\langle {\left[{\hat H, \hat A} \right]} \right\rangle + \left\langle {\frac{{\partial \hat A}}{{\partial t}}} \right\rangle \end{eqnarray}$ (1)

Assuming that $\hat A$ does not depend explicitly on time (which is typically the case), the second term vanishes and the expectation value of the $\left[{\hat H, \hat A} \right]$ is generally nonzero. Therefore, it satisfies the generalized Heiseinberg uncertainty relation:

 $\begin{eqnarray} \sigma _H^2 \sigma _A^2 \ge \left| {\frac{1}{{2i}}\left\langle {\left[{\hat H, \hat A} \right]} \right\rangle } \right|^2 \hspace{-0.15cm}& =&\hspace{-0.15cm}\left| {\frac{1}{{2i}}\frac{\hbar }{i}\frac{{\textrm{d}\left\langle {\hat A} \right\rangle }}{{\textrm{d}t}}} \right|^2 \nonumber\\ &=&\hspace{-0.15cm} \frac{{\hbar ^2 }}{4}\left| {\frac{{\textrm{d}\left\langle {\hat A} \right\rangle }}{{\textrm{d}t}}} \right|^2 \end{eqnarray}$ (2)

Notably, since the Hamiltonian is the energy operator, $\sigma _H$ corresponds to the uncertainty in energy.

Taking square roots now gives the relation:

 $\begin{eqnarray} \sigma _H \sigma _A \ge \frac{\hbar }{2}\left| {\frac{{\textrm{d}\left\langle {\hat A} \right\rangle }}{{\textrm{d}t}}} \right| \end{eqnarray}$ (3)

Define $\sigma _t$ by the relation:

 $\begin{eqnarray} \begin{array}{l} \displaystyle{\sigma _A = \left| {\frac{{\textrm{d}\left\langle {\hat A} \right\rangle }}{{\textrm{d}t}}} \right|\sigma _t } \\ \displaystyle{\frac{{\sigma _A }}{{|\textrm{d}\left\langle {\hat A} \right\rangle /\textrm{d}t|}} \approx \frac{{\sigma _A }}{{|\sigma _A /\sigma _t |}} = \sigma _t} \\ \end{array} \end{eqnarray}$ (4)

In other words, $\sigma _t$ corresponds to the time which takes the measured value of $\left\langle {\hat A} \right\rangle$ or shift by $\sigma _A$, for any operator $\hat A$. Given this definition, substituting Eq. (4) into Eq. (3) yields the relation:

 $\begin{eqnarray} \sigma _H \sigma _t \ge \frac{\hbar }{2} \end{eqnarray}$ (5)
B. Wave-packet collapse

Ultrashort pulse excitation creates a linear superposition of the eigenstates within the spectrum of the laser rather than a single eigenstate of the system excited by continuous working laser. The linear superposition is such as to create wave packet. The superposition principle states that a system is in all possible states at the same time, until it is measured [27]. In fact, such states are very fragile in the presence of dissipation, and would rapidly collapse to eigenstates after measurement, thus destroying the original configuration leaving no unusual interference features. This is also known as the collapse of the wave packet: when the system is excited with a laser pulse expressed as a pulse electric field $\vec E(t)$=$\vec E_0 \exp (-t^2 /\alpha \tau ^2) \cos (\omega t)$, in the optical frequency, $\alpha$=$2\pi ^2 /\ln 2$, $\tau$ is the pulse width. The resulting wave function of the system after excitation can be written as a superposition of the eigenstates of the system:

 $\begin{eqnarray} \psi = \sum\limits_{n = 1} {a_n \varphi _n } \end{eqnarray}$ (6)

where $a_n$ are the coefficients (or amplitudes) of the eigenstates $\varphi _n$. These coefficients, whose squared moduli represent the populations or occupation probabilities of the individual eigenstates after the pulse is applied, can be determined from the following first-order perturbation-theory formula:

 $\begin{eqnarray} a_n = A\left\langle {\varphi _n } \right|\left. {\chi _0 } \right\rangle \int\limits_{ - \infty }\hspace{-0.15cm}&&\hspace{-0.15cm} \exp [-i(\omega _n-\omega _0 )t] \times\nonumber\\&& \exp ( - t^2 /\alpha \tau ^2 )\cos (\omega t)\textrm{d}t \end{eqnarray}$ (7)
 $\begin{eqnarray}a_n \hspace{-0.15cm}&=&\hspace{-0.15cm} A\left\langle {\varphi _n } \right|\left. {\chi _0 } \right\rangle \exp [-(\omega _n-\omega _0 )^2 \alpha \tau ^2 /4] \end{eqnarray}$ (8)

$A$=$\vec \mu \cdot \vec E$ is the transition dipole moment, $\left\langle {\varphi _n } \right|\left. {\chi _0 } \right\rangle$ is the Franck-Condon factor, $h\omega _n$ is the eigen energy value of the individual excited state $\varphi _n$, and $h\omega _0$ is that of the ground state. Under the case of very short pulse excitation where $\tau$$\rightarrow0, one has  \begin{eqnarray} \psi = A\sum\limits_{n = 1}^{} {\left\langle {\varphi _n } \right|\left. {\chi _0 } \right\rangle \varphi _n } = \sum\limits_{n = 1} {c_n \varphi _n } \end{eqnarray} (9) where  \begin{eqnarray} c_n = A\left\langle {\varphi _n } \right|\left. {\chi _0 } \right\rangle \end{eqnarray} (10) When considering the time-evolution, one has  \begin{eqnarray} \psi (t) = \sum\limits_{n = 1}^{} {c_n \exp \left( { - iE_n t/\hbar } \right)\varphi _n } \end{eqnarray} (11) This wave packet is a coherent superposition of excited eigenstates. Over time, these states dephase with respect to one another. Ⅲ. HOMOGENEOUS AND INHOMOGENEOUS SPECTRAL BROADENING CAUSED BY ELECTRONIC DEPHASING The following describes the concepts of line broadening and optical dephasing pertaining to an ensemble comprising of two-level systems [28]. Considering an ensemble of identical oscillators that do not interact with the surrounding environment or each other, the system is free to evolve over time and will be correlated for long times, limited only by the natural lifetime \tau_e which is approximately 10^{-8} s (typically). In the time domain (for time\leq$$\tau_e$) the system will evolve as an undamped sinusoid which would approach an infinitesimally narrow line shape in the frequency domain after Fourier transformation (FIG. 1, the left panel). If the same ensemble is coupled to a weakly interacting local environment, the local environment will cause the system to depart from a perfect sinusoidal behavior as each two-level system interacts differently with the surrounding environment. This causes each oscillator in the ensemble to become out of phase, a process called optical dephasing. The system now begins to lose the long time correlation and decays to zero as $t $$\to$$\tau _e$. Thus, the process of optical dephasing of the oscillators in the ensemble causes damping in the time domain and broadens the absorption and fluorescence spectra. The two major types of dephasing or broadening mechanisms for molecules in solution are homogeneous and inhomogeneous broadening (FIG. 1). Two typical cases of oscillators interacting with a liquid solvent will now be considered. In the first scenario the surroundings of each oscillator change rapidly to produce a similar environment for each oscillator. The time evolution of each oscillator appears more stochastic in this case than in the natural lifetime limit due to the interaction between the oscillators and their surroundings changing rapidly with time. The optical dephasing in this case is said to be homogeneous and typically produces an exponential decay in the time-dependent correlation function $C(t)$ which results in a Lorentzian form of the broadened absorption and fluorescence line shape (FIG. 1, the middle panel). Dynamic interactions between a two-level system and its environment can be inelastic, in which an exchange of energy occurs, or elastic, where no energy is exchanged but the frequency is modulated in some way. The optical dephasing from inelastic collisions which contributes to the decay of the time-dependent correlation function is referred to as energy relaxation and is typically characterized by an exponential decay rate of $\displaystyle{\frac{1}{{2T_1 }}}$, where $T_1$ is the population relaxation time. The term pure dephasing is used to describe the effect of elastic interactions on the decay of the correlation function, or equivalently the broadening of the spectrum. The decay of the time-dependent correlation function due to pure dephasing is often characterized by an exponential decay rate of 1/$T_2^*$, where $T_2^*$ is the pure dephasing time. The decay rates of energy relaxation and pure dephasing are combined to produce a total homogeneous dephasing rate,

 FIG. 1 Schematic diagram illustrating the homogeneous and inhomogeneous spectral broadening in time and frequency domain for three hypothetical two-level systems. Left panel: a system of non-interacting harmonic oscillators; middle panel: a homogenously broadened system; right panel: an inhomogenously broadened system. $\nu_i$($t$): sample coordinate trajectories; $M$($t$): time-dependent correlation function; $I_a$($\omega$): absorption spectra [29].
 $\begin{eqnarray} \frac{1}{{T_2 }} = \frac{1}{{2T_1 }} + \frac{1}{{T_2^* }} \end{eqnarray}$ (12)

In the second hypothetical scenario each oscillator has a slightly different local environment, which does not change significantly over the course of time. In this case, the interaction between the inhomogeneous environments and the oscillators results in a shift in the energy separation away from the unperturbed energy $\hbar \omega _{{eg}}^0$, each frequency shift is different given that the ensemble of oscillators experiences different static environments. As the oscillators evolve over time they still exhibit sinusoidal behavior but oscillate at different frequencies due to their inhomogeneous surroundings. Assuming the distribution of frequencies is Gaussian, the time-dependent correlation function will have a Gaussian decay which results in a Gaussian broadened absorption and fluorescence lineshape when transformed (FIG. 1, the right panel). This hypothetical scenario is a limiting case of an inhomogeneously broadened system.

An absorption or fluorescence spectrum of a sample in a real solution may have both homogeneous and inhomogeneous characteristics. One could imagine a hypothetical two-level system between the two limiting cases in which a distribution of inhomogeneous environments arises from extremely slow solvent motions while dynamical interactions from faster solvent motions occur within each environment. Unfortunately, the linear response of system does not distinguish between dephasing processes of different time scales. The advantage of some nonlinear techniques is their ability to selectively eliminate inhomogeneous contributions and extract only the part for which dynamical (homogeneous) interactions dominate.

Ⅳ. PRINCIPLE OF 3PPE

The molecular system for nonlinear optical measurement is not a molecule of eigenstates, but a statistical ensemble of the inhomogeneous system, the temporal evolution of the ensemble can be described by the density function theory. Due to the ensemble inhomogeneity, any optical information retrieval occurs in the form of a photon echo [30], where photon echo is one of the first examples of an optical analogue of NMR.

For an ensemble of two-level systems characterized by a ground state $\left| g \right\rangle$ and an excited state $\left| e \right\rangle$, time independent eigen functions of the Hamiltonian with eigen energies being $E_g$ and $E_e$ superposition states, the wave function superposition state is

 $\begin{eqnarray}\left| {\psi (t)} \right\rangle = c_g e^{ - i\omega _g t} \left| g \right\rangle + c_e e^{ - i\omega _e t} \left| e \right\rangle \end{eqnarray}$ (13)

It should be noted that the observables corresponding to operators that do not commute with the Hamiltonian, such as the dipole operator, will oscillate over time with frequency $\omega$=$\omega _e $$-$$\omega _g$. This phase factor is a unique feature of quantum mechanics. We anticipate here that the quantum beating in such observables is a direct consequence of the superposition character [31]. To generalize this treatment for ensembles containing a very large number of particles, the common approach is to introduce the density operator defined as $\rho$=$\left| \psi \right\rangle \left\langle \psi \right|$, so that:

 $\begin{eqnarray} \rho &=& \sum\limits_{ij} {\left\langle {c_i (t)c_j^* (t)} \right\rangle } \left| {\varphi _i } \right\rangle \left\langle {\varphi _j } \right| \nonumber\\ &= &\sum\limits_{ij} {\rho _{ij} \left| {\varphi _i } \right\rangle \left\langle {\varphi _j } \right|} \end{eqnarray}$ (14)

In the matrix representation of the operator, each diagonal element $\rho _{{ii}}$ ("population") represents the average probability of finding a molecule in state $i$ at time $t$, whereas the off-diagonal elements $\rho _{{ij}}$ ("coherences") reflect the phase coherence of the ensemble and are complex numbers carrying phase information. Obviously, the density matrix reflects the temporal evolution of the wave functions which is directly related to the time-dependent Schrödinger equation. The time-dependent Schrödinger equation $\hat H\Psi$=$i\hbar \displaystyle{\frac{{\textrm{d}\Psi }}{{\textrm{d}t}}}$ with a time-dependent wave function $\Psi (r, t)$=$\displaystyle{\sum {c_n (t)u_n (r)}}$ consisting of a temporal part $c_n (t)$ and a spatial part $u_n (r)$ is equivalent to the density matrix expression called Liouville-von Neumann equation:

 $\begin{eqnarray} \frac{{\partial \rho _{nm} }}{{\partial t}} \hspace{-0.15cm}&=& \hspace{-0.15cm}\frac{i}{\hbar }\left({\sum\limits_k {\rho _{nk} H_{km} - H_{nk} \rho _{km} } } \right) \nonumber\\ &=&\hspace{-0.15cm} \frac{i}{\hbar }\left[{\rho _n H_m-H_n \rho _m } \right] \nonumber\\ &=&\hspace{-0.15cm}\frac{i}{\hbar }\left[{\hat \rho, \hat H} \right] \end{eqnarray}$ (15)

where

 $\begin{eqnarray}\rho _{nm} = c_n c_m^* \end{eqnarray}$ (16)
A. Evolution of density matrix interpreted as photo echoes [32-34]

Considering an optical transition in a two-level system, e.g., the ground state and an excited state of a dye molecule having two different energies $E_g$ and $E_e$, are involved in the interaction. The external force is $F$=$e$$\cdot$${\bf{E(r, t)}}$, where $e$ is the charge of the electron, $\textbf{r}$ is the position vector in the three dimensional space, and $\textbf{E}$ is the electric field. The electric field is: ${\bf{E(r, t)}}$=${\bf{E}}_{\bf{0}} \cos (\omega t)$.

Using $F$=$-\nabla V$ and the electric dipole moment is $\mu$=$-e \cdot {\bf{r}}$, the external potential can be expressed as:

 $\begin{eqnarray} V \hspace{-0.15cm}&= &\hspace{-0.15cm} - \mu \cdot {\bf{E(r, t)}} \end{eqnarray}$ (17)
 $\begin{eqnarray} \mu _{gg} \hspace{-0.15cm} &=&\hspace{-0.15cm} \mu _{ee} = 0 \end{eqnarray}$ (18)
 $\begin{eqnarray} \mu _{ge} \hspace{-0.15cm}&=&\hspace{-0.15cm}\ \mu _{eg} = \mu \end{eqnarray}$ (19)

When molecules interact with electrical field ${\bf{E(r, t)}}$, the ensemble average of $\left\langle \mu \right\rangle$ would be

 $\begin{eqnarray} \left\langle \mu \right\rangle &=& \textrm{tr}(\rho \mu ) \nonumber\\ &=& \rho _{ge} \mu _{eg} + \rho _{eg} \mu _{ge} + \rho _{gg} \mu _{gg} + \rho _{ee} \mu _{ee} \nonumber\\ &=& \mu (\rho _{ge} + \rho _{eg} ) \end{eqnarray}$ (20)

Then consider the radiation field whose frequency fulfills $\omega _{eg}$=$({E_e}-{E_g})/\hbar$. As a result, the density matrix is reduced to a 2$\times$2 matrix with the elements $\rho _{gg}$, $\rho _{ge}$, $\rho _{eg}$ and $\rho _{ee}$.

Introducing Rabi frequency

 $\begin{eqnarray} \Omega \hspace{-0.15cm}&=& \hspace{-0.15cm}\frac{V}{\hbar }\nonumber\\ \hspace{-0.15cm}&=& \hspace{-0.15cm}\frac{{ - {\bf{ \pmb{\mathit{ μ}} }} \cdot {\bf{E}}(r, t)}}{\hbar } \end{eqnarray}$ (21)
 $\begin{eqnarray} H(t) \hspace{-0.15cm}&=& \hspace{-0.15cm}H_0 + V(t) \end{eqnarray}$ (22)
 $\begin{eqnarray} \frac{{\textrm{d}\rho _{gg} }}{{\textrm{d}t}}\hspace{-0.15cm}& = &\hspace{-0.15cm} - \frac{i}{\hbar }(V_{ge} \rho _{eg} - V_{eg} \rho _{ge} ) \end{eqnarray}$ (23)
 $\begin{eqnarray} \frac{{\textrm{d}\rho _{ee} }}{{\textrm{d}t}}\hspace{-0.15cm}& =&\hspace{-0.15cm} - \frac{i}{\hbar }(V_{eg} \rho _{ge} - V_{ge} \rho _{eg} ) \end{eqnarray}$ (24)
 $\begin{eqnarray} \frac{{\textrm{d}\rho _{ge} }}{{\textrm{d}t}} \hspace{-0.15cm}&=&\hspace{-0.15cm} - \frac{i}{\hbar }(E_g - E_e )\rho _{ge} - \frac{i}{\hbar }V_{ge} (\rho _{ee} - \rho _{gg} ) \nonumber\\ &=& i\omega _{eg} \rho _{ge} - i\Omega (\rho _{ee} - \rho _{gg} ) \end{eqnarray}$ (25)
 $\begin{eqnarray} \frac{{\textrm{d}\rho _{eg} }}{{\textrm{d}t}} &=& - \frac{i}{\hbar }(E_e - E_g )\rho _{eg} - \frac{i}{\hbar }V_{eg} (\rho _{gg} - \rho _{ee} ) \nonumber\\ &=& - i\omega _{eg} \rho _e g + i\Omega (\rho _{ee} - \rho _{gg} ) \end{eqnarray}$ (26)

Apparently the solutions of $\rho _{ge} (t)$ and $\rho _{eg} (t)$ would contain a time varying factor of the type: $\exp (\pm i\omega _{eg} t)$. What we want is a time-independent, or at least slowly varying solution when the frequency of the perturbing field is close or equal to the transition frequency, i.e., $\omega$$\approx$$\omega _{eg}$. To realize this, the rotating wave approximation (RWA) is employed by multiplying $\rho _{ge} (t)$ and $\rho _{eg} (t)$ with the complex conjugate of their time depending factor:

 $\begin{eqnarray} \hat \rho _{gg} &=& \rho _{gg} \end{eqnarray}$ (27)
 $\begin{eqnarray} \hat \rho _{ee} &= &\rho _{ee} \end{eqnarray}$ (28)
 $\begin{eqnarray} \tilde \rho _{ge}&\equiv&\textrm{e}^{ - i\omega t} \cdot \rho _{ge} (t) \end{eqnarray}$ (29)
 ${\tilde \rho _{eg}}{\rm{ }} \equiv {\rm{ }}{{\rm{e}}^{i\omega t}} \cdot {\rho _{{e_g}}}(t)$ (30)

Then we have

 $\begin{eqnarray} \frac{{\textrm{d}\tilde \rho _{ge} }}{{\textrm{d}t}} = - i(\omega - \omega _{eg} ) \cdot \tilde \rho _{ge} - i\textrm{e}^{ - i\omega t} \cdot \Omega \cdot (\rho _{ee} - \rho _{gg} ) \end{eqnarray}$ (31)

Rewrite the cosine form of the electrical field in the Rabi frequency into a complex form:

 $\begin{eqnarray} \Omega (t) \hspace{-0.15cm}&=&\hspace{-0.15cm} \frac{1}{2} \cdot \left| \Omega \right| \cdot (\textrm{e}^{i\omega t} +\textrm{ e}^{ - i\omega t} ) \end{eqnarray}$ (32)
 $\begin{eqnarray} \tilde \Omega (t) \hspace{-0.15cm}&=&\hspace{-0.15cm} \Omega (t) \cdot \textrm{e}^{ - i\omega t} = \frac{1}{2} \cdot \left| \Omega \right| \cdot (1 + \textrm{e}^{ - 2i\omega t} ) \end{eqnarray}$ (33)

$\tilde \Omega (t)$ contains a time-independent term and a term oscillating on twice the field frequency, 2$\omega$. The effect of the later term is expected to be small since it oscillates much faster than any other frequency in the system. We thus neglect the time-dependent part, denoting as $\tilde \Omega (t)$=$\tilde \Omega ^* (t) \cong \displaystyle{\frac{1}{2}}\Omega$. This approximation is called the rotating wave approximation. By introducing the frequency detuning as $\Delta$=$\omega -\omega _{eg}$, we have

 $\begin{eqnarray} \frac{{\textrm{d}\tilde \rho _{ge} }}{{\textrm{d}t}} = - i\Delta \cdot \tilde \rho _{ge} - \frac{i}{2}\Omega \cdot (\rho _{ee} - \rho _{gg} ) \end{eqnarray}$ (34)
 $\begin{array}{l} \frac{{\rm{d}}}{{{\rm{d}}t}}\left[{\begin{array}{*{20}{c}} {{{\tilde \rho }_{ge}}}\\ {{{\tilde \rho }_{eg}}}\\ {{{\tilde \rho }_{gg}}}\\ {{{\tilde \rho }_{ee}}} \end{array}} \right] = - i\left[{\begin{array}{*{20}{c}} \Delta &0&{-\tilde \Omega (t)}&{\tilde \Omega (t)}\\ 0&{-\Delta }&{{{\tilde \Omega }^*}(t)}&{-{{\tilde \Omega }^*}(t)}\\ { - {{\tilde \Omega }^*}(t)}&{\tilde \Omega (t)}&0&0\\ {{{\tilde \Omega }^*}(t)}&{ - \tilde \Omega (t)}&0&0 \end{array}} \right]\left[{\begin{array}{*{20}{c}} {{{\tilde \rho }_{ge}}}\\ {{{\tilde \rho }_{eg}}}\\ {{{\tilde \rho }_{gg}}}\\ {{{\tilde \rho }_{ee}}} \end{array}} \right]\\ = - i\left[{\begin{array}{*{20}{c}} \Delta &0&{-\frac{1}{2}\Omega }&{\frac{1}{2}\Omega }\\ 0&{-\Delta }&{\frac{1}{2}\Omega }&{-\frac{1}{2}\Omega }\\ { - \frac{1}{2}\Omega }&{\frac{1}{2}\Omega }&0&0\\ {\frac{1}{2}\Omega }&{ - \frac{1}{2}\Omega }&0&0 \end{array}} \right]\left[{\begin{array}{*{20}{c}} {{{\tilde \rho }_{ge}}}\\ {{{\tilde \rho }_{eg}}}\\ {{{\tilde \rho }_{gg}}}\\ {{{\tilde \rho }_{ee}}} \end{array}} \right] \end{array}$ (35)

The resulting Liouville-von Neumann equation can now be recast in the form

 $\frac{{{\rm{d}}\tilde \rho }}{{{\rm{d}}t}}{\rm{ }} =- \frac{i}{\hbar }\left[{{{\hat H}_{eff, }}\hat \rho } \right]$ (36)
 ${\hat H_{eff}} = \left[{\begin{array}{*{20}{c}} {\hbar \Delta }&{\frac{\hbar }{2}\Omega }\\ {\frac{\hbar }{2}\Omega }&0 \end{array}} \right]$ (37)

In this model, no relaxation is introduced so far. Now we introduce relaxation by two phenomenological constants. The population of the excited state decays to the ground state at a decay rate written as 1/$T_1$ and the decay rate of the coherence as $1/T_2^*$, this relaxation happens due to phonons in the molecule. $T_2^*$ is called the irreversible loss off coherence time. $T_2$ is a combination of the two relaxation times above according to Eq.(12), which is called as the coherence decay rate. Then we have

 $\begin{eqnarray} \begin{array}{l} \displaystyle{\frac{{\textrm{d}\tilde \rho _{ee} }}{{\textrm{d}t}} = - \frac{i}{2}\Omega \cdot (\tilde \rho _{eg} - \tilde \rho _{ge} ) - \tilde \rho _{ee} /T_1 } \\[2ex] \displaystyle{\frac{{\textrm{d}\tilde \rho _{gg} }}{{\textrm{d}t}} = - \frac{i}{2}\Omega \cdot (\tilde \rho _{ge} - \tilde \rho _{eg} ) + \tilde \rho _{ee} /T_1} \\[2ex] \displaystyle{\frac{{\textrm{d}\tilde \rho _{eg} }}{{\textrm{d}t}} = i\Delta \cdot \tilde \rho _{eg} + \frac{i}{2}\Omega \cdot (\tilde \rho _{ee} - \tilde \rho _{gg} ) - \tilde \rho _{eg} /T_2 } \\[2ex] \displaystyle{\frac{{\textrm{d}\tilde \rho _{ge} }}{{\textrm{d}t}} = - i\Delta \cdot \tilde \rho _{ge} - \frac{i}{2}\Omega \cdot (\tilde \rho _{ee} - \tilde \rho _{gg} ) - \tilde \rho _{ge} /T_2 } \\ \end{array} \end{eqnarray}$ (38)

The term $\tilde \rho _{gg}$, $\tilde \rho _{ee}$, $\tilde \rho _{ge}$, and $\tilde \rho _{eg}$ of the density matrix describe the molecular or atomic state and its evolution when a electromagnetic field is applied.

B. Bloch vector

The Bloch vector formalism consists in replacing the density matrix elements by three real components $(u, v, w)$ with a population conservation equation $\rho _{gg}$+$\rho _{ee}$=1. The components are defined as

 $\begin{eqnarray} \begin{array}{l} \displaystyle{ \tilde \rho _{gg} = \frac{{1 - w}}{2}} \\[2ex] \displaystyle{ \tilde \rho _{ee} = \frac{{1 + w}}{2}} \\[2ex] \displaystyle{\tilde \rho _{eg} = \frac{{u - iv}}{2}} \\[2ex] \displaystyle{ \tilde \rho _{ge} = \frac{{u + iv}}{2}} \\[2ex] \end{array} \end{eqnarray}$ (39)

or

 $\begin{eqnarray}\begin{array}{l} w = \tilde \rho _{ee} - \tilde \rho _{gg} \\ u = \tilde \rho _{eg} + \tilde \rho _{ge} \\ v = i(\tilde \rho _{eg} - \tilde \rho _{ge} ) \\ \end{array} \end{eqnarray}$ (40)

Notice that $u$, $v$, and $w$ are real. $w$ represents the population difference and $u$ and $v$ are the real and imaginary parts of the off-diagonal elements in the density matrix, and thus represents the real and imaginary parts of the optical coherence. Also $u$ and $v$ are sometimes called the in-phase and put-of-phase components of Bloch vector.

The Bloch vector equations of motion thus can be written as

 $\begin{eqnarray} \begin{array}{l} \dot u = - u/T_2 + \Delta v + \Omega w \\ \dot v = - \Delta u - v/T_2 - \Omega w \\ \dot w = \Omega (v - u) - (1 + w)/T_1 \\ \end{array} \end{eqnarray}$ (41)
 $\begin{eqnarray} \left[{\begin{array}{*{20}c} {\dot u} \\ {\dot v} \\ {\dot w} \\ \end{array}} \right] = \left[{\begin{array}{*{20}c} 0&\Delta &\Omega \\ {-\Delta }&0&{-\Omega } \\ {-\Omega }&\Omega &0 \\ \end{array}} \right]\left[{\begin{array}{*{20}c} u \\ v \\ w \\ \end{array}} \right] - \left[{\begin{array}{*{20}c} {u/T_2 } \\ {v/T_2 } \\ {(1 + w)/T_1 } \\ \end{array}} \right] \end{eqnarray}$ (42)

If the relaxation process is disregarded, the optical Bloch equation has the same form as that on NMR by setting

 $\begin{eqnarray} \vec S(t) = \left[{\begin{array}{*{20}c} u \\ v \\ w \\ \end{array}} \right] \end{eqnarray}$ (43)
 $\begin{eqnarray} \vec \Omega (t) = \left[{\begin{array}{*{20}c} \Omega \\ \Omega \\ {-\Delta } \\ \end{array}} \right] \end{eqnarray}$ (44)

And the coherent evolution part can be represented by

 $\begin{eqnarray} \frac{{\textrm{d}\vec S(t)}}{{\textrm{d}t}} = \vec \Omega (t) \times \vec S(t) \end{eqnarray}$ (45)
C. Bloch diagram

Neglecting the relaxation processes, the length of the Bloch vector is thus $\left| {\bf{S}} \right| = \sqrt {u^2 + v^2 + w^2 }$.

The whole Bloch diagram is rotating around the $w$-axis at the field frequency $\omega$, where the RWA is introduced. When the whole population is in the ground level, the arrow points to $w$=-1. How do the relaxation processes alter the direction and length of the Bloch vector? The factor 1/$T_1$ relaxes the component in the $w$-direction towards the value $w$=-1. The factor 1/$T_2$ strives to let $u$$\to0 and v$$\to$0. The latter leads to a decrease of the length of the Bloch vector.

The detuning factor $\Delta$ will causes $\textbf{S}$ to start to smear out like a fan in the ($u$, $v$)-plane. This happens since all the different molecules in the ensemble have slightly different excitation energies, due to individual surroundings interaction with the matrix, etc. Because the Bloch diagram is rotating around the $w$-axis on the field frequency $\omega$, a molecule excited with a slightly different frequency, $\omega _m$, will have its own Bloch vector rotating on the detuning frequency $\Delta$=$\omega _m -\omega$ with respect to the coordinate system.

An electric field turns the Bloch vector around the $u$-axis by a certain number of degrees depending on its integrated intensity. A $\pi$/2-pulse causes the Bloch vector to make a $\pi$/2 degree turn around the $u$-axis, while a $\pi$-pulse turns the Bloch vector $\pi$ degree around the $u$-axis as shown in FIG. 2(c, d) [35].

 FIG. 2 Bloch diagrams in a rotation frame with $w$ as the rotating frequency. (a) The Bloch vector, $\textbf{S}$, pointing out with a length of unit in the space ($u$, $v$, and $w$). (b) The Bloch vector starts to smear out in the ($u$, $v$)-plane, because of the slightly different excitation energies of the molecules in the ensemble. (c) and (d) How the Bloch vector is turned when a $\pi$/2-pulse, respectively, a $\pi$-pulse is interacting with the ensemble of molecules, all in the ground level.

In photo echo experiments, the $\pi$-or $\pi$/2-pulse excitation is realized by the sequenced laser pulses separated by a time delay $\tau$. FIG. 3 shows how Bloch vectors change in a two-photon echo process. The initial state is at $w$=-1 with ${\bf{S}}$ vector pointing at the south pole, and the first $\pi$/2 pulse excitation coherently turn the ${\bf{S}}$ vector around $u$ axis by $\pi$/2 to overlap with $v$ axis, then the decoherence occurs by smearing out ${\bf{S}}$ for each member of the inhomogeneous ensemble from their initial overlapped sate at a rate proportional to the frequency detuning $\Delta$=$\omega _m-\omega$ within the $uv$ plane owing to the frequency detuning during the interval $\tau$. Then the second $\pi$-pulse comes, which flips the smeared vectors ${\bf{S}}$ around $u$ axis by $\pi$ to overlap with $-v$ axis while keeping the precessing direction of each individual ${\bf{S}}$ same as that before flipping. It should be noted that after the $\pi$-pulse excitation, the smearing motion becomes reversed as a converging process, which is exactly the reverse process of the decoherence. It is expected, after an evolution time same as the interval $\tau$, the coherence is resumed and a stimulated emission occurs as the photo echo. Photo echo is a physical vision of the time evolution of the density matrix elements, however in reality, the direction in which the echo is emitted, $\vec k_e$, is determined by the wave vectors of the two excitation pulses, $\vec k_1$ and $\vec k_2$, according to the following equation based on the momentum conservation in the four-wave mixing of nonlinear optics as shown in FIG. 4 [36]:

 FIG. 3 How the Bloch vectors change according to the excitation pulses separated in time by τ in 2PPE. The blue arrows indicate the direction of precessing for each individual vector.
 FIG. 4 Phase-matched optical configuration in 3PPE. (a) The echo is emitted in the above shown direction. The direction is given by the two excitation pulses wave vectors. (b) Various phase-matched stimulated echo beams on the projection plane [34].
 $\begin{eqnarray} \vec k_e = 2\vec k_2 - \vec k_1 \end{eqnarray}$ (46)

If the temporal order of the two writing beams is changed, the echo is emitted in the direction implied as a dashed arrow above, which provides a convenient way to acquire the rephasing and non-rephaisng signals in the 2DES.

Although in the time sequence, there appear two excitation pulses called two-pulse photo echo. However, in terms of nonlinear wave mixing, it consists of three beams, i.e., one beam 1 and two degenerated beam 2 having the same propagating direction and a zero respective delay time. Therefore, by nature, 2PPE is a simplified version of three-photo echo emission.

The three-pulse photon echo (3PPE) experiment involves three pulses with time sequence as $\tau$, $T$, and $t$ as shown in FIG. 3, where $T$=0. Like 2PPE, the first $\pi$/2 pulse again creates a macroscopic polarization, then at delay time $\tau$ is a $\pi$/2-pulse rather than a $\pi$-pulse. This second $\pi$/2-pulse creates population in either the ground or excited states in which the coherence is stored [37]. A third $\pi$/2-pulse then interacts with the sample at delay time $T$ and causes the system to rephase the inhomogeneous component of the polarization decay to produce a photon echo at time $t$. The echo signal is measured in the phase matching direction $\vec k_e$=$\vec k_3 + \vec k_2 -\vec k_1$ and is maximum for three $\pi$/2-pulses [28].

D. 3PPE in view of density matrix

FIG. 5 shows the stimulated photon echo pulse sequence along with the 2$\times$2 density matrix for a two-level system [14]. The sequence has three time intervals, $\tau$, $T$, and $t$, where due to interactions with the three fields $\textbf{E}_1$, $\textbf{E}_2$, and $\textbf{E}_3$, the system evolves alternately in coherences ($\rho _{jk}$) and populations ($\rho _{jj}$). An "echo" is the phase conjugation between the two coherence evolution periods. During $\tau$ the density is $\rho_{ge}$ and during $t$ it is $\rho_{eg}$, thus effectively reversing the sign of time in the field-free Green function responsible for the system's evolution during the second coherence interval. The lower set of density matrices shown during $\tau$ and $t$ indicate an alternative pathway where the system returns to the ground state following the second excitation pulse. Both the excited state and ground state paths are required to evaluate the system response to the applied fields. The fact that echo signal $\rho_{ge}$ is a phase-conjugated term underlies the principle for 3PPE as a quantum data storage technique in time sequence.

 FIG. 5 Illustration of the evolution for the density matrix elements of two-level system in 3PPE. During $\tau$, the system evolves in a coherence. The system is stored in a population on the ground or excited state during $T$, and the final interaction induces a second coherence, which radiates a stimulated photo echo signal in the phase match direction [14].

In summary, the conventional interpretation of a photon echo experiment represents the state of the system using a pseudo spin vector on the Bloch sphere. By transforming to the rotating frame, each member of the inhomogeneous ensemble precesses away from the origin at a rate proportional to the difference between the specific frequency $\omega_m$ and the frame rotation frequency $\omega$, $i.e.$, $\Delta$. When $T$=0, the second and third pulses flip the vectors leading to the phase conjugation. The echo is emitted when the vectors precess freely back to the origin, hence the gradual increase and decrease of the echo signal can be detected. Choosing $T$$>0 permits incoherent relaxation to occur before the refocusing pulse, providing access to lifetime measurements and observation of spectral diffusion processes due to microscopic site randomization. Indeed, it is precisely the T-dependence that provides much of the chemical information in condensed phase photon echo experiments. The homogenous and inhomogeneous decay process of the system in response to respective one-pulse, two-pulses and three-pulses excitation is depicted schematically in FIG. 6 [28].  FIG. 6 Schematic diagram of coherent techniques involving one, two and three pulses. The electric field pulses (blue solid peaks) excite the sample. The resulting coherent signal is shown for homogeneously broadened transitions (red dashed line) and for inhomogenously broadened transitions (green dotted and dashed line). (a) A single excitation pulse {\bf{E}}_1 with a wave vector \vec k_1. Homogeneously broadened transitions lead to an exponential decay with dephasing time of T_2^*, while interference of an inhomogenously distribution of the resonant frequencies causes a faster decay. (b) Two excitation pulses {\bf{E}}_1 and {\bf{E}}_2 with wave vectors of \vec k_1 and \vec k_2 gives an echo in the direction 2\vec k_2 -\vec k_1. For inhomogeneously broadened transitions, the echo is emitted at t=\tau. The echo amplitude measured as a function of \tau allows determination of the homogeneously dephasing time constant T_2^*. (c) When three excitation pulses {\bf{E}}_1, {\bf{E}}_2, and {\bf{E}}_3 are applied, the coherent signal emitted in the direction \vec k_3+\vec k_2 -\vec k_1 allows determination of the dephasing time T_2^* and the energy relaxation time T_1 [38]. E. 3PPE as time-domain data storage 3PPE has been shown as potential technique for time-domain data storage [39] and word by word logic processing [40] possibly used in quantum computation. The idea of photon-echo memory is to perform spectral hole-burning memory [37] in the time domain by utilizing the temporal interference effect of two optical pulses. From the viewpoint of memory, the three excitation pulses are renamed as WRITE, DATA and READ pulses and, rather than being a delta-function pulse, the DATA pulse carries some kind of temporal information as shown in FIG. 7. Calculation shows that this echo signal simply replicates the shape of the DATA pulse. In terms of memory, the WRITE pulse as a trigger and the DATA pulse with temporal information complete the storage of the data into the medium. The readout process starts by applying a trigger called a READ pulse, and the echo then appears as the OUTPUT DATA. The DATA can be understood in spectral region, since pulse 1 and pulse 2 would create spectral interference which provides spectral selectivity by tuning the relative delay \tau.  FIG. 7 Pulse sequence and expected signal in the stimulated photon echo (upper panel) and the simulated photon-echo memory (lower panel) [39]. Ⅴ. THIRD ORDER NONLINEAR POLARIZABILITY THEROY FOR 3PPE 3PPE is one of implementations of 2DES [41]. In experiment, three ultrashort pulses are arranged in a unique temporal sequence to excite the sample, creating a 3rd order polarization P^{(3)}(\tau, T, t) which can be expressed as the convolution of the nonlinear response function R^{(3)} (\tau, T, t) with the respective three field envelope {\bf{E}}_i (\vec k_i, t):  \begin{array}{l} {P^{(3)}}(\tau, T, t) \propto \int\limits_0^\infty {{\rm{d}}\tau \int\limits_0^\infty {{\rm{d}}T} } \int\limits_0^\infty {{\rm{d}}t} {R^{(3)}}(\tau, T, t){E_1}({{\vec k}_1}, \tau )\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{E_2}({{\vec k}_2}, T){E_3}({{\vec k}_2}, t) \end{array} (47) where the third order nonlinear response function is expressed as a function of the dipole moment and the density matrix of the unperturbed system \rho (-\infty):  \begin{gathered} {R^{(3)}}({t_3}, {t_2}, {t_1}) = {\left( {- \frac{i}{\hbar }} \right)^n}\left\langle {\mu ({t_3} + {t_2} + {t_1})[\mu ({t_2} + {t_1}), } \right. \hfill \\ \left. {\left[{\mu ({t_1}), [\mu ({t_0}), \rho (-\infty )]]]} \right.} \right\rangle \hfill \\ \end{gathered} (48) Apparently, P^{(3)} (\tau, T, t) contains three interacting electric fields, when the laser pulses E_1(t), E_2(t), and E_3(t) are shorter than the time separation between them (semi-impulsive limit), they do not overlap in time which gives rise to a strict time ordering. In that case, the first interaction \mu(t_0) originates from the pulse E_1(t), \mu_1(t) from E_2(t), and so on. This reduces the number of terms in P^{(3)} (\tau, T, t) to 32 corresponding to all possible contributions from three incoming electric fields to polarization. Each of these terms can be associated with a distinct time evolution of the density matrix, known as Liouville pathways [34]. In spectroscopy these pathways are usually graphically visualized by means of diagrams capable of highlighting the temporal sequence of the field interactions and the transitions promoted in the systems by such interactions, while the most famous are the double-sided Feynman diagrams. The fourth homologous laser pulse serves as a local oscillator to detect the emitted signal by means of optical heterodyne detection as shown in FIG. 8 [42]. The detected signal by the spectrometer S(\tau, T, \omega _t) is related to the 3rd polarization P^{(3)} (\tau, T, t) via one-dimensional Fourier transformation on the detection period t at the phase-matched condition \vec k_s=\vec k_2$$-\vec k_1$+$\vec k_3$.

 FIG. 8 Pulse sequence in 3PPE. LO: local oscillator for heterodyne detection.
 $\begin{eqnarray} \tilde S(\tau, T, \omega _t ) &=& i \cdot sig(\omega _t )\int_{ - \infty }^\infty {P^{(3)} (\tau, T, t)\exp (i\omega _t t)\textrm{d}t} \nonumber\\ &=&i \cdot sig(\omega _t )\tilde P^{(3)} (\tau, T, \omega _t ) \end{eqnarray}$ (49)

where the third-order polarization is related to the measured signal field by $\tilde P^{(3)} (\tau, T, \omega _t)$$\propto$$\tilde E(\tau, T, \omega _t)/\omega_t$. Here the spectrometer acts as a hardware of performing Fourier transformation along the detection time $t$ resulting in a probed spectrum spanned along $\omega_t$, leaving the pump frequency to be resolved from the time space $\tau$. Further performing the 2D-FT of the experimentally measured nonlinear optical response function along the evolution time $\tau$, one has the 2D-FT spectrum at a fixed $T$ [43-45]:

 $\begin{eqnarray} \int_{ - \infty }^\infty {P^{(3)} (\tau, T, \omega _t )\exp (i\omega _\tau t)\textrm{d}t} {\rm{ = }}\tilde P^{(3)} (\omega _\tau, T, \omega _t ) \end{eqnarray}$ (50)

The 2D spectrum at a particular $T$ involves frequency dimensions $\omega_\tau$ and $\omega_t$, corresponding to the evolution time $\tau$ (coherence time) and the detection time $t$ respectively. In 2D-FT spectroscopy, two different phase-matched signals are typically recorded: the rephasing and non-rephasing signals in the $k_2$-$k_1$+$k_3$ and $k_1$-$k_2$+$k_3$ directions, respectively. When summed, these signals yield the absorptive 2D-FT spectrum, free from broadening refractive contributions. Another signal, with the phase-matching direction $k_1$+$k_2$-$k_3$, yields information about double-quantum coherences [43]. Typical 2DES is shown in FIG. 9.

 FIG. 9 Schematic illustration of the information contained in absorptive two-dimensional Fourier transform spectra. The observed spectral range shown in the unshaded region is determined by the excitation and detection laser bandwidths. (a) The $T$=0 correlation spectrum reveals homogeneous and inhomogeneous line widths, excitonic coupling, and excited state absorption (ESA) features. (b) At $T$$>0, the broadening of peaks in the anti-diagonal direction reflects spectral diffusion. The growth of cross peaks indicates energy transfer. The emergence of an entirely new peak at later waiting times, such as the B'D cross peak in panel b, represents the formation of a new product species, populated upon excitation of B', that absorbs at D (cited from Ref.[41]). The time resolution of dynamics over the waiting time, in both pump-probe and 2D experiments, is dependent on the duration of the pump pulse, with shorter pulses giving higher time resolution. The commensurately larger bandwidth of a shorter pulse leads to a lower spectral resolution as restricted by the uncertainty principle. 2D-FT spectroscopy use the Fourier transform methodology by converting time to the frequency, thus the restriction imposed by the uncertainty principle in 2D spectroscopy is circumvented [43], and the result has both high temporal and spectral resolution, limited only by the signal-to-noise ratio. However, the wave-packet collapse to the eigenstates as reflected in both excitation and detection frequency axis at any fixed waiting time, i.e., coherence cannot be informed from the corresponding 2D spectra. Ⅵ. ELECTRONIC COHERENCE (QUANTUM BEAT) ALONG WAITING TIME AND DOUBLE-SIDED FEYNMAN DIAGRAM A. Double-sided Feynman diagram [31] Double-sided Feynman diagrams describe the evolution of the density matrix in the presence of several coherent fields in time sequence [46, 47] by depicting interaction sequences between external perturbations such as electromagnetic fields and ensemble of molecules. One advantage of using Feynman diagrams is to physically capture the most important interaction paths without referring to the actual form of either the density matrix or the Hamiltonians describing the ensemble. The drawing of double-sided Feynman diagrams derives from taking separately each term of the commutator expansion (Eq.(47) and Eq.(48)) to create the corresponding diagram, using specific rules: (ⅰ) two vertical lines represent the time evolution of the ket (left line) and bra (right line) of the density matrix, with time running from the bottom to the top. (ⅱ) The interactions between the density matrix and the electric field are indicated by arrows acting either on the bra or ket side pointing towards or away from the system corresponding to absorption or emission. (ⅲ) The ket and bra sides of the diagram are the complex conjugate of each other and the last emission is always from the ket side, by convention. (iv) The last interaction must end in a population state. (v) An arrow pointing to the right represents an electric field with rm{e}^{-i\omega t + ikr}(+k), an arrow pointing to the left represents an electric field with rm{e}^{ + i\omega t -ikr}(-k). The emitted light, i.e. the last interaction, has a frequency and wavevector which is the sum of the input frequencies and wavevectors, in accord to the energy conservation and phase matching condition respectively. Considering a molecular system with \left| 0 \right\rangle, \left| \alpha \right\rangle, and \left| \beta \right\rangle as the electronic ground, first and second electronic states respectively, with further assumption that the energy difference between \left| \alpha \right\rangle and \left| \beta \right\rangle is smaller than the spectral bandwidth of the excitation laser pulses, i.e., \left| \alpha \right| and \left| \beta \right| can be simultaneously excited. The first and the second electric field interactions create an absorption resulting in a population state. Therefore, the possible population states can be \left| \alpha \right\rangle \left\langle \alpha \right|, \left| \beta \right\rangle \left\langle \beta \right|, and \left| \alpha \right\rangle \left\langle \beta \right| together with their conjugated complexes. Obviously, \left| \alpha \right\rangle \left\langle \alpha \right| and \left| \beta \right\rangle \left\langle \beta \right| are the excited eigenstates while \left| \alpha \right\rangle \left\langle \beta \right| is the electronic coherent state, their corresponding double-sided Feynman diagrams are shown in FIG. 10. In measuring of the echo signal by scanning the waiting time T, static population decay kinetics are observed for either \left| \alpha \right\rangle \left\langle \alpha \right| or \left| \beta \right\rangle \left\langle \beta \right| only. However, for the coherent excited state, except for the static population relaxation, there is an additional coherent oscillation term as \exp[\pm(\omega _\alpha - \omega _\beta)T] defining the quantum beating signal overriding on the population relaxation envelope. An example of observation of the electronic quantum coherence is illustrated in FIG. 11 and FIG. 12 respectively. Therefore, whether a specific Liouville pathway is represented by the Feynman diagram leads to coherent oscillation can be easily determined by the generated population state after interaction with the first two pulses. Consequently, in 2DES the coherence is observed at the energy difference between the two eigenstates in the excited wave packet along the waiting time T axis, and the restriction arising from the collapse of the wave packet is thus circumvented.  FIG. 10 Feynman diagrams. (a) Diagrammatic representation of the interaction of an electromagnetic field with a material system, the interactions at t_1 and t_2 represent absorption and the interactions at t_3 and t_4 emissions at the ket and bra side respectively of the density matrix for a two-level system with \left| 0 \right\rangle and \left| 1 \right\rangle, time increases from bottom to top. The very last interaction (usually emission) corresponds to the coherently induced polarization in the ensemble. (b), (c) and (d) Feynman diagrams for excitonic system with \left| 0 \right\rangle, \left| \alpha \right\rangle and \left| \beta \right\rangle as the ground, first, and second electronic excited state respectively showing non-oscillation and oscillation properties.  FIG. 11 Double-sided Feynman diagrams describing the dominant contributions to electronic beats in 2DES spectra at diagonal and off-diagonal positions. Not all Feynman diagrams are included, the focus being on those important for the electronic beats. (a) and (b) for rephasing pathways while (c) and (d) for non-rephasing pathways (with modification to match the convention of Feynman diagram) [5].  FIG. 12 Typical example of 2DES showing the quantum beat signals. (a) The 2DES spectrum depicts the signal amplitude on an arcsinh scale (color scale, arbitrary units) plotted as the excitation frequency \omega_\tau and detection frequency \omega_t at a population delay time of T=100 fs for the antenna complex PC45 recorded at room temperature. (b) Amplitude of the cross-peaks at positions indicated in the 2D map with open squares [(\omega_\tau, \omega_t)=(2.185, 2.06) eV in black and (\omega_\tau, \omega_t)=(2.06, 2.185) eV in red] as a function of time T. The dashed lines interpolate the data points (solid circles). The solid line is a fit to a sum of damped sine functions revealing the quantum beats along T. For further information, see Refs.[6, 48]. B. Electronic versus vibronic coherence Two-dimensional spectroscopy has the capability to reveal the quantum coherence as the oscillatory behavior of the excitation dynamics of molecular systems. However, the situation is complicated regarding the exact origin of the observed coherent phenomena, i.e., what is actually being observed: excitonic or vibrational wave-packet motion or evidence of quantum transport. Especially, in some molecules and their aggregates, electronic transitions are coupled to various intra-and intermolecular vibrational modes, with the magnitudes of the resonant couplings J in excitonic aggregates being in the same range as those of vibrational energies. Thus, vibronic and excitonic systems show considerable spectroscopic similarities, and the presence of electronic and/or vibrational beats in the 2DES spectrum is expected. Indeed, similar spectral beats originating entirely from a high-energy vibrational wave-packet motion have been studied theoretically for weak electron-phonon coupling [49] and have been observed experimentally [50]. Number of groups have contributed to addressing this important issue in both theory and experiments [49-52]. Two generic model systems which exhibit distinct internal coherent dynamics have been compared [53]. One model system of pure electronic transition which shows similar spectroscopic properties but has completely different coherent internal dynamics without vibrations, is an excitonic dimer (ED) [54, 55]. It consists of two two-level chromophores (sites) with identical transition energies \varepsilon. The two sites are coupled by the inter-site resonance coupling J. As a result, the ED has one ground state \left| g \right\rangle, two single-exciton states \left| {e_1 } \right\rangle and \left| {e_2 } \right\rangle with energies \varepsilon _{e_1, e_2 } = \varepsilon$$ \mp$$J, respectively, and a single double-exciton state |f$$\rangle$ with energy $\varepsilon _f$=$2\varepsilon $$-$$ \Delta _{\textrm{binding}}$, where $\Delta _{\textrm{binding}}$ is the bi-exciton binding energy. The lowest transition frequency is denoted as $\omega _{eg}$=$\displaystyle{\frac{{\varepsilon -J}}{\hbar }}$, the transition frequency difference between $\left| {e_1} \right\rangle$ and $\left| {e_2 } \right\rangle$ is $\Delta \omega$=$\displaystyle{\frac{{2J}}{h}}$, the corresponding energy diagram, and 2DES pattern constructed from all the possible Liouville paths revealing oscillation (colored) and non-oscillation (black) in rephasing and non-rephasing 2DES are shown in FIG. 13 (a), (b), and (c) respectively.

 FIG. 13 Fourier maps and the Feynman diagrams for the excitonic dimer. (a) Energy levels. (b) and (c) Fourier maps for rephasing and non-rephasing configuration respectively. (d) Feynman diagrams leading to oscillation. $\omega$: positive frequency; $-\omega$: negative frequency; GB: ground-state bleaching; SE: stimulated emission. ESA: excited-state absorption. Black symbols for non-oscillating pathways, colored symbols for oscillating pathways.

The other model is the vibronic system of an isolated molecular electronic excitation represented by two electronic states, $\left| g \right\rangle$ and $\left| e \right\rangle$, which are coupled to a one-dimensional nuclear coordinate $q$ with a vibrational frequency energy of $h\nu$, assuming that $\nu$ keeps the same for all the vibrational excited states either in the electronic ground or excited state with its magnitude equal to that of $\Delta\omega$. This model is denoted as a displaced oscillator (DO) system (FIG. 14). Its energy diagram and 2DES pattern constructed from all the possible Liouville paths revealing oscillation and non-oscillation in 2DES are shown in FIG. 14 (a), (b) and (c) respectively.

 FIG. 14 Fourier map and the Feynman diagrams for the displaced oscillator model with coupling of a single vibrational mode. (a) Energy levels for displaced oscillator; (b) Fourier map for rephasing configuration. (c) Feynman diagrams leading to oscillation in rephasing configuration. $\nu$: positive frequency, $-\nu$: negative frequency, GB: ground-state bleaching, SE: stimulated emission. Black symbols for non-oscillating pathways, colored symbols for oscillating pathways.

To reveal oscillatory contributions in the ED and DO systems, all contributions into either oscillatory or static can be compared, shown in FIG. 13(b) and 14(b). As a function of $T$, the DO system has 8 oscillatory and 8 static configurations, which organize into six peaks, while the ED system has only 4 oscillatory and 8 static contributions which give four peaks. The net result is that the diagonal peaks in the rephasing and cross-peaks in the non-rephasing signals are non-oscillatory in the ED, while all peaks except the upper diagonal peak in rephasing signals are oscillatory in the DO. Thus significant differences in oscillatory peaks between ED and DO systems are found. The typical oscillatory pattern of 2DES spectrum appearing in the rephasing configuration of DO model with a single vibrational mode exhibits as a "chair", and is commonly referred as the "chair" pattern [56, 48].

Generally, the 2DES spectra are presented as the slices of pump and probe frequency projection at a fixed waiting time, and the coherence is given as oscillating amplitude at a certain point on the projection. To study the coherent oscillations, it is more straightforward to view the 2DES spectra along the beating frequency $\Delta\omega_T$ instead of $T$ by performing further Fourier transformation along $T$, leading to the so-called oscillation Fourier maps. They are constructed by applying the Fourier transformation on the real part of the rephasing spectrum $S_R (\omega _\tau, T, \omega _t)$ with respect to the delay time $T$ [48, 51, 57]: \begin{eqnarray} A(\omega _\tau, \Delta \omega _T, \omega _t) = \int_0^\infty {\textrm{d}T\textrm{e}^{i\Delta \omega _T T} {\mathop{\rm Re}\nolimits} } S_R (\omega _\tau, T, \omega _t) \end{eqnarray} In this way, at a given oscillation, all the contributions from the non-oscillation pathways and those of oscillation at a frequency other than the given one would be filtered, providing a powerful method to look into the coherent system having multiple beating frequencies. It should be noted that the Fourier maps for oscillation/non-oscillation contributions in 2DES spectra are similar whether it is projected at a fixed waiting $T$ or at a fixed beating frequency $\Delta\omega_T$, while only the latter can be used to trace the contribution from a given single beating frequency in the 2DES.

C. Single-vibrational mode versus multi-vibrational mode coupling in vibronic transition

Single vibrational-mode coupling to the vibronic transition is the simplest model. When multi-mode coupling is considered, the Hamiltonian as well as the corresponding sets of wave functions for a single-vibrational mode displaced oscillator coupled to the electronic states [58-60] needs to be extended to the multi-mode case [61]. Considering a multi-mode coupled displaced oscillator in a chromophore molecule with only two electronic states, $i.e$., the ground state $\left| g \right\rangle$ and the excited state $\left| e \right\rangle$, coupled to a collection of harmonic oscillators within the adiabatic approximation, the total Hamiltonian can be expressed as

 $\begin{eqnarray} H = H_g \left| g \right\rangle \left\langle g \right| + H_e \left| e \right\rangle \left\langle e \right| \end{eqnarray}$ (52)

where

 $\begin{eqnarray} H_g \hspace{-0.15cm}&=&\hspace{-0.15cm} \sum\limits_i {\left( {\frac{{p_i^2 }}{{2\mu _{gi} }} + \frac{1}{2}\mu _{gi} \omega _{gi}^2 q_i^2 } \right)} \end{eqnarray}$ (53)
 $\begin{eqnarray} H_e \hspace{-0.15cm}& =&\hspace{-0.15cm} h\omega _{eg} + \sum\limits_i {\left[{\frac{{p_i^2 }}{{2\mu _{gi} }} + \frac{1}{2}\mu _{ei} \omega _{ei}^2 (q_i-d_i )^2 } \right]} \end{eqnarray}$ (54)

Here $i$ represents the $i$th vibrational mode, $p_i$ and $q_i$ are the corresponding momentum and the position operator respectively. $\mu _{gi}$ ($\mu _{ei}$) is the reduced mass of the ground (excited) state, $\omega _{gi}$ ($\omega _{ei}$) is the frequency of the oscillator on the ground (excited) state, $d_i$ is the displacement of the origin on the excited-state potential relative to the ground state along the coordinate $q_i$, and $\omega_{eg}$ is the energy difference between the minima of the two electronic states. An assumption is made that the nuclear potential surfaces are well described as harmonic and that the excited energy surface is allowed to differ in curvature, giving rise to that the values of $\omega_{ei}$ and $\omega_{gi}$ can be different [62, 63]. The oscillator strength for the vibronic transition is determined by the Frank-Condon factor while each vibrational mode contributes with a weight given by its Huang-Rhys factor in terms of difference of the mode displacement of the initial and final nuclear configurations, $i.e$., the extent for the wavefuction overlap of nuclear part between the initial and the final states [64-66].

The simplest model of multi-mode coupling in displaced oscillator is the two vibtrational modes system with their vibrational frequencies denoted as $\nu$ and $\Omega$ respectively, which is schematically shown in FIG. 15(a), while the expected Fourier map for a given oscillation frequency $\nu$ with a coupled mode $\Omega$ is in FIG. 15(b) and the corresponding four Feynman diagrams leading to the mode coupling is in FIG. 15(c). The "chair" pattern indicated by the blue ellipses in FIG. 15(b) reflects the single vibrational mode coupling in DO model at the observing frequency $\nu$, $i.e$., the Fourier map is the projection at fixed value of $\nu$. Except for the expected chair pattern for the single mode, there are also two-mode coupled patterns appearing at the location of the rephasing frequency at $\omega _{eg} -\nu + \Omega$ excited at $\omega _{eg}$, $\omega _{eg}$+$\nu$ and $\omega _{eg}$+$\Omega$; rephasing frequency at $\omega _{eg}+\nu-\Omega$ excited at $\omega_0$; rephasing frequency at $\omega _{eg}-\nu$ excited at $\omega _{eg}+\Omega$; rephasing frequency at $\omega _{eg}-\Omega$ excited at $\omega _{eg}+\nu$. These indicate that the observed beating frequency $\nu$ can be excited at a vibrational frequency corresponding either to $\omega _{eg} + \nu$ or to $\omega _{eg} + \Omega$. The fact suggests that the vibrational frequencies $\nu$ and $\Omega$ are coupled to each other as shown in the red ellipses.

 FIG. 15 (a) Multi-vibrational mode-coupled displaced oscillator model with electronic transition energy $\omega _{eg}$ and two vibrational modes having a specific transition frequency of $\nu$ and $\Omega$ respectively; (b) Rephasing Fourier map for the vibrational mode-coupled displaced model at the oscillation frequency $\Delta \omega _T$=$\nu$ displaying the uncoupled pattern (blue ellipses) and the two-vibrational mode coupled pattern (red ellipses). The squares represent the beating appearing in the ground-state bleaching (GB) signals while the circles in the stimulated emission (SE) signals. The filled symbols for the beating with +$\nu$ frequency while the open ones for the $-\nu$ frequency. (c) Feynman diagrams of seven typical coupling cases having a beating frequency of $\pm\nu$. The subscript indicates ladder of the mode, and the superscript indicates the vibrational mode.

Obviously, since the detected oscillation signals originate from all the possible Liouville pathways including the excitonic and single-mode or multi-mode vibrationally coupled coherences, new methods either of apparatus or data analysis are called for clear distinguishing the electronic coherence from those of vibronic. It has been shown that the polarization modulated [67, 68] 2DES can be an effective way toward this purpose as successfully demonstrated in a recent experiment, where the oscillation maps extracted from the double-cross polarized schemes (relative polarizations of $\pi/4$, $-\pi/4$, $\pi/2$, and 0 for laser pulses 1 to 4) strongly suppresses signals of intramolecular vibrational coherences produced by the Franck-Condon excitation of vibrational wave packets [69]. On the other hand, 2D single-time-domain exciton perturbation spectroscopy (2D-STEPS) and 2D time-resolved interexciton perturbation spectroscopy (2D-TRIPS) are the new approaches to extract information from 2DES spectra. These time-domain optical spectroscopies probe vibrational motion by examining correlated dynamics of each exciton's spectral motion. 2D-STEPS measures how excited-state nuclear motions systematically perturb individually excitonic transitions after excitation, whereas 2D-TRIPS probes how nuclear motions spawned by exciting one exciton affect other unoccupied excitonic states [70].

Ⅶ. CONCLUSION

Quantum coherence manifests itself not only as a concept but also as a reality in many aspects. An emerging field of quantum computation involves coherence control and the quantum information storage, while 3PPE provides possible quantum information storage in time sequence. While in natural and artificial photosynthetic systems, it also has been suggested that coherent energy transport could be more efficient than classical transport. Development of advanced theories would provide more applicable feasibility while the new advances in the experimental methods and designs of novel quantum coherent systems would speed up the application of quantum coherence. Though quantum coherence was considered fragile, recent evidence of coherence in chemical and biological systems suggests that the phenomena are robust and can survive in the face of disorder and noise, suggesting that coherence can be used in complex chemical systems [71].

Ⅷ. ACKNOWLEDGMENTS

This work was supported by the National Natural Science Foundation of China (No.21227003, No.21433014, No.11721404). Mr. Xuan Leng and Rui-dan Zhu are kindly acknowledged for proof-reading and Figure preparation.

Reference
 [1] R. Lambert W., M. Felker P., and H. Zewail A., J. Chem. Phys. 75 , 5958 (1981). DOI:10.1063/1.442052 [2] H. Zewail A., J. Phys. Chem. 97 , 12427 (1993). DOI:10.1021/j100150a001 [3] Katsuki H., Chiba H., Girard B., Meier C., and Ohmori K., Science 311 , 1589 (2006). DOI:10.1126/science.1121240 [4] G. D. Reid and K. Wynne, Ultrafast Laser Technology and Spectroscopy, in Encyclopedia of Analytical Chemistry, Chichester: John Wiley & Sons Ltd., (2000). [5] S. Engel G., R. Calhoun T., L. Read E., K. Ahn T., Mančal T., C. Cheng Y., E. Blankenship R., and R. Fleming G., Nature 446 , 782 (2007). DOI:10.1038/nature05678 [6] Collini E., Y. Wong C., E. Wilk K., M. G. Curmi P., Brumer P., and D. Scholes G., Nature 463 , 644 (2010). DOI:10.1038/nature08811 [7] Panitchayangkoon G., Hayes D., A. Fransted K., R. Caram J., Harel E., Wen J., E. Blankenship R., and S. Engel G., Proc. Natl. Acad. Sci. USA 107 , 12766 (2010). DOI:10.1073/pnas.1005484107 [8] Romero E., and I. Novoderezhkin V., R. van Grondelle, Nature 543 , 355 (2017). [9] X. Weng Y., Physics (in Chinese) 39 , 331 (2010). [10] X. Weng Y., Physics (in Chinese) 36 , 817 (2007). [11] Collini E., and D. Scholes G., Science 323 , 369 (2009). DOI:10.1126/science.1164016 [12] Hao K., Moody G., Wu F., Dass C.K., Xu L., H. Chen C., Sun L., Y. Li M., J. Li L., H. MacDonald A., and Li X., Nature Phys. 12 , 677 (2016). DOI:10.1038/nphys3674 [13] L. Cowan M., P. Ogilvie J., and J. D. Miller R., Chem. Phys. Lett. 386 , 184 (2004). DOI:10.1016/j.cplett.2004.01.027 [14] P. Ogilvie J., J. Kubarych K., E. Arimondo in:, and Berman P.R., C. C. Lin (Eds.), Adv. At. Mol. Opt. Phys. 57 , 249 (2009). DOI:10.1016/S1049-250X(09)57005-X [15] M. G. Faeder S., and M. Jonas D., J. Phys. Chem. A 103 , 10489 (1999). [16] H. Zhang T., N. Borca C., Q. Li X., and T. Cundiff S., Optics Express 13 , 7432 (2005). DOI:10.1364/OPEX.13.007432 [17] D. Hybl J., W. Albrecht A., M. Gallagher Faeder S., and M. Jonas D., Chem. Phys. Lett. 297 , 307 (1998). DOI:10.1016/S0009-2614(98)01140-3 [18] D. Hybl J., A. Ferro A., and M. Jonas D., J. Chem. Phys. 115 , 6606 (2001). DOI:10.1063/1.1398579 [19] Brixner T., Stenger J., M. Vaswani H., Cho M., E. Blankenship R., and R. Fleming G., Nature 434 , 625 (2005). DOI:10.1038/nature03429 [20] I. Prokhorenko V., Halpin A., and J. D. Miller R., Opt. Express 17 , 9764 (2009). DOI:10.1364/OE.17.009764 [21] F. Tekavec P., A. Myers J., L. M. Lewis K., and P. Ogilvie J., Opt. Lett. 34 , 1390 (2009). DOI:10.1364/OL.34.001390 [22] T. Cundiff S., D. Bristow A., Siemens M., Li H., Moody G., Karaiskaj D., Dai X., and Zhang T., IEEE J. Sel. Top. Quant. 18 , 318 (2012). DOI:10.1109/JSTQE.2011.2123876 [23] Cho M., Chem. Rev. 108 , 1331 (2008). DOI:10.1021/cr078377b [24] Sen D., Current Science 107 , 203 (2014). [25] Fano U., Rev. Modern Phys. 29 , 74 (1957). DOI:10.1103/RevModPhys.29.74 [26] G. Dantus M., and Gross P., Encyclopedia of Applied Physics, 22 , 431 (1998). [27] J. Griffiths D., Introduction to Quantum Mechanics, Cambridge University Press (2016). [28] C. N. Lincoln, Ph. D Thesis, Swinburne University of Technology, (2007). [29] B. Williams R., Ph. D thesis, Cornell University (2001). [30] S. V. Poltavtsev, M. Salewski, Y. V. Kapitonov, I. A. Yugova, I. A. Akimov, C. Schneider, M. Kamp, S. Hoefling, D. R. Yakovlev, A. V. Kavokin, and M. Bayer, Phys. Rev. B 93, 121304(R) (2016). [31] Collini E., Chem. Soc. Rev. 42 , 4932 (2013). DOI:10.1039/c3cs35444j [32] Mukamel S., Principles of Nonlinear Optical Spectroscopy, Oxford University Press (1999). [33] Nilsson C., Lund Reports in Atomic Physic (1997). [34] Hamm P., and Zanni M., Concepts and Methods of 2D Infrared Spectroscopy, Cambridge University Press (2011). [35] H. Eberly J., Optical Resonance and Two-level Atoms, Dover (1987). [36] Polli D., Rivalta I., Nenov A., Weingart O., Garavelli M., and Cerullo G., Photochem. Photobiol. Sci. 14 , 213 (2015). DOI:10.1039/C4PP00370E [37] T. Fourkas J., Annu. Rev. Phys. Chem. 53 , 17 (2002). DOI:10.1146/annurev.physchem.53.082001.144216 [38] Zinth W., and Kaiser W., Ultrashort Laser Pulses, Springer (1988). [39] Mitsunaga M., Optical and Quantum Electronics 24 , 1137 (1992). DOI:10.1007/BF00620310 [40] Kroll S., and Elman U., Opt. Lett. 18 , 1834 (1993). DOI:10.1364/OL.18.001834 [41] D. Fuller F., P. Ogilvie J., and M. A. Johnson in:, T. J. Martinez (Eds.), Annu. Rev. Phys. Chem. 66 , 667 (2015). DOI:10.1146/annurev-physchem-040513-103623 [42] H. Cho M., Brixner T., Stiopkin I., Vaswani H., and R. Fleming G., J. Chin. Chem. Soc. 53 , 15 (2006). DOI:10.1002/jccs.v53.1 [43] Jonas D.M., Annu. Rev. Phys. Chem. 54 , 425 (2003). DOI:10.1146/annurev.physchem.54.011002.103907 [44] Brixner T., Mančal T., V. Stiopkin I., and R. Fleming G., J. Chem. Phys. 121 , 4221 (2004). DOI:10.1063/1.1776112 [45] X. Weng Y., and Chen H., Ultrafast Spectroscopy: Principles and Techniques (in Chinese), Beijing: Chemical Industry Press (2013). [46] Wynne K., Muller M., Brandt D., and D. W. Vanvoorst J., Chem. Phys. 125 , 211 (1988). DOI:10.1016/0301-0104(88)87075-7 [47] D. Fainberg B., and Narbaev V., J. Chem. Phys. 113 , 8113 (2000). DOI:10.1063/1.1315998 [48] B. Turner D., Dinshaw R., K. Lee K., S. Belsley M., E. Wilk K., M. G. Curmi P., and D. Scholes G., Phys. Chem. Chem. Phys. 14 , 4857 (2012). DOI:10.1039/c2cp23670b [49] Nemeth A., Milota F., Mancal T., Lukes V., Hauer J., F. Kauffmann H., and Sperling J., J. Chem. Phys. 132 , 184514 (2010). DOI:10.1063/1.3404404 [50] Christensson N., F. Kauffmann H., Pullerits T., and Mancal T., J. Phys. Chem. B 116 , 7449 (2012). DOI:10.1021/jp304649c [51] Butkus V., Zigmantas D., Abramavicius D., and Valkunas L., Chem. Phys. Lett. 587 , 93 (2013). DOI:10.1016/j.cplett.2013.09.043 [52] Mancal T., Nemeth A., Milota F., Lukes V., F. Kauffmann H., and Sperling J., J. Chem. Phys. 132 , 184515 (2010). DOI:10.1063/1.3404405 [53] Butkus V., Zigmantas D., Valkunas L., and Abramavicius D., Chem. Phys. Lett. 545 , 40 (2012). DOI:10.1016/j.cplett.2012.07.014 [54] Sakurai A., and Tanimura Y., J. Phys. Chem. A 115 , 4009 (2011). DOI:10.1021/jp1095618 [55] Xu J., Xu R., Abramavicius D., Zhang H., and Yan Y., Chin. J. Chem. Phys. 24 , 497 (2011). DOI:10.1088/1674-0068/24/05/497-506 [56] S. Senlik S., R. Policht V., and P. Ogivie J., J. Phys. Chem. Lett. 6 , 2413 (2015). DOI:10.1021/acs.jpclett.5b00861 [57] R. Calhoun T., S. Ginsberg N., S. Schlau-Cohen G., C. Cheng Y., Ballottari M., Bassi R., and R. Fleming G., J. Phys. Chem. B 113 , 16291 (2009). DOI:10.1021/jp908300c [58] Butkus V., Valkunas L., and Abramavicius D., J. Chem. Phys. 140 , 034306 (2014). DOI:10.1063/1.4861466 [59] Basinskaite E., Butkus V., Abramavicius D., and Valkunas L., Photosyn. Res. 121 , 95 (2014). DOI:10.1007/s11120-014-0002-z [60] P. Wang Y., Wang Z., and X. Weng Y., Chin. Sci. Bull. (Chinese Version) 57 , 2895 (2012). DOI:10.1360/972012-587 [61] F. Fidler A., and S. Engel G., J. Phys. Chem. A 117 , 9444 (2013). DOI:10.1021/jp311713x [62] Yue S., Wang Z., Leng X., D. Zhu R., L. Chen H., and X. Weng Y., Chem. Phys. Lett. 683 , 591 (2017). DOI:10.1016/j.cplett.2017.03.029 [63] Yue S., Wang Z., C. He X., B. Zhu G., and X. Weng Y., Chin. J. Chem. Phys. 28 , 509 (2015). DOI:10.1063/1674-0068/28/cjcp1506136 [64] Huang K., and Rhys A., Proc. R. Soc. Lond. A 204 , 406 (1950). DOI:10.1098/rspa.1950.0184 [65] Pieper J., Voigt J., and J. Small G., J. Phys. Chem. B 103 , 2319 (1999). DOI:10.1021/jp984460e [66] Ratsep M., Pajusalu M., and Freiberg A., Chem. Phys. Lett. 479 , 140 (2009). DOI:10.1016/j.cplett.2009.07.094 [67] S. Schlau-Cohen G., Ishizaki A., R. Calhoun T., S. Ginsberg N., Ballottari M., Bassi R., and R. Fleming G., Nature Chem. 4 , 389 (2012). DOI:10.1038/nchem.1303 [68] T. Zanni M., H. Ge N., S. Kim Y., and M. Hochstrasser R., Proc. Natl. Acad. Sci. USA 98 , 11265 (2001). DOI:10.1073/pnas.201412998 [69] Palecek D., Edlund P., Westenhoff S., and Zigmantas D., Sci. Advances 3 , e1603141 (2017). [70] S. Rolczynski B., Zheng H., P. Singh V., Navotnaya P., R. Ginzburg A., R. Caram J., Ashraf K., T. Gardiner A., H. Yeh S., Kais S., J. Cogdell R., and S. Engel G., Chem. 4 , 138 (2018). DOI:10.1016/j.chempr.2017.12.009 [71] D. Scholes G., R. Fleming G., Chen L.X., AspuruGuzik A., Buchleitner A., F. Coker D., S. Engel G., van Grondelle R., Ishizaki A., Jonas D.M., S. Lundeen J., K. McCusker J., Mukamel S., P. Ogilvie J., Olaya-Castro A., A. Ratner M., C. Spano F., B. Whaley K., and Zhu X., Nature 543 , 647 (2017). DOI:10.1038/nature21425