Sum-over-states vs quasiparticle pictures of coherent correlation spectroscopy of excitons in semiconductors; femtosecond analogues of multidimensional NMR

TL;DR

This paper compares the sum-over-states and quasiparticle approaches to understanding 2D correlation spectroscopy of excitons in semiconductors, providing a formal connection between the two methods using Green's functions.

Contribution

It derives Green's function expressions linking the sum-over-states and quasiparticle pictures for 2DCS signals in semiconductors, clarifying their relationship.

Findings

Established formal connection between the two theoretical pictures

Linked the approaches to Frenkel and vibrational excitons

Provided Green's function framework for 2DCS analysis

Abstract

Two-dimensional correlation spectroscopy (2DCS) based on the nonlinear optical response of excitons to sequences of ultrafast pulses, has the potential to provide some unique insights into carrier dynamics in semiconductors. The most prominent feature of 2DCS, cross peaks, can best be understood using a sum-over-states picture involving the many-body eigenstates. However, the optical response of semiconductors is usually calculated by solving truncated equations of motion for dynamical variables, which result in a quasiparticle picture. In this work we derive Green's function expressions for the four wave mixing signals generated in various phase-matching directions and use them to establish the connection between the two pictures. The formal connection with Frenkel excitons (hard-core bosons) and vibrational excitons (soft-core bosons) is pointed out.

Figures

Equations

$$E\left(\mathbf{r},t\right)=\sum _{j=1}^3{E}_j\left(\mathbf{r},t\right)=E^{+}(\mathbf{r},t)+E^{-}(\mathbf{r},t),$$

$$E^{+}(\mathbf{r},t)=\sum _{j=1}^3{\mathcal{E}}_j^{+}(t-{\tau }_j)e^{-i{\omega }_{j}t}{e}^{i{\mathbf{k}}_j\mathbf{r}},$$

$$E^{-}(\mathbf{r},t)=\sum _{j=1}^3{\mathcal{E}}_j^{-}(t-{\tau }_j)e^{i{\omega }_{j}t}{e}^{-i{\mathbf{k}}_j\mathbf{r}}.$$

$$\mathbf{P}(\mathbf{r},{\tau }_4)={\iiint }_{-\infty }^{\infty }\mathrm{d}{\tau }_3\mathrm{d}{\tau }_2\mathrm{d}{\tau }_1{\mathbf{S}}^{(SOS)}({\tau }_4,{\tau }_3,{\tau }_2,{\tau }_1)\mathbf{E}(\mathbf{r},{\tau }_3)\mathbf{E}(\mathbf{r},{\tau }_2)\mathbf{E}(\mathbf{r},{\tau }_1),$$

$${\mathbf{S}}^{(SOS)}({\tau }_4,{\tau }_3,{\tau }_2,{\tau }_1)$$

$$-\theta ({\tau }_{43})\theta ({\tau }_{42})\theta ({\tau }_{21})⟨\hat{\mathbf{\mu }}({\tau }_3)\hat{\mathbf{\mu }}({\tau }_4)\hat{\mathbf{\mu }}({\tau }_2)\hat{\mathbf{\mu }}({\tau }_1)⟩$$

$$+\theta ({\tau }_{42})\theta ({\tau }_{23})\theta ({\tau }_{41})⟨\hat{\mathbf{\mu }}({\tau }_3)\hat{\mathbf{\mu }}({\tau }_2)\hat{\mathbf{\mu }}({\tau }_4)\hat{\mathbf{\mu }}({\tau }_1)⟩$$

$$-\theta ({\tau }_{41})\theta ({\tau }_{12})\theta ({\tau }_{23})⟨\hat{\mathbf{\mu }}({\tau }_3)\hat{\mathbf{\mu }}({\tau }_2)\hat{\mathbf{\mu }}({\tau }_1)\hat{\mathbf{\mu }}({\tau }_4)⟩].$$

$$⟨\mathbf{\mu }({\tau }_4)\mathbf{\mu }({\tau }_3)\mathbf{\mu }({\tau }_2)\mathbf{\mu }({\tau }_1)⟩$$

$$\displaystyle=\sum_{a_{1},a_{2},a_{3}}\bm{\mu}_{ga_{3}}\bm{\mu}_{a_{3}a_{2}}\bm{\mu}_{a_{2}a_{1}}\bm{\mu}_{a_{1}g}e^{-i\big{[}(\varepsilon_{a_{3}}-\varepsilon_{g})\tau_{4}+(\varepsilon_{a_{2}}-\varepsilon_{a_{3}})\tau_{3}+(\varepsilon_{a_{1}}-\varepsilon_{a_{2}})\tau_{2}+(\varepsilon_{g}-\varepsilon_{a_{1}})\tau_{1}\big{]}}~{}.$$

$${\hat{\mathbf{\mu }}}^{+}$$

$${\hat{\mathbf{\mu }}}^{-}$$

$$H_I\left(t\right)=-{\hat{\mathbf{\mu }}}^{+}{\mathbf{E}}^{+}(\mathbf{r},\tau )-{\hat{\mathbf{\mu }}}^{-}{\mathbf{E}}^{-}(\mathbf{r},\tau )$$

$$⟨\hat{\mathbf{\mu }}\hat{\mathbf{\mu }}\hat{\mathbf{\mu }}\hat{\mathbf{\mu }}⟩=⟨{\hat{\mathbf{\mu }}}^{-}{\hat{\mathbf{\mu }}}^{+}{\hat{\mathbf{\mu }}}^{-}{\hat{\mathbf{\mu }}}^{+}⟩+⟨{\hat{\mathbf{\mu }}}^{-}{\hat{\mathbf{\mu }}}^{-}{\hat{\mathbf{\mu }}}^{+}{\hat{\mathbf{\mu }}}^{+}⟩.$$

$${\mathbf{S}}^{(SOS)}({\tau }_4,{\tau }_3,{\tau }_2,{\tau }_1)$$

$$+\theta ({\tau }_{43})\theta ({\tau }_{32})\theta ({\tau }_{21})⟨{\hat{\mathbf{\mu }}}^{-}({\tau }_4){\hat{\mathbf{\mu }}}^{+}({\tau }_3){\hat{\mathbf{\mu }}}^{-}({\tau }_2){\hat{\mathbf{\mu }}}^{+}({\tau }_1)⟩$$

$$-\theta ({\tau }_{43})\theta ({\tau }_{42})\theta ({\tau }_{21})⟨{\hat{\mathbf{\mu }}}^{-}({\tau }_3){\hat{\mathbf{\mu }}}^{-}({\tau }_4){\hat{\mathbf{\mu }}}^{+}({\tau }_2){\hat{\mathbf{\mu }}}^{+}({\tau }_1)⟩$$

$$+\theta ({\tau }_{42})\theta ({\tau }_{23})\theta ({\tau }_{41})⟨{\hat{\mathbf{\mu }}}^{-}({\tau }_3){\hat{\mathbf{\mu }}}^{+}({\tau }_2){\hat{\mathbf{\mu }}}^{-}({\tau }_4){\hat{\mathbf{\mu }}}^{+}({\tau }_1)⟩]$$

$$+c.c.$$

$$\theta ({\tau }_{43})\theta ({\tau }_{42})\theta ({\tau }_{21})=\theta ({\tau }_{42})\theta ({\tau }_{21})\theta ({\tau }_{13})+\theta ({\tau }_{42})\theta ({\tau }_{23})\theta ({\tau }_{31})+\theta ({\tau }_{43})\theta ({\tau }_{32})\theta ({\tau }_{21}).$$

$$\theta ({\tau }_{23})\theta ({\tau }_{42})\theta ({\tau }_{41})=\theta ({\tau }_{42})\theta ({\tau }_{21})\theta ({\tau }_{13})+\theta ({\tau }_{42})\theta ({\tau }_{23})\theta ({\tau }_{31})+\theta ({\tau }_{41})\theta ({\tau }_{12})\theta ({\tau }_{23})$$

$${\mathbf{S}}^{(SOS)}({\tau }_4,{\tau }_3,{\tau }_2,{\tau }_1)$$

$$\displaystyle+\theta(\tau_{43})\theta(\tau_{32})\theta(\tau_{21})\big{\langle}\hat{\bm{\mu}}^{-}(\tau_{4})\hat{\bm{\mu}}^{+}(\tau_{3})\hat{\bm{\mu}}^{-}(\tau_{2})\hat{\bm{\mu}}^{+}(\tau_{1})\big{\rangle}$$

$$\displaystyle-\theta(\tau_{43})\theta(\tau_{32})\theta(\tau_{21})\big{\langle}\hat{\bm{\mu}}^{-}(\tau_{3})\hat{\bm{\mu}}^{-}(\tau_{4})\hat{\bm{\mu}}^{+}(\tau_{2})\hat{\bm{\mu}}^{+}(\tau_{1})\big{\rangle}$$

$$\displaystyle-\theta(\tau_{42})\theta(\tau_{23})\theta(\tau_{31})\big{\langle}\hat{\bm{\mu}}^{-}(\tau_{3})\hat{\bm{\mu}}^{-}(\tau_{4})\hat{\bm{\mu}}^{+}(\tau_{2})\hat{\bm{\mu}}^{+}(\tau_{1})\big{\rangle}$$

$$\displaystyle-\theta(\tau_{42})\theta(\tau_{21})\theta(\tau_{13})\big{\langle}\hat{\bm{\mu}}^{-}(\tau_{3})\hat{\bm{\mu}}^{-}(\tau_{4})\hat{\bm{\mu}}^{+}(\tau_{2})\hat{\bm{\mu}}^{+}(\tau_{1})\big{\rangle}$$

$$\displaystyle+\theta(\tau_{41})\theta(\tau_{12})\theta(\tau_{23})\big{\langle}\hat{\bm{\mu}}^{-}(\tau_{3})\hat{\bm{\mu}}^{+}(\tau_{2})\hat{\bm{\mu}}^{-}(\tau_{4})\hat{\bm{\mu}}^{+}(\tau_{1})\big{\rangle}$$

$$\displaystyle+\theta(\tau_{42})\theta(\tau_{23})\theta(\tau_{31})\big{\langle}\hat{\bm{\mu}}^{-}(\tau_{3})\hat{\bm{\mu}}^{+}(\tau_{2})\hat{\bm{\mu}}^{-}(\tau_{4})\hat{\bm{\mu}}^{+}(\tau_{1})\big{\rangle}$$

$$\displaystyle\left.+\theta(\tau_{42})\theta(\tau_{21})\theta(\tau_{13})\big{\langle}\hat{\bm{\mu}}^{-}(\tau_{3})\hat{\bm{\mu}}^{+}(\tau_{2})\hat{\bm{\mu}}^{-}(\tau_{4})\hat{\bm{\mu}}^{+}(\tau_{1})\big{\rangle}\right]$$

$$+c.c.$$

$$\mathbf{P}(\mathbf{r},{\tau }_4)={\iiint }_0^{+\infty }d{t}_{3}d{t}_{2}d{t}_1{\mathbf{S}}^{(3)}(t_3,t_2,t_1){\mathbf{E}}_3(\mathbf{r},{\tau }_4-t_3){\mathbf{E}}_2(\mathbf{r},{\tau }_4-t_2-t_3){\mathbf{E}}_1(\mathbf{r},{\tau }_4-t_1-t_2-t_3),$$

$$\mathbf{P}(\mathbf{r},{\tau }_4)={\mathbf{S}}^{(SOS)}(t_3,t_2,t_1){\mathcal{E}}_3^{{\lambda }_3}({\tau }_4-t_3){\mathcal{E}}_2^{{\lambda }_2}({\tau }_4-t_3-t_2){\mathcal{E}}_1^{{\lambda }_1}({\tau }_4-t_3-t_2-t_1)\times$$

$$e^{i({\lambda }_1{\mathbf{k}}_1+{\lambda }_2{\mathbf{k}}_2+{\lambda }_3{\mathbf{k}}_3)\mathbf{r}}{e}^{-i({\lambda }_1{\omega }_1+{\lambda }_2{\omega }_2+{\lambda }_3{\omega }_3){\tau }_4}{e}^{i({\lambda }_1{\omega }_1+{\lambda }_2{\omega }_2+{\lambda }_3{\omega }_3)t_3}{e}^{i({\lambda }_1{\omega }_1+{\lambda }_2{\omega }_2)t_2}{e}^{i{\lambda }_1{\omega }_1{t}_1}$$

$$\mathbf{P}\left(\mathbf{r},{\tau }_4\right)=\sum _{s=1}^4\mathbf{P}\left({\mathbf{k}}_s,{\omega }_s\right)e^{i{\mathbf{k}}_s\mathbf{r}-i{\omega }_s{\tau }_4}+c.c.$$

$$\mathbf{P}\left({\mathbf{k}}_s,{\omega }_s\right)={\mathbf{S}}_s^{(SOS)}\left(t_3,t_2,t_1\right){\mathcal{E}}_3^{{\lambda }_3}{\mathcal{E}}_2^{{\lambda }_2}{\mathcal{E}}_1^{{\lambda }_1}.$$

$${\mathbf{S}}^{(SOS)}(t_3,t_2,t_1)={\mathbf{S}}_I^{(SOS)}(t_3,t_2,t_1)+{\mathbf{S}}_{II}^{(SOS)}(t_3,t_2,t_1)+{\mathbf{S}}_{III}^{(SOS)}(t_3,t_2,t_1).$$

$${\mathbf{S}}_I^{(SOS)}$$

$$+⟨{\hat{\mathbf{\mu }}}^{-}(0){\hat{\mathbf{\mu }}}^{+}(t_1){\hat{\mathbf{\mu }}}^{-}(t_1+t_2+t_3){\hat{\mathbf{\mu }}}^{+}(t_1+t_2)⟩$$

$$-⟨{\hat{\mathbf{\mu }}}^{-}(0){\hat{\mathbf{\mu }}}^{-}(t_1+t_2+t_3){\hat{\mathbf{\mu }}}^{+}(t_1+t_2){\hat{\mathbf{\mu }}}^{+}(t_1)⟩].$$

$${\mathbf{S}}_{II}^{(SOS)}$$

Full text

I Introduction

Exciton models are widely used to describe the linear and nonlinear optical properties of many types of systems, including semiconductor nanostructures (quantum wells, dots and wires), molecular aggregates and crystals,Davydov (1962); Mukamel (1993); Chernyak et al. (1998a) as well as vibrations in proteins.Zhuang et al. (2005); Mukamel and Abramavicius (2004) In semiconductors, nonlinear optical experiments reveal a wealth of interesting phenomena.Chemla and Shah (2001); Kaindl et al. (2001); Shih et al. (2005); Shah (1999); Meier et al. (2006) For instance, such experiments provide information about many-exciton states such as biexcitons, their interactions, relaxation and dissociation.Chemla (2001); Adachi et al. (2004); Rice et al. (2005); Danckwerts et al. (2006)

The introduction of multidimensional techniques had revolutionized NMR in the seventiesAue et al. (1976) and established it as a powerful tool for studying complex systems and identifying specific structural and dynamical correlations.Ernst et al. (1989) In such experiments the system is subjected to a sequence of well separated pulses. Correlation plots of the signals vs. two (or more) time delay periods then provide multidimensional spectroscopic windows into the system. The correlated dynamics of spins carefully prepared by the pulse sequence is very sensitive to their interactions. Analysis of these correlation plots then provides a powerful probe for molecular geometries and dynamical correlations. These techniques were recently extended to the infrared and the visible regime and were shown to be very useful for Frenkel excitons in molecular systems.Mukamel (2000); Tanimura and Mukamel (1993); Mukamel and Abramavicius (2004); Zhang et al. (1999) There are some differences between the optical and the NMR techniques. NMR uses strong saturating fields whereas optical pulses are most effective in the weak field regime. NMR signals are essentially isotropic in space whereas coherent optical signals are generated in well defined (phase-matching) directions. These differences were explored in detail in Refs. Scheurer and Mukamel, 2002a; Scheurer and Mukamel, 2001, 2002b. Nevertheless the NMR and optical techniques are conceptually similar and many ideas of pulse sequences developed in NMR may be adopted in the optical regime, where the millisecond NMR time-scale is pushed to the femtosecond regime. The same ideas may be extended to study interband and intersuband excitations in semiconductors.Yang et al. (2006); Kuznetsova et al. (2007); Erementchouk et al. (2006); Li et al. (2006); Borca et al. (2005) Multidimensional analysis of the nonlinear optical response of semiconductors to sequences of femtosecond pulses could provide a novel probe for many-body interactions. In a recent workYang et al. (2006) on semiconductor Quantum Wells, 2D correlation spectra from three 3rd order optical techniques have been calculated. The unique character of 2D spectroscopy allowed to easily recognize and classify features due to different types of biexcitons. Such features are sometimes difficult to separate in the usual one-dimensional mode of displaying non-linear spectra, due to the strong line broadening and the highly congested exciton spectra.

Two types of approaches have been traditionally used towards modeling the nonlinear optical response of excitonic systems. The first is based on the many-body eigenstates obtained by exact diagonalization of the Hamiltonian.Mukamel (1995) Sum-over-states (SOS) expressions can then be derived for the nonlinear response functions and optical signals. This method is practical in many applications to electronic and vibrational Frenkel excitons in moleculesBrixner et al. (2005); Fang et al. (2004); Abramavicius et al. (2004) and allows clear identification and classification of possible single- and multi-photon resonances. Calculating the eigenstates is a serious computational bottleneck in extended structures. For an $N$ site tight-binding Frenkel-exciton model the number of single and two-exciton states scales as $\sim N$ and $\sim N^{2}$ respectively. For Wannier excitons in semiconductors these scalings are $\sim N^{2}$ and $\sim N^{4}$, making the simulations prohibitively expensive. This is why the approach is not widely used for electron-hole excitations in semiconductors. Instead, one adopts a second strategy, which describes the response in terms of quasiparticles (QP), and the many-particle eigenstates are never calculated. Spano and Mukamel (1989); Leegwater and Mukamel (1992); Mukamel (1995); Chernyak and Mukamel (1996); Spano and Mukamel (1991); Mukamel (1993) Calculations are performed by solving equations of motion for microscopic coherences, which are coupled to other dynamical variables. Even for a simple system such as a single semiconductor quantum well, solving the equations numerically to create a 2D map of a nonlinear response function is computationally expensive,Yang et al. (2006) since these equations must be solved repeatedly for different pulse delays. Only after obtaining the optical signal on a 2D time grid, a Fourier transform can be performed to get the 2DCS. Apart from direct, numerical solutions of equations of motion Koch et al. (1999); Weiser et al. (2000) there exist other theoretical approaches to exciton correlation effects, such as memory kernel representation Östreich et al. (1995); Axt et al. (2001) or Coupled Cluster Expansion for doped semiconductors.Primozich et al. (2000); Shahbazyan et al. (2000)

In this paper we derive closed expressions for 2DCS of semiconductors by solving the Nonlinear Exciton Equations (NEE)Axt and Mukamel (1998); Chernyak et al. (1998a) for the third order response. Both time-ordered and non-ordered forms of the response function which represent time and frequency domain techniques, respectively, are derived. Our QP expressions for the response are given in terms of the single exciton Green’s function and the exciton scattering matrix. The SOS response functions, in contrast, are expressed in terms of many-exciton eigenstates. Even though the response functions calculated using both techniques must be identical, the relation between the two pictures is not obvious. The expressions look very different and it is not possible to see their equivalence by a simple inspection. The SOS expressions contain large terms, which grow with system size and have opposite signs, thus they almost cancel. This complicates their numerical implementation. In contrast these cancellations are built-in from the outset in the QP approach, which uses a harmonic reference system. The nonlinearities are then attributed to exciton-exciton scattering which is absent in the harmonic reference system. The second goal of this paper is to show precisely how the two pictures of many-body correlations are connected. We write down the SOS expressions using the Keldysh loop and then derive the QP expressions directly from the SOS ones. This provides a time-domain interpretation for the interference effects. The SOS and the QP expressions provide complementary views into the origin of features seen in 2D spectrograms.

In Sec. II we present the SOS expressions for the third order response obtained from time-dependent perturbation theory. Their QP counterparts are derived in Sec. III. We use the method developed in Refs. Chernyak and Mukamel, 1996; Chernyak et al., 1998a to transform the Hamiltonian to a form typical for interacting oscillators. The starting many-electron Hamiltonian can be written in an ab-initio,Oszwałdowski et al. (2005) tight-binding Sieh et al. (1999a) or a $\mathbf{k}\cdot\mathbf{p}$ basis. One of the key results of this paper, i.e., the equivalence of the SOS and QP pictures is proven in Sec. IV. In Sec. V we derive closed expressions for 2D correlation signals. The QP approach provides a unified description for electron-hole excitations in semiconductors as well as to Frenkel excitons in molecular aggregates (Paulions) and anharmonic vibrations (bosons), which are described by the same general Hamiltonian. QP formulae for nonlinear response have been derived previously along similar lines for Frenkel excitons. This connection is shown in Appendix F. In the last Section (VI) we discuss the results.

II Sum-over-states expressions for the time-ordered nonlinear response

We consider a 4 wave-mixing experiment performed with three femtosecond laser pulses (Fig. 1). The optical electric field is:

[TABLE]

[TABLE]

[TABLE]

The $j$-th pulse is centered at $\tau_{j}$, has an envelope $\mathcal{E}_{j}(t-\tau_{j})$, carrier frequency $\omega_{j}$, and wavevector $\mathbf{k}_{j}$. $E^{+}$ ($E^{-}$) denotes the positive (negative) frequency part of the field, and $\mathcal{E}_{j}^{-}=\left(\mathcal{E}_{j}^{+}\right)^{\ast}$. The induced polarization in the system is recorded as a function of time-delays between pulses.

Assuming the dipole interaction with the optical field $\hat{H}_{I}={\hat{\bm{\mu}}}\cdot\bm{E}(\bm{r},\tau)$, where ${\hat{\bm{\mu}}}$ is the dipole operator, the third-order contribution to the system’s polarization can be written as

[TABLE]

where the response function $\bm{S}^{(SOS)}$, which connects the induced polarization with the laser field envelopes, is given by (throughout this paper we set $\hbar=1$):

[TABLE]

We shall use double-sided Feynman diagrams to represent the time ordering of various interactions.Mukamel (2005) The four terms in Eq. (5) are represented by diagrams a, b, c, d shown on Figure 2. These diagrams should be read starting at the bottom left and proceeding along the loop, clockwise, as indicated by the arrows. The $\tau_{i}$ variables are ordered on the Keldysh-Schwinger loop, but not necessarily in real (physical) time. $\tau_{i}$ in diagrams (a) and (d) are also ordered in real time. This is not the case for diagrams (b) and (c): in (b) $\tau_{3}$ can come either before or after $\tau_{1}$ and $\tau_{2}$, whereas in (c) $\tau_{1}$ can come either before or after $\tau_{3}$ and $\tau_{2}$.

If the eigenstates $|a\rangle$ and eigenvalues $\varepsilon_{a}$ of the system are known, Eq. (5) may be expanded in terms of the corresponding matrix elements:

[TABLE]

So far we considered a general multilevel system. We next turn to the response of excitons, where the energy levels form manifolds, classified by the number of excitons: the ground state $\left(g\right)$, single exciton $\left(e\right)$, two-exciton $\left(f\right)$ (or biexciton), etc. (Fig. 3). We shall assume that the dipole operator can only create and annihilate a single exciton at a time. Only the single and the two-exciton states then contribute to the third order signals. We further partition the dipole operator as $\hat{\bm{\mu}}=\hat{\bm{\mu}}^{+}+\hat{\bm{\mu}}^{-}$, where $\hat{\bm{\mu}}^{+}$ is the positive frequency part which induces upward $g$ to $e$ and $e$ to $f$ transitions, while its Hermitian conjugate $\hat{\bm{\mu}}^{-}$ (the negative frequency part) induces the opposite transitions. We thus write

[TABLE]

Invoking the rotating-wave approximation (RWA), we neglect all terms where at least one of the transitions is not in resonance with one of the incident carrier frequencies. The system-field interaction term then becomes

[TABLE]

Each correlation function in Eq. (5) will split into $2^{4}=16$ terms upon substituting $\hat{\bm{\mu}}=\hat{\bm{\mu}}^{+}+\hat{\bm{\mu}}^{-}.$ Assuming that the system is initially in the ground state, only two of these contributions are non-zero

[TABLE]

Substitution of Eq. (7) into Eq. (5) gives

[TABLE]

The four terms represented by the diagrams in Fig. 4 were obtained by taking $\hat{\bm{\mu}}(\tau_{4})=\hat{\bm{\mu}}^{-}(\tau_{4})$ for the last interaction, $\hat{\bm{\mu}}(\tau_{4})=\hat{\bm{\mu}}^{+}(\tau_{4})$ gives the complex conjugates. Hereafter left/right direction of the arrows corresponds to $\hat{\bm{\mu}}^{-}/\hat{\bm{\mu}}^{+}$ in Eq.(8). Note that time-reversal symmetry implies $\left\langle\hat{\bm{\mu}}^{-}(\tau_{4})\hat{\bm{\mu}}^{-}(\tau_{3})\hat{\bm{\mu}}^{-}(\tau_{2})\hat{\bm{\mu}}^{-}(\tau_{1})\right\rangle^{\ast}=\left\langle\hat{\bm{\mu}}^{+}(\tau_{1})\hat{\bm{\mu}}^{+}(\tau_{2})\hat{\bm{\mu}}^{+}(\tau_{3})\hat{\bm{\mu}}^{+}(\tau_{4})\right\rangle$. If the pulse envelopes are much shorter than their delays, the system is forced to interact sequentially first with pulse $\bm{k}_{1}$, then $\bm{k}_{2}$ and finally $\bm{k}_{3}$. This means that in the integral of Eq. (4) one must replace $E(\bm{r},\tau_{j})$ with one of the $\mathcal{E}_{j}$, depending on the time-ordering of the integration variables in real (physical) time. We note that the first and the second terms in Eq. (8) impose a full time ordering of the integration variables while the third and the fourth terms do not. Term (b) is only partially time ordered. Depending on the position of $\tau_{3}$ relative to the $\tau_{1}<\tau_{2}<\tau_{4}$ sequence, the diagram can be separated into three fully time ordered terms: $\tau_{3}<\tau_{1}$, $\tau_{1}<\tau_{3}<\tau_{2}$ or $\tau_{2}<\tau_{3}<\tau_{4}$. Formally we do that by separating the product of step functions as follows:

[TABLE]

Using this relation, diagram (b) of Fig. 4 is split into (b3), (b2) and (b1) as shown in the first line of Fig. 5. The interactions on the l.h.s. of this diagrammatic equation are ordered on the loop. On the other hand, the arrows in the open, double-sided diagrams on the r.h.s. are ordered in real (physical) time. All diagrams on the r.h.s. are obtained from (b) by moving the arrows while preserving their order along the loop (but not in physical time!). Similarly we write for term (c)

[TABLE]

and the diagram is split into (c2), (c3), (c1). $\bm{S}^{(3)}$ can now be recast in the fully time-ordered form

[TABLE]

The labels on the right correspond to the various diagrams shown in Figs. (4) and (5).

Once split into fully time-ordered contributions, it is convenient to change the integration variables in Eq. (4) from $\tau_{4},\tau_{3},\tau_{2},\tau_{1}$ that label the actual interaction times with the fields, to the three delays $t_{3},t_{2},t_{1}$ between successive interactions. Note that the correlation functions are invariant to time translation $\left\langle\hat{\bm{\mu}}(\tau_{4}-\tau)\hat{\bm{\mu}}(\tau_{3}-\tau)\hat{\bm{\mu}}(\tau_{2}-\tau)\hat{\bm{\mu}}(\tau_{1}-\tau)\right\rangle=\left\langle\hat{\bm{\mu}}(\tau_{4})\hat{\bm{\mu}}(\tau_{3})\hat{\bm{\mu}}(\tau_{2})\hat{\bm{\mu}}(\tau_{1})\right\rangle$. Eq. (4) thus assumes the form

[TABLE]

where $t_{1}=\tau_{2}-\tau_{1},$ $t_{2}=\tau_{3}-\tau_{2},$ $t_{3}=\tau_{4}-\tau_{3}.$ In the impulsive limit, where all pulses are shorter than all system’s response time scales, we can substitute Eqs. (1)-(3) in Eq. (10) and eliminate the time integrations. This gives

[TABLE]

The polarization is created along 8 possible directions $\bm{k}_{s}=\lambda_{1}\bm{k}_{1}+\lambda_{2}\bm{k}_{2}+\lambda_{3}\bm{k}_{3}$ with $\lambda_{i}=\pm 1$

[TABLE]

where

[TABLE]

$\bm{k}_{1}+\bm{k}_{2}+\bm{k}_{3}$ vanishes for the assumed dipole selection rules in our model. Since $\bm{P}\left(-\bm{k}_{s},-\omega_{s}\right)=\bm{P}^{\ast}\left(\bm{k}_{s},\omega_{s}\right)$, we are left with three independent combinations $\bm{k}_{I}\equiv-\bm{k}_{1}+\bm{k}_{2}+\bm{k}_{3}$, $\bm{k}_{II}\equiv+\bm{k}_{1}-\bm{k}_{2}+\bm{k}_{3}$, $\bm{k}_{III}\equiv+\bm{k}_{1}+\bm{k}_{2}-\bm{k}_{3}$:

[TABLE]

We can classify the diagrams in Fig. 5 according to the directions of the arrows: arrow pointing to the right (left) represents $+\bm{k}$ ($-\bm{k}$), arrows are read from the bottom up on either side. We obtain for $\bm{k}_{I}$ (Fig. 6)

[TABLE]

For the $\bm{k}_{II}$ technique we similarly have (Fig. 7):

[TABLE]

Finally $\bm{k}_{III}$ is given by (Fig. 8):

[TABLE]

Each term is labelled according to Eq. (9). Eqs. (13-15) can be used to express the third order SOS response in terms of transition dipoles, system frequencies and dephasing rates (see App. E and Sec. V).

In the next section we employ the EOM approach to derive the alternative QP expressions for these signals. These will then be connected with the current SOS expressions in Section IV.

III Quasiparticle expressions for Wannier excitons in semiconductors

Interband transitions in semiconductors may be described by the two-band many-electron Hamiltonian:Axt and Mukamel (1998); Haug and Koch (2004)

[TABLE]

with the single-particle part

[TABLE]

where $c^{\dagger}$ create electrons and $d^{\dagger}$ create holes. The Coulomb interaction is:

[TABLE]

while

[TABLE]

is the dipole interaction with light, and the optical electric field $E$ will be treated as a scalar for simplicity. $\hat{H}_{T}$ can describe both bulk and low-dimensional semiconductor systems. All the steps in this Section are independent of the single-electron basis used. $\hat{H}_{0}$ would be diagonal in the basis of the system’s single-particle eigenstates, i.e., $t_{m1,n1}^{\left(i\right)}=\varepsilon_{m_{1}}^{\left(i\right)}\delta_{m_{1}n_{1}}$. In this paper we focus on the coherent response and we neglect coupling with phonons, which would result in additional, relevant dynamical variables and new contributions to the response function. Chernyak et al. (1998a); Mukamel (1993) The SOS and QP pictures should be equivalent also when dephasing is included. In that case, however, the theory becomes more complicated. For the sake of simplicity and transparency we restrict the following analysis to the coherent response, where we do not include phonons explicitly. Dephasing effects, necessary for a realistic description, will be simply introduced by adding imaginary parts to excitonic frequencies.

To introduce the exciton representation we define electron-hole operators:Chernyak and Mukamel (1996)

[TABLE]

where we have employed shorthand notation for pairs of indices: $m\equiv\left(m_{1},m_{2}\right)$. Using these operators we construct an effective Hamiltonian $\hat{H}$ (see App. A):

[TABLE]

The Hamiltonians $\hat{H}$ and $\hat{H}_{T}$ are equivalent in the single and double excitations subspace, which is relevant for the response to third order in $E$.Mukamel (1995) This transformation from fermion to exciton variables is crucial for our approach, since it allows us to view the electronic degrees of freedom as a system of coupled oscillators. The parameters of the transformed Hamiltonian $\hat{H}$ are given by:

[TABLE]

The commutation relations for the $\hat{B}$ operators can be obtained using the elementary fermion anticommutators: $\left[c_{m_{1}}^{\dagger},c_{k_{1}}\right]_{+}=\delta_{m_{1},k_{1}}$. Within the subspace of $\left|0\right\rangle$ and $\hat{B}_{i}^{\dagger}\left|0\right\rangle$ states (i.e., the ground state and single excitations), we get Chernyak and Mukamel (1996)

[TABLE]

where $\delta_{mn}=\delta_{m_{1}n_{1}}\delta_{m_{2}n_{2}}$ and

[TABLE]

Eqs. (19) and (20) are obtained in a similar way to (17) and (18).  Terms with additional $\hat{B}_{i}^{\dagger}\hat{B}_{j}$ pairs (e.g. $\hat{B}^{\dagger}\hat{B}^{\dagger}\hat{B}\hat{B}$) are neglected in (19), because they would introduce corrections higher than $O\left(E^{3}\right)$ to the nonlinear response. Note the symmetry $\mathcal{P}_{mnpq}=\mathcal{P}_{mnqp}$.

Using Eqs. (17) and (19) we obtain the nonlinear exciton equations (see Appendix B) for single-exciton variables $\left\langle\hat{B}_{m}\right\rangle$:Spano and Mukamel (1991); Mukamel (1993); Chernyak et al. (1998b); Mukamel (2000)

[TABLE]

where $V$ is given by

[TABLE]

Here $Y_{mn}\equiv\left\langle\hat{B}_{m}\hat{B}_{n}\right\rangle$ are two-exciton variables. The Heisenberg equations give:

[TABLE]

Calculating the optical response by numerical integration of these equationsSieh et al. (1999b); Yang et al. (2006) is straightforward but numerically expensive. An alternative, more tractable approach, which further provides a better insight into the nature of the response, is to integrate the equations formally using one-exciton Green’s functions $G\left(t\right)$ and exciton scattering matrix $\Gamma\left(t\right)$. The scattering matrix depends on quasiparticle statistics through the $\mathcal{P}$ matrix (Eqs. 49, 50) as well as on exciton-exciton coupling. This results in closed quasiparticle expressions for the 3rd order contributions $S_{I}$, $S_{II}$ and $S_{III}$ to the response function (for details see Appendix C and Ref. Abramavicius and Mukamel, 2005)

[TABLE]

where $\tau_{43}=\tau_{4}-\tau_{3},$ etc. and $G_{mn}\left(t\right)=-i\theta\left(t\right)\left[\exp\left(-iht\right)\right]_{mn}$.

The response functions for the other phase-matching directions can be derived along the same lines. We get

[TABLE]

and:

[TABLE]

Just as in the SOS case, time translation symmetry implies that these response functions only depend on the three pulse delays $t_{3},t_{2},t_{1}$. Eqs. (24-26) will be used next to connect the QP and the SOS pictures.

IV Connecting the sum-over-states and the quasiparticle pictures

We first recast Eqs. (13-15) using Green’s functions (in all expressions $t_{1}>0$, $t_{2}>0$, $t_{3}>0$):

[TABLE]

Here $\hat{G}(t)\equiv-i\theta(t)\exp(-i\hat{H}t)$ and $\hat{G}^{{\dagger}}(t)\equiv+i\theta(t)\exp(i\hat{H}t)$ represent the retarded and the advanced Green’s function respectively; $\mathbb{\hat{G}}$, $\hat{G}$ and $\mathcal{\hat{G}}$ describe the evolution within the ground-state, single-exciton and double-exciton blocks of the Hamiltonian (Eq. 17) respectively. We also set the ground state energy $\varepsilon_{g}$ to zero.

Our goal is to show the equivalence of the QP and SOS pictures by deriving Eqs. (24-26) from Eqs. (27-29). To that end we adopt a harmonic reference system of noninteracting quasiparticle and expand the SOS response in anharmonicities. Harmonic oscillators are linear, and their nonlinear response vanishes identically Spano and Mukamel (1989); Kuhn et al. (1996); Mukamel (1993, 1995), as can be easily seen from the Heisenberg equations of motion. This means that the various Liouville space pathways for all nonlinear response function interfere destructively. Exploiting this property in the following derivation, we show that the quasiparticle physical picture has built-in cancellations in the reference harmonic system.

We shall use the Dyson equation for the two particle Green’s function, also known as the Bethe-Salpeter equationLeegwater and Mukamel (1992)

[TABLE]

or in the time domain:

[TABLE]

$\mathcal{\hat{G}}_{0}$ is taken to be the Green’s function of a doubly excited, harmonic system. It can be factorized into the product of a single-exciton Green’s functions

[TABLE]

The exciton scattering matrix $\Gamma$ is defined by Eq. (31).

Let us start with the $\bm{S}_{I}$ technique and show the equivalence of Eq. (27) to (24). The second and third terms of Eq. (27) (diagrams (c1) and (c3) in Fig. 6) are purely harmonic, independent on the quasiparticle interactions. This is a direct result of the ordering of $\bm{\mu}^{\pm}$, whereby the system only evolves in the ground and first excited state. Exciton-exciton interactions influence the evolution only in the second excited manifold. The first term in Eq. (31), i.e. $\mathcal{\hat{G}}_{0}$, represents harmonic evolution in the two-exciton manifold. Thus the first term of Eq. (27) with $\mathcal{\hat{G}}$ replaced by $\mathcal{\hat{G}}_{0}$ must cancel the other two terms, because the nonlinear response of a harmonic system vanishes. Substituting the second term from Eq. (31) in Eq. (27) we obtain a single term for $\bm{S}_{I}$:

[TABLE]

The equivalence of the Eqs. (32) and (24) can be directly seen using the diagrams shown in Fig. (9). In these diagrams the scattering matrix $\Gamma$ is represented by dashed regions. Note that $G\left(t\right)G^{\dagger}\left(t\right)=\theta\left(t\right)\exp(-iht)\exp(iht)=\theta\left(t\right)$. The QP diagram in Fig. (9) is obtained from the SOS one by changing the integration variables $\tau^{\prime}=t_{3}-\tau_{s}^{\prime}$ and $\tau^{\prime\prime}=t_{3}-\tau_{s}^{\prime\prime}.$ This completes the derivation of the QP expression for $S_{I}$ (Eq. 24) starting from the SOS expression (Eq. 27).

$S_{II}$ can be calculated similarly. By combining Eqs. (28) and (31) the same type of cancellation of harmonic terms yields

[TABLE]

Eq. (33) is identical to Eq. (25) as illustrated in Fig. (10).

We finally turn to $\bm{S}_{III}$, (Eq. 29). Using again the Bethe Salpeter equation (31) and the fact that terms that only depend on $\mathcal{\hat{G}}_{0}$ must cancel (harmonic reference), we get

[TABLE]

The equivalence of QP (Eq. 26) and SOS (Eq. 29) expressions can be shown as follows: the two terms in Eq. (34) are labeled (SOSa) and (SOSb). The term (SOSb) can further be split into two terms (SOSb1) and (SOSb2), the first corresponding to $\tau^{\prime}<t_{2}$, the second to $\tau^{\prime}>t_{2}$ (Fig. 11). (SOSb1) is identical to (SOSa), but with opposite sign coming from $\hat{G}^{\dagger}(t_{3})$. Only the second term (SOSb2) remains, and it is equivalent to the (QP) diagram. We thus obtained Eq. (26) from Eq. (34).

V 2D correlation signals

2D signals are displayed as correlation plots obtained by the double Fourier transforms of the various signals.Mukamel (2000) We shall denote the frequencies conjugate to the pulse delay times $t_{1},t_{2}$ and $t_{3}$ by $\Omega_{1},$ $\Omega_{2}$ and $\Omega_{3}$. Starting with Eq. (11), and deleting some inessential factors, we obtain the induced polarization, which depends parametrically on the delay times $t_{1}$, $t_{2}$ and $t_{3}$:

[TABLE]

Specifying the three possible signals by a proper choice of $\lambda$ factors we obtain:

[TABLE]

The 2DCS for $\bm{P}_{I}$ and $\bm{P}_{II}$ is defined as

[TABLE]

For the SOS picture we use the expansions in eigenstates given by Eqs. (59), (60) and (61). The QP expressions for $\bm{P}_{I}^{\left(QP\right)}$, $\bm{P}_{II}^{\left(QP\right)}$ and $\bm{P}_{III}^{\left(QP\right)}$ are obtained along the lines presented in App. C. Dephasing is introduced phenomenologically by adding a decay rate $\gamma$ to the Green’s functions. We thus obtain $\bm{P}_{I}$

[TABLE]

[TABLE]

The Green’s function Fourier transform is defined as $G(\omega)=\int dt\exp(i\omega t)G(t)$ [and $G(t)=\int\frac{d\omega}{2\pi}\exp(-i\omega t)G(\omega)$]. We have

[TABLE]

Eq. (40) is obtained by transforming $\mathcal{G}_{0\,kljr}\left(t\right)$ to the single-exciton basis and performing the Fourier transform. We also define

[TABLE]

Similarly we obtain for $\bm{P}_{II}$:

[TABLE]

[TABLE]

The $\bm{P}_{III}$ 2DCS signal is defined as

[TABLE]

This yields:

[TABLE]

[TABLE]

Both $\bm{P}^{(SOS)}$ and $\bm{P}^{\left(QP\right)}$ depend on the single-exciton energies. However, the SOS expressions contain two-exciton eigenenergies ($\varepsilon_{f}$) explicitly, while the QP counterparts contain the scattering matrix $\Gamma$ instead. The equivalence of the two representations has been established in Sec. IV. Eqs. (37-45) constitute our final expressions for the various 2DCS signals. In this form they may be readily used in numerical simulations. The SOS expressions (Eqs. 37, 41 and 44) were recently used to survey the various possible resonances and cross-peaks in 2DCS of semiconductors.Yang et al. (2006)

VI Discussion

The quasiparticle representation is obtained using the Heisenberg equations for the exciton oscillator variables. These equations form an infinite hierarchy involving successively higher numbers of excitons.Abramavicius and Mukamel (2006a); Chernyak et al. (1998a) The hierarchy may be truncated, depending on the observable of interest. For instance, the absorption originates from single-exciton creation/annihilation. Only single-exciton variables should then be considered, and exciton-exciton interaction terms may be neglected. The nonlinear response depends on the exciton interactions, thus single- and double-exciton variables need to be treated explicitly. The two coupled NEE equations (21, 23) describe the third order response. These equations are exact in the absence of dephasing. When dephasing is included by adding linear coupling to a phonon bath, two additional variables $\left\langle\hat{B}^{\dagger}\hat{B}\right\rangle$ and $\left\langle\hat{B}^{\dagger}\hat{B}\hat{B}\right\rangle$ must be included in the NEE to describe the third order response. Chernyak et al. (1998a) Without dephasing these may be factorized as $\left\langle\hat{B}^{\dagger}\hat{B}\right\rangle=\left\langle\hat{B}^{\dagger}\right\rangle\left\langle\hat{B}\right\rangle$ and $\left\langle\hat{B}^{\dagger}\hat{B}\hat{B}\right\rangle=\left\langle\hat{B}^{\dagger}\right\rangle\left\langle\hat{B}\hat{B}\right\rangle$. We then recover the coherent limit considered in this article. For some techniques the present equations provide a good approximation even in the presence of dephasing. $\left\langle\hat{B}^{\dagger}\hat{B}\right\rangle$ describes incoherent exciton transport and is only relevant during $t_{2}$, while $\left\langle\hat{B}^{\dagger}\hat{B}\hat{B}\right\rangle$ is generated during $t_{3}$. It describes the optical coherence between one-exciton and two-exciton manifolds, which are represented by $\left\langle\hat{B}\right\rangle$ and $\left\langle\hat{B}\hat{B}\right\rangle$.

The quasiparticle approach avoids the explicit calculation of multiple exciton states: their influence is represented by the scattering matrix, which can be calculated provided the exciton interactions are known. We have shown how the quasiparticle expressions for the various third order techniques, ordinarily derived by solving equations of motion, can be obtained directly from the sum-over-states expressions by employing the Bethe-Salpether equation. These expressions explicitly contain the two-exciton Green’s functions and have many interfering terms with large cancellations,Spano and Mukamel (1989) which complicate their numerical implementation. In the QP picture, on the other hand, these interference effects are built-in, considerably simplifying the expressions for the nonlinear response.Spano and Mukamel (1991)

The interpretation of 2DCS signals using the SOS expressions is straightforward.Yang et al. (2006) In the $\mathbf{k}_{I}$ technique one-exciton coherences are observed during $t_{1}$, and the coherences between excitons and biexcitons are observed during $t_{3}$. Thus the 2DCS shows peaks along $\Omega_{1}$ and $\Omega_{3}$ corresponding to these resonances. $\mathbf{k}_{III}$ shows biexciton resonances along the $\Omega_{2}$ axis, this technique is known in NMR as double-quantum coherence. We have established the connection between the SOS and the QP pictures by using time-ordering on the Keldysh-Schwinger loop, which only maintains partial time ordering in real (physical) time.

The Hamiltonian $\hat{H}$ (Eq. 17) can describe several microscopic models other than the Wannier excitons considered here ($\hat{H}_{T}$). Vibrational excitations (soft-core bosons) and Frenkel excitons (hard-core bosons, Paulions) in molecules can be mapped into the same model. Zhang et al. (1999); Abramavicius and Mukamel (2006a); Kuhn et al. (1996) The equations of motion for these other systems are similar, but not identical, because of the different commutation relations (QP statistics). Eq. (19) provides a unified description for all of these systems, by specifying the proper commutation rules:Chernyak and Mukamel (1996); Chernyak et al. (1998a) for bosons $\mathcal{P}_{mnpq}=0$ and for Paulions $\mathcal{P}_{mnpq}=\delta_{mn}\delta_{mq}\delta_{np}$. These expressions for $\mathcal{P}$ may be substituted into our final expressions for the response functions, where they only affect the exciton scattering matrix, which in the frequency domain reads (App. D)

[TABLE]

where $V$ is given in Eq. (22), $\mathcal{G}_{0}$ is the free two–exciton Green’s function (App. B) and $\mathbb{I}$ is the tetradic identity matrix. The nonlinearity of the system depends on QP interactions as well as non-boson statistics; both enter through $\Gamma$. For noninteracting bosons, where $U\equiv\mathcal{P}\equiv 0$, $\Gamma$ vanishes and so does the nonlinear response. In Appendix F we present $\Gamma$ for bosons and Paulions.

We have used the symmetry $P_{mnpq}=P_{mnqp}$ in our derivation. Since the boson commutation relations are simpler than for Fermi or Pauli operators, a considerable effort has been devoted to mapping the original problem with complicated commutation relations into a boson picture.Ilinskaia and Ilinski (1996); Agranovich and Toshich (1967) The resulting boson Hamiltonian contains additional interactions which compensate for the statistics. For instance, the Frenkel exciton Hamiltonian for Paulions may be mapped into an anharmonic Hamiltonian of bosons with quartic couplings. BosonizationChernyak and Mukamel (1996) is very convenient for describing exciton scattering: the response functions derived for bosons can be applied for arbitrary operators, provided we modify the Hamiltonian and express it in terms of boson operators.

Acknowledgements.

We wish to thank Dr Igor V. Schweigert for valuable discussions. This research was supported by the National Science Foundation Grant no. CHE-0446555 and the National Institutes of Health 2RO1-GM59230-05.

Appendix A Exciton representation of the two-band Hamiltonian for fermions

By construction, the Hamiltonians $H$ (Eq. (17)) and $\hat{H}_{T}$ (Eq. (16)) are equivalent only in the physically relevant space of single and double excitations. This is sufficient to calculate the response to third order in the field $E\left(t\right)$. $\hat{H}$ may be constructed using the following rules:

• •

since the Hamiltonian (16) conserves the number of excitons, it should only contain products with equal number of $\hat{B}^{\dagger}$ and $\hat{B}$ operators (except for the $H_{I}$ term, which does change the number of excitons)

• •

a term $\hat{B}_{a_{1}}^{\dagger}\hat{B}_{a_{2}}^{\dagger}\ldots\hat{B}_{a_{p}}^{\dagger}\hat{B}_{b_{1}}\hat{B}_{b_{2}}\ldots\hat{B}_{b_{p}}$ gives zero when acting on states with less than $p$ excitations and only affects manifolds with $p$ excitations and higher.

The parameters of $\hat{H}$ can be obtained as follows. First we note that no constant term $k$ should be added to (17), since it would yield: $\left\langle 0\left|k\right|0\right\rangle\neq 0$, while $\left\langle 0\left|H_{T}\right|0\right\rangle=0$. The $\sum_{m1,n2,k2,l1}W_{m_{1}n_{2}l_{1}k_{2}}c_{m_{1}}^{\dagger}d_{n_{2}}^{\dagger}d_{k_{2}}c_{l_{1}}$ term of $H_{T}$ can be written directly as $\sum_{m_{1},m_{2},n_{1},n_{2}}W_{m_{1}m_{2}n_{1}n_{2}}c_{m_{1}}^{\dagger}d_{m_{2}}^{\dagger}d_{n2}c_{n1}=\sum_{m,n}W_{mn}\hat{B}_{m}^{\dagger}\hat{B}_{n}$. Also the term describing the interaction with light can be obtained directly. Using the second rule given above we immediately see that no terms higher than $\hat{B}_{a_{1}}^{\dagger}\hat{B}_{a_{2}}^{\dagger}\hat{B}_{b_{1}}\hat{B}_{b_{2}}$ are necessary in the sub-space defined by functions $\left|0\right\rangle$, $\hat{B}_{i}^{\dagger}\left|0\right\rangle$ and $\hat{B}_{i}^{\dagger}\hat{B}_{j}^{\dagger}\left|0\right\rangle$. We thus obtain the form given in (17). We next calculate, in this sub-space, matrix elements of $\hat{H}$, and compare to matrix elements of $\hat{H}_{T}$. In this way a one-to-one correspondence of the parameters of $\hat{H}$ and $\hat{H}_{T}$ can be established.

Additional terms must be included in $\hat{H}$ in order to describe higher order response functions. This can be done using the same rules.

Appendix B The Nonlinear Exciton Equations

The Heisenberg equation of motion (NEE) for the Hamiltonian (17) reads:

[TABLE]

Here we invoked RWA and used the notation of Eq. (1). Employing (20) we see that $\mathcal{P}_{nmpq}=\mathcal{P}_{nmqp}$, so the last two terms in Eq. (46) can be recast as: $2\bm{E}^{\bm{+}}\sum_{mpq}\mathcal{P}_{nmpq}\left\langle\hat{B}_{m}^{\dagger}\hat{B}_{q}\right\rangle\bm{\mu}_{p}^{\ast}$. We now make the following factorization:

[TABLE]

which is exact for pure states when dephasing is neglectedLeegwater and Mukamel (1992) and is a good approximation in the absence of incoherent exciton transport. Eq. (46) then yields the Eqs. (21) and (23), where

[TABLE]

We next expand the EOMs in orders of $E$. Using $B_{m}^{\left(1\right)}$ for $\left\langle\hat{B}_{m}\right\rangle^{\left(1\right)}$ we obtain:

[TABLE]

The Green’s function (tetradic matrix) for $Y^{\left(2\right)}$ is $\mathcal{G}\left(t\right)_{mnkl}=-i\theta\left(t\right)\left[\exp\left(-ih^{\left(Y\right)}t\right)\right]_{mnkl}$, thus

[TABLE]

We also define the zero-order tetradic Green’s function $\mathcal{G}_{0}$ for $Y^{\left(2\right)}$ for the case $V=0$, i.e. $\mathcal{G}_{0\,mnkl}\left(t\right)=-i\theta\left(t\right)\left[\exp\left(-i\bar{h}t\right)\right]_{mnkl}$, it will be used later. $B_{m}^{\left(3\right)}$ is given as

[TABLE]

This expression can be simplified using the symmetry $\mathcal{G}_{klfg}=\mathcal{G}_{klgf}$. At this point we introduce the tetradic exciton scattering matrix $\Gamma$ defined as:

[TABLE]

which in time domain can be written as (see App. D):

[TABLE]

where the tetradic identity matrix $\mathbb{I}_{fgjr}=\delta_{fj}\delta_{gr}.$ Since $\mathcal{G}\left(t-\tau\right)\sim\theta\left(t-\tau\right)$ is retarded, $\Gamma$ must be retarded as well, i.e., $\Gamma\left(t-\tau_{1}\right)\sim\theta\left(t-\tau_{1}\right)$. To proceed further we take advantage of the factorization:

[TABLE]

which can be easily shown in the single-exciton eigenbasis. After a rearrangement of terms we obtain:

[TABLE]

The 3rd order polarization is

[TABLE]

where we used $\Gamma_{n_{4}^{\prime}n_{3}^{\prime}n_{1}^{\prime}n_{2}^{\prime}}\left(t\right)=\Gamma_{n_{4}^{\prime}n_{3}^{\prime}n_{2}^{\prime}n_{1}^{\prime}}\left(t\right)$.

The above expression is finite only for $\tau_{2}>\tau_{1}.$ Hence there are 3 possible intervals for $\tau_{3},$ that define three contributions to the third order response function $\bm{S}^{\left(QP\right)}$

[TABLE]

This definition of $\bm{S}_{I}^{\left(QP\right)}$, $\bm{S}_{II}^{\left(QP\right)}$and $\bm{S}_{III}^{\left(QP\right)}$ is used in Appendix (C) to obtain the Eqs. (24-26).

Appendix C Response functions of quasiparticles

For calculating each of the contributions to Eq. (52) we need to switch to a different set of time-ordered variables. For $\bm{S}_{I}^{\left(QP\right)}$ we set: $\tau_{1}\rightarrow\tau_{2},\tau_{2}\rightarrow\tau_{3},\tau_{3}\rightarrow\tau_{1}$:

[TABLE]

where $\tau_{4}>\tau_{3}>\tau_{2}>\tau_{1}$.

Substituting $\tau_{s}^{\prime}=\tau_{4}-\tau^{\prime\prime}$ and $\tau_{s}^{\prime\prime}=\tau_{4}-\tau^{\prime}$ and exchanging the dummy indices $n_{1}\rightarrow n_{2},n_{2}\rightarrow n_{3},n_{3}\rightarrow n_{1}$ (same for primed indices) we obtain Eq. (24). Integration limits for $\tau_{s}^{\prime}$ have been limited by $G_{n_{4}n_{4}^{\prime}}\left(\tau_{s}^{\prime}\right)$ and $\Gamma\left(\tau_{s}^{\prime\prime}-\tau_{s}^{\prime}\right)$, while for $\tau_{s}^{\prime\prime}$ by $G_{n_{3}^{\prime}n_{3}}\left(\tau_{43}-\tau_{s}^{\prime\prime}\right)$. Eqs. (25) and (26) are obtained in a similar way.

Eq. (24) can be simplified considerably by performing the double time-integrations analytically. We first express the exciton Green’s function $G\left(\tau\right)$ and $\Gamma\left(\tau\right)$ in the one-exciton basis $\psi_{e}$, defined by:

[TABLE]

where $h_{mn}$ is given by Eq. (18). The energies $\varepsilon_{e}$ define the lowest optically-excited manifold of the system, i.e., single excitons. In this basis we express the time and frequency-domain one-exciton Green’s functions:

[TABLE]

where we introduce dephasing via $\varepsilon\rightarrow\varepsilon-i\gamma_{e}$. The tetradic exciton scattering matrix is given by

[TABLE]

and the transformed dipole matrix elements $\bm{\mu}_{eg}$ are $\bm{\mu}_{e}=\sum_{m}\bm{\mu}_{m}\psi_{em}$. Using these quantities we express the $\bm{S}_{I}^{\left(QP\right)}$ in the single-exciton basis:

[TABLE]

We next introduce $\Gamma\left(t\right)=\int\frac{d\omega}{2\pi}e^{-i\omega t}\Gamma\left(\omega\right)$. Since the response function depends only on the pulse delays and not on the absolute times, we denote $\bm{S}_{I}(\tau_{4},\tau_{3},\tau_{2},\tau_{1})=\bm{S}_{I}(t_{3},t_{2},t_{1})$, where $t_{3}=\tau_{43}$, $t_{2}=\tau_{32}$ and $t_{1}=\tau_{21}$. We perform a Fourier transform with respect to the first and last arguments. We thus obtain

[TABLE]

The $\omega$ integration can be performed by noting that

[TABLE]

which is obtained from Eq. (49) by noting that $\mathcal{G}\left(\omega\right)$ has poles only at two-exciton energies, and that $2\varepsilon-2i\gamma$ is a good approximation for two-exciton energy and dephasing rate. Hence, if we close the Cauchy integration path in the positive half-plane, there will be only a single pole at $\omega=\Omega_{3}+\varepsilon_{e_{1}}+i\gamma_{e_{1}}$ as seen from (39). This finally gives

[TABLE]

To account for carrier frequencies $\omega_{1}$, $\omega_{2}$ and $\omega_{3}$ appearing in the polarization (Eq. 35) at this stage, we can perform the substitution $\Omega_{1}\rightarrow\Omega_{1}+\omega_{1}$ and $\Omega_{3}\rightarrow-\Omega_{1}+\omega_{1}+\omega_{2}+\omega_{3}$. In this way we obtain Eq. (38). Eqs. (42) and (45) are derived similarly.

Appendix D The exciton scattering-matrix

In order to use equations (24-26) for calculating the quasiparticle response function, we should calculate the scattering matrix $\Gamma$. We first write $\mathcal{G}\left(\omega\right)$ and $\mathcal{G}_{0}\left(\omega\right)$ in an operator form

[TABLE]

where $\bar{h}$ is defined in (48). The Dyson equation then reads

[TABLE]

which can be recast in the form

[TABLE]

Using Eq. (49), we obtain:

[TABLE]

which results in the final expression for $\Gamma$

[TABLE]

The l.h.s. of Eq. (49) can be expressed as a convolution:

[TABLE]

The r.h.s. can be written as

[TABLE]

note that $\mathcal{P}$ is independent on $\tau^{\prime}$ or $\omega$. Since l.h.s.=r.h.s. for any $\omega$, we must have:

[TABLE]

Substituting $\tau^{\prime}\rightarrow t-\tau$ and $t_{1}\rightarrow t-\tau_{1}$ we obtain Eq. (50).

Appendix E SOS expressions for third order techniques.

Upon expansion in the eigenstates for the exciton level scheme shown in Fig. 3 we get

[TABLE]

Expanding Eqs. (13 - 15) in the eigenstates, we obtain the sum-over-states expressions for the third-order response functions:

[TABLE]

[TABLE]

[TABLE]

where $I_{e}(t)$, defined in Eq. (54), is the Green’s function in the single-exciton eigenstate basis. Eqs. (37,41,44) immediately follow by substituting Eqs. (59,60,61) in Eq. (36) and (43).

Appendix F Quasiparticle picture for soft-core and hard-core bosons

In this Appendix we apply our QP expressions to two other types of quasiparticles with different statistics. These two examples demonstrate the generality of our approach discussed briefly in Sec. VI.

We first consider the Hamiltonian of a system of coupled anharmonic oscillators (soft-core bosons):

[TABLE]

where $\hat{B}_{m}^{+}$ and $\hat{B}_{n}$ are boson creation and annihilation operators with commutation $[B_{m},B_{n}^{+}]=\delta_{mn}$ $h_{mm}$ is the fundamental transition energy of the $m$th oscillator, while $h_{mn}$ is the coupling between the $m$th and the $n$th oscillators. $U_{mnkl}$ is the anharmonic coupling. This Hamiltonian has been used to describe infrared nonlinear spectra of proteins.Abramavicius and Mukamel (2006a, 2005)

For this model the scattering matrix can be obtained from (58) by putting $\mathcal{P}=0.$ In the site representation it reads:

[TABLE]

here $\Gamma$ is a tetradic matrix, $V=2U$ and $\mathcal{G}_{0}\left(\omega\right)$ is defined in Eq. (30).

We next turn to electronic excitations in molecular aggregates or crystals with weakly interacting molecules. These are described using the Frenkel Exciton Hamiltonian. If the excited-state absorption frequency of each molecule is well separated from the ground state absorption, the excitations can be modelled as coupled two-level systems.Leegwater and Mukamel (1992); Juzeliunas and Knoester (2000) The Hamiltonian is

[TABLE]

The nonlinearities are now hidden in the statistics of exciton creation $(\hat{B}_{m}^{+})$ and annihilation $(\hat{B}_{n})$ operators. These are bosonic for different oscillators (units) and fermionic for the same oscillator. Their Pauli commutation relation is $[\hat{B}_{m},\hat{B}_{n}^{+}]=\delta_{mn}\left(1-2\hat{B}_{n}^{\dagger}\hat{B}_{n}\right)$. The commutation relation ensures that two excitations are not allowed to reside on the same site (hard-core bosons). The scattering matrix in this case is given by:

[TABLE]

and

[TABLE]

This form of the exciton scattering matrix was recently successfully applied to study molecular chirality induced signals in molecules.Abramavicius and Mukamel (2006b) It can be obtained from Eq. (58) in the limit $U=0$. All QP-statistics effects (Paulion commutation relations) are included in the relation $V_{nmpq}=-2\sum_{l}\mathcal{P}_{nmlp}h_{lq}$ (Eq. 22) and $\mathcal{P}_{nmlp}=\delta_{nm}\delta_{nl}\delta_{mp}$.