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?