Study of resonant inelastic light scattering in Keldysh-Schwinger functional integral formalism

The scattering cross section of the resonant inelastic light scattering is represented as a correlation function in the Keldysh-Schwinger functional integral formalism. The functional integral approach enables us to compute the cross section in the Feynman diagram perturbation theory where many-body effects can be fully incorporated. This approach is applied to the one G-phonon Raman scattering of graphene, and the result is shown to agree with the one previously obtained by the conventional Fermi golden rule formula. Also, this approach is generalized to the systems in non-equilibrium conditions.


Introduction
The inelastic light scattering (ILS) is a very important experimental method for the investigation of the physical properties of condensed matter. In particular, the resonant inelastic X-ray scattering (often called RIXS) has developed very rapidly recently, mainly owing to the advent of the high power X-ray light sources [1]. RIXS has proven to be a versatile tool for the study of elementary excitations of strongly correlated electron materials, such as magnons [2] and orbitons, where many-body interactions play crucial roles [1]. Thus it is very desirable to have a theoretical scheme for the interpretation of the experimental ILS data which can incorporate many-body effects precisely.
In general, the transition rate w ILS for ILS can be obtained by Fermi golden rule which is applied up to the second order with respect to the interaction Hamiltonian H ′ between light and matter: where |I , |n , and |F denote the initial states, the (many-body) intermediate eigenstates, and the final states of the light-matter system, respectively. E I,n are the energy eigenvalues of the states and p I is the probability for the initial states ( I p I = 1).Ĥ ′ NR is the part of the HamiltonianĤ ′ quadratic in the photon vector potential operator A, which contributes to the nonresonant ILS [3].Ĥ ′ R is linear in A, and it is responsible for the resonant ILS [3] [the second term of (1)]. The nonresonant transition amplitude often dominates over the resonant one, but in special cases where E I is close to E n , the resonant transition amplitude can be significant. After some manipulations involving the expectation value of the photon operators, (1) can be shown to be essentially the Kramers-Heisenberg formula of ILS. In this paper, for simplicity, we will consider only the caseĤ ′ R = − 1 c d r A( r) · J e ( r), where J e ( r) is the electric current operator.
It is well known that the scattering cross-section of nonresonant contribution can be expressed in the form of correlation function of stress tensor operators [3,4] by using van Hove-Placzek procedure [5,6]. By employing the fluctuation-dissipation theorem the correlation function can be directly linked to the retarded Green's function (e.g., see p.151 of [4]). Then the retarded Green's function can be computed systematically (including full many-body effects) in imaginary time Matsubara formalism [4,7] by using analytic continuation technique. Thus as far as the nonresonant contribution to ILS is concerned we have a firmly established theoretical scheme where full manybody physics can be incorporated precisely. Henceforth we will assume that resonance condition is satisfied and will focus on the resonant contributions only.
The first goal of this paper is to express the resonant contribution to ILS in the form of correlation function just like the case of the nonresonant contribution mentioned above. To achieve this goal, in principle, we have to carry out the sum over the intermediate states |n of the second term of (1). In general, very little is known about the states |n , and we have to employ some kinds of assumptions or approximations about the states. Correlation functions are defined without reference to intermediate states, so in our approach the explicit sum over the intermediate states |n should be avoided. Then we have to go back to the time-dependent perturbation theory involving time evolution operator. Recall that the time-dependent perturbation theory requires time-ordering of operators [8], namely T(X(t)Ŷ (t ′ )) = θ(t − t ′ )X(t)Ŷ (t ′ ) ± θ(t ′ − t)Ŷ (t ′ )X(t), (2) where minus sign applies when bothX(t) andŶ (t ′ ) are fermion operators and T is the time-ordering symbol [9]. θ(t) is the unit step function. Now as can be seen in (1), we have to take the complex conjugate of transition amplitude to obtain the transition rate. In quantum mechanics the complex conjugate of a transition amplitude of a process is related to a transition amplitude of the time-reversed process, and this observation necessitates the introduction of anti-timeordering:T whereT is the anti-time-ordering symbol. From the above discussion it is evident that the scattering cross section of ILS (generically all transition probability) requires not only forward time evolution but also backward time evolution in the time representation of the time-dependent perturbation theory. Therefore, we are naturally led to the concept of closed time contour introduced by Keldysh and Schwinger (KS) [10,11].
The KS formulation of time-dependent perturbation theory has been employed in the earlier works on the X-ray spectra of metals [12,13,14]. However, these works did not fully disentangle the time orderings of involved operators, so that they were limited to the time-representation of Feynman diagrams or the time-ordered Feynman diagrams [see, for example, the figure 2 of [13] and the equation (11) of [14]]. If the time orderings are fully disentangled then all time-integrals can be done from −∞ to ∞, and all Feynman diagrams can be calculated in energy-momentum representation.
The second goal of this paper is to elaborate on the formalism developed in the references [12,13,14], so that the correlation function obtained in the first goal can be recast in the KS functional integral formalism, which will disentangle the time orderings completely. The correlation functions in KS functional integral formalism can be computed in the standard Feynman diagram method of quantum field theory [15] without being limited to time-ordered diagrams, etc.
The third goal of this paper is to generalize the KS functional integral approach to the systems in non-equilibrium conditions. The KS methods were originally developed for the study of non-equilibrium phenomena. As such, the KS functional integral formalism can be readily generalized for the study of ILS in non-equilibrium situation.
The main results of this paper are (21), (23), (47), (74), and (75). This paper is organized as follows. In section 2 we derive the correlation function corresponding to the scattering cross-section of the resonant ILS. In section 3 the correlation function is represented in KS functional integral formalism. The key step is the disentanglement of the time orderings of electric current operators. In section 4 we compare our approach with the Kramers-Heisenberg formula and point out differences. In section 5 we apply our method to the case of G-phonon Raman scattering of graphene and compare the results with the one obtained by conventional Fermi golden rule method. In section 6 we generalize the KS functional integral formalism of ILS to the non-equilibrium situations. We close this paper with discussions and summary in section 7.

Scattering cross-section of ILS and correlation function
First we recall that the scattering cross-section of ILS is defined by the transition rate (induced by light-matter interaction) divided by the flux of incoming photons: where M FI is the transition amplitude from an initial state |I to a final state |F in time interval [t I , t F ] and V is the volume of space for the box quantization of photons. Our goal in this section is to express (4) in the form of correlation function. The key quantity is the transition amplitude M FI (t F , t I ) which is most naturally presented in the framework of time-dependent perturbation theory [8].
The transition amplitude M FI (t F , t I ) is defined by where U Htot (t F , t I ) is the time evolution operator governed by the total Hamiltonian H tot of light-matter system from time t I to t F : The Next we employ the interaction picture and treat H ′ perturbatively. In the interaction picture the time evolution operator U Htot (t F , t I ) becomes (see, for example, pp.722-724 of [16]) where UĤ 0 (t F , t I ) is the time evolution operator governed by the Hamiltonian H 0 [ see (7)], and where H ′ H0 denotes the Hamiltonian H ′ in the Heisenberg picture with respect to H 0 . Now plugging (8) into (5) and expanding (9) up to the second order in H ′ H0 , we obtain the second order contribution which is responsible for the resonant ILS: Then the transition probability can be expressed as Note the appearance of the anti-time orderingT in (11).
The next key step is to compute the matrix elements of (11) with respect to the photon states (only). For that purpose, the initial |I and the final |F states are represented by the tensor product of photon Hilbert space |k i/f ,ǫ i/f and matter Hilbert space |i/f : where k i,f andǫ i,f are the wavenumber vector and the polarization vector of incoming and outgoing photon, respectively. Since the matter Hamiltonian H m commutes with the photon Hamiltonian H p which is non-interacting [see (A.1)], the action of H p and its time evolution operator U Hp in (11) on photon Hilbert space can be evaluated exactly. After evaluating the action of U Hp , we can define an operatorDk i ,k f and its Hermitian conjugateD † ki,k f which act only on the matter Hilbert space.
wherek i,f denotes k i,f andǫ i,f collectively. Then the transition probability (11) can be recast as We note that the operator Dk The initial photon state |k i ,ǫ i can be taken to a pure state, so that the probability distribution for the initial state p I pertains only to the initial matter Hilbert space. The probability distribution p i for the initial matter Hilbert space can be specified by a density matrixρ m ( Tr m indicates trace over matter Hilbert space): The scattering cross-section of resonant ILS now takes the form The sum over the final matter Hilbert space of (17) can be done using the closure relation of matter Hilbert space f |f f| = I ( I is an identity operator). From (15) we obtain The sum over the initial matter Hilbert space |i of (17) can done by using (16).
Equation ( Substituting (A.6) and (A.7) into (19) we arrive at where the correlation function For matter systems at equilibrium the correlation function C ki,k f (t 3 , t 4 , t 1 , t 2 ) possesses the time translation invariance which is manifest in the following form: Plugging (22) into (20), changing variable t i → t i + (t F + t I )/2, and carrying out the time integrals in the limit where the energy conserving delta function is substituted by (21) and (23) constitute the first goal of this paper.

Keldysh-Schwinger functional integral formulation of correlation function
As we have seen in section 2 the computation of scattering cross-section has been reduced to the computation of the correlation function (21), which we rephrase as for the sake of more convenient notations. From now on we will drop the subscript "m" denoting the matter sector. In this section the correlation function (25) will be represented in KS functional integral formulation, and in this representation all time-orderings will be disentangled completely. The literature on the method of KS closed time contour are very vast, and we cite only [17,18] which appear to be more useful for condensed matter physicists in author's opinion. We will also need the method of coherent state path integral, and the detailed treatments can be found in [19]. In the following many well-known detailed steps (for which readers are referred to the references cited above) will be omitted and we will focus on the disentangling the time orderings of (25).
The key concept of the KS method is the time evolution along closed time contour. A characteristic property of this closed time contour which follows from the unitarity of time evolution operator U. We start by resolving the time orderings of (25) explicitly: where the correlation functions Π (1,2,3,4) ({τ j }) are given by Next we represent each correlation function Π (1,2,3,4) in Heisenberg picture can be expressed in terms of the time evolution operator. For A(τ i ), and similarly for other operators. Substituting (30) into Π (1) of (29) and employing the composition property of time evolution operator U(t 1 , where the turning point of the closed time contour is marked with a vertical bar. From now on we will take a limit t F → ∞, t I → −∞. The time evolution to the left of the bar of (31) is backward in time, while that to the right of the bar is forward in time.
Thus in (31) we have time evolution along the closed time contour, and the operators are inserted along the contour.
To develop the functional integral representation, the time evolution operators of (31) are to be expressed by the Trotter product formula [19] in accordance with the closed time contour of figure 1. For t > t ′ ( the time evolution in forward branch) And for t < t ′ (the time evolution in backward branch) we note Then between each infinitesimal time evolution operator e ±i∆t H/ , we insert the following closure relation of coherent states (I is an identity operator): where {φ ±,k } and {φ * ±,k } denote the collective coherent state variables at the k-th time slice along the forward (+) and the backward (-) branch. Recall that the ket coherent state is an eigenstate of annihilation operator, while the bra coherent state is an eigenstate of creation operator: This property makes the computation of matrix elements of any normal-ordered operators very easy: Note that the symbol A in the right hand side of (36) is not an operator but a function of coherent state variables φ * , φ ′ . All operators appearing in correlation function (31) are assumed to be normal-ordered. Employing (34), (35), and (36) it is straightforward to develop functional integral representation of the correlation function (31). Below we exhibit a few details of the computations near operator insertion points. Consider, for example, the insertion of the operator C in (31). Clearly the time slicing can be done in such a way that the set of time slices include the time slice at τ 3 (or arbitrarily close to it). Let |φ +,k be the ket of the time slice corresponding to the time τ 3 (recall that the operatorĈ is located in the forward branch). Then the relevant matrix elements are Since the time interval ∆t is infinitesimally small, we can approximate e −i∆t H/ ≈ 1 − i∆t H/ and obtain Inserting a closure relation (34) into (38), (38) becomes The operator C =Ĉ(ĉ † ,ĉ) and the Hamiltonian H can be assumed to be general polynomials of the creationĉ † and the annihilation operatorĉ. Now ξ, ξ * integrals of (39) can be done exactly. This is because the operator C and H are assumed to be polynomials ofĉ andĉ † , so that the most general integrals involving ξ and ξ * are of the following form: (n, m are non-negative integers) From (40), it follows that (since polynomial form is preserved through integral) In the continuum limit N → ∞, the difference between φ * k+1 and φ * k in C(φ * +,k+1 , φ +,k ) and H(φ * +,k+1 , φ +,k ) can be ignored. Then (41) proves that the operator C(τ ) which is inserted in the forward branch maps to C(φ * + (τ ), φ + (τ )) in KS coherent functional integral formulation: In an entirely similar way it can be shown that Combining the above results we arrive at the following KS functional integral representation of the correlation function (31): where S ± is the well-known classical action of functional integral defined on the forward (+) and backward(-) branch. δS +− is infinitesimally small terms which couple the forward and the backward branch functional integral variables. For the subtle roles played by δS +− , see [17] and [20]. We emphasize that the A, B, C, D of (44) are functions of ordinary variables (or anti-commuting Grassman numbers for fermion operators) and that they are not operators, so that the relative order between them does not matter. Now recalling the definitions of the correlation functions Π (1,2,3,4) ({τ i }) from (29), we find that all of Π (1,2,3,4) ({τ i }) are given by the same expression, namely (44), in the KS functional integral. But they are defined in the separate regions of time variables {τ i } as signified by the products of step functions in (28). Noting a trivial identity Applying the result (46) to the correlation function of ILS (21) we arrive at which is the second main result of this paper. (47) is the four-current correlation function which can be computed directly by using the standard functional integral methods such as Feynman diagram perturbation theory or semi-classical approximations.

Comparison with Kramers-Heisenberg formula
Let us compare the results (21) and (23) of this paper with Kramers-Heisenberg formula which is essentially Fermi golden rule (1). At this point we emphasize that in deriving (21) and (23) no approximations have been made regarding any intermediate states.
Below we will point out that Kramers-Heisenberg formula was obtained by ignoring certain terms which are attributed to the sudden turning on of perturbation [8]. This implies that our results (21) and (23) contain more contributions than Kramers-Heisenberg formula, so that we have to select a subset of Feynman diagrams contributing to (21) and (23) which indeed correspond to Kramers-Heisenberg formula.
We start by evaluating the time integrals of (10) by introducing intermediate states |n . For simplicity put t I = 0, and resolve time-ordering explicitly.
In the limit t F − t I → ∞, only the first term in the parenthesis of (48) yields a contribution which is proportional to t F − t I and it implements the energy conservation δ(E F − E I ). This contribution leads to the Fermi golden rule (1) and to the Kramers-Heisenberg formula. In the same limit the second term does not give substantial contribution unless E F ∼ E n , and this constraint greatly suppresses the number of possible intermediate states |n , resulting in negligible contribution compared to the first one [8]. The constraint E F ∼ E n can be enforced by introducing δ(E F − E n ) in the sum over intermediate states. This delta function can be naturally understood as a spectral function which is the imaginary part of retarded Green's function. Therefore, in our formalism, we have to ignore the Feynman diagrams which have more spectral functions than those which correspond to Kramers-Heisenberg formula. This assertion will be explicitly demonstrated in section 5. In transition probability, the ignored terms (which is to be squared) will allow processes which are not constrained by the energy conservation E F = E I . These processes will turn out to be strongly suppressed by more factors of spectral function. An explicit example will be given in section 5. So this is the price we have to pay: we have obtained a correlation function which does not depend on intermediate states, but we have to select appropriate Feynman diagrams. Fortunately, the selection of diagrams will turn out to be more or less obvious and straightforward.
We also note that the ignored terms can be eliminated by adiabatic turning on of perturbations [8], but in our formalism this method turns out to be very cumbersome, so it is not pursued further.

Application: single G-phonon Raman intensity of graphene
In this section we apply KS functional integral formalism to the case of the single G-phonon Raman intensity of graphene and compare the results with those obtained by conventional approach [22] based on Fermi golden rule.
For explicit calculations we need more elaborations in KS formalism. More details can be found in references [17,20]. For fermions let us use the notation (ψ ± , ψ * ± ) for coherent state variables instead of (φ ± , φ * ± ), and reserve the notation (φ ± , φ * ± ) for bosons. It turns out that in practical calculations it is more convenient to take a linear combination of coherent state variables (often called Keldysh basis) in the following way [17,20]. For bosons, the Keldysh basis is given by ( "cl"=classical, "q"= quantum, see [17] ) For fermions, a different combination is chosen.
Then Green's functions of bosons (such as phonon) are given by (a, b = cl, q ) where the superscripts "K,R,A" indicate Keldysh, retarded, and advanced Green's functions, respectively. On the other hand. Green's functions of fermions are given by (a, b = 1, 2) At equilibrium, the Keldysh Green's function is determined by retarded and advanced Green's functions owing to the fluctuation-dissipation theorem [17,18].
Since the retarded and the advanced Green's functions in frequency space are complex conjugate of each other, (53) implies that the Keldysh Green's function is purely imaginary (or anti-Hermitian). We recall the difference between the retarded and the advanced Green's function in energy space is 2i times the spectral function, which means that in equilibrium the Keldysh Green's function is always proportional to the spectral function. Now we turn to the case of graphene [for a review, see [21]]. The graphene has honeycomb lattice structure with two sublattices [denoted by A, B in figure 2 ], so that all the electron Green's functions become 2x2 matrix in sublattice space (the spin will be neglected for simplicity). More concretely, where t k = −t 0 i=1,2,3 e −ik·di with d i being vectors connecting nearest neighbors and t 0 is a hopping amplitude. I 2 is the identity matrix in sublattice space. The retarded Green's function of G-phonon ( u in figure 2) is given by (M is the mass of carbon atom) where ω q is the G-phonon dispersion. The Einstein phonon approximation ω q ≈ ω ph will be assumed. The electric currents of (47) of graphene [22] can be found in terms of coherent state variables using (50): where v k = ∂t k ∂k , Ψ = (ψ 1A , ψ 1B , ψ 2A , ψ 2B ) t , and For the Raman scattering of graphene, the photon momentum can be neglected, so that q ≈ 0. Finally, the (simplified) interaction between electron and G-phonon in action form is given by (N is the number of lattice sites) where γ cl = I 2 (identity matrix) and γ q = σ x (Pauli matrix) act on the Keldysh ] is the phonon coordinate in Keldysh basis and λ is the electron-phonon coupling constant. The detailed forms of functions v k and f k,q do not concern us in this paper. Now we are ready to compute the single G-phonon Raman intensity graphene in KS functional integral formalism. The single phonon Raman intensity in the lowest order perturbation theory is described by the Feynman diagram of figure 3. The direct application of Feynman rules yields the following result for the correlation function (47):C where (M a −− is for the loop in the left side of Fig. 3, and M b ++ is for the loop in the right) The trace "Tr" is over both Keldysh and sublattice spaces. The evaluation of the correlation function proceeds as follows: (1) take the trace over Keldysh space (2) then using the (anti) Hermiticity properties of matrices of Green's functions it is easily shown that M cl −− = (M cl ++ ) * and M q −− = −(M q ++ ) * (3) sort out all terms with more than two Keldysh Green's function G K or spectral functions; these terms belong to the ignored terms discussed in section 4. Direct integration indeed shows that there is no resonant behaviour (roughly two delta functions eliminate all singular behaviours). (4) after the sorting, a relation M cl −− = −M q −− can be shown by using the Hermiticity of matrices and the property dǫG R (ǫ)G R (ǫ ± ω i/f )G R (ǫ ± ω f/i ) = 0 (also for the advanced Green's function) which follows from the analyticity of the retarded Green's functions (5) finally using the property [ −− can be shown to be purely imaginary. Then the correlation function boils down tõ where (53) is employed in the second line and the explicit form of phonon Greens' function (55) is used in the third line. The delta function δ(ω i − ω f − ω ph ) of (62) corresponds to the Stokes process of phonon emission which dominates over the anti-Stokes process of phonon absorption in the low temperature limit due to the thermal occupation number factor coth[(ω i − ω f )/2k B T ] + 1. Now the computation of M cl −− is rather straightforward (but tedious). Substituting (62) into (23), the integral over the final photon momentum k f is done by the delta function of (62), leaving behind an integral over the solid angle of k f , which defines a differential cross section. The result obtained in this way is identical with the equations (7), (12a) of [22] [F 0 α 3 part of (12b,c,d)]. More precisely, [22] also includes additional phonon-photon coupling specific to graphene which is irrelevant in our discussion. The Feynman diagram of figure 4 depicts a process for the correlation function (47) which is allowed by Feynman rule but does not satisfy the correct energy conservation, as can be seen clearly in the phonon energy ω i + ω f . This Feynman diagram corresponds to the squared ignored terms discussed in section 4. Direct computation shows that all terms contributing to this Feynman diagrams contain at least two factors of spectral functions or Keldysh Green function, so they are strongly suppressed and do not contribute to the resonant scattering. A simple criterion for this type of diagrams is that they have E + and E − vertices in the same loop (compare with figure 3).

Inelastic light scattering for systems in nonequilibrium conditions
KS functional integral approach for ILS can be readily generalized to the systems in non-equilibrium conditions. In fact, the study of non-equilibrium phenomena was the original motivation for the introduction of KS method. We will assume that the nonequilibrium conditions are created by certain driving forces for the matter systems. Specifically we will consider the semiconductors of direct transition type such as GaAs, where the driving forces are pump pulse laser field which can be treated classically. This pump pulse field prepares the initial matter state |i at time t I for ILS, and then the probe photon of ILS is incident [23].
In the presence of the pumping field, the time translation invariance is certainly violated, so that the definition of differential cross-section (4) is modified in the following way: The initial matter state |i is prepared from the remote past (say −∞) state [which is denoted as |i −∞ ] by the action of the pumping field and the interactions in matter systems. The pumping field is assumed to be turned off before the probe photon comes in, H pump (t) = 0 for t > t I . Thus we have where The distribution of the remote past states are descrbied by density matrix: Now the equation (18) is modified to: Then the differential cross-section takes the form The expressions for the operators Dk i ,k f and D † ki,k f [ the equations (A.6,A.7) ] remain unchanged since the pumping field is turned off after t I . Substituting (A.6,A.7) into (68) and defining a correlation functioñ we arrive at The equation (70) is the differential cross section of ILS for systems in nonequilibrium conditions in the form of correlation function.
Next The difference with the equilibrium case (25) lies in the presence of the pumping factor U Hm+ Hpump . Next we follow the same steps as (27,28,29). Employing the composition property of the time evolution operator we obtaiñ Since the pumping field is turned off in the time interval [t I , t F ], the subscripts denoting the Hamiltonian governing time evolution are, in fact, redundant. Namely we can simply take the total Hamiltonian to be H m + H pump with the understanding that the pumping field vanishes for t > t I . Then again by the composition property of evolution operator, the equation (72) reduces tõ where we have inserted an identity U(t F , +∞) U(+∞, t F ) = I. Now the equation (73) is of the same form as that of (31), so that they have the identical form of the following functional integral representation: where the actions S ± now include the driving pumping field term.
With the equations (74,75), ILS in non-equilibrium situation can be computed by using the standard functional integral method. Next let us apply the above results to the case of ILS of the semiconductor pumped by intense laser pulse. The pumping Hamiltonian in the dipole approximation of optical interband transition [24] is given by [ for simplicity two band approximation (one conduction and one valence band ) is made and spin degrees of freedom is ignored, see Chapter 5 of [24] ] where E(t) is the electric field of the pumping pulse laser and d cv is the dipole matrix element.ĉ c,v is the electron annihilation operator for conduction and valence band, respectively. The electric current operator in dipole approximation is given by [w k,q is the optical transition matrix element] In terms of the coherent state variables in Keldysh basis (77) becomes J e,± (q) = 1 2 where Ψ = (ψ 1c , ψ 1v , ψ 2c , ψ 2v ) t , and Then the action for the pumping term is given by We can consider only the intraband electron-phonon interaction owing to band gap.
The action corresponding to (81) is Now everything boils down to the computation of the following correlation function from which the differential cross-section can be obtained via (70).
The correlation function for the one-phonon Raman scattering [corresponding to (59) in equilibrium case with similar notations] is given bỹ whose Feynman diagram is given in figure 5. The electron Green's function satisfy the following Dyson equation which includes the pumping term [ G (85) follows from the standard Feynman rules. The phonon subsystem is also driven into non-equilibrium through the phonon self-energy. Its influence on ILS can be seen in higer order Feynman diagrams [ see figure 6 (b) ]. To best of our knowledge, the ILS in non-equilibrium condition has not been presented in the form of the equations (84,85) before. The detailed investigation of (84,85) is beyond the scope of this paper and is left for future studies.

Discussions
The KS functional integral approach to ILS possesses several advantages over the conventional approach based on Fermi golden rule [22] even in the lowest order calculations. First, the finite temperature effect is naturally incorporated as can be seen in (62), which then allows the discussion of both Stokes and anti-Stokes processes on equal footing. Second, the energy conservation for one-phonon process has been imposed by hand in [22], while in our approach it is a natural consequence of Feynman rule.
True advantages of KS functional approach lie in the study of higher order manybody effects. In terms of Feynman diagrams the many-body effects are often described by the self-energy corrections and the vertex corrections. For example, the Feynman diagram of the lowest order electron self-energy correction to ILS is given by the diagram (a) of figure 6, and the diagram (b) describes the the lowest order phonon self-energy correction. A lowest order vertex correction to ILS is described in figure  7. With phonon self-energy corrections included, the phonon Green's function in (59) can be taken to be fully interacting Green's function, so that the phonon frequency shift and the phonon broadening effect (both stemming from phonon self-energy) for ILS can be precisely addressed.
Both self-energy and vertex correction are known to possess the logarithmic corrections [25], which implies that the inclusion of many-body effects are crucial for the proper interpretation of experimental ILS data. Our formalism provides a natural framework for the investigation of such higher order many effects. The study of the above many-body effects based on Kramers-Heisenberg formula would be exceedingly difficult, if not impossible. In summary, we have expressed the scattering cross section of the resonant ILS in the form of correlation function both in equilibrium and non-equilibrium conditions, which then is recast in the framework of KS functional integral. The correlation function in KS functional integral can be computed by the Feynman diagram perturbation theory which permits the systematic study of many-body effects. Now we are ready to compute the operator Dk i ,k f [ see equations (13) and (14) for definitions ] explicitly. Plugging (A.3) into (A.2) and using the result (A.5) it is straightforward to obtain (ω i,f = c|k i,f |) where x = ( r, t) and k = (k, ω k ) are 4-vectors, and k · x = ω k t − k · r. Similarly, we can obtain The spatial Fourier transform of the electric current operator is defined as J a e (k j , t j ) = d r j e −ikj · rj J a e ( r j , t j ). (A.8)