The following identity is used in Peskin & Schroeder's book Eq.(19.43), page 660:
$$\int\frac{d^4k}{(2\pi)^4}\,\frac{1}{(k^2)^2}e^{ik\cdot\epsilon}=\frac{i}{(4\pi)^2}\log\frac{1}{\epsilon^2},\quad \epsilon\rightarrow 0$$
I can't figure out why it holds. Could someone provide a method to prove this? Many thanks in advace.
That's equivalent simply to $c\int dx/x$. Switch to the Euclidean spacetime, $k_0=ik_4$ where $(k_1,\dots k_4)$ is $k_E$; i.e. analytically continue in $k_0$ (Wick rotation). The integral is $$\int \frac{i\cdot d^4 k_E}{(2\pi)^4} \frac{1}{(k_E^2)^2} \exp(ik\cdot \epsilon)$$ So it's proportional to the Fourier transform of $1/k_E^4$. The original function is $SO(4)$ symmetric, so the Fourier transform must be symmetric as well and depend on $\epsilon^2$ only. Dimensional analysis implies that the result is dimensionless i.e. it must be a combination of a constant and $\ln(\epsilon^2)$. The logarithm is there with a nonzero coefficient so the constant only determines how to take the logarithm: it should properly be written as $\ln(\epsilon^2/\epsilon_0^2)$ for some constant $\epsilon_0$ with the same dimension.
The only remaining unknown is the coefficient and one gets $4\pi^2$ from the remaining integral. It's a sort of waste of resources to compute this special integral; it's better to compute the more general integrals in appendix A.4, see especially formulae (A.44)-(A.49) on page 807, which I won't copy here because that's why Peskin and Schroeder wrote the textbook.
I will give another approach to this identity. First, we notice that
$$\int\frac{d^4k}{(2\pi)^4}\,\frac{1}{(k^2-m^2)^2}e^{ik\cdot\epsilon}=-i\frac{\partial}{\partial m^2}D_F(x)\big|_{x=\epsilon}$$
For space-like vector $\epsilon^2=-r^2<0$, we have
$$D_F(x)=\frac{m}{4\pi^2r}K_1(mr)$$
whose derivation refers to Weinberg's book vol. 1, page 202. For $r\rightarrow 0$, the following expansion holds
$$ K_1(mr)=\frac{1}{mr}+\frac{mr}{2}\log\frac{mr}{2}$$
With this conditions, we finally obtain
$$\int\frac{d^4k}{(2\pi)^4}\,\frac{1}{(k^2-m^2)^2}e^{ik\cdot\epsilon}=\frac{i}{16\pi^2}\log\frac{1}{\epsilon^2}$$
Further another approach.
After Wick rotation($k^0=ik_E^0,\,\epsilon^0=i\epsilon_E^0$) the integral is
$$ I_1 \equiv \int \frac{i\cdot d^4 k_E}{(2\pi)^4} \frac{1}{(k_E^2)^2} \exp(-ik_E\cdot \epsilon_E) . $$
In (19.43) we want the quntity
$$ \frac{\partial}{\partial \epsilon^{\gamma}}I_1. $$
Note that
$$ \frac{1}{k_E^2}=\int_0^{\infty}du e^{-k_E^2 u}, $$
$$ \frac{1}{(k_E^2)^2}=\int_0^{\infty}du\int_0^{\infty}dv e^{-k_E^2 (u+v)}, $$
$$ I_1=\frac{i}{(2\pi)^4} \int d^4k_E \int_0^{\infty}du\int_0^{\infty}dv e^{-(u+v)k_E^2-i\epsilon_E \cdot k_E}. $$
$ \displaystyle \int_{-\infty}^{\infty} dk_E^i \exp[-(u+v)(k_E^i)^2-i\epsilon_E^i k_E^i]
= \sqrt{\frac{\pi}{u+v}} \exp\left[ -\frac{(\epsilon_E^i)^2}{4(u+v)} \right] $
\begin{alignat}{2}
\therefore I_1&=&& \frac{i}{(2\pi)^4} \int_0^{\infty}du\int_0^{\infty}dv \frac{\pi^2}{(u+v)^2} \exp\left[-\frac{\epsilon_E^2}{4(u+v)} \right] \\
&=&& \frac{i}{16\pi^2}I_2 \left(\frac{\epsilon_E^2}{4}\right)
\end{alignat}
where $\displaystyle I_2(x) \equiv \int_0^{\infty}du\int_0^{\infty}dv\, \frac{1}{(u+v)^2}\exp\left(-\frac{x}{u+v}\right) $. The calculation of $I_2$ is here.
\begin{alignat}{2}
I_1 &=&& \frac{i}{16\pi^2}\left( -\log\left(\frac{\epsilon_E^2}{4}\right) +\gamma-1+\lim_{M\to \infty}\log M +\cal{O}(\epsilon_E^2) \right) \\
&=&& \frac{i}{16\pi^2}\left( -\log\left(-\frac{\epsilon^2}{4}\right) +\gamma-1+\lim_{M\to \infty}\log M +\cal{O}(\epsilon^2) \right)
\end{alignat}
After $\epsilon \to 0$, we have
$$ \frac{\partial}{\partial \epsilon^{\gamma}}I_1 = \frac{i}{16\pi^2}\frac{\partial}{\partial \epsilon^{\gamma}}\log \frac{1}{\epsilon^2}. $$
Here is yet another solution, which probably is not a physicist's way of thinking.
After Wick rotation, we may work on Euclidean space. The function $f(k_E) = |k_E|^{-4}$ is not square-integrable on $\mathbb{R}^4$, however, so its Fourier transform does not exist in ordinary sense. A moment of thought suggests that it can be realized as distribution on the space
$$ \mathcal{A} := \{ \varphi \in \mathcal{S}(\mathbb{R}^4) : \textstyle \int_{\mathbb{R}^4} \varphi(\epsilon) \, \mathrm{d}^4\epsilon = 0 \}, $$
where $\mathcal{S}(\mathbb{R}^4)$ is the Schwarz space. Then computing
$$I(\epsilon) = i \check{f}(\epsilon) = \frac{i}{(2\pi)^4} \int_{\mathbb{R}^4} |k_E|^{-4} e^{i\epsilon \cdot k_E} \, \mathrm{d}^4 k_E$$
in distribution sense amounts to identifying the following pairing
$$ \langle I, \varphi \rangle = \langle \check{f}, \varphi \rangle = \langle f, \check{\varphi} \rangle, \qquad \forall \varphi \in \mathcal{A}. $$
Since $\check{\varphi}(0) = 0$ and $\check{\varphi}$ has rapid decay near infinity, the pairing $\langle f, \check{\varphi} \rangle$ is realized as Lebesgue integral. Then by the Fubini's theorem,
\begin{align*} \langle f, \check{\varphi} \rangle &= i \int_{\mathbb{R}^4} \frac{1}{|k_E|^4} \check{\varphi}(k_E) \, \mathrm{d}^4 k_E = i \int_{\mathbb{R}^4} \bigg( \int_{0}^{\infty} t \mathrm{e}^{-|k_E|^2 t} \, \mathrm{d}t \bigg) \check{\varphi}(k_E) \, \mathrm{d}^4 k_E, \\ &= i \int_{0}^{\infty} t \bigg( \int_{\mathbb{R}^4} \mathrm{e}^{-t |k_E|^2} \check{\varphi}(k_E) \, \mathrm{d}^4 k_E \bigg) \mathrm{d}t. \end{align*}
Using $\langle \mathrm{e}^{-t|\cdot|^2}, \check{\varphi} \rangle = \langle (\mathrm{e}^{-t|\cdot|^2})^{\vee}, \varphi \rangle$ and the formula $\int_{\mathbb{R}} \mathrm{e}^{-tx^2}\mathrm{e}^{ix\epsilon} \, dx = \sqrt{\frac{\pi}{t}} \mathrm{e}^{-\epsilon^2/4t}$ for $t > 0$, we have
\begin{align*} \langle f, \check{\varphi} \rangle &= i \int_{0}^{\infty} t \bigg( \int_{\mathbb{R}^4} \frac{1}{(4\pi t)^2} \mathrm{e}^{-\frac{|\epsilon|^2}{4t}} \varphi(\epsilon) \, \mathrm{d}^4 \epsilon \bigg) \mathrm{d}t \\ &= \frac{i}{(4\pi)^2} \int_{0}^{\infty} \frac{1}{t} \bigg( \int_{\mathbb{R}^4} \mathrm{e}^{-\frac{|\epsilon|^2}{4t}} \varphi(\epsilon) \, \mathrm{d}^4 \epsilon \bigg) \mathrm{d}t. \end{align*}
We want to finalize the computation by switching the order of integration, but the Fubini's theorem is not applicable in this case and even the heuristic computation produces a divergent integral. Thankfully, using the fact that $\int_{\mathbb{R}} \varphi(\epsilon) \, \mathrm{d}^4\epsilon = 0$, we can regularize the inner integral so that the Fubini's theorem works:
\begin{align*} \langle f, \check{\varphi} \rangle &= \frac{i}{(4\pi)^2} \int_{0}^{\infty} \frac{1}{t} \bigg( \int_{\mathbb{R}^4} \big( \mathrm{e}^{-\frac{|\epsilon|^2}{4t}} - \mathbf{1}_{ \{ t \geq 1 \} } \big) \varphi(\epsilon) \, \mathrm{d}^4 \epsilon \bigg) \mathrm{d}t \\ &= \frac{i}{(4\pi)^2} \int_{\mathbb{R}^4} \bigg( \int_{0}^{\infty} \frac{1}{t}\big( \mathrm{e}^{-\frac{|\epsilon|^2}{4t}} - \mathbf{1}_{ \{ t \geq 1 \} } \big) \, \mathrm{d}t \bigg) \varphi(\epsilon) \, \mathrm{d}^4 \epsilon. \end{align*}
Now the inner integral can be computed using the substitution $u = |\epsilon|^2/4t$ as follows
\begin{align*} \int_{0}^{\infty} \frac{1}{t}\big( \mathrm{e}^{-\frac{|\epsilon|^2}{4t}} - \mathbf{1}_{ \{ t \geq 1 \} } \big) \, \mathrm{d}t &= \int_{0}^{\infty} \frac{1}{u} \big( \mathrm{e}^{-u} - \mathbf{1}_{ \{ u \leq |\epsilon|^2/4 \} } \big) \, \mathrm{d}u = \log\frac{4}{|\epsilon|^2} -\gamma. \end{align*}
Here, $\gamma$ is the Euler-Mascheroni constant. Therefore it follows that
$$ I(\epsilon) = \frac{i}{(4\pi)^2} \bigg( \log\frac{4}{|\epsilon|^2} -\gamma \bigg) = \frac{i}{(4\pi)^2} \log\frac{1}{\epsilon^2}. $$
The last equality follows from the fact that constants as distribution on $\mathcal{A}$ is equal to zero, i.e., for any constant $c$ we have $\langle c, \varphi \rangle = 0$ for all $\varphi \in \mathcal{A}$.
