Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Phys Rev A. Author manuscript; available in PMC 2010 October 27.
Published in final edited form as:
PMCID: PMC2964889

Nonlinear optical spectroscopy of single, few, and many molecules; nonequilibrium Green’s function QED approach


Nonlinear optical signals from an assembly of N noninteracting particles consist of an incoherent and a coherent component, whose magnitudes scale ~ N and ~ N(N − 1), respectively. A unified microscopic description of both types of signals is developed using a quantum electrodynamical (QED) treatment of the optical fields. Closed nonequilibrium Green’s function expressions are derived that incorporate both stimulated and spontaneous processes. General (n + 1)-wave mixing experiments are discussed as an example of spontaneously generated signals. When performed on a single particle, such signals cannot be expressed in terms of the nth order polarization, as predicted by the semiclassical theory. Stimulated processes are shown to be purely incoherent in nature. Within the QED framework, heterodyne-detected wave mixing signals are simply viewed as incoherent stimulated emission, whereas homodyne signals are generated by coherent spontaneous emission.


Nonlinear optical processes in bulk materials are traditionally classified as either coherent or incoherent [1]–[4]. This classification refers to the relation of the macroscopic signal intensity of the sample to the nonlinear response of the individual molecules. Processes where the detected signal is obtained by simply adding up the contributions of the individual particles in the sample, are termed incoherent. Such signals hence scale as the number of particles N in the interaction volume. With coherent processes, on the other hand, the contributions of the constituent particles add up on the level of amplitudes rather than intensities, and the detected signal is proportional to N2.

Coherent optical signals are usually calculated semiclassically by adopting a two step procedure [2] – [4]. First a (linear or nonlinear) polarization P(n) induced in the system by interactions with classical external fields is calculated microscopically. In the second step, P(n) is used as a source in the macroscopic Maxwell equations to generate the signal. Frequency-domain signals are given by the absolute square of an amplitude related to the susceptibilities χ(n) which, in turn, are proportional to N. This gives rise to the ~ N2 dependence of coherent signals [5],[6]. Examples are Rayleigh scattering (n = 1), sum frequency generation (n = 2), four-wave mixing (n = 3). In contrast, signals such as Raman, fluorescence, two photon fluorescence etc. are incoherent and may not be calculated by the above semiclassical approach. Instead, they require a quantum description of the field to account for spontaneous emission.

Treating both types of processes by a unified approach is of fundamental interest [1], particularly when they coexist. For example, two photon fluorescence and second harmonic generation have been observed in the same sample [7]. There is a considerable current effort to perform nonlinear optical measurements on small samples and even single nano particles or molecules [8]–[11]. Recent experiments include 4-wave mixing spectroscopy of single semiconductor quantum dots [12] or nonlinear optical measurements on single electrons in F centers[13]. CARS microscopy [14] is being extended to small interaction volumes. The traditional macroscopic formulation does not apply for small samples and a quantum electrodynamical (QED) formulation is called for.

A QED treatment is essential for predicting observables related to the quantum nature of the field. A powerful quantum master equation approach has been successfully used in quantum optics to describe fluctuations of the laser field, photon statistics, entanglement and squeezing [15]–[25]. In this paper we consider simple observables that are usually calculated by the semiclassical approach. However a quantum description of the fields allows a more precise definition of the signal, provides new insights and shows the limitations of the semiclassical approach. Our final results are expressed in terms of multipoint correlation functions of the material system. These can be evaluated for more complicated models than the Bloch-type equations ordinarily used in quantum optics. Details of molecular complexity, coupling to a solvent with arbitrary time scale and excitonic effects may be readily described by these expressions [3].

The present approach is based on the many-body non-equilibrium Green’s functions (NEGF) technique [26] – [32]. It provides a single microscopic definition of the signal that includes both spontaneous and stimulated emission and generates both incoherent and coherent signals. We shall use a superoperator in Liouville space representation that is commonly applied to study electrical conduction in quantum junctions [33]–[35]. In Sec. II we give a general definition of the signal (Eq. (11)) and express it in terms of Green’s functions of the system and the field. The formalism will be demonstrated by first considering the signal generated by a single molecule. Applications are made to two types of incoherent experiments, spontaneous light emission (SLE) (Sec. III) and to pump-probe spectroscopy (Sec. IV). A purely diagrammatic derivation of the final expressions for the SLE signal using the technique of Keldysh Schwinger loops [30] is developed. We present rules that directly translate the diagrams to partially-time-ordered expressions for the signal, both in the time-and in the frequency-domain. This representation, has been shown to be most useful for the study of frequency-domain experiments within the semiclassical approximation [36]. We extend it to quantum fields, and show that it yields more compact expressions than the fully-time-ordered double-sided Feynman diagrams.

In Sec. V, we turn to the response of an assembly of noninteracting particles. The expressions of Sec. II are generalized by simply replacing the dipole operator of a single molecule by that of an assembly of particles. The signal is separated into an incoherent and a coherent part. The latter, which is described by semiclassical optical susceptibilities, is dominant in macroscopic samples. However for small N both have comparable magnitudes and coexist. The coherent signal does not exist in the limit of a single molecule. This indicates a breakdown of the semiclassical theory which does predict such a signal. Finally, in Sec. VI we show that within QED heterodyne detection emerges as an incoherent stimulated emission process in the detection mode. This is in contrast with semiclassical treatment where we first generate a coherent macroscopic signal, interfere it with a local oscillator and then obtain the signal. All of these steps are avoided by the microscopic QED definition given here. Sec. VII provides a summary of our results.


We first consider a single molecule driven by an optical field E(t). This is all we need for calculating incoherent signals. The extension to coherent signals from molecular ensembles will be done in Sec. V.

The total Hamiltonian is given by


where Ĥ0 represents the free molecule, and


is the radiation field Hamiltonian. In Eq. (2) a^s(a^s) denotes the destruction (creation) operator for the sth field mode and a^s(a^s) satisfy the boson commutation relation [a^s,a^s]=δss. For clarity, we shall omit the unit vector describing the field polarization of a particular mode, and treat the field as a scalar. The generalization to include the full tensorial expressions is straightforward.

We shall treat the emitted photon modes with frequency ωs quantum mechanically. An external file that holds a picture, illustration, etc.
Object name is nihms212364ig1.jpg(An external file that holds a picture, illustration, etc.
Object name is nihms212364ig2.jpg) are the positive (negative) frequency components of the quantized electric field




and Ω is the quantization volume.

In the rotating wave approximation (RWA) the molecule-field interaction is given by


where we have partitioned the dipole operator [mu] as [mu] = V + V. Here V(V) are the creation (annihilation) operators for excitations. The time-dependence in Eq. (5) is in the interaction picture with respect to ĤF. Consequently, the time dependence of the total (field+molecule) density operator [rho with circumflex]T(t) is governed by the Hamiltonian


We assume that the field is initially in a coherent state,




is the normalization such that left angle bracketΨFFright angle bracket = 1. In Eq. (7) αs is the eigenvalue of âs, âsFright angle bracket = αsF right angle bracket, and |0right angle bracket is the vacuum state. The expectation value of the field is then




is the field amplitude at space-point r.

We shall define the time and frequency resolved signal S(t) as the rate of change of the photon occupation, i.e.




denotes the photon number operator. This definition applies to both spontaneous and stimulated process, as will be demonstrated later. In Eq. (12) the sum extends over all modes of interest in the experiment under consideration.

We shall adopt of the following convention for the various ensemble averages, one of which has been used in Eq. (11). In Eq. (11), (…)D,T stands for a trace with respect to the density operator of the molecule where Ĥint includes all field modes (“total” system). According to the definition of the signal, Eqs. (11) and (12), the field modes can be partitioned into detected (or signal) modes and the incoming modes. For later use, we introduce the average (…)D, which denotes a trace with respect to the density operator of the molecule calculated in the interaction picture where Ĥint only includes the incoming modes. Averages with respect to a noninteracting system, where Ĥint=0, will be denoted by left angle bracketright angle bracket. This corresponds to the limit t → −∞ (adiabatic switching).

The expectation value in Eq. (11) may be conveniently evaluated by switching to the Heisenberg picture, using the basic identity


and the canonical commutation relations for a^s(a^s). In Eq. (13) we denote operators in the Heisenberg picture by a subscript “H”.

Evaluating the commutator in Eq. (13) and transforming back to the original interaction picure, we finally obtain for the signal, Eq. (11),


Taking into account the time dependence of the total (system + field) density operator, which is determined by the Hamiltonian of Eq. (6), the trace on the rhs of Eq. (14) can be computed perturbatively in the field E(r, t) by switching to the interaction picture with respect to the Hamiltonian Ĥ0. We shall compute this trace using superoperators in Liouville space. To this end we associate with each operator  in Hilbert space, a left (L) and a right (R) operation in Liouville space, by defining


where X denotes an arbitrary Hilbert space operator. Furthermore, we introduce the following linear combinations of L/R operations, which will be referred to as +/− operations


A key tool in the following manipulations is the time-ordering operator in Liouville space, An external file that holds a picture, illustration, etc.
Object name is nihms212364ig3.jpg; when acting on a product of superoperators, it reorders them such that their time arguments increase from right to left. In the interaction picture Eq. (14) can be expressed as


where An external file that holds a picture, illustration, etc.
Object name is nihms212364ig4.jpg is the superoperator corresponding to Ĥint, Eq. (5).

Note that Ĥint contains two interaction terms, one with the incoming field modes H^int and the other with the signal modes H^int. As before, the explicit time dependence represents the interaction picture


where X = L, R and An external file that holds a picture, illustration, etc.
Object name is nihms212364ig5.jpg (An external file that holds a picture, illustration, etc.
Object name is nihms212364ig6.jpg) is the superoperator corresponding to ĤF(Ĥ0).

To lowest order in c H^int Eq. (17) gives


The superoperator expression corresponding to Eq. (5) gives


Substituting in Eq. (19) and factorizing the correlation functions into products of field and system parts, we get


where we have introduced the field and system nonequilibrium Green’s functions, dXY and DXY,


In Eq. (21) we assume that the system is in its ground state for t → −∞, hence the following system correlation functions vanish


From the basic definitions of the Green’s functions we have


Substituting Eq. (24) in (21) gives


Using Eqs. (7) and (22), the Green’s functions for the field are given by


The signal, Eq. (25), then becomes


In Eq. (28), the term with δss represents spontaneous emission, whereas the terms proportional to asas correspond to stimulated processes. Note also that the Green’s functions DXY contain the density matrix of the driven system which involves all interactions with the incoming field.

Equation (28) can also be written as


or equivalently


In Eq. (30) we made use of the symmetry relation


to arrive at an expression, where interactions at time t always occur from the left (i.e. on the ket). This choice will be adopted in the diagrammatic representation, which will be introduced in Sec. III.

Equations (29) or equivalently (30) constitute our basic NEGF expression for incoherent optical signals. The first term in the brackets describes the creation of excitations in the system by absorbing photons from the field and the second term represents the reverse process. Both processes are stimulated. The signal is thus given by the net photon flux. The last term represents spontaneous emission. In the coming two sections we shall apply these expressions to compute different signals.


To describe spontaneous light emission (SLE), which is an incoherent process, we assume that the scattered field is initially in its vacuum state and is generated by spontaneous emission. Thus An external file that holds a picture, illustration, etc.
Object name is nihms212364ig7.jpg(r, t) = 0 and Eq. (30), reduces to


where it is furthermore assumed that the scattered field has only one mode ωs thus dropping the sum over s (signal modes). By integrating over modes we can replace 1/Ω by the density of modes.

Expanding the Green’s function DLR to second order in the incoming field yields


In order to go any further, we need to expand the four-point correlation function of Eq. (33) in L/R operations. A systematic, diagrammatic technique for selecting the terms in Eq. (33) that contribute within the RWA, will be introduced next.

We start with Eq. (33), where the correlation function originates from expanding the Green’s function DLR(t, τ) in Eq. (32) to second order in the incoming field. Hence in what follows, we associate the time variables τ1 and τ2 with the incoming field, whereas the variables t and τ correspond to the signal field. We also note that thanks to the time ordering operator, Eq. (33) is symmetric with respect to interchanging the dummy variables τ1 and τ2.

In order to derive an explicit expression for this correlation function in terms of superoperators, we first note that the interactions at times t and τ both represent photon emission, one on the ket at time t, the other on the bra at time τ. Since the system is assumed to be initially in its ground state, both emission events have to be preceded by two absorptions on either side. Although this implies a certain order of interactions on both the ket and the bra individually, their relative ordering in physical time remains a priori undetermined. This partial time ordering can be expressed diagrammatically by arranging the interactions on a loop. We adopt the following general rules for constructing and reading the diagrams:

  1. Time runs along the loop clockwise from bottom left to bottom right.
  2. The left strand of the loop represents the ket, the right corresponds to the bra.
  3. Each interaction with a field mode is represented by a wavy line on either the right (R-operators) or the left (L-operators).
  4. The field is indicated by dressing the wavy lines with arrows, where an arrow pointing to the right represents the field annihilation operator An external file that holds a picture, illustration, etc.
Object name is nihms212364ig1.jpg(r, t), which involves the term ei(ks·rωst) (see Eq. (4)). Conversely, an arrow pointing to the left corresponds to the field creation operator An external file that holds a picture, illustration, etc.
Object name is nihms212364ig2.jpg(r, t), being associated with ei(ks·rωst). This is made explicit by adding the wavevectors ±ks to the arrows.
  5. Within the RWA (Eq. (5)), each interaction with An external file that holds a picture, illustration, etc.
Object name is nihms212364ig1.jpg(r, t) is accompanied by applying the operator V, which leads to excitation of the state represented by ket and deexcitation of the state represented by the bra, respectively. Arrows pointing “inwards” (i.e. pointing to the right on the ket and to the left on the bra) consequently cause absorption of a photon by exciting the system, whereas arrows pointing “outwards” (i.e. pointing to the left on the bra and to the right on the ket) represent deexciting the system by photon emission.
  6. The interaction at the observation time t, is fixed and is always the last. As a convention, it is chosen to occur from the left. This can always be achieved by a reflection of all interactions through the center line between the ket and the bra, which corresponds to taking the complex conjugate of the original correlation function.
  7. Interactions within each strand are time-ordered, but interactions on different strands are not. Each loop can be further decomposed into several fully-time-ordered diagrams (double sided Feynman diagrams). These can be generated from the loop by simply shifting the arrows along each strand, thus changing their position relative to the interactions on the other strand. Each of these relative positions then gives rise to a particular fully time-ordered diagram.
  8. The overall sign of the correlation function is given by (−1)NR, where NR stands for the number of interactions from the right, at times τ1 and τ2.

We note that loop diagrams drawn according to the rules presented above, lead to double sided Feynman diagrams that follow the standard conventions employed within the semiclassical theory of nonlinear optics [3]. Using these rules, the SLE is represented by the single loop diagram displayed in Fig. 1. We denote the incoming field as ±k1 and the signal field as ±k2. The loop translates Eq. (33) into the following expression for the signal,

FIG. 1
Loop diagram for SLE. k1 is the incoming field, and k2 is the signal field. Note that the interactions are time ordered within each strand, but not between strands. s1, s2, s3 are the delay times between the interactions along the loop.


where an additional prefactor of 2 was added to take into account the symmetry with respect to the dummy variables τ1 and τ2.

If desired, this can further be decomposed into fully time-ordered terms using the double sided diagrams of Fig. 2. As indicated above we can now e.g. shift the arrow in Fig. 1 corresponding to time τ2 along the ket, thus obtaining 3 possible relative positions with respect to the interactions at times τ and τ1. The loop diagram consequently splits into 3 double sided Feynman diagrams. These are depicted in Fig. 2.

FIG. 2
Time-ordered Feynman diagrams of SLE, generated by shifting the arrows of Fig. 1 along each strand thus changing their relative time ordering. Each of the possible three relative positions then gives one fully time ordered diagram (double sided Feynman ...

The resulting fully-time-ordered expression for the signal reads


Here, we have made use of Eq. (7) to evaluate the correlation functions for the incoming field, since


Note that, in Eq. (35) full time-ordering is expressed by adding a product of step functions to each of the contributing terms. We apply the short hand notation θ(τ1τ2) = θ(τ1τ2).

Next, we calculate the signal for a frequency-domain experiment with stationary beams. This will also illustrate how to extend the diagrammatic rules to frequency-domain experiments. We adopt the following convention for the Fourier transform of a function f


When the incoming field is stationary, the integrals in Eq. (34) can be evaluated directly, which by application of Eq. (37) leads to a multiple Fourier transform of the partially ordered 4-point system correlation function in Eq. (34). To simplify this calculation, it is advantageous to switch to a new set of time variables (“s-variables”), s1:= tτ2, s2:= tτ, s3:= ττ1, which represent the time intervals of interactions along the loop. This choice becomes intuitive when looking at the diagram, Fig. 1.

We can now show that the correlation function in Eq. (35) depends on the s-variables in a simple way. To this end, we transform the system correlation function in Eq. (34) to a form that only involves L-superoperators, hence


In terms of the s-variables, Eq. (38) becomes


where we define the Liouville space propagator


Note that, when associating all interactions with the ket, one needs to propergate twice backwards in time. This amounts to applying the retarded propagator An external file that holds a picture, illustration, etc.
Object name is nihms212364ig8.jpg twice.

Using Eq. (39) the integration in Eq. (34), leads to a multiple Fourier transform of the correlation function in terms of the loop (s)- variables. This is readily calculated, giving the following expression for the frequency-domain SLE signal


Since all interactions in Eq. (41) are “left”-operations that act on the ket, it is possible to express the final result for the signal in Hilbert space, giving


Upon conversion of Eq. (41), the frequency-domain Liouville space propagators,


are replaced by Hilbert-space propagators


Here, ωg is the material frequency of the ground state, which accounts for the free evolution of the bra. In both Eqs. (43) and (44) the infinitesimal η > 0 arises from causality and guarantees the convergence of the Fourier transform.

Equation (42) can be recast in the form


which by using the level scheme of Fig. 3(a) yields the Kramers-Heisenberg formula

FIG. 3
The three level systems with sequential dipole couplings used in the derivation of the Kramers Heisenberg relation for the SLE (panel (a); Eq. (47)) and the pump-probe signal (panel (b)).




In Eqs. (46) and (47) ωca = ωcωa denotes the transition frequency between the levels c and a and we have added a phenomenological dephasing rate Γca.

By proceeding along the same line for an arbitrary loop diagram, one can establish the following rules, that allow to translate a given diagram into its frequency-domain expression. These complement the rules given earlier for time-domain expressions.

  1. In the frequency-domain the loop translates into an alternating product of interactions (arrows) and periods of free evolutions (vertical solid lines) along the loop.
  2. Since the loop time goes clockwise along the loop, periods of free evolution on the left amount to propagating forward in real time (Ĝ), whereas evolution on the right corresponds to backward propergation (Ĝ).
  3. Each Ĝ adds the multiplicative factor i, whereas each Ĝ results in a multiplication by (−i).
  4. Frequency arguments of each propagator are cumulative, i.e. they are given by the sum of all “earlier” interactions along the loop. The ground state frequency ωg must be added to all Green’s functions.

Equation (42) can be immediately generated from Fig. 4 by applying these rules.

FIG. 4
The eight loop diagrams for the pump-probe signal (Eq. (53)). Note that diagram (a) coincides with the SLE (Fig. 1).


The simplest nonlinear optical technique involves two fields, k1 (the pump) and k2 (the probe). The signal is defined as the difference of the probe transmitted intensity with and without the pump. We assume that the probe intensity is high so that spontaneous emission is negligible compared to stimulated emission. The second term in Eqs. (29) and (30), which describes spontaneous emission, is thus neglected. Using Eq. (29) and expanding to second order in the incoming field, the signal then reduces to


Using the identity


we can express the pump-probe signal in terms of +/− operators as


This is equivalent to the classical expression for the signal in terms of χ(3). This becomes even more apparent if Eq. (50) is rewritten in the form


which is equivalent to Eq. (50), if only terms proportional to the field intensities are taken into account. The details are given in appendix A.

The forgoing example of SLE illustrated the use of Keldysh-Schwinger loops to compute the signal. In particular, by using this partially time ordered diagrammatic representation, we were able to reduce the number of diagrams from three to one. The merits of this compact notation will become even more apparent for the pump-probe signal, as will be shown next.

We start with an expression for the signal using Eq. (30)


where in analogy to Eq. (48), we neglect the spontaneous emission term.

We note that the second term in Eq. (52) coincides with the SLE correlation function, Eq. (33). What remains therefore is to resolve the first term in terms of L/R operations. For simplicity we assume a three level ladder with sequential transition dipole moments as shown in Fig. 3(b).

We again start our analysis by looking at the two interactions at times t and τ. Contrary to the second term in Eq. (52), both of these correspond to absorptions rather than emissions. For the dipole selection rules shown in Fig. 3(b), this leaves two possibilities for placing the remaining interactions along the loop: We either have three interactions on one side, or two interactions on either side (see Fig. 4).

For the former case and the model of Fig. 3(b), these can only correspond to two absorbtions and one emission. By applying the same argument used for SLE, the earliest of these three interactions must necessarily be an absorption. To proceed further, one then needs to distinguish whether these interactions occur on the ket or the bra. Since the last interaction is fixed to time t, the former case correspondingly allows only one possible loop diagram, which is displayed in panel (b) of Fig. 4.

When placing the three integrations on the right strand, one has to allow for all possible orderings of their associated times τ, τ1 and τ2. At this point we reiterate that when drawing a loop diagram, even though the relative time ordering of the ket and the bra interactions is unspecified, the temporal order within each strand is fixed. Hence having three interactions with the bra, generates 4 loop diagrams, which arise upon performing permutions of the arrows on the right side. Panels (e) – (h) of Fig. 4 show these diagrams. Note that in this case the loop diagrams are actually fully time ordered.

We now turn to the diagrams with two left and two right interactions. Since the last interaction occurs on the ket at time t, for a system starting off in the ground state, the other interaction on the ket also must be an absorption. Finally, we note that when calculating the correlation functions Eq. (52), one needs to take a trace in the end, which restricts both interactions on the bra to be absorptions as well. Arguing along this line leads to the two loop diagrams shown in panels (c) and (d) of Fig. 4. The two loops again take into account the two possible permutations of placing interactions on the bra.

In summary, the pump-probe signal can be represented by eight loop diagrams, seven stem from the first correlation function in Eq. (52), whereas the remaining one is SLE-type. The corresponding expression for the signal derived using the time-domain rules is given in appendix B.

For the sake of completness, we shall compare the QED signal with the semiclassical result, where the pump-probe signal is calculated as a third order response in the direction +k1k1 + k2. in Note that the diagrams in panels (b) – (h) in Fig. 4, correspond to +k1k1k2. This difference can however easily be resolved by taking the complex conjugate of these diagrams, which leaves the expression for the signal in Eq. (B1) invariant. Pictorially, since this amounts to reflecting the arrows through the center line between the ket and bra, we recover the classical combination of wavevectors +k1k1 + k2. As was the case for SLE, we can generate the corresponding fully-time ordered diagrams from the loops displayed in Fig. 4. The resulting 16 double-sided Feynman diagrams are summarized in Fig. 5. In addition, as has been shown in Sec. III for SLE, frequency-domain expressions follow naturally from the loop diagrams, since the field permutations are already built in. Applying the rules given in Sec. III to the diagrams for the pump-probe (Fig. 4), we can immediately write down the frequency-domain signal

FIG. 5
The Double sided Feynman diagrams resulting from Fig. 4 upon breaking the loops into time ordered contributions. Panel (a) – (d) are the time ordered diagrams corresponding to the loops shown in panel (a) – (d) of Fig. 4, respectively. ...

The eight terms correspond respectively to the eight diagrams in Fig. 4.


So far we have focused on incoherent processes. We could thus consider a single molecule and simply multiply the signal by N in the end. To describe coherent signals we consider V in Eq. (14) as the dipole moment of a collection of noninteracting molecules at positions rα, i.e.


where Vα now denotes the dipole operator of a single molecule. Subsequent interactions of the detecting field can take place with the same (α = β) or with different (αβ) molecules. We shall describe an (n + 1)-wave mixing experiment, where n incoming modes generate a signal in a new direction. When the signal mode (index s) is initially in its vacuum state, the spontaneous emission signal can be described by Eq. (29) modified to include Eq. (54). We then get


where we have further added a time-integration to describe a frequency-domain experiment [37]. Equation (55) is a general expression for spontaneously generated signals. The α = β and αβ terms in Eq. (55) give rise to incoherent (SI) and coherent (SC) signals, respectively. We thus write


where the incoherent signal is given by


For uncorrelated particles, we have [Vα, Vβ] = 0. The αβ sum may thus be factorized into α and β,


This gives rise to the coherent term in Eq. (56),


For this we have also made use of


in Eq. (57).

In Eq. (59) we assumed that the sample is much smaller than the optical wavelength or that we have exact phase matching Δk = 0, where Δk is the difference between the wavevector of the signal mode and the sum of wavevectors of the n incoming fields. More generally, we should replace the N(N − 1) factor in Eq. (56) by


For macroscopic samples, the number of terms with αβ in the double sum of Eq. (61) will by far exceed the terms with equal indices. Evaluating this double sum in the continuous limit and letting the sample volume go to infinity, gives the standard phase matching condition Fk) → δk), which yields a directed signal.

For small samples both coherent and incoherent terms need to be considered. It is interesting to note that for a single molecule we never get SC, we only have SI. Nonlinear susceptibilities, even though they are calculated for single molecules, do not represent a spontaneous wave-mixing experiment performed on a single molecule. Using the semiclassical approach one may conclude that coherent signals are possible even from a single molecule. The QED approach shows that this is not possible.

Below we compare the expressions for SI(n) and P(n) obtained by expanding Eqs. (57) and (59) for n = 1, 2, 3, where n denotes the number of incoming modes. We demonstrate that the incoherent signal for n incoming modes, relates to a 2(n + 1)-point correlation function in the dipole moment, compared to an (n + 1)-point correlation function contained in P(n).

For n = 1


gives Rayleigh scattering, whereas


is responsible for SLE (Raman and Fluorescence) as discussed in Sec. III.

For n = 2


describes 3-wave mixing e.g. second harmonic generation [41]–[43], and


may represent many possible processes, depending on the level scheme and off-resonant detunings, e.g. two photon fluorescence [38]–[40].

For n = 3


represents 4-wave mixing e.g. third harmonic generation and


can represent many possible processes, e.g. three photon fluorescence.

A nanocrystal is made out of many unit cells. Its dipole operator is given by [mu] = Σα[mu]α, where [mu]α represents an electron-hole pair on unit α. It thus behaves as a collection of many molecules and it can show spontaneous coherent response [8]–[9]. N is thus large even though we have a single particle.


When detecting nonlinear optical signals in macroscopic samples, one can either measure the intensity (homodyne detection) or the magnitude and phase of the signal field (heterodyne detection). The later technique uses an intense external field An external file that holds a picture, illustration, etc.
Object name is nihms212364ig9.jpg, usually referred to as “local oscillator”, which interferes with the created nonlinear signal in the same direction, but is assumed not to interact with the molecule. Using the standard semiclassical approach to nonlinear optics, the quantity detected in a heterodyne measurement is [3]


where Ps is the polarization in the signal direction ks. With the QED formalism developed here, we can greatly simplify the derivation of Eq. (62) and show that it emerges naturally as an incoherent rather than a coherent signal.

Equation (55) describes the signal generated from an assembly of particles by spontaneous emission. Applying Eq. (54) to Eq. (29), we can give a similar expression for the stimulated signal


In analogy to Sec. (V), we can now split Eq. (63) into a coherent (αβ) and an incoherent (α = β) contribution. Since the dipole operators of different, uncorrelated particles αβ commute, the coherent contribution vanishes identically. The present formalism hence yields the general result that stimulated processes are always incoherent in nature. In particular this applies to heterodyne experiments where the local oscillator serves as the stimulating field. It is hence sufficient to look only at a single particle and use Eq. (14) derived in Sec. II.

Equation (14) is reminiscent of the expression for heterodyne detection (Eq. (62)), if one only looks at one particular mode “s” in the definition of the signal (see Eqs. (11) and (12)), thus


Here, by adding a subscript α, we make explicit that Eq. (64) is the signal generated by one molecule at position rα. For an assembly of noninteracting particles, this has then to be summed over rα. The only difference of Eq. (64) as compared to Eq. (62), is that it contains a field operator An external file that holds a picture, illustration, etc.
Object name is nihms212364ig1.jpg, rather than a classical field An external file that holds a picture, illustration, etc.
Object name is nihms212364ig7.jpg.

Heterodyne detection can be viewed as an incoherent stimulated emission process in the detected mode. Since Eq. (64) is already first order in the mode s, in a subsequent perturbative expansion of the density operator the detecting mode will only contribute to higher order. Hence, to first order in the coupling with the detecting mode, the time-dependence of the density operator is only given by all other modes s′ ≠ s. The field operator therfore acts directly on the state of the system (molecule + field) at t → −∞, which when assuming a coherent state of the field in the same limit (see Eq. (7)) yields


In Eq. (65) we additionally took the complex conjugate within the imaginary part.

The expectation value in Eq. (65) can be further expanded in the incoming modes. For n modes this gives the classical nth order polarization


Consequently, integrating over all times, we obtain the following generic result for an (n+1)-wave mixing process (i.e. 1 detected mode “s” + n incoming modes “s′ ≠ s”)


where we explicitly display the spatial dependence of the fields, e.g.


In Eq. (67) we have added the sum over rα to make the transition from the signal generated by single molecule to one of several molecules (incoherent process).

Under exact phase-matching conditions (Δk = 0), or if the sample is much smaller than the optical wavelength, we have fk) = N. Rather than creating a coherent signal, propagating it and interfering as done in the semiclassical approach, we can simply describe it microscopically as an incoherent signal.

Equation (67), hence allows us to interpret heterodyne detection as follows: All fields (the n fields (ωs, ks) as well as the detected field (ωs, ks)) interact with the molecule. This leads to stimulated emission of photons of frequency ωs. A sum over all molecular positions has to be performed (Eq. (67)). Going to the continuum limit, which is justified for a macroscopic sample, this sum can be replaced by an integral over the sample volume. Extending this integral over the entire real space, we obtain the phase-matching condition fk) → δk). Therefore, in a macroscopic sample the signal is finite only when the stimulating field An external file that holds a picture, illustration, etc.
Object name is nihms212364ig10.jpg(r, t) fulfills this condition.


A microscopic QED treatment of nonlinear optical processes induced by a weak quantized field is developed in this paper. In the traditional, widely-used approaches to nonlinear optical response, the matter-field interaction is treated semiclassically, hence the signal field is obtained by simultaneously solving macroscopic (Maxwell-) and microscopic (quantum Liouville-) equations. Here we treat the entire process as a single event.

Based on a microscopic definition of the spontaneous and stimulated signal, general expressions are derived for the nonlinear optical signal involving non-equilibrium Green’s functions in the incoming field. Depending on the experiment under consideration, these can then be expanded order by order in the incoming field modes. Application is made to spontaneous light emission and pump-probe spectroscopy. Coherent and incoherent processes involving both spontaneous and stimulated emission are treated using a unified framework.

A diagrammatic derivation of the contributing terms within the RWA using the technique of Keldysh Schwinger loop is developed. This respresentation reflects the partially time-ordered nature of the Green’s functions and hence yields more compact expressions than the well established fully time ordered double sided Feynman diagrams. The loop diagrams are particularly useful for frequency-domain measurements where the bookeeping of time ordering is not necessary anymore. If needed for time-domain experiments, the fully time ordered expressions can be obtained from the loop diagrams using a simple prescription. Rules for constructing and reading these diagrams are laid out, making it possible to intuitively derive the expressions for the signals.

The practical merits of this partially time-ordered approach are illustrated for the pump-probe technique where the number of diagrams is reduced from 16 to 8. This example also demonstrates an interesting fundamental merit of the QED treatment of nonlinear optical processes, even though in pump-probe spectroscopy quantum effects of the field (i.e. terms stemming from spontaneous emission) are usually negligible: Semiclassically, the pump-probe signal is obtained by calculating a third order response induced by two fields in the direction k1k1 +k2. The resulting signal field is then obtained using heterodyne detection, where mode 2 acts as its own local oscillator field. Even though the final expression for the signal is identical to a QED result, the semiclassical approach creats an artificial asymmetry between the local oscillator (mode 2), which by definition does not interact with the molecule, and the remaining fields (mode 1 and 2), which give rise to the nonlinear polarization. This asymmetry, which is eliminated here, is a direct consequence of the limitations imposed by separately treating matter and field within the semiclassical approach.

Finally, using the QED formalism, we are able to treat both coherent and incoherent processes in a unified way. This distinction comes about by calculating the signal for an assembly of molecules rather than for a single particle. The signal then splits into two parts, a coherent term scaling ~ N(N − 1) and an incoherent term ~N, where N denotes the number of particles. As the number of particles in the sample is increased, the coherent term becomes more pronounced. For experiments carried out with only a few molecules however, both coherent and incoherent terms will generally make comparable contributions to the signal. In single molecule experiments only the incoherent term survives.

Since incoherent signals may not be calculated semiclassically, this leads to a striking consequence. For an (n +1)-wave mixing processes carried out on a single molecule we show that the signal is not given by an (n + 1)-point response function of the dipole operator related to P(n), as predicted by the semiclassical approach. It is rather given by a different 2(n + 1)-point combination of correlation functions. This doubling arises since the signal may not be recast as an amplitude square and we must calculate the signal itself, not an amplitude. The present formalism allows a microscopic calculation of nonlinear optical experiments on single molecules which have only become feasible recently.

Finally, we address this apparent limitation of the semiclassical description from a more general viewpoint and show that heterodyne detection for an (n+1)-wave mixing experiment can be simply viewed as an incoherent stimulated emission process in the detected mode. In contrast, homodyne-detected n-wave mixing is a coherent spontaneous emission process. For a macroscopic sample, by summing over the molecules in the interaction volume, we recover the phase-matching condition, i.e. the heterodyne-signal is finite only if the stimulating mode (i.e. the local oscillator) matches this condition.


We wish to thank Prof. Sunney Xie for useful discussions. The support of the National Science Foundation (Grant No. CHE-0446555, CBC-0533162), NIRT (Grant No. EEC 0303389)) and the National Institutes of Health (GM59230) is gratefully acknowledged.


In this appendix we illustrate the equivalence of the QED expression for the pump-probe signal with the classical signal expressed in terms of χ(3). Starting from Eq. (50) and applying the definition of Eq. (16), we can write


Making use of Eq. (A1), the correlation function in Eq. (50) can be recast in terms of correlations of the system and the field. Calculating the field correlations, only 6 of the possible 16 terms are nonzero, giving


where the system is assumed to be in its ground state for t → −∞, hence


Note that since Eq. (50) is symmetric with respect to the variables τ1 and τ2, Eq. (A2) simplifies further and only includes three distinct terms upon integration, i.e.


For a classical incoming field, the correlation function in Eq. (A2) is equal to only the first term on the rhs, where


For a classical field the second term in Eq. (A5) will be negligible.

The last two terms in Eq. (A4) are due to the quantum character of the incoming field, since


Hence, neglecting the last two terms in Eq. (A4) and making use of Eq. (A5), we obtain Eq. (51).


In this appendix we give the expressions for the pump-probe signal in the time-domain obtained by applying the rules summarized in Sec. III to the loop diagrams in Fig. 4. This complements the frequency-domain expression given in Eq. (53). Applying the time-domain rules for the loop diagrams in Fig. 4, we obtain


where we kept the order of the diagrams in Fig. 4. Similar to Eq. (34), an additional factor 2 accounts for the symmetry with respect to τ1 and τ2. This can again be broken into fully time-ordered terms using the double sided Feynman diagrams displayed in Fig. 5. We then get


Comparison of Eqs. (B1) and (B2) illustrates the tremendous simplification achieved by the more compact partially-time-ordered approach. Employing the diagrammatic Keldysh-Schwinger loop technique, the number of diagrams is reduced by half (16 vs. 8 diagrams).


PACS numbers: 42.65-k, 12.20.-m, 42.62. Fi, 42.50.-p, 42.50.Ct


1. Andrews D, Allcock P. Optical Harmonics in Molecular Systems; Quantum Electrodynamical Theory. Wiley-VCH; Weinheim: 2002.
2. Bloembergen N. Nonlinear optics. Benjamin; New York: 1965.
3. Mukamel S. Principles of Nonlinear Optical Spectroscopy. Oxford University Press; New York: 1995.
4. Scully M, Zubairy M. Quantum Optics. Cambridge University Press; Cambridge: 1997.
5. Hanamura E, Mukamel S. J Opt Soc Am B. 1986;3:1124–1126.
6. Mukamel S. Adv Chem Phys. 1988;70(Part I):165–230.
7. Bidault S, Brasselet S, Zyss J. Optics Letters. 2004;29:1242. [PubMed]
8. Van Dijk MA, Lippitz M, Orrit M. Accounts of chemical research. 2005;38:594–601. [PubMed]
9. Van Dijk MA, Lippitz M, Orrit M. Phys Rev Lett. 2005;95:267406. [PubMed]
10. Sanchez EJ, Novotny L, Holtom GR, Xie XS. Phys Chem A. 1997;101:7019.
11. Muskens OL, Del Fatti N, Vallee F. Nano Letters. 2006;6:552. [PubMed]
12. Langbein W, et al. J Phys¿ Condens Matter. 2007;19:295203. [PubMed]
13. Lukin MD, et al. Science. 2007;316:1312. [PubMed]
14. Potma EO, Xie XS. Chem Phys Chem. 2005;6:77. [PubMed]
15. Haroche S, Raimond JM. Exploring the Quantum: Atoms, Cavities, and Photons. Oxford University Press; Oxford: 2006.
16. Schleich WP. Quantum Optics in Phase Space. Wiley-VCH; Berlin: 2001.
17. Scully M, Lamb WE., Jr Physical Review. 1967;159:159.
18. Sarget M, III, Holm DA, Zubairy MS. Physical Review A. 1985;31:3112. [PubMed]
19. Stenholm S, Holm DA, Sarget M., III Physical Review A. 1985;31:3124. [PubMed]
20. Reid MD, Walls DF. Physical Review A. 1986;34:4929. [PubMed]
21. Agarwal GS. Physical Review A. 1986;34:4055. [PubMed]
22. Raymer MG, Walmsley IA. Progress in Optics. 1990;28:181.
23. Gorbachev VN, Zhiliba AI. J Opt B. 2000;2:764.
24. Panteleev AA, Rerikh VlK, Starostin AN. Journal of Experimental and Theoretical Physics. 2000;90:50.
25. Macovei M, Li GX. Physical Review A. 2007;76:023818.
26. Negele J, Orland H. Quantum Many Particle Physics. Westview Press; 1998.
27. Rammer J. Quantum field theory of non-equilibrium states. Cambridge University Press; 2007.
28. Keldysh LV. Sov Phys JETP. 1965;20:1018.
29. Gorkov LP, Dzyaloshinski IE, Abrikosov AA. Methods of quantum field theory. Prentice Hall; 1963.
30. Mills R. Propagators for many-particle systems; an elementary treatment. New York: Gordon and Breach; 1969.
31. Datta S. Electronic Transport in Mesoscopic Systems. Cambridge University Press; Cambridge: 1995.
32. Haug H, Jauho A-P. Quantum Kinetics in Transport and Optics of Semiconductors. Springer-Verlag; Berlin, Heidelberg: 1996.
33. Harbola U, Maddox JB, Mukamel S. Phys Rev B. 2006;73:205404.
34. Harbola U, Mukamel S. The Journal of Chemical Physics. 2006;124:044106. [PubMed]
35. Harbola U, Maddox JB, Mukamel S. Phys Rev B. 2006;73:075211.
36. Mukamel S. submitted to Phys. Rev. A.
37. Chernyak V, Wang N, Mukamel S. Physics Reports. 1995;263:213–309.
38. Callis PR. Annu Rev Phys Chem. 1997;48:271. [PubMed]
39. Pons T, Mertz J. Optics Communications. 2006;258:203.
40. Rehms AA, Callis PR. Chemical Physics Letters. 1993;208:276.
41. Campagnola PJ, Loew LM. Nature Biotechnology. 2003;21:1356. [PubMed]
42. Moreaux L, Pons T, Dambrin V, Blanchard-Desce M, Mertz J. Optics Letters. 2003;28:625. [PubMed]
43. Naujok RR, Paul HJ, Corn RM. The Journal of Physical Chemistry. 1996;100:10497.