\section*{Appendix}

The purpose of this Appendix is to give an explicit calculation of
the general integral (\ref{integral}).
Performing multiplication, one can get
$$
G(y_1,y_2,\epsilon)=\int
d^dx\varepsilon^{\rho+\sigma} 
\frac{\varepsilon^2+(x-y_1)_k\gamma^k(x-y_2)_j\gamma^j-\varepsilon(y_1-y_2)_k\gamma^k\gamma_0}
{(\varepsilon^2+(\vec{x}-\vec{y}_1)^2)^{d/2+1/2+\sigma}
(\varepsilon^2+(\vec{x}-\vec{y}_2)^2)^{d/2+1/2+\rho}}
\times$$ $$
\times\left(\delta_{mn}+4\frac{(x-y_1)_i(x-y_2)^i(x-y_1)_m(x-y_2)_n}
{(\varepsilon^2+(\vec{x}-\vec{y}_1)^2)
(\varepsilon^2+(\vec{x}-\vec{y}_2)^2)}\right.$$ \begin{equation}\left.
-2\frac{(x-y_1)_m(x-y_1)_n}
{\varepsilon^2+(\vec{x}-\vec{y}_1)^2}-2\frac{(x-y_2)_m(x-y_2)_n}
{\varepsilon^2+(\vec{x}-\vec{y}_2)^2}\right).
\label{fullintegral}
\end{equation}

We will consider (\ref{fullintegral}) as a sum of different terms
with the similar structure, it will unify our investigation.
Also the succeeding analysis will be simplified if we rewrite each term
applying the Fourier transformation.
We use the following Fourier transformation formula:
\begin{equation}
\label{fur1}
2^{\alpha}\Gamma\left(\frac d2 +\alpha\right)\frac{\varepsilon^{\alpha}}{(\varepsilon^2+(\vec{x}-\vec{y})^2)^{d/2+\alpha}}=
\mbox{const}\cdot
\int
d^dpe^{i\vec{p}(\vec{x}-\vec{y})}p^{\alpha}K_{\alpha}(\varepsilon
p)
\end{equation}
where $K_{\alpha}$ is a modified Bessel function with index $\alpha$ and $p=|\vec{p}|$.
The contraction of two propagators in (\ref{fullintegral}) turns into the
product of their Fourier images. The remaining integral could be found
explicitly in \cite{Prud}, so we get
\begin{eqnarray}
&&\int e^{-i\vec{p}(\vec{y_1}-\vec{y_2})}
p^{\mu+\nu}K_{\mu}(\varepsilon p)K_{\nu}(\varepsilon p)
d^dp=\nonumber\\
&&=\frac{\mbox{const}}{\varepsilon^{d+\mu+\nu}}
\frac{2^{d/2+\mu+\nu-2}\Gamma\left(d/2+\mu+\nu\right)
      \Gamma\left(d/2+\nu\right)
      \Gamma\left(d/2+\mu\right)}
     {\Gamma\left(d+\mu+\nu\right)}
      \times\nonumber\\
&&_3F_2\left(d/2+\mu+\nu,
           d/2+\nu,
           d/2+\mu;
           \frac{d+\mu+\nu}2,
           \frac{d+\mu+\nu+1}2;
           -\frac{|\vec y_1-\vec y_2|^2}{4\varepsilon^2}\right).
\label{genint}
\end{eqnarray}
All factors like $p_m$ appearing after differentiation
of (\ref{fur1})
could be obtained by applying the operator $i\frac{\partial}{\partial
y_m}$ to (\ref{genint}).  The circle of convergence for the
hypergeometric function $_3F_2(z)$ is $|z|<1$.  So, in order to find
the limit $\varepsilon\rightarrow 0$ we have to apply the formula of
analytic continuation, known in the theory of the hypergeometric
functions $$ _{3}F_2(a_1,a_2,a_3;b_1,b_2;z)=
\frac{\Gamma(b_1)\Gamma(b_2)}{\Gamma(a_1)\Gamma(a_2)\Gamma(a_{3})}
\left[\Gamma(a_1)(-z)^{-a_1}\frac{\Gamma(a_2-a_1)\Gamma(a_3-a_1)}
{\Gamma(b_1-a_1)\Gamma(b_2-a_1)}\times\right.$$ $$
{}_3F_2(1+a_1-b_1,1+a_1-b_2,a_1;1+a_1-a_2,1+a_1-a_3;1/z)$$
$$+\Gamma(a_2)(-z)^{-a_2}\frac{\Gamma(a_1-a_2)\Gamma(a_3-a_2)}
{\Gamma(b_1-a_2)\Gamma(b_2-a_2)}\times$$ $$
\times{}_3F_2(1+a_2-b_1,1+a_2-b_2,a_2;1+a_2-a_1,1+a_2-a_3;1/z)$$ $$
+\Gamma(a_3)(-z)^{-a_3}\frac{\Gamma(a_1-a_3)\Gamma(a_2-a_3)}
{\Gamma(b_1-a_3)\Gamma(b_2-a_3)}\times$$ \begin{equation}
\left.\times\frac{}{}{}_3F_2(1+a_3-b_1,1+a_3-b_2,a_3;1+a_3-a_1,1+a_3-a_2;1/z)\right].
\label{invert}
\end{equation}
Taking into account, that $F(0)=1$ we find three terms with the different
behavior as $\varepsilon$ tends to zero. In our concrete case of the integral
(\ref{fullintegral}) we have $\mu=\sigma+1/2-n$, $\nu=\rho+1/2-m$,
where $m$, $n$ are integers, which are different for each term in
 (\ref{fullintegral}). So we find that three mentioned types of
behavior as $\varepsilon\rightarrow 0$
are
\begin{equation}
\varepsilon^{\sigma+\rho+d+1-n-m},\qquad
\varepsilon^{\sigma-\rho-n+m},\qquad \varepsilon^{\rho-\sigma+n-m}.
\end{equation}
Note that the powers in the second and the third types are inverse
to each other, so if the difference between $\rho$ and $\sigma$  is not
integer we have a singular contribution.

