I think there is no problem (at least, no mathematical problem) in defining time evolution of the Dirac equation by way of a Hamiltonian.
The Dirac equation can be written in the form
<br />
\partial_0 \psi(t,x) = -i (\gamma^0 m + \gamma^0\gamma^j\partial_j) \psi(t,x)<br />
where psi is a 4-component spinor, and sum over j=1,2,3.
If we take Fourier transforms in the 3 spatial dimensions, that is, define
<br />
\hat{\psi}_t(k) = \int d^3x \; \psi(t,x) \; e^{-ikx} \qquad<br />
\psi(t,x) = \int \frac{d^3k}{(2\pi)^{3}} \; \hat{\psi}_t(k) \; e^{ikx}<br />
then the Dirac equation becomes
<br />
\partial_0 \hat{\psi}_t(k) &= -i (\gamma^0 m + \gamma^0\gamma^j ik_j) \hat{\psi}_t(k).<br />
Put
<br />
H = (\gamma^0 m + \gamma^0\gamma^j ik_j), <br />
and this can be solved; the solution is
<br />
\hat{\psi}_t(k) &= \exp( -iHt) \hat{\psi}_0(k).<br />
where \psi_0 are the initial conditions. It's possible to evaluate the exponential in closed form, because H^2 simplifies using the result H^2 = E^2 = (k^2 + m^2).
<br />
\exp( -iHt ) = \left( \cos( Et ) - \frac{i \sin( Et )}{E}H \right)<br />
Given that \hat{\psi}_t is now given explicitly for all t, this implies (by inverse Fourier transform) that \psi(t,x) is constrained for all t.
As a related comment, I thought it was quite interesting to note that evolution in Fourier space is simply pointwise multiplication by the matrix exp( -iHt ). This means that evolution in real space is convolution with the Fourier transform of exp( -iHt ), which is precisely what we mean by a propagator. I *think* that if you explicitly evaluate the Fourier transform of exp( -iHt ), you get to the retarded Green's function, but I got stuck with the math.
Dave