\documentstyle [12pt] {article}


% chamlin, 2/20/92. comment out next 2 lines and select eqsecnum in REVTeX runs
\renewcommand{\theequation}{\arabic{section}.\arabic{equation}}
\makeatletter\@addtoreset{equation}{section}\makeatother

% --- Definitions for equation numberings (section #.equation #)
% --- Use \mysection rather than \section
% --- Use \myappendix rather than \appendix
% --- Use \myappsection rather than \section (for appendix sections only!)
% --- End of definitions



\begin{document}

\rightline{TAUP 1944-92}
\rightline{February, 1992} \vskip .5in
%\baselineskip=21pt

{%\baselineskip=14pt
\centerline{\large \bf FERMION PAIR PRODUCTION IN A STRONG ELECTRIC FIELD}
\vskip .5in
\centerline{Y. Kluger, J. M. Eisenberg, and B. Svetitsky}
\vskip .25in
{\it
\centerline{School of Physics and Astronomy}
\centerline{Raymond and Beverly Sackler Faculty}
\centerline{of Exact Sciences}
\centerline{Tel Aviv University}
\centerline{69978 Tel Aviv, Israel}  }
\vskip .4in
\centerline{F. Cooper and E. Mottola}
\vskip .25in
{\it
\centerline{Theoretical Division}
\centerline{Los Alamos National Laboratory}
\centerline{Los Alamos, New Mexico 87545 USA}  }

\vskip 0.5in
\noindent {\it ABSTRACT.\/}
The initial-value problem for the quantum back-reaction
in spinor QED is formulated and solved in the semiclassical mean field approximation, for a homogeneous but time-dependent electric field $E(t)$.
We apply the method of adiabatic regularization to the Dirac equation in
order to renormalize the expectation value of the current and derive a
finite coupled set of ordinary differential equations for the
time evolution of the system.
We solve this system in (1+1) dimensions numerically and compare
the solution to a simple model based on a relativistic Boltzmann-Vlasov
equation, with a particle
production source term inferred from the Schwinger particle creation rate
and a Pauli-blocking factor. This model reproduces very well the time
behavior of the electric field and the creation rate of $e^+e^-$ pairs
of the semiclassical calculation.   

} % end of \baselineskip 14pt
\newpage

\section {Introduction}
 The rate of creation of pairs of charged particles in a static
 and homogeneous external electric field was computed long ago [1--6].
This process  has been used
extensively in color-flux-tube models to describe multiparticle production 
in hadronic collisions [6--13]. A strong color-electric field is assumed to 
be formed  between receding 
hadronic sources---quarks or hadrons---and
quark-antiquark and gluon pairs
emerge in the presence of the field by tunneling.
In models for the production of a quark-gluon
plasma the source for the field is the fragmenting nuclei.

In solving a dynamical problem with a strong initial  electric field the effect 
of the produced particles on the electric field (the back-reaction)
should be taken into consideration [11--13]. 
The construction of the usual models for quark-gluon plasma creation has often been based on oversimplifications.  One is the modification of the Schwinger expression for the time independent rate of pair production 
so that it becomes time dependent through the time variation of the 
electric field (even though a fixed external electric field is inherent in Schwinger's derivation);  moreover the transverse and
longitudinal momentum distributions of the produced particles are chosen
according to results of a WKB calculation [12] . 
It is not clear whether this choice is entirely appropriate.

In order to evaluate the dynamic rate of pair production self-consistently,
we propose to study the  quantum back-reaction of 
spin-$\frac {1}{2}$ fields in a spatially homogeneous classical electric field 
through the semiclassical Maxwell equations.  
The back-reaction problem for quantum fields has been developed particularly
in the study of quantum fields in curved space [14--19]. 
In formulating the back-reaction problem,  divergences appear in the 
expectation value of the conserved currents, and  a renormalization procedure
which guarantees well-defined finite equations is required.
Adiabatic regularization is a useful approach which enables one to 
dispose of these infinities in a way that is consistent with  conventional 
renormalization
for a variety of spatially homogeneous problems [14--20], and is easy to 
implement in a practical numerical procedure. 
This method has been successfully applied to the study of the back-reaction
problems in which the matter field is a scalar field [17--20].

In Ref.~[20] we found that results based on calculations in  
semiclassical scalar electrodynamics are very similar to those obtained
from a model based on a phenomenological relativistic Boltzmann-Vlasov equation.
It is of interest to see if the same results obtain for fermions which,
unlike bosons, possess no classical limit, and for which Pauli blocking, 
rather than Bose enhancement, implies a quite different consequence from
the consideration of quantum statistics.

To apply adiabatic regularization to the case of fermions, we again express 
the Fourier components of the field operators in a WKB-like form. 
This enables us to isolate the ultraviolet divergences via an adiabatic 
expansion and to perform a mode-by-mode subtraction. 
In Section 2 we derive the coupled
equations for the fields in the semiclassical limit of QED, and in
Section 3 we discuss the adiabatic regularization procedure.
In Section 4 we present our numerical results in (1+1) dimensions 
and compare them with the phenomenological Boltzmann-Vlasov model.




\newpage
\section {QED in the semiclassical mean field limit}

The lagrangian density for electrodynamics is 
\begin{equation}
{\cal L }=\bar {\psi} i\gamma^{\mu}
(\partial_{\mu} +ieA_{\mu})\psi -m\bar {\psi}\psi-
\frac {1}{4}F_{\mu \nu}F^{\mu \nu},
\end{equation}
where the metric convention is taken as $(+---)$.
For the $\gamma$ matrices we use the  convention of Bjorken and Drell [21]
\begin{equation}
\begin{array}{lr}
\gamma^{0}=\left(\begin{array}{cc}
            I & 0\\
            0 & -I \end{array} \right)            

            & \gamma^{i}=\left(\begin{array}{cc}
             0 & \sigma^{i}\\
            -\sigma^{i} & 0 \end{array} \right),

\end{array}
\quad i=1,2,3
\end{equation}
where $I$ is the identity matrix and $\sigma^{i}$ are the Pauli matrices.

In the mean field approximation we quantize only the Dirac field, while the
electromagnetic field is treated classically. This approximation may be
derived formally as the leading term in the large $N$ limit of QED, where
$N$ is the number of charged matter fields [19]. 
The resulting coupled field equations read
\begin{eqnarray}
(i\gamma^{\mu}\partial_{\mu}-e\gamma^{\mu}A_{\mu}-m)\psi (x) 
= 0 \\ 
\partial_{\mu} F^{\mu \nu }=\langle j^{\nu} \rangle =
\frac{e}{2}\langle [\bar {\psi},\gamma ^{\nu}\psi]\rangle 
\end{eqnarray}
where the expectation value is with respect to the initial state of
the spinor field. The commutator in the electric 
current guarantees a zero
expectation value for any charge conjugation eigenstate. 
Expressing the solution of the Dirac equation as
\begin{equation}
\psi (x)=(i\gamma^{\mu}\partial_{\mu}-e\gamma^{\mu}A_{\mu}+m)\phi (x) ,
\end{equation}
and inserting (2.5) into (2.3), it follows that $\phi$ satisfies 
the quadratic Dirac equation 
\begin{equation} 
\left [ (i\partial_{\mu} -eA_{\mu})^2-\frac{e}{2}\sigma^{\mu \nu}F_{\mu \nu} 
-m^2 \right ]\phi (x)=0 ,
\end{equation}
where $\phi $ is a four component spinor.
Here we consider the case where the electric field is spatially 
homogeneous so that
the field strength $F^{\mu \nu}$  depends only on time.  Owing to  homogeneity
the semiclassical Maxwell equations (2.4) allow only configurations
where $\langle j^0 \rangle =0$. We take the electric field to be in the  
direction of the
$z$ axis, and we choose a gauge such that only $A\equiv A^3 (t)$ 
is nonvanishing.    
Then the second order Dirac equation becomes
\begin{equation} 
\left [ \Box +e^2 A^2 (t)+2iA(t)\partial_3-ie\partial_0 A(t)\gamma^0 \gamma^3
+m^2 \right ]\phi (x)=0.
\end{equation}


Spatial homogeneity implies that there exist solutions of the form
\begin{equation}
\phi_{{\bf{k}}s}(x) = e^{i{\bf k \cdot x}} f_{{\bf{k}}s}(t) \chi_s,
\end{equation} 
where 
\begin{equation}
\begin{array}{lr}
\chi_1=\left(\begin{array}{cc}
            1 \\
            0 \\
            1 \\
            0 \end{array} \right),            

            & \chi_2=\left(\begin{array}{cc}
             0 \\
             1 \\
             0 \\
            -1 \end{array} \right),
\end{array}
\end{equation}
\begin{equation}
\begin{array}{lr} 

              \chi_3=\left(\begin{array}{cc}
            1 \\
            0 \\
           -1 \\
            0 \end{array} \right),            

            & \chi_4=\left(\begin{array}{cc}
             0 \\
             1 \\
             0 \\
             1 \end{array} \right). 

\end{array}
\end{equation}
These spinors are chosen to be eigenvectors of $\alpha^3=\gamma^0\gamma^3$
in the representation (2.2) for the $\gamma$ matrices. The eigenvalues of 
$\alpha^3$ are 
$\lambda_{s}=1$ for $s=1,2$ and  $\lambda_{s}=-1$ for $s=3,4$.
These spinors satisfy the normalization and completeness
conditions,
\begin{eqnarray}
\sum_{\alpha =1}^{4} (\chi^{\dag}_r)^{\alpha}(\chi_s)_{\alpha} &=& 2 \delta_{rs}
\nonumber \\
\sum_{r =1}^{4} (\chi^{\dag}_r)^{\alpha}(\chi_r)_{\beta} &=& 
2 \delta^{\alpha}_{\ \beta}\ .
\end{eqnarray}
Substituting (2.8) into (2.7) it follows that the mode functions 
$f_{{\bf{k}}s}(t)$  
satisfy 
\begin{equation}
\frac{d^2 f_{{\bf{k}}s}(t)}{dt^2} +
\left [{\omega _{\bf{k}}^2 (t)-i\lambda_{s}e\frac{dA}{dt}}\right ]
f_{{\bf{k}}s}(t)=0 ,
\end{equation} 
with
\begin{equation}
\omega^2 _{\bf{k}} (t) \equiv p_{3}^{2} +{\bf{k}}_{\perp}^2+m^2  \quad 
,\quad  
{\bf{k}}_{\perp}^2\equiv k_1^2+k_2^2 ,\quad p_i \equiv k^i-eA^i(t) \quad.
\end{equation}


Equations (2.12)  are second order differential equations, and therefore
for each $s$ there are two independent solutions. 
Let $f^{+}_{{\bf{k}}s}$ and $f^{-}_{{\bf{k}}s}$ be the two independent
solutions of (2.12),
which become positive- and negative-frequency  solutions 
in the absence of the electric field.
Clearly at the moment we have eight different solutions for the second-order
equation  (2.6), namely $f^{\pm}_{{\bf{k}}s}$ for $s=1,2,3,4$.
However, the Dirac equation (2.3) has only four independent solutions . If we
restrict ourselves to solutions which belong to the set $s=1,2$ or
to the set $s=3,4$ we shall see that from each set one can construct 
a linearly independent set of solutions of the Dirac equation.
The form introduced in (2.8) allows us to write $\psi$ as
\begin{equation}
\psi_{{\bf{k}}s}(x)\equiv \hat{D}\phi_{{\bf{k}}s} (x)
\equiv (i\gamma^{0}\partial_{0}+{\gamma}^i k_i
-e\gamma^{3} A_3+m)\phi_{{\bf{k}}s} (x).
\end{equation}
Explicitly, the two sets of independent solutions of the Dirac equation
may be taken to be
 \begin{eqnarray}
\psi^{\pm}_{{\bf{k}}s}= 
e^{i{\bf{k\cdot x}}}\hat{D}f^{\pm}_{{\bf{k}}s}\chi_s, \quad s=1,2
\end{eqnarray}
and
\begin{eqnarray}
\psi^{\pm}_{{\bf{k}}s}= 
e^{i{\bf{k\cdot x}}}\hat{D}f^{\pm}_{{\bf{k}}s}\chi_s, \quad s=3,4 \ .
\end{eqnarray}
Using Eqs.~(2.14)--(2.16) we find (for a given ${\bf{k}}$)
\begin{eqnarray}
{\psi^{\pm}_{r}}^{\dagger} {\psi^{\pm}_{s}} &=&
{\chi_{r}}^{\dagger}{f_{r}^{\ast \pm }}{\hat{D}}^{\dagger}
\hat{D}f_{s}^{\pm}\chi_{s}  \\  \nonumber
&=&2 \left \{ \omega^2 _{\bf{k}} (t) {f_{r}^{\ast \pm }} 
f_{s}^{\pm}+
{{\dot{f}}_{r}^{\ast \pm }}
\dot{f}_{s}^{\pm}-i\lambda_{s}
p_3 \left [{f_{r}^{\ast \pm }}
\dot{f}_{s}^{\pm}-
{\dot{f}_{r}^{\ast \pm }} f_{s}^{\pm}\right ] \right \}\delta_{rs} 
\end{eqnarray}
where either $r,s =1,2$ or $r,s=3,4$. An exactly analogous formula may be
derived for ${\psi^{\pm}_{r}}^{\dagger} {\psi^{\mp}_{s}}$.
By differentiating these expressions with respect to time and by using
Eq.~(2.12), it can be readily verified that these inner products are
time-independent.
As $t \rightarrow -\infty$ there is no
interaction between the fermion field and  the electromagnetic field, and    
we can choose two
independent plane wave solutions for Eq.~(2.12),
\begin{eqnarray}
\lim_{t \rightarrow -\infty}f^{\pm}_{{\bf{k}}s}=
c_{s} e^{\mp i\omega_{\bf{k}}t}, 
\end{eqnarray}
where $c_s$ are constants.
Insertion of these free solutions in the relation for
${\psi^{\pm}_{r}}^{\dagger} {\psi^{\mp}_{s}}$ yields immediately 

\begin{equation}
 {\psi^{\pm}_{r}}^{\dagger}{\psi^{\mp}_{s}}=0 \ .
\end{equation}
Since this result is time-independent, it is valid at any time, and 
each set $\psi^{\pm}_{s}$, with $s=1,2$ or $s=3,4$,
is a complete set of 
linearly independent solutions of the Dirac equation.
Note that these complete
systems are not identical, and orthonormality conditions holds 
for each set separately.
 In principle, we need only one of these sets in order to expand the 
field operator $\Psi$ in terms of single-particle solutions.
In order to ensure that with our initial conditions the Dirac current 
vanishes at $t=0$ it is advantageous to use {\it both}
sets in our calculations, as will be shown in the next section.   


Now we can construct the quantized spinor 
field operator in the form
\begin{eqnarray}
\Psi (x)&=&\int [d{\bf{k}}]\sum_{s =1,2}[b_{s}({\bf{k}})
{\psi^{+}_{{\bf{k}}s}}
+d_{s}^{\dagger}({\bf{-k}}){\psi^{-}_{{\bf{k}}s}} ]  
\\ \nonumber
&= &
\int [d{\bf{k}}]\sum_{s =3,4}[b_{s}({\bf{k}})
{\psi^{+}_{{\bf{k}}s}}
+d_{s}^{\dagger}({\bf{-k}}){\psi^{-}_{{\bf{k}}s}} ],
\end{eqnarray}
where the two lines show the field expressed in terms of the two alternative 
bases. 
Here $[d{\bf{k}}]={d^3{\bf{k}}}/{(2\pi)^3}$. The  fermion fields obey
canonical anticommutation relations: 
$\{ \Psi_{\alpha} (t,{\bf{x}}),\Psi^{\dagger} _{\beta} (t,{\bf{y}})\}=
\delta^3 ({\bf{x}}-{\bf{y}})\delta_{\alpha \beta}$.
The creation and annihilation operators of each set
($r,s=1,2$ or $r,s=3,4$) will obey the standard anticommutation relations,
\begin{eqnarray}
 \{ b_{r}({\bf{k}}),b_{s}^{\dagger}({\bf{q}}) \}=
 \{ d_{r}({\bf{k}}),d_{s}^{\dagger}({\bf{q}}) \}=
(2\pi )^3\delta^3({\bf{k}}-{\bf{q}})\delta_{rs}  ,
\end{eqnarray}
if we impose the normalization condition
\begin{eqnarray}
 {\psi^{\pm}_{r}}^{\dagger}{\psi^{\pm}_{s}}= \delta_{rs},
\end{eqnarray} 
which fixes the normalization of the mode functions.

We turn now to the calculation of the expectation value of the electric current.
For the sake of simplicity we choose the initial state to be the vacuum.
Using the anticommutation relations (2.21) we find
\begin{eqnarray}
\langle 0\vert j^3\vert 0 \rangle &=&\frac{e}{2}
\langle 0\vert [\bar{\Psi},\gamma^3\Psi]\vert 0\rangle \\ \nonumber
&=&\frac{e}{2}\int [d{\bf{k}}]\sum_{s=1,2}\left \{
-{\psi^+_{{\bf{k}}s}}^{\dagger}\gamma^0\gamma^3
{\psi^+_{{\bf{k}}s}}
+{\psi^-_{{\bf{k}}s}}^{\dagger}\gamma^0\gamma^3
{\psi^-_{{\bf{k}}s}} \right \} \ .
\end{eqnarray}
Alternatively,
\begin{eqnarray}
\langle 0\vert j^3\vert 0\rangle =\frac{e}{2}\int [d{\bf{k}}]
\sum_{s=3,4}\left \{-{\psi^+_{{\bf{k}}s}}^{\dagger}\gamma^0\gamma^3
{\psi^+_{{\bf{k}}s}}
+{\psi^-_{{\bf{k}}s}}^{\dagger}\gamma^0\gamma^3
{\psi^-_{{\bf{k}}s}}\right  \} 
\end{eqnarray}
Averaging the two expressions,
\begin{eqnarray}
\langle 0\vert j^3\vert 0 \rangle =\frac{e}{4}\int [d{\bf{k}}]
\sum_{s=1}^{4}\left \{-{\psi^+_{{\bf{k}}s}}^{\dagger}\gamma^0\gamma^3
{\psi^+_{{\bf{k}}s}}
+{\psi^-_{{\bf{k}}s}}^{\dagger}\gamma^0\gamma^3
{\psi^-_{{\bf{k}}s}}\right \}. 
\end{eqnarray}
This form will be useful 
when we turn to the adiabatic regularization in the next section.
The other components of the current are zero since the electric field 
is in the $z$ direction.
Using (2.14) we find
\begin{equation}
\gamma ^{0}\gamma ^{3} {\psi^{\pm}_{{\bf{k}}s}}
= \lambda_{s} \left [-i\gamma^{0}\partial_{0} +{\bf{\gamma}}^{\perp} 
{\bf{k}}_{\perp}
-(k_3-eA_3)\gamma^{3}+m \right ]{\phi^{\pm}_{{\bf{k}}s}} \quad , 
\end{equation}
and thus (2.17) and (2.22) give
\begin{eqnarray}
&&{\psi^{\pm}_{{\bf{k}}s}}^{\dagger}\gamma^0\gamma^3 
{\psi^{\pm}_{{\bf{k}}s}}=
\chi_{s}^{\dagger}f_{{\bf{k}}s}^{\ast \pm }{\hat{D}}^{\dagger}
\gamma^0\gamma^3 \hat{D} f_{{\bf{k}}s}^{\pm }\chi_{s}
 \nonumber \\
&=&2\lambda_{s} \left \{ (k^2_{\perp}+m^2-p_3{^2})
f_{{\bf{k}}s}^{\ast \pm }f_{{\bf{k}}s}^{\pm}-
\dot{f}_{{\bf{k}}s}^{\ast \pm }
\dot{f}_{{\bf{k}}s}^{\pm}+i \lambda_{s}
p_3(f_{{\bf{k}}s}^{\ast \pm }
\dot{f}_{{\bf{k}}s}^{\pm}-
\dot{f}_{{\bf{k}}s}^{\ast \pm } f_{{\bf{k}}s}^{\pm}) \right \} 
\nonumber \\
&=&\lambda_{s}\left 
[4(k^2_{\perp}+m^2)\vert f_{{\bf{k}}s}^{\pm}\vert ^2 -1 \right ]  ,
\end{eqnarray}
where the index $\alpha$ is not summed over. 
Inserting (2.27) into (2.24) gives
\begin{equation}
\langle 0\vert j^3\vert 0\rangle =e \sum_{s=1}^{4}\int [d{\bf{k}}](k^2_{\perp}+m
^2)
\lambda_{s} \left (\vert f_{{\bf{k}}s}^{-}\vert ^2 - 
\vert f_{{\bf{k}}s}^{+}\vert ^2\right ).
\end{equation}
From (2.17) and (2.22) it can be shown [22] that
\begin{equation}
2(k^2_{\perp} +m^2)\left (\vert f_{{\bf{k}}s}^{+}\vert ^2 + 
\vert f_{{\bf{k}}s}^{-}\vert ^2\right )=1 .
\end{equation}
Eq.~(2.29) then gives the current as
\begin{equation}
\langle 0\vert j^3\vert 0\rangle 
=-2e \sum_{s=1}^{4}\int [d{\bf{k}}](k^2_{\perp} +m^2) 
\lambda_{s}\vert f_{{\bf{k}}s}^{+}\vert ^2 .
\end{equation}

\newpage
\section{Adiabatic regularization}
The difficulty in solving the coupled semiclassical equations (2.3) and (2.4) 
originates in the fact that the expectation value of the current (2.30) 
diverges in the interacting theory.
 This infinity is related to charge renormalization.
 The adiabatic regularization approach is suitable for identifying and 
removing the ultraviolet divergences [19]. 
In order to isolate the ultraviolet behavior of the current integrand 
by adiabatic expansion, we need to express the mode 
equations (2.12) in a WKB-like form. The generic
problem is to find a suitable parameterization of the
solution of the differential equation  
$\ddot{u}(t) + \epsilon(t) u(t) = 0$, where in the present case
$\epsilon$ is the complex quantity in square brackets in (2.12). 
Such a parameterization was found in
[23], namely,

 \begin{eqnarray}
f_{{\bf{k}}s}^+ (t) = N_{{\bf{k}}s}
 \frac {1}{\sqrt{2\Omega_{{\bf{k}}s}}} \exp\left \{ \int_{0}^{t} 
\left ( -i\Omega_{{\bf{k}}s} (t^{\prime}) 
-\lambda_{s}  \frac {e\dot{A}(t^{\prime})}{2\Omega_{{\bf{k}}s} (t^{\prime})}
\right )
dt^{\prime}\right \}  , 
\end{eqnarray}
where $N_{{\bf{k}}s}$ are normalization constants and $\Omega_{{\bf{k}}s}$
is a real generalized frequency. (The second solution
$f_{{\bf{k}}s}^-$  for the mode equation can be found by using its Wronskian. 
When choosing the form (3.1) for $f_{{\bf{k}}s}^+$,
the second solution  does not 
 have a simple form, and for this
reason we expressed the current  (2.30) in terms of the positive-frequency
solutions only.)
By substituting (3.1) into (2.12) we obtain 
the WKB-like equation for $\Omega_{{\bf{k}}s}$,
\begin{eqnarray}
\Omega^2_{{\bf{k}}s}(t)= 
-\frac{\ddot{\Omega}_{{\bf{k}}s}}{2\Omega_{{\bf{k}}s}}
+\frac{3}{4}
\frac{\dot{\Omega}^2_{{\bf{k}}s}}{ \Omega^2_{{\bf{k}}s}}
+\omega^2_{\bf{k}}+\left( \frac{e\dot{A}}{2\Omega_{{\bf{k}}s}}\right )^2
-\lambda_{s} \frac{e\ddot{A}}{ 2\Omega_{{\bf{k}}s}} 
+\lambda_{s} \frac{e\dot{A}\dot{\Omega}_{{\bf{k}}s}}{ \Omega^2_{{\bf{k}}s}}  .
\end{eqnarray}
As in the bosonic case, the equation for $\Omega$ is a second order 
non-linear differential equation.

This equation enables us to study the large-momentum behavior of the
solutions.
It will be shown that an adiabatic expansion of (3.2)  to second order is needed to identify and isolate the divergences in the current (2.30).
The adiabatic expansion, an expansion in powers of $1/\omega_{\bf{k}}$ 
for large ${\bf{k}}$, is obtained by successive iteration: Inserting 
the zeroth order solution  $\Omega^{(0)}_{{\bf{k}}s}=\omega_{\bf{k}}$ 
into the right-hand side of (3.2) one obtains $\Omega^2$ up to second order;
inserting this value to the right-hand side the fourth order is obtained, 
and so forth.
It is not difficult to see that higher-order adiabatic approximations 
contain terms of higher order in $1/\omega_{\bf{k}}$. 

Noting that
\begin{equation}
\ddot{\omega}_{\bf{k}}=\frac {-e\ddot{A}p_3}{\omega_{\bf{k}}} 
+O(1/ \omega_{\bf{k}}) ,
\end{equation} 
we have up to second order
\begin{equation}
\Omega_{{\bf{k}}s}=
\omega_{\bf{k}}-e\ddot{A}\left (\lambda_{s} \omega_{\bf{k}}
- p_3 \right ) /4\omega^3_{\bf{k}}
+O(1/ \omega^3_{\bf{k}}),
\end{equation}
We shall refer to
\begin{equation}
\Omega^{(2),min}_{{\bf{k}}s}\equiv
\left [\omega_{\bf{k}}-e\ddot{A}\left (\lambda_{s} \omega_{\bf{k}}
- p_3 \right ) /4\omega^3_{\bf{k}}\right ] 
\end{equation}
as the minimal second order adiabatic approximation.
 Using the ansatz (3.1) the current reads
\begin{eqnarray}
\langle 0\vert j^3\vert 0 \rangle 
= -2e \sum_{s=1}^{4}\int [d{\bf{k}}](k^{2}_{\perp}+m^2) 
\lambda_{s} \left [  \frac{\vert N_{{\bf{k}}s}\vert ^2}{2\Omega_{{\bf{k}}s}}
\exp\left \{-\lambda_{s}\int \frac {e\dot{A}(t^{\prime})}
{\Omega_{{\bf{k}}s}(t^{\prime})}dt^{\prime}\right \} \right ]
\end{eqnarray}

Eqs.~(2.17) and (2.22) determine the normalization
constants, and we define
\begin{eqnarray}
\Gamma_{{s}}({\bf{k}})&\equiv &\frac{\vert N_{{\bf{k}}s}\vert ^2}
{\Omega_{{\bf{k}}s}}
 \exp\left \{ -\lambda_{s} \int \frac{e\dot{A}(t^{\prime})}
{\Omega_{{\bf{k}}s}(t^{\prime})} dt^{\prime}\right \} \nonumber \\
&=&\left [
\omega^2_{\bf{k}}+\left ( \frac{-\lambda_{s} e\dot{A}}{ 2\Omega_{{\bf{k}}s}}-
\frac{\dot{\Omega}_{{\bf{k}}s}}{2\Omega_{{\bf{k}}s}}\right ) ^2 
+\Omega^2_{{\bf{k}}s}
-\lambda_{s} 2p_3\Omega_{{\bf{k}}s}    \right ] ^{-1}
\end{eqnarray}
to represent the constants of the square brackets in (3.6).
With the identity
\begin{equation}
k^2_{\perp}+m^2=[\omega_{\bf{k}}+p_3][\omega_{\bf{k}}-p_3]  
\end{equation}
we obtain
\begin{eqnarray}
\ddot{A}&=&\langle 0\vert j^3\vert 0\rangle  \\ \nonumber 
&=&e\int {[d{\bf{k}}]}
\left \{  \frac{\omega_{\bf{k}}-p_3}
{\omega_{\bf{k}}+e\ddot{A}(\omega_{\bf{k}}+p_3)
/ 4\omega^3_{\bf{k}}} 
 -  \frac{\omega_{\bf{k}}+p_3}
{\omega_{\bf{k}}-e\ddot{A}(\omega_{\bf{k}}- p_3)
/ 4\omega^3_{\bf{k}}} 
 + O(\omega^{-3}_{\bf{k}})\right  \}  . 
\end{eqnarray}
 At large-momentum, we approximate
\begin{eqnarray}
\frac {1}{1\pm e\ddot{A}[\omega_{\bf{k}} \pm  p_3]
/4\omega^4_{\bf{k}}}
\simeq
1\mp \frac{e\ddot{A}}{4\omega^4_{\bf{k}}}[\omega_{\bf{k}}\pm p_3] .
\end{eqnarray}
 After we perform the angular integrations and drop
terms which are odd functions of $p_3$, the Maxwell
equation becomes
\begin{eqnarray}
\ddot{A}&=&\langle 0\vert j^3\vert 0\rangle  ={-e^2\ddot{A}} \int[d{\bf{k}}] 
\left ( \frac{1}{2\omega^3_{\bf{k}}} - \frac{p_3^2}{2\omega^5_{\bf{k}}}\right)
 + ({\it finite\; part})\\ \nonumber
&=& {-e^2\ddot{A}}\,\delta e^2 + ({\it finite\; part})
\end{eqnarray}
where 
\begin{equation}
\delta e^2 \equiv \int [d{\bf{k}}] 
\left ( \frac{1}{2\omega^3_{\bf{k}}} - \frac{p_3^2}{2\omega^5_{\bf{k}}}\right)
=\frac{1}{4\pi^2} \int_{0}^{\infty}dk\left[ 
\frac{k^2}{(k^2+m^2)^{\frac{3}{2}}}
- \frac{k^4}{3(k^2+m^2)^{\frac{5}{2}} } \right ].
\end{equation} 

The current in (3.11) diverges logarithmically, as expected, with the same 
divergence as the vacuum polarization $\Pi (q^2=0)$.
We define as usual
\begin{equation}
e^2_{R}=e^2(1+e^2\delta e^2)^{-1}\equiv Ze^2 ,\quad \quad A_R=Z^{-1/2}A \quad ,
\end{equation}
so that $eA=e_RA_R$. We can also write $Z=(1-e^2_R\delta e^2)$.
Multiplying (3.6) by $Ze/e_R$  we obtain
\begin{eqnarray}
\ddot{A}_R-e^2_R\ddot{A}_R\delta e^2=e_R\sum_{s=1}^{4}\int [d{\bf{k}}]
 (k^2_{\perp}+m^2)  
(-\lambda_{s})\Gamma_{s}({\bf{k}}) .
\end{eqnarray}
Upon subtracting the last term on the left-hand side,(3.14) becomes
\begin{equation}
\ddot{A}_R=e_R\int [d{\bf{k}}]\left [(k^2_{\perp}+m^2) 
 \sum_{s=1}^{4}(-\lambda_{s})\Gamma_{s}({\bf{k}})
+{e_R^2\ddot{A_R}} \left ( \frac{1}{2\omega^3_{\bf{k}}} - 
\frac{p_3^2}{2\omega^5_{\bf{k}}}\right) \right],
\end{equation}
so the expression for $\ddot{A}_R$ is finite (see  Eq.~(3.11)). 
Since this is so, and since $\Gamma_{{s}}$ depends on $e$ and $A$ 
through the product
$eA=e_R A_R$ only, the $R$ subscripts will be omitted from now on.

Consider now the difference between the exact expression (3.15) 
and its adiabatic 
approximation. Examining (3.9), (3.10), and (3.14), we can  write
\begin{eqnarray}
(k^2_{\perp}+m^2)  \sum_{s=1}^{4}(-\lambda_{s}) \Gamma_{s}({\bf{k}})  
=-\frac{2p_3}{\omega_{\bf{k}}}-\frac{e\ddot{A}}{2\omega^5_{\bf{k}}}
(\omega^2_{\bf{k}}-\pi^2_3) + R_{\bf{k}}(t) \quad .
\end{eqnarray}
$R_{\bf{k}}(t)$ is the difference between the exact current
integrand and its minimal adiabatic approximation. At large momentum
this minimal adiabatic approximation has to match the exact integrand 
up to terms
that fall off as $O(1/\omega^3_{\bf{k}})$, so the remainder 
$R_{\bf{k}}(t)$ falls off faster. 
Substituting (3.16) into (3.15) and using (3.9)--(3.11) 
the finite Maxwell equation reads
\begin{equation}
\ddot{A}=e\int [d{\bf{k}}]R_{\bf{k}}(t)  .
\end{equation}
Superficially it seems that the second derivative of $A$ appears only on the 
left-hand side of (3.17), but in fact the subsidiary condition (3.16)
defining $R_{\bf{k}}$ is an intrinsic part of (3.17).

We are interested in solving an initial value problem where the initial 
conditions for the Maxwell equation are given by

\begin{equation}
\dot{A}(t=0) = -E_0 ,\quad  A(t=0)=0   ,
\end{equation}
 and where the initial state is the adiabatic vacuum, selected by
 matching the exact solutions in (3.1) to their adiabatic approximation, 
viz.,
\begin{eqnarray}
\Omega _{{\bf{k}}s}(t=0)&=&\omega_{{\bf{k}}}(t=0) \\ \nonumber
\dot{\Omega} _{{\bf{k}}s}(t=0)&=&\dot{\omega}_{{\bf{k}}}(t=0)\ .
\end{eqnarray}

Non-vacuum initial conditions may be handled in a manner analogous to
the bosonic case, by adding to the current expectation value non-zero
particle number densities, without changing the initial conditions (3.19).
As in the bosonic case the initial conditions are not completely arbitrary,
because the asymptotic form of the adiabatic expansion (3.16) in which
$\ddot A$ and $R_{\bf{k}}$ appear must be consistent with the finite
renormalized Maxwell equation (3.17). By substituting (3.19) into (3.7)
and (3.16) we find that $\ddot{A}(0) = 0$, but

\begin{equation}
R_{\bf{k}}(0)=\frac{ 2p_3 }{\omega_{\bf{k}} }
\left \{ 1 - \frac{1+e^2E^2_0/4\omega^4_{\bf{k}}}
{ \left[ (1+e^2E^2_0/ 8\omega^4_{\bf{k}})^2- (e^2E^2_0p_3
/ 8\omega^5_{\bf{k}})^2 
\right ] } \right \}  
\end{equation}
is not zero. However, the integration over $\bf k$ in (3.17) is zero
by the charge conjugation symmetry $p_3 \rightarrow -p_3$.
Hence we find $\ddot{A}(0) = 0$ in the initial value of the Maxwell equation,
and the initial conditions (3.19) are consistent with the renormalization
requirements.
We emphasize that choosing $R_{\bf{k}}(0)=0$ is {\it not} consistent with
the initial conditions (3.19).

Given such a set of consistent initial conditions,
one can evolve the set of back-reaction equations (3.2), (3.16), and (3.17).
Even though the results we present below are limited to the simpler
(1+1)-dimensional case, we do have some experience with a full (3+1)-dimensional
calculation for scalar electrodynamics [24]. The following
technical comments should apply to the fermion problem as well.
By stepping  the back-reaction equations forward in time, one arrives at
$\{ A,\dot{A},\Omega_{{\bf{k}}s},\dot{\Omega}_{{\bf{k}}s} \}$
at each time.
In principle, $\ddot{A}$ can be calculated without using (3.17):
One can examine the asymptotic behavior of $\Gamma_s ({\bf{k}})$ at
large momentum and extract both $R_{\bf{k}}$ and $\ddot{A}$.
For instance, one can parametrize the left-hand side of (3.16) as
\begin{equation}
-\frac{2p_3}{\omega_{\bf{k}}}-\frac{e\ddot{A}}{2\omega^5_{\bf{k}}}
(\omega^2_{\bf{k}}-\pi^2_3) + \frac{\alpha}{\omega^4_{\bf{k}}}
+\frac{\beta}{\omega^5_{\bf{k}}}+ ...   \quad ,
\end{equation} 
and directly extract the value of $\ddot{A}$. The fitted value of 
$\ddot{A}$ depends, however, on the number of terms used to expand the
above series.
In addition, at some momentum a fit to a finite series
in powers of $\omega^{-1}_{\bf{k}}$ breaks down.
Instead, we propose an iterative scheme making use of (3.17).
We take $R_{\bf{k}}$ to be zero at an extremely large momentum
as a trial value, so $\ddot{A}$ is extracted from (3.16) automatically.
Using this value for $\ddot{A}$ we use (3.16) to
extract $R_{\bf{k}}$ for each ${\bf{k}}$ 
up to this very large momentum.
Then, substituting $R_{\bf{k}}$ in (3.17) we get a new and slightly
different value for $\ddot{A}$. This procedure may be iterated
until convergence of the sequence of $\ddot{A}$, $R_{\bf{k}}$ is reached.
Thereafter the next time step is taken.

\newpage
\section{Semiclassical QED in (1+1) dimensions}
For the  sake of simplicity we solve semiclassical QED in (1+1) 
dimensions for a system 
that initially is taken to be in an adiabatic vacuum state, with given
initial electric field and zero initial electric
current.
We follow the same steps as in the (3+1)-dimensional case.
Let us indicate the modifications  for (1+1) dimensions.
The $\gamma$ matrices are given by
\begin{equation}
\begin{array}{lr}
\gamma^{0}=\left(\begin{array}{cc}
            1 & 0\\
            0 & -1 \end{array} \right) ,           

            & \gamma^{1}=\left(\begin{array}{cc}
             0 & 1\\
            -1 & 0 \end{array} \right),

\end{array}
\end{equation} 
and $\gamma^1$ plays the role that $\gamma^3$ played
in the (3+1) dimensional case. 
For a spatially homogeneous electric field, we choose a
gauge  where $A^0=0$ and define $A^1=A(t)$, 
so that the second order Dirac equation becomes
\begin{equation} 
\left [ \Box +e^2 A^2 (t)+2iA(t)\partial_1-ie\partial_0 A(t)\gamma^0 \gamma^1
+m^2 \right ]\phi (x)=0.
\end{equation}
The Dirac equation in two dimensions has two independent solutions.
Either
\begin{equation}
 f_{k1}^+ \chi_1 \quad      
 f_{k1}^- \chi_1 \quad 
\end{equation}
or
\begin{equation}
 f_{k2}^+ \chi_2 \quad     
 f_{k2}^- \chi_2 
\end{equation}
may be taken as the basis set of independent solutions.
Here $\gamma^0 \gamma^1 \chi_{{s}}=\lambda_{s} \chi_{{s}}$ 
with $\lambda_1=1$ and $\lambda_2=-1$, and $\omega^2_k$ is now given by
\begin{equation}
\omega^2 _k (t) \equiv p_{1}^{2} +m^2  \quad ,\quad 
p_1 \equiv k^1-eA^1(t) \quad.
\end{equation} 
$\chi_{{s}}$ are given by
\begin{equation}
\begin{array}{lr}
\chi_1=\left(\begin{array}{cc}
            1 \\
            1 \end{array} \right),            

            & \chi_2=\left(\begin{array}{cc}
             1 \\
            -1 \end{array} \right).
\end{array}
\end{equation}

The current expectation value is
\begin{equation}
\ddot{A}=\langle 0\vert j^1\vert 0\rangle =2e\int \frac{dk}{2\pi}m^2
\{ \vert f_{k2}^{+}\vert ^2 - \vert f_{k1}^{+}\vert ^2   \}.
\end{equation}
Notice that in two dimensions there is no spin, and thus there are half as many
terms as in (2.30)
after summation.


Now that we have the matter field equation (2.12)
coupled to the Maxwell equation (4.7), we can solve the system numerically.
However, the parameter $e$ in these equations is not yet
renormalized . 
This renormalization can be done in the same way as  
in Section 3.  
In two dimensional QED the charge renormalization is finite and we find
that $\delta e^2=(6\pi m^2)^{-1}$. Therefore in the  renormalized
Maxwell equation, $\ddot{A}$ can be isolated [in contrast to (3.15)], 
and we obtain
\begin{eqnarray}
\ddot{A}_R&=&\frac{e_R}{1-e^2_R \delta e^2} \int  \frac{dk}{2\pi} m^2 
\left \{\sum_{s=1,2}\left(-\lambda_{s} \Gamma_{s}(k)\right ) 
+2\frac{p_1}{\omega_k} \right \}
\\ 
\Gamma_{{s}}(k)  &\equiv &
\left [\omega^2_{k}+\left ( \frac{-\lambda_{s} e_R\dot{A}_R}
{ 2\Omega_{ks}}-
\frac{\dot{\Omega}_{ks}}{2\Omega_{ks}}\right ) ^2 
+\Omega^2_{ks}
-\lambda_{s} 2p_1 \Omega_{ks}    \right] ^{-1}. 
\end{eqnarray}  
The last term in the braces in (4.8) does not contribute to the integral  
but was included for numerical purposes.
The set of equations (3.2) and (4.8)  with the
initial conditions  (3.18) and (3.20) defines the numerical back-reaction
problem.

We show in Figures~1~and~2  the time evolution of the scaled electric field
$\tilde{E}\equiv eE/m^2$ and the induced current $\tilde{j}\equiv ej/m^3$ 
as functions of $\tau\equiv mt$. Stability was attained for these results, 
as well
as for those presented below, with the time step $d\tau=5\times 10^{-4}$ and a
momentum grid with $d\tilde{k}\equiv dk/m =0.003$, forcing very long running 
times.  
With these strong initial electric fields we find that the induced current 
increases rapidly and becomes saturated at a constant value for some time,
then it decreases until it saturates at a higher constant value in the
opposite direction, after which clear plasma oscillations are seen.


The saturation in the first period is easy to understand. In a classical
kinetic picture, we have have $j=2en\langle v \rangle $, where $n$
is the density of particles (or antiparticles) and $v$ is their velocity.
$j$ saturates as $v$ is driven to the speed of light by the strong electric 
field; then $E$ changes its sign and the existing particles, together with
the newly
created particles, are accelerated in the
opposite direction until again they almost reach the speed of light. 
Due to the additional produced particles the absolute value of $j$ in
the second plateau is larger than in the first one. 

The envelope of the electric field amplitude decreases 
substantially only in the first few oscillations and remains almost constant	
at later times, i.e., the pair production happens essentially
in the first stages of the evolution.
In subsequent oscillations $n$ is larger and $E$ is weaker, 
and therefore the frequency of oscillations
increases. For quantitative comparisons to classical plasma oscillation
it should be noted that as long as the particles are relativistic
the frequency of plasma oscillations depends not only on $n$ but
also on the strength of $E$; the weaker the field the higher the frequency,
in contrast to the nonrelativistic case where the plasma
frequency does not depend on the amplitude of $E$.

As time  progresses, there develops a highly oscillatory behavior
(as a function of $k$) for the integrand in (4.8). 
However, we expect that the particle creation rate decreases for
weaker electric fields, and that the particle number becomes
approximately conserved at late times. Let us {\it define} the
particle number density in phase space by expanding the exact
fermionic field operator $\Psi$ in terms of the lowest, zeroth-order
mode functions, $f_k^{(0)}(t)$. Then the corresponding
Fock-space creation and destruction operators, $a_k^{(0)}$~and~$b_k^{(0)}$,
become time-dependent, and 

\begin{equation}
n(k;t)\equiv\langle 0|a^{(0)\dagger}_{k}(t) a^{(0)}_{k}(t) |0 \rangle
\end{equation}
may be computed by a Bogolyubov transformation (see Appendix A).
The result is shown in Figures~3(a) and 4(a). 
Clearly, the occupation number does not exceed one
in accordance with the Pauli principle.

In modelling the creation of a quark-gluon plasma in heavy-ion collisions, 
a flux-tube model has been used, employing a relativistic
Boltzmann-Vlasov equation with a particle source term.  
As in our calculation, one studies the evolution of a system with a homogeneous
(color) electric field, which initially
contains no particles [10--13] .
Let us compare the results of this phenomenological model
to those of our semiclassical analysis.
The relativistic kinetic equation in the presence of a homogeneous 
electric field is 
\begin{equation}
\frac{\partial f}{\partial t}+eE\frac{\partial f}{\partial p} =
\frac{dN}{dt\,dx\,dp}=-{\vert  eE(t)\vert  }
\log \left[1-\exp \left(-\frac {\pi m^2}{\vert  eE(t) \vert }\right)\right]
\, \delta (p),
\end{equation}
where
$f(p,t)$ is the ($x$-independent) classical phase space distribution,
expressed as a function of the {\it kinetic} momentum $p$, and
the right-hand side is the fermion pair-production 
rate in (1+1) dimensions.
We assume that the
particles are produced at rest [20], i.e., the source term is proportional to
$\delta (p)$. (The singularity on the right hand side of (4.11) for $ m=0 $
is spurious in that there can be no pair production from a classical homogeneous
electric field for massless fermions in one spatial dimension.  This
selection rule is not, however, built into the usual WKB arguments
[6] used to determine tunneling rates based on energy balance alone.)  
Initially, $f(p,0)=0$. Eq.~(4.11) may be solved using the characteristics
$\frac{dp}{dt}=eE$, giving \begin {eqnarray} f(p,t)=-\int_{0}^{t}dt'\,{\vert
eE(t')\vert } \,\log\left[1-\exp\left(-\frac{\pi m^2}{\vert eE(t')\vert
}\right)\right]\, \delta\bigl(p-eA(t')+eA(t)\bigr). \nonumber \\
\end{eqnarray}
The $\delta$ function allows us to perform the integration, and hence
\begin{eqnarray}
f(p,t)=-\sum_{i} \log \left[1-\exp\left(-\frac{\pi m^2}
{\vert eE(t_i)\vert }
\right)\right]\quad ,
\end{eqnarray}
where the $t_i$'s fulfill $p+eA(t)-eA(t_i)=0$ and $t_i<t$.  

The field equation for $A$ is
\begin{equation}
\frac{d^2 A}{dt^2}=j_{total}=j_{cond}+j_{pol}
\end{equation}
where the conduction current is
\begin{equation}
j_{cond}=2e\int\frac{dp}{2\pi}\,\frac{p}{\epsilon_p}\,f(p,t)  \quad,
\end{equation}
with $\epsilon_p\equiv\sqrt {p^2 + m^2}$, and the polarization current is
[12]
\begin{equation}
j_{pol}=\frac{2}{E}\int \frac{dp}{2\pi}\,\epsilon_p\frac{dN}{dt\,dx\,dp} \quad.
\end{equation} 
(The factors of 2 in (4.15) and (4.16) account for the contributions of 
the \break antiparticles.) 
Inserting (4.12) 
into (4.15) reduces the system to a single equation, 
\begin{eqnarray}
\frac{d^2 \tilde{A}}{d\tau ^2}&=&-\frac{e^2}{\pi m^2}\int_{0}^{\tau}\,d\tau'
\frac{\tilde{A}(\tau')-\tilde{A}(\tau)}{\sqrt{\left[\tilde{A}(\tau')-\tilde{A}
(\tau)\right]^2+1}}
\left\vert  \tilde{E}(\tau')\right\vert  
\log \left[1-\exp \left(-\frac {\pi}{\left\vert \tilde{E}(\tau')
\right\vert }\right)\right] \nonumber \\
&& -\frac{e^2}{\pi m^2}\,{\rm sign}\left(\tilde{E}(\tau)\right)
\log \left[1-\exp\left(-\frac {\pi}{\left\vert \tilde{E}(\tau)\right\vert }
\right)\right] 
\end{eqnarray}
in terms of the dimensionless variables $\tilde{A}\equiv eA/m$,
$\tilde{E}$, and $\tau$.


The time evolutions of $\tilde{E}$ and $\tilde\jmath$ are shown in the dashed
curves of Figures~1(a) and~2(a).
In the former, we see that for an initial field $\tilde{E}_{\tau=0}=1$ 
there is good quantitative
agreement between the results obtained with the two very different methods. 
The oscillations are {\it slower} and the electric fields decay
more {\it slowly} in the semiclassical calculation than in the
Boltzmann-Vlasov model.
For $\tilde{E}_{\tau=0}=4$ the plasma frequencies turn out to very
different in the two calculations.

The kinetic theory can be improved by use of a source
term that takes  Pauli blocking into account.
We replace the right-hand side of (4.11) with (see appendix B)
\begin{equation}
\frac{dN}{dt\,dx\,dp}=-\left ( 1-2f(p,t) \right ) {\vert eE(t)\vert }
\log \left[1-\exp \left(-\frac{\pi m^2}{\vert eE(t)\vert }\right)\right]\,
\delta(p)\quad.
\end{equation}
With this source term the agreement between the kinetic theory
and the quantum theory is even more striking, as demonstrated 
by the dashed-dotted lines in Figures~1(b) and~2(b).
The improvement is especially dramatic in the strong-field case,
Fig.~2.

The amplitude of the electric field approaches
a limiting value after a few oscillations, meaning that 
thereafter the production
of particles is negligible. In the boson case an analogous effect 
is seen when Bose enhancement [20] is considered, but it sets in somewhat later 
than for fermions.
The constant amplitude reflects the absence of pair creation from virtual
photons and the exponentially small spontaneous
pair creation at this stage of the evolution. The fact 
that the electric field reaches its
limiting value  more quickly for fermions than for bosons may be due to the
difficulty of producing more fermions once the low-momentum states have been
occupied.
Our semiclassical description is equivalent [19]
to QED with a large number $N$ of flavors. In the large $N$ limit any processes
involving virtual photons are excluded. To include pair production via a
virtual photon one must go to the next order in a $1/N$ expansion.
A systematic $1/N$ expansion for the pair
production problem can be obtained by using Schwinger's closed time loop
generating functional [25]. 

The distribution function $f(p,t)$, measured after the electric field has
reached an almost constant amplitude, may be compared to the quantum theory's
$n(k)$ after the latter is smoothed, as shown in Figures~3 and~4.
The curves have a relative
displacement due to the slightly different value of $A$.
This comparison displays the necessity of including the Pauli blocking
term in the phenomenological model. In the absence of this term the
occupation number exceeds one if the
initial electric field is strong enough as can be seen in Fig.~4(b).
This is to be expected since the source term in (4.10)
is very large for a strong electric field and it is the $(1-2f)$ term 
that suppresses the violation of the Pauli exclusion principle.

Because of the striking similarity between the results 
in Figures~3(c) and~4(c),
one can use the kinetic theory model
to explain detailed features of the particle distribution as has often been
assumed [9-13] in the past.
Past treatments have not, however, included the
Pauli blocking term which is crucial for strong fields.
A direct
derivation of the Boltzmann-Vlasov equation  through some chain of
approximations beginning with a field formulation would be of great interest in
elucidating the physical mechanisms involved here. A step in this direction
using the Wigner distribution function has been taken in Ref.~[26]. 

\newpage
\appendix
\section {Appendix A}


To obtain the number of particles per unit volume 
of phase space we consider 
\begin{eqnarray}
n({\bf{k}};t)=
\sum_{s=1,2} \langle 0\vert b^{(0)\dagger}_{s}({\bf{k}};t)b^{(0)}_{s}({\bf{k}};t) \vert 0 \rangle\ ,
\nonumber \\ \end{eqnarray}
where
we expand the field operator $\Psi (x)$ in terms of the adiabatic
order zero mode functions, viz.,
\begin{eqnarray}
\Psi (x)&=&\int [d{\bf{k}}]\sum_{r =1,2}\left[b^{(0)}_{r}({\bf{k}};t)
{u_{r,{\bf{k}}}}(t)e^{-i\int\omega_{\bf{k}} dt}
+{d^{(0)}_{r}}^{\dagger}({\bf{-k}};t){v_{r,{-\bf{k}}}}(t)
e^{i\int\omega_{\bf{k}}dt}\right] e^{i{\bf{k}}\cdot{\bf{x}}} ,\nonumber \\
\end{eqnarray}
having defined
\begin{equation}
\begin{array}{lr}

{u_{r,{\bf{k}}}}=\left(\begin{array}{cc}
               \sqrt{\frac{(\omega_{\bf{k}} +m)}{2\omega_{\bf{k}}}}
\varphi^{r}\\
\frac{{\vec{\sigma} \cdot \vec{\pi}}}
         {\sqrt{ 2\omega_{\bf{k}}
(\omega_{\bf{k}} +m)}}\varphi^{r}\end{array}\right),

            & {v_{r,{-\bf{k}}}}=\left(\begin{array}{cc}
            \frac{{-\vec{\sigma} \cdot \vec{\pi}}}
         {\sqrt{ 2\omega_{\bf{k}}(\omega_{\bf{k}}+m)}}\varphi^{r} \\
   \sqrt{\frac{(\omega_{\bf{k}}+m)}{2\omega_{\bf{k}}}}\varphi^{r}\end{array}
\right).            
\end{array}
\end{equation}
Here $\varphi^{r}$ are the spinors
\begin{equation}
\begin{array}{lr}
\varphi^1=\left(\begin{array}{cc}
            1 \\
            0 \end{array} \right),            

            & \varphi^2=\left(\begin{array}{cc}
             0 \\
             1 \end{array} \right).
\end{array}
\end{equation}
These ${u_{r,{\bf{k}}}(t)}$ and ${v_{r,{\bf{k}}}(t)}$  have been chosen to 
satisfy
\begin{equation}
{u_{r,{\bf{k}}}}^{\dagger}{u_{r',{\bf{k}}}}=\delta_{rr'} \quad
{v_{r,{\bf{k}}}}^{\dagger}{v_{r',{\bf{k}}}}=\delta_{rr'} \quad
{u_{r,{\bf{k}}}}^{\dagger}{v_{r',{\bf{k}}}}=0.
\end{equation}
They are {\it not} the $t\rightarrow \infty$ limits of $\psi_{s}^{\pm}$.
On the other hand, the field operator has the expansion (2.20), 
\begin{eqnarray}
\Psi (x)&=&\int [d{\bf{k}}]\sum_{s=1,2}[b_{s}({\bf{k}})
\psi^+_{\bf{k}s}
+d_{s}^{\dagger}({\bf{-k}})\psi^-_{\bf{k}s} ] \ ,
\end{eqnarray}
in terms of the time-independent operators $b_{s}({\bf{k}})$ and 
$d_{s}^{\dagger}({\bf{-k}})$.
The operators  $b^{(0)}_r({\bf{k}};t)$ and ${d^{(0)\dagger}_r ({\bf{k}};t)}$
are related to the operators
$b_s ({\bf{k}})$ and ${d_s^{\dagger} ({\bf{-k}})}$
by the following time-dependent Bogolyubov transformation:
\begin{eqnarray}
b^{(0)}_{r}({\bf{k}};t)&=&\sum_{s=1,2}
\alpha^{s}_{{\bf{k}}r}(t)b_{s}({\bf{k}})
+\beta^{s}_{{\bf{k}}r}(t){d_{s}}^{\dagger}({\bf{-k}}) \nonumber \\
{d^{(0)\dagger}_{r}}(-{\bf{k}};t)&=&\sum_{s=1,2}
\beta^{\ast s}_{{\bf{k}}r}(t)b_{s}({\bf{k}})+
\alpha^{\ast s}_{{\bf{k}} r}(t){d_{s}}^{\dagger}({\bf{-k}}).
\end{eqnarray}
Since we expand the field using two different bases, (A.2) and (A.6), 
which are distinct even without interaction, the Bogolyubov transformation
(A.7) has a matrix form.
From the canonical anticommutation relations it follows 
that 
\begin{equation}
\sum_{r=1,2}(\vert \alpha^{s}_{{\bf{k}} r}\vert ^2 + 
\vert \beta^{s}_{{\bf{k}} r}\vert ^2)=1,
\end{equation}
independent of time.

Substituting (A.7) into (A.2) we obtain
\begin{eqnarray}
\Psi (x)&=&\int [d{\bf{k}}]
\sum_{s =1,2}\sum_{r=1,2}\left[b_{s}({\bf{k}})
\left(\alpha^{s}_{{\bf{k}} r}{u_{r,{\bf{k}}}}e^{-i\int\omega_{\bf{k}} dt}
+\beta^{\ast s}_{{\bf{k}} r} {v_{r,{\bf{-k}}}}
e^{i\int\omega_{\bf{k}} dt} \right) \right. \nonumber \\
&& +{d_{s}^{\dagger}({\bf{-k}})}
\left.\left(\beta^{s}_{{\bf{k}} r}{u_{r,{\bf{k}}}}e^{-i\int\omega_{\bf{k}} dt}
+\alpha^{\ast s}_{{\bf{k}} r} {v_{r,{\bf{-k}}}}
e^{i\int\omega_{\bf{k}} dt}\right)\right] e^{i{\bf{k}}\cdot{\bf{x}}}.
\end{eqnarray}
We then identify $\psi^+_{{\bf{k}}s}$ and
$\psi^-_{{\bf{k}}s}$ as
\begin{eqnarray}
\psi^+_{{\bf{k}}s}&=&\sum_{r=1,2}
\left(\alpha^{s}_{{\bf{k}} r}{u_{r,{\bf{k}}}}
e^{-i\int\omega_{\bf{k}}dt}
+\beta^{\ast s}_{{\bf{k}} r} {v_{r,{-\bf{k}}}}
e^{i\int\omega_{\bf{k}} dt} \right) e^{i{\bf{k}}\cdot{\bf{x}}} \nonumber \\
\psi^-_{{\bf{k}}s}&=&\sum_{r=1,2}
\left(\beta^{s}_{{\bf{k}} r}
{u_{r,{\bf{k}}}}e^{-i\int\omega_{\bf{k}} dt}
+\alpha^{\ast s}_{{\bf{k}} r} {v_{r,{-\bf{k}}}}
e^{i\int\omega_{\bf{k}} dt} \right)e^{i{\bf{k}}\cdot{\bf{x}}}
\end{eqnarray}

The number of particles (or antiparticles) produced per unit phase space volume
at a given momentum is then given by
\begin{equation} 
n({\bf{k}};t)=\sum_{r=1,2}\langle 0_{in}\vert {b^{(0){\dagger}}_r({\bf{k}};t)}
b^{(0)}_r ({\bf{k}};t)\vert 0_{in}\rangle =\sum_{s=1,2}\sum_{r=1,2}
\vert \beta^{s}_{{\bf{k}} r}(t)\vert ^2.
\end{equation}
We now calculate $\vert \beta^{s}_{{\bf{k}} r}\vert ^2$ in terms of the 
solutions of the quadratic Dirac equation $f_{{\bf{k}}s}^+$.
Multiplying the first equation in (A.10) by ${v_{r,{\bf{k}}}}^{\dagger}$
we find
\begin{equation}
\sum_{r=1,2}\vert {v_{r,{\bf{k}}}}^{\dagger}
\psi^+_{{\bf{k}}s}\vert ^2
=\sum_{r=1,2}\vert \beta^{s}_{{\bf{k}} r}\vert ^2.
\end{equation}
Using (2.12) and (2.14), $\psi^+_{{\bf{k}}s}$ can
be written as
\begin{equation}
\psi^+_{{\bf{k}}s}\equiv \hat{D}f_{{\bf{k}}s}^+\chi_{s}=
 (i\gamma^{0}\partial_{0}-\gamma^{i}p_{i}+m)f_{{\bf{k}}s}^+\chi_{s}.
\end{equation}
We may use the Dirac representation of the $\gamma$ matrices to
recast $\psi^+_{{\bf{k}}s}$ explicitly, 
\begin{equation}
\begin{array}{lr}
\hat{D}=\left(\begin{array}{cccc}
  (i\partial_0+m)  & 0 & -p_3   & (-p_1+ip_2) \\
  0 & (i\partial_0+m) & (-p_1-ip_2) & p_3  \\
  p_3  & (p_1-ip_2) & (-i\partial_0+m) & 0 \\
  (p_1+ip_2)  & -p_3  & 0 & (-i\partial_0+m)   
\end{array} \right) ,           
\end{array}
\end{equation}
and it follows that
\begin{equation}
\begin{array}{lr} 

\psi^+_{{\bf{k}}1}\equiv \hat{D}f_{{\bf{k}}1}^+\chi_{1}=
\left(\begin{array}{c}
            i\dot{f}_{{\bf{k}}1}^+ +(m-p_3)f_{{\bf{k}}1}^+\\
            -(p_1+ip_2) f_{{\bf{k}}1}^+\\
            - i\dot{f}_{{\bf{k}}1}^+ +(m+p_3)f_{{\bf{k}}1}^+\\
             (p_1+ip_2) f_{{\bf{k}}1}^+\end{array} \right), 
\end{array}
\end{equation}
\begin{equation}
\begin{array}{lr} 

 \psi^+_{{\bf{k}}2}\equiv \hat{D}f_{{\bf{k}}2}^+\chi_{2}
=\left(\begin{array}{c}
             (p_1-ip_2) f_{{\bf{k}}2}^+ \\
             i\dot{f}_{{\bf{k}}2}^+ +(m-p_3)f_{{\bf{k}}2}^+ \\
             (p_1-ip_2) f_{{\bf{k}}2}^+ \\
              i\dot{f}_{{\bf{k}}2}^+ +(-m-p_3)f_{{\bf{k}}2}^+ 
\end{array} \right). 
\end{array}
\end{equation}
Let us also write ${u_{r,{\bf{k}}}}$ and
 ${v_{r,{-\bf{k}}}}^{\dagger}$ in their explicit form
\begin{eqnarray}
{u_{1,{\bf{k}}}}^{\dagger}&=&(\omega_{\bf{k}} +m,\quad 0,\quad p_3,
\quad p_1 -ip_2)/
\sqrt{2\omega_{\bf{k}}(\omega_{\bf{k}} + m)}   \nonumber \\
{u_{2,{\bf{k}}}}^{\dagger}&=&( 0,\quad \omega_{\bf{k}} +m,\quad p_1+ip_2,
\quad -p_3)/
\sqrt{2\omega_{\bf{k}}(\omega_{\bf{k}} + m)} \\ 
{v_{1,{-\bf{k}}}}^{\dagger}&=&(-p_3,\quad -p_1 +ip_2,\quad\omega_{\bf{k}} 
+m
,\quad 0)/
\sqrt{2\omega_{\bf{k}}(\omega_{\bf{k}} + m)}  \nonumber \\
{v_{2,{-\bf{k}}}}^{\dagger}&=&(-p_1 -ip_2,\quad +p_3,\quad 0,
\quad \omega_{\bf{k}} +m)/
\sqrt{2\omega_{\bf{k}}(\omega_{\bf{k}} + m)} \nonumber
\end{eqnarray}
Thus, we obtain
\begin{eqnarray}
\vert \beta^{\ast 1}_{\bf{k}1} \vert
&=&\vert {v_{1,{-\bf{k}}}}^{\dagger}\psi^+_{{\bf{k}}1} \vert
= \left \vert \frac{(\omega_{\bf{k}} +m+p_3)\left [\omega_{\bf{k}} 
f_{{\bf{k}}1}^+ -i\dot{f}_{{\bf{k}}1}^+\right ]}
{\sqrt{2\omega_{\bf{k}}(\omega_{\bf{k}} + m)}}  \right \vert  \nonumber \\
\vert \beta^{\ast 1}_{\bf{k}2} \vert
&=&\vert {v_{2,{-\bf{k}}}}^{\dagger}\psi^+_{{\bf{k}}1} \vert =
\left  \vert \frac{(p_1+ip_2)\left [\omega_{\bf{k}} f_{{\bf{k}}1}^+ 
-i\dot{f}_{{\bf{k}}1}^+ \right ]}
{\sqrt{2\omega_{\bf{k}}(\omega_{\bf{k}} + m)}} \right  \vert \\ 
\vert \alpha^{1}_{\bf{k}1} \vert 
&=&\vert {u_{1,{\bf{k}}}}^{\dagger}\psi^+_{{\bf{k}}1} \vert =
\left  \vert \frac{(\omega_{\bf{k}} +m-p_3)
\left [\omega_{\bf{k}} f_{{\bf{k}}1}^+ +i\dot{f}_{{\bf{k}}1}^+ \right ]}
{\sqrt{2\omega_{\bf{k}}(\omega_{\bf{k}} + m)}}\right  \vert  \nonumber \\
\vert \alpha^{1}_{\bf{k}2} \vert
&=&\vert {u_{2,{\bf{k}}}}^{\dagger}\psi^+_{{\bf{k}}1} \vert =
\left  \vert \frac{-(p_1+ip_2)
\left [\omega_{\bf{k}} f_{{\bf{k}}1}^+ +i\dot{f}_{{\bf{k}}1}^+ \right ]}
{\sqrt{2\omega_{\bf{k}}(\omega_{\bf{k}} + m)}}\right   \vert \nonumber
\end{eqnarray}
To see the difference between the two different bases,
we observe that in the free case $f^{+}_{s}=c_{s}e^{-i\omega_{\bf{k}} t}$,
(A.15) and (A.16) are different from $u_{1,{\bf{k}}}$ and $u_{2,{\bf{k}}}$.

Using (A.18), (2.17), and (2.22) it is straightforward to show that equation
(A.8) is satisfied.
Finally we find that the number of particles produced 
is given by
\begin{eqnarray}
n({\bf{k}};t)=\sum_{\stackrel{s=1,2}{or\quad s=3,4}}\frac{(\omega_{\bf{k}}+\lambda
_{s} p_3)}
{\omega_{\bf{k}}} 
\left[\omega^{2}_{\bf{k}} \vert f_{{\bf{k}}s}^+\vert ^2
+ \vert \dot{f}_{{\bf{k}}s}^+ \vert ^2-i\omega_{\bf{k}}
({f}_{{\bf{k}}s}^{\ast +} \dot{f}_{{\bf{k}}s}^+
- \dot{f}_{{\bf{k}}s}^{\ast +} {f}_{{\bf{k}}s}^+) \right ]  \nonumber \\
=\sum_{\stackrel{s=1,2}{or\quad s=3,4}}
\frac{(\omega_{\bf{k}}+\lambda_{s} p_3)}{2\omega_{\bf{k}}}  
\frac{\left [ \omega^2_{\bf{k}}+\left ( \frac{-\lambda_{s} e\dot{A}}
{ 2\Omega_{{\bf{k}}s}}-
\frac{\dot{\Omega}_{{\bf{k}}s}}{2\Omega_{{\bf{k}}s}}\right ) ^2 
+\Omega^2_{{\bf{k}}s}
- 2\omega_{\bf{k}}\Omega_{{\bf{k}}s}    \right]}
{\left [
\omega^2_{\bf{k}}+\left ( \frac{-\lambda_{s} e\dot{A}}{ 2\Omega_{{\bf{k}}s}}-
\frac{\dot{\Omega}_{{\bf{k}}s}}{2\Omega_{{\bf{k}}s}}\right ) ^2 
+\Omega^2_{{\bf{k}}s}
-\lambda_{s} 2p_3\Omega_{{\bf{k}}s}    \right]}  \nonumber \\
\end{eqnarray}
Averaging these two summations, we find
\begin{eqnarray}
n({\bf{k}};t)=\sum_{s=1}^{4}\frac{(\omega_{\bf{k}} +\lambda_{s} p_3)}
{2\omega_{\bf{k}}} 
\left[\omega^{2}_{\bf{k}} \vert f_{{\bf{k}}s}^+\vert ^2
+ \vert \dot{f}_{{\bf{k}}s}^+ \vert ^2-i\omega_{\bf{k}}
({f}_{{\bf{k}}s}^{\ast +} \dot{f}_{{\bf{k}}s}^+
- \dot{f}_{{\bf{k}}s}^{\ast +} {f}_{{\bf{k}}s}^+) \right ]  \nonumber \\
=\sum_{s=1}^{4}
\frac{(\omega_{\bf{k}} +\lambda_{s} p_3)}{4\omega_{\bf{k}}}  
\frac{\left [ \omega^2_{\bf{k}}+\left ( \frac{-\lambda_{s} e\dot{A}}
{ 2\Omega_{{\bf{k}}s}}-
\frac{\dot{\Omega}_{{\bf{k}}s}}{2\Omega_{{\bf{k}}s}}\right ) ^2 
+\Omega^2_{{\bf{k}}s}
- 2\omega_{\bf{k}}\Omega_{{\bf{k}}s}    \right]}
{\left [
\omega^2_{\bf{k}}+\left ( \frac{-\lambda_{s} e\dot{A}}{ 2\Omega_{{\bf{k}}s}}-
\frac{\dot{\Omega}_{{\bf{k}}s}}{2\Omega_{{\bf{k}}s}}\right ) ^2 
+\Omega^2_{{\bf{k}}s}
-\lambda_{s} 2p_3\Omega_{{\bf{k}}s}    \right]} .
\end{eqnarray}
In (1+1) dimensions the number of produced particles per unit length at momentum
 ${\bf{k}}$
is obtained  by replacing $p_3$ by $p_1$ in Eq.~(A.20) and dividing
 by 2. This is the quantity which is computed from our numerical solution of the equations and plotted in Figures 3(a) and 4(a).


\newpage
\section {Appendix B}
   

We derive here the Pauli blocking correction factor 
of Eq.~(4.18).
The operators  $b_r({\bf{k}};t+\Delta t)$ and 
${d^{\dagger}_r (-{\bf{k}};t+\Delta t)}$
 are related to the operators
$b_r ({\bf{k}};t)$ and ${d_r^{\dagger} ({\bf{-k}};t)}$
by the following Bogolyubov transformation:
\begin{eqnarray}
b_{r}({\bf{k}};t+\Delta t)&=&
\alpha_{{\bf{k}}r}(t+\Delta t)\, b_{r}({\bf{k}};t)
+\beta_{{\bf{k}}r}(t+\Delta t)\, d^{\dagger}_{r}({\bf{-k}};t)\\ \nonumber
{d^{\dagger}_{-r}}(-{\bf{k}};t+\Delta t)&=&
\beta^{\ast}_{{\bf{k}}r}(t+\Delta t)\, b_{r}({\bf{k}};t)+
\alpha^{\ast}_{{\bf{k}} r}(t+\Delta t)\, d_{r}^{\dagger}({\bf{-k}};t).
\end{eqnarray}
From the time-independence of the anticommutation relations it follows 
that 
\begin{equation}
\left \vert \alpha_{{\bf{k}} r}(t+\Delta t)\right \vert ^2 + 
\left \vert \beta_{{\bf{k}} r}(t+\Delta t)\right \vert ^2=1.
\end{equation}

Let us denote the density of particles (antiparticles) in momentum ${\bf{k}}$ 
and spin state $r$ at time $t$ by $n^{\pm}_{{\bf{k}}r}(t)$. We emphasize
that the wavenumber $\bf k$ refers to the {\it canonical} momentum
of the produced particle(s), which remains constant in time. We shall
denote by $\bf p$ the corresponding time-dependent {\it kinetic} momentum
[see (2.13)].
Using the Bogolyubov transformation (B.1) the number of particles at time 
$t+\Delta t$ is given by
\begin{eqnarray} 
&&n({\bf{k}};t+\Delta t) \nonumber \\
&=&\sum_{r=1,2}
\langle n^{+}_{{\bf{k}}1}(t)n^{+}_{{\bf{k}}2}(t) n^{-}_{{\bf{k}}1}(t)
n^{-}_{{\bf{k}}2}(t)\vert 
{b^{{\dagger}}_r({\bf{k}};t+\Delta t)} b_r ({\bf{k}};t+\Delta t)
\vert n^{+}_{{\bf{k}}1}(t)n^{+}_{{\bf{k}}2}(t) n^{-}_{{\bf{k}}1}(t)
n^{-}_{{\bf{k}}2}(t)\rangle \nonumber \\
&=&\sum_{r=1,2} \left \{
\vert \alpha_{{\bf{k}} r}(t+\Delta t)\vert ^2 n^{+}_{{\bf{k}}r}(t)+
\vert \beta_{{\bf{k}} r}(t+\Delta t)\vert ^2 
\left (1-n^{-}_{{\bf{k}}r}(t)\right )\right \} \\ 
&=& \sum_{r=1,2}\left [n^{+}_{{\bf{k}}r}(t)+
\vert \beta_{{\bf{k}} r}(t+\Delta t)\vert ^2\left 
(1-n^{+}_{{\bf{k}}r}(t)-n^{-}_{{\bf{k}}r}(t) \right )\right ] \nonumber .
\end{eqnarray} 
Due to the homogeneity of the problem, at each moment the number of 
particles is equal to the number of  antiparticles  with opposite quantum
number $r$, i.e.,
$n^{+}_{{\bf{k}}r}(t)=n^{-}_{{\bf{k}}r}(t)$. The probability to create 
particles in each of the two spin states is equal, so that 
$\vert \beta_{{\bf{k}} 1}(t)\vert ^2=\vert \beta_{{\bf{k}} 2}(t)\vert ^2$. 
Introducing $n({\bf{k}};t)=n^{+}_{{\bf{k}}1}(t)+n^{+}_{{\bf{k}}2}(t)$ 
in (B.3) we obtain
\begin{equation}
n({\bf{k}};t+\Delta t)=n({\bf{k}};t)+2\vert \beta_{{\bf{k}}}(t+\Delta t)\vert ^2
[1-n({\bf{k}};t)].
\end{equation}
Finally we write the rate of pair production in (3+1) dimensions as
\begin{equation}
\frac{\Delta n({\bf{k}};t)}{\Delta t}=2R_{\bf{k}}[1-n({\bf{k}};t)]\ ,
\end{equation}
where $R_{\bf{k}}={\vert \beta_{\bf{k}}(t+\Delta t)\vert ^2}/{\Delta t}$
is the rate of 
pair production  when $n({\bf{k}},t) = 0$.
Schwinger's result for pair-production assumed that no particles are
initially present. We will assume as a reasonable ansatz that even in the
presence of particles $R_{\bf{k}}$  is given by Schwinger's form for pair production from
the vacuum.  Thus we assume
\begin{equation}
R_{\bf{k}}=-{\vert  eE(t)\vert  }
\log \left[1-\exp \left(-\frac {\pi(k^2_{\perp}+ m^2)}
{\vert  eE(t) \vert }\right)\right]
\, \delta (k_{\parallel}-eA(t)),
\end{equation}
where the particles are assumed to be created at rest, {\it i. e.},
with kinetic momentum $p_{\parallel} = k_{\parallel}-eA(t) = 0$.

The effect of Pauli blocking is incorporated by the factor multiplying this
form of $R_{\bf{k}}$ in (B.5). 
In (1+1) dimensions there is no spin, and therefore
the analogous result has the form
\begin{equation}
n(k;t+\Delta t)=n(k;t)+\vert \beta_{k}(t+\Delta t)\vert ^2[1-2n(k;t)],
\end{equation}
and we can write the rate of pair production as
\begin{equation}
\frac{\Delta n(k;t)}{\Delta t}=R_k[1-2n(k;t)],
\end{equation}
where $R_k={\vert \beta_{k}(t+\Delta t)\vert ^2}/{\Delta t}$ .
Now we can identify the canonical phase space density $n(k;t)$ with the
Boltzmann density function in terms of the kinetic momentum $p$, viz.,
\begin{equation}
n(k;t) =  f(p,t) = f\bigl(k-eA(t),t\bigr)\ .
\end{equation}

We conclude that in the phenomenological Boltzmann-Vlasov approach,
Pauli suppression is properly taken into account by  
the source term given in Eq.~(B.8). This is precisely what we have used in
(4.18).

We note that in [20] the normalization of the source term was 
smaller by $2\pi$  than in (4.11) and at the same time
the definition of the current was larger by $2\pi$ than
in Eqs.~(4.15) and (4.16). When one
takes the Pauli-blocking term into account
in Eq.~(4.18) (or the boson enhancement factor in Eq.~(17) of Ref.~[20]), 
one must use the normalization convention of this work. 
The numerical results presented in [20] reflect correct normalization.
\newpage

\centerline{\bf ACKNOWLEDGMENTS}
\vskip\baselineskip
We wish to thank Dr. I.~Paziashvili for drawing our attention
to Refs.~[14]--[16], and for very helpful discussions concerning
them and this work in general.
This work was partially supported by the German-Israel Foundation.
Further support was provided by
the Ne'eman Chair in Theoretical Nuclear Physics at Tel Aviv University.
The work of B. S. was supported
by a Wolfson Research Award administered by the Israel Academy of
Sciences and Humanities. 
Y.K. thanks the Theoretical Division of Los Alamos National Laboratory
for their hospitality. 
J.M.E. thanks the Institute for Theoretical Physics
of the University of Frankfurt, and its director, Professor Walter
Greiner, for their hospitality, and
the Alexander von Humboldt-Stiftung for partial support of this work.


\newpage
\centerline{\bf REFERENCES} {\frenchspacing
\begin{enumerate}
\item F.~Sauter, Z. Phys. {\bf 69}, 742 (1931).
\item J.~Schwinger, Phys. Rev. Lett. {\bf 82}, 664 (1951).
\item N.~B.~Narozhny and A.~I.~Nikishov, Yad. Fiz. {\bf11}, 1072 (1970)
[Sov. J. Nucl. Phys. {\bf11}, 596 (1970)].
\item M.~S.~Marinov and V.~S.~Popov, Fortsch. Phys. {\bf25}, 373 (1977).
\item E.~Brezin and C.~Itzykson, Phys. Rev. D {\bf2}, 1191 (1970).
\item A.~Casher, H.~Neuberger, and S.~Nussinov, Phys. Rev. D 
	{\bf28}, 179 (1979).
\item B.~Andersson, G.~Gustafson, G.~Ingelman, and T.~Sjostrand,
       Phys. Rep. {\bf 97} (1983) 31.
\item T.~S.~Bir\'{o}, H.~B.~Nielsen, and J.~Knoll, Nucl. Phys. {\bf B245},
	449 (1984).

\item A.~K.~Kerman, T.~Matsui, and B.~Svetitsky, Phys. Rev. Lett. {\bf 56},
	219 (1986).
\item K.~Kajantie and T.~Matsui, Phys. Lett {\bf164B}, 373 (1985).
\item A.~Bia\l as and W.~Czy\.z, Acta Phys. Pol. B {\bf17}, 635 (1986).
\item G.~Gatoff, A.~K.~Kerman, and T.~Matsui, Phys. Rev. D {\bf36},
	114 (1987).
\item  A.~Bia\l as, W.~Czy\.z, A.~Dyrek, and W.~Florkowski, Nucl. Phys.
{\bf B296}, 611 (1988).
\item A.~A.~Grib, S.~G.~Mamayev and V.~M.~Mostepanenko, 
{\it Quantum Effects in Strong External Fields}
(Atomizdat, Moscow, 1980). (In Russian.)
\item  V.~M.~Mostepanenko, Yad. Fiz. {\bf 30}, 208 (1979)
[Sov. J. Nucl. Phys. {\bf 30}, 107 (1979)]. 
\item S.~G.~Mamaev and N.N.~Trunov, Yad. Fiz. {\bf 30}, 1301 (1979)
[Sov. J. Nucl. Phys. {\bf 30}, 677 (1979)]. 
\item  W.~M.~Suen and P.~R.~Anderson, Phys. Rev. D {\bf 35}, 2940 (1987). 
\item  W.~M.~Suen, Phys. Rev. D {\bf 40}, 315 (1989).
\item F.~Cooper and E.~Mottola, Phys. Rev. D {\bf40}, 456 (1989).
\item Y.~Kluger, J.~M.~Eisenberg, B.~Svetitsky, F.~Cooper, and E.~Mottola,
Phys. Rev. Lett. {\bf 67}, 2427 (1991). 
\item J.~D.~Bjorken and S.~D.~Drell, {\it Relativistic Quantum Mechanics}
\break(McGraw-Hill, New York, 1964).
\item A.~A.~Grib, V.~M.~Mostepanenko, and V.~M.~Frolov, Teor. Mat. Fiz.
{\bf 13}, 377 (1972) 
[Theor. Math. Phys. {\bf 13}, 1207 (1972)]. 
\item P.~C.~Waterman, Am. J. Phys. {\bf 41}, 373 (1973).
\item Y.~Kluger, J.~M.~Eisenberg, B.~Svetitsky, F.~Cooper, and E.~Mottola,
in preparation.
\item J.~Schwinger, J. Math. Phys. {\bf 2}, 407 (1961).
\item I.~Bialynicki-Birula, P.~Gornicki, and J.~Rafelski, Phys. Rev. D {\bf 44}, 1825 (1991).

\end{enumerate}}
\newpage
\centerline{\bf FIGURE CAPTIONS}
\vskip\baselineskip
\parindent0pt
\parskip14pt
FIGURE 1.
(a) Time evolution of the scaled electric field $\tilde E$ and current
$\tilde \jmath$,
with initial value $\tilde E=1.0$ and coupling $e^2/m^2=0.1$.
Solid line is semiclassical QED, and dashed line is
Boltzmann-Vlasov model.
(b) Same as (a), but the dashed-dotted curve is
the Boltzmann-Vlasov model with the Pauli blocking correction,
Eq.~(4.18).

FIGURE 2.
As in Fig.~1, but for initial field $\tilde E=4.0$.

FIGURE 3. (a) Momentum distribution of produced pairs, for the evolution shown
in Fig.~1, at time $\tau=600.$  The abscissa is the scaled kinetic momentum
$\tilde p\equiv\tilde k-\tilde A$, with $\tilde k\equiv k/m$.
(b) Data of (a) after binning (histogram), compared with 
Boltzmann-Vlasov model (curve).
(c) Data of (a) after binning (histogram), compared with 
Boltzmann-Vlasov model with Pauli blocking 
(curve).

FIGURE 4. The same as Fig.~3, but for the evolution shown
in Fig.~2 (i.e., initial field $\tilde E=4.0$), at time $\tau=200.$  

\end{document}
As this paper was being written, the work of Grib and coworkers
[18--20] came to our attention.  Because of differences between our
renormalization method and theirs, as well as for the sake of
completeness, we present our version in some detail.  We go beyond
Refs. [18--20] in presenting results and a
consistent physical interpretation thereof; moreover, the comparison
with the Boltzmann-Vlasov model is novel and should be of interest to
phenomenologists.





