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
[tex]
\partial_0 \psi(t,x) = -i (\gamma^0 m + \gamma^0\gamma^j\partial_j) \psi(t,x)[/tex]
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
[tex]
\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}[/tex]
then the Dirac equation becomes
[tex]
\partial_0 \hat{\psi}_t(k) &= -i (\gamma^0 m + \gamma^0\gamma^j ik_j) \hat{\psi}_t(k).[/tex]
Put
[tex]
H = (\gamma^0 m + \gamma^0\gamma^j ik_j), [/tex]
and this can be solved; the solution is
[tex]
\hat{\psi}_t(k) &= \exp( -iHt) \hat{\psi}_0(k).[/tex]
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).
[tex]
\exp( -iHt ) = \left( \cos( Et ) - \frac{i \sin( Et )}{E}H \right)[/tex]
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