quantum field theory - Propagator for Dirac equation in real space

I'm interested in the retarded propagator for a free massless Dirac fermion, i.e. solutions $ψ$ to the inhomogeneous PDE

$$(∂_t- \nabla·\vec σ) ψ(x,t) = f(x,t)$$

with boundary conditions

$$\quad ψ(x,t) \to 0 \text{ for } t \to -∞$$

where $\vec σ = (σ_1,σ_2,σ_3)^T$ are the three Pauli matrices. (The boundary conditions can be even more restrictive, I just want the solution to decay sufficiently quickly at infinity so that it becomes unique and has a well-defined Fourier transform.)

Now, solving the Dirac equation is a standard exercise in virtually every QFT book, but all the books I've looked at only consider the Fourier transform of the propagator.

However, I am interested in the real space formula for the retarded propagator

Using the retarded propagator for the wave equation in $3+1$ dimensions, we can write

$$ψ(x,t) = (∂_t + \nabla·\vec σ)(∂_t^2 - \nabla^2)^{-1} f(x,t)$$

$$= (∂_t + \nabla·\vec σ) \frac1{4π·\text{something}}∫d^3x'dt' \frac1{|x-x'|}\delta(|x-x'|-|t-t'|) f(x',t')$$

but this formula strikes me as seriously weird: carrying out the differentiation with respect to $x$ and $t$ will differentiate the $\delta$-function in the integral, which means that the solution depends on the derivatives of the function $f$. This goes against my intuition that a linear first-order PDE should depend on the initial values directly, and not on their time and space derivatives!

Is there a reference where I can find a discussion of the retarded propagator of the (massless) Dirac equation in real space?

Answer

The real space propagator for the massive Dirac fermion in $3+1$ dimensions is calculated in R. Feynman's book Quantum Electrodynamics (Lecture 17, page 84 in the edition linked to).

The result is very much as indicated in the question: first solve the wave equation, then differentiate the solution with the Dirac operator again. In particular, Feynman calculates the propagator of the Klein-Gordon equation in real space:

$$I_+(t,x) = ∫ \frac{d^4p}{(2π)^4} \frac{\exp[-i(p\cdot x)]}{p^2 - m^2 + i\varepsilon} = -(4π)^{-1} \delta(s^2) + \frac{m}{8πs}H_1^{(2)}(ms)$$

Here, $s = +(t^2-x^2)^{1/2}$ for $t>|x|$ and $s = -i(x^2-t^2)^{1/2}$ for $t < |x|$. Moreover, $\delta(s^2)$ is a delta function and $H^{(2)}_1(ms)$ is a Hankel function. Then, you have to differentiate the delta function, indeed.

However, note that to be physically meaningful, the propagator $G(x_2,t_2;x_1,t_1)$ for the Dirac equation should only take into account the positive energy states for the retarded time frame $t_2 - t_1 > 0$, while the advanced portion $t_2 - t_1 < 0$ should only take into account the negative energy eigenstates (holes).