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?

 A: maybe you already thought about it, but I'll ask anyway: why don't you try to tackle directly the differential operator: $\partial_t -\vec{\sigma}\cdot\nabla$ ? It's a first order PDE so you would get something of the kind:
$$
\psi(x,t) = \psi_0(x,t) + \int G(x,x',t,t')f(x',t')\, dx'dt' \qquad \text{(1)}
$$
Where $\psi_0(x,t)$ is a solution of the homogeneous equation and $G(x,x',t,t')$ is the solution of
$$
(\partial_t -\vec{\sigma}\cdot\nabla)G(x,x',t,t') = \delta(x-x')\delta(t-t')
$$
Now, it can be shown that equation 1 can be cast into the form (see Ref.):
$$
\psi(x,t) = \psi_0(x,t) + \int dx' \int_{-\infty}^t G_1(x,x',t,t')f(x',t')\, dt' \qquad \text{(2)}
$$
where $G_1(x,x',t,t')$ is an object that satisfy:
$$
(\partial_t -\vec{\sigma}\cdot\nabla)G_1(x,x',t,t') = 0
$$
and therefore should be less clumsy to solve.
Reference: Mathematical Method of Classical and Quantum Physics, F. W. Byron and R. W. Fuller. Chapter 7: Time-dependent Green's functions: First Order
A: 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).
