D'Alembertian Green's function from Laplacian Green's function by Wick rotation Following this review paper (in particular eq.(14)), I am trying to understand how to obtain the Green's for the D'Alembert operator from the kernel of the Laplace operator by ''going to imaginary time''.
I shall denote by $x=(x^0,x^i)$, with $i=1,2,\ldots,n-1$, the usual Minkowski coordinates for $\mathbb R^{1,n-1}$ and by $x_E=(x_E^0, x^i)$ the cooresponding coordinates for the Euclidean space $\mathbb R^{n}$. The D'Alembert and Laplace operators are then defined respectively as 
$$
\Box = \left(\frac{\partial}{\partial x^0}\right)^2-\sum_{i=1}^{n-1}\left(\frac{\partial}{\partial x^i}\right)^2\,,\qquad
\Delta = \left(\frac{\partial}{\partial x_E^0}\right)^2+\sum_{i=1}^{n-1}\left(\frac{\partial}{\partial x^i}\right)^2\,.
$$
The heuristic observation underlying the Wick rotation trick is that $\Delta = - \Box$ if we formally identify $x_E^0=ix^0$, namely
$$
\Delta \varphi(x^0_E,x^1,\ldots,x^{n-1})=-\Box \varphi(ix^0,x^2,\ldots,x^{n-1}).
$$ 
In fact, $-x^2 = x_E^2$ under this identification, where $x^2 = (x^0)^2-(x^1)^2-\cdots- (x^{n-1})^2$ and $x_E^2$ is the Euclidean squared norm.
Then starting from the Laplace Green's function (I employ the Einstein summation convention)
$$
\Delta \left[(x_E^0)^2+(x^i)^2 \right]^{1-n/2}=\frac{2\pi^{n/2}(2-n)}{\Gamma(n/2)}\delta(x_E^0,x^1,\ldots,x^{n-1})
$$
one might expect that
$$
\Box \left[-(x^0)^2+(x^i)^2 \right]^{1-n/2}=\frac{2i\pi^{n/2}(2-n)}{\Gamma(n/2)}\delta(x^0,x^1,\ldots,x^{n-1})\,,
$$
which of course doesn't make much sense as the left-hand side is real and the right-hand side is purely imaginary...
In fact, the paper claims that the right answer should be (eq.(14))
$$
\Box \left(\lim_{\epsilon\to0^+} \mathrm{Im}  \left[-(x^0-i\epsilon)^2+(x^i)^2 \right]^{1-n/2}\right) = \frac{(n-2)\pi^{n/2}}{\Gamma(n/2)}\delta(x)\,.
$$
For instance for $n=4$, restricting to negative $x^0$ we correctly get
$$
\lim_{\epsilon\to0^+}\mathrm{Im}\, \frac{1}{-x^2-i\epsilon}=\delta(x^2)\implies \Box \delta(x^2)=2\pi \delta(x)\,.
$$
I think the $i\epsilon$ should arise studying the singularities one has to bypass in the Wick rotation, which are to be defined according to the boundary conditions: for instance the retarded Green's function is given in four dimensions by the formal integral
$$
G_{\text{ret}}(x)=\lim_{\epsilon\to0}\int \frac{d^{4}k}{(2\pi)^4}\frac{e^{-ik^0 x^0+i\mathbf k \cdot \mathbf x}}{-(k^0+i\epsilon)^2+|\mathbf k|^2}=\theta(x^0)\int \frac{d^{3}k}{(2\pi)^{3}}\frac{\sin(|\mathbf k | x^0)}{ |\mathbf k|}e^{i\mathbf k \cdot \mathbf x}\,,
$$ 
so that the support of $G_{\text{ret}}$ lies in $x^0>0$, but I cannot make the connection to the above formula.
Can someone lend me a hand? 
 A: Let us consider the region $\Omega=\{z\in\mathbb C: \mathrm{Re}(z)>0\}\smallsetminus\{z=x^0\in\mathbb R: x^0\ge |\mathbf x|\}$ of the complex $z$ plane, for $|\mathbf x|=r>0$. In this region, we define
$$\boxed{
f_n(z, \mathbf x)=\frac{\Gamma(n/2)}{(2-n)2\pi^{n/2}}(-z^2+|\mathbf x|^2)^{1-n/2}\,.
}
$$
Choosing the branch lines as $\{z=x^0\in\mathbb R: |x^0|\ge |\mathbf x|\}$ we see that $f_n(z)$ is holomorphic in $\Omega$.
Since by construction $f_n(z)$ reduces to the Green's function of the Laplace operator on the imaginary axis, setting $z=ix^0_E + \eta$, as $\eta\to0^+$ we have
$$
[(\partial_0^E)^2+(\partial_i)^2]f_n(ix^0_E, \mathbf x)=\delta(x^0_E, \mathbf x)\,.
$$
We analytically continue the left-hand side equation to the whole region $\Omega$, where it has no singularities and satisfies
$$
[-\partial_z^2+(\partial_i)^2] f_n(z,\mathbf x) = 0\,.
$$
On the other hand, we may now integrate $\varphi(z,\mathbf x)[-\partial_z^2+(\partial_i)^2]f_n(z,\mathbf x)$ along the contour in the figure

yielding
$$
- \lim_{\epsilon\to0^+}\int d^{n-1}x\int_0^{+\infty}\!\!\! \varphi(x^0,\mathbf x)\, \Box \big[ f_n(x^0+i\epsilon,\mathbf x) - f_n(x^0-i\epsilon,\mathbf x) \big]
dx^0
\\
= 
\int d^{n-1}x\int_{-\infty}^{+\infty}\!\!\!\varphi(ix^0_E,\mathbf x) \delta(x^0_E,\mathbf x)\, idx^0_E =i{\varphi(0)}\,.
$$
Thus, 
$$\boxed{
-2\Box\left(\mathrm{Im}\lim_{\epsilon\to0^+}f_n(x^0+i\epsilon)\right)=\delta(x^0,\mathbf x)\,.
}$$
which is what we wanted to prove. Note in particular that this is restricted to $x^0>0$, so no sign ambiguity arises when we expand $(x^0+i\epsilon)^2$. For instance, in dimension $n=4$,
$$
-2f_4(x+i\epsilon)=\frac{1}{2\pi^2}\frac{1}{|\mathbf x|^2-(x^0+i\epsilon)^2}=\frac{1}{2\pi^2}\frac{1}{-x^2-i\epsilon}=\frac{1}{2\pi^2}\mathrm{PV}\frac{1}{-x^2}+\frac{i}{2\pi}\delta(x^2)\,,
$$
as $\epsilon\to0^+$, and we correctly retrieve the retarded Green's function
$$
\theta(x^0)\frac{\delta(x^2)}{2\pi}\,.
$$
