Spinor field normalisation from poles in the propagator In the theory of free scalar bosons (KG field) it is a basic result that the propagator $\Delta(p)$ has poles at $p^2=m^2$, with residue $1$ (or any other constant, depending on conventions). Thinking of $p^2$ as a variable is acceptable because it can be shown that the propagator is a function of $p^2$, even in the interacting case.
On the other hand, in the theory of free spinor fields (Dirac field), we can show that the propagator is not a function of $p^2$, but of $\not p$, so we can't use the $p^2$ trick. Conventions aside, the propagator is
$$
S(p)=\frac{\not p+m}{p^2-m^2+i \epsilon}
$$
It is clear that we cannot think of poles as $p^2\to m^2$, because $S(p)$ depends on $p$ in a non-so-trivial way. The usual solution is to think of $S$ as a function of $\not p$, as if $\not p$ was a number instead of a matrix, thus writing
$$
S(\not p)=\frac{1}{\not p-m+i\epsilon}
$$
and, again, we find a pole at $\not p=m$, with residue $1$. I think this is not mathematics. This is just nonsensical to me (may be I'm too skeptic, and we can give a meaning to this last formula).
In an interacting theory, we define the field strength normalisation $Z$ as the residue of the propagator at the poles:
$$
\Delta(p^2)=\frac{Z}{p^2-m^2+i\epsilon}+\int_{M_\text{thresh}^2}^\infty \mathrm d\mu^2\rho(\mu^2)\frac{1}{p^2-\mu^2+i\epsilon}
$$
$$
S(\not p)=\frac{Z}{\not p-m+i\epsilon}+\int_{M_\text{thresh}^2}^\infty \mathrm d\mu^2\frac{\not p\rho_1(\mu^2)+\mu \rho_2(\mu^2)}{\not p^2-\mu^2+i\epsilon}
$$
Now, I don't understand what's the proper definition of $Z$ as a residue. In what sense is it a residue? What's the variable? I just can't accept it is a residue as $\not p\to m$. Do we really take this "think of $\not p$ as a variable'' seriously? Is it possible to formalise these matters, by thinking of a residue for an actual variable (as $p_0$)?
Perhaps the $\not p$ trick is just that: a trick, which may simplify calculations, but such that it can be shown that the actual result, found by a more standard procedure, is the same. Is this the case? If so, what is the correct procedure?
(I would really appreciate if the answers don't assume that we can boost to the rest frame of the particle, as I'd like it to be as general as posible. I want to take into account the possibility of $m=0$, so please don't boost into $k=(m,\boldsymbol 0)$ if it's not really necessary)
 A: First of all, stating that the Klein-Gordon propagator is a function of $p^2$ is stating that it is a function of $p^\mu$ but does not depend on Lorentz transforms of the vector $p^\mu$, so it must be in fact a function in the form $S(p^2(p^\mu))$. However, that the propagator has a pole at $p^2=m^2$ is only true if we understand $p^2$ as an independent c-number. 
In the full $p^\mu$ space it is harder to naturally define what this statement means because $p^\mu p^\nu \eta_{\mu \nu}$ fails as a "radial" coordinate on the $p^2=0$ light-cone. (Note that we do integrate in light-like and space-like $p^\mu$ space upon propagation). The only way to give meaning to the statement of a pole is to understand $p^2$ as a locally (but not globally) valid coordinate in the $p^\mu$ space as given by $p^\mu=\sqrt{p^2} n^\mu,\,n^\mu n_\mu=1$, where $n^\mu$ is characterized by three "angular" coordinates. In the sense of this local coordinate and the complex continuation through it, the propagator has a pole of residue $1$ (in your phase convention) at $p^2=m^2$ same for every value of $n^\mu$.

Now let's turn to the fermionic case. The matrix $\not p=p^\mu \gamma_\mu$ is Lorentz-covariant and so is the matrix 
$$S(\not p) = \frac{\not p + m}{p^2-m^2}$$ 
The notation $\frac{1}{\not p - m}$ only means the inverse of the matrix in the denominator; it indeed is the only sense it can be interpreted in. The propagator is, similarly to the scalar propagator, covariant with respect to Lorentz coordinate transformations. But only covariant, so we cannot really say it is a function of just $p^2$. $\not p$ contains the same information as $p^\mu$ so it is interchangeable to say that $S$ is a function of $\not p$ or $p^\mu$.
Now let us use the local coordinates $p^\mu=\sqrt{p^2} n^\mu,\,n^\mu n_\mu=1$. You can verify yourself that once again the propagator has a pole in the very sense the scalar propagator had, albeit here the value of the residue is $m(1+\not n)=m+\not p_{on-shell}$, i.e. it is dependent on the "angular" coordinates in the $p^\mu$ space.
The renormalization factor $Z$ is then simply the extra factor to this residue.

The case of massless particles is more difficult because a light-like four-momentum $k^\mu$ corresponds to the same physical state as $\lambda k^\mu$ with any $\lambda>0$ and $k^2$ does nothing to discriminate. One must then choose a certain coordinate such that $k^2$ is a "good" coordinate at least in the vicinity of the light-cone, and then the pole is given as the pole of $k^2$ in this sense. (An example would be classical 4D spherical coordinates in the $k^\mu$ space with one of the north poles at $k^\mu=(k^0,0,0,0)$ and $k^2$ taking the role of the $\theta$ coordinate.)
One must also be careful that the value of the residue generally depends on the gauge employed (the $m \to 0$ limit gives you the propagator only in a specific gauge), but otherwise $Z$ can be read of from the residue in the very same manner as in the massive-particle case.
A: Very good question, OP! Your scepticism is most certainly justified. The good news is, someone already has addressed your concerns. You can find the answer in Ticciati's Quantum Field Theory for Mathematicians, section 10.13.
Long story short: given
\begin{equation}
S(\not p)=\frac{1}{\not p-m-\Sigma(\not p)+i\epsilon}
\end{equation}
you can always parametrise the matrix $\Sigma$ as
\begin{equation}
\Sigma(\not p)=a(p^2)1+b(p^2)\not p
\end{equation}
for a pair of scalar functions $a,b$. This general expression is the result of Lorentz- and parity-invariance. But you already know that.
With this, you can rationalise the denominator into
\begin{equation}
S(\not p)=\frac{i(\not p+\alpha)}{(1-b)(p^2-\alpha^2)+i\epsilon}
\end{equation}
where
\begin{equation}
\alpha(p^2)\overset{\mathrm{def}}=\frac{m+a(p^2)}{1-b(p^2)}
\end{equation}
The required pole at $p^2=m^2$ implies
\begin{equation}
{\color{red}{\alpha(m^2)=m}},\tag1
\end{equation}
that is,
\begin{equation}
S(\not p)=\frac{1}{(1-b)(1-2m\alpha')}\frac{i(\not p+m)}{p^2-m^2+i\epsilon}+\mathcal O(1)
\end{equation}
If the residue is to be equal to that of a normalised field, that is,
\begin{equation}
S(\not p)=\frac{i(\not p+m)}{p^2-m^2+i\epsilon}+\mathcal O(1)
\end{equation}
then we must have
\begin{equation}
{\color{red}{(1-b(m^2))(1-2m\alpha'(m^2))=1}}\tag2
\end{equation}
So far so good: we have rigorously characterised the normalisation conditions of a bispinor field.
The key point is the following: if we introduce the formal complex variable $\not p\in\mathbb C$, we can formally write $(1),(2)$ as
\begin{equation}
\begin{aligned}
\Sigma(m)&=0\\
\Sigma'(m)&=0
\end{aligned}
\end{equation}
as can be checked by a straightforward computation. This justifies the introduction of the formal complex variable $\not p$: the formal manipulations of this variable are equivalent to the more correct procedure of introducing the pair of scalar functions $a,b$, and the end result is the same. All's well that ends well I guess.
A: Consider the simple case of time-like momentum in 2+1 dimensions so we can go to a frame where $p = (\omega,0,0)$.  Choose $\gamma^t = \sigma_z$ to be the third Pauli matrix, $\sigma_z = diag(1,-1)$.  Then the propagator for the electron can be written
$$
S(p) = {1 \over \omega^2 - m^2} (\omega \sigma_3 - m) \ .
$$
In this simple case, it's pretty clear what's going on.  $S(p)$ is a 2x2 (diagonal) matrix with eigenvalues $1/(\omega+m)$ and $-1/(\omega-m)$.  In other words (and I believe more generally) the eigenvalues of $S(p)$ have poles when $p^2 = m^2$ with residue $\pm 1$.
