Hot electron relaxation: An exactly solvable model and improved quantum kinetic equations

V. Meden, C. Wöhler, J. Fricke, and K. Schönhammer

Institut für Theoretische Physik der Universität Göttingen, Bunsenstr.9, D-37073 Göttingen, Germany (March 1, 1995)

Abstract

We present a one-dimensional microscopic model of interacting electrons coupled to phonons which allows the exact calculation of the time evolution of the momentum distribution functions starting from a nonequilibrium initial state. The model shows relaxation features known from approximate calculations for more realistic three-dimensional models. The exact results are compared with approximate ones resulting from the use of quantum kinetic equations which avoid the assumption of completed successive scattering processes implicit in simple Boltzmann equations. To renormalize the electronic but not the phonon propagation in the non-Markovian equations yields no improvement uniform in time compared to simpler approximations. Due to the presence of the Fermi sea only the additional renormalization of the phonon propagation leads to a drastic improvement in comparison with the exact results.

PACS numbers: 71.38.+i, 78.47.+p, 72.10.Bg

Typeset using REVTEX

I. INTRODUCTION

The relaxation of hot carriers has traditionally been described using a Boltzmann equation, where the collision term is calculated assuming successive completed electron-electron and electron-phonon collisions with energy conservation in each scattering process. From the standard derivation of Fermi’s Golden Rule in quantum mechanics textbooks it is obvious that these assumption can certainly not be valid in the short time regime after the nonequilibrium state has been established. The proper inclusion of quantum mechanical coherence leads to quantum kinetic equations of non-Markovian nature. Quantum kinetic equations can be derived by breaking the hierarchy of the equations of motion for the relevant reduced density matrices [1] or by using nonequilibrium Green’s functions within the Keldysh formalism [2,3], which generalizes equilibrium many body techniques, like the use of Feynman diagrams, by introducing time-ordering on a contour running along the real time axis and backwards. As in equilibrium an approximate description of the one-particle Green’s function involves an approximation for the selfenergy. In order to arrive at quantum kinetic equations unfortunately additional approximations have to be made. To get rid of the two-time electronic nonequilibrium correlation functions G^{>} and G^{<} the generalized Kadanoff-Baym ansatz (GKBA) is usually introduced and the remaining retarded and advanced nonequilibrium Green’s functions are replaced by the corresponding equilibrium ones [3]. Recently Schoeller has presented an illuminating diagrammatic analysis to clarify the approximation involved in the GKBA [4]. The phonons are usually treated as a thermal bath [3].

While there are very few microscopic models of interacting electrons and phonons which allow an exact calculation of equilibrium properties the situation concerning exact results for the nonequilibrium behavior of microscopic models is even worse. Due to the uncontrolled approximations involved in the derivation of quantum kinetic equations, which can lead e. g. to negative probabilities, it is highly desirable to obtain exact results for models which show relaxation behavior. As a simplified model to study non-Markovian relaxation, the Jaynes-Cumming model, was used recently to compare results from approximate kinetic equations with exact ones [5]. Due to the simple two level nature of the model no real relaxation occurs.

Interacting fermions in one-dimension were studied intensively after the pioneering work of Tomonaga [6], who showed that the exact eigenstates can be obtained for a long range interaction between the fermions. The calculation of expectation values, like the occupation numbers in momentum space, can elegantly be performed using the technique of the bosonization of the electronic field operator [7]. If the electrons are linearly coupled to phonons the model can still be solved exactly, when the coupling involves long wavelength phonons only [8]. In this paper we show that the method used to calculate exact equilibrium properties can be extended to obtain exact results also for nonequilibrium expectation values, like the momentum distribution for an initial state with hot electrons.

In Sec. II we present the model of one-dimensional electrons coupled to the phonons and exactly calculate various nonequilibrium expectation values. This model shows all the qualitative features known from approximate quantum kinetic treatments of more realistic three-dimensional models [9].

In Sec. III quantum kinetic equations are derived for the electron and phonon momentum distribution in our model closely following the scheme mentioned above. For the small electron-phonon coupling used in most of the figures a simple quasi-particle form for the retarded and advanced equilibrium electron Green’s function agrees very well with the corresponding exact equilibrium functions. We study the equations with unperturbed as well as renormalized phonon propagation. The resulting quantum kinetic equations are non-Markovian. The numerical solution of these equations as well as their Markovian versions and the comparison with the exact results is discussed in Sec. IV. It turns out that none of the approximate kinetic equations with the unperturbed phonon propagation leads to an overall satisfactory description for all time regimes studied. The inclusion of a renormalized phonon propagation in the non-Markovian equation leads to a drastic improvement in the agreement with the exact results.

We summarize our results in Sec. V. The details of the exact calculation of the nonequilibrium electron momentum distribution are presented in an Appendix.

II. THE MODEL AND ITS SOLUTION

To study the relaxation of excited electrons in a microscopic model we start with the Tomonaga-Luttinger (T-L) model, which correctly describes the low energy physics of interacting electrons in one dimension [6,7]. The characteristic feature of the model is an energy dispersion which is linearized around the Fermi points ±k_{F}. This leads to two independent branches of fermions. The dispersion and the assumption of only low momentum scattering processes (no 2k_{F} -backscattering) makes it possible to calculate all many-body eigenstates and energies. Recently there has been progress in the understanding of the equilibrium electronic spectral properties of the model, by exact calculations of the one-particle spectral function [10,11]. Contrary to interacting electrons in three dimensions the low energy physics of this one-dimensional model of interacting electrons cannot be described by the Fermi liquid (FL) scenario. The momentum resolved one-electron spectral function at k = k_{F} shows a power law singularity with an interaction dependent exponent and no delta peak occurs as in the FL-case, the momentum distribution function shows a nonanalytic behavior at k_{F} but no discontinuity, and spin and charge degrees of freedom decouple (spin-charge separation) [10,11]. All these phenomena are characteristic for so called Luttinger liquids (LL) [7].

The key to the exact solution for the eigenstates of the T-L model is the bosonization of the Hamiltonian. The calculation of arbitrary equilibrium expectation values also requires the bosonization of the electronic field operator ψ^{†}(x) [7].

In this paper we show that this bosonization of the field operator is also the essential ingredient for the theoretical description of nonequilibrium phenomena within this model. As in the following we study relaxation processes due to the coupling to phonons, we simplify the electronic model in the present publication as much as possible in order to reduce the numerical effort. The complete model will be considered in a forthcoming publication [12]. We consider only one branch of spinless electrons with the dispersion ε_{k} = v_{F} (k−k_{F}) (k_{F} = 0 in the following). The dispersion is assumed to be linear down to k → −∞, with the Dirac sea of negative energy states filled in the groundstate and the energy measured relative to the groundstate. Because the one branch model is a nontrivial model of interacting electrons the one-particle spectral function for momenta different from k_{F} shows a nontrivial structure [10] but the k = k_{F} spectral function consists of a delta peak of weight one and the momentum distribution function is identical to the noninteracting one, as the interacting groundstate is given by the noninteracting one. In the “g-ology” classification [13] this is the g_{4}-model for which the anomalous dimension is zero.

The equilibrium properties of the T-L model can still be calculated exactly if an electronphonon coupling of the Holstein type is introduced which only involves long wavelength phonons [8,14]. As we discuss here only a single electronic branch (“chiral fermions”) we introduce a corresponding simplification for the phonons and in the form of the electronphonon coupling. The full T-L model with the usual Holstein type coupling can be treated with similar methods as those presented below [12].

The Hamiltonian for a finite system of length L and phonon dispersion ω_{q} is

H = H_{e} + H_{p} + H_{ee} + H_{ep} ,  (1)

with

\begin{aligned}H_{e} & = & ∑_{k} ε_{k} (a_{k}^{†}a_{k} − ⟨a_{k}^{†}a_{k}⟩_{0}) , H_{p} = ∑_{q>0} ω_{q}B_{q}^{†}B_{q} , \\ H_{ee} & = & \frac{1}{2L} ∑_{q≠0} ṽ(q)ρ_{q}ρ_{−q} , \\ H_{ep} & = & (\frac{2π}{L} )^{1/2} ∑_{{k \atop q>0}} g(q) (a_{k+q}^{†}a_{k}B_{q} + B_{q}^{†}a_{k}^{†}a_{k+q}) , \end{aligned}  (2)

where a_{k} (B_{q}) is the annihilation operator of an electron (phonon) and < ... >_{0} is the expectation value in the noninteracting groundstate. Here ρ_{q} denotes the Fourier transform of the density of electrons

ρ_{q} = ∑_{k} a_{k}^{†}a_{k+q} ,  (3)

for q ≠ 0. For q = 0 we define the electron number operator relative to the groundstate as

N̂ = ∑_{k} (a_{k}^{†}a_{k} − ⟨a_{k}^{†}a_{k}⟩_{0}) .  (4)

ṽ(q) is the Fourier transform of the two-particle electron-electron interaction and g(q) the coupling function between electrons and phonons. With a proper normalization the densities ρ_{q} obey Bose commutation relations [15]. If one defines

b_{q} = (\frac{2π}{qL})^{1/2} ρ_{q} , b_{q}^{†} = (\frac{2π}{qL})^{1/2} ρ_{−q} ,  (5)

for q > 0, the commutation relations read

[b_{q}, b_{q′}^{†}] = δ_{q,q′} , [b_{q}, b_{q′}] = 0 .  (6)

It is possible to write H_{ee} bilinearly in b_{q} (b_{q}^{†}) and also to express the electron-hole excitation in H_{ep} linearly in the boson operators b_{q} (b_{q}^{†}). For one-dimensional electrons with a linear dispersion it is also possible to write the kinetic energy of the electrons as [16,7]

H_{e} = ∑_{q>0} ε_{q}b_{q}^{†}b_{q} + c(N̂) ,  (7)

where c(N̂) only contains the electronic particle number operator. In the following this term can be neglected because we are working in a sector of the Hilbert space with constant particle number. Now we can formulate Eq. (2) in terms of the bosonic electron densities and the phonons as

H = H̃_{e} + H̃_{p} + H̃_{ep}  (8)

with

\begin{aligned}H̃_{e} & = & ∑_{q>0} ε̃_{q}b_{q}^{†}b_{q} , H̃_{p} = ∑_{q>0} ω_{q}B_{q}^{†}B_{q} , \\ H̃_{ep} & = & ∑_{q>0} \sqrt{q}g(q) (b_{q}^{†}B_{q} + B_{q}^{†}b_{q}) , \end{aligned}  (9)

and ε̃_{q} = v_{F} q (1 + \frac{ṽ(q)}{2πv_{F}} ). As one can see from Eq. (9) the groundstate of the model is still the tensor product of the Fermi sea with the phonon vacuum. Therefore all groundstate expectation values like the momentum distribution are given by the noninteracting ones. The anomalous dimension is zero as in the g_{4}-model without the electron-phonon coupling. In the model with two branches an additional canonical transformation is necessary to write the Hamiltonian in a form similar to Eq. (9). In the full model the groundstate contains electron-phonon correlations [12].

Using a canonical transformation the problem of the coupled bosonic modes Eq. (9) can be diagonalized [8]. The transformation is given by

\begin{aligned}b_{q} & = & α_{q}c_{q} − β_{q}s_{q} , \\ B_{q} & = & α_{q}s_{q} + β_{q}c_{q} , \end{aligned}  (10)

with

c_{q} := cos (γ_{q}/2) , s_{q} := sin (γ_{q}/2)  (11)

and

tan (γ_{q}) = \frac{2g(q)\sqrt{q}}{ε̃_{q} − ω_{q}} .  (12)

In the new bosonic operators α_{q} and β_{q} the Hamiltonian reads

H = ∑_{q>0} [λ_{+}(q)α_{q}^{†}α_{q} + λ_{−}(q)β_{q}^{†}β_{q}] ,  (13)

with the energies

λ_{±}(q) := \frac{1}{2} \{ε̃_{q} + ω_{q} ± [(ε̃_{q} − ω_{q})^{2} + 4g^{2}(q)q]^{1/2}\} .  (14)

In order to obtain a stable ground state the boson energies λ_{±}(q) have to be positive [8,17]. This leads to a restriction of the coupling strength g(q) that can be used. With Eqs. (13) and (14) all eigenstates and energies are determined in terms of two bosonic modes.

To calculate electronic properties like the one-particle Green’s function and the momentum distribution function we use the bosonization of the electronic field operator

ψ^{†}(x) = \frac{1}{\sqrt{L}} ∑_{k} e^{−ikx}a_{k}^{†} .  (15)

One has the operator identity [7]

ψ^{†}(x) = \frac{e^{−ixπ/L}}{\sqrt{L}} e^{−iΦ^{†}(x)}U^{†}e^{−iΦ(x)} ,  (16)

with

Φ(x) = \frac{π}{L}N̂x − i ∑_{q>0} e^{iqx} (\frac{2π}{Lq} )^{1/2} b_{q} ,  (17)

where U^{†} denotes a fermionic raising operator which commutes with the b_{q} and maps the N-electron groundstate to the (N + 1)-electron groundstate.

To discuss relaxation in its elementary form an additional electron is created in the groundstate |ϕ_{0} > at time t = 0. The resulting many-particle state a_{k_{0}}^{†}|ϕ_{0} > is not an eigenstate of H and the excited electron will relax. We are interested in the evolution of the momentum distribution functions

f(k, k_{0}, t) :=< ϕ_{0}|a_{k_{0}}a_{k}^{†}(t)a_{k}(t)a_{k_{0}}^{†}|ϕ_{0} >  (18)

and

h(q, k_{0}, t) :=< ϕ_{0}|a_{k_{0}}B_{q}^{†}(t)B_{q}(t)a_{k_{0}}^{†}|ϕ_{0} > ,  (19)

as well as the evolution of the distribution of the excitation energy between the electron and phonon system

< H̃_{i} > (t) :=< ϕ_{0}|a_{k_{0}} H̃_{i}(t)a_{k_{0}}^{†}|ϕ_{0} >  (20)

for i = e, p, ep. The nonequilibrium phonon distribution function h(q, k_{0}, t) and the energy contributions < H̃_{i} > (t) can be calculated without the bosonization of the electronic field operator. Using Eqs. (10) and (13) one obtains

B_{q}(t) = c_{q}s_{q} (e^{−iλ_{+}(q)t} − e^{−iλ_{−}(q)t}) b_{q}
+ (s_{q}^{2}e^{−iλ_{+}(q)t} + c_{q}^{2}e^{−iλ_{−}(q)t}) B_{q} .  (21)

As B_{q}a_{k_{0}}^{†}|ϕ_{0} >= 0 the second term does not contribute in the calculation of h(q, k_{0}, t). With

[b_{q}, a_{k_{0}}^{†}] = (\frac{2π}{qL})^{1/2} a_{k_{0}−q}^{†} the application of B_{q}(t) to the initial state simplifies to

B_{q}(t)a_{k_{0}}^{†}|ϕ_{0} >= (\frac{2π}{qL})^{1/2} c_{q}s_{q} (e^{−iλ_{+}(q)t} − e^{−iλ_{−}(q)t}) a_{k_{0}−q}^{†}|ϕ_{0} > .  (22)

Taking the scalar product of this state with itself yields

h(q, k_{0}, t) = 2\frac{2π}{qL}c_{q}^{2}s_{q}^{2} (1 − cos \{[λ_{+}(q) − λ_{−}(q)] t\})
×Θ(k_{0} − q) .  (23)

The application of b_{q}(t) to the initial state can be calculated similarly

b_{q}(t)a_{k_{0}}^{†}|ϕ_{0} >= (\frac{2π}{qL})^{1/2} (c_{q}^{2}e^{−iλ_{+}(q)t} + s_{q}^{2}e^{−iλ_{−}(q)t}) a_{k_{0}−q}^{†}|ϕ_{0} > ,  (24)

and yields

< ϕ_{0}|a_{k_{0}}b_{q}^{†}(t)b_{q}(t)a_{k_{0}}^{†}|ϕ_{0} > = Θ(k_{0} − q)\frac{2π}{qL}
× (c_{q}^{4} + s_{q}^{4} + 2c_{q}^{2}s_{q}^{2} cos {[λ_{+}(q) − λ_{−}(q)] t}) .  (25)

The conservation law < B_{q}^{†}B_{q} + b_{q}^{†}b_{q} > (t) = const., which follows from the special form of the electron-phonon interaction H̃_{ep} in Eq. (9), is fulfilled. We obtain for < H̃_{e} > (t) and < H̃_{p} > (t)

< H̃_{e} > (t) = ∑_{q>0} ε_{q} < ϕ_{0}|a_{k_{0}}b_{q}^{†}(t)b_{q}(t)a_{k_{0}}^{†}|ϕ_{0} >  (26)

and

< H̃_{p} > (t) = ∑_{q>0} ω_{q}h(q, k_{0}, t) .  (27)

The calculation of

< H̃_{ep} > (t) = ∑_{q>0} \sqrt{q}g(q)2Re \{< ϕ_{0}|a_{k_{0}}B_{q}^{†}(t)b_{q}(t)a_{k_{0}}^{†}|ϕ_{0} >\}  (28)

involves

Re \{< ϕ_{0}|a_{k_{0}}B_{q}^{†}(t)b_{q}(t)a_{k_{0}}^{†}|ϕ_{0} >\} = \frac{2π}{qL}c_{q}s_{q} (c_{q}^{2} − s_{q}^{2})
× (1 − cos {[λ_{+}(q) − λ_{−}(q)] t}) Θ(k_{0} − q) .  (29)

Using the explicit form of the < H̃_{i} > (t) it is easy to show that the expectation value of the total Hamiltonian has the time independent value v_{F} k_{0}. The results for h(q, k_{0}, t) and < H̃_{i} > (t) presented above are very simple as bosonic operators in the Heisenberg picture can easily be calculated as in Eq. (21). The fermionic operators a_{k}(t) in the Heisenberg picture, on the other hand, can only be obtained using the bosonization of the field operator ψ(x). As the function Φ(x) in the exponents in Eq. (16) is linear in the boson operators b_{q}, the main tool in the calculation of f(k, k_{0}, t) is the repeated use of the Baker-Hausdorff formula. Details of the calculation are given in the Appendix. In order to be able to perform the various x-integrations, that occur due to the expression of a_{k} in terms of ψ(x), we use an iterative approach we developed earlier for the calculation of one-particle spectral functions [10]. This leads to

f(k_{n}, k_{n_{0}}, t) = ∑_{r=max(n,0)}^{n_{0}} ∑_{s=0}^{n_{0}−r} ∑_{l=0}^{min(r+s,r−n)} a_{r+s−l}^{(n_{c})}(t)
× (a_{r}^{(n_{c})}(t))^{∗} b_{l}^{(n_{c})}(t) (b_{s}^{(n_{c})}(t))^{∗} ,  (30)

for n_{0} ≥ n > 0, where q_{c} = \frac{2π}{L}n_{c} is a cutoff for the momentum transfer due to electron-electron and electron-phonon interaction (and the maximum phonon momentum) and k_{n} = \frac{2π}{L}n. The time dependent coefficients a_{n}^{(m)} and b_{n}^{(m)} are determined iteratively for m > 1, l ∈ IN_{0}, and i = 0, ..., m − 1 by

\begin{aligned}a_{lm+i}^{(m)}(t) & = & ∑_{j=0}^{l} \frac{[u_{m}(t)/m]^{j}}{j!} a_{m(l−j)+i}^{(m−1)}(t) , \\ b_{lm+i}^{(m)}(t) & = & ∑_{j=0}^{l} \frac{[−u_{m}(t)/m]^{j}}{j!} b_{m(l−j)+i}^{(m−1)}(t) , \end{aligned}  (31)

with the starting values (m = 1)

\begin{aligned}a_{l}^{(1)}(t) & = & ∑_{k=0}^{l}[e(t)]^{l−k} \frac{[u_{1}(t)]^{k}}{k!} , l ∈ IN_{0} , \\ b_{l}^{(1)}(t) & = & \frac{[−u_{1}(t)]^{l}}{l!} − \frac{[−u_{1}(t)]^{l−1}}{(l − 1)!} e(t) , l ∈ IN , \\ b_{0}^{(1)}(t) & = & 1 . \end{aligned}  (32)

The functions u_{l}(t) and e(t) are defined as

\begin{aligned}u_{l}(t) & := & c_{q_{l}}^{2}e^{−iλ_{+}(q_{l})t} + s_{q_{l}}^{2}e^{−iλ_{−}(q_{l})t} − e^{−iv_{F} q_{l}t} , \\ e(t) & := & e^{−iv_{F} t2π/L} . \end{aligned}  (33)

Performing the iteration on a computer leads to the exact numerical result for f for large but finite systems. As a test of our program calculating f, we determined the energy in the electronic system for a given time via the electron momentum distribution, and compared the result with Eq. (26) calculated for the same system size. We found both calculations to agree within our numerical accuracy. As an additional test we used the electron number conservation.

As can be seen from Eqs. (26), (27), and (28) it is possible to carry out the thermodynamic limit in the calculation of the energy expectation values by replacing the sum by an integral. At a time t_{r} the energies < H̃_{i} > (t) for a finite system start to deviate appreciably from the values for the infinite system due to recurrence effects always present in finite systems. It is important to know t_{r} because we calculate f for large but finite systems.

From now on we restrict ourselves to the relaxation due to the coupling of the electrons to an optical phonon mode and set ṽ(q) ≡ 0 and ω_{q} ≡ ω_{0}. We also have to specify the coupling function g(q). In this paper we use two forms. The figures shown in this section are calculated for a constant coupling g(q) = gΘ(q_{c} − q). We also performed all the calculations for a Fröhlich type coupling g(q) = \frac{g}{q+q_{0}} Θ(q_{c} −q), with a low momentum cutoff q_{0} in order to avoid divergencies in our one-dimensional model. To guarantee λ_{−}(q) ≥ 0 for all q we have to fulfill a relation between q_{0} and g. There are only some minor differences between the results with the two forms of the coupling, which we point out later. In Sec. IV we also show figures for the Fröhlich coupling. The constant coupling function leads to a dimensionless coupling constant Γ := g^{2}/(v_{F}^{2} q_{c}), where v_{F} q_{c} is related to the bandwith in a more realistic model.

Fig. 1 shows the expectation values of the various energies as a function of time for L → ∞. For t = 0 the electron system contains all the excitation energy. Afterwards a part of the energy is transferred to the phonon system and for large times all < H̃_{i} > (t) tend to an asymptotic value in an oscillatory way. Using the stationary phase method it is easy to determine the asymptotic behavior of Eqs. (26), (27), and (28). For couplings not too strong (2Γ < ω_{0}/(v_{F} q_{c})) and k_{0}/q_{c} not too small (ω_{0}/(v_{F} q_{c}) − k_{0}/q_{c} < 2Γ), such that the stationary point lies in the integration interval [0, k_{0}], the oscillation frequency is given by the minimum of the mode energy difference λ_{+}(q) − λ_{−}(q) and the oscillation falls off as 1/\sqrt{t}. The asymptotic values are obtained by neglecting the cosine-terms in Eqs. (26), (27), and (28). Increasing the coupling strength the asymptotic energy in the phonon system increases as expected. The same behavior also occurs for a Fröhlich coupling, but with different bounds for the coupling strength and excitation momentum. If in the model with a Fröhlich coupling one chooses a q_{0} which leads to a nearly smooth crossover at q_{c}, i. e. g(q_{c}) is very small compared to g(0), there is an additional modulation of < H̃_{p} > (t) and < H̃_{ep} > (t) with frequency ω_{0} and a small amplitude. This is because in contrast to the constant coupling the cosine-prefactors in (27) and (28) c_{q}^{2}s_{q}^{2}/q and c_{q}s_{q}(c_{q}^{2} −s_{q}^{2})/\sqrt{q} have their maximum at q → 0. Therefore the weight of the oscillation with the frequency λ_{+}(q → 0) − λ_{−}(q → 0) = ω_{0} is enhanced. For < H̃_{e} > (t) the prefactor is given by c_{q}^{2}s_{q}^{2} and as c_{q}^{2} ∼ q and s_{q}^{2} ∼ 1 for q → 0, for all coupling functions which are finite at q = 0, this additional oscillation is not seen in < H̃_{e} > (t). By a variation of q_{0} we could discuss the continuous transition from a Fröhlich coupling to a constant coupling.

In Fig. 2 we present the phonon distribution h as a function of the momentum q/q_{c} and the time for a finite system. The system size is determined by the dimensionless noninteracting energy level spacing ν := \frac{v_{F} 2π/L}{v_{F} q_{c}}. As one can read off Eq. (23) h oscillates for a given q/q_{c} with the mode energy difference at the given wave vector. In a calculation of h in a Boltzmann approximation, with momentum and energy conservation in each elementary scattering process (see below), there would only be phonons at momentum q_{B} := ω_{0}/v_{F} (here q_{B}/q_{c} = 0.2). As Fig. 2 shows other momenta have an appreciable occupation in the exact solution. The occupation number for each wave vector shows undamped oscillations. Nevertheless for t → ∞ < H̃_{p} > (t) tends to a time independent value due to interference in the thermodynamic limit. For the Fröhlich coupling h looks almost the same but as discussed above the weight at small momenta is enhanced compared to the constant coupling.

Fig. 3 shows the evolution of the electron distribution f for k_{0}/q_{c} = 1.5. The distribution is plotted for momenta down to k/q_{c} = 0.8 only. The slight spreading of the Fermi sea for the times used is not shown. To present the behavior of f in one figure the occupation number at k_{0}/q_{c} is scaled by a factor 1/5. For the system size chosen (v_{F} q_{c})t_{r} is approximately 500, i. e. much larger than the largest time shown. In Fig. 3 one can see pronounced features at wave vectors k = k_{0} − nq_{B} with n ∈ IN. An approximate treatment using a Boltzmann equation leads to nonzero values of f only at these momenta. Contrary to that, at small t the exact distribution f shows a broad maximum around k_{0} − q_{B}, which sharpens up later, but never reaches the sharpness of the initial distribution at k_{0}. The same holds for higher replicas which build up as the time increases. In addition to the sharpening the weight in the peak increases until it reaches a maximum value, and decreases later on. For even larger times more and more weight is distributed evenly not restricted to the replication regions. For increasing coupling strength the replication tendency becomes more and more irrelevant, i. e. the weight is distributed more evenly over the momenta. For the Fröhlich coupling one difference to the behavior described above has to be pointed out. For small t the distribution of the weight over the k-states is a copy of the q-dependence of g(q) and only for larger t we find the replication structure discussed above.

III. KINETIC EQUATIONS

In this section we discuss the approximate description of the relaxation process using kinetic equations. In the derivation of the kinetic equations for our model Eq. (2), we restrict ourselves from the beginning to the relaxation due to the electron-phonon interaction with an optical phonon mode. Apart from the fact that our model contains phonons with only positive momenta the derivation of quantum kinetic equations for nonequilibrium Green’s functions within the Keldysh formalism using a GKBA is very similar to the derivations presented in Refs. [2] and [3]. The initial state a_{k_{0}}^{†}|ϕ_{0} > introduced in Sec. II is an eigenstate of the noninteracting Hamiltonian and allows a Wick’s theorem, i. e. no initial correlations [4] occur. From the Keldysh formalism one obtains for the derivative of the electron momentum distribution n_{k}(t)

\frac{d}{dt}n_{k}(t) = − ∫_{0}^{t} dt′ \{Σ_{k}^{>}(t, t′)G_{k}^{<}(t′, t) − Σ_{k}^{<}(t, t′)G_{k}^{>}(t′, t)

+G_{k}^{<}(t, t′)Σ_{k}^{>}(t′, t) − G_{k}^{>}(t, t′)Σ_{k}^{<}(t′, t)} ,  (34)

where Σ_{k}^{>} and Σ_{k}^{<} denote the nonequilibrium selfenergies and the electron correlation functions G_{k}^{<}(t, t′) and G_{k}^{>}(t, t′) are defined as

\begin{aligned}iG_{k}^{<}(t, t′) & = & − ⟨a_{k}^{†}(t′)a_{k}(t)⟩ , \\ iG_{k}^{>}(t, t′) & = & ⟨a_{k}(t)a_{k}^{†}(t′)⟩ \end{aligned}  (35)

where < ... > denotes the expectation value with the initial state. As in [2] and [3] we neglect vertex corrections in the equations for the nonequilibrium selfenergies

\begin{aligned}Σ_{k}^{>}(t, t′) & = & \frac{2π}{L} ∑_{q>0} g^{2}(q) \{iD_{q}^{>}(t, t′)G_{k−q}^{>}(t, t′) + iD_{q}^{<}(t′, t)G_{k+q}^{>}(t, t′)\} , \\ Σ_{k}^{<}(t, t′) & = & \frac{2π}{L} ∑_{q>0} g^{2}(q) \{iD_{q}^{<}(t, t′)G_{k−q}^{<}(t, t′) + iD_{q}^{>}(t′, t)G_{k+q}^{<}(t, t′)\} , \end{aligned}  (36)

with the phonon correlation functions D_{q}^{<}(t, t′) and D_{q}^{>}(t, t′) given by

\begin{aligned}iD_{q}^{<}(t, t′) & = & ⟨B_{q}^{†}(t′)B_{q}(t)⟩ , \\ iD_{q}^{>}(t, t′) & = & ⟨B_{q}(t)B_{q}^{†}(t′)⟩ . \end{aligned}  (37)

We neglect vertex corrections in Eq. (36) although it is known that already in equilibrium the Migdal theorem [18] does not hold in one dimension [14]. In the diagrams presenting vertex corrections there is no additional factor ω_{0}/B, where B denotes the bandwidth (here B ∼ v_{F} q_{c}), compared to other higher order diagrams. Therefore for the equilibrium theory all diagrams of a given order in the coupling lead to contributions of the same order even if ω_{0}/B is small. To improve the quality of the derived kinetic equations we therefore choose small couplings so that vertex corrections are small due to the small coupling constant instead of the factor ω_{0}/B present in three dimensions. To neglect vertex corrections in the equilibrium selfenergy is not a serious approximation for the parameters used in the following. In Fig. 4 we compare the exact equilibrium spectral function ρ^{>}(k_{0}, ω) (solid line) with the one calculated in second order perturbation theory for the selfenergy (dotted line) for the constant coupling function described in Sec. II. As parameters we choose the values used in the comparison of the exact nonequilibrium electron distribution and the one determined from the approximate kinetic equations in the next section. The exact spectral function is calculated with a method we presented in [14] and [10]. In the exact calculation we have to add a small Lorentzian broadening, which is also added in the approximation using the second order selfenergy. For the small coupling constant chosen the agreement between both results is very good and this also holds for all other momenta. Therefore all higher order terms in the selfenergy, including vertex corrections, are unimportant in equilibrium and it seems plausible to assume that this also holds in nonequilibrium. The same statement holds for a Fröhlich coupling if the coupling constant is chosen small enough.

We use the GKBA in order to replace the electron correlation functions by the momentum distribution function and the retarded and advanced Green’s functions [2,4]

\begin{aligned}G_{k}^{>}(t, t′) & = & \{G_{k}^{r}(t, t′)(1 − n_{k}(t′)) \\ & & −(1 − n_{k}(t))G_{k}^{a}(t, t′)\} , \\ G_{k}^{<}(t, t′) & = & − \{G_{k}^{r}(t, t′)n_{k}(t′) − n_{k}(t)G_{k}^{a}(t, t′)\} . \end{aligned}  (38)

with

\begin{aligned}G_{k}^{r}(t, t′) & = & −iΘ(t − t′) ⟨\{a_{k}(t), a_{k}^{†}(t′)\}⟩ , \\ G_{k}^{a}(t, t′) & = & iΘ(t′ − t) ⟨\{a_{k}(t), a_{k}^{†}(t′)\}⟩ . \end{aligned}  (39)

In a further approximation the remaining Green’s functions are replaced by their equilibrium versions. Here we choose for the equilibrium Green’s function a damped quasi-particle approximation [3]

Gk^{{r \atop a}} (t, t′) = {− \atop +}iΘ({+ \atop −}(t − t′))e^{{− \atop +} (γ_{k} {+ \atop −} iE_{k})(t−t′)} ,  (40)

where E_{k} = ε_{k} + Re{Σ^{(2)}(k, E_{k})} ≈ ε_{k} + Re{Σ^{(2)}(k, ε_{k})} and γ_{k} = |Im{Σ^{(2)}(k, E_{k})}| ≈

|Im{Σ^{(2)}(k, ε_{k})}| are given by the second order selfenergy Σ^{(2)}. For the constant coupling Σ^{(2)}(k, ε_{k}) is given by

\frac{Im \{Σ^{(2)}(k, ε_{k})\}}{v_{F} q_{c}} = \{ \begin{matrix}−πΓ : & ω_{0} < ε_{k} \\ 0 : & else \end{matrix}  (41)

and

\frac{Re \{Σ^{(2)}(k, ε_{k})\}}{v_{F} q_{c}} = \{ \begin{matrix}−Γ log |\frac{ω_{0}}{ω_{0}−ε_{k}} | : & ε_{k} < v_{F} q_{c} \\ −Γ log |\frac{ω_{0}}{ω_{0}−v_{F} q_{c}} | : & ε_{k} > v_{F} q_{c} \end{matrix}  (42)

for k > 0 and ω_{0} < v_{F} q_{c}. Similar expressions can be given for the Fröhlich coupling. The dashed line in Fig. 4 shows the spectral function for the constant coupling in this approximation, which is a Lorentzian around E_{k}. The width of the curve is given by the sum of γ_{k} and the artificial broadening. The Lorentzian agrees very well with the exact solution. For weak couplings the quasi-particle approximation holds for all momenta, excluding a small momentum range around q_{B}, where the exact equilibrium spectral function changes its character [19,14]. As we will see in the next section this momentum range does not contribute in the relaxation for the time regimes studied. Therefore the approximation of the equilibrium retarded and advanced electron Green’s functions by damped quasi-particles should be very good in the derivation of kinetic equations if the coupling is not too strong. The same holds for the Fröhlich coupling.

In a last step in the derivation we have to find an expression for the phonon correlation functions D^{<} and D^{>}. It is a special feature of our model that these functions can be easily calculated exactly using Eq. (21). It turns out to be more consistent to use again the GKBA in its usual form

\begin{aligned}D_{q}^{>}(t, t′) & = & D_{q}^{r}(t, t′)(1 + N_{q}(t′)) − (1 + N_{q}(t))D_{k}^{a}(t, t′) , \\ D_{q}^{<}(t, t′) & = & D_{q}^{r}(t, t′)N_{q}(t′) − N_{q}(t)D_{q}^{a}(t, t′) . \end{aligned}  (43)

with

\begin{aligned}D_{q}^{r}(t, t′) & = & −iΘ(t − t′) ⟨[B_{q}(t), B_{q}^{†}(t′)]⟩ , \\ D_{q}^{a}(t, t′) & = & iΘ(t′ − t) ⟨[B_{q}(t), B_{q}^{†}(t′)]⟩ \end{aligned}  (44)

and the phonon momentum distribution function N_{q}(t). In the literature the phonons are usually treated as a thermal bath at temperature T, i. e. N_{q}(t) is replaced by the time independent thermal distribution N_{q}^{th}, and the nonequilibrium phonon propagation is replaced by the equilibrium noninteracting one [3,9]

(Dq^{{r \atop a}} )^{(0)} (t, t′) = {− \atop +}iΘ({+ \atop −}(t − t′))e^{−iω_{0}(t−t′)} .  (45)

In order to have a more consistent comparison with the exact solution in this paper we treat the phonon momentum distribution to be time dependent. As the presence of the Fermi sea leads to an important renormalization of the phonons the equilibrium phonon propagation used in Eq. (43) is not assumed to be the unperturbed one. We calculate the second order equilibrium phonon selfenergy Π^{(2)}(q, ω)

Re \{Π^{(2)}(q, ω)\} = \frac{g^{2}(q)}{v_{F}} \frac{v_{F} q}{ω − v_{F} q}  (46)

and

Im \{Π^{(2)}(q, ω)\} = −π \frac{g^{2}(q)}{v_{F}}v_{F} qδ(ω − v_{F} q)  (47)

and obtain for the corresponding phonon propagator (Dq^{{r \atop a}} )^{(2)} (t, t′)

(Dq^{{r \atop a}} )^{(2)} (t, t′) = {− \atop +}iΘ({+ \atop −}(t − t′)) (s_{q}^{2}e^{−iλ_{+}(q)(t−t′)} + c_{q}^{2}e^{−iλ_{−}(q)(t−t′)}) ,  (48)

which is already the exact result for the equilibrium phonon propagation. The fact that the RPA-like resummation of all electron bubbles leads to the exact phonon Green’s function is well known for our model [20]. Generally the second order phonon selfenergy will of course not lead to the exact result, but we expect a similar importance of this improved description of the phonon propagation also for more realistic models.

With these approximations the kinetic equations for the electron distribution reads

\frac{d}{dt}n_{k}(t) = −2\frac{2π}{L} ∑_{q>0} g^{2}(q) ∫_{0}^{t} dt′
× ([s_{q}^{2}e^{−(γ_{k−q}+γ_{k})(t−t′)} cos [(λ_{+}(q) + E_{k−q} − E_{k})(t − t′)]
× {n_{k}(1 − n_{k−q})(N_{q} + 1) − n_{k−q}(1 − n_{k})N_{q}}_{t′}
+s_{q}^{2}e^{−(γ_{k}+γ_{k+q})(t−t′)} cos [(λ_{+}(q) + E_{k} − E_{k+q})(t − t′)]
× {n_{k}(1 − n_{k+q})N_{q} − n_{k+q}(1 − n_{k})(N_{q} + 1)}_{t′}]
+ [s_{q}^{2}, λ_{+}(q) → c_{q}^{2}, λ_{−}(q)]) .  (49)

Following corresponding steps and approximations as above we obtain as the kinetic equations for the phonon distribution

\frac{d}{dt}N_{q}(t) = −2\frac{2π}{L} g^{2}(q) ∑_{k} ∫_{0}^{t} dt′
× ([s_{q}^{2}e^{−(γ_{k−q}+γ_{k})(t−t′)} cos [(λ_{+}(q) + E_{k−q} − E_{k})(t − t′)]
× {−n_{k}(1 − n_{k−q})(1 + N_{q}) + n_{k−q}(1 − n_{k})N_{q}}_{t′}]
+ [s_{q}^{2}, λ_{+}(q) → c_{q}^{2}, λ_{−}(q)]) .  (50)

For γ_{k} = 0, E_{k} = ε_{k}, s_{q}^{2} = 1, c_{q}^{2} = 0, and λ_{+}(q) = ω_{0} these non-Markovian equations with memory (49) and (50) could also be derived by the decoupling of the equations of motion [1].

From the non-Markovian equations (49) and (50) simplified Markovian versions follow by changing the time dependence of the momentum distribution functions in the integral from t′ to t. The time integration can then be performed. Setting the upper integration boundary equal to infinity (completed collision approximation) yields the Boltzmann approximation. In the next section we also discuss the results obtained by these additional approximations. For the non-Markovian approximation of the electron distribution function Eq. (49) the damping factor γ_{k} is very important, in order to supress the well known deficiency of momentum regimes with negative spectral weight [3,5]. For γ_{k} = 0 this occurs already at small times and the results become totally unphysical for larger times. The neglect of the damping in Eq. (50) does not lead to such drastic effects.

The solution of Eqs. (49) and (50) with the initial conditions used in Sec. II leads to an approximation for the exact distributions described in Eqs. (18) and (19). We solved the equations numerically by standard routines. For the non-Markovian case we mapped the system of integro-differential equations to a larger system of differential equations by introducing new functions. As a test for the numerical stability of the solution we always checked the particle number by summing n_{k} over k. As one can see from Eq. (49) particle number conservation holds as in the exact solution.

IV. COMPARISON OF THE EXACT SOLUTION AND THE APPROXIMATIONS

In this section we compare the exact results for the distribution functions described in Sec. II with the results calculated with help of the kinetic equations presented in Sec. III. In the first subsection we compare the exact results with results from the kinetic equations (49) and (50) where we replace the phonon propagation by the unperturbed one. We obtain the corresponding kinetic equations by replacing s_{q}^{2} → 1, c_{q}^{2} → 0, and λ_{+}(q) → ω_{0} in Eqs. (49) and (50).

A. Unperturbed phonon propagation

In addition to this approximation for the phonon propagation the phonons are usually treated as a thermal bath [3,9]. Here we also solve Eq. (50) to determine N_{q}(t). As for the initial condition studied the phonon occupation numbers are very small for all times, there are only some minor differences in n_{k}(t) if the phonons are treated as time dependent or as a thermal bath at T = 0 (see below). In Figs. 5-7 we compare the exact solution for the electron momentum distribution with the Markovian approximations with γ_{k} = 0, E_{k} = ε_{k} and γ_{k} ≠ 0, E_{k} and the non-Markovian approximation with γ_{k} ≠ 0, E_{k} for the model with a constant coupling. Within the accuracy of the plotting there is no difference, if we replace the quasi-particle dispersion E_{k} by the noninteracting one ε_{k}, because the small difference E_{k} − ε_{k} cancels out for allmost all relevant momenta in the expression E_{k−q} − E_{k} because E_{k} − ε_{k} is constant for k > q_{c}.

For the small coupling constant chosen and not too large times (v_{F} q_{c})t < 150, the qualitative overall agreement between the exact solution and these approximations is reasonable.

All the plotted approximations show the replication behavior with the first broad and later sharpened peaks around k_{0} − nq_{B}. We want to point out, that a similar kind of replication, sharpening, and shifting of weight through the replicas has been observed in approximate calculations for three-dimensional models [9]. This shows that our one-dimensional model contains essential parts of the physics of relaxation also present in more realistic threedimensional models. The other approximations mentioned in Sec. III lead to a qualitatively wrong behavior already at small times. A Boltzmann approximation with γ_{k} = 0 (energy conservation in each elementary scattering process) and E_{k} = ε_{k} leads to a sharp replication of the weight at t = 0 at the momenta k = k_{0} − nq_{B} with n ∈ IN contrary to the exact solution. A Boltzmann approximation with γ_{k} ≠ 0 (energy conservation with uncertainty γ_{k} instead of exact energy conservation in each scattering process) leads to a broadened replication, but not to the sharpening shown at later times by the exact solution. As demonstrated in Sec. II for even larger times there is a spreading of the weight around all momenta. This cannot be produced by a Boltzmann approximation. The Boltzmann approximations therefore strongly deviates from the exact distribution for larger times. As mentioned above the non-Markovian approximation with γ_{k} = 0 leads to unphysical negative values for the momentum distribution already at very small times.

For times earlier than those shown in Figs. 5-7 the Markovian approximation with γ_{k} = 0 works surprisingly well. As it tends towards the incorrect Boltzmann limit for large times its agreement with the exact solutions deteriorates from Fig. 5 to Fig. 7. The Markovian approximation with γ_{k} ≠ 0 provides a better approximation in the time range shown in the figures. It is surprising that the non-Markovian approximation including quasi-particle corrections is the worst approximation for the time (v_{F} q_{c})t = 60 in Fig. 5 and even for larger times it is not the overall best. The main deficiencies are the fast growing of the first replica and the occurence of small negative occupation numbers for larger times. As discussed below a satisfying description is only obtained by the additional inclusion of the renormalized phonon propagation.

It is remarkable, that already for the small coupling shown in the figures all the approximations show weight at momenta larger than k_{0} around k_{0} + nq_{B}. In contrast the exact solution never shows weight at momenta larger than k_{0}.

Fig. 8 shows the depletion of the initially excited momentum state k_{0}. The coupling strength is chosen to be ten times larger than in Figs. 5-7. Apart from the Boltzmann approximation with γ_{k} = 0 all the approximations agree with the exact solution for small times (v_{F} q_{c})t < 10. In the Boltzmann approximation we find the typical exponential decay at small times. The best agreement is given by the Markovian approximation with γ_{k} = 0 and the non-Markovian approximation with γ_{k} ≠ 0. The Markovian approximation with γ_{k} ≠ 0 falls off too slowly at intermediate times.

A comparison of the exact electron distribution and the approximate one for the Fröhlich coupling and not too strong couplings leads to similar results as described above. We also do not find a systematic improvement going from the simple Markovian approximation to the more complicated non-Markovian approximation. For intermediate values of the coupling strength the use of the perturbative value for γ_{k} can lead to negative occupation numbers. Because of the damping the influence of these unphysical negative values on the distribution for larger times is not as bad as for the non-Markovian approximation with γ_{k} = 0. Also for γ_{k}-values different from the perturbative one we found such unphysical negative values.

In Fig. 9 we show the depletion of the occupation at k_{0} for the Fröhlich coupling. In contrast to the constant coupling case the decay is superposed by an oscillation with frequency ω_{0}. The reason for this oscillation is the same as for the ω_{0}-oscillation in the energy discussed in Sec. II. Due to the different q → 0 prefactor behavior the weight of the ω_{0}-oscillation is enhanced compared to the constant coupling. The approximations also show this oscillatory decay.

We repeated all the approximate calculations for thermal phonons where we replaced N_{q}(t) in Eq. (49) by the Bose distribution for the given temperature. Because the phonon numbers in the consistent calculations shown in the above figures are very small, for most of the momenta the electron distribution is very similar and there are only two essential differences. At T = 0 because of N_{q}^{th} ≡ 0 there is no weight at momenta larger than k_{0} and the Fermi sea is inert.

The non-Markovian approximation using electronic quasi-particles but unperturbed phonon propagation is not overall better than the corresponding Markovian version, which is numerically much easier to handle. This holds for the time dependent and the thermal description of the phonon distribution. In the next subsection we discuss the importance of the additional renormalization of the phonon propagation.

B. Renormalized phonon propagation

Using degenerate perturbation theory it follows for a general model of electrons coupled to optical phonons, that near the crossing of the phonon mode energy ω_{0} and the excitation energies of the electron system (in our model described by the boson mode ε_{q}) there is an energy splitting and the phonon energy is renormalized. In our model ω_{0} → λ_{±}(q_{B}) which leads to a gap of phonon frequencies around ω_{0}. This shows that the simple Boltzmann picture of energy and momentum conserving collisions of noninteracting electrons and phonons cannot be correct. If the renormalization of the phonon frequency is not taken into account, as in the approximations discussed in the preceding subsection, too much weight is produced near the momenta k_{0} − nq_{B} for larger times.

In Fig. 10 we compare the exact momentum distribution with the one calculated with the non-Markovian approximations with unperturbed equilibrium phonon propagation and with the equilibrium phonon propagation given in Eq. (48) for the same parameters as in Fig. 7. The agreement with the exact result is drastically improved using the renormalized phonon propagation. Especially the first replica, which is too large using the unrenormalized propagation, fits very well. For smaller times a similar improvement is obtained. It also holds for times much larger as shown Figs. 11 and 12. The tendency of distributing the spectral weight more evenly around all momenta for larger times is incorporated in the kinetic equations only if the renormalized phonon propagation is used. Surprisingly this approximation describes well even the detailed features around k_{0} − q_{B} of the exact solution. The non-Markovian approximation with renormalized phonon propagation shows an overall good agreement with the exact result also for larger electron-phonon coupling.

The same improvement occurs if the Fröhlich interaction is used. To demonstrate this we compare in Fig. 13 the depletion of the initially excited momentum state for the exact solution and the two non-Markovian approximations for the same parameters as in Fig. 9.

As for the approximations discussed in subsection IV A we repeated all calculations with a thermal phonon distribution instead of the time dependent one. The differences between both treatments are again minor.

Using the renormalized phonon propagation also improves the quality of the Markovian version of the kinetic equations, but not as drastically as for the non-Markovian approximation. Therefore also for more realistic models non-Markovian equations with renormalized electron and phonon propagation should be used, if a partially filled band is present.

V. SUMMARY

We have presented a model which allows to study exactly the relaxation of electrons from a nonequilibrium configuration due to the interaction with each other and with phonons. This microscopic model shows real electron relaxation with qualitative features known from approximate calculations in more realistic models. We have compared the exact results for the distribution functions with the results from approximate kinetic equations.

The basic ingredients in the derivation of the kinetic equations are the GKBA and the replacement of nonequilibrium retarded and advanced Green’s functions by their equilibrium versions. As another approximation in the derivation we neglected the vertex correction in the nonequilibrium selfenergy. We found that for the coupling strengths chosen this approximation leads to very good results for the electronic equilibrium selfenergy. Therefore it seems plausible that this approximation is reasonable also in the nonequilibrium case. The approximation of the retarded and advanced equilibrium electron Green’s functions by damped quasi-particles was found to be very good for the relevant momenta and the weak couplings chosen. Using the second order selfenergy for the equilibrium phonon propagation already leads to the exact result in our model. The uncontrolled approximations in the derivation of the kinetic equations (49) and (50) are therefore the GKBA and the subsequent replacement of the retarded and advanced nonequilibrium Green’s functions by the equilibrium ones.

The additional approximation to use the unperturbed phonon propagation after the GKBA has been introduced was shown to be responsible for the shortcomings of the corresponding approximate non-Markovian quantum kinetic equation. Only the additional renormalization of the phonon propagation led to improved kinetic equations that allow to describe in a satisfactory way even detailed features of the exact solution. We expect an improved description of the phonon propagation to be important also in more realistic models, if a partially filled band is present.

ACKNOWLEDGMENTS

Three of the authors (V. M., C. W., and J. F. ) are grateful to the Deutsche Forschungsgemeinschaft (SFB 345 “Festkörper weit weg vom Gleichgewicht”) for financial support.

APPENDIX A

In this Appendix we present the derivation of Eq. (30). By inserting the inverse of Eq.

(15) into Eq. (18) we obtain

f(k, k_{0}, t) = \frac{1}{L^{2}} ∫_{−L/2}^{L/2} dx_{1}dx_{2}dx_{3}dx_{4}e^{−ik_{0}(x_{1}−x_{4})}e^{−ik(x_{3}−x_{2})}I(x_{1}, x_{2}, x_{3}, x_{4}, t)  (A1)

with

I(x_{1}, x_{2}, x_{3}, x_{4}, t) :=< ϕ_{0}|ψ(x_{1})ψ^{†}(x_{2}, t)ψ(x_{3}, t)ψ^{†}(x_{4})|ϕ_{0} > .  (A2)

Using the bosonization of the field operators Eq. (16) this leads to

I(x_{1}, x_{2}, x_{3}, x_{4}, t) = \frac{1}{L^{2}} ⟨0 |exp \{− ∑_{q>0} \sqrt{\frac{2π}{Lq}} e^{−iqx_{1}}b_{q}^{†}\} exp \{∑_{q>0} \sqrt{\frac{2π}{Lq}} e^{iqx_{1}}b_{q}\}

×e^{iHt} exp \{∑_{q>0} \sqrt{\frac{2π}{Lq}} e^{−iqx_{2}}b_{q}^{†}\} exp \{− ∑_{q>0} \sqrt{\frac{2π}{Lq}} e^{iqx_{2}}b_{q}\}e^{−iHt}

×e^{iHt} exp \{− ∑_{q>0} \sqrt{\frac{2π}{Lq}} e^{−iqx_{3}}b_{q}^{†}\} exp \{∑_{q>0} \sqrt{\frac{2π}{Lq}} e^{iqx_{3}}b_{q}\}e^{−iHt}

× exp \{∑_{q>0} \sqrt{\frac{2π}{Lq}} e^{−iqx_{4}}b_{q}^{†}\} exp \{− ∑_{q>0} \sqrt{\frac{2π}{Lq}} e^{iqx_{4}}b_{q}\}| 0⟩ .  (A3)

We now express the b_{q} (b_{q}^{†}) in terms of α_{q} (α_{q}^{†}) and β_{q} (β_{q}^{†}) and normal order afterwards. By using the Baker-Hausdorff formula several times we end up with

I(x_{1}, x_{2}, x_{3}, x_{4}, t) = \frac{1}{L^{2}} e^{F (x_{1},x_{2},x_{3},x_{4},t)} ,  (A4)

where

F (x_{1}, x_{2}, x_{3}, x_{4}, t) = \frac{2π}{L} ∑_{q>0} \frac{1}{q} \{e^{iq(x_{1}−x_{4}+i0)} + e^{iq(x_{2}−x_{3}+i0)}
+ (e^{iq(x_{3}−x_{4}+i0)} − e^{iq(x_{2}−x_{4}+i0)}) (c_{q}^{2}e^{−iλ_{+}(q)t} + s_{q}^{2}e^{−iλ_{−}(q)t})
+ (e^{iq(x_{1}−x_{2}+i0)} − e^{iq(x_{1}−x_{3}+i0)}) (c_{q}^{2}e^{iλ_{+}(q)t} + s_{q}^{2}e^{iλ_{−}(q)t})} .  (A5)

We now assume that the electron-electron and the electron-phonon interaction are cutoff at q_{c} = \frac{2π}{L}n_{c}. Therefore for q > q_{c} we have s_{q} = 0, c_{q} = 1, and λ_{+}(q) = v_{F} q. Using this and the relation

exp \{ ∑_{m=1}^{∞} \frac{1}{m}e^{i \frac{2π}{L} m(x+i0)}\} = (1 − e^{i \frac{2π}{L} (x+i0)})^{−1}  (A6)

we obtain

I(x_{1}, x_{2}, x_{3}, x_{4}, t) = \frac{1}{L^{2}} (1 − e^{i \frac{2π}{L} (x_{2}−x_{4}−v_{F} t+i0)}) (1 − e^{i \frac{2π}{L} (x_{1}−x_{3}+v_{F} t+i0)})

× (1 − e^{i \frac{2π}{L} (x_{3}−x_{4}−v_{F} t+i0)})^{−1} (1 − e^{i \frac{2π}{L} (x_{1}−x_{2}+v_{F} t+i0)})^{−1}

× (1 − e^{i \frac{2π}{L} (x_{2}−x_{3}+i0)})^{−1} (1 − e^{i \frac{2π}{L} (x_{1}−x_{4}+i0)})^{−1}

× ∏_{m=1}^{nc} exp \{\frac{1}{m} [(e^{i \frac{2π}{L} m(x_{3}−x_{4}+i0)} − e^{i \frac{2π}{L} m(x_{2}−x_{4}+i0)}) u_{m}(t)

+ (e^{i \frac{2π}{L} m(x_{1}−x_{2}+i0)} − e^{i \frac{2π}{L} m(x_{1}−x_{3}+i0)}) u_{m}^{∗}(t)]\} ,  (A7)

with u_{m}(t) defined in Eq. (33). By expanding the exponential functions and the denominators in power series in e^{i \frac{2π}{L} (x_{k}−x_{l})} and ordering the powers we can write

I(x_{1}, x_{2}, x_{3}, x_{4}, t) = \frac{1}{L^{2}} ∑_{j,m,p=0}^{∞} ∑_{l,r,s=0}^{∞} e^{i \frac{2π}{L} j(x_{2}−x_{3})}e^{i \frac{2π}{L} m(x_{1}−x_{4})}e^{i \frac{2π}{L} p(x_{3}−x_{4})}
×e^{i \frac{2π}{L} r(x_{1}−x_{2})}e^{i \frac{2π}{L} l(x_{2}−x_{4})}e^{i \frac{2π}{L} s(x_{1}−x_{3})}
×a_{p}^{(n_{c})}(t)(a_{r}^{(n_{c})}(t))^{∗}b_{l}^{(n_{c})}(t)(b_{s}^{(n_{c})}(t))^{∗} ,  (A8)

with the coefficients a_{m}^{(n_{c})}(t) and b_{m}^{(n_{c})}(t) determined iteratively as given by (31) and (32). The integrations over x_{1}, x_{2}, x_{3}, and x_{4} in (A1) can then be performed and lead to Kronecker deltas which reduces the number of summation variables from six to three. This leads to Eq. (30).

REFERENCES

[1] see e. g. R. Zimmermann, Phys. Stat. Sol. 159, 317 (1990).

[2] P. Lipavský et al., Phys. Rev. B 34, 6933 (1986).

[3] H. Haug, Phys. Stat. Sol. 173, 139 (1992); L. Bányai et al., Phys. Stat. Sol. 173, 149 (1992); D. B. Tran Thoai and H. Haug, Phys. Rev. B 47, 3574 (1993).

[4] H. Schoeller, Ann. Phys. (NY) 229, 273 (1994).

[5] R. Zimmermann and J. Wauer, J. Lum. 58, 271 (1994).

[6] S. Tomonaga, Prog. Theor. Phys. 5, 544 (1950); J.M. Luttinger, J. Math. Phys. 4, 1154 (1963).

[7] F. D. M. Haldane, J. Phys. C 14, 2585 (1981).

[8] G. Wentzel, Phys. Rev. 83, 168 (1951).

[9] J. Schilp, T. Kuhn, and G. Mahler, Phys. Rev. B 50, 5435 (1994).

[10] V. Meden and K. Schönhammer, Phys. Rev. B 46, 15753 (1992); K. Schönhammer and V. Meden, Phys. Rev. B 47, 16205 (1993) and Phys. Rev. B 48, 11521 (1993).

[11] J. Voit, Phys. Rev. B 47, 6740 (1993).

[12] C. Wöhler, V. Meden, and K. Schönhammer, unpublished.

[13] J. Sólyom, Adv. Phys. 28, 201 (1979).

[14] V. Meden, K. Schönhammer, and O. Gunnarsson, Phys. Rev. B 50, 11179 (1994).

[15] D. C. Mattis and E. H. Lieb, J. Math. Phys. 6, 304 (1965).

[16] R. de L. Kronig, Physica 2, 968 (1935).

[17] D. Loss and Th. Martin, Phys. Rev. B 50, 12160 (1994).

[18] A. B. Migdal, Zh. Eksp. Teor. Fiz. 34, 1438 (1958) [Sov. Phys. JETP 7, 996 (1958)].

[19] S. Engelsberg and J. R. Schrieffer, Phys. Rev. 131, 993 (1963).

[20] M. Apostol and I. Baldea, Phys. Lett. A 88, 73 (1982).

FIGURES

FIG. 1. Evolution of the distribution of the excitation energy for L → ∞. The solid line presents the excitation energy in the electronic system, the dotted line the excitation energy in the phonon system, the dashed line the expectation value of the coupling part of the Hamiltonian, and the dashed-dotted line the constant sum of the three contributions. The parameters are Γ = 0.01, ω_{0}/(v_{F} q_{c}) = 0.2, and k_{0}/q_{c} = 1.5.

FIG. 2. Phonon distribution function h for Γ = 0.001, ω_{0}/(v_{F} q_{c}) = 0.2, and k_{0}/q_{c} = 1.5 for a finite system with ν ≡ 2π/(q_{c}L) = 1/200.

FIG. 3. Electron distribution function f for Γ = 0.001, ω_{0}/(v_{F} q_{c}) = 0.2, and k_{0}/q_{c} = 1.5 for a finite system with ν = 1/80. The occupation at k_{0}/q_{c} is scaled by a factor of 1/5.

FIG. 4. Electron spectral function for a finite system as a function of ω. The solid line presents the exact result, the dotted line the second order selfenergy result, and the dashed line the quasi-particle approximation used in the kinetic equations. The parameters are as in Fig. 3 and the momentum is chosen to be k/q_{c} = 1.5. The artificial broadening (see text) is 0.003 in units of (v_{F} q_{c}).

FIG. 5. Comparison of the exact and approximate electron distribution function for the constant coupling and a given time t. The parameters are as in Fig. 3 and (v_{F} q_{c})t = 60. The solid line presents the exact result, the dotted line the Markovian approximation with unrenormalized electron and phonon propagation, the dashed line the Markovian approximation, and the dashed-dotted line the non-Markovian approximation both with renormalized electron but unrenormalized phonon propagation. The occupation at k_{0}/q_{c} is scaled by a factor of 1/10.

FIG. 6. The same as in Fig. 5, but for (v_{F} q_{c})t = 120.

FIG. 7. The same as in Fig. 5, but for (v_{F} q_{c})t = 200.

FIG. 8. Comparison of the exact and approximate depletion of the initially excited momentum state for the constant coupling. The solid line presents the exact result, the dotted lines the Markovian and the Boltzmann approximations with unrenormalized electron and phonon propagation, the dashed line the Markovian approximation, and the dashed-dotted line the non-Markovian approximation both with renormalized electron but unrenormalized phonon propagation. The parameters are: Γ = 0.01, ω_{0}/(v_{F} q_{c}) = 0.2, k_{0}/q_{c} = 1.5, and ν = 1/80.

FIG. 9. Comparison of the exact and approximate depletion of the initially excited momentum state for the Fröhlich coupling. The solid line presents the exact result, the dotted line the Markovian approximation with unrenormalized electron and phonon propagation, the dashed line the Markovian approximation, and the dashed-dotted line the non-Markovian approximation both with renormalized electron but unrenormalized phonon propagation. The parameters are: g^{2}/(v_{F}^{2} q_{c}^{3}) = 0.000019, q_{0}/q_{c} = 0.0215, ω_{0}/(v_{F} q_{c}) = 0.2, k_{0}/q_{c} = 1.5, and ν = 1/80.

FIG. 10. Comparison of the exact and approximate electron distribution function for the constant coupling. The parameters are as in Fig. 3 and (v_{F} q_{c})t = 200. The solid line presents the exact result, the dotted line the non-Markovian approximation with renormalized electron and renormalized phonon propagation (RPP), and the dashed line the non-Markovian approximation with renormalized electron but unrenormalized phonon propagation (UPP). The occupation at k_{0}/q_{c} is scaled by a factor of 2/5.

FIG. 11. The same as in Fig. 10, but for (v_{F} q_{c})t = 300.

FIG. 12. The same as in Fig. 10, but for (v_{F} q_{c})t = 400.

FIG. 13. Comparison of the exact and approximate depletion of the initially excited momentum state for the Fröhlich coupling. The parameters are as in Fig. 9. The solid line presents the exact result, the dotted line the non-Markovian approximation with renormalized electron and renormalized phonon propagation (RPP), and the dashed line the non-Markovian approximation with renormalized electron but unrenormalized phonon propagation (UPP).