The analytical result for free massless fermion propagator For massless fermion, the free propagator in quantum field theory
is
\begin{eqnarray*}
 &  & \langle0|T\psi(x)\bar{\psi}(y)|0\rangle=\int\frac{d^{4}k}{(2\pi)^{4}}\frac{i\gamma\cdot k}{k^{2}+i\epsilon}e^{-ik\cdot(x-y)}.
\end{eqnarray*}
In Peskin & Schroeder's book, An introduction to quantum field theory
(edition 1995, page 660, formula 19.40), they obtained the analytical
result for this propagator,
\begin{eqnarray*}
 &  & \int\frac{d^{4}k}{(2\pi)^{4}}\frac{i\gamma\cdot k}{k^{2}+i\epsilon}e^{-ik\cdot(x-y)}=-\frac{i}{2\pi^{2}}\frac{\gamma\cdot(x-y)}{(x-y)^{4}} .\tag{19.40}
\end{eqnarray*}
Question: Is this analytical result right? Actually I don't know
how to obtain it.
 A: Yes it is correct. The derivation in P&S is straightforward but I will expand on it a bit. The key observation is that
\begin{equation}
\int\frac{d^4k}{(2\pi)^4}e^{-ik\cdot(y-z)}\frac{i\gamma^{\mu}k_{\mu}}{k^2+i\epsilon}
=-\gamma^{\mu}\partial_{\mu}\int\frac{d^4k}{(2\pi)^4}\frac{1}{k^2+i\epsilon}e^{-ik\cdot(y-z)},
\end{equation}
where the integral on the right hand side is the Feynman propagator for a 
massless scalar. Performing the $k$ integrals to get to position space yields
\begin{equation}
\int\frac{d^4k}{(2\pi)^4}\frac{1}{k^2+i\epsilon}e^{-ik\cdot(y-z)}=\frac{i}{4\pi^2}\frac{1}{(y-z)^2-i\epsilon}.
\end{equation}
If you aren't sure about this last step, it is easier to consider the massive case and then take the limit as $m\rightarrow 0$ at the end. Schwinger parameters are also helpful for proving this identity. 
Once you have transformed to position space you simply act with $-\gamma^{\mu}\partial_{\mu}$ to arrive at the final expression.
A: As alluded to in the other answer here, the integral can basically be evaluated in several ways. One of them is the one that OP has himself followed. (PS - OP, congratulations on completing that feat!)
Let me present here a way of computing this using Schwinger parameterization. We will use
$$
\frac{1}{a} = \int_0^\infty d\tau e^{- \tau a}  ~, a > 0~. 
$$
We want to compute 
$$
I = \int \frac{d^4k}{(2\pi)^4} \frac{i  e^{- i k \cdot x}}{k^2+ i \epsilon}  = \int \frac{dk^0}{2\pi} \frac{d^3k}{(2\pi)^3} \frac{i  e^{- i k^0 t + i \vec{k} \cdot \vec{x} }}{(k^0)^2 - \vec{k}^2+ i \epsilon} 
$$
Do a Wick rotation $k^0 \to i k^0_E$, $t \to - i t_E$. Then, 
$$
I  =   \int \frac{dk_E^0}{2\pi} \frac{d^3k}{(2\pi)^3} \frac{ e^{ k_E^0 t + i \vec{k} \cdot \vec{x} }}{(k_E^0)^2 + \vec{k}^2+ i \epsilon} =   \int   \frac{d^4k}{(2\pi)^3} \frac{ e^{ - i k_E^0 t_E + i \vec{k} \cdot \vec{x} }}{ k^2 }
$$
where in the last equation, we now have a Euclidean $k^2$ that is always positive over the range of integration. Now, we may use the Schwinger parameterization so that
$$
I = \int_0^\infty d\tau \int \frac{d^4k}{(2\pi)^4}  e^{- k^2 \tau - i k_E^0 t_E + i \vec{k} \cdot \vec{x}} 
$$
Now, we can do the integral of $k$ quite easily since $k^2 = \sum_i k_i^2$. This gives
$$
I =  \int_0^\infty d\tau \frac{1}{(4\pi)^2} \frac{e^{ - \frac{1}{4\tau} ( t_E^2 + \vec{x}^2 ) } }{\tau^2} 
$$
Now to perform the integral over $\tau$, define new integration variable $y = \frac{1}{4\tau} $. Then
$$
I =  \int_0^\infty dy \frac{1}{(2\pi)^2}  e^{  -  y ( t_E^2 + \vec{x}^2 ) } 
$$
This last integral is again the Schwinger parameter one. It converges and is nice so we compute it and find
$$
I = \frac{1}{4\pi^2 ( t_E^2 + \vec{x}^2 ) } = - \frac{1}{4\pi^2 ( t^2 - \vec{x}^2 ) } = - \frac{1}{4\pi^2 x^2 } ~. 
$$
where in the last step, we have performed the inverse Wick rotation to go back to Lorentzian time. 
A: In the following, I will carefully deal with Wick rotation. In the
end, I have found that I was confused.
The integration is
\begin{eqnarray*}
 &  & \int\frac{d^{4}k}{(2\pi)^{4}}\frac{1}{k^{2}+i\epsilon}e^{-ik\cdot x}\\
 & = & \frac{1}{(2\pi)^{4}}\int d^{3}ke^{i\mathbf{k}\cdot\mathbf{x}}\int_{-\infty}^{\infty}dk_{0}\frac{1}{k_{0}^{2}-(E_{k}-i\epsilon)^{2}}e^{-ik_{0}t}\\
 & \equiv & \frac{1}{(2\pi)^{4}}\int d^{3}ke^{i\mathbf{k}\cdot\mathbf{x}}\times\mathrm{I}
\end{eqnarray*}
with
\begin{eqnarray*}
\mathrm{I} & = & \int_{-\infty}^{\infty}dk_{0}\frac{1}{k_{0}^{2}-a^{2}}e^{-ik_{0}t}\\
a & = & E_{k}-i\epsilon=\sqrt{m^{2}+\mathbf{k}^{2}}-i\epsilon
\end{eqnarray*}
Now we will use Wick rotation to calculate $\mathrm{I}$. Note that
$\pm a$ are two singularities of the integrand. Consider following
contour. The radii of coutours $l_{5},l_{6}$ are both $R$ and $R\rightarrow\infty$.

According to contour integral theorem, we can see
\begin{eqnarray*}
\mathrm{I} & = & \int_{-\infty}^{\infty}dk_{0}\frac{1}{k_{0}^{2}-a^{2}}e^{-ik_{0}t}\\
 & = & \int_{l_{1}}dz\frac{1}{z^{2}-a^{2}}e^{-izt}+\int_{l_{2}}dz\frac{1}{z^{2}-a^{2}}e^{-izt}\\
 & = & \bigg(\int_{l_{5}}dz\frac{1}{z^{2}-a^{2}}e^{-izt}+\int_{l_{3}}dz\frac{1}{z^{2}-a^{2}}e^{-izt}\bigg)+\bigg(\int_{l_{4}}dz\frac{1}{z^{2}-a^{2}}e^{-izt}+\int_{l_{6}}dz\frac{1}{z^{2}-a^{2}}e^{-izt}\bigg)\\
 &  & \bigg[\text{note: set }z=ik_{E}^{0}\text{ in }l_{3},l_{4}\text{ and combine }l_{5},l_{6}\bigg]\\
 & = & (-i)\int_{-\infty}^{\infty}dk_{E}^{0}\frac{1}{(k_{E}^{0})^{2}+a^{2}}e^{tk_{E}^{0}}+\int_{l_{6}}dz\frac{1}{z^{2}-a^{2}}(e^{-izt}+e^{izt}),\ \bigg[\text{set }z=Re^{i\phi}\text{ in }l_{6}\bigg]\\
 & = & (-i)\int_{-\infty}^{\infty}dk_{E}^{0}\frac{1}{(k_{E}^{0})^{2}+a^{2}}e^{tk_{E}^{0}}\\
 &  & -iR\int_{0}^{\frac{\pi}{2}}d\phi e^{i\phi}\frac{1}{R^{2}e^{2i\phi}-a^{2}}(e^{-itR\cos\phi+tR\sin\phi}+e^{itR\cos\phi-tR\sin\phi})\\
 & \equiv & (-i)\int_{-\infty}^{\infty}dk_{E}^{0}\frac{1}{(k_{E}^{0})^{2}+a^{2}}e^{tk_{E}^{0}}+\mathrm{II}
\end{eqnarray*}
with
\begin{eqnarray*}
\mathrm{II} & = & -iR\int_{0}^{\frac{\pi}{2}}d\phi e^{i\phi}\frac{1}{R^{2}e^{2i\phi}-a^{2}}(e^{-itR\cos\phi+tR\sin\phi}+e^{itR\cos\phi-tR\sin\phi})
\end{eqnarray*}
Actually, I do not know how to prove $\mathrm{II}=0$ as $R\to\infty$.
But if $\mathrm{II}\neq0$ as $R\to\infty$, then we can not simply
obtain
\begin{eqnarray*}
\int_{-\infty}^{\infty}dk_{0}\frac{1}{k_{0}^{2}-a^{2}}e^{-ik_{0}t} & = & (-i)\int_{-\infty}^{\infty}dk_{E}^{0}\frac{1}{(k_{E}^{0})^{2}+a^{2}}e^{tk_{E}^{0}}
\end{eqnarray*}
I am just confused at this point.
A: The calculation of the propagator in four dimensions is as follows. 
\begin{eqnarray*}
\int\frac{d^4 k}{(2\pi)^4}e^{-ik\cdot (x-y)}\frac{1}{k^2} 
&=&  i\int \frac{d^4 k_E}{(2\pi)^4}e^{ik_E\cdot (x_E-y_E)}\frac{1}{-k_E^2}  \\
&=& \frac{-i}{(2\pi)^4} \left( \int_0^{2\pi}d\theta_3 \int_0^{\pi}d\theta_2 \sin \theta_2 \right) \int_0^{\infty} dk_E k_E^3 \frac{1}{k_E^2} \int_0^{\pi}d\theta_1 \sin^2 \theta_1 e^{ik_E | x_E-y_E | \cos \theta_1} \\
&=& \frac{-i4\pi}{(2\pi)^4}  \int_0^{\infty} dk_E k_E \int_0^{\pi}d\theta_1 \frac{1-\cos 2\theta_1}{2} e^{ik_E | x_E-y_E | \cos \theta_1} \\
&=& \frac{-i}{4\pi^3} \frac{1}{| x_E-y_E |^2} \int_0^{\infty} ds s (\frac{\pi}{2} J_0(s)- \frac{\pi i^2}{2} J_2(s)) 
\end{eqnarray*}
where $s\equiv k_E\| x_E-y_E \| $, and $J_n(s)$'s are bessel functions and I made use of Hansen-Bessel Formula.
\begin{eqnarray*}
&=& \frac{-i}{4\pi^3} \frac{1}{| x_E-y_E |^2}  \int_0^{\infty} ds s \frac{\pi}{2} \frac{2}{s} J_1(s) \\
&=& -\frac{i}{4\pi^2} \frac{1}{| x_E-y_E|^2} \int_0^{\infty} ds \, J_1(s) \\
&=& -\frac{i}{4\pi^2} \frac{1}{| x_E-y_E |^2} \\
&=& \frac{i}{4\pi^2} \frac{1}{(x-y)^2}
\end{eqnarray*}
