Sum-over-states vs quasiparticle pictures of coherent correlation
  spectroscopy of excitons in semiconductors; femtosecond analogues of
  multidimensional NMR
Sum-over-states vs quasiparticle pictures of coherent correlation
spectroscopy of excitons in semiconductors; femtosecond
analogues of multidimensional NMR
Shaul Mukamel, Rafal Oszwaldowski, Darius Abramavicius
Chemistry Department, University of California,
Irvine, CA 92697-2025, United States
(Dated: August 28, 2021)
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.
PACS numbers: 78.47.+p,71.35.-y
http://arxiv.org/abs/0704.0636v1
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,1,2,3 as well as vibrations in proteins.4,5 In semi-
conductors, nonlinear optical experiments reveal a wealth of interesting phenomena.6,7,8,9,10
For instance, such experiments provide information about many-exciton states such as biex-
citons, their interactions, relaxation and dissociation.11,12,13,14
The introduction of multidimensional techniques had revolutionized NMR in the
seventies15 and established it as a powerful tool for studying complex systems and iden-
tifying specific structural and dynamical correlations.16 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.5,17,18,19 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. 20,21,22. 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 inter-
suband excitations in semiconductors.23,24,25,26,27 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 work23 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 op-
tical response of excitonic systems. The first is based on the many-body eigenstates obtained
by exact diagonalization of the Hamiltonian.28 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 molecules29,30,31 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 ∼ N and ∼ N2 respectively. For Wannier excitons in semiconductors these
scalings are ∼ N2 and ∼ N4, 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.2,28,32,33,34,35 Calculations are per-
formed 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 quan-
tum well, solving the equations numerically to create a 2D map of a nonlinear response
function is computationally expensive,23 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 motion36,37 there exist other theoretical approaches to exciton correlation ef-
fects, such as memory kernel representation38,39 or Coupled Cluster Expansion for doped
semiconductors.40,41
In this paper we derive closed expressions for 2DCS of semiconductors by solving the
Nonlinear Exciton Equations (NEE)3,42 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. 3,34 to transform the Hamiltonian to a form typical
for interacting oscillators. The starting many-electron Hamiltonian can be written in an
ab-initio,43 tight-binding44 or a k · 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 NONLIN-
EAR RESPONSE
We consider a 4 wave-mixing experiment performed with three femtosecond laser pulses
(Fig. 1). The optical electric field is:
E (r, t) =
Ej (r, t) = E
+(r, t) + E−(r, t), (1)
E+(r, t) =
E+j (t− τj)e
−iωjteikjr, (2)
E−(r, t) =
E−j (t− τj)e
iωjte−ikjr. (3)
The j-th pulse is centered at τj, has an envelope Ej(t − τj), carrier frequency ωj , and
wavevector kj . E
+ (E−) denotes the positive (negative) frequency part of the field, and
E−j =
. The induced polarization in the system is recorded as a function of time-
delays between pulses.
Assuming the dipole interaction with the optical field ĤI = µ̂ · E(r, τ), where µ̂ is the
dipole operator, the third-order contribution to the system’s polarization can be written as
P (r, τ4) =
∫∫∫ ∞
dτ3dτ2dτ1S
(SOS)(τ4, τ3, τ2, τ1)E(r, τ3)E(r, τ2)E(r, τ1), (4)
where the response function S(SOS), which connects the induced polarization with the laser
field envelopes, is given by (throughout this paper we set ~ = 1):
(SOS)(τ4, τ3, τ2, τ1) = i
3 [θ(τ43)θ(τ32)θ(τ21) 〈µ̂(τ4)µ̂(τ3)µ̂(τ2)µ̂(τ1)〉 (5)
− θ(τ43)θ(τ42)θ(τ21) 〈µ̂(τ3)µ̂(τ4)µ̂(τ2)µ̂(τ1)〉
+ θ(τ42)θ(τ23)θ(τ41) 〈µ̂(τ3)µ̂(τ2)µ̂(τ4)µ̂(τ1)〉
−θ(τ41)θ(τ12)θ(τ23) 〈µ̂(τ3)µ̂(τ2)µ̂(τ1)µ̂(τ4)〉] .
We shall use double-sided Feynman diagrams to represent the time ordering of various
interactions.45 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 τi variables are ordered on the Keldysh-
Schwinger loop, but not necessarily in real (physical) time. τi in diagrams (a) and (d) are
also ordered in real time. This is not the case for diagrams (b) and (c): in (b) τ3 can come
either before or after τ1 and τ2, whereas in (c) τ1 can come either before or after τ3 and τ2.
If the eigenstates |a〉 and eigenvalues εa of the system are known, Eq. (5) may be expanded
in terms of the corresponding matrix elements:
〈µ(τ4)µ(τ3)µ(τ2)µ(τ1)〉 (6)
a1,a2,a3
µga3µa3a2µa2a1µa1ge
(εa3−εg)τ4+(εa2−εa3)τ3+(εa1−εa2)τ2+(εg−εa1)τ1
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 (g), single exciton (e), two-exciton (f) (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 µ̂ = µ̂+ + µ̂−, where µ̂+ is the positive frequency part
which induces upward g to e and e to f transitions, while its Hermitian conjugate µ̂− (the
negative frequency part) induces the opposite transitions. We thus write
εν>εν′
µνν′ |ν〉 〈ν
εν<εν′
µνν′ |ν〉 〈ν
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
HI (t) = −µ̂
+(r, τ)− µ̂−E−(r, τ)
Each correlation function in Eq. (5) will split into 24 = 16 terms upon substituting
µ̂ = µ̂+ + µ̂−. Assuming that the system is initially in the ground state, only two of these
contributions are non-zero
〈µ̂µ̂µ̂µ̂〉 =
. (7)
Substitution of Eq. (7) into Eq. (5) gives
(SOS)(τ4, τ3, τ2, τ1) = i
θ(τ43)θ(τ32)θ(τ21)
−(τ4)µ̂
−(τ3)µ̂
+(τ2)µ̂
+(τ1)
(a1) (8)
+ θ(τ43)θ(τ32)θ(τ21)
−(τ4)µ̂
+(τ3)µ̂
−(τ2)µ̂
+(τ1)
− θ(τ43)θ(τ42)θ(τ21)
−(τ3)µ̂
−(τ4)µ̂
+(τ2)µ̂
+(τ1)
+θ(τ42)θ(τ23)θ(τ41)
−(τ3)µ̂
+(τ2)µ̂
−(τ4)µ̂
+(τ1)
+ c.c.
The four terms represented by the diagrams in Fig. 4 were obtained by taking µ̂(τ4) =
µ̂−(τ4) for the last interaction, µ̂(τ4) = µ̂
+(τ4) gives the complex conjugates. Hereafter
left/right direction of the arrows corresponds to µ̂−/µ̂+ in Eq.(8). Note that time-reversal
symmetry implies 〈µ̂−(τ4)µ̂
−(τ3)µ̂
−(τ2)µ̂
−(τ1)〉
= 〈µ̂+(τ1)µ̂
+(τ2)µ̂
+(τ3)µ̂
+(τ4)〉. If the
pulse envelopes are much shorter than their delays, the system is forced to interact se-
quentially first with pulse k1, then k2 and finally k3. This means that in the integral of
Eq. (4) one must replace E(r, τj) with one of the Ej, 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
τ3 relative to the τ1 < τ2 < τ4 sequence, the diagram can be separated into three fully time
ordered terms: τ3 < τ1, τ1 < τ3 < τ2 or τ2 < τ3 < τ4. Formally we do that by separating the
product of step functions as follows:
θ(τ43)θ(τ42)θ(τ21) = θ(τ42)θ(τ21)θ(τ13) + θ(τ42)θ(τ23)θ(τ31) + θ(τ43)θ(τ32)θ(τ21).
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)
θ(τ23)θ(τ42)θ(τ41) = θ(τ42)θ(τ21)θ(τ13) + θ(τ42)θ(τ23)θ(τ31) + θ(τ41)θ(τ12)θ(τ23)
and the diagram is split into (c2), (c3), (c1). S(3) can now be recast in the fully time-ordered
(SOS)(τ4, τ3, τ2, τ1) = i
θ(τ43)θ(τ32)θ(τ21)
−(τ4)µ̂
−(τ3)µ̂
+(τ2)µ̂
+(τ1)
(a1) (9)
+ θ(τ43)θ(τ32)θ(τ21)
−(τ4)µ̂
+(τ3)µ̂
−(τ2)µ̂
+(τ1)
− θ(τ43)θ(τ32)θ(τ21)
−(τ3)µ̂
−(τ4)µ̂
+(τ2)µ̂
+(τ1)
− θ(τ42)θ(τ23)θ(τ31)
−(τ3)µ̂
−(τ4)µ̂
+(τ2)µ̂
+(τ1)
− θ(τ42)θ(τ21)θ(τ13)
−(τ3)µ̂
−(τ4)µ̂
+(τ2)µ̂
+(τ1)
+ θ(τ41)θ(τ12)θ(τ23)
−(τ3)µ̂
+(τ2)µ̂
−(τ4)µ̂
+(τ1)
+ θ(τ42)θ(τ23)θ(τ31)
−(τ3)µ̂
+(τ2)µ̂
−(τ4)µ̂
+(τ1)
+θ(τ42)θ(τ21)θ(τ13)
−(τ3)µ̂
+(τ2)µ̂
−(τ4)µ̂
+(τ1)
+ c.c.
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 integra-
tion variables in Eq. (4) from τ4, τ3, τ2, τ1 that label the actual interaction times with the
fields, to the three delays t3, t2, t1 between successive interactions. Note that the correla-
tion functions are invariant to time translation 〈µ̂(τ4 − τ)µ̂(τ3 − τ)µ̂(τ2 − τ)µ̂(τ1 − τ)〉 =
〈µ̂(τ4)µ̂(τ3)µ̂(τ2)µ̂(τ1)〉. Eq. (4) thus assumes the form
P (r, τ4) =
∫∫∫ +∞
dt3dt2dt1S
(3)(t3, t2, t1)E3(r, τ4−t3)E2(r, τ4−t2−t3)E1(r, τ4−t1−t2−t3),
where t1 = τ2 − τ1, t2 = τ3 − τ2, t3 = τ4 − τ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
P (r, τ4) = S
(SOS)(t3, t2, t1)E
3 (τ4 − t3)E
2 (τ4 − t3 − t2)E
1 (τ4 − t3 − t2 − t1)× (11)
ei(λ1k1+λ2k2+λ3k3)re−i(λ1ω1+λ2ω2+λ3ω3)τ4ei(λ1ω1+λ2ω2+λ3ω3)t3ei(λ1ω1+λ2ω2)t2eiλ1ω1t1
The polarization is created along 8 possible directions ks = λ1k1+λ2k2+λ3k3 with λi = ±1
P (r, τ4) =
P (ks, ωs) e
iksr−iωsτ4 + c.c. (12)
where
P (ks, ωs) = S
(SOS)
s (t3, t2, t1) E
k1 + k2 + k3 vanishes for the assumed dipole selection rules in our model. Since
P (−ks,−ωs) = P
∗ (ks, ωs), we are left with three independent combinations kI ≡
−k1 + k2 + k3, kII ≡ +k1 − k2 + k3, kIII ≡ +k1 + k2 − k3:
(SOS)(t3, t2, t1) = S
(SOS)
I (t3, t2, t1) + S
(SOS)
II (t3, t2, t1) + S
(SOS)
III (t3, t2, t1).
We can classify the diagrams in Fig. 5 according to the directions of the arrows: arrow
pointing to the right (left) represents +k (−k), arrows are read from the bottom up on
either side. We obtain for kI (Fig. 6)
(SOS)
I = i
3θ(t1)θ(t2)θ(t3)
−(0)µ̂+(t1 + t2)µ̂
−(t1 + t2 + t3)µ̂
+(t1)
(c3) (13)
−(0)µ̂+(t1)µ̂
−(t1 + t2 + t3)µ̂
+(t1 + t2)
−(0)µ̂−(t1 + t2 + t3)µ̂
+(t1 + t2)µ̂
+(t1)
. (b3)
For the kII technique we similarly have (Fig. 7):
(SOS)
II = i
3θ(t1)θ(t2)θ(t3)
−(t1 + t2 + t3)µ̂
+(t1 + t2)µ̂
−(t1)µ̂
(a2) (14)
−(t1)µ̂
+(t1 + t2)µ̂
−(t1 + t2 + t3)µ̂
−(t1)µ̂
−(t1 + t2 + t3)µ̂
+(t1 + t2)µ̂
Finally kIII is given by (Fig. 8):
(SOS)
III = i
3θ(t1)θ(t2)θ(t3)
−(t1 + t2 + t3)µ̂
−(t1 + t2)µ̂
+(t1)µ̂
(a1) (15)
−(t1 + t2)µ̂
−(t1 + t2 + t3)µ̂
+(t1)µ̂
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 expres-
sions for these signals. These will then be connected with the current SOS expressions in
Section IV.
III. QUASIPARTICLE EXPRESSIONS FOR WANNIER EXCITONS IN SEMI-
CONDUCTORS
Interband transitions in semiconductors may be described by the two-band many-electron
Hamiltonian:42,46
ĤT = Ĥ0 + ĤC + ĤI , (16)
with the single-particle part
Ĥ0 =
m1,n1
t(1)m1,n1c
cn1 +
m2,n2
t(2)m2,n2d
dn2 ,
where c† create electrons and d† create holes. The Coulomb interaction is:
ĤC =
m1,n1,k1,l1
m1n1k1l1
c†m1c
ck1cl1 +
m2,n2,k2,l2
m2n2k2l2
d†m2d
dk2dl2
m1,n2,k2,l1
Wm1n2l1k2c
d†n2dk2cl1 ,
while
ĤI = −
m1,m2
+ (t)µ∗m1m2c
d†m2 +E
− (t)µm1m2dm2cm1
is the dipole interaction with light, and the optical electric field E will be treated as a scalar
for simplicity. 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. Ĥ0 would
be diagonal in the basis of the system’s single-particle eigenstates, i.e., t
m1,n1 = ε
m1δm1n1.
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.2,3 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:34
B̂m ≡ dm2cm1 , B̂
m ≡ c
d†m2 ,
where we have employed shorthand notation for pairs of indices: m ≡ (m1, m2). Using these
operators we construct an effective Hamiltonian Ĥ (see App. A):
hmnB̂
mB̂n +
UmnklB̂
nB̂kB̂l −
+ (t)µ∗mB̂
− (t)µmB̂m
. (17)
The Hamiltonians Ĥ and ĤT are equivalent in the single and double excitations subspace,
which is relevant for the response to third order in E.28 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 Ĥ are given by:
hmn = t
m1,n1
δm2n2 + t
m2,n2
δm1n1 −Wm1m2n1n2 , (18)
Umnkl = −
m1,k1
δm2k2δn1l1δn2l2 + t
m2,k2
δm1k1δn1l1δn2l2+
n1,l1
δm1k1δm2k2δn2l2 + t
n2,l2
δm1k1δm2k2δn1l1
m1n1k1l1
δm2k2δn2l2 + V
m2n2k2l2
δm1k1δn1l1
The commutation relations for the B̂ operators can be obtained using the elementary fermion
anticommutators:
c†m1 , ck1
= δm1,k1. Within the subspace of |0〉 and B̂
i |0〉 states (i.e., the
ground state and single excitations), we get34
B̂m, B̂
= δmn − 2
PmnpqB̂
pB̂q , (19)
where δmn = δm1n1δm2n2 and
Pmnpq =
δm1q1δp1n1δm2p2δn2q2 +
δm2q2δp2n2δm1p1δn1q1. (20)
Eqs. (19) and (20) are obtained in a similar way to (17) and (18). Terms with additional
i B̂j pairs (e.g. B̂
†B̂†B̂B̂) are neglected in (19), because they would introduce corrections
higher than O (E3) to the nonlinear response. Note the symmetry Pmnpq = Pmnqp.
Using Eqs. (17) and (19) we obtain the nonlinear exciton equations (see Appendix B) for
single-exciton variables
:2,17,35,47
− µ∗mE
+ (t) +
Vmnkl
B̂kB̂l
+ 2E+ (t)
Pmnpq
where V is given by
Vnmpq = 2Unmpq − 2
Pnmlphlq − 2
PnmklUklpq. (22)
Here Ymn ≡
B̂mB̂n
are two-exciton variables. The Heisenberg equations give:
mn,kl
Ykl −E
+ (t)
+ 2E+ (t)
Pmnkl
Calculating the optical response by numerical integration of these equations23,48 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 (t) and exciton scattering matrix Γ (t). The
scattering matrix depends on quasiparticle statistics through the P matrix (Eqs. B4, B5) as
well as on exciton-exciton coupling. This results in closed quasiparticle expressions for the
3rd order contributions SI , SII and SIII to the response function (for details see Appendix
C and Ref. 49)
(QP )
I (τ4, τ3, τ2, τ1) = (24)
− 2θ (τ43) θ (τ32) θ (τ21)
n4...n1
dτ ′′s
τ ′′s
dτ ′s
(τ ′′s − τ
s)Gn4n′4 (τ
s)Gn′3n3 (τ43 − τ
s )Gn′2n2 (τ42 − τ
(τ41 − τ
where τ43 = τ4 − τ3, etc. and Gmn (t) = −iθ (t) [exp (−iht)]mn.
The response functions for the other phase-matching directions can be derived along the
same lines. We get
(QP )
II (τ4, τ3, τ2, τ1) = (25)
− 2θ (τ43) θ (τ32) θ (τ21)
n4...n1
dτ ′′s
τ ′′s
dτ ′s
...n′
(τ ′′s − τ
s)Gn4n′4(τ
s)Gn′3n3(τ43 − τ
(τ42 − τ
s)Gn′1n1(τ41 − τ
(QP )
III (τ4, τ3, τ2, τ1) = (26)
− 2θ (τ43) θ (τ32) θ (τ21)
n4...n1
µn4µn3µ
dτ ′′s
τ ′′s
dτ ′s
...n′
(τ ′′s − τ
s)Gn4n′4(τ
(τ43 − τ
s)Gn′2n2(τ42 − τ
s )Gn′1n1(τ41 − τ
Just as in the SOS case, time translation symmetry implies that these response functions
only depend on the three pulse delays t3, t2, t1. 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 t1 > 0, t2 > 0,
t3 > 0):
(SOS)
−Ĝ†(t1 + t2 + t3)µ̂
−Ĝ(t3)µ̂
+Ĝ(t2)µ̂
(b3) (27)
−Ĝ†(t1 + t2)µ̂
†(t3)µ̂
−Ĝ(t2 + t3)µ̂
−Ĝ† (t1) µ̂
†(t2 + t3)µ̂
−Ĝ(t3)µ̂
, (c1)
(SOS)
II =−
−Ĝ(t3)µ̂
Ĝ(t2)µ̂
−Ĝ(t1)µ̂
(a2) (28)
−Ĝ†(t2 + t3)µ̂
−Ĝ(t3)µ̂
+Ĝ(t1 + t2)µ̂
−Ĝ†(t2)µ̂
†(t3)µ̂
−Ĝ(t1 + t2 + t3)µ̂
, (c2)
(SOS)
III =−
−Ĝ(t3)µ̂
−Ĝ(t2)µ̂
+Ĝ(t1)µ̂
(a1) (29)
−Ĝ†(t3)µ̂
−Ĝ(t2 + t3)µ̂
+Ĝ(t1)µ̂
, (b1)
Here Ĝ(t) ≡ −iθ(t) exp(−iĤt) and Ĝ†(t) ≡ +iθ(t) exp(iĤt) represent the retarded and
the advanced Green’s function respectively; Ĝ, Ĝ and Ĝ describe the evolution within the
ground-state, single-exciton and double-exciton blocks of the Hamiltonian (Eq. 17) respec-
tively. We also set the ground state energy ε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 identically2,28,32,50, 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 equation33
Ĝ (ω) = Ĝ0 (ω) + Ĝ0 (ω) Γ (ω) Ĝ0 (ω) , (30)
or in the time domain:
Ĝ(τ) = Ĝ0(τ) +
∫ τ ′
dτ ′′Ĝ0 (τ − τ
′) Γ (τ ′ − τ ′′) Ĝ0(τ
′′). (31)
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
Ĝ0(τ)n1n2n3n4 = iĜn1n3(τ)Ĝn2n4(τ).
The exciton scattering matrix Γ is defined by Eq. (31).
Let us start with the SI 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 µ±,
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. Ĝ0, represents harmonic evolution in the two-exciton manifold. Thus the
first term of Eq. (27) with Ĝ replaced by 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 SI :
(SOS)
I = −θ (t3) θ (t2) θ (t1)
∫ τ ′
dτ ′′× (32)
−Ĝ†(t1 + t2 + t3)µ̂
−Ĝ0(t3 − τ
′)Γ(τ ′ − τ ′′)Ĝ0(τ
′′)µ̂+Ĝ(t2)µ̂
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 Γ is represented by dashed regions. Note
that G (t)G† (t) = θ (t) exp(−iht) exp(iht) = θ (t). The QP diagram in Fig. (9) is obtained
from the SOS one by changing the integration variables τ ′ = t3 − τ
s and τ
′′ = t3 − τ
This completes the derivation of the QP expression for SI (Eq. 24) starting from the SOS
expression (Eq. 27).
SII can be calculated similarly. By combining Eqs. (28) and (31) the same type of
cancellation of harmonic terms yields
(SOS)
II = −θ (t3) θ (t2) θ (t1)
∫ τ ′
dτ ′′× (33)
−Ĝ† (t2 + t3) µ̂
−Ĝ0 (t3 − τ
′) Γ (τ ′ − τ ′′) Ĝ0 (τ
′′) µ̂+Ĝ(t1 + t2)µ̂
Eq. (33) is identical to Eq. (25) as illustrated in Fig. (10).
We finally turn to SIII , (Eq. 29). Using again the Bethe Salpeter equation (31) and the
fact that terms that only depend on Ĝ0 must cancel (harmonic reference), we get
(SOS)
III = −θ (t3) θ (t2) θ (t1)× (34)
∫ τ ′
dτ ′′
−Ĝ(t3)µ̂
−Ĝ0(t2 − τ
′)Γ(τ ′ − τ ′′)Ĝ0(τ
′′)µ̂+Ĝ(t1)µ̂
∫ t2+t3
∫ τ ′
dτ ′′
−Ĝ†(t3)µ̂
−Ĝ0(t2 + t3 − τ
′)Γ(τ ′ − τ ′′)Ĝ0(τ
′′)µ̂+Ĝ(t1)µ̂
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 τ ′ < t2, the second
to τ ′ > t2 (Fig. 11). (SOSb1) is identical to (SOSa), but with opposite sign coming from
Ĝ†(t3). 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.17 We shall denote the frequencies conjugate to the pulse delay times
t1, t2 and t3 by Ω1, Ω2 and Ω3. Starting with Eq. (11), and deleting some inessential factors,
we obtain the induced polarization, which depends parametrically on the delay times t1, t2
and t3:
Ps(t3, t2, t1) = Ss(t3, t2, t1)e
i(λ1ω1+λ2ω2+λ3ω3)t3ei(λ1ω1+λ2ω2)t2eiλ1ω1t1 . (35)
Specifying the three possible signals by a proper choice of λ factors we obtain:
PI(t3, t2, t1) = SI(t3, t2, t1)e
i(−ω1+ω2+ω3)t3ei(−ω1+ω2)t2e−iω1t1 ,
PII(t3, t2, t1) = SII(t3, t2, t1)e
i(ω1−ω2+ω3)t3ei(ω1−ω2)t2eiω1t1 ,
PIII(t3, t2, t1) = SIII(t3, t2, t1)e
i(ω1+ω2−ω3)t3ei(ω1+ω2)t2eiω1t1 .
The 2DCS for PI and PII is defined as
Pα(Ω3, t2,Ω1) ≡
dt1Pα(t3, t2, t1) exp {iΩ3t3 + iΩ1t1} , α = I, II (36)
For the SOS picture we use the expansions in eigenstates given by Eqs. (E1), (E2) and (E3).
The QP expressions for P
(QP )
I , P
(QP )
II and P
(QP )
III are obtained along the lines presented
in App. C. Dephasing is introduced phenomenologically by adding a decay rate γ to the
Green’s functions. We thus obtain PI
(SOS)
I (Ω3, t2,Ω1) = (37)
e2,e1
µge1µ
µge2I
(−Ω1 + ω1)I
g (t2)Ig(t2)Ie1(Ω3 − ω1 + ω2 + ω3)
e2,e1
µge1µ
µge2I
(−Ω1 + ω1)I
(t2)Ie1(t2)Ie1(Ω3 − ω1 + ω2 + ω3)
e2,e1f
µe2fµ
µge2I
(−Ω1 + ω1)I
(t2)Ie1(t2)Ffe2(Ω3 − ω1 + ω2 + ω3),
(QP )
I (Ω3, t2,Ω1) = (38)
e1,e2,e3,e4
(t2)Ie2(t2)I
(−Ω1 − ω1)Ie4(Ω3 − ω1 + ω2 + ω3)
× Γe4e1e3e2(Ω3 − ω1 + ω2 + ω3 + εe1 + iγe1)G0 e3e2(Ω3 − ω1 + ω2 + ω3 + εe1 + iγe1).
The Green’s function Fourier transform is defined as G(ω) =
dt exp(iωt)G(t) [and G(t) =
exp(−iωt)G(ω)]. We have
Ie(ω) ≡
Ĝ (ω)
= (ω − εe + iγe)
−1, (39)
G0 e2e1(ω) ≡
Ĝ0(ω)
ω − εe2 − εe1 + i (γe2 + γe1)
. (40)
Eq. (40) is obtained by transforming G0 kljr (t) to the single-exciton basis and performing the
Fourier transform. We also define
Fab(t) ≡ −iθ(t) exp (i (εb − εa) t− (γa + γb) t) ,
Fab(ω) = (ω − εa + εb + iγa + iγb)
Similarly we obtain for PII :
(SOS)
II (Ω3, t2,Ω1) = (41)
e2,e1
µge2µge1µ
Ie1(Ω1 + ω1)I
g (t2)Ig(t2)Ie2(Ω3 + ω1 − ω2 + ω3)
e2,e1
µge1µ
µge2µ
Ie1(Ω1 + ω1)I
(t2)Ie1(t2)Ie1(Ω3 + ω1 − ω2 + ω3)
e2,e1,f
fµge2µ
Ie1(Ω1 + ω1)I
(t2)Ie1(t2)Ffe2(Ω3 + ω1 − ω2 + ω3),
(QP )
II (Ω3, t2,Ω1) = (42)
e4..e1
I∗e2(t2)Ie1(t2)Ie1(Ω1 + ω1)Ie4(Ω3 + ω1 − ω2 + ω3)
× Γe4e2e3e1(Ω3 + ω1 − ω2 + ω3 + εe2 + iγe2)G0 e3e1(Ω3 + ω1 − ω2 + ω3 + εe2 + iγe2).
The PIII 2DCS signal is defined as
PIII(Ω3,Ω2, t1) ≡
dt2PIII(t3, t2, t1) exp {iΩ3t3 + iΩ2t2} . (43)
This yields:
(SOS)
III (Ω3,Ω2, t1) = (44)
e2,e1,f
µge1µe1fµ
Ie2(t1)If(Ω2 + ω1 + ω2)Ie1(Ω3 + ω1 + ω2 − ω3)
e2,e1,f
µge1µe1fµ
Ie2(t1)If(Ω2 + ω1 + ω2)Ffe1(Ω3 + ω1 + ω2 − ω3),
(QP )
III (Ω3,Ω2, t1) = (45)
e4...e1
µe4µe3µ
Ie1(t1)Ie4(Ω3 + ω1 + ω2 − ω3)I
(Ω2 − Ω3 + ω3)×
[Γe4e3e2e1(Ω2 + ω1 + ω2)G0 e2e1(Ω2 + ω1 + ω2)
−Γe4e3e2e1(Ω3 + ω1 + ω2 − ω3 + εe3 + iγe3)G0 e2e1(Ω3 + ω1 + ω2 − ω3 + εe3 + iγe3)] .
Both P (SOS) and P (QP ) depend on the single-exciton energies. However, the SOS expres-
sions contain two-exciton eigenenergies (εf) explicitly, while the QP counterparts contain the
scattering matrix Γ 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.23
VI. DISCUSSION
The quasiparticle representation is obtained using the Heisenberg equations for the ex-
citon oscillator variables. These equations form an infinite hierarchy involving successively
higher numbers of excitons.3,51 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 ab-
sence of dephasing. When dephasing is included by adding linear coupling to a phonon bath,
two additional variables
B̂†B̂
B̂†B̂B̂
must be included in the NEE to describe the
third order response.3 Without dephasing these may be factorized as
B̂†B̂
B̂†B̂B̂
. We then recover the coherent limit considered in this arti-
cle. For some techniques the present equations provide a good approximation even in the
presence of dephasing.
B̂†B̂
describes incoherent exciton transport and is only relevant
during t2, while
B̂†B̂B̂
is generated during t3. It describes the optical coherence between
one-exciton and two-exciton manifolds, which are represented by
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 exci-
ton 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 inter-
fering terms with large cancellations,32 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.35
The interpretation of 2DCS signals using the SOS expressions is straightforward.23 In
the kI technique one-exciton coherences are observed during t1, and the coherences between
excitons and biexcitons are observed during t3. Thus the 2DCS shows peaks along Ω1 and
Ω3 corresponding to these resonances. kIII shows biexciton resonances along the Ω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 Ĥ (Eq. 17) can describe several microscopic models other than the Wan-
nier excitons considered here (ĤT ). Vibrational excitations (soft-core bosons) and Frenkel
excitons (hard-core bosons, Paulions) in molecules can be mapped into the same model.19,50,51
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:3,34 for bosons Pmnpq = 0
and for Paulions Pmnpq = δmnδmqδnp. These expressions for 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)
Γ(ω) = (I− V G0(ω))
V G0(ω) (I− P)G
0 (ω)− PG
0 (ω),
where V is given in Eq. (22), G0 is the free two–exciton Green’s function (App. B) and I
is the tetradic identity matrix. The nonlinearity of the system depends on QP interactions
as well as non-boson statistics; both enter through Γ. For noninteracting bosons, where
U ≡ P ≡ 0, Γ vanishes and so does the nonlinear response. In Appendix F we present Γ for
bosons and Paulions.
We have used the symmetry Pmnpq = Pmnqp in our derivation. Since the boson commuta-
tion relations are simpler than for Fermi or Pauli operators, a considerable effort has been de-
voted to mapping the original problem with complicated commutation relations into a boson
picture.52,53 The resulting boson Hamiltonian contains additional interactions which com-
pensate for the statistics. For instance, the Frenkel exciton Hamiltonian for Paulions may be
mapped into an anharmonic Hamiltonian of bosons with quartic couplings. Bosonization34 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.
Acknowledgments
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 HAMIL-
TONIAN FOR FERMIONS
By construction, the Hamiltonians H (Eq. (17)) and 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 (t). 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 B̂† and B̂ operators (except for the HI term, which
does change the number of excitons)
• a term B̂†a1B̂
. . . B̂†apB̂b1B̂b2 . . . B̂bp gives zero when acting on states with less than p
excitations and only affects manifolds with p excitations and higher.
The parameters of Ĥ can be obtained as follows. First we note that no con-
stant term k should be added to (17), since it would yield: 〈0 |k| 0〉 6= 0, while
〈0 |HT | 0〉 = 0. The
m1,n2,k2,l1Wm1n2l1k2c
d†n2dk2cl1 term of HT can be written directly
m1,m2,n1,n2
Wm1m2n1n2c
d†m2dn2cn1 =
WmnB̂
mB̂n. Also the term describing the
interaction with light can be obtained directly. Using the second rule given above we imme-
diately see that no terms higher than B̂†a1B̂
B̂b1B̂b2 are necessary in the sub-space defined
by functions |0〉, B̂
i |0〉 and B̂
j |0〉. We thus obtain the form given in (17). We next
calculate, in this sub-space, matrix elements of Ĥ , and compare to matrix elements of ĤT .
In this way a one-to-one correspondence of the parameters of Ĥ and ĤT can be established.
Additional terms must be included in Ĥ in order to describe higher order response func-
tions. 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:
− µ∗nE
Vnmpq
B̂†mB̂pB̂q
Pnmpq
B̂†mB̂p
B̂†mB̂q
Here we invoked RWA and used the notation of Eq. (1). Employing (20) we see that Pnmpq =
Pnmqp, so the last two terms in Eq. (B1) can be recast as: 2E
Pnmpq
B̂†mB̂q
We now make the following factorization:
B̂†mB̂p
B̂†mB̂pB̂q
B̂pB̂q
, (B2)
which is exact for pure states when dephasing is neglected33 and is a good approximation
in the absence of incoherent exciton transport. Eq. (B1) then yields the Eqs. (21) and (23),
where
mn,kl
= δmkhnl + δnlhmk + Vmnkl ≡ h̄ + V. (B3)
We next expand the EOMs in orders of E. Using B
m for
we obtain:
n − µ
+ (t) ,
mn,kl
−E+ (t)
B(1)n µ
+ 2E+ (t)
PmnklB
VmnklB
+ 2E+ (t)
PmnpqB
The Green’s function (tetradic matrix) for Y (2) is G (t)
= −iθ (t)
−ih(Y )t
Y (2)mn (t) =
Gmnkl (t− τ)E
+ (τ)
(τ)µ∗k +B
(τ)µ∗l
PklpqB
p (τ)µ
Gmnkl (t− τ)×
Gla (τ − τ
′)µ∗k +Gka (τ − τ
′)µ∗l − 2
PklpqGpa (τ − τ
′)µ∗q
+ (τ)E+ (τ ′) .
We also define the zero-order tetradic Green’s function G0 for Y
(2) for the case V = 0, i.e.
G0mnkl (t) = −iθ (t)
−ih̄t
, it will be used later. B
m is given as
B(3)m (t
Gmn (t
′ − t)×
VnpklB
p (t) Y
kl (t) + 2E
+ (t)
PnkpqB
k (t)B
q (t)µ
This expression can be simplified using the symmetry Gklfg = Gklgf . At this point we
introduce the tetradic exciton scattering matrix Γ defined as:
Γ (ω)G0 (ω) = V G (ω) (I−P)−P, (B4)
which in time domain can be written as (see App. D):
V G (t− τ) (I− P) = Pδ (t− τ) +
dτ1Γ (t− τ1)G0 (τ1 − τ) , (B5)
where the tetradic identity matrix Ifgjr = δfjδgr. Since G (t− τ) ∼ θ (t− τ) is retarded, Γ
must be retarded as well, i.e., Γ (t− τ1) ∼ θ (t− τ1). To proceed further we take advantage
of the factorization:
G0 kljr (t) = iGkj (t)Glr (t) , (B6)
which can be easily shown in the single-exciton eigenbasis. After a rearrangement of terms
we obtain:
n4 (τ4) = −2
∫ ∫ ∫ ∫ ∫
dτ ′dτ2dτ1dτ
′′dτ3
n1,n2,n3
Gn4n′4 (τ4 − τ
′′) Γn′
(τ ′′ − τ ′)G∗n′
(τ ′′ − τ3)×
n2 (τ
′ − τ2)Gn′
n1 (τ
′ − τ1) θ (τ2 − τ1)µ
+ (τ2)E
+ (τ1)E
− (τ3) .
The 3rd order polarization is
P (3) (τ4) =
+ µ∗n4B
∫ ∫ ∫
dτ3dτ2dτ1
n4,n3,n2,n1
µn4µn3µ
θ (τ2 − τ1)
dτ ′′
(τ ′′ − τ ′)×
Gn4n′4 (τ4 − τ
′′)G∗n′
(τ ′′ − τ3)Gn′
n2 (τ
′ − τ2)Gn′
n1 (τ
′ − τ1)×
− (τ3)E
+ (τ2)E
+ (τ1) + complex conjugate,
where we used Γn′
(t) = Γn′
The above expression is finite only for τ2 > τ1. Hence there are 3 possible intervals for τ3,
that define three contributions to the third order response function S(QP )
P (τ4) =
dτ1θ (τ2 − τ1)× (B7)
(QP )
I dτ3 +
(QP )
II dτ3 +
(QP )
III dτ3
− (τ3)E
+ (τ2)E
+ (τ1) + c.c.
This definition of S
(QP )
I , S
(QP )
II and S
(QP )
III 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. (B7) we need to switch to a different set
of time-ordered variables. For S
(QP )
I we set: τ1 → τ2, τ2 → τ3, τ3 → τ1:
(QP )
I (τ4, τ3, τ2, τ1) = −2
n4,n3,n2,n1
µn4µn3µ
dτ ′′
(τ ′′ − τ ′)Gn4n′4 (τ4 − τ
′′)G∗n′
(τ ′′ − τ1)Gn′
n2 (τ
′ − τ3)Gn′
n1 (τ
′ − τ2) ,
where τ4 > τ3 > τ2 > τ1.
Substituting τ ′s = τ4 − τ
′′ and τ ′′s = τ4 − τ
′ and exchanging the dummy indices n1 →
n2, n2 → n3, n3 → n1 (same for primed indices) we obtain Eq. (24). Integration limits for τ
have been limited by Gn4n′4 (τ
s) and Γ (τ
s − τ
s), while for τ
s by Gn′3n3 (τ43 − τ
s ). Eqs. (25)
and (26) are obtained in a similar way.
Eq. (24) can be simplified considerably by performing the double time-integrations ana-
lytically. We first express the exciton Green’s function G (τ) and Γ (τ) in the one-exciton
basis ψe, defined by:
hmnψen = εeψem, (C1)
where hmn is given by Eq. (18). The energies ε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:
Gmn (τ) =
ψemIe (τ)ψ
en ⇒ Ie (τ) = −iθ (τ) exp (−iεet) , (C2)
where we introduce dephasing via ε → ε − iγe. The tetradic exciton scattering matrix is
given by
Γm4m3m2m1 (τ) =
e1...e4
ψe4m4ψe3m3Γe4e3e2e1 (τ)ψ
ψ∗e1m1 .
and the transformed dipole matrix elements µeg are µe =
µmψem. Using these quantities
we express the S
(QP )
I in the single-exciton basis:
(QP )
I (τ4, τ3, τ2, τ1) = −2θ (τ43) θ (τ32) θ (τ21)
e1...e4
µe3 (C3)
∫ τ43
dτ ′′s
∫ τ ′′s
dτ ′sΓe4e3e2e1(τ
s − τ
× Ie4(τ
s)Ie2(τ43 − τ
s )Ie1(τ42 − τ
(τ41 − τ
We next introduce Γ (t) =
e−iωtΓ (ω). Since the response function depends only on the
pulse delays and not on the absolute times, we denote SI(τ4, τ3, τ2, τ1) = SI(t3, t2, t1), where
t3 = τ43, t2 = τ32 and t1 = τ21. We perform a Fourier transform with respect to the first and
last arguments. We thus obtain
(QP )
I (Ω3, t2,Ω1) = −2iθ (t2)
e1...e4
(−Ω1) I
(t2) Ie
(t2) Ie
(Ω3)× (C4)
(ω)G0 e
(ω) I∗e1 (ω − Ω3) ,
The ω integration can be performed by noting that
Γ(ω)G0(ω) ∼
ω − 2ε+ 2iγ
which is obtained from Eq. (B4) by noting that G (ω) has poles only at two-exciton energies,
and that 2ε−2iγ 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 ω = Ω3 + εe1 + iγe1 as seen from (39). This finally gives
(QP )
I (Ω3, t2,Ω1) = −2
e4...e1
(t2)Ie2(t2)I
(−Ω1)Ie4(Ω3) (C5)
× Γe4e1e3e2(Ω3 + εe1 + iγ1)G0 e3e2(Ω3 + εe1 + iγ1),
To account for carrier frequencies ω1, ω2 and ω3 appearing in the polarization (Eq. 35) at
this stage, we can perform the substitution Ω1 → Ω1 + ω1 and Ω3 → −Ω1 + ω1 + ω2 + ω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 Γ. We first write G (ω) and G0 (ω) in an operator
G0 (ω) =
ω − h̄
G (ω) =
ω − h(Y )
ω − h̄− V
where h̄ is defined in (B3). The Dyson equation then reads
G = G0 + G0V G (ω) = G0 + G0V G0 + G0V G0V G0 + . . . ,
which can be recast in the form
V G = (1− V G0)
V G0.
Using Eq. (B4), we obtain:
ΓG0 = (I− V G0)
V G0 (I− P)− P,
which results in the final expression for Γ
Γ = (I− V G0)
V G0 (I− P)G
0 −PG
0 (D1)
The l.h.s. of Eq. (B4) can be expressed as a convolution:
dt1Γ (t1)G0 (τ
′ − t1) exp (iωτ
The r.h.s. can be written as
dτ ′G (τ ′) (I−P) exp (iωτ ′)−
dτ ′Pδ (τ ′) exp (iωτ ′) ,
note that P is independent on τ ′ or ω. Since l.h.s.=r.h.s. for any ω, we must have:
dt1Γ (t1)G0 (τ
′ − t1) = V G (τ
′) (I−P)− Pδ (τ ′) .
Substituting τ ′ → t− τ and t1 → t− τ1 we obtain Eq. (B5).
APPENDIX E: SOS EXPRESSIONS FOR THIRD ORDER TECHNIQUES.
Upon expansion in the eigenstates for the exciton level scheme shown in Fig. 3 we get
Expanding Eqs. (13 - 15) in the eigenstates, we obtain the sum-over-states expressions
for the third-order response functions:
(SOS)
I (t3, t2, t1) = i
3θ(t3)θ(t2)θ(t1)
µge′µe′gµgeµegI
e′ (t1) Ie (t3) (E1)
+ i3θ(t3)θ(t2)θ(t1)
µge′µe′gµgeµegI
e′ (t2 + t1) Ie (t2 + t3)
− i3θ(t3)θ(t2)θ(t1)
µge′µe′fµfeµegI
e′ (t1 + t2 + t3) If (t3) Ie (t2) ,
(SOS)
II (t3, t2, t1) = i
3θ(t3)θ(t2)θ(t1)
µge′µe′gµgeµegIe′ (t3) Ie (t1) (E2)
+ i3θ(t3)θ(t2)θ(t1)
µge′µe′gµgeµegI
e′ (t2) Ie (t1 + t2 + t3)
− i3θ(t3)θ(t2)θ(t1)
µge′µe′fµfeµegI
e′ (t2 + t3) If (t3) Ie (t1 + t2) ,
(SOS)
III (t3, t2, t1) = i
3θ(t3)θ(t2)θ(t1)
µge′µe′fµfeµegIe′ (t3) If (t2) Ie (t1) (E3)
− i3θ(t3)θ(t2)θ(t1)
µge′µe′fµfeµegI
e′ (t3) If (t2 + t3) Ie (t1) ,
where Ie(t), defined in Eq. (C2), is the Green’s function in the single-exciton eigenstate
basis. Eqs. (37,41,44) immediately follow by substituting Eqs. (E1,E2,E3) 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):
hmnB̂
mB̂n +
UmnklB̂
nB̂kB̂l,
where B̂+m and B̂n are boson creation and annihilation operators with commutation
[Bm, B
n ] = δmn hmm is the fundamental transition energy of the mth oscillator, while hmn
is the coupling between the mth and the nth oscillators. Umnkl is the anharmonic coupling.
This Hamiltonian has been used to describe infrared nonlinear spectra of proteins.49,51
For this model the scattering matrix can be obtained from (D1) by putting P = 0. In the
site representation it reads:
Γ = (I− V G0)
here Γ is a tetradic matrix, V = 2U and G0 (ω) 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.33,54 The Hamilto-
nian is
hmnB̂
mB̂n.
The nonlinearities are now hidden in the statistics of exciton creation (B̂+m) and annihilation
(B̂n) operators. These are bosonic for different oscillators (units) and fermionic for the
same oscillator. Their Pauli commutation relation is [B̂m, B̂
n ] = δmn
1− 2B̂†nB̂n
. 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:
Γmnkl = δmnδklΓ̄mn ,
Γ̄ = −Ḡ (ω)
Ḡmn (ω) = δmm1δnn1G0mm1nn1 (ω) .
This form of the exciton scattering matrix was recently successfully applied to study molec-
ular chirality induced signals in molecules.55 It can be obtained from Eq. (D1) in the limit
U = 0. All QP-statistics effects (Paulion commutation relations) are included in the relation
Vnmpq = −2
Pnmlphlq (Eq. 22) and Pnmlp = δnmδnlδmp.
1 A. Davydov, A Theory of Molecular Excitons (McGraw-Hill: New York, 1962).
2 S. Mukamel, in Molecular Nonlinear Optics, Materials, Physics, and Devices, J. Zyss ed. (Else-
vier, Amsterdam) (1993).
3 V. Chernyak, W. M. Zhang, and S. Mukamel, J. Chem. Phys. 109, 9587 (1998).
4 W. Zhuang, D. Abramavicius, and S. Mukamel, Proc. Nat. Acad. Sci. USA 102, 7443 (2005).
5 S. Mukamel and D. Abramavicius, Chem. Rev. 104, 2073 (2004).
6 D. S. Chemla and J. Shah, Nature 411, 549 (2001).
7 R. A. Kaindl, M. Wurm, K. Reimann, M. Woerner, T. Elsaesser, C. Miesner, K. Brunner, and
G. Abstreiter, Phys. Rev. Lett. 86, 1122 (2001).
8 T. Shih, K. Reimann, M. Woerner, T. Elsaesser, I. Waldmüller, A. Knorr, R. Hey, and K. H.
Ploog, Phys. Rev. B 72, 195338 (2005).
9 J. Shah, Ultrafast Spectroscopy of Semiconductors and Semiconductor Nanostructures, Springer
Series in Solid-State Sciences (Springer, 1999), 2nd ed.
10 T. Meier, P. Thomas, and S. W. Koch, Coherent Semiconductor Optics: From Basic Concepts
to Nanostructure Applications (Springer, 2006), 1st ed.
11 D. S. Chemla, C. R. Acad. Sci. Paris, t.2, Serie IV, Solides, fluides: proprietes electroniques et
optiques/ Solids, fluids: electronic and optical properties pp. 1427–1438 (2001).
12 S. Adachi, K. Hazu, T. Sota, S. F. Chichibu, G. Cantwell, D. B. Eason, D. C. Reynolds, and
C. W. Litton, Semicond. Sci. Technol. 19, S276 (2004).
13 J. H. Rice, J. W. Robinson, J. H. Na, K. H. Lee, R. A. Taylor, D. P. Williams, E. P. O’Reilly,
A. D. Andreev, Y. Arakawa, and S. Yasin, Nanotechnology 16, 1477 (2005).
14 J. Danckwerts, K. J. Ahn, J. Förstner, and A. Knorr, Phys. Rev. B 73, 165318 (2006).
15 W. P. Aue, E. Bartholdi, and R. R. Ernst, Journal Of Chemical Physics 64, 2229 (1976).
16 R. R. Ernst, G. Bodenhausen, and A. Wokaun, Principles of Nuclear magnetic Resonance in
One and Two Dimensions (Oxford Science Publications, 1989).
17 S. Mukamel, Annu. Rev. Phys. Chem. 51, 691 (2000).
18 Y. Tanimura and S. Mukamel, J. Chem. Phys. 99, 9496 (1993).
19 W. M. Zhang, V. Chernyak, and S. Mukamel, J. Chem. Phys. 110, 5011 (1999).
20 C. Scheurer and S. Mukamel, Bull. Chem. Soc. Jpn. 75, 989999 (2002).
21 C. Scheurer and S. Mukamel, J. Chem. Phys. 115, 4889 (2001).
22 C. Scheurer and S. Mukamel, J. Chem. Phys. 116, 6803 (2002).
23 L. Yang, I. V. Schweigert, S. T. Cundiff, and S. Mukamel, accepted by Phys. Rev. B,
cond-mat/0701424 (2006).
24 I. Kuznetsova, P. Thomas, T. Meier, T. Zhang, X. Li, R. P. Mirin, and S. T. Cundiff, accepted
for publication in Solid State Communications, cond-mat/0702107 (2007).
25 M. Erementchouk, M. N. Leuenberger, and L. J. Sham, cond-mat/0611196 (2006).
26 X. Li, T. Zhang, C. N. Borca, and S. T. Cundiff, Physical Review Letters 96, 057406 (2006).
27 C. N. Borca, T. Zhang, X. Li, and S. T. Cundiff, Chemical Physics Letters 416, 311315 (2005).
28 S. Mukamel, Principles of Nonlinear Optical Spectrscopy (Oxford University Press, New York,
1995).
29 T. Brixner, J. Stenger, H. M. Vaswani, M. Cho, R. E. Blankenship, and G. R. Fleming, Nature
434, 625 (2005).
30 C. Fang, J. Wang, Y. S. Kim, A. K. Charnley, W. Barber-Armstrong, A. B. Smith III, S. M.
Decatur, and R. M. Hochstrasser, J. Phys. Chem. B 108, 10415 (2004).
31 D. Abramavicius, W. Zhuang, and S. Mukamel, J. Phys. Chem. B 108, 18034 (2004).
32 F. C. Spano and S. Mukamel, Phys. Rev. A 40, 5783 (1989).
33 J. A. Leegwater and S. Mukamel, Phys. Rev. A 46, 452 (1992).
34 V. Chernyak and S. Mukamel, J. Opt. Soc. Am. B 13, 1302 (1996).
35 F. C. Spano and S. Mukamel, Phys. Rev. Lett. 66, 1197 (1991).
36 S. W. Koch, C. Sieh, T. Meier, F. Jahnke, A. Knorr, P. Brick, M. Hubner, C. Ell, J. Prineas,
G. Khitrova, et al., Journal of Luminescence 83-84, 1 (1999).
37 S. Weiser, T. Meier, J. Möbius, A. Euteneuer, E. J. Mayer, W. Stolz, M. Hofmann, W. W.
Rühle, P. Thomas, and S. W. Koch, Phys. Rev. B 61, 13088 (2000).
38 T. Östreich, K. Schönhammer, and L. J. Sham, Phys. Rev. Lett. 74, 4698 (1995).
39 V. M. Axt, S. R. Bolton, U. Neukirch, L. J. Sham, and D. S. Chemla, Phys. Rev. B 63, 115303
(2001).
40 N. Primozich, T. V. Shahbazyan, I. E. Perakis, and D. S. Chemla, Phys. Rev. B 61, 2041 (2000).
41 T. V. Shahbazyan, N. Primozich, I. E. Perakis, and D. S. Chemla, Phys. Rev. Lett. 84, 2006
http://arxiv.org/abs/cond-mat/0701424
http://arxiv.org/abs/cond-mat/0702107
(2000).
42 V. M. Axt and S. Mukamel, Rev. Mod. Phys. 70 (1), 145 (1998).
43 R. Oszwa ldowski, M. Reichelt, T. Meier, S. W. Koch, and M. Rohlfing, Phys. Rev. B 71, 235324
(2005).
44 C. Sieh, T. Meier, F. Jahnke, A. Knorr, S. W. Koch, P. Brick, M. Hübner, C. Ell, J. Prineas,
G. Khitrova, et al., Phys. Rev. Lett. 82, 3112 (1999).
45 S. Mukamel, Phys. Rev. B 72, 235110 (2005).
46 H. Haug and S. W. Koch, Quantum Theory of the Optical and Electronic Properties of Semi-
conductors (World Scientific, Singapore, 2004), 4th ed.
47 V. Chernyak, S. Yokojima, T. Meier, and S. Mukamel, Phys. Rev. B 58, 4496 (1998).
48 C. Sieh, T. Meier, A. Knorr, F. Jahnke, P. Thomas, and S. Koch, Eur. Phys. J. B 11, 407
(1999).
49 D. Abramavicius and S. Mukamel, Chem. Phys. 318, 50 (2005).
50 O. Kuhn, V. Chernyak, and S. Mukamel, J. Chem. Phys. 105, 8586 (1996).
51 D. Abramavicius and S. Mukamel, J. Chem. Phys. 124, 034113 (2006).
52 A. V. Ilinskaia and K. N. Ilinski, J. Phys. A: Math. Gen. 29, L23 (1996).
53 V. M. Agranovich and B. S. Toshich, Zh. Eksp. Teor. Fiz. 53, 149 (1967), [Sov. Phys. JETP
26, 104 (1968)].
54 G. Juzeliunas and J. Knoester, J. Chem. Phys. 112, 2325 (2000).
55 D. Abramavicius and S. Mukamel, Journal of Physics B: Atomic, Molecular and Optical Physics
39, 5051 (2006).
FIG. 1: The sequence of light pulses in a time-domain Four Wave Mixing Experiment: the pulses
are centered at times τ1, τ2, τ3, while the delays are t1, t2 and t3. (The latter are sometimes
denoted as τ , T and t.) The signal is generated in the kS direction.
(a) (b) (c) (d)
FIG. 2: Diagrams representing the four partially time-ordered terms contributing to the third-order
polarization (Eq. 5).
FIG. 3: Energy levels of the exciton model.
(a1) (a2) (b) (c)
FIG. 4: Loop diagrams representing the partially time-ordered terms contributing to the third-
order polarization [Eq. (8)] within the rotating wave approximation. Arrows pointing to the right
represent µ+ and arrows pointing to the left µ−. The diagrams are obtained by adding arrows to
the interactions in Fig. (2).
(b2)(b1) (b3)
(c2) (c3)(c1)
FIG. 5: Diagrams representing the fully time-ordered terms contributing to the third-order polar-
ization within the rotating wave approximation (Eq. (9)). τj repesent the interaction times with
the various fields. Arrows pointing to the right (left) represent µ+ (µ−). Time variables in loop
diagrams (b) and (c) on the left are ordered in the loop. The other open diagrams are fully ordered
in physical time.
'e'egg
FIG. 6: Feynman diagrams for the kI technique (Eq. 13).
g gg g g g
g g g g
FIG. 7: Feynman diagrams for the kII technique (Eq. 14).
(a1) 
g gg g
g g 'e'e
FIG. 8: Feynman diagrams for the kIII technique (Eq. 15).
xx Ws”
(SOS)
FIG. 9: Loop diagrams showing the equivalence of the SI expressions in the QP (Eq. 24) and SOS
(Eq. 32) pictures. Dotted, single and double lines show ground, single exciton and double exciton
states’ evolution respectively. Dashed region represents scattering matrix.
xx Ws”
(SOS)
FIG. 10: Loop diagrams showing the order of time variables in the QP (Eq. 25) and SOS (Eq. 33)
expressions for SII .
FIG. 11: Loop diagrams showing the order of time variables in the QP (Eq. 26) and SOS (Eq. 34)
expressions for SIII . (SOSb1) cancels with the (SOSa) term in (Eq. 34) (not shown). The remaining
diagram (SOSb2) is identical to the (QP) diagram with a simple change of time variables.
	Introduction
	Sum-over-states expressions for the time-ordered nonlinear response 
	Quasiparticle expressions for Wannier excitons in semiconductors 
	Connecting the sum-over-states and the quasiparticle pictures
	2D correlation signals
	Discussion
	Acknowledgments
	Exciton representation of the two-band Hamiltonian for fermions 
	The Nonlinear Exciton Equations
	Response functions of quasiparticles
	The exciton scattering-matrix
	SOS expressions for third order techniques.
	Quasiparticle picture for soft-core and hard-core bosons
	References
