So I start with the Hamiltonian operator for a Dirac field:
$$\hat{H} = \gamma.\partial + i \gamma_5 m$$
I define the transition Green's function so that it acts like the time evolution operator:
$$\int G(x-y,t_2-t_1)\psi(y)dy^3 = e^{i(t_2-t_1)\hat{H}}\psi(x)$$
Then converting to momentum coordinates we get:
$$G(x-y,t_2-t_1)= \int e^{i(t_2-t_1)(\gamma.k + \gamma_5 m) + ik.(x-y)} dk^3 $$
So this seems to satisfy the Dirac equation: $(\partial_t + \gamma.\partial +i \gamma_5 m)G(x-y,t_2-t_1)=0$. But for the Feynman propagator the RHS should be a delta function not zero. So I must be doing something wrong. Also it doesn't look much like the Feynamn propagator for a fermion.
(I may have missed out some factor's of $i$ some places).
Edit: This method seems to work for the non-relativistic case with:
$$\hat{H} = \frac{1}{2m}\partial^2$$
in which case I get:
$$G(x-y,t_2-t_1) = \int e^{i(t_2-t_1)(|k|^2/2m)+ik.(x-y)}dk^3 = e^{-im|x-y|^2/(t_2-t_1)} \left(\frac{im}{t_2-t_1}\right)^{3/2}$$
So I must just be going wrong for the Dirac case. Can you see where?
Edit 2: I think that maybe it's a case of taking the right limits to get the delta functions? e.g. when t is small in $e^{-x^2/t}/\sqrt{t}$ it becomes a delta function.