%	$Id: FinalLaser.tex,v 1.5 1998/01/11 21:40:23 lvov Exp lvov $	
\documentstyle[12pt]{article}
\oddsidemargin 1cm
\evensidemargin 1cm
\textwidth 16cm
\topmargin 0pt
\headheight 0pt
\headsep 0pt
\baselineskip 40pt
\newcounter{newcnt}
\renewcommand{\theequation}{\thesection.\arabic{equation}}
\begin{document}
\def \F {\Phi}
\def \n  {^{(0)}}
\def \Im {\rm Im}
\def \a {\alpha}
\def \b {\beta}
\def \e {\epsilon}
\def \o {\omega}
\def \O {\Omega}
\def \d {\delta}
\def \D {\Delta}
\def \k {\kappa}
\def \g {\gamma} 
\def \G {\Gamma}  
\def \s {\sigma}
\def \S {\Sigma}
\def \l {\lambda}
\def \L {\Lambda}
\def \t {\tau}
\def \p {\pi}
\def \m {\mu}
\def \BE {\begin{equation}}
\def \EE {\end{equation}}
\def \BEA {\begin{eqnarray}}
\def \EEA {\end{eqnarray}}
\def \CR {\nonumber \\}
\def \P {{{\cal{P}}}}
\def \Q {{{{\cal Q}}}}
\def \R {{{\cal{R}}}}
\def \sgn {{\rm{sgn}}}
\def \r{ \rho}
\def \NN{\CR}
\def \bbox {\bf }
\begin{center}
        { \Large Quantum Weak Turbulence with Applications \\  to 
        Semiconductor Lasers.}
\vskip 2cm  {Y.V. Lvov$^{1,2}$, R.  Binder$^3$ and 
        A.C. Newell$^{1,4}$ }\vskip 1cm 
{\bf Physica D, {\bf 121}, pp. 317 - 343,(1998).}
\end{center}
{\small

$^1$ Department of Mathematics, The University of Arizona, Tucson,
85721 Arizona

$^2$ Department of Physics,     The University of Arizona, Tucson, 85721 Arizona

$^3$ Optical Sciences Center, The University of Arizona,
Tucson, 85721 Arizona 

$^4$ Department of Mathematics, University of Warwick, Coventry,
CV47AL, UK 
}
\abstract { Based on a model Hamiltonian appropriate for
the description of fermionic systems such as semiconductor lasers, we
describe a natural asymptotic closure of the BBGKY hierarchy in
complete analogy with that derived for classical weak turbulence. The
main features of the interaction Hamiltonian are the inclusion of full
Fermi statistics containing Pauli blocking and a simple,
phenomenological, uniformly weak two particle interaction potential
equivalent to the static screening approximation. We find a new class
of solutions to the quantum kinetic equation which are analogous to
the Kolmogorov spectra of hydrodynamics and classical weak
turbulence. They involve finite fluxes of particles and energy {in
momentum space} and are particularly relevant for describing the
behavior of systems containing sources and sinks. We make a prima
facie case that these finite flux solutions can be important in the
context of semiconductor lasers and show how they might be used to
enhance laser performance.}



\section{Introduction and General Discussion.} 

%{\huge ack grant AFOSR 94-1-0144DEF Final Report due on March 15 
%            Close to Final Copy Feb 15 }
\setcounter{equation}{0}
\subsection {Background.}
Gaining a good understanding of the relaxation processes of many
particle (or many wave) systems is a difficult task. A successful
approach requires many ingredients. First, one must make
approximations concerning the relative strength and uniformity (in
momentum space) of the nonlinear coupling. Second, one needs to
understand how the infinite hierarchy of moment equations can
effectively decouple and give rise to a closed kinetic Boltzmann
equation for the redistribution of particles (waves) and energy in
momentum space.  One of the goals of this paper is to show that one
can indeed, as one can do for classical systems, derive the quantum
kinetic Boltzmann equation in a self consistent manner without
resorting to a priori statistical hypotheses or cumulant {discard}
assumptions.  Third, systems are rarely in isolation and frequently
involve sources (some instability or forcing that injects particles
and energy, often at specific locations in momentum space) and sinks
(regions of absorption of particles and energy).  Moreover, and in
analogy with hydrodynamic turbulence and classical wave systems, the
presence of sources and sinks dramatically changes the nature of the
equilibria reached by the quantum kinetic equation. A second and major
goal of this paper is to describe these equilibria and  show how 
they can affect the
output of semiconductor lasers.

\input epsf




Fermionic quantum systems have been studied intensively for more than
five decades.  Many theoretical approaches (see, for example,
\cite{b1}-\cite{kuznetsov}) have been developed but the common feature
of all is the derivation of the quantum mechanical Boltzmann kinetic
equation (henceforth referred to as the quantum kinetic equation or
QKE) which describes the evolution and relaxation of the particle
number, $<a_k^\dagger a_k>$, due to collisions.  It is the analogue of
the classical kinetic equation (KE) of wave turbulence (see, e.g., \cite{ZLF}) and
its natural quantum extension to boson gases.
The QKE accounts for phase space blocking effects (Pauli exclusion
principle) and is equivalent to what one would obtain by treating the
scattering cross-section using the second Born approximation.
Examples of variations of this equation and its derivation include the
Lenard-Balescu equation which accounts for dynamical screening effects
(see, e.g., \cite{b7}-\cite{b8}).  Other extensions include
generalized scattering cross sections, such as the exchange effects
(crossed diagrams), as well as T-matrix effects (see, e.g.,
\cite{b2},\cite{b9}).  T-matrix approaches are especially important in
systems that allow for bound states.  Further generalizations of the
quantum Boltzmann equation are the various forms of the Kadanoff-Baym
equations. The most general Kadanoff-Baym equations are two-time
equations which describe charge-carrier correlations consistently with
relaxation dynamics.  The Markov approximation, which is the
lowest-order gradient expansion with respect to macroscopic times,
yields the familiar form of the Lenard-Balescu equation, and the
additional static screening and small momentum transfer approximation
yields the Landau kinetic equation.  In contrast to the slow
relaxation dynamics for which the Kadanoff-Baym gradient expansion is
applicable, recent investigations have addressed the issue of
ultrafast relaxation and the related problems of memory effects as
contained in the Kadanoff-Baym equations (see, e.g., \cite{koehler}
and all references therein) and in the generalized Kadanoff-Baym
equations of Lipavski, Spicka, and Velicky (see, e.g.,
\cite{Lipavsky}-\cite{kuznetsov}).

In all of the derivations, the principal obstacle to be overcome is
the closure problem. Due to the nonlinear character of the quantum
mechanical Coulomb interaction Hamiltonian, the time evolution of the
expectation values of two operator products such as $<a_k^\dagger
a_k>$ are determined by four operator expectation values
$<{a^\dagger_{k_1}a^\dagger_{k_2}a_{k_3}a_{k_4}}>$.  The problem
compounds. The time evolution of the N operator product expectation
value is determined by (N+2) operator expectation values and one is
left with an infinite set of moment operator equations known as the
BBGKY hierarchy. The closure problem is to find a self consistent
approximation of this infinite hierarchy which reduces the infinite
set of coupled equations to an infinite set of equations which are
essentially decoupled. To do this, one needs to make approximations
and these almost always involve the introduction of small parameters.

In what follows we introduce one such small parameter, the relative
strength $\e, 0<\e\ll 1$, of the coupling coefficient $T_{k k_1,k_2
k_3}$ in the system Hamiltonian ($V=$ d-dimensional system volume)
\BEA H  &=&  \frac{V}{(2\pi)^d}
\int d {\bf k} \hbar \o_k
a_k^\dagger a_k  \CR &&+ 
 \frac{V^3}{(2 \pi)^{3d}}
\frac{1}{2} \int d {\bf k_1} d{\bf k_2} d {\bf k_3}
d {\bf k_4} T_{k k_1, k_2 k_3}a^\dagger_{k} a^\dagger_{k_1} a_{k_2}
a_{k_3} \delta ({\bbox k}+{\bbox k_1} -{\bbox k_2} -{\bbox k_3}) .\NN
\label{Ham}\EEA By relative, we mean the ratio of $T_{k k_1,k_2 k_3}$
to $\o_k$.  In other words, the linear response of the Hamiltonian is
dominant over short times.  This is not a trivial approximation
because in order to develop a consistent asymptotic closure we need
$T_{k k_1, k_2 k_3} $ to be small compared to $\o_k $ {\it uniformly}
in $k$. There are many important examples where this is not the
case. For instance, for the classical nonlinear Schroedinger system
(see, e.g., \cite{Optical}), which describes the weak turbulence of
optical waves of diffraction in a nonlinear medium, $\o_k\propto k^2$
and $T_{k k_1, k_2 k_3}=a/(4 \p^2)$ where $a$ is a positive
(negative) constant for the focusing (defocusing) case and it is clear
that the ratio of $T_{ k k_1, k_2 k_3}$ to $\o_k$ increases
dramatically as $k\to 0 $. This has real physical consequences. The
asymptotic closure for the classical case, analogous to one we are
about to derive for fermions, whose central feature is the classical
kinetic equation describing resonant four-wave energy exchange, is
only valid when the particles and energy are at finite $k$ values. But
the dynamics of the classical kinetic equation are such that most of
the energy travels to a sink at large $k$ and most of the particles
and some energy travels to low $k$. Unless there is a strong damping
near ${\bbox k}=0$, the flow of particles and energy towards ${\bbox
k}=0$ will trigger collapsing filaments in the focusing case or build
condensates in the defocusing case. In either case, unless there is
strong damping near ${\bbox k}=0$, the weak turbulence theory of four
wave interactions breaks down. Likewise in plasmas and and highly
excited semiconductors, while the frequency $\o_k$ is $\a k^2$, the
bare coupling coefficient $T_{ k k_1, k_2 k_3}$ behaves as $1/k^2$
reflecting Coulomb interactions. Again, even though the carriers may
be initially excited at ${\bbox k}$ values much greater than zero, the
natural dynamics of either the classical or quantum kinetic equations
will lead to particle deposit at small ${\bbox k}$. However, in both
of these cases, there is another effect which modifies the
$T_{kk_1,k_2k_3}$ and brings the theory closer to the nonlinear
Schroedinger case.  The physical origin of this effect is screening
and it is manifested as a weakened potential at long distances or
small ${\bbox k}$ values. The modification is a renormalization of the
coupling coefficient $T_{kk_1,k_2k_3}$ near ${\bbox k}=0$ which
effectively cancels the singularity, or equivalently replaces $k^{-2}$
by $(k^2+\kappa^2)^{-1}$, where $\kappa$ is the inverse screening
length. This is the small $k$ limit of the dynamically screened
effective interaction in the Lenard-Balescu approach (see, e.g.,
\cite{b7}). In this paper, we do not take issue \footnote{ The strong
analogy with the nonlinear Schroedinger equation would suggest that
the question of the self consistency of the quantum weak turbulence
theory without small $k$ damping might be revisited.} with the self
consistency of this approximation because we will be introducing a
sink at small $k$ values (greater than $\kappa$) which remove the
potential danger in the nonuniformity of the ratio of $T_{k k_1,k_2
k_3}$ to $\o_k$. This is not at all artificial. The sink is, in the
semiconductor laser context, nothing other than the lasing output.  In
what follows, then, we assume that in the $k$ regions of interest, the
ratio of $T_{k k_1,k_2 k_3}$ to $\o_k$ is uniformly small.

There are essentially two reasons for the successful closure of the
hierarchy over long times (asymptotic closure).  First, to leading
order in $\e$, the cumulants (the non Gaussian part) corresponding to
the expectation values of products of $N$ operators ($N>2$) play no
role in the long time behavior of cumulants of order $r<N$. This is a
consequence of phase mixing and due to the nondegenerate nature of the
dispersion relation $\o_k=\o(|{\bbox k}|)$. In fact, we will find that
these (zeroth order in $\e$) cumulants slowly decay. The net result is
that the initial state (subject to certain smoothness conditions)
plays no role in the long time dynamics. Second, and more important,
the higher order (in $\e$) contributions to these $N$-operator
cumulants, which are generated directly by the nonlinearities in the
system, are dominated by products of lower order cumulants. Some of
these dominant terms are supported only on certain well defined low
dimensional (resonant) manifolds in momentum space. These are
responsible for the redistribution of particle number among the
different momentum states and for the slow decay in time of the
leading order approximation to the cumulants of order $N>2$.  The
other dominant terms are nonlocal and lead to a nonlinear frequency
renormalization of $\o_k$.

As a result, we are led to a very simple and natural asymptotic closure of
the BBGKY hierarchy.

\subsection{Brief discussion of
principal results} Here we list and briefly discuss the main and new
results of this article.  \begin{itemize}\item {\it 1.} We
systematically derive evolution equations over long times $\e^{-2}$,
of the order of the inverse square of the coupling strength, for the
{\it leading order} approximations of {\it all} $N^{th}$ order
cumulants $N\ge 2$. The first of these is the QKE, a closed equation
for the particle number, and is consistent with what one obtains using
the Hartree-Fock approximation
\footnote{ Here, we use the term Hartree-Fock approximation
 in the commonly
used sense, which is not identical to the strict meaning
of a quantum system with a Hartree-Fock wavefunction. In the
latter, scattering and correlations are excluded.}
{ when one takes proper account of
the frequency renormalization {\it to order $\e^2$.}} The remaining
equations for the leading order behavior of the higher order cumulants
decouple and can be solved by a simple renormalization of the
frequency which depends only on particle number density. The first
correction is the well known Hartree-Fock self energy correction
whereas the second correction is new and moreover complex. The sign of
the imaginary part of the latter is positive definite as the particle
number approaches an equilibrium solution of the QKE. As a result, the
leading order contribution to the cumulants of order $N>2$ (which
contain information on the initial state) slowly decays.  Not so the
higher order corrections.  We obtain specific expressions for these and
show that they contain terms with cumulants of higher order, products
of cumulants of order $N$ with products of particle number and
products simply containing particle number alone. Only the latter two
types of terms exert a long time influence and it is this fact that
leads us to a natural asymptotic closure of the hierarchy.
\item{\it 2.} The connections with classical systems and boson systems
are discussed. By analogy with results obtained in the case of the
classical kinetic equation, we introduce a new four parameter family
of steady (equilibrium) solutions for the QKE. The four parameters are
temperature $T$, chemical potential $\mu$, particle number flow rate
$Q$ and energy density flow rate $P$.  When $P=Q=0$, the equilibrium
solution is the well known Fermi-Dirac distribution. For either $P$ or
$Q$ nonzero, the distributions are the analogues of the Kolmogorov
distributions of classical wave turbulence, in which particles
(energy) flow from sources at intermediate momentum scales to sinks at
low (high) momentum values.

These solutions have not been considered in the fermionic context
before. It turns out in fact, in contrast to the pure Kolmogorov
solutions (for which $Q=0,\ \ P>0$, or $P=0, \ \ Q>0$) that the
relevant equilibrium solutions of the QKE which describe (a) the
finite flux of particles and a little energy between intermediate
momenta at which the system is pumped and a sink at low momenta and
(b) the finite flux of energy and a few particles between the
intermediate momenta and the energy dissipation sink at high momenta,
involve special relations between $Q$ and $P$.  These solutions are
also new for the cases of classical wave turbulence and boson systems.
\item{\it 3.} Again by analogy with the classical case, we introduce a
simple, local approximation to the collision integral in the QKE which
gives excellent qualitative results. This is known as the differential
approximation (see, e.g., \cite{Optical},\cite{HASS},\cite{Balk}) and is strictly
valid only when the principal transfer of energy and particles is
between close neighbors in momentum space. Its most important feature,
however, is that it allows us explore both analytically and
numerically the behavior of the collision integral and thereby gain a
clear qualitative picture about the nature of its steady (equilibrium)
solutions.
\item {\it 4.}  We briefly examine the relevance of these results in
 the context of semiconductor lasers using a simplified model which
 assumes that any relaxation due to carrier-phonon interactions is
 negligible with respect to carrier-carrier interactions. We suggest
 that, because of Pauli blocking which results in a reduced pumping
 efficiency into those momentum states that are already filled, it may
 be advantageous to pump these lasers in a fairly narrow momentum
 window whose corresponding frequency is greater than the lasing
 frequency and allow the finite flux steady (equilibrium) solutions to
 carry particles and energy back to the lasing frequency via collision
 processes.  Indeed we show that, at least for the model we use, this
 strategy is more efficient than pumping the system across a wide
 range of momenta as it is currently done.  Indeed, in several
 {numerical} experiments, for equal amounts of energy pumping, we find
 that the laser turns on when the finite flux equilibrium is excited
 but fails to turn on under broadband pumping. When the pumping is
 increased so that the laser turns on under the latter conditions, the
 output power is significantly smaller than that emitted when the
 finite flux solution is operative.
\end{itemize}

\section{Systematic Derivation of the 
Kinetic Equation and Evolution equations for Higher Order Cumulants}
\setcounter{equation}{0}
\subsection{Basic Definitions and Evolution Equation}
We start from the Hamiltonian (\ref{Ham}) of a spatially homogeneous
system of particles with binary interactions.  Here, $e_k=\hbar \o_k$
is the energy level of momentum state ${\bbox k}$ (for example, in
semiconductors, the parabolic band approximation is given by $\o_k =\a
{\bbox k}^2$) where ${\bbox k}$ is a d-dimensional wave vector, and
$a_k ^\dagger , a_k$ are fermionic creation/annihilation operators
fulfilling the anticommutation relations, \BE a_i a_j^\dagger +
a_j^\dagger a_i = \d_{ij} \ \ \ \ \ \ \ a_i a_j + a_j a_i
=0. \label{commutation}\EE We include the size of the system in the
definition of the interaction matrix element $T_{12,34}$.  We
introduce for convenience the short hand notation: 
${\bbox k_1}\equiv 1$, 
$\frac{V}{(2\pi)^d} d {\bbox k_1} \equiv d1$, and 
$\frac{(2\pi)^d}{V} \delta ({\bbox k_1}+{\bbox
k_2}-{\bbox k_3}- {\bbox k_4}) \equiv \d^{12}_{34}$, $\hbar=1$.  The
Hamiltonian now reads $$H= \int d1 \o_1 a_1^\dagger a_1 + \frac{1}{2}
\int d1234 T_{12,34}a^\dagger_{1} a^\dagger_{2} a_{3} a_{4}
\delta^{12}_{34}.$$ If one interchanges the indices $1$ and $2$ or $3$
and $4$ in the above expression and uses the fact that the Hamiltonian
is Hermitian, the following properties hold: $$ T_{12,34}
=-T_{21,34}=T_{21,43}=T^*_{43,21}.$$ In the Heisenberg picture, the
equations of motion are \BE\dot a_k = i[H,a]_-\label{Heizenberg}\EE
which give \BE \dot a_k^\dagger = i \o_k a_k^\dagger + i \int
T_{32,10} a^\dagger_3 a^\dagger_2 a_1 \delta_{32}^{10} d123 \EE and
\BE \dot a_k = -i \o_k a_k - i \int T_{01,23} a^\dagger_1 a_2 a_3
\delta_{32}^{10} d123.
\label{eqsmotion}\EE From the Heisenberg equations of motion, one can now 
derive the BBGKY hierarchy of equations for the normal ordered 
expectations values. The first three
are, \BEA \frac{d}{d t} \r_k &=& 2 \Im\int T_{01,23} 
N_{0123}\d_{32}^{10}d123, \\ \frac{d}{d t} N_{1' 2' 3' 4'} &=& i 
 \Delta^{1'2'}_{3'4'} N_{1'2'3'4'}+i \int d123 \left( T_{3 2, 1 1'} (N_{3 
2 3' 4'} \d^{1}_{ 2'} - N_{3 2 2' 1 3' 4'})\d^{3 2}_{ 1 1'} + \right.\CR&& 
T_{3 2, 1 2'}N_{1' 3 2 1 3' 4'}\d^{3 2}_{ 1 2'}- \CR&&\left.  T_{3' 1, 2 
3} N_{1' 2' 1 2 3 4'}\d^{3 2}_{ 1 3'} - T_{4' 1, 2 3 } (N_{1' 2' 2 
3}\d^{3'}_{ 1} - N_{1' 2' 1 3'23})\d^{3 2}_{ 1 4'}\right),\CR&& 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \  {\bbox k_1'}+{\bbox k_2'}-{\bbox 
k_3'}- {\bbox k_4'}=0,\\ \label{N6timeev}
\frac{d}{d 
t}N_{1'2'3'4'5'6'}&=&i \Delta^{1'2'3'}_{4'5'6'}N_{1'2'3'4'5'6'}\cr &+&i \int 
d123 \left[ {T_{32,11'}\left( 
\d^{1}_{2'}N_{323'4'5'6'}-\d^{1}_{3'} N_{322'4'5'6'}+ 
        N_{322'3'14'5'6'}\right)\d^{32}_{11'}} \right.\cr
&&+T_{32,12'}\left( \d^{1}_{3'} N_{1'324'5'6'}-N_{1'323'14'5'6'}
\right) \d^{32}_{12'}\cr
&&+T_{32,13'} N_{1'2'3214'5'6'}\d^{32}_{13'}-T_{4'1,23}
N_{1'2'3'1235'6'}\d^{4'1}_{23}\cr
&&-T_{5'1,23}\left( \d^{1}_{4'}
N_{1'2'3'236'}-N_{1'2'3'14'236'}\right)\d^{5'1}_{23}\cr
&&-T_{6'1,23}\left(\d^{1}_{5'}N_{1'2'3'4'23} -\d^{1}_{4'} N_{1'2'3'5'23}+
        N_{1'2'3'14'5'23}\right) \d^{6'1}_{23}\left.\right], \CR &&
\ \ \ \ \ \ \   \ \ \ \ \ \ \ \ \ 
\ \ {\bbox k_1'}+{\bbox k_2'}+{\bbox k_3'}- {\bbox k_4'}-{\bbox 
k_5'}+{\bbox k_6'}=0, \EEA where 
\BEA\Delta^{12}_{34}&\equiv&(\o_{k_1}+\o_{k_2}-\o_{k_3}-\o_{k_4}),\CR
\Delta^{123}_{456}&\equiv&(\o_{k_1}+\o_{k_2}+\o_{k_3}-\o_{k_4}-
\o_{k_5}-\o_{k_6}).
\label{deltas}
\EEA
Here, the expectation value is taken with respect to an arbitrary
initial state $\Phi$, i.e., 
\BEA 
 \r(k_1)\d^{1}_{2}&=&<\Phi|a_1^\dagger a_2|\Phi>\cr
 N_{1234}&=&   <\Phi|a_1^\dagger a_2^\dagger a_3 a_4|\Phi>\cr
 N_{123456}&=& <\Phi|a_1^\dagger a_2^\dagger a_3^\dagger a_4 a_5 a_6|\Phi> \cr
 N_{12345678}&=& <\Phi|a_1^\dagger 
 a_2^\dagger a_3^\dagger a_4^\dagger a_5a_6a_7a_8|\Phi> \label{average}\EEA 
In the definitions  of these $2 m \ (m=1,2,...) $ order 
expectation values, the first $m$ indices correspond to creation
operators and the last 
$m$ indices correspond to annihilation operators. The number of 
creation and number of annihilation operators are equal to each other  
because the Hamiltonian (\ref{Ham}) conserves number of particles. The 
fact that the right hand sides of (\ref{average}) are zero on 
${\bbox k_1}+{\bbox k_2}-{\bbox k_3}-{\bbox 
k_4}=0$,\ \ \ ${\bbox k_1}+{\bbox k_2}+{\bbox k_3}-{\bbox k_4}-{\bbox 
k_5}-{\bbox k_6}=0$ respectively is a direct consequence of the spatial 
homogeneity of the system. This means that the $N^{th}$ order moment of 
the spatially dependent  field operators $\psi^\dagger({\bbox 
x_j}), \psi({\bbox x_j})$, the generalized Fourier transforms of the 
creation and annihilation operators
$a^\dagger({\bbox k_j}),\ \  a({\bbox k_j})$, depends only on the relative 
spacing; i.e. on the differences of the coordinates ${\bbox x_j},\ 
j=1..N$.
\subsection{Cumulants and their evolution.} 
Let us define the cumulant of the $N$ product of spatially dependent
field operators $\psi({\bbox x_j}),\ \ j=1,2...N$ to be the moment of
order $N$ from which the appropriate combinations of products of lower
order moments are subtracted so that the resulting expression has the
property that it decays to zero as the separations ${\bbox x_j}-{\bbox
x_i}$ become large.  The Fourier transforms of these cumulants are
therefore well defined ordinary (as opposed to generalized) functions
and these are the objects $\r(k), \P_{1234}, ...$ with which we deal.
Moreover, if the statistics are exactly Gaussian (namely, Hartree-Fock
like), then all cumulants of order $N$, $N>2$, are zero. Because these
weakly interacting fermionic systems relax to a near Gaussian state,
the cumulants are the most convenient dependent variables.


For the two and three particle functions we define the fourth and
sixth order cumulants ${\P_{1234}}$ and $\Q_{123456}$ respectively by
\BE\label{facN4} N_{1234}= \r_1 \r_2 (\delta^{2}_{3}\delta^{1}_{4}-
\delta^{2}_{4}\delta^{1}_{3})+\P_{1234}\d^{12}_{34}\label{dN4def}\EE
and \BEA N_{123456}= \r_1 \r_2 \r_3 \cdot \left( \right.  \d^{3}_{4}
(\d^{2}_{5}\d^{1}_{6}-\d^{1}_{5}\d^{2}_{6})+ 
\d^{3}_{5}(\d^{2}_{6}\d^{1}_{4}- \d^{2}_{4}\d^{1}_{6})+ \d^{3}_{6}
(\d^{1}_{5}\d^{2}_{4}- \d^{1}_{4}\d^{2}_{5}) \left.  \right)+\CR
\r_3[+\P_{1256}\d^{3}_{4}\d^{12}_{56}-\P_{1246}\d^{3}_{5}\d^{12}_{46}
+ \P_{1245}\d^{3}_{6}\d^{12}_{45}]+\cr
\r_2[-\P_{1356}\d^{2}_{4}\d^{13}_{56}+\P_{1346}\d^{2}_{5}\d^{13}_{46}
- \P_{1345}\d^{2}_{6}\d^{13}_{45}]+\cr
\r_1[+\P_{2356}\d^{1}_{4}\d^{23}_{56}-\P_{2346}\d^{1}_{5}\d^{23}_{46}
+ \P_{2345}\d^{1}_{6}\d^{23}_{45}]\cr+ \Q_{123456}\d^{123}_{456}
\NN
\label{dN6def} \EEA The expressions are analogues 
to what one would obtain in the classical case.  There, the symbols
$a^\dagger$ and $a$ are complex numbers (or {\it c-numbers}, as
opposed to operators in the quantum case) and one defines the
cumulants by expanding the fourth order moments into all possible
decompositions; namely \BE \left<a_1^*a_2^* a_3
a_4\right>=\{a_1^*a_3\}\{a_2^*a_4\}\d^1_3\d^2_4+
\{a_1^*a_4\}\{a_2^*a_3\}\d^1_4\d^2_3+
\{a_1^*a_2^*a_3a_4\}\d^{12}_{34}\EE where angle brackets denote
moments and curly brackets denote the corresponding cumulants.  In the
classical case, it is consistent, but not necessary (one can keep the
other terms and discover they play no essential role), to set all
correlations such as $\{a^*_1a^*_2\}$ and $\{a_3a_4\}$ equal to zero.
In the quantum case, one also decomposes $<\Phi|a_1^\dagger
a_2^\dagger a_3 a_4|\Phi>$ into products of all possible
decompositions. Again, it is consistent, but not necessary, to set 
terms such as $\{a^\dagger_1a^\dagger_2\}$ and $\{a_3a_4\}$ equal
to zero. The resulting decomposition should be consistent with the
anticommutation relations ({\ref{commutation}}), from which follows
that $N_{1234}=-N_{2134}$. Therefore, certain terms (for example, 
$\r_1\r_2\delta^{2}_{4}\delta^{1}_{3}$) are  negative in
({\ref{dN4def}}-\ref{dN6def}).

A general algorithm for the decomposition of the $N^{\rm th}$ order
expectation value is given in the Appendix.

Having defined the higher order cumulants, we can now write down the
evolution equations for the cumulant hierarchy. For the purpose of
deriving the QKE, it is sufficient to consider only the equations for
$\r_k$ and $\P_{1234}$. To obtain the frequency corrections to order
$\e^2$ ($\e, 0<\e\ll1$, is a measure of the strength of the coupling
coefficient), we need to consider contributions coming from the
equation for $\Q_{123456}$.  In carrying out the analysis on $\r_k$,
$\P_{1234}$ and $\Q_{123456}$, one finds, just as in the classical
case, that certain patterns emerge which allow one to identify the
terms in the equations for the cumulants of arbitrary high order that
gives rise to long time effects. Taking account of these terms gives
the expansions {(\ref{oddvaa},\ref{oddvab}) } which will be discussed
in the next section. In this section we only write down equations for
$\r_k$ and $\P_{1234}$. They are \BE \frac{d}{dt}\r_{k_0}= 2 \Im\int
T_{01,23}\P_{0123}\d^{01}_{23}d123 ,\label{rotimeevol}\EE where ${\rm
Im}$ denotes the imaginary part and the symbol zero denotes ${\bbox
k}$, and \BEA \frac{d}{dt}\P_{1'2'3'4'}&=& i
\tilde\Delta^{1'2'}_{3'4'} \P_{1'2'3'4'}- (\dot \r_{1'}\r_{2'} +
\r_{1'}\dot \r_{2'})
(\d^{2'}_{3'}\d^{1'}_{4'}-\d^{1'}_{3'}\d^{2'}_{4'}) \cr &+&
\underline{2 i T_{4'3',2'1'}\r_{3'}\r_{4'}(1-\r_{2'}-\r_{1'})+ 2 i
T_{4'3',2'1'}\r_{1'}\r_{2'}(\r_{3'}+\r_{4'}-1)} \cr + i \int d123
\left[\right.  && T_{32,11'}\P_{323'4'}\d^{1}_{2'}\d^{32}_{11'} -
T_{4'1,23}\P_{1'2'23}\d^{1}_{3'}\d^{32}_{14'} \cr -T_{32,11'} \cdot
\left( \right.  &&\r_2'(+\P_{323'4'} \d^{2'}_{1}-\P_{3214'}
\d^{2'}_{3'} + \P_{3213'} \d^{2'}_{4'})+\cr &&\r_2(
\P_{32'14'}\d^{2}_{3'} - \P_{32'13'}\d^{2}_{4'})+\cr &&\r_3(
-\P_{22'14'}\d^{3}_{3'} + \P_{22'13'}\d^{3}_{4'})
+\underline{\underline{\Q_{322'13'4'}}} \left.\right) \d^{32}_{11'}\cr
+ T_{32,12'}\cdot \left( \right.  &&\r_2(-\P_{1'314'}\d^{2}_{3'} +
\P_{1'313'}\d^{2}_{4'})+\cr &&\r_3(+\P_{1'214'}\d^{3}_{3'} -
\P_{1'213'}\d^{3}_{4'})+\cr
&&\r_1'(+\P_{323'4'}\d^{1'}_{1}-\P_{3214'}\d^{1'}_{3'} +
\P_{3213'}\d^{1'}_{4'}) +\underline{\underline{\Q_{1'3213'4'}}}
\left.\right) \d^{32}_{12'}\cr - T_{3'1,23} \left( \right.  &&\r_1( +
\P_{1'2'23}\d^{1}_{4'})+\cr
&&\r_2'(-\P_{1'134'}\d^{2'}_{2}+\P_{1'124'}\d^{2'}_{3} -
\P_{1'123}\d^{2'}_{4'})+\cr
&&\r_1'(+\P_{2'134'}\d^{1'}_{2}-\P_{2'124'}\d^{1'}_{3} +
\P_{2'123}\d^{1'}_{4'}) + \underline{\underline{\Q_{1'2'1234'}}}
\left.\right) \d^{32}_{13'}+ \cr T_{4'1,23}\left(\right.
&&\r_1(+\P_{1'2'23}\d^{1}_{3})+ \cr
&&\r_2'(-\P_{1'123}\d^{2'}_{3'}+\P_{1'13'3}\d^{2'}_{2} -
\P_{1'13'2}\d^{2'}_{3})+\cr
&&\r_1'(+\P_{2'123}\d^{1'}_{3'}-\P_{2'13'3}\d^{1'}_{2} +
\P_{2'13'2}\d^{1'}_{3}) +\underline{\underline{\Q_{1'2'13'23}}} \left.
\left.  \right)\d^{32}_{14'} \right],\NN\label{PPP} \EEA on ${\bbox
k_{1'}}+{\bbox k_{2'}}={\bbox k_{3'}}+{\bbox k_{4'}}$, and the
Hartree-Fock self-energy is
$$\tilde \D^{1'2'}_{3'4'}=\o_{1'}+\o_{2'}-\o_{3'}-\o_{4'}+ 2 \int d1
\left((T_{1'1,11'}+ T_{2'1,12'}-T_{3'1,13'}-T_{4'1,14'})n_1\right).$$

At this stage, it is worthwhile pointing out precisely those terms,
underlined in  (\ref{PPP})  that 
give rise to the various long term effects:
\begin{enumerate}
\item The terms that gives rise to particle number transfer in the 
QKE are
$${2 i T_{4'3',2'1'}\r_{3'}\r_{4'}(1-\r_{2'}-\r_{1'})+ 2 i
  T_{4'3',2'1'}\r_{1'}\r_{2'}(\r_{3'}+\r_{4'}-1)}.$$ The reason is
  that when one solves for $\P_{1234}$, one obtains an expression
  which contains this term multiplied by
$$A_t(\D^{1'2'}_{3'4'})=\int\limits_0^t d \t \exp[i \D^{1'2'}_{3'4'}t]
=(i \D^{1'2'}_{3'4'})^{-1}(\exp{(\D^{1'2'}_{3'4'}t)}-1).$$
In the long time limit, 
$$\lim_{t\to\infty} A_t(x)=\pi \sgn t \cdot \d(x)+i P(\frac{1}{x}).$$
Under the operator $\Im$ in (\ref{rotimeevol}), the delta-function is
counted twice, and the principal value term cancels. The observant
reader will notice that the QKE can be effectively derived by simply
ignoring all terms in the equation for $\P_{1'2'3'4'}$ (except
$i\D^{1'2'}_{3'4'}\P_{1'2'3'4'}$) proportional to cumulants of order
greater than two. In the literature, this is called the Hartree-Fock
approximation.  What we show in this paper is that, for the magnitude
$\e$ of the coupling coefficient uniformly (in ${\bbox k}$) small, the
Hartree-Fock approximation is indeed self consistent when one takes
proper account of the frequency renormalization {\it to order $\e^2$.}
\item
The order $\e$ renormalizations to the frequency comes from the
decomposition of sixth order moments such as
$N_{1'2'13'23}\d_{11'}^{23}$ in (\ref{N6timeev}). These give rise to
terms in the equation for $\P_{1'2'3'4'}$ which are proportional to
$\P_{1'2'3'4'}$ itself. Indeed one obtains one such contribution from
each of the sixth order moments in (\ref{N6timeev}) leading to an
expressions in the equation for $\P_{1'2'3'4'}$ equal to
$$2 i \int d1 (T_{1'1,11'}+T_{2'1,12'}-T_{3'1,13'}-T_{4'1,14'})
\r_1.$$ When added to the frequency factor
$i(\o_{1'}+\o_{2'}-\o_3-\o_4)$, we obtain the term denoted by $\tilde
\D_{3'4'}^{1'2'}$ in (\ref{PPP}). It is not too difficult to see that,
in the equation for every cumulant $Q_N$, there is a term proportional
to $$i Q_N \sum\limits_{i=1}^N(\o_j+2\int T_{j1,1j}\r_1 d 1)$$ and
this gives rise to the first contribution in the renormalization of
the frequency (\ref{great}).  In the literature, $\tilde
\D_{3'4'}^{1'2'}$ is called Hartree-Fock self energy for fourth order
averages.  \item The order $\e^2$ terms in the renormalization to the
frequency arise from the terms \BEA i\int &&\left(
T_{32,11'}Q_{322'13'4'}\d^{32}_{11'}+
T_{32,12'}Q_{1'3213'4'}\d^{32}_{12'}\right.\NN &&\left.  +
T_{3',23}\Q_{1'2'1234'}+ T_{3'1,23}Q_{1'2'13'23}\d^{32}_{13'}\right) d
123 \NN\EEA containing the sixth order cumulants in (\ref{PPP}). The
equations for the sixth order cumulant contains, in addition to terms
proportional to a product of lower order particle number densities,
terms proportional to $\P_{1'2'3'4'}$ with a factor containing
$\r_k(1-\r_k)$.
\item All other terms are integrals which contain highly oscillatory
factors which, because of the Riemann-Lebesgue lemma, contribute
nothing in the long time limit.
\end{enumerate} 
\subsection{Asymptotic Expansions and Closure} We 
take advantage of the small parameter, the strength $\e$ of the
coupling coefficient and make the formal substitution $T_{12,34}\to\e
T_{12,34}$. We expand all cumulants $Q_{N}$ in an asymptotic expansion
\BE Q_2\equiv\r_k=n_k+\e Q_{2}^{(1)}+\e^2
Q_{2}^{(2)}+...\label{oddvaa} \EE and \BE Q_{N}(k_1,k_2,k_3...k_N;t)=
Q_N^{(0)}+\e Q_N^{(1)}+\e^2 Q_N^{(2)}+ ...\ \ \ \ \ \ \ \ N>2
\label{oddvab} \EE Because of resonant interactions, these asymptotic
expansions will be nonuniform in time.  Namely, terms proportional to
$\e t$, $\e^2 t$, etc., will appear in $Q_{N}$.  These terms are
removed and the asymptotic expansions {(\ref{oddvaa},\ref{oddvab}) }
rendered well ordered by
allowing corrections to the time dependence of the leading order
cumulants $n_k$, $Q^{(0)}_N, N>2$. We make the ansatz \BEA
\frac{\partial n_k}{\partial t} &=& \e^2 F_2^{(2)}+...\label{qw1}\\
\frac{\partial Q^{(0)}_{1,2,..N}}{\partial t} &=& i(\o_1+\o_2+... -
\o_{N})Q^{(0)}_{1,2,..N}+ \e F_{1,2,..N}^{(1)}+\e^2
F_{1,2,..N}^{(2)}+...  \CR\label{qw2} \EEA and choose $F_2^{(2)}, \ \
F_N^{(1)},\ \ F_N^{(2)}\ ... \ \ \ N>2$ in order to keep
{(\ref{oddvaa},\ref{oddvab}) }
asymptotically uniform for times $\o_{k_0}t=O(\e^{-2})$.  Equations
(\ref{qw1}) and (\ref{qw2}) describe the long time behavior of the
system.

This method goes by many names, averaging, multiple time scales etc.,
familiar to nonlinear physicists (see, e.g.,
\cite{b3},\cite{Newell}-\cite{Ben}).  It will turn out
that $F_2^{(2)}$ depends only on $n_k$ itself leading to a closed
equation for the particle number, the QKE.  It will also turn out that
$F_N^{(1)},\ \ F_N^{(2)},\ \ N>2$ are simply products of ${Q^{(0)}_N}$
with a function of $n_k$, which is the symmetric sum of $N$
components.  This means that the resulting equations (\ref{qw2}) are
easily solved by simply renormalizing the frequency.

The quantum kinetic equation is given by
\BEA\label{KE}
\frac{\partial} {\partial t}n_{k}\equiv\e^2F_2^{(2)}=
 (4\e^2  \int |T_{01',2'3'}|^2\d^{01'}_{2'3'}d1'2'3'\times
(\pi \sgn(t)\d(\tilde\Delta^{01'}_{2'3'})) \cr
\times(  n_{2'}n_{3'}(1-n_{1'}-n_{k})+
   n_{k}n_{1'}(n_{2'}+n_{3'}-1)))                      .
\EEA  
where the Hartree Fock self-energy is
$$\tilde \D^{1'2'}_{3'4'}=\o_{1'}+\o_{2'}-\o_{3'}-\o_{4'}+ 2 \int d1
\left((T_{1'1,11'}+ T_{2'1,12'}-T_{3'1,13'}-T_{4'1,14'})n_1\right).$$
From the form
of (\ref{KE}), it is clear that number density is redistributed by binary
particle collisions which satisfy momentum and energy conservation. In
particular, exchange of particle number (momenta, energy) is
associated with particles whose momenta and energies lie on the
resonant manifold defined to a good approximation by \BE {\bbox
k_1}+{\bbox k_2}-{\bbox k_3}-{\bbox k_4}=0,\ \ \ \ \ \
\o_{k_1}+\o_{k_2}-\o_{ k_3}-\o_{k_4}=0. \label{man}\EE

The evolution equation for $Q_{1'2'...N}^{(0)}$
can be written
\BEA 
\frac{\partial \ln{Q_{12..N}^{(0)}}}{\partial t}&=&
i \left(\sum_{i=1}^{N/2}\O_{i}-\sum_{i=N/2+1}^{N}\O^*_{i} \right),\NN
\O_{k'}&=&\o_{k'}+2\e\int d1 n_1 T_{1'1,11'}\CR&&+2 \e^2
\int d123 \left( n_1+n_2 n_3-n_1
 n_3- n_1 n_2 \right) 
T^2_{k'123}\d^{k'1}_{23}\CR&&
\ \ \ \ \ \ \ \ \ \ \times
\left(P(\frac{1}{\tilde\Delta^{k'1}_{23}})+
	i \pi \sgn(t)\d( \tilde\Delta^{k'1}_{23})\right) 
\CR\label{great}
\EEA 
which can be interpreted as a complex frequency renormalization.

We can calculate the sign of $\Im\O$ once $n_k$ reaches its 
steady (equilibrium) state. 
For steady state $\dot n_k=0$ we rewrite the QKE 
(\ref{KE}) as \BEA\label{KEstatic} 
\int 
 |T_{01',2'3'}|^2\d^{01'}_{2'3'}d1'2'3'\times (\pi 
\sgn(t)\d(\tilde\Delta^{01'}_{2'3'})) \cr \times(  
n_{2'}n_{3'}-n_{1'}n_{2'}-n_{1'}n_{3'}+n_{1'} )\CR
=\CR  \frac{1}{n_k}\int 
   |T_{01',2'3'}|^2\d^{01'}_{2'3'}d1'2'3'\times (\pi 
\sgn(t)\d(\tilde\Delta^{01'}_{2'3'})) \cr \times(  
n_{2'}n_{3'}(1-n_{1'}))\ge 0.\EEA
The LHS of the above equation is the imaginary part of $\O_k$. 
Observe that, because $\Im\O_k\ge 0$, the leading order 
approximation to the $N^{th}$ ($N\ge 2$) order cumulant decays  with 
time. This means that the memory of (smooth) initial states is gradually 
forgotten.
 
We want to make two very important points which are often overlooked.
While the leading order contributions (which at $t=0$ is the initial
state multiplied by an oscillatory factor) to the $N^{th}$ order
cumulants for $N>2$ play no role in the long time behavior of the
system and indeed slowly decay, higher order (in $\e$) contributions
do not disappear in the long time limit. The system retains a weakly
non Gaussian character which is responsible for and essential for
particle number and energy transfer. For example, in the long time
limit, the order $\e$ fourth order cumulant has the quasi stationary
contribution (the terms with higher order cumulants asymptote to zero
by means of phase mixing and the Riemann-Lebesgue lemma) \BE\nonumber
\P_{1'2'3'4'}^{(1)}(t)= 2i
T_{4'3',2'1'}(n_{3'}n_{4'}(1-n_{2'}-n_{1'})+
n_{1'}n_{2'}(n_{3'}+n_{4'}-1)) \times A_t(\tilde\Delta^{1'2'}_{3'4'})
\nonumber\EE where
$$A_t(x)=\int\limits_0^t d \t \exp[i x \t].$$ Because
$\lim_{t\to\infty} A_t(x)=\pi \sgn t \cdot \d(x)+i P(\frac{1}{x}),$ in
the limit $t\to\infty$, \BEA\nonumber \P_{1'2'3'4'}^{(1)}(t)\to 2i
T_{4'3',2'1'}(n_{3'}n_{4'}(1-n_{2'}-n_{1'})+
n_{1'}n_{2'}(n_{3'}+n_{4'}-1)) \nonumber\cr\times\left( \pi \sgn t
\cdot \d(\tilde\Delta^{1'2'}_{3'4'})+i P(\frac{1}
{\tilde\Delta^{1'2'}_{3'4'}}) \right) .\nonumber\EEA Thus, in the long
time limit, the order $\e$ contribution $\P_{1234}^{(1)}$ to the
fourth order cumulant is not smooth but is given by a sum of
generalized functions represented by the Dirac delta function and the
Cauchy Principal value. We may therefore legitimately ask: in what
sense is the asymptotic series {(\ref{oddvaa},\ref{oddvab}) } well
ordered if it contains terms
which are products of powers of $\e$ with generalized functions? The
answer is that to analyze properly the asymptotic behavior of the
system, we must always revert to physical space and look at the
corresponding asymptotic expansion for the cumulants
$\{\psi^\dagger({\bbox r})\psi^\dagger({\bbox r'})\psi({\bbox
r''})\psi({\bbox r'''})\}$ connected with fourth order expectation
values of the field operators. These objects, which decay to zero as
the separations ${\bbox r'}-{\bbox r},\ {\bbox r''}-{\bbox r},\ \
{\bbox r'''}-{\bbox r}$ tend to infinity, are simply the Fourier
transform of $\P_{kk'k''k'''}$. We can also show that no terms worse
than single delta functions occur (or products of delta functions which
have their support on different resonant submanifolds) at later powers
of $\e$ so that the resulting asymptotic expansion for the spatial
cumulants is indeed well ordered.

The second important point concerns the reversibility or rather the
retracebility of solutions of (\ref{KE},\ref{great}). In the
derivation of (\ref{KE},\ref{great}), we assumed that the initial
cumulants were sufficiently smooth so that integrals over momentum
space of multiplications of the initial values of $Q_N$ by
$\exp(-i(\o_1+.... - \o_N)t)$ tend to zero in the asymptotic
limit. However, it is clear from (\ref{great}) that the regenerated
cumulants have terms of higher order in $\e$ which are not smooth and
indeed have their (singular) support precisely on the resonant
manifold which is the exponent of the oscillatory exponential. What
would happen, then, if one were to redo the initial value problem from
a later time $t_1=O(\e^{-2})$, either positive or negative, after
which the fourth order cumulant had developed a nonsmooth part? On the
surface, it would seem that the $\sgn$ term in (\ref{great}) would be
$\sgn(t-t_1)$ so that, at every time $t_1$, there would be a
discontinuity in the slope of $n_k(t)$. But that is not the case. If
one accounts for the nonsmooth behavior (\ref{great}) in the new
initial value for $\P_{1234}^{(1)}$, then one gets additional terms in
(\ref{great}) which give exactly the same collision integral but with
the factor $ (\sgn t-\sgn(t-t_1))$.  Adding the two contributions, we
find the QKE is exactly the same as the one derived beginning at
$t=0$.  It is not that the point $t=0$ is so special. Rather, there is
a range of times $t$, $-\e^{-2}\ll\o_{k_0}t\ll\e^{-2}$ such that, if
one begins anywhere within this range, an initially smooth
distribution stays smooth. But once the limit $t\to\infty$, $\e^2 t $
finite, is taken, an irreversibility and nonsmoothness in the
cumulants is introduced.

In a very real sense, then, the infinite dimensional Hamiltonian
system acts as if there is an attracting manifold (an inertial or
generalized center manifold in the modern vernacular) in its phase
space to which the system relaxes as $\o_{k_0}t\to O(\e^{-2})$ (in
either time direction) on which the slow dynamics is given by the
closure equations (\ref{KE}),(\ref{great}).  On this attracting
manifold, the higher order cumulants are essentially slaved to the
particle number density and their frequencies are renormalized by
contributions which also depends on particle number density. The
attenuation in this case is due to losses to the heat bath consisting
of all momenta which do not lie on a resonance manifold associated
with ${\bbox k}$.

\section{Analysis of the Kinetic Equation.} \setcounter{equation}{0} 
\subsection{Conservation Laws, Thermodynamic and Finite-Flux 
solutions of the Kinetic Equation.}

The collision integral in (\ref{KE}) has the following constants of
motion \BE {\cal N}= \frac{V}{(2\pi)^d} \int d {\bbox k} n_k, \ \ \ {\cal
{\bbox P}}= \frac{V}{(2\pi)^d} \int \hbar{\bbox k} d {\bbox k} n_k, \ \ \ \  {\cal
E}= \frac{V}{(2\pi)^d} \int \frac{\hbar^2|\vec k|^2}{2m} d {\bbox k} n_k,
\label{ConsLaws}\EE which can be identified as number of particles,
momentum and kinetic energy. In a spatially homogeneous systems,
${\bbox P}=0$ so that the only relevant constants of motion are $\cal
N$ and $\cal E$.

In this article we will be dealing with the isotropic case only, and,
for simplicity, neglect the spin degree of freedom.  Therefore we
simplify the collision integral by averaging it over all
angles. First, we change variables from particles momentum $\bbox k$
to the particle kinetic energy \BE\e_k=\hbar \o_k=\frac{\hbar^2 k^2}{2
m},\label{dispersion}\EE where $m$ is the coefficient of
proportionality between $\o_k$ and $\hbar k^2/2$, and can be
associated with (effective) mass of the interacting particles.

We introduce  
$n_{\o}=n({\bf k}(\o))$ and rewrite the kinetic equation as 
\BEA \dot n_\o 
= \frac{1}{\O_0 k^{d-1} (dk/d\o)}\times&& 
\int\int\limits_{\o_i>0,i=1,2,3}\int{\cal K}(\o,\o_1,\o_2,\o_3) S_{\o \o_1 
\o_2 \o_3}\CR&& \ \ \ \ \ \ \ \ \ \ \ \times\d(\o+\o_1-\o_2-\o_3) 
d\o_1d\o_2d\o_3 \NN \label{KEfreq}\EEA where $S_{\o \o_1 \o_2 \o_3}$ is 
the angle-averaged potential,  \BE S_{\o \o_1 \o_2 \o_3}=4\pi\O_0(k k_1 
k_2 k_3)^{(d-1)} \frac{d k}{d\o} \frac{d k_1}{d\o_1}\frac{d 
k_2}{d\o_2}\frac{d k_3}{d\o_3}\cdot <|T_{kk_1,k_2k_3}|^2 \d({\bf k}+{\bf 
k_1}-{\bf k_2}-{\bf k_3})>.\label{potential}\EE
The 
 brackets $<..>$ denote averages  over unit spheres 
(including the $\frac{V}{(2\pi)^d}$ factors)
in 
${\bbox k},{\bbox k_1},{\bbox k_2},{\bbox k_3}$ space (i.e., we have 
integrated over all angular contributions) and $\O_0$ is the surface 
area of the unit sphere in $d$ dimensions.  Although in
our case $d=3$, we want to stress with the above notation 
the results can easily be extended to other dimensions,
for example to $d=2$ which is appropriate for the description of
semiconductor quantum wells. 
$\cal K$ is the kernel of 
kinetic equation, which for the quantum (fermionic), quantum (bosonic) and 
classical cases is  given respectively by:  \BEA{\cal K^{\rm 
fermionic}}(\o,\o_1,\o_2,\o_3)&=& n_{\o_2} n_{\o_3} (1-n_{\o_1})(1-n_{\o}) 
- n_{\o} n_{\o_1} (1-n_{\o_2})(1-n_{\o_3}),\NN {\cal K}^{\rm bosonic} 
(\o,\o_1,\o_2,\o_3)&=& (n_{\o_2} n_{\o_3}(n_{\o}+n_{\o_1}+1)-n_{\o} 
n_{\o_1} (n_{\o_2}+n_{\o_3}+1)), \NN {\cal K^{\rm 
classical}}({\o,\o_1,\o_2,\o_3})&=& (n_{\o_2} 
n_{\o_3}(n_{\o}+n_{\o_1})-n_{\o} n_{\o_1} 
  (n_{\o_2}+n_{\o_3})).\label{KernalS}\EEA   
We then introduce  the particle density 
per frequency  $$N_\o={\O_0 k^{d-1} (dk/d\o)}n_\o$$ so that $\int 
 N_\o d\o=\int n_k d{\bbox k}$ and \BEA \dot N_\o = 
\int\limits_{\o_1,\o_2,\o_3>0}\int\int {\cal 
K}(\o,\o_1,\o_2,\o_3)\times 
S_{\o \o_1 \o_2 \o_3}\d(\o+\o_1-\o_2-\o_3) d\o_1d\o_2d\o_3.\CR 
\label{KEfreq2}\EEA 
One class of steady (equilibrium) solutions of the KE corresponds to
the thermal equilibrium. For fermionic systems, it is given by the
Fermi Dirac (FD) distribution, and for bosonic systems it is given by
the Bose-Einstein distribution, \BE
n_{k}=\frac{1}{\exp{(\b(\e_k-\mu))}\pm 1}\label{FD}\EE where the plus
sign corresponds to the fermion case and the minus sign to the bosonic
case.  $\mu$ is the chemical potential and $\b$ is the inverse
temperature in energy units.  The classical analogue of the quantum
thermal equilibrium distribution is given by the Rayleigh-Jeans
distribution \BE n(k)=T/(\mu + \e_k) \label{RG}.\EE It is easy to
check that the solutions (\ref{FD},\ref{RG}) make the integrand in
(\ref{KEfreq},\ref{KEfreq2}) exactly zero in all three cases.

However, the thermodynamic equilibrium is not the most general {
steady (equilibrium)} solution of the kinetic equation and indeed in
some cases has little relevance. The solutions we are most interested
in are those which describe the steady state reached {\it between}
ranges of frequencies where particles and energy are added to or
removed from the system. These regions, where there is no pumping or
dumping, are called "windows of transparency" or "inertial
ranges". In particular, we have in mind the following
situation. Particles and energy are added to the system in a narrow
range of intermediate frequencies about $\o_0$. Particles and energy
are drained from the system in a range of frequencies about
$\o<\o_L<\o_0$ and for $\o>\o_R>\o_0$.  Because of conservation of
energy and particles in the inertial ranges between $\o_L$ and $\o_0$
and between $\o_0$ and $\o_R$ where there is no pumping or damping and
because the relations between particle number $N_\o$ and energy
density $E_\o=\o N_\o$, we will find that a net flux of energy to the
higher frequencies must be accompanied by a net flux of particles to
lower frequencies {as it might be expected by analogy with classical
wave turbulence}. The presence of sources and sinks drives the system
away from the thermodynamic equilibrium. Therefore, in the windows of
transparency, $\o_L<\o<\o_0$ and $\o_0<\o<\o_R$, the system can also
relax to equilibrium distributions corresponding to a finite flux of
particles and energy flowing through these windows from the sources to
the sinks. These are the new solutions of the QKE. The number of such
finite flux solutions corresponds to the number of conserved densities
(here two, $n_k$ and $\o_k n_k$, or $N_\o $ and $\o N_\o$) of the QKE.

To demonstrate the existence of such solutions, we rewrite the KE in
the following form: 
\BEA \dot N_\o& =&
\frac{\partial^2}{\partial \o^2} {\cal W}[n_\o],\NN {\cal
W}[n_\o]&=&\int\int\limits_\O\int 
(\o+\o_1-\o_2-\o_3) \times
{\cal K}_{\o_2+\o_3-\o,\o_2,\o_3}S_{\o_2+\o_3-\o,\o_2,\o_3}  d\o_1 d\o_2
d\o_3,\NN  
 \label{FluxTheorem}\EEA  
where the integration is over the region 
$\ \O \{\o_1,\o_2,\o_3>0\ \ \o_1<\o_2+\o_3<\o+\o_1\}$
{This expression can be checked
by direct differentiation.}

The relevant equation kinetic equation, which includes the presence of sources and
sinks is \BE \frac{\partial N_\o}{\partial t}=
\frac{\partial^2}{\partial\o^2}{\cal W}[n_\o]+F_\o-D_\o,
\label{N310}\EE where we think of $F_\o$ as having its support near
$\o=\o_0$ and $D_\o$ its support below $\o_L$ and above $\o_R$. We ask
if (\ref{N310}) leads to a steady (equilibrium) solutions in the
{transparency } regions $\o_L<\o<\o_0$ and $\o_0<\o<\o_R$ where
$F_\o=D_\o=0$ corresponding to a finite flux of particles and energy
across these windows. One can readily associate the quantities \BEA
&&Q=\frac{\partial {\cal W}}{\partial \o},\ \frac{\partial
N_\o}{\partial t}=\frac{\partial Q}{\partial \o}, \NN {\rm and }\ \ \
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ &&\NN &&P={\cal W} -\o\frac{\partial
{\cal W}}{\partial \o},\ \frac{\partial \o N_\o}{\partial
t}=-\frac{\partial P}{\partial \o},
\label{Fluxdef}
\EEA with the fluxes $P$ and $Q$ of particles $N_\o$ and energy
$E_\o=\o N_\o$.  $Q$ and $P$ are taken positive if leftward and
rightward flowing respectively. In the windows of transparency, we
look for solutions for which $N_\o$ is constant in time and then
(\ref{FluxTheorem}) integrates to \BE {\cal W}=Q\o+P\label{N313}\EE
where $Q$ and $P$ are constants.  The two parameter family of
thermodynamic solutions, parametrized by $T$ and $\mu$ is given by
solving the homogeneous equation ${\cal W}=0$ for which
$P=Q=0$. Therefore, the thermodynamic solutions carry no fluxes of
particles or energy.  A more general\footnote{ On purpose we do not
say "the most general solutions", as there may be some hidden
symmetries of the QKE which generate more conserved quantities and
thus more general solutions, corresponding to fluxes of those
quantities. The statement that the general KE with $M$ conserved
quantities has "a most general solution" depending on $2M$ parameters,
is also not proven yet for the general collision integral.  Remember
also that we are considering the isotropic case. This is, of course,
an idealization, and generally "drift" solutions may be of equal
importance.} steady (equilibrium) solution to (\ref{FluxTheorem})
therefore is the four-parameter family \BE n_\o=n_\o(T,\mu,P,Q). \EE

We are particularly interested in the solutions for which $Q_0$ particles 
per unit time and $\o_0 Q_0$ units of energy per 
unit time are fed to the system in a narrow frequency 
window about $\o=\o_0$. We will assume that the flux of particles passing 
through the left (right) window $\o_L<\o<\o_0$ ($\o_0<\o<\o_R$) is $Q_L$ 
($Q_R$) and the flux of energy though the  right (left) window is $P_R$ 
($P_L$). We will also assume that the sinks consume all the particles and 
energy that reach them. Then (see Figure 1)

\BEA Q_L-Q_R=Q_0\NN P_R-P_L=\o_0 Q_0 \NN P_L=-\o_L Q_L\NN 
Q_R=-\frac{1}{\o_R} P_R.\label{PQrela}\EEA

\

\


%picture1
\noindent
\noindent
\noindent
\setlength{\unitlength}{.8cm}
\begin{picture}(16,4)(0,1)
\put(0,1.5){\line(1,0){16}} \put(8,1.3){\line(0,1){.4}}
\put(8,1){\makebox(0,0){$\omega_0$}} \put(8,2.2){\vector(0,1){1.6}}
\put(8,0.5){\makebox(0,0){IN}} \put(7.6,4.3){$Q_0,\omega_0 Q_0$}
\put(1.25,1.3){\line(0,1){.4}} \put(1.25,1){\makebox(0,0){$\omega_L$}}
\put(14.6,1.3){\line(0,1){.4}} \put(14.6,1){\makebox(0,0){$\omega_R$}}
\put(7,1.2){\line(0,1){2.8}} \put(9,1.2){\line(0,1){2.8}}
\put(2.5,1.2){\line(0,1){2.8}} \put(1,2){\vector(0,-1){.8}}
\put(1,0.5){\makebox(0,0){OUT}} \put(15,2){\vector(0,-1){.8}}
\put(15.2,0.5){\makebox(0,0){OUT}} \put(10.2,2){\vector(-1,0){1}}
\put(10.5,1.9){$Q_R$} \put(9.2,2.8){\vector(1,0){1}}
\put(10.5,2.7){$P_R$} \put(6.8,2){\vector(-1,0){1}} \put(5,1.9){$Q_L$}
\put(5.8,2.8){\vector(1,0){1}} \put(5,2.7){$P_L$}
\noindent
\noindent
\noindent
\end{picture}

\vskip 1cm
{\begin{centerline}{\large \bf Figure 1}\end{centerline}} 

\vskip 3cm The first two relations in (\ref{PQrela}) express
conservation of particles and energy. The second two express the fact
that, in order to maintain equilibrium, the rate of particle
destruction at $\o_R$ is the rate of energy destroyed there divided by
the energy per particle.  Likewise the amount of energy destroyed at
$\o_L$ (which absorption, in the context of application discussed in
the section 5, will be due to semiconductor lasing) must be $\o_L$
times the number of particles absorbed there.  Solving (\ref{PQrela})
we obtain \BEA Q_L&=&Q_0 (\o_R-\o_0)/(\o_R-\o_L),\CR
Q_R&=&-Q_0(\o_0-\o_L)/(\o_R-\o_L),\CR P_R&=&Q_0 \o_R
(\o_0-\o_L)/(\o_R-\o_L)\CR P_L&=& -\o_L Q_0
(\o_R-\o_0)/(\o_R-\o_L).\label{PQSol}\EEA {We see that for}
$\o_L\ll\o_0\ll\o_R$, $|Q_R|\ll Q_L\simeq Q_0, |P_L|\ll P_R\simeq \o_0
Q_0$ so that the solutions are almost pure Kolmogorov in the sense
that almost all energy flows to $\o_R$ and almost all particles flow
to $\o_L$.  However, it is important to stress that in the left
window, (\ref{N313}) becomes \BE {\cal W}=Q_L(\o-\o_L)\label{N318}\EE
and in the right window \BE {\cal
W}=P_R(1-\frac{\o}{\o_R})\label{N319}\EE so that the right hand sides
vanish if the frequency approaches the sink value. This aids
convergence in these windows.

Solutions to (\ref{N318}), (\ref{N319}) have not been investigated
even in the classical case. In the classical case, Zakharov (see,
e.g., \cite{Z68a},\cite{Z68b}) had found the pure Kolmogorov solutions
$T=\mu=P=0$, $T=\mu=Q=0$ which turns out to have power law behavior
$n_{\o}=c\o^{-x}$. Likewise, in the bosonic case, several authors have
attempted to find power law solutions which essentially balance the
quadratic terms in ${\cal W}^{\rm bosonic}$ with a finite energy
flux. However, in the differential approximation, there are no power
law solutions.

In many cases it may be that whenever $\o_R\gg\o_0$, $\o_L$ may 
not be all that much smaller than $\o_0$. In particular, in order to 
exploit these solutions in the context of semiconductor lasers, it is 
advantageous to have $\o_0$ close enough to $\o_L$ to minimize energy 
losses (the ratio 
$P_R/|P_L|=\o_R(\o_0-\o_L)/( \o_L(\o_R-\o_0))$) but 
far enough away to facilitate  pumping unimpeded  by 
Pauli blocking.  We return to this application after introducing an 
enormous simplification for $\cal W$ which gives a very good qualitative 
description of the collision integral. 

\section{Differential Kinetic Equation} \setcounter{equation}{0} 
\subsection{Derivation of Differential Quantum Kinetic Equation}

To simplify our analysis, we use an approximation to the kinetic
equation known as the differential approximation (see, e.g.,
\cite{Optical},\cite {HASS},\cite{Balk}). The differential quantum kinetic
equation (DQKE) gives qualitatively correct behavior in general, but
is strictly valid only when the particle and energy transfer happens
primarily between neighbors in momentum space.  It is easy to verify
that it has a four-parameter-family of a steady (equilibrium)
solutions and it is easy to identify two of these parameters as the flux of
particle number and energy respectively and the other two as temperature
and chemical potential.  The analytical expressions for the fluxes can
be calculated so that, for any given distribution, the corresponding
fluxes may be easily computed numerically.  The steady (equilibrium)
solutions can be found analytically in various limits.  The DQKE is
much more suitable for numerical experiments than the full collision
integral, simply because it is easier and faster to compute
derivatives than integrals of the collision type.

{We now demonstrate the derivation of the differential approximation
in the fermion case. This result is new. The results for the cases of
classical waves and bosonic systems are given in (\ref{DKE}).}
Assume that $S_{\o \o_1 \o_2 \o_3}$ is dominated by its contribution
from the region $\o\simeq\o_1\simeq\o_2\simeq\o_3$.  We call such a
coupling coefficient ``strongly diagonal''. Then obviously the
integrand of the QKE deviates significantly from zero in the same
region.  The first derivation of the differential kinetic equation
proceeds as follows.  Multiply both sides of (\ref{KEfreq2})  by a
sufficiently smooth function $\F(\o)$ and integrate with respect to
$\o$: \BEA \int \dot N_\o \F(\o)d\o = &&\int\int\int {\cal
K}(\o,\o_1,\o_2,\o_3)\NN &&\times S_{\o,\o_1,\o_2,\o_3}
\d(\o+\o_1-\o_2-\o_3) \F(\o)d \o d\o_1d\o_2d\o_3\NN
\label{symm2}\EEA
Symmetrise the RHS of Eqn. (\ref{symm2}) to get  
\BEA                           
\int \dot N_\o \F(\o)d\o = 
\int\int\int
{\cal K}(\o,\o_1,\o_2,\o_3)\times S_{\o,\o_1,\o_2,\o_3}
\d(\o+\o_1-\o_2-\o_3) \NN\times\frac{1}{4}\left(\F(\o)+
\F(\o_1)-\F(\o_2)-\F(\o_3)\right)
d \o d\o_1d\o_2d\o_3
\label{symm}\EEA
To do this we use the symmetries of the kernel of collision 
integral responsible for particle number and 
energy  conservation.  Now make a change of variables 
$\o_i=\o+\D_i,\ \ i=1,2,3$ and expand $\F$ in the Taylor 
series with respect to $\D_i$ around $\o$. 
The first nonzero term in the expansion contains the second 
derivative $\F''(\o)$:
 \BEA 
\F(\o)+\F(\o+\D_1)-\F(\o+\D_2)&-&\F(\o+\D_1-\D_2)\NN && =  
(\D_1-\D_2)\D_2\F''(\o) +O(\D^2)\NN&&
\label{expansion}\EEA
We also expand $n_{\o_i}=n_{\o}+\D_i n_{\o}'+\D_i^2/2 n''_{\o}$ in the
kernel ${\cal K^{\rm fermionic}}(\o,\o_1,\o_2,\o_3)$ of the kinetic
equation, \BEA {\cal K^{\rm
fermionic}}(\o,\o+\D_1,\o+\D_2,\o+\D_1-\D_2)\NN=(\D_1-\D_2)\D_2
(n_\o'^2(1-2n_\o)+n_\o n''_\o(n_\o-1))\NN\label{expansion2}\EEA We
then substitute (\ref{expansion}) and (\ref{expansion2}) to
(\ref{symm}), integrate by parts the $\F''(\o)$ term to get \BEA \int
d \o\F(\o)[{\dot{ N_\o}}&-&\frac{\partial^2}{\partial\omega^2}
(n_{\o}'^2(1-2n_\o)+n_\o n{''}_\o(n_\o-1)) \NN&&\times\int d \D_1 d
\D_2 S_{\o,\o+\D_1,\o+\D_2,\o+\D_1-\D_2}(\D_2(\D_1-\D_2))^2]=0\NN
\label{rez}\EEA
Using the arbitrariness of $\F$ we finally get \BEA \dot n_\o = -
\frac{1}{\O_0 k^{d-1} (dk/d\o)}
\frac{1}{4}\frac{\partial^2}{\partial\o^2}\left[ \left( { n_\o^4
\frac{\partial^2}{\partial\o^2}(\frac{1}{n_\o})+n_\o^2
\frac{\partial^2}{\partial\o^2}(\ln (n_\o))}\right)  \right.
\nonumber \cr \times\left. \int d \D_1 d \D_2
S_{\o,\o+\D_1,\o+\D_2,\o+\D_1-\D_2}(\D_2(\D_1-\D_2))^2
\right]\nonumber.\EEA Finally, we assume that $S_{\o \o_1 \o_2 \o_3}$
is a homogeneous function of its arguments of degree $\g$: \BE
S_{\e\o,\e\o_1,\e\o_2,\e\o_3}=\e^\g S_{\o \o_1 \o_2 \o_3}.
\label{scaleini}\EE
 We define $\o \d_i=\D_i$, and rewrite the DQKE as \BEA \dot n_\o = - \frac{1}{\O_0 
k^{d-1} (dk/d\o)} \frac{\partial^2}{\partial\o^2} \left[\left( { n_\o^4 
\frac{\partial^2}{\partial\o^2}(\frac{1}{n_\o})+n_\o^2 
        \frac{\partial^2}{\partial\o^2}(\ln (n_\o))} \right)\times 
  I\times\o^s\right] \label{DKEfin1} \EEA where $I$ is the interaction 
strength 
$$ I = {\frac{1}{4}\int d \d_1 d \d_2
S_{1,1+\d_1,1+\d_2,1+\d_1-\d_2}(\d_2(\d_1-\d_2))^2 },$$ and $s=\g+6.$

An alternative derivation of the DQKE can be given by applying the
Zakharov transformation (see, e.g., \cite{Z68b}) directly to the QKE.
The Zakharov transformation is a conformal change of variables which
reveals the symmetry of original collision integral. It transforms
certain regions of the integration domain and, in the classical case,
makes the transformed collision integral have a zero integrand for
certain power law distributions. The Zakharov-transformed KE takes the
form \BEA \dot n_\o =\frac{1}{\O_0 k^{d-1} (dk/d\o)}\times
\frac{1}{4}\times
\int\limits_{\D_1}\int ({\cal
K}(\o,\o_1,\o_2,\o_3)+(\frac{\o}{\o_2})^{\g+3} {\cal
K}(\o,\frac{\o\o_3}{\o_2},\frac{\o^2}{\o_2},\frac{\o\o_1}{\o_2})\CR
+(\frac{\o}{\o_1})^{\g+3} {\cal K}(\o,\frac{\o^2}{\o_1},\frac{\o
\o_2}{\o_1},\frac{\o\o_3}{\o_1})+ (\frac{\o}{\o_3})^{\g+3}{\cal
K}(\o,\frac{\o\o_2}{\o_3}, \frac{\o\o_1}{\o_3},\frac{\o^2}{\o_3})) \CR
\times S_{\o \o_1 \o_2 \o_3}\d(\o+\o_1-\o_2-\o_3) d\o_1d\o_2d\o_3\CR
\label{zakh}\EEA
We then expand the RHS of the above equation in powers of $\D's$.  The
first nonvanishing term is of order $\D^4$ and can be represented as
second order derivative  with respect to $\o$ of 
$\left[\left( { n_\o^4
\frac{\partial^2}{\partial\o^2}(\frac{1}{n_\o})+n_\o^2
\frac{\partial^2}{\partial\o^2}(\ln (n_\o))} \right)\times
I\times\o^s\right]$.  The resulting DQKE is given in  (\ref{DKEfin1}).


\subsection{Solutions and properties of the DQKE}
Let us now rewrite the DQKE in the form: \BEA \dot N_\o&=&
\frac{\partial^2}{\partial\o^2} {\cal W}[n_\o],\CR {\cal W}^{\rm
fermionic} [n_\o]&=&-I\left( { n_\o^4
\frac{\partial^2}{\partial\o^2}(\frac{1}{n_\o})+n_\o^2
\frac{\partial^2}{\partial\o^2}(\ln (n_\o))}\right)\times\omega^s, \NN
{\cal W}^{\rm bosonic} [n_\o]&=&I\left( { n_\o^4
\frac{\partial^2}{\partial\o^2}(\frac{1}{n_\o})-n_\o^2
\frac{\partial^2}{\partial\o^2}(\ln (n_\o))}\right)\times\omega^s, \NN
{\cal W}^{\rm classical} [n_\o]&=&I \left( { n_\o^4
\frac{\partial^2}{\partial\o^2} (\frac{1}{n_\o})}
\right)\times\omega^s. \NN
\label{DKE}\EEA 
We can now use (\ref{Fluxdef}) to calculate the fluxes $P$ and $Q$ in
terms of $n_\o$ and its derivatives.  We concentrate on the fermionic
case. There, \BEA Q=\frac{\partial {\cal W}}{\partial\o}&=&I s
\o^{s-1}\left(-n'^2(2n-1)-n n''(1-n)\right)\CR &&\ \ \ + I
\o^s\left(-2n'^3+n'n''(1-2n)+n n'''(n-1) \right)\CR
P=\left({\cal{W}}-\o\frac{\partial{\cal W}}{\partial \o}\right)&=&\ I
\o^s(1-s)\left(-n'^2(2n-1)-n n''(1-n)\right) \CR && \ \ \ -I \o^{s+1}
\left(-2n'^3+n'n''(1-2n)+n n'''(n-1) \right).\CR\label{fluxes}\EEA Let
us make a change of variables $n=1/(G+1)$ and $n=1/(e^m+1)$.  $n, G, m$
are functions of $\o$ and $t$, ``dot'' is used to denote
differentiation with respect to time, and ``prime'' with respect to
$\o$.  Then \BEA \dot G &=& (1+G)^2 \frac{I}{\O_0 k^{d-1} (dk/d\o)}
\frac{\partial^2}{\partial \o^2}\left(\o^s\frac{G'^2-G G''}{(1+G)^4}
\right)\NN {\rm or} \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ &&\CR \dot m_\o
&=&- \frac{I}{\O_0 k^{d-1} (dk/d\o)} \cosh^2 (\frac{m}{2})
\frac{\partial^2}{\partial\o^2}\o^s \frac{m''}{4 \cosh^4(\frac{m}{2})}
\label{DKEsech}\EEA Since the stationary DQKE is a fourth order ODE,
its solutions will have four free parameters.  Indeed, assume a steady
(equilibrium) state and integrate (\ref{DKE},\ref{DKEsech}) twice to
get \BEA {\cal W}[n_\o]&=&Q \o + P,\CR \ \ \ {\rm or}\ \ \ \ \ \ \ &&
\CR G&=&\exp{(-\frac{\partial^{-2}}{\partial \o^2}
\frac{(Q\o+P)(G+1)^4}{I\o^s G^2})},\CR {\rm or} \ \ \ \ \ \ \ \ \ \ \
&& \CR m''&=& \frac{Q \o + P}{I\o^s}\cosh
^4(\frac{m}{2}),\label{DKEstat}\EEA where $Q$ and $P$ are fluxes of
particle number and energy.  For $P=Q=0$, (\ref{DKEstat}) trivially
gives
$$n(\o)=\frac{1}{e^{\frac{1}{T}( \o -\mu)}+1},$$
the Fermi-Dirac distribution function.
Therefore we observe that the Fermi-Dirac distribution 
function corresponds to a zero flux solution of the kinetic 
equation, {consistent with our findings of the previous section}.

\subsection{Numerical Results.}

We now investigate numerically the DQKE to gain 
an intuitive understanding of its properties. 

We begin by {studying} the time independent solutions of
(\ref{DKE}) and solve (\ref{N318}) as an initial value problem with
${\cal W}$ given by ${\cal W}^{\rm fermionic}$ (see (\ref{DKE})), \BE
Q_L(\o-\o_L)={\cal W}^{\rm fermionic}\equiv -I\left( { n_\o^4
\frac{\partial^2}{\partial\o^2}(\frac{1}{n_\o})+n_\o^2
\frac{\partial^2}{\partial\o^2} (\ln
(n_\o))}\right)\times\omega^s. \label{LeftWindow}\EE We start from
$\o=0.1$ and take the initial conditions to be
$n_{\o=0.1}=f_{\o=0.1},\ \ n'_{\o=0.1}=f'_{\o=0.1}$, where $f_{\o}$ is
a conveniently chosen Fermi-Dirac distribution function
$f_\o=1/(\exp{(7\o-4)}+1)$.  The results are shown on Figure 2a for
the range of $10^6 Q_L=0,3,6,9$.  For $Q_L=0$ we recover the
Fermi-Dirac distribution function $f_\o$. However, for nonzero $Q$ the
distribution deviates from the Fermi-Dirac function, and the bigger
the flux, the bigger the deviation. For the same initial conditions,
the graph of the number density $n_\o$ for the positive finite flux
distribution lies below that for the Fermi-Dirac
distribution. Further, observe that, for the Fermi-Dirac solution,
the chemical potential $\mu$ is simply the distance to the point of
the inflection of $n_\o$ and the temperature is the width of the
exponential decay region. For finite positive fluxes, the distance to
the point of inflection (which we call $\mu_{Q}$) is less than
$\mu_{Q=0}=\mu$. Likewise, $T_{Q}$, the width of the $n_\o$
distribution for a finite flux is less than $T_0=T$. As fluxes become
larger, the finite flux equilibrium distribution becomes more and more
like a Heavyside function with rapidly decreasing $\mu_{Q}$. This
means that $n_\o$ is very small for all $\o>\mu_Q$ and thus the system
can be pumped in this region without losses due to the Pauli
blocking. We exploit this in the next section where the application
to semiconductor lasers is considered.

In semiconductor lasers, it is the finite temperature effect of
broadening the Fermi-Dirac distribution that contributes to
inefficiency. If one could operate at $T=0$, then one could simply
choose the chemical potential, related to the total carrier number (see,
e.g., \cite{Koch}), so that the distribution cuts off immediately
after the lasing frequency. However, the finite temperature broadens
the distribution and means that one has to pump momentum values which
play no role in the lasing process. The effect of the finite flux is
to make $T$ effectively smaller.

We next solve (\ref{N319}) for the steady state solutions in the
larger momentum region with ${\cal W}^{\rm fermionic}$ : \BE {\cal
W}=P_R(1-\o/\o_R)={\cal W}^{\rm fermionic}\equiv -I\left( { n_\o^4
\frac{\partial^2}{\partial\o^2}(\frac{1}{n_\o})+n_\o^2
\frac{\partial^2}{\partial\o^2} (\ln (n_\o))}\right)\times\omega^s
\label{RightWindow}
\EE We start from $\o=0.3$ with initial conditions
$n_{\o=0.3}=f_{\o=0.3},\ \ n^\prime_{\o=0.3}=f^\prime_{\o=0.3}$.  The
results are shown on Figure 2b for the range of
$10^4P_R=0,3,6,9$. Again, we recover the Fermi-Dirac function for zero
flux, and observe a similar deviation from thermodynamical equilibrium
for nonzero fluxes, namely the effective chemical potential and
temperature diminishes with increasing of the flux value $P_R$.


We then consider the time evolution of the distribution function as
given by the DQKE. The fundamental property of the kinetic equation
that any distribution function relaxes to its thermodynamical
equilibrium value in the absence of forcing (pumping/damping) is also
true for the DQKE, as illustrated on Figure 3(a). There an initial
distribution function, shown by a thin solid line, relaxes to the FD
function, shown by a thick solid line, through several intermediate
states shown by dashed lines. Since there is no forcing to the system,
we take "fluxless" boundary conditions $P=Q=0$ in (\ref{fluxes}) on
the boundaries, so that no particles or energy cross the boundaries.
The distribution relaxes to the FD distribution roughly by the time
$\t_{\rm relax}$, which can be estimated as
$\o_{max}^{(2-s)}/(I(n+n^2))$, where $\o_{max}$ is the frequency where
distribution approaches zero value.  To check that the final
distribution is indeed FD, we calculate $\ln({1/n_\o-1})$ and verify
that it is a linear function.  The total particle number $N$ and
energy $E$ are conserved in our numerical runs to an accuracy of
$10^{-5}$.

We then address the question of what is the steady (equilibrium)
solution when the system has some external forcing. To model the
forcing, we specify some positive flux of $Q$ on the boundaries, and
wait for the distribution to reach a new equilibrium, a hybrid state
with a constant flux of $Q$, a zero $P$ flux, energy $E$ and particle
number $N$ (see Figure 3(b)). The more the flux of number of particles
is, the more the final distribution is bent in the manner of Figure 2
and according to (\ref{DKEstat}). The total particle number $N$ and
energy $E$ are the same for all curves on Figure 3(b).

\section{Application to Semiconductor Lasers}
We will now  investigate the relevance of finite flux solutions
in the context of semiconductor lasers.  
In many ways, semiconductor lasers are similar to two level lasers in
that the coherent light output is associated with the in-phase
transitions of an electron from a higher to lower energy state. In
semiconductors, the lower energy state is the valence band, from which
sea electrons are removed leaving behind positively charged holes. The
higher energy state is the conduction band.  The quantum of energy
released corresponds to an excited electron in the conduction band
combining with a hole directly below (because the emitted photon has
negligible momentum) in the lower band.

However, there are several important ways in which the semiconductor
laser differs from and is more complicated than the traditional
two-level laser model.  First, in order for there to be lasing,
there must be both an electron and hole available at the same
momentum (and spin) value whereas, in the traditional two level laser,
the ground state is always available for an excited electron. As a
result, the excitation levels of both electrons and holes must be
above a certain level. Second, there is a continuum of transition
energies parameterized by the electron momentum ${\bbox k}$ and the
laser output is a weighted sum of contributions from polarizations
corresponding to electron-hole pairs at each momentum value. In this
feature, the semiconductor laser resembles an inhomogeneosly broadened
two level laser. Third, electrons and holes interact with each other
via Coulomb forces. Although this interaction is screened by the
presence of many electrons and holes, it is nonetheless sufficiently
strong to lead to a nonlinear coupling between electrons and holes at
different momenta. It is this effect that gives rise to the coupling
coefficient $T_{kk_1,k_2k_3}$ of the earlier sections of this
paper. The net effect of these collisions is a redistribution of
carriers (the common name for both electrons and holes) across the
momentum spectrum. In fact it is the fastest ($\approx 100$ fs.)
process (for electric field pulses of duration greater than
picoseconds) and because of this, the gas of carriers essentially
relaxes to a distribution corresponding to an equilibrium of this
collision process.  This equilibrium state is commonly taken to be
that of thermodynamic equilibrium for fermion gases, the Fermi-Dirac
distribution characterized by two parameters, the chemical potential
$\mu$ and temperature $T$, slightly modified by the presence of
broadband pumping and damping. However, as we have shown, in
situations where there is applied forcing and damping and in
particular where these processes take place in separate regions of
momentum space, it is the finite flux equilibrium which is more
relevant.
 
In the semiconductor laser the applied forcing is usually the
electrical pumping process. The low-energy sink is the actual lasing
process, i.e. stimulated emission at the laser frequency.  The
second sink has contributions from a variety of processes. One
contribution is due to that fact that some of the charge carriers
with high kinetic energies can leave the optically active region
and, therefore, contribute to the electrial pumping current without
contributing to the light amplification. Other processes acting as
sinks are less well localized at high energies. They are distributed
over an extended range of momentum values. Examples are
non-radiative recombination of electron-hole pairs mediated by
impurities, dislocations, interface roughness, etc. In addition,
Auger processes contribute significantly to the damping (i.e., loss
of charge carriers).

Although it is beyond the scope of this paper, we should like to
mention that, for a complete description of relaxation and
thermalization processes in semiconductor lasers, one would also have
to take electron-phonon scattering into consideration.  This
interaction insures that the temperature of the electron-hole plasma
is driven towards the lattice temperature. The main electron-phonon
interactions involve longitudinal optical and acoustic phonons. The
former couple via the Fr\"{o}hlich interaction to the charge carriers
and typical thermalization times are almost as short as those due to
carrier-carrier interaction in semiconductor lasers (within a factor
of about 5).  Much slower (three to four orders of magnitude) is the
deformation potential coupling with acoustic phonons. Making the
assumption that the electron-electron interaction dominates over
electron-phonon interactions, we proceed now and investigate the role
of equilibrium distributions other than Fermi-Dirac distributions in
laser performance.  As we have shown in previous sections, there are
finite flux equilibria, for which there is a finite and constant flux
of carriers and energy across a given spectral window. It is the aim
of this work to suggest that these finite flux equilibria are more
relevant to situations in which energy and carriers are added in one
region of the spectrum, redistributed via collision processes to
another region where they are absorbed.  Moreover, it may be
advantageous to pump the laser in this way because such a strategy may
partially overcome the deleterious effects of Pauli blocking.  
In conventional diode lasers, pumping is a process in which charge
carriers are injected into the depletion layer region in a p-n or
p-i-n structure (see, for example, \cite{bhattacharya.94}). If the
active layer is bulk-like, this process is based on a regular
drift-diffusion current, whereas in quantum-well lasers there is the
additional process of carrier capture into the quantum well by means
of inelastic scattering processes.  In the following we restrict
ourselves to the simplest model for injection pumping that neglects
the intrinsically anisotropic aspect of injection pumping but includes
the basic features of Pauli blocking effects in the pump
process
\cite{Koch}.  Within this model, the rate of change of the carrier
distribution is proportional to $\L_k (1 - n_k)$, where the pump
coefficient $\L_k$ is taken as a Fermi function with a given density
modeling the incoming equilibrated carriers, $n_k$ is the actual
carrier distribution in the active region, and $(1 - n_k)$ takes into
account the Pauli blocking effects. This means that only non-occupied
states can be filled by the pump current.  Since $\L_k$ is an
function that extends over a large range of k-values, we call this
pump model ``broad band pumping.''  In contrast to broad band pumping
one can also pump a laser locally in momentum space. However, usually
such a local pump process is narrow-band optical pumping, but this is
not commonly used to pump semiconductor lasers. There are, however, also
electrical pumping schemes available, which are based on tunneling
processes, and which, in principle, can allow for selective and
localized pumping and damping (see, e.g., \cite{zhang-etal.96}).

In the following we will examine the laser process and discuss, in
particular, the influence of the pumping process and its relation to
the equilibrium distribution function in stationary laser
operation. We base our numerical solutions on a greatly simplified
laser model. We assume that the distribution functions for electrons
and holes are identical (in other words, we assume electrons and holes
to have identical effective masses); we model the cavity losses by a
simple phenomenological loss term in the propagation equation for the
light field amplitude; we assume ideal single-mode operation; we make
the rotating wave approximation in the equation for the distribution
functions and the optical polarization function $p_k$; we neglect all
electron-hole Coulomb correlations (the so-called Coulomb enhancement,
see, e.g., \cite{Koch} ); we neglect bandgap energy renormalization;
and we neglect, as mentioned above, electron-phonon interaction.  In
spite of the approximations made, our model still captures the basic
processes in a semiconductor laser. The equations of motion (a form of
the semiconductor Maxwell-Bloch equations, see, e.g.,
\cite{Koch},\cite{Koch1}) read: \BEA \frac{\partial e}{\partial t}&=&
i \frac{\O}{2\e_0} \frac{V}{(2\pi)^d}
\int\mu_k p_k d{\bbox k} - \g_E e ,\label{two}\\
\frac{\partial p_k}{\partial t}&=&(i\O-i\tilde\o_k-\g_P)p_k - \frac{i
\mu_k}{2\hbar}(2n_k-1)e,\label{three}\\ \frac{\partial n_k}{\partial
t}&=&\L_k(1-n_k) -\g_k n_k +\left( \frac{\partial n_k}{\partial t}
\right)_{\rm collision} -\frac{i}{2\hbar} \left( \mu_k p_k e^*-\m_k
p_k^* e\right). \label{four}\EEA Here $e(t)$ and $p_k(t)$ are the
electric field and polarization (at momentum $\bbox k$) envelopes of
the carrier wave $\exp{(-i\O t + i K z)}$ where $\O$ is the cavity
frequency (we assume single mode operation only) and $n_k$ is the
distribution function for electrons and holes. The constants $\g_E=6\times
10^{10} /sec$, $\g_P=10^{13}/sec$ model electric field losses and
polarization decay (dephasing), $\e_0$ is the dielectric constant,
$\mu_k$ is the weighting accorded to different $\bbox k$ momentum
values and is modeled by $\mu_k=\mu_{k=0}/(1+\e_k/\e_{\rm gap}),\ \
\mu_{k=o}=3/10^{10} M e$, $e$ is the electron charge,
$\g_k=10^{10}/sec$ represent nonradiative carrier damping. In
(\ref{four}), ${\L_k}$ is the pumping due to the injection current
(taken to be between $0.001 ps^{-1}$ and $0.1 ps^{-1}$) 
and in
(\ref{three}) $\hbar\tilde\o_k=\e_{\rm gap}+\e_{e,k}+\e_{h,k}$.  We
further assume that all fields are isotropic and make a convenient
transformation from $k$ ($=|{\bbox k}|$) to $\o$ via the dispersion
relation $\o=\o({\bbox k})$ defining the carrier density
$n_\o=n({\bbox k}(\o))$ and approximate the collision term $\left(
\frac{\partial n_k}{\partial t}\right)_{\rm collision}= \left(
\frac{\partial n_{k(\o)}}{\partial t}\right)_{\rm collision}$ by the
differential kinetic expression (\ref{DKEfin1}): \BEA {\left(
\frac{\partial n_{k(\o)}}{\partial t}\right)_{\rm collision} } &=& -
\frac{1}{\O_0 k^{d-1} (dk/d\o)} \NN&&\times\frac{\partial^2}{\partial\o^2}
\left[\left( { n_\o^4
\frac{\partial^2}{\partial\o^2}(\frac{1}{n_\o})+n_\o^2
\frac{\partial^2}{\partial\o^2}(\ln (n_\o))} \right)\times
I\times\o^s\right]. \NN\label{six} \EEA We choose the value of the
constant $I$ to ensure that a solution of (\ref{six}) relaxes in a time
of 100 fs to its equilibrium value and $s$ is taken to be $7$.

We now compare the laser efficiencies in two numerical experiments in
which we arrange to: (i) Pump broadly a across wide range of momenta, so
that the effective carrier distribution equilibrium has zero (or
small\footnote{ Even for broad band pumping the small amount (much
smaller than in local band pumping) of $P$ and $Q$ fluxes are also
excited, because the form of the pumping $f_k$ gets effectively
multiplied by the Pauli blocking $(1-n_k)$ term. But these fluxes are
much smaller and relatively local in $k$ space.})  flux. We take the
pump profile to be given by the Fermi-Dirac distribution. (ii) Pump
carriers and energy into a narrow band of frequencies about $\o_0$ and
simulate this by specifying carrier and energy flux rates $Q_L$ and
$P_L=-\o_L Q_L$ at the boundary $\o=\o_0$.  $P_L$ is chosen so that the
energy absorbed by the laser is consistent with the number of carriers
absorbed there.  $\o=\o_L<\o_0$ is the frequency at which the system
lases.  We compare only cases in which the total amount of energy
supplied is the same. Because of the distribution of the supply, the
particle number in the broad band pumping has to be higher.

In the first numerical experiment, we show that for a very small
amount of pumping the laser operates for the narrow band pumping, while it
fails to operate in the broad band pumping case (i.e.  the threshold
pumping value for narrow band pumping is much lower than in the broad
band case).  The carriers supplied through the pumping process get
totally absorbed by the global damping $-\g_k n_k$ in the
broad band pumping case. There  no lasing occurs.  In the local band
pumping case, 
for the same amount of energy and lower amount of particle supplies,
the  laser operates. This is the qualitative
difference between two cases.


The results of this numerical experiment are presented in Figure 4.
The narrow-band pumped laser switches on and generates a nonzero
output power. We pump in the narrow region around $\o_0\simeq 200 {\rm
meV}$ and we model this by specifying the boundary conditions at
$\o_0$ to correspond to carrier and energy flux rates $Q_L$ and
$P_L=-\o_L Q_L$ respectively.  The initial value of the distribution
function (shown by a thin line), taken to be just below the lasing
threshold, builds up because of the influx of particles and energy
from the right boundary (dashed lines) until the laser switches
on. The final (steady) distribution function is shown by a thick solid
line and corresponds to a flux of particles and energy from the right
boundary (where we add particles and energy) to the left boundary,
where the system lases (Figure 4a).  The output power as a function of
time is also shown (Figure 4b). Time is measured in units of
relaxation times $\t_{\rm relax}=100fs$. In the contrast, the
broad-band pumped laser fails to switch on for such weak pumping
because of pumping inefficiency due to Pauli blocking.

We then increase the level of pumping, to a point where the laser
turns on for both the broad band and narrow band pumping cases, and
examine the output power in both cases. It turns out that for the same
amount of energy pumped{\footnote {Much more carriers are pumped in
the case of broad band pumping because of the pump distribution.
Because there are fewer carriers in the narrow band pumping case,the
"effective" chemical potential also decreases, which increases
efficiency.  }}, and for almost the same amount of carriers pumped the
output power in the narrow band case is significantly higher than in
the broad band case. The results are presented in Figure 5. We first
pump broadly, so that the effective carrier distribution has (almost)
zero flux.  The initial distribution function (thin line, Figure 5a)
builds up because of a global pumping (dashed lines) until the laser
switches on. The final (steady) distribution function is shown by the
thick solid line.  The output power as a function of time is also
shown (Figure 5b).

If we pump in the narrow region around $\o_0\simeq 200 {\rm
meV}$ and we model this by specifying the carrier and energy flux rates
$Q_L$ and $P_L=-\o_L Q_L$, then the initial distribution function
(thin line) builds up because of an influx of particles and energy from
the right boundary (dashed lines) until the laser switches on. The final
(steady) distribution function is shown by a thick solid line and
corresponds to a flux of particles and energy from the right boundary
(where we add particles and energy) to the left boundary, where the
system lases (Figure 5c).  The output power as a function of time is
also shown (Figure 5d).

We observe that the output power is an order of magnitude bigger in the
case of narrow-band pumping for the same amount of energy
influx, at least in our  model.  This can be explained
qualitatively by noting that, in the case of broad-band pumping, most
of the particles are injected at momenta where the distribution function
is roughly between 1/2 and 1 (the condition necessary for lasing) and
thus Pauli blocking is significant. In contrast, when one pumps
at  high momentum values, where there are
almost no particles, Pauli blocking is negligible, so that for the same
amount of pumping, more carriers are able to reach active zone of the
laser and contribute to inversion.

These results certainly suggest that the possibility of using narrow
band pumping and the resulting finite flux equilibrium of the QKE is
an option which is worth exploring further.

\section{Conclusions}
\setcounter{equation}{0} 

We have shown that, for weak coupling, the full BBGKY hierarchy of
operator equations has a natural asymptotic closure consistent with
Hartree-Fock factorization and complex frequency renormalization.

We have demonstrated that the QKE has a new class of steady
(equilibrium) solutions which carry constant fluxes of particles and
energy across momentum space. We have explored these solutions by
using the differential approximation to the collision integral.

We have then applied these ideas to semiconductor lasers and shown
that the output power is much higher when the finite flux equilibrium
is excited.  While we do not claim that, if we take all effects into
account, the advantages of narrow-band pumping will necessary
 remain, we do suggest that the evidence suggests that it is a
possibility worth pursuing further.

\section{Acknowledgment}
The authors wish to thank the Arizona Center for Mathematical
Sciences (ACMS) for support. ACMS is sponsored by AFOSR contract
F49620-97-1-0002 under the University Research Initiative Program
at the University of Arizona, Department of Mathematics.

\section{Figure Captions} 
\begin{itemize} 
\item {\bf Figure 1.\\} This pictures explains the setup for
input-output fluxes.  Carriers and energy are added at $\o_0$ at rates
$Q_0$ and $\o_0 Q_0$ respectively.  Energy and some carriers are
dissipated at $\o_R>\o_0$ (an idealization) and carriers and some
energy are absorbed by the laser at $\o_L$. Finite flux stationary
solutions are realized in the windows $(\o_L,\o_0)$ and $(\o_0,\o_R)$
although in practice there will be some losses through both of these
regions.
\item {\bf Figure 2\\} {\bf (a)} Numerical solution of
(\ref{LeftWindow}). Initial conditions are $n_{\o=0.1}=f_{\o=0.1}$,
$n'_{\o=0.1}=f'_{\o=0.1}$, where $f_\o=1/(\exp{(7\o-4)}+1)$. The
results are shown for $10^6 Q_L=0,3,6,9$.  The dot size indicates the
value of the flux, with larger dots corresponding to larger flux.
\\{\bf (b)} Same as {\bf(a)} but with (\ref{RightWindow}) and
$n_{\o=0.3}=f_{\o=0.3},\ \ n^\prime_{\o=0.3}=f^\prime_{\o=0.3}$.  The
results are shown for $10^4P_R=0,3,6,9$.

\item {\bf Figure 3.\\} {\bf (a)} Time evolution of the distribution
function as described by (\ref{DKE}), with boundary conditions $P=Q=0$
at both ends.  The initial distribution (thin line) relaxes to a
Fermi-Dirac state (thick line).  Several intermediate states are shown
by long-dashed and short-dashed lines.\\ {\bf (b)} Same as {\bf
(a)}, but for $P=0, Q=Q_1>Q_0$ at both ends.  The initial distribution
function (thick line) relaxes to finite-$Q$-equilibria as shown by
the long-dashed line (only the final steady state is shown). 
Also shown are the results for 
the boundary conditions $P=0, Q=Q_1>Q_0$ at both ends (short-dashed
line) and for  $P=0, Q=Q_2>Q_1$ at both ends (dotted line).
\item {\bf Figure 4 }\\ 
Time dependence of distributions {\bf (a)} and laser intensity {\bf
(b)} according to (\ref{two}-\ref{six}). Local-band pumping is modeled
by specifying the carrier and energy flux rates $Q_L$ and $P_L=-\o_L
Q_L$ at the right ($\o_0\simeq 200 {\rm meV}$) boundary.  The initial
distribution function (shown by a thin line in this and consequent
figures) evolves through a number of intermediate states shown by
dashed lines, until the laser switches on. The final (steady)
distribution function is shown by a thick solid line (in this and
subsequent figures). Time is given in units of the
relaxation time $\t_{\rm relax}=100fs$ on this and subsequent figures.

The laser fails to switch on in the broad-band pumping case for these
levels of pumping.

\item {\bf Figure 5.}\\ 
Same as in Figure 4 but for increased pump rate. \\{\bf(a)} Broad-band
pumping {\bf(b)} Narrow-band pumping with $\o_0\simeq 200{\rm meV}$.
\end{itemize}
\section{Appendix. Diagrams}
\setcounter{equation}{0} In this section we present simple diagrams
which can help one to visualize the definition of cumulants.  These
diagrams are only meant as an aid to visualize the correlation
contribution and should not be confused with other diagrams such as
Feymann diagrams.  The following pictures are intended to illustrate
 the ways to
factorize $N_{1234}$, $N_{123456}$ and $N_{12345678}$.  Let us denote
by square boxes operator averages according to (\ref{average}), so that two
connected arrows present $\r$ (2 operator expectation value), and

\setlength{\unitlength}{.22cm}
\begin{picture}(54,6)(0,0)
\put(1,3){\makebox(0,0){$N_{1234}=$}}
\put(6,2){\vector(1,0){6}}
\put(6,4){\vector(1,0){6}}
\put(12,0){\line(1,0){6}}
\put(12,6){\line(1,0){6}}
\put(12,0){\line(0,1){6}}
\put(18,0){\line(0,1){6}}
\put(18,2){\vector(1,0){6}}
\put(18,4){\vector(1,0){6}}
\put(7,1){\makebox(0,0){$3$}}
\put(7,5){\makebox(0,0){$4$}}
\put(23,5){\makebox(0,0){$1$}}
\put(23,1){\makebox(0,0){$2$}}
\end{picture}
\\Let us present cumulants in the form of vertices, with incoming arrows
representing the arguments of annihilation operators, and the outgoing
the arguments of creation operators. Then the second order cumulant
(which is the same as the two operator average) is represented by two
arrows:

\begin{picture}(54,6)(0,0)
\put(31,3){\makebox(0,0){$\r_1\d^1_2=$}}
%\put(36,0){\vector(1,1){3}}
\put(36,6){\vector(1,-1){3}}
\put(39.25,3){\circle*{.5}}
\put(39.25,3){\vector(1,1){3}}
%\put(39.25,3){\vector(1,-1){3}}
\put(36,1.5){\makebox(0,0){$ $}}
\put(36,4.5){\makebox(0,0){$2$}}
\put(42,4.4){\makebox(0,0){$1$}}
\put(42,1.5){\makebox(0,0){$ $}}
\end{picture}
\\ The fourth order cumulant is  represented by four arrows:\\
\begin{picture}(54,6)(0,0)
\put(31,3){\makebox(0,0){$\P_{1234}=$}}
\put(36,0){\vector(1,1){3}}
\put(36,6){\vector(1,-1){3}}
\put(39.25,3){\circle*{.5}}
\put(39.25,3){\vector(1,1){3}}
\put(39.25,3){\vector(1,-1){3}}
\put(36,1.5){\makebox(0,0){$3$}}
\put(36,4.5){\makebox(0,0){$4$}}
\put(42,4.4){\makebox(0,0){$1$}}
\put(42,1.5){\makebox(0,0){$2$}}
\end{picture}
%Picture for $N_{<8>}$ factorization should contain one extra term:$
%n_.n_.\P_{<4>}$
\\so that the definition of fourth order cumulant $\P_{1234}$ is 
\BE\label{facN4d} N_{1234}= \r_1 \r_2 (\delta^{2}_{3}\delta^{1}_{4}-
\delta^{2}_{4}\delta^{1}_{3})+\P_{1234}\d^{12}_{34}.\label{dN4def2}\EE
This partition can be represented graphically as

\begin{picture}(54,6)(10,8)
\put(0,2){\vector(1,0){6}}
\put(0,4){\vector(1,0){6}}
\put(6,0){\line(1,0){6}}
\put(6,6){\line(1,0){6}}
\put(6,0){\line(0,1){6}}
\put(12,0){\line(0,1){6}}
\put(12,2){\vector(1,0){6}}
\put(12,4){\vector(1,0){6}}
\put(1,1){\makebox(0,0){$1$}}
\put(1,5){\makebox(0,0){$2$}}
\put(17,5){\makebox(0,0){$3$}}
\put(17,1){\makebox(0,0){$4$}}

\put(20,3){\makebox(0,0){$= -$}}

\put(23,-1){\vector(1,1){3}}
\put(23,7){\vector(1,-1){3}}
\put(26,4){\vector(1,1){3}}
\put(26,2){\vector(1,-1){3}}
\put(23,0.5){\makebox(0,0){$2$}}
\put(23,5.5){\makebox(0,0){$1$}}
\put(29,5.5){\makebox(0,0){$3$}}
\put(29,0.5){\makebox(0,0){$4$}}

\put(32,3){\makebox(0,0){$+$}}

\put(35,-1){\vector(1,1){3}}
\put(35,7){\vector(1,-1){3}}
\put(38,4){\vector(1,1){3}}
\put(38,2){\vector(1,-1){3}}
\put(35,0.5){\makebox(0,0){$2$}}
\put(35,5.5){\makebox(0,0){$1$}}
\put(41,5.5){\makebox(0,0){$4$}}
\put(41,0.5){\makebox(0,0){$3$}}

\put(43,3){\makebox(0,0){$+$}}

\put(46,0){\vector(1,1){3}}
\put(46,6){\vector(1,-1){3}}
\put(49.25,3){\circle*{.5}}
\put(49.25,3){\vector(1,1){3}}
\put(49.25,3){\vector(1,-1){3}}
\put(46,1.5){\makebox(0,0){$2$}}
\put(46,4.5){\makebox(0,0){$1$}}
\put(52,4.4){\makebox(0,0){$3$}}
\put(52,1.5){\makebox(0,0){$4$}}
\end{picture}
\vskip 2cm 
\noindent
 Because of the commutation relations (\ref{commutation}), if we
interchange two indices corresponding to two creation or two
annihilation operators, the average should change its sign, for
example: $N_{1234}=-N_{2134}$. The definitions of cumulants should not
contradict this property, so each product of lower order cumulants
should be either positive or negative, depending upon whether it
corresponds to an odd or even permutation. This explains the negative
sign in from of $\r_1 \r_2\delta^{2}_{4}\delta^{1}_{3}$ term in
(\ref{dN4def2}). Similarly, the definition of the sixth order cumulant
$\Q_{123456}$ is \BEA N_{123456}= \r_1 \r_2 \r_3 \cdot \left( \right.
\d^{3}_{4} (\d^{2}_{5}\d^{1}_{6}-\d^{1}_{5}\d^{2}_{6})+
\d^{3}_{5}(\d^{2}_{6}\d^{1}_{4}- \d^{2}_{4}\d^{1}_{6})+ \d^{3}_{6}
(\d^{1}_{5}\d^{2}_{4}- \d^{1}_{4}\d^{2}_{5}) \left.  \right)+\CR
\r_3[+\P_{1256}\d^{3}_{4}\d^{12}_{56}-\P_{1246}\d^{3}_{5}\d^{12}_{46}
+ \P_{1245}\d^{3}_{6}\d^{12}_{45}]+\cr
\r_2[-\P_{1356}\d^{2}_{4}\d^{13}_{56}+\P_{1346}\d^{2}_{5}\d^{13}_{46}
- \P_{1345}\d^{2}_{6}\d^{13}_{45}]+\cr
\r_1[+\P_{2356}\d^{1}_{4}\d^{23}_{56}-\P_{2346}\d^{1}_{5}\d^{23}_{46}
+ \P_{2345}\d^{1}_{6}\d^{23}_{45}]\cr+ \Q_{123456}\d^{123}_{456}.
\label{dN6defd} \EEA
This partition can be represented as

\begin{picture}(54,8)(13,16)

\put(0,2){\vector(1,0){6}}
\put(0,4){\vector(1,0){6}}
\put(0,6){\vector(1,0){6}}
\put(6,0){\line(1,0){6}}
\put(6,8){\line(1,0){6}}
\put(6,0){\line(0,1){8}}
\put(12,0){\line(0,1){8}}
\put(12,2){\vector(1,0){6}}
\put(12,4){\vector(1,0){6}}
\put(12,6){\vector(1,0){6}}
\put(1,3){\makebox(0,0){$1$}}
\put(1,5){\makebox(0,0){$2$}}
\put(1,7){\makebox(0,0){$3$}}
\put(17,7){\makebox(0,0){$4$}}
\put(17,5){\makebox(0,0){$5$}}
\put(17,3){\makebox(0,0){$6$}}


\put(19,4){\makebox(0,0){\tiny $\ =-$}}

%1
\put(21,0){\vector(1,1){2}}
\put(21,8){\vector(1,-1){2}}
\put(21,4){\vector(1,0){2}}
\put(23,4){\vector(1,0){2}}
\put(23,6){\vector(1,1){2}}
\put(23,2){\vector(1,-1){2}}
\put(21,1){\makebox(0,0){$3$}}
\put(21,3){\makebox(0,0){$2$}}
\put(21,7){\makebox(0,0){$1$}}
\put(25,7){\makebox(0,0){$4$}}
\put(25,3){\makebox(0,0){$5$}}
\put(25,1){\makebox(0,0){$6$}}
\put(27,4){\makebox(0,0){$+$}}

%2
\put(29,0){\vector(1,1){2}}
\put(29,8){\vector(1,-1){2}}
\put(29,4){\vector(1,0){2}}
\put(31,4){\vector(1,0){2}}
\put(31,6){\vector(1,1){2}}
\put(31,2){\vector(1,-1){2}}
\put(29,1){\makebox(0,0){$3$}}
\put(29,3){\makebox(0,0){$2$}}
\put(29,7){\makebox(0,0){$1$}}
\put(33,7){\makebox(0,0){$4$}}
\put(33,3){\makebox(0,0){$6$}}
\put(33,1){\makebox(0,0){$5$}}
\put(35,4){\makebox(0,0){$+$}}

%3
\put(37,0){\vector(1,1){2}}
\put(37,8){\vector(1,-1){2}}
\put(37,4){\vector(1,0){2}}
\put(39,4){\vector(1,0){2}}
\put(39,6){\vector(1,1){2}}
\put(39,2){\vector(1,-1){2}}
\put(37,1){\makebox(0,0){$3$}}
\put(37,3){\makebox(0,0){$2$}}
\put(37,7){\makebox(0,0){$1$}}
\put(41,7){\makebox(0,0){$5$}}
\put(41,3){\makebox(0,0){$4$}}
\put(41,1){\makebox(0,0){$6$}}
\put(43,4){\makebox(0,0){$-$}}

%4
\put(45,0){\vector(1,1){2}}
\put(45,8){\vector(1,-1){2}}
\put(45,4){\vector(1,0){2}}
\put(47,4){\vector(1,0){2}}
\put(47,6){\vector(1,1){2}}
\put(47,2){\vector(1,-1){2}}
\put(45,1){\makebox(0,0){$3$}}
\put(45,3){\makebox(0,0){$2$}}
\put(45,7){\makebox(0,0){$1$}}
\put(49,7){\makebox(0,0){$5$}}
\put(49,3){\makebox(0,0){$6$}}
\put(49,1){\makebox(0,0){$4$}}
\put(51,4){\makebox(0,0){$-$}}

%5
\put(53,0){\vector(1,1){2}}
\put(53,8){\vector(1,-1){2}}
\put(53,4){\vector(1,0){2}}
\put(55,4){\vector(1,0){2}}
\put(55,6){\vector(1,1){2}}
\put(55,2){\vector(1,-1){2}}
\put(53,1){\makebox(0,0){$3$}}
\put(53,3){\makebox(0,0){$2$}}
\put(53,7){\makebox(0,0){$1$}}
\put(57,7){\makebox(0,0){$6$}}
\put(57,3){\makebox(0,0){$4$}}
\put(57,1){\makebox(0,0){$5$}}
\put(59,4){\makebox(0,0){$+$}}

%6
\put(61,0){\vector(1,1){2}}
\put(61,8){\vector(1,-1){2}}
\put(61,4){\vector(1,0){2}}
\put(63,4){\vector(1,0){2}}
\put(63,6){\vector(1,1){2}}
\put(63,2){\vector(1,-1){2}}
\put(61,1){\makebox(0,0){$3$}}
\put(61,3){\makebox(0,0){$2$}}
\put(61,7){\makebox(0,0){$1$}}
\put(65,7){\makebox(0,0){$6$}}
\put(65,3){\makebox(0,0){$5$}}
\put(65,1){\makebox(0,0){$4$}}
\put(67,4){\makebox(0,0){$+$}}
\end{picture}

\begin{picture}(54,8)(13,22)
\put(19,4){\makebox(0,0){$+$}}

%1
\put(21,8){\vector(1,-1){2}}
\put(23,6){\vector(1,1){2}}
\put(21,0){\vector(1,1){2}}
\put(21,4){\vector(1,-1){2}}
\put(23.25,2){\circle*{.5}}
\put(23.25,2){\vector(1,1){2}}
\put(23.25,2){\vector(1,-1){2}}
\put(21,1){\makebox(0,0){$3$}}
\put(21,3){\makebox(0,0){$2$}}
\put(21,7){\makebox(0,0){$1$}}
\put(25,7){\makebox(0,0){$4$}}
\put(25.3,3){\makebox(0,0){$5$}}
\put(25.3,1){\makebox(0,0){$6$}}
\put(27,4){\makebox(0,0){$-$}}

%2
\put(29,8){\vector(1,-1){2}}
\put(31,6){\vector(1,1){2}}
\put(29,0){\vector(1,1){2}}
\put(29,4){\vector(1,-1){2}}
\put(31.25,2){\circle*{.5}}
\put(31.25,2){\vector(1,1){2}}
\put(31.25,2){\vector(1,-1){2}}
\put(29,1){\makebox(0,0){$3$}}
\put(29,3){\makebox(0,0){$2$}}
\put(29,7){\makebox(0,0){$1$}}
\put(33,7){\makebox(0,0){$5$}}
\put(33.3,3){\makebox(0,0){$4$}}
\put(33.3,1){\makebox(0,0){$6$}}
\put(35,4){\makebox(0,0){$+$}}

%3
\put(37,8){\vector(1,-1){2}}
\put(39,6){\vector(1,1){2}}
\put(37,0){\vector(1,1){2}}
\put(37,4){\vector(1,-1){2}}
\put(39.25,2){\circle*{.5}}
\put(39.25,2){\vector(1,1){2}}
\put(39.25,2){\vector(1,-1){2}}
\put(37,1){\makebox(0,0){$3$}}
\put(37,3){\makebox(0,0){$2$}}
\put(37,7){\makebox(0,0){$1$}}
\put(41,7){\makebox(0,0){$6$}}
\put(41.3,3){\makebox(0,0){$4$}}
\put(41.3,1){\makebox(0,0){$5$}}
\put(43,4){\makebox(0,0){$-$}}

%4
\put(45,8){\vector(1,-1){2}}
\put(47,6){\vector(1,1){2}}
\put(45,0){\vector(1,1){2}}
\put(45,4){\vector(1,-1){2}}
\put(47.25,2){\circle*{.5}}
\put(47.25,2){\vector(1,1){2}}
\put(47.25,2){\vector(1,-1){2}}
\put(45,1){\makebox(0,0){$3$}}
\put(45,3){\makebox(0,0){$1$}}
\put(45,7){\makebox(0,0){$2$}}
\put(49,7){\makebox(0,0){$4$}}
\put(49.3,3){\makebox(0,0){$5$}}
\put(49.3,1){\makebox(0,0){$6$}}
\put(51,4){\makebox(0,0){$+$}}

%5
\put(53,8){\vector(1,-1){2}}
\put(55,6){\vector(1,1){2}}
\put(53,0){\vector(1,1){2}}
\put(53,4){\vector(1,-1){2}}
\put(55.25,2){\circle*{.5}}
\put(55.25,2){\vector(1,1){2}}
\put(55.25,2){\vector(1,-1){2}}
\put(53,1){\makebox(0,0){$3$}}
\put(53,3){\makebox(0,0){$1$}}
\put(53,7){\makebox(0,0){$2$}}
\put(57,7){\makebox(0,0){$5$}}
\put(57.3,3){\makebox(0,0){$4$}}
\put(57.3,1){\makebox(0,0){$6$}}
\put(59,4){\makebox(0,0){$-$}}
\end{picture}


\begin{picture}(54,8)(13,28)
\put(19,4){\makebox(0,0){$-$}}
%1
\put(21,8){\vector(1,-1){2}}
\put(23,6){\vector(1,1){2}}
\put(21,0){\vector(1,1){2}}
\put(21,4){\vector(1,-1){2}}
\put(23.25,2){\circle*{.5}}
\put(23.25,2){\vector(1,1){2}}
\put(23.25,2){\vector(1,-1){2}}
\put(21,1){\makebox(0,0){$3$}}
\put(21,3){\makebox(0,0){$1$}}
\put(21,7){\makebox(0,0){$2$}}
\put(25,7){\makebox(0,0){$6$}}
\put(25.3,3){\makebox(0,0){$4$}}
\put(25.3,1){\makebox(0,0){$5$}}
\put(27,4){\makebox(0,0){$+$}}

%2
\put(29,8){\vector(1,-1){2}}
\put(31,6){\vector(1,1){2}}
\put(29,0){\vector(1,1){2}}
\put(29,4){\vector(1,-1){2}}
\put(31.25,2){\circle*{.5}}
\put(31.25,2){\vector(1,1){2}}
\put(31.25,2){\vector(1,-1){2}}
\put(29,1){\makebox(0,0){$2$}}
\put(29,3){\makebox(0,0){$1$}}
\put(29,7){\makebox(0,0){$3$}}
\put(33,7){\makebox(0,0){$4$}}
\put(33.3,3){\makebox(0,0){$5$}}
\put(33.3,1){\makebox(0,0){$6$}}
\put(35,4){\makebox(0,0){$-$}}


%3
\put(37,8){\vector(1,-1){2}}
\put(39,6){\vector(1,1){2}}
\put(37,0){\vector(1,1){2}}
\put(37,4){\vector(1,-1){2}}
\put(39.25,2){\circle*{.5}}
\put(39.25,2){\vector(1,1){2}}
\put(39.25,2){\vector(1,-1){2}}
\put(37,1){\makebox(0,0){$2$}}
\put(37,3){\makebox(0,0){$1$}}
\put(37,7){\makebox(0,0){$3$}}
\put(41,7){\makebox(0,0){$5$}}
\put(41.3,3){\makebox(0,0){$4$}}
\put(41.3,1){\makebox(0,0){$6$}}
\put(43,4){\makebox(0,0){$+$}}

%4
\put(45,8){\vector(1,-1){2}}
\put(47,6){\vector(1,1){2}}
\put(45,0){\vector(1,1){2}}
\put(45,4){\vector(1,-1){2}}
\put(47.25,2){\circle*{.5}}
\put(47.25,2){\vector(1,1){2}}
\put(47.25,2){\vector(1,-1){2}}
\put(45,1){\makebox(0,0){$2$}}
\put(45,3){\makebox(0,0){$1$}}
\put(45,7){\makebox(0,0){$3$}}
\put(49,7){\makebox(0,0){$6$}}
\put(49.3,3){\makebox(0,0){$4$}}
\put(49.3,1){\makebox(0,0){$5$}}
\put(51,4){\makebox(0,0){$+$}}

%5
\put(53,4){\vector(1,0){3}}
\put(56,4){\vector(1,0){3}}
\put(53,1){\vector(1,1){3}}
\put(53,7){\vector(1,-1){3}}
\put(56.25,4){\circle*{.5}}
\put(56.25,4){\vector(1,1){3}}
\put(56.25,4){\vector(1,-1){3}}
\put(53,0){\makebox(0,0){$2$}}
\put(53,3){\makebox(0,0){$1$}}
\put(54,7){\makebox(0,0){$3$}}
\put(58,7){\makebox(0,0){$6$}}
\put(59.3,3){\makebox(0,0){$4$}}
\put(59.3,0){\makebox(0,0){$5$}}
\end{picture}
\vskip 6.5  cm Again, because of the commutation relations
(\ref{commutation}), if we interchange two indices corresponding to
two creation or two annihilation operators, the average should change
its sign, for example: $N_{123456}=-N_{213456}$. The picture below
illustrates a simple algorithm of counting the parity of
permutation by counting the number of crossing between lines
connecting different arguments. In the example below, one sees that
the parity of $n_1\d^1_5\P_{2346}$ term in the expansion of $N_{123456}$ is
odd (because of the odd number of crossings), so the product is
negative.

\setlength{\unitlength}{1.0cm}
\begin{picture}(54,3)(0,0)
%\put(39.25,3){\vector(1,-1){3}}
%\put(42,1.5){\makebox(0,0){$2$}}
\put(1,-0.3){\makebox(0,0){$1$}}
\put(2,-0.3){\makebox(0,0){$2$}}
\put(3,-0.3){\makebox(0,0){$3$}}
\put(4,-0.3){\makebox(0,0){$4$}}
\put(5,-0.3){\makebox(0,0){$5$}}
\put(6,-0.3){\makebox(0,0){$6$}}
\put(3,2){\line(-1,-2){1}}
\put(3,2){\line( 3,-2){3}}
\put(3,2){\line( 1,-2){1}}
\put(3,2){\line( 0,-2){2}}
\put(1,1){\line( 0,-1){1}}
\put(1,1){\line( 3, 0){4}}
\put(5,1){\line( 0,-1){1}}
\end{picture}
\setlength{\unitlength}{.22cm}
\vskip 3cm
\noindent
In the  same manner, the definition of $\Q_{12345678}$
\BEA N_{12345678}&=&\r_1\r_2\r_3\r_4
(\d^{4}_{5}\d^{3}_{6}\d^{2}_{7}\d^{1}_{8}+...)\cr
&&\r_4\Q_{123678}\d^{4}_{5}+...\ + \r_1 \r_2 \P_{3456}\d^{1}_{8}\d^{2}_{7}+..
.+
\P_{3456} \P_{1278}+ ...\ \cr
& &+\R_{12345678}\label{dN8def}
\EEA
can be presented as
\
\vskip 1 cm
%factorization of N(12345678)
\setlength{\unitlength}{.22cm}
\begin{picture}(54,10)(0,0)
\put(0,2){\vector(1,0){6}}
\put(0,4){\vector(1,0){6}}
\put(0,6){\vector(1,0){6}}
\put(0,8){\vector(1,0){6}}
\put(6,0){\line(1,0){6}}
\put(6,10){\line(1,0){6}}
\put(6,0){\line(0,1){10}}
\put(12,0){\line(0,1){10}}
\put(12,2){\vector(1,0){6}}
\put(12,4){\vector(1,0){6}}
\put(12,6){\vector(1,0){6}}
\put(12,8){\vector(1,0){6}}
\put(1,3){\makebox(0,0){$4$}}
\put(1,5){\makebox(0,0){$3$}}
\put(1,7){\makebox(0,0){$2$}}
\put(1,9){\makebox(0,0){$1$}}
\put(17,9){\makebox(0,0){$5$}}
\put(17,7){\makebox(0,0){$6$}}
\put(17,5){\makebox(0,0){$7$}}
\put(17,3){\makebox(0,0){$8$}}
\put(19,5){\makebox(0,0){$=$}}
%1
\put(21,0){\vector(1,1){2}}
\put(21,10){\vector(1,-1){2}}
\put(21,4){\vector(1,0){2}}
\put(23,4){\vector(1,0){2}}
\put(21,6){\vector(1,0){2}}
\put(23,6){\vector(1,0){2}}
\put(23,8){\vector(1,1){2}}
\put(23,2){\vector(1,-1){2}}
\put(27,5){\makebox(0,0){$+$}}
%2
\put(32,7){\vector(1,1){3}}
\put(29,10){\vector(1,-1){3}}
\put(29,3){\vector(1,0){3}}
\put(32,3){\vector(1,0){3}}
\put(29,0){\vector(1,1){3}}
\put(29,6){\vector(1,-1){3}}
\put(32.25,3){\circle*{.5}}
\put(32.25,3){\vector(1,1){3}}
\put(32.25,3){\vector(1,-1){3}}
\put(37,5){\makebox(0,0){$+$}}
%3

\put(39,0){\vector(1,1){2}}
\put(39,4){\vector(1,-1){2}}
\put(41.25,2){\circle*{.5}}
\put(41.25,2){\vector(1,1){2}}
\put(41.25,2){\vector(1,-1){2}}

\put(39,5){\vector(1,1){2}}
\put(41,7){\vector(1,-1){2}}

\put(39,10){\vector(1,-1){2}}
\put(41,8){\vector(1,1){2}}

%\put(39,6){\vector(1,1){2}}
%\put(39,10){\vector(1,-1){2}}
%\put(41.25,8){\circle*{.5}}
%\put(41.25,8){\vector(1,1){2}}
%\put(41.25,8){\vector(1,-1){2}}

\put(45,5){\makebox(0,0){$+$}}

%4
\put(47,0){\vector(1,1){2}}
\put(47,4){\vector(1,-1){2}}
\put(49.25,2){\circle*{.5}}
\put(49.25,2){\vector(1,1){2}}
\put(49.25,2){\vector(1,-1){2}}
\put(47,6){\vector(1,1){2}}
\put(47,10){\vector(1,-1){2}}
\put(49.25,8){\circle*{.5}}
\put(49.25,8){\vector(1,1){2}}
\put(49.25,8){\vector(1,-1){2}}
\put(53,5){\makebox(0,0){$+$}}
%5
\put(55,0){\line(1,1){10}}
\put(55,10){\line(1,-1){10}}
\put(55,2){\line(5,3){10}}
\put(55,8){\line(5,-3){10}}
\put(60,5){\circle*{.5}}
\put(55,11){\makebox(0,0){$1$}} \put(55,7){\makebox(0,0){$2$}}
\put(55,3){\makebox(0,0){$3$}} \put(55,-1){\makebox(0,0){$4$}}
\put(65,11){\makebox(0,0){$5$}} \put(65,7){\makebox(0,0){$6$}}
\put(65,3){\makebox(0,0){$7$}} \put(65,-1){\makebox(0,0){$8$}}
\end{picture} 
%***************************************************************
\vskip 1cm
Because of the large amount of terms in this case, we show only
schematically the
factorization of $N_{1'2'3'4'5'6'7'8'}$. One has to choose {\it all}
possible permutations of indices, corresponding to different
topologies putting ``annihilation'' arguments to the incoming and
``creation'' arguments to the outgoing arrows.
\begin{thebibliography}{100} 
\bibitem{b1}
L. D. Landau and E. M. Lifshitz, {\em Physical Kinetics}
 (Pergamon Press, New York, 1981).
\bibitem{b2}
L. P. Kadanoff and G. Baym, {\em Quantum Statistical Mechanics}
  (Addison-Wesley, New York, 1989).
\bibitem{b3}
P. P. J. M. Schram, {\em Kinetic Theory of Gases and Plasmas}
(Kluwer Academic Press, Dordrecht, 1991).
 \bibitem{b4}
D. F. Dubois,  in {\em Lectures in Theoretical Physics}, edited by W. E.
  Brittin (Gordon and Breach, New York, 1967), Vol. IX c, pp.\ 469--620.
\bibitem{b5}
Y. L. Klimontovich, D. Kremp, and W. D. Kraeft, Advances of Chemical
Phys. {\bf  68},  175  (1987).
\bibitem{b6}
P. Danielewicz, Ann. of Physics {\bf 152},  239  (1984).
\bibitem{b7} A. Balescu, Physics of Fluids, {\bf 3},52, (1960).
\bibitem{b8} A. Lenard, Annals of Physics (N.Y.), {\bf 3}, 390, (1960).
\bibitem{b9} W. D. Kraeft, D. Kremp, W. Ebeling, and G. R{\"{o}}pke,
{\em Quantum statistics of charged particle systems} (Akademie-Verlag,
Berlin, 1986).
\bibitem{Kadanoff} L.P.  Kadanoff and G.  
Baym, "Quantum Statistical Mechanics, Green's Function methods in 
equilibrium problems.
\bibitem {koehler} H.S.K{\"o}hler, Phys. Rev. E {\bf 53}, 3145, (1996).
\bibitem{Lipavsky} P.Lipavski, V.Spica, Phys Rev B.
{\bf 52}, 14615,(1995), \ \ \ \ Phys. Rev. Lett. {\bf 73}, 3439, (1994) 
\bibitem{elsayed} K.E. Sayed, L.Banyai and H.Haug, Phys.Rev. B {\bf
50} 1541 (1994).
\bibitem{kuznetsov} A.V.Kuznetsov, Phys. Rev. B {\bf 44} 8721 (1991).
\bibitem{ZLF}V.E.  Zakharov, V.S.  L'vov 
and G.Falkovich, "Kolmogorov Spectra of Turbulence", Springer-Verlag, 1992.
\bibitem{Optical} S.  Dyachenko, A.C. Newell, A.  Pushkarev, V.E.  
Zakharov, Physica D {\bf 57}, 96-160, (1992).
\bibitem{HASS}
S. Hasselmann, K. Hasselmann, J.H. Allender and T.P. Barnett 
J.Phys.Oceonography, {\bf 15}, 1378, (1985).
\bibitem{Balk} A.M.Balk, J. of Fluid Mechanics, {\bf 315}, 139, (1996).
\bibitem{Newell} B.  J . Benney and A.C.  Newell, Studies in Appl.  
Math.  {\bf 48} (1) 29 (1969).  
\bibitem{Newell1} A.C. Newell, Review of 
Geophysics, {\bf 6 }, (1968a), 1.
\bibitem{Ben} D.J. Benney and P.Saffman, Proc Royal. Soc, A(1966), 289, 
301-320.
\bibitem{Z68a} V.E.Zakharov, Zh. Priklad. Tech. Fiz, 1, 35 (1965).
J.Appl.Mech. Tech. Phys, {\bf 1}, 22, (1967) 
\bibitem{Z68b} V.E.Zakharov Zh.Eksp.Teor.Fiz. {\bf 51},688,(1966) 
Sov.Phys.JETP {\bf 24} (1967).
\bibitem{Koch} W.W. Chow, S.W. Koch, 
M.  Sargent, "Semiconductor Laser Physics", Springer-Verlag, 1994.
\bibitem{bhattacharya.94} P.Bhattacharya {\em Semiconductor
Optoelectronic Devices}(Prentice Hall, Englewood Cliff, 1994)
\bibitem{zhang-etal.96} X.Zhang, Y.Yuan, P.Bhattacharya, Applied
Physics Letters {\bf 69},2309 (1996).
\bibitem{Koch1} H.Haug and S.W. Koch  ``Quantum theory of the
optical and electronic properties of the semiconductors''. 2nd edition, World
Scientific, Singapore 1993.

%\bibitem {Balk1} Balk, Pushkarev {private 
%communications-1991}.  
%\bibitem{Binder1} M.  Linberg, Y.Z.Hu, R.Binder 
%and S.W.Koch, Phys Rev B (50) 1994 

\end{thebibliography} 
\section{Figures}
\begin{figure}
\epsfxsize=17cm 
\epsffile{FinalPictures.eps.0}\epsfxsize=17cm 
\end{figure}
\epsffile{FinalPictures.eps.1}\epsfxsize=17cm 

\epsffile{FinalPictures.eps.2}\epsfxsize=17cm 

\epsffile{FinalPictures.eps.3}\epsfxsize=17cm 

\epsffile{FinalPictures.eps.4}\epsfxsize=17cm 

\epsffile{FinalPictures.eps.5}\epsfxsize=17cm 

\epsffile{FinalPictures.eps.6}\epsfxsize=17cm 

\epsffile{FinalPictures.eps.7}\epsfxsize=17cm 

\epsffile{FinalPictures.eps.8}\epsfxsize=17cm 

\epsffile{FinalPictures.eps.9}


\end{document}

