Nonlinear response theory on the Keldysh contour

The correction to the dipole moment to linear order in the perturbation is given by

$$\begin{equation} \langle \hat{\mu }^{(1)}(t) \rangle = - \frac{1}{{\rm i}\hbar } \int _c {\rm d}\tau _1 {\rm Tr}\lbrace {\rm T_c}[ \hat{\mu }(\tau _1) \hat{\mu }(t) ] \hat{\varrho }(t_0) \rbrace E(\tau _1). \end{equation} \tag{ 11 }$$

The contour time ordering of the two dipole operators, $\hat{\mu }(\tau _1)$ and $\hat{\mu }(t)$, is written out explicitly using Heaviside step functions

$$\begin{eqnarray} &&{\rm T_c}[ \hat{\mu }(\tau _1) \hat{\mu }(t) ] \,{=}\, \vartheta (t \,{-}\, \tau _1) \hat{\mu }(t) \hat{\mu }(\tau _1) \,{+}\, \vartheta (\tau _1 - t) \hat{\mu }(\tau _1) \hat{\mu }(t), \nonumber\\ \end{eqnarray} \tag{ 12 }$$

and the linear response term now splits into two:

$$\begin{eqnarray} &&\fl \langle \hat{\mu }^{(1)}(t) \rangle = - \frac{1}{{\rm i}\hbar } \int _{t_0}^{t} {\rm d}\tau _1 {\rm Tr}\lbrace \hat{\mu }(t) \hat{\mu }(\tau _1) \hat{\varrho }(t_0) \rbrace E(\tau _1) \end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray} &&\fl \hphantom{\langle \hat{\mu }^{(1)}(t) \rangle =} - \frac{1}{{\rm i}\hbar } \int _{t}^{t_0} \! {\rm d}\tau _1 {\rm Tr}\lbrace \hat{\mu }(\tau _1) \hat{\mu }(t) \hat{\varrho }(t_0) \rbrace E(\tau _1). \end{eqnarray} \tag{ 14 }$$

The two loop diagrams contributing to linear response (representing equations (13) and (14)) are depicted in figure 2. (Strictly speaking, we identify each diagram with a Liouville space pathway; see below.) Physical time runs on the real axis, where we focus on the interval between times t0 and t. The contour introduced via the contour-ordered exponential comprises a closed loop in the complex time plane (see figure 1). We take contour 'time' to begin at t0 and to propagate clockwise around the contour. The contour crosses the real axis at the observation time t and propagates back to t0. The first term, equation (13), has τ1 before t on the contour; thus, the scattering event occurs on the upper branch of the loop contour and is represented by the diagram to the left in figure 2. Conversely, term two has the scattering occurring at a point τ1 occurring after t, and τ1 is thus marked on the lower branch in the right diagram in figure 2. To project contour time onto physical time, we follow, in spirit, the strategy laid out by Langreth [4]. Under the general assumption that the integrand has no poles within the loop diagram, we invoke Cauchy's integral theorem to continuously deform the contour. Having pinned the contour to the real axis at times t0 and t, we tighten the loop contour. On the upper branch of the contour, contour time morphs continuously into physical time and we obtain the limits ∫cdτ → ∫tt0dτ on the time integrals. For the lower time contour, we obtain ∫cdτ → ∫t0tdτ. Flipping the limits of the time integration on the lower branch at the cost of a sign ($\int _{t}^{t_0} {\rm d}\tau = - \int _{t_0}^{t} {\rm d}\tau$), we could obtain Kubo's expression from equations (13) and (14):

$$\begin{eqnarray} &&\fl \langle \hat{\mu }^{(1)}(t) \rangle = - \frac{1}{{\rm i}\hbar } \int _{t_0}^{t} {\rm d}\tau _1 {\rm Tr}\lbrace [ \hat{\mu }(t), \hat{\mu }(\tau _1)] \hat{\varrho }(t_0) \rbrace E(\tau _1). \end{eqnarray} \tag{ 15 }$$

We note that by working with the loop contour we avoid the need for commutators. Disregarding correlations in the initial state we let t₀ → −∞, also we extend the integration to infinity by introducing a Heaviside step function. Substituting the time variable t₁ = t − τ₁,

$$\begin{eqnarray} &&\fl \langle \hat{\mu }^{(1)}(t) \rangle = - \frac{1}{{\rm i}\hbar } \int \! {\rm d}t_1 \vartheta (t_1) {\rm Tr}\lbrace \hat{\mu }(t) \hat{\mu }(t - t_1) \hat{\varrho }(-\infty ) \rbrace \nonumber\\ &&\hphantom{ \langle \hat{\mu }^{(1)}(t) \rangle =} \times E(t - t_1) \end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray} &&\fl \hphantom{ \langle \hat{\mu }^{(1)}(t) \rangle =} + \frac{1}{{\rm i}\hbar } \int \! {\rm d}t_1 \vartheta (t_1) {\rm Tr}\lbrace \hat{\mu }(t - t_1) \hat{\mu }(t) \hat{\varrho }(-\infty ) \rbrace \nonumber\\ &&\hphantom{ \langle \hat{\mu }^{(1)}(t) \rangle =} \times E(t - t_1). \end{eqnarray} \tag{ 17 }$$

Define the linear response function, S⁽¹⁾(t₁), indirectly, via the convolution integral

$$\begin{equation} \langle \hat{\mu }^{(1)}(t) \rangle = \int \!{\rm d}t_1 S^{(1)}(t_1)E(t-t_1), \end{equation} \tag{ 18 }$$

which gives us the expression

$$\begin{eqnarray} \fl S^{(1)}(t_1) &=& \frac{{\rm i}}{\hbar } \vartheta (t_1) \left [{\rm Tr}\lbrace \hat{\mu }(t) \hat{\mu }(t - t_1) \hat{\varrho }(-\infty ) \rbrace \right. \nonumber\\ &&\left. - {\rm Tr}\lbrace \hat{\mu }(t - t_1) \hat{\mu }(t) \hat{\varrho }(-\infty ) \rbrace \right ] . \end{eqnarray} \tag{ 19 }$$

The dependence on t is only apparent, since when disregarding correlations in the initial state, we secured time translation invariance, and we can choose t − t₁ = 0 to obtain

$$\begin{eqnarray} \fl S^{(1)}(t_1) &=& \frac{{\rm i}}{\hbar } \vartheta (t_1) [{\rm Tr}\lbrace \hat{\mu }(t_1) \hat{\mu }(0) \hat{\varrho }(-\infty ) \rbrace \nonumber \\ && - {\rm Tr}\lbrace \hat{\mu }(0) \hat{\mu }(t_1) \hat{\varrho }(-\infty ) \rbrace ]. \end{eqnarray} \tag{ 20 }$$

The Heaviside function explicitly ensures causality. The two terms are complex conjugated. Defining $J(t_1) = {\rm Tr}\lbrace \hat{\mu }(t_1) \hat{\mu }(0) \hat{\varrho }(-\infty ) \rbrace ,$ we can write

$$\begin{equation} S^{(1)}(t_1) = \frac{{\rm i}}{\hbar } \vartheta (t_1) [J(t_1) - J^*(t_1)] . \end{equation} \tag{ 21 }$$

Complex conjugation corresponds to the reversal of contour time; for the loop diagrams this amounts to reflection in the real axis.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The second-order correction to the dipole moment is given by

$$\begin{eqnarray} \fl \langle \hat{\mu }^{(2)}(t) \rangle &=& \frac{1}{2} \left( \frac{1}{{\rm i}\hbar } \right)^2 \int _c \! {\rm d}\tau _2 \int _c \! {\rm d}\tau _1 \nonumber\\ &&\times \,{\rm Tr}\lbrace {\rm T_c}[\hat{\mu }(\tau _2) \hat{\mu }(\tau _1) \hat{\mu }(t)]\hat{\varrho }(t_0) \rbrace E(\tau _2) E(\tau _1) . \nonumber\\ \end{eqnarray} \tag{ 22 }$$

Expansion of the contour ordering of n + 1 operators gives (n + 1)! terms (in this case 3! = 6):

$$\begin{eqnarray} &&{\rm T_c} [\hat{\mu }(\tau _2) \hat{\mu }(\tau _1) \hat{\mu }(t)] = \vartheta (t - \tau _2) \vartheta (\tau _2 - \tau _1)\ \hat{\mu }(t) \hat{\mu }(\tau _2) \hat{\mu }(\tau _1)\nonumber \\ \end{eqnarray} \tag{ 23 }$$

$$\begin{equation} \;\;+ \vartheta (t - \tau _1) \vartheta (\tau _1 - \tau _2)\ \hat{\mu }(t) \hat{\mu }(\tau _1) \hat{\mu }(\tau _2) \end{equation} \tag{ 24 }$$

$$\begin{equation} + \vartheta (\tau _2 - t) \vartheta (t - \tau _1)\ \hat{\mu }(\tau _2) \hat{\mu }(t) \hat{\mu }(\tau _1) \end{equation} \tag{ 25 }$$

$$\begin{equation} + \vartheta (\tau _1 - t) \vartheta (t - \tau _2)\ \hat{\mu }(\tau _1) \hat{\mu }(t) \hat{\mu }(\tau _2) \end{equation} \tag{ 26 }$$

$$\begin{equation} \;\;+ \vartheta (\tau _2 - \tau _1) \vartheta (\tau _1 - t)\ \hat{\mu }(\tau _2) \hat{\mu }(\tau _1) \hat{\mu }(t) \end{equation} \tag{ 27 }$$

$$\begin{equation} \;\;\;+ \vartheta (\tau _1 - \tau _2) \vartheta (\tau _2 - t)\ \hat{\mu }(\tau _1) \hat{\mu }(\tau _2) \hat{\mu }(t) . \end{equation} \tag{ 28 }$$

We group these terms according to the position of $\hat{\mu }(t)$. Consider first the n! (in this case 2! = 2) terms where t occurs first on the contour (we will refer to these terms as III):

$$\begin{eqnarray} &&\fl \frac{1}{2} \left( \frac{1}{{\rm i}\hbar } \right)^2 \int _c \! {\rm d}\tau _2 \int _c \! {\rm d}\tau _1 \vartheta (\tau _2 - \tau _1) \vartheta (\tau _1 - t) \nonumber\\ &&\times {\rm Tr}\lbrace \hat{\mu }(\tau _2) \hat{\mu }(\tau _1) \hat{\mu }(t) \hat{\varrho }(t_0) \rbrace E(\tau _2) E(\tau _1) \end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray} &&\fl + \frac{1}{2} \left( \frac{1}{{\rm i}\hbar } \right)^2 \int _c \! {\rm d}\tau _2 \int _c \! {\rm d}\tau _1 \vartheta (\tau _1 - \tau _2) \vartheta (\tau _2 - t) ) \nonumber\\ &&\times {\rm Tr}\lbrace \hat{\mu }(\tau _1) \hat{\mu }(\tau _2) \hat{\mu }(t) \hat{\varrho }(t_0) \rbrace E(\tau _2) E(\tau _1) . \end{eqnarray} \tag{ 30 }$$

When we tighten the loop contour the contour 'times' τ₁ and τ₂ lie on the lower, anti-time-ordered branch of the Keldysh contour.

$$\begin{eqnarray} &&\fl =\frac{1}{2} \left( \frac{1}{{\rm i}\hbar } \right)^2 \int _{t}^{t_0} \! {\rm d}\tau _2 \int _{t}^{t_0} \! {\rm d}\tau _1 \vartheta (\tau _1 - \tau _2) \nonumber\\ &&\times {\rm Tr}\lbrace \hat{\mu }(\tau _2) \hat{\mu }(\tau _1) \hat{\mu }(t) \hat{\varrho }(t_0) \rbrace E(\tau _2) E(\tau _1) \end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray} &&\fl + \frac{1}{2} \left( \frac{1}{{\rm i}\hbar } \right)^2 \int _{t}^{t_0} \! {\rm d}\tau _2 \int _{t}^{t_0} \! {\rm d}\tau _1 \vartheta (\tau _2 - \tau _1)\nonumber\\ &&\times {\rm Tr}\lbrace \hat{\mu }(\tau _1) \hat{\mu }(\tau _2) \hat{\mu }(t) \hat{\varrho }(t_0) \rbrace E(\tau _2) E(\tau _1) . \end{eqnarray} \tag{ 32 }$$

For each of the k (here 2) contour 'times' on the lower branch, we flip time, i.e. the limits of the 'time' integration, at the cost of a sign. We obtain an overall factor of ( − 1)^k (here 1). The n! terms now differ only by the permutation of the integration variables. We choose one ordering and eliminate the factor of $\frac{1}{n!}$. To connect to experiments, we identify the specific ordering of the dummy time variables with ordering in physical time:

$$\begin{eqnarray} &&\fl = \left( \frac{1}{{\rm i}\hbar } \right)^2 \int _{t_0}^{t} \! {\rm d}\tau _2 \int _{t_0}^{t} \! {\rm d}\tau _1 \vartheta (\tau _2 - \tau _1) \nonumber\\ &&\times {\rm Tr}\lbrace \hat{\mu }(\tau _1) \hat{\mu }(\tau _2) \hat{\mu }(t) \hat{\varrho }(t_0) \rbrace E(\tau _2) E(\tau _1) . \end{eqnarray} \tag{ 33 }$$

We let the reference time t₀ go to the infinite past $\hat{\varrho }(t_0) \rightarrow \hat{\varrho }(-\infty )$ and introduce a second step function ϑ(t − τ₂):

$$\begin{eqnarray} &&\fl \left( \frac{1}{{\rm i}\hbar } \right)^2 \int \! {\rm d}\tau _2 \int \! {\rm d}\tau _1 \vartheta (t - \tau _2) \vartheta (\tau _2 - \tau _1)\nonumber\\ &&\times {\rm Tr}\lbrace \hat{\mu }(\tau _1) \hat{\mu }(\tau _2) \hat{\mu }(t) \hat{\varrho }(-\infty ) \rbrace E(\tau _2) E(\tau _1) . \end{eqnarray} \tag{ 34 }$$

Now we make the substitution t_i = τ_{i + 1} − τ_i (τ_{n + 1} ≡ t) of the time variables τ_i. The new variables t_i refer to the time intervals between scattering events. We note that $\frac{\partial t_i}{\partial \tau _i} = -1$; hence, the Jacobian of the variable substitution equals unity, |J_n| = |( − 1)ⁿ| = 1:

$$\begin{eqnarray} &&\fl \left( \frac{1}{{\rm i}\hbar } \right)^2 \int \! {\rm d} t_2 \int \! {\rm d} t_1 \vartheta (t_2) \vartheta (t_1) {\rm Tr}\lbrace \hat{\mu }(t - t_2 - t_1) \hat{\mu }(t - t_2) \nonumber\\ &&\times \hat{\mu }(t) \hat{\varrho }(-\infty ) \rbrace E(t - t_2) E(t - t_2 - t_1) . \end{eqnarray} \tag{ 35 }$$

Define the second-order response function indirectly by

$$\begin{eqnarray} &&\fl \langle \hat{\mu }^{(2)}(t) \rangle = \int \!{\rm d}t_2 \int \!{\rm d}t_1 S^{(2)}(t_2,t_1) E(t - t_2) E(t - t_2 - t_1) . \nonumber\\ \end{eqnarray} \tag{ 36 }$$

So far, we have obtained the contribution from the term III:

$$\begin{eqnarray} &&\fl S^{(2)}_{{\rm III}}(t_2,t_1) = C^{(2)} {\rm Tr}\lbrace \hat{\mu }(t - t_2 - t_1) \hat{\mu }(t - t_2) \hat{\mu }(t) \hat{\varrho }(-\infty ) \rbrace \nonumber\\ \end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray} &&\fl \hphantom{S^{(2)}_{{\rm III}}(t_2,t_1)}= C^{(2)}{\rm Tr}\lbrace \hat{\mu }(0) \hat{\mu }(t_1) \hat{\mu }(t_2 + t_1) \hat{\varrho }(-\infty ) \rbrace , \end{eqnarray} \tag{ 38 }$$

where $C^{(n)} = \left( \frac{{\rm i}}{\hbar } \right)^n \prod _{i=1}^{n}\vartheta (t_n)$. Next consider the n! terms with one time before t (we refer to them as II)

$$\begin{eqnarray} &&\fl \frac{1}{2} \left( \frac{1}{{\rm i}\hbar } \right)^2 \int _c \! {\rm d}\tau _2 \int _c \! {\rm d}\tau _1 \vartheta (\tau _2 - t) \vartheta (t - \tau _1)\nonumber \\ &&\times {\rm Tr}\lbrace \hat{\mu }(\tau _2) \hat{\mu }(t) \hat{\mu }(\tau _1) \hat{\varrho }(t_0) \rbrace E(\tau _2) E(\tau _1) \end{eqnarray} \tag{ 39 }$$

$$\begin{eqnarray} &&\fl +\frac{1}{2} \left( \frac{1}{{\rm i}\hbar } \right)^2 \int _c \! {\rm d}\tau _2 \int _c \! {\rm d}\tau _1 \vartheta (\tau _1 - t) \vartheta (t - \tau _2)\nonumber \\ &&\times {\rm Tr}\lbrace \hat{\mu }(\tau _1) \hat{\mu }(t) \hat{\mu }(\tau _2) \hat{\varrho }(t_0) \rbrace E(\tau _2) E(\tau _1) \end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray} &&\fl =\frac{1}{2} \left( \frac{1}{{\rm i}\hbar } \right)^2 \int _{t}^{t_0} \! {\rm d}\tau _2 \int _{t_0}^{t} \! {\rm d}\tau _1 {\rm Tr}\lbrace \hat{\mu }(\tau _2) \hat{\mu }(t) \hat{\mu }(\tau _1) \hat{\varrho }(t_0) \rbrace\nonumber \\ &&\times E(\tau _2) E(\tau _1) \end{eqnarray} \tag{ 41 }$$

$$\begin{eqnarray} &&\fl +\frac{1}{2} \left( \frac{1}{{\rm i}\hbar } \right)^2 \int _{t_0}^{t} \! {\rm d}\tau _2 \int _t^{t_0} \! {\rm d}\tau _1 {\rm Tr}\lbrace \hat{\mu }(\tau _1) \hat{\mu }(t) \hat{\mu }(\tau _2) \hat{\varrho }(t_0) \rbrace \nonumber \\ &&\times E(\tau _2) E(\tau _1) \end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray} &&\fl = - \left( \frac{1}{{\rm i}\hbar } \right)^2 \int _{t_0}^{t} \! {\rm d}\tau _2 \int _{t_0}^{t} \! {\rm d}\tau _1 {\rm Tr}\lbrace \hat{\mu }(\tau _2) \hat{\mu }(t) \hat{\mu }(\tau _1) \hat{\varrho }(t_0) \rbrace\nonumber \\ &&\times E(\tau _2) E(\tau _1) . \end{eqnarray} \tag{ 43 }$$

We arrive at a single term with no relative order of the scattering times τ₁ and τ₂. We split this term in two by introducing two Heaviside functions: 1 = ϑ(τ₂ − τ₁) + ϑ(τ₁ − τ₂). Subsequently, we swap the two dummy time variables in the second term to obtain two terms (IIA and IIB):

$$\begin{eqnarray} &&\fl - \left( \frac{1}{{\rm i}\hbar } \right)^2 \int \! {\rm d}\tau _2 \int \! {\rm d}\tau _1 \vartheta (t - \tau _2) \vartheta (\tau _2 - \tau _1) \nonumber \\ &&\times {\rm Tr}\lbrace \hat{\mu }(\tau _2) \hat{\mu }(t) \hat{\mu }(\tau _1) \hat{\varrho }(-\infty ) \rbrace E(\tau _2) E(\tau _1) \end{eqnarray} \tag{ 44 }$$

$$\begin{eqnarray} &&\fl - \left( \frac{1}{{\rm i}\hbar } \right)^2 \int \! {\rm d}\tau _2 \int \! {\rm d}\tau _1 \vartheta (t - \tau _2) \vartheta (\tau _2 - \tau _1)\nonumber \\ &&\times {\rm Tr}\lbrace \hat{\mu }(\tau _1) \hat{\mu }(t) \hat{\mu }(\tau _2) \hat{\varrho }(-\infty ) \rbrace E(\tau _2) E(\tau _1) . \end{eqnarray} \tag{ 45 }$$

Making time substitutions as above and writing only the contributions to the response function, we have

$$\begin{equation} S^{(2)}_{{\rm IIA}}(t_2,t_1) = - C^{(2)} {\rm Tr}\lbrace \hat{\mu }(t_1) \hat{\mu }(t_2 + t_1) \hat{\mu }(0) \hat{\varrho }(-\infty ) \rbrace \end{equation} \tag{ 46 }$$

$$\begin{equation} S^{(2)}_{{\rm IIB}}(t_2,t_1) = - C^{(2)} {\rm Tr}\lbrace \hat{\mu }(0) \hat{\mu }(t_2 + t_1) \hat{\mu }(t_1) \hat{\varrho }(-\infty ) \rbrace . \end{equation} \tag{ 47 }$$

Defining two of the Liouville space pathways as

$$\begin{equation} Q_1(t_2,t_1) = \hspace{8.5359pt} {\rm Tr}\lbrace \hat{\mu }(t_2 + t_1) \hat{\mu }(t_1) \hat{\mu }(0) \hat{\varrho }(-\infty ) \rbrace , \end{equation} \tag{ 48 }$$

$$\begin{equation} Q_2(t_2,t_1) = - {\rm Tr}\lbrace \hat{\mu }(t_1) \hat{\mu }(t_2 + t_1) \hat{\mu }(0) \hat{\varrho }(-\infty ) \rbrace , \end{equation} \tag{ 49 }$$

the response function can be written as

$$\begin{equation} S^{(2)}(t_2,t_1) = C^{(2)} \sum _{i = 1}^{2}[Q_i(t_2,t_1) + Q_i^*(t_2,t_1) ]. \end{equation} \tag{ 50 }$$

Each of the contributions is given by

$$\begin{equation} S^{(2)}_{{\rm I}}(t_2,t_1) = C^{(2)} Q_1(t_2,t_1), \end{equation} \tag{ 51 }$$

$$\begin{equation} S^{(2)}_{{\rm IIA}}(t_2,t_1) = C^{(2)} Q_2(t_2,t_1), \end{equation} \tag{ 52 }$$

$$\begin{equation} S^{(2)}_{{\rm IIB}}(t_2,t_1) = C^{(2)} Q^*_2(t_2,t_1), \end{equation} \tag{ 53 }$$

$$\begin{equation} S^{(2)}_{{\rm III}}(t_2,t_1) = C^{(2)} Q^*_1(t_2,t_1). \end{equation} \tag{ 54 }$$

The four Liouville space pathways are shown in figure 3. We started out with (n + 1)! terms, which were grouped into (n + 1) groups of terms that turned out to be identical, eliminating the factor of n! from the contour-ordered exponential. To have a specific time ordering, we split each term having k times on the lower branch into $\left({n}\atop{k}\right)$ individual terms, for a grand total of $\sum _{k=0}^n \left({n}\atop{k}\right) = 2^n$ terms, a result we could have anticipated since we distribute n times on two branches of the contour.

Zoom In Zoom Out Reset image size

Many coherent spectroscopy experiments performed today relate to the third-order response function, S⁽³⁾(t₃, t₂, t₁), and we shall briefly go through its derivation. Specifically, we shall focus on how to stack the decks when shuffling two time orderings, one on each branch of the Keldysh contour, into a single time ordering in physical time. We shall demonstrate, by explicit derivation, which Liouville space pathway is contributed by each loop diagram:

$$\begin{eqnarray} &&\fl \langle \hat{\mu }^{(3)}(t) \rangle = -\frac{1}{3!} \frac{1}{({\rm i}\hbar )^3} \int _c \! {\rm d}\tau _3 \int _c \! {\rm d}\tau _2 \int _c \! {\rm d}\tau _1 {\rm Tr}\lbrace {\rm T_c} [\hat{\mu }(\tau _3) \hat{\mu }(\tau _2)\nonumber\\ &&\times \hat{\mu }(\tau _1) \hat{\mu }(t)] \hat{\varrho }(t_0) \rbrace E(\tau _3) E(\tau _2) E(\tau _1) . \end{eqnarray} \tag{ 55 }$$

A full expansion of the contour ordering gives (3 + 1)! = 24 terms, of which we write four representative terms (that we will refer to as terms I–IV):

$$\begin{eqnarray} &&\fl {\rm T_c} [\hat{\mu }(\tau _3) \hat{\mu }(\tau _2) \hat{\mu }(\tau _1) \hat{\mu }(t)] = \vartheta (t - \tau _3) \vartheta (\tau _3 - \tau _2)\nonumber\\ &&\times \vartheta (\tau _2 - \tau _1) \hat{\mu }(t) \hat{\mu }(\tau _3) \hat{\mu }(\tau _2) \hat{\mu }(\tau _1) + \cdots \end{eqnarray} \tag{ 56 }$$

$$\begin{eqnarray} &&\fl {+} \vartheta (\tau _3 {-} t) \vartheta (t {-} \tau _2) \vartheta (\tau _2 {-} \tau _1)\ \hat{\mu }(\tau _3) \hat{\mu }(t) \hat{\mu }(\tau _2) \hat{\mu }(\tau _1) {+} \cdots\nonumber\\ \end{eqnarray} \tag{ 57 }$$

$$\begin{eqnarray} &&\fl {+} \vartheta (\tau _3 {-} \tau _2) \vartheta (\tau _2 {-} t) \vartheta (t {-} \tau _1)\ \hat{\mu }(\tau _3) \hat{\mu }(\tau _2)\hat{\mu }(t) \hat{\mu }(\tau _1) {+} \cdots \nonumber\\ \end{eqnarray} \tag{ 58 }$$

$$\begin{eqnarray} &&\fl + \vartheta (\tau _3 {-} \tau _2) \vartheta (\tau _2 {-} \tau _1) \vartheta (\tau _1 {-} t)\ \hat{\mu }(\tau _3) \hat{\mu }(\tau _2) \hat{\mu }(\tau _1)\hat{\mu }(t) {+} \cdots . \nonumber \\ \end{eqnarray} \tag{ 59 }$$

Consider term III, where one time occurs before t, and thus lies on the upper branch of the loop diagram. The remaining two times lie on the lower contour. Deforming to the Keldysh contour, we obtain

$$\begin{eqnarray} &&{-}\frac{1}{3!} \frac{1}{({\rm i}\hbar )^3} \int _{t}^{t_0} \! {\rm d}\tau _3 \int _{t}^{t_0} \! {\rm d}\tau _2 \int _{t_0}^t \! {\rm d}\tau _1 \vartheta (\tau _3 {-} \tau _2) \vartheta (\tau _2 {-} t) \vartheta (t {-} \tau _1) \nonumber \\ \end{eqnarray} \tag{ 60 }$$

$$\begin{equation} \times \ {\rm Tr}\lbrace \hat{\mu }(\tau _3) \hat{\mu }(\tau _2) \hat{\mu }(t) \hat{\mu }(\tau _1) \hat{\varrho }(t_0) \rbrace E(\tau _3) E(\tau _2) E(\tau _1) \end{equation} \tag{ 61 }$$

$$\begin{eqnarray} &&\fl = -\frac{1}{3!} \frac{1}{({\rm i}\hbar )^3} \int _{t_0}^t \! {\rm d}\tau _3 \int _{t_0}^t \! {\rm d}\tau _2 \int _{t_0}^t \! {\rm d}\tau _1 \vartheta (\tau _2 - \tau _3) \end{eqnarray} \tag{ 62 }$$

$$\begin{equation} \times \ {\rm Tr}\lbrace \hat{\mu }(\tau _3) \hat{\mu }(\tau _2) \hat{\mu }(t) \hat{\mu }(\tau _1) \hat{\varrho }(t_0) \rbrace E(\tau _3) E(\tau _2) E(\tau _1) . \end{equation} \tag{ 63 }$$

Note the sign change in the argument of the remaining Heaviside function. This is a consequence of contour ordering being deformed to anti-time ordering, and the subsequent time flip. We now have three times on the real time axis, but only two of these are ordered. To relate to physical time, full time ordering is desired. We split the term into three terms, where the time from the upper branch occurs before, in between, or after the two ordered times from the lower branch (we shall refer to these as IIIA, IIIB and IIIC):

$$\begin{eqnarray} &&\fl \vartheta (\tau _2 - \tau _3) = \vartheta (\tau _2 - \tau _3) \vartheta (\tau _3 - \tau _1) + \vartheta (\tau _2 - \tau _1)\vartheta (\tau _1 - \tau _3)\nonumber\\ &&+ \vartheta (\tau _1 - \tau _2) \vartheta (\tau _2 - \tau _3) . \end{eqnarray} \tag{ 64 }$$

As discussed above, a term with k times on the lower branch can be split into $\left({n}\atop{k}\right)$ terms with specific time ordering. Define the third-order response function indirectly by

$$\begin{eqnarray} &&\fl \langle \hat{\mu }^{(3)}(t) \rangle = \int \!{\rm d}t_3 \int \!{\rm d}t_2 \int \!{\rm d}t_1 S^{(3)}(t_3,t_2,t_1)\nonumber\\ &&\times E(t - t_3) E(t - t_3 - t_2) E(t - t_3 - t_2 - t_1). \end{eqnarray} \tag{ 65 }$$

Each term, for example IIIA, is accompanied by five identical terms, thus eliminating the factor $\frac{1}{3!}$. We substitute times to obtain

$$\begin{eqnarray} &&\fl S^{(3)}_{{\rm IIIA}}(t_3,t_2,t_1)\nonumber \\ &&= C^{(3)} {\rm Tr}\lbrace \hat{\mu }(t_1) \hat{\mu }(t_2 + t_1) \hat{\mu }(t_3 + t_2 + t_1) \hat{\mu }(0) \hat{\varrho }(-\infty ) \rbrace . \nonumber\\ \end{eqnarray} \tag{ 66 }$$

Defining four Liouville space pathways as

$$\begin{eqnarray} &&\fl R_1(t_3,t_2,t_1) {=} {\rm Tr}\lbrace \hat{\mu }(t_1) \hat{\mu }(t_2 {+} t_1) \hat{\mu }(t_3 {+} t_2 {+} t_1) \hat{\mu }(0) \hat{\varrho }(-\infty ) \rbrace , \nonumber\\ \end{eqnarray} \tag{ 67 }$$

$$\begin{eqnarray} &&\fl R_2(t_3,t_2,t_1) \hspace{-0.6pt}\,{=}\, \hspace{-0.6pt}{\rm Tr}\lbrace \hat{\mu }(0) \hat{\mu }(t_2\hspace{-0.6pt} \,{+}\,\hspace{-0.6pt} t_1) \hat{\mu }(t_3 \hspace{-0.6pt}\,{+}\,\hspace{-0.6pt} t_2 \hspace{-0.6pt}\,{+}\,\hspace{-0.6pt} t_1) \hat{\mu }(t_1) \hat{\varrho }(-\infty ) \rbrace , \nonumber\\ \end{eqnarray} \tag{ 68 }$$

$$\begin{eqnarray} &&\fl R_3(t_3,t_2,t_1)\hspace{-0.6pt}\, {=}\,\hspace{-0.6pt} {\rm Tr}\lbrace \hat{\mu }(0) \hat{\mu }(t_1) \hat{\mu }(t_3 \hspace{-0.6pt}\,{+}\,\hspace{-0.6pt} t_2 \hspace{-0.6pt}\,{+}\,\hspace{-0.6pt} t_1) \hat{\mu }(t_2 \hspace{-0.6pt}\,{+}\,\hspace{-0.6pt} t_1) \hat{\varrho }(-\infty ) \rbrace ,\nonumber\\ \end{eqnarray} \tag{ 69 }$$

$$\begin{eqnarray} &&\fl R_4(t_3,t_2,t_1)\hspace{-0.6pt}\, {=}\,\hspace{-0.6pt} {\rm Tr}\lbrace \hat{\mu }(t_3\hspace{-0.6pt} \,{+}\,\hspace{-0.6pt} t_2 \hspace{-0.6pt}\,{+}\,\hspace{-0.6pt} t_1) \hat{\mu }(t_2 \hspace{-0.6pt}\,{+}\,\hspace{-0.6pt} t_1) \hat{\mu }(t_1) \hat{\mu }(0) \hat{\varrho }(-\infty ) \rbrace , \nonumber\\ \end{eqnarray} \tag{ 70 }$$

the response function is given by

$$\begin{equation} \fl S^{(3)}(t_3, t_2, t_1) = C^{(3)} \sum _{i = 1}^{4}[R_i(t_3,t_2,t_1) - R_i^*(t_3,t_2,t_1) ]. \end{equation} \tag{ 71 }$$

The contribution from the various terms are as follows:

$$\begin{equation} S^{(3)}_{{\rm IIIA}}(t_3,t_2,t_1) = C^{(3)} R_1(t_3,t_2,t_1) \end{equation} \tag{ 72 }$$

$$\begin{equation} S^{(3)}_{{\rm IIIB}}(t_3,t_2,t_1) = C^{(3)} R_2(t_3,t_2,t_1) \end{equation} \tag{ 73 }$$

$$\begin{equation} S^{(3)}_{{\rm IIIC}}(t_3,t_2,t_1) = C^{(3)} R_3(t_3,t_2,t_1) \end{equation} \tag{ 74 }$$

$$\begin{equation} S^{(3)}_{\rm I}(t_3,t_2,t_1) = C^{(3)} R_4(t_3,t_2,t_1) . \end{equation} \tag{ 75 }$$

The four remaining terms are related to the first by complex conjugation:

$$\begin{equation} S^{(3)}_{{\rm IIA}}(t_3,t_2,t_1) = - C^{(3)} R^*_1(t_3,t_2,t_1), \end{equation} \tag{ 76 }$$

$$\begin{equation} S^{(3)}_{{\rm IIB}}(t_3,t_2,t_1) = - C^{(3)} R^*_2(t_3,t_2,t_1), \end{equation} \tag{ 77 }$$

$$\begin{equation} S^{(3)}_{{\rm IIC}}(t_3,t_2,t_1) = - C^{(3)} R^*_3(t_3,t_2,t_1), \end{equation} \tag{ 78 }$$

$$\begin{equation} S^{(3)}_{{\rm IV}}(t_3,t_2,t_1) = - C^{(3)} R^*_4(t_3,t_2,t_1) . \end{equation} \tag{ 79 }$$

The eight Liouville space pathways contributing to the third-order response function are depicted in figure 4. With the strategy laid out in these sections, one could now derive expressions for the 2ⁿ Liouville space pathways contributing to the nth-order response function, and subsequently evaluate the corresponding correction to the dipole moment operator:

$$\begin{equation} \langle \hat{\mu }^{(n)}(t) \rangle \hspace{-0.6pt}=\hspace{-0.6pt} \prod _{i=1}^n \int \! {\rm d}t_i S^{(n)}(t_n,\ldots ,t_1) \prod _{i=1}^n E\left(t\hspace{-0.6pt} -\hspace{-0.6pt} \sum _{j=i}^n t_j\right)\hspace{-0.6pt} . \ \end{equation} \tag{ 80 }$$

Zoom In Zoom Out Reset image size

As the avid reader has noted, the loop diagrams relate to the double-sided Feynman diagrams encountered in nonlinear optical spectroscopy. Imagine the loop diagram in figure 5 to be rotated at 90°. Note that the information about the phases of the electrical field implied by the corresponding arrows are not given in the loop diagram.

Zoom In Zoom Out Reset image size

As a prelude, we shall consider transitions in a purely electronic system and obtain Fermi's golden rule as a linear response function. Consider a two-level system, H₀ = ε_a|a〉〈a| + ε_b|b〉〈b|, initially prepared in the state |a〉, $\hat{\varrho }(-\infty ) = |a \rangle \langle a|$. At time t₀, a coupling between the states, V = V_ab|a〉〈b| + V_ba|b〉〈a|, is introduced and we monitor the rate of transition $\hat{k}$, calculated using the Liouville equation:

$$\begin{equation} \hat{k} = -\frac{{\rm i}}{\hbar } [H, |a \rangle \langle a|] \end{equation} \tag{ 81 }$$

$$\begin{equation} = -\frac{{\rm i}}{\hbar } (V_{ba} |b \rangle \langle a| - V_{ab} |a \rangle \langle b|). \end{equation} \tag{ 82 }$$

Now, evaluate the transition rate to first order in the coupling (using expression (17)):

$$\begin{eqnarray} &&\fl \langle \hat{k}^{(1)} \rangle = \frac{1}{{\rm i}\hbar } \int \! {\rm d}t_1 \vartheta (t_1) [ {\rm Tr}\lbrace \hat{k}(t_1)\hat{V}(0) \hat{\varrho }(-\infty ) \rbrace \nonumber\\ &&\quad\;\; - {\rm Tr}\lbrace \hat{V}(0) \hat{k}(t_1) \hat{\varrho }(-\infty ) \rbrace ] \end{eqnarray} \tag{ 83 }$$

$$\begin{eqnarray} &&\quad\! = \frac{1}{{\rm i}\hbar } \int \! {\rm d}t_1 \vartheta (t_1) [ \langle a | \hat{k}(t_1)\hat{V}(0) | a \rangle - \langle a | \hat{V}(0)\hat{k}(t_1) | a \rangle ] \nonumber\\ \end{eqnarray} \tag{ 84 }$$

$$\begin{eqnarray} &&= \frac{1}{{\rm i}\hbar } \int \! {\rm d}t_1 \vartheta (t_1) [ V_{ba} \langle a | {\rm e}^{{\rm i}H_0 t_1} \hat{k} {\rm e}^{-{\rm i}H_0 t_1} | b \rangle \nonumber\\ &&- V_{ab} \langle b | {\rm e}^{{\rm i}H_0 t_1} \hat{k} {\rm e}^{-{\rm i}H_0 t_1} | a \rangle ] \hspace{42pt} \end{eqnarray} \tag{ 85 }$$

$$\begin{equation} = \frac{1}{{\rm i}\hbar ^2} \int \! {\rm d}t_1 \vartheta (t_1) |V_{ab}|^2 [ {\rm i} {\rm e}^{{\rm i} (\varepsilon _a - \varepsilon _b) t_1} + {\rm i} {\rm e}^{{\rm i} (\varepsilon _b - \varepsilon _a) t_1} ] \end{equation} \tag{ 86 }$$

$$\begin{equation} = \frac{{\rm i}}{\hbar } |V_{ab}|^2 \left[ \frac{1}{\varepsilon _a - \varepsilon _b +{\rm i}\eta } + \frac{1}{\varepsilon _b - \varepsilon _a +{\rm i}\eta }\right] \end{equation} \tag{ 87 }$$

$$\begin{equation} = \frac{2}{\hbar } |V_{ab}|^2 \frac{\eta }{(\varepsilon _b - \varepsilon _a)^2 + \eta ^2} \end{equation} \tag{ 88 }$$

$$\begin{equation} = \frac{2\pi }{\hbar } |V_{ab}|^2 \delta (\varepsilon _b - \varepsilon _a) . \end{equation} \tag{ 89 }$$

We recognize Fermi's golden rule and note that the rate is independent of time t, since the perturbation is not associated with an external field.

Mathematically, the derivation of the Marcus equation is a slight generalization of the derivation of Fermi's golden rule. Considering the two-level system coupled to a single bosonic mode,

$$\begin{eqnarray} &&\fl H = [\varepsilon _a + v_a(c^{\dagger } + c) ] |a \rangle \langle a| + [\varepsilon _b + v_b(c^{\dagger } + c) ] |b \rangle \langle b| \nonumber \\ &&+ V_{ab} |a \rangle \langle b| + V_{ba} |b \rangle \langle a| + \hbar \omega c^{\dagger }c . \end{eqnarray} \tag{ 90 }$$

Shift the bosonic mode ($a = c + \frac{v_a}{\hbar \omega }$ and $a^{\dagger } = c^{\dagger } + \frac{v_a}{\hbar \omega }$) to obtain

$$\begin{eqnarray} &&\fl H = \left[\varepsilon _a - \frac{v^2_a}{\hbar \omega } \right] |a \rangle \langle a| + \left[\varepsilon _b - \frac{v^2_b}{\hbar \omega } + \frac{(v_b - v_a)^2}{\hbar \omega } \right. \nonumber \\ &&\!\!\!\!\!\!\!\!\!\left. + (v_b - v_a)(a^{\dagger } + a)\right] |b \rangle \langle b| + V_{ab} |a \rangle \langle b| + V_{ba} |b \rangle \langle a| + \hbar \omega a^{\dagger }a \nonumber \\ \end{eqnarray} \tag{ 91 }$$

$$\begin{eqnarray} &&= \tilde{\varepsilon }_a |a \rangle \langle a| + \left[\tilde{\varepsilon }_b + \frac{\Delta v^2}{\hbar \omega } + \Delta v (a^{\dagger } + a)\right] |b \rangle \langle b| \nonumber \\ &&+ V_{ab} |a \rangle \langle b| + V_{ba} |b \rangle \langle a| + \hbar \omega a^{\dagger }a \hspace{48pt} \end{eqnarray} \tag{ 92 }$$

$$\begin{equation} = \tilde{\varepsilon }_a |a \rangle \langle a| + \tilde{\tilde{\varepsilon }}_b |b \rangle \langle b| + V_{ab} |a \rangle \langle b| + V_{ba} |b \rangle \langle a| + \hbar \omega a^{\dagger }a , \end{equation} \tag{ 93 }$$

where we have defined $\tilde{\varepsilon }_a = \varepsilon _a - \frac{v^2_a}{\hbar \omega }$, Δv = v_b − v_a and $\tilde{\tilde{\varepsilon }}_b = \tilde{\varepsilon }_b + \frac{\Delta v^2}{\hbar \omega } + \Delta v (a^{\dagger } + a)$. We use the transition operator of equation (82) and choose $\hat{\varrho }(-\infty ) = |a \rangle \langle a| \hat{\varrho }^{{\rm bath}}$ with $\hat{\varrho }^{{\rm bath}} = {\rm e}^{-\beta \hbar \omega a^{\dagger }a} (1 - {\rm e}^{-\beta \hbar \omega })$. The transition rate to linear order in the coupling is given by

$$\begin{eqnarray} \fl \langle \hat{k}^{(1)} \rangle &=& \frac{1}{{\rm i}\hbar } \int {\rm d}t_1 \vartheta (t_1) [{\rm Tr}\lbrace \hat{k}(t_1) \hat{V}(0) \hat{\varrho }(-\infty )\rbrace \nonumber\\ &&- {\rm Tr}\lbrace \hat{V}(0) \hat{k}(t_1) \hat{\varrho }(-\infty )\rbrace ] \end{eqnarray} \tag{ 94 }$$

$$\begin{eqnarray} \fl\hphantom{\langle \hat{k}^{(1)} \rangle} &=& \frac{1}{{\rm i}\hbar } \int \! {\rm d}t_1 \vartheta (t_1) [ {\rm Tr_{bath}}\lbrace \langle a | \hat{k}(t_1)\hat{V}(0) | a \rangle \hat{\varrho }^{{\rm bath}} \rbrace \nonumber\\ &&- {\rm Tr_{bath}}\lbrace \langle a | \hat{V}(0)\hat{k}(t_1) | a \rangle \hat{\varrho }^{{\rm bath}} \rbrace ] \end{eqnarray} \tag{ 95 }$$

$$\begin{eqnarray} \fl\hphantom{\langle \hat{k}^{(1)} \rangle}&=& \frac{1}{{\rm i}\hbar } \int \! {\rm d}t_1 \vartheta (t_1) [ V_{ba} {\rm Tr_{bath}}\lbrace \langle a | {\rm e}^{{\rm i}H_0 t_1} \hat{k} {\rm e}^{-{\rm i}H_0 t_1} | b \rangle \hat{\varrho }^{{\rm bath}} \rbrace \nonumber\\ &&- V_{ab}{\rm Tr_{bath}}\lbrace \langle b | {\rm e}^{{\rm i}H_0 t_1} \hat{k} {\rm e}^{-{\rm i}H_0 t_1} | a \rangle \hat{\varrho }^{{\rm bath}}\rbrace ] \end{eqnarray} \tag{ 96 }$$

$$\begin{eqnarray} \fl \hphantom{\langle \hat{k}^{(1)} \rangle}&=& \frac{|V_{ab}|^2}{\hbar ^2} \int \! {\rm d}t_1 \vartheta (t_1) \langle {\rm e}^{{\rm i}(\tilde{\varepsilon }_a + \hbar \omega a^{\dagger }a )t_1} {\rm e}^{-{\rm i}(\tilde{\tilde{\varepsilon }}_b + \hbar \omega a^{\dagger }a )t_1} \nonumber\\ &&+ {\rm e}^{{\rm i}(\tilde{\tilde{\varepsilon }}_b + \hbar \omega a^{\dagger }a )t_1} {\rm e}^{-{\rm i}(\tilde{\varepsilon }_a + \hbar \omega a^{\dagger }a )t_1} \rangle \end{eqnarray} \tag{ 97 }$$

$$\begin{eqnarray} &&\fl\hphantom{\langle \hat{k}^{(1)} \rangle} = \frac{|V_{ab}|^2}{\hbar ^2} \int \! {\rm d}t_1 \langle {\rm e}^{{\rm i}(\tilde{\tilde{\varepsilon }}_b + \hbar \omega a^{\dagger }a )t_1} {\rm e}^{-{\rm i}(\tilde{\varepsilon }_a + \hbar \omega a^{\dagger }a )t_1} \rangle . \end{eqnarray} \tag{ 98 }$$

With the merging of the two terms, we are back on the beaten track and the final steps are found in multiple textbooks [15, 6]. The integrand factorizes as

$$\begin{equation} \langle {\rm e}^{{\rm i}(\tilde{\tilde{\varepsilon }}_b + \hbar \omega a^{\dagger }a )t_1} {\rm e}^{-{\rm i}(\tilde{\varepsilon }_a + \hbar \omega a^{\dagger }a )t_1} \rangle = {\rm e}^{{\rm i}(\tilde{\varepsilon }_b - \tilde{\varepsilon }_a)t_1} C_{ba}(t_1) . \end{equation} \tag{ 99 }$$

The correlation function C_ba(t₁) evaluates to

$$\begin{equation} C_{ba}(t_1) = {\rm e}^{- \frac{\Delta v^2}{(\hbar \omega )^2}(2n(\omega )+1) + \frac{\Delta v^2}{(\hbar \omega )^2} ((n(\omega ) +1){\rm e}^{-{\rm i}\omega t_1} + n(\omega ) {\rm e}^{{\rm i}\omega t_1})}, \end{equation} \tag{ 100 }$$

where n(ω) = (e^βℏω − 1)⁻¹ is the Bose distribution function. In the limit where the correlation function is short-lived, $\left( \frac{\Delta v}{\hbar \omega } \right)^2 n(\omega ) \gg \ 1$, the exponential functions in the exponent can be expanded to second order in ωt₁ and the integral becomes a Gaussian integral. We obtain the celebrated Marcus equation

$$\begin{equation} \langle \hat{k}^{(1)} \rangle = \frac{|V_{ab}|^2}{\hbar ^2} \sqrt{\frac{\pi }{a}} {\rm e}^{- \frac{(\tilde{\varepsilon }_b - \tilde{\varepsilon }_a - \frac{\Delta v^2}{\hbar \omega } )^2}{4a\hbar ^2}}, \end{equation} \tag{ 101 }$$

where $a = \frac{1}{2}(2n(\omega ) +1)\frac{\Delta v^2}{\hbar ^2}$. The Keldysh contour derivation given here turned out to be no more than a response function derivation. The methodology used, however, applies also to higher order corrections.