This is a follow-up question to a post I made on Physics SE. TLDR; Essentially I have derived an apparent contradiction, starting from a simple integral. It feels like I'm missing something quite elementary and trivial, but I cannot seem to resolve it.
Consider the following integral, whose approximate form is commonly finds in physics.
$$f(x)=\int_{-\infty}^{+\infty}\frac{dz}{\left(-z^2 + m^2 - i\epsilon\right )^2 \left(-(z-x)^2+ M^2 -i \epsilon\right)} \tag{1}$$
$m,M>0$ are positive real numbers, $x$ is a positive real number (which we may allow to be complex), and $\epsilon>0$ is an infinitesimal small number which shifts the poles away from the real axis. In the complex $z$ plane, the integral looks as follows:
The blue dots correspond to the double-poles at $z=\pm (m - i\epsilon)$, and red dots correspond to the single-poles at $z=x\pm (M-i\epsilon )$. The green line represents the integration contour. Notice that the pole at $z=x-M+i\epsilon$ resides in the $1^{\textrm{st}}$ quadrant.
Now let's apply the technique of Feynman parameterization, which combines both the rational terms in the integrand into a single rational function.
$$\begin{align} f(x)&=\int_{-\infty}^{+\infty}dz \int_0^1 du_1 \int_0^{1-u_1} du_2 \frac{2!}{\left[ (1-u_1-u_2)\left( -(z-x)^2+ M^2 -i \epsilon \right)+u_2 \left( -z^2 + m^2 - i\epsilon\right) +u_3 \left( -z^2 + m^2 - i\epsilon \right) \right]^3} \\ &= \int_0^1 du_1 \int_0^{1-u_1} du_2 \int_{-\infty}^{+\infty}dz \frac{2}{\left[-(z-(1-u_1-u_2)x)^2-(u_1+u_2)(1-u_1-u_2)x^2+(1-u_1-u_2)M^2+(u_1+u_2)m^2-i\epsilon\right]^3} \\ &= \int_0^1 du_1 \int_0^{1-u_1} du_2 \int_{-\infty}^{+\infty}dz \frac{2}{\left[-z^2+g(u_1,u_2,x)^2-i\epsilon\right]^3} \tag{2} \end{align}$$
In the final expression above, note that no matter if $g^2$ is positive or negative, the triple-poles will always lie in the $2^{\textrm{nd}}$ and $4^{\textrm{th}}$ quadrants.
Since no poles $1^{\textrm{st}}$ and $3^{\textrm{rd}}$ quadrants, we may rotate the contour counter-clockwise by $90^\circ$, also known as Wick-Rotation. We now have the following expression:
$$f(x)=\int_0^1 du_1 \int_0^{1-u_1} du_2 \int_{-\infty}^{+\infty} i dz \, \frac{1}{\left[z^2+g(u_1,u_2,x)^2 -i\epsilon\right]^3} \tag{3}$$
Let's now reverse all the previous algebraic steps, but using the Wick-rotated expression. This will finally give us (compare with (1)):
$$f(x)=\int_{-\infty}^{+\infty}\frac{d(iz)}{\left(z^2 + m^2 -i\epsilon\right )^2 \left(-(iz-x)^2+ M^2 -i\epsilon\right)} \tag{4}$$
But this is none other than the Wick-rotated version of the original integral! This is explicitly telling us that:
$$\int_{\Gamma_1} I(z) dz= \int_{\Gamma_2} I(z) dz \tag{5}$$
where $\Gamma_1$ is the contour along the real axis, $\Gamma_2$ is the contour along the imaginary axis, and $I(z)$ is the integrand in (1). Note that this is expression was derived by @QMechanic in the PSE post mentioned previously. The essence of this derivation is found in most introductory textbooks on Quantum Field Theory.
On the other hand, if we tried to smoothly rotate the original contour by $90^\circ$, we would unavoidably pick up the residue of the pole in the $1^{\textrm{st}}$ contour.
$$\int_{\Gamma_1} I(z) dz= 2\pi i \,\textrm{Res} \left(I(z);z=x-M-i\epsilon \right)+ \int_{\Gamma_2} I(z) dz \tag{6}$$
Below you can see an illustration of this intuitive equation.
Combining (5) and (6), we see that the residue must be zero. But we can calculate this residue and show it isn't zero, a contradiction!
$$\textrm{Res} \left(I(z);z=x-M-i\epsilon \right)=\frac{1}{2M\left(-(x-M)^2+m^2\right)}\overset{?}{\neq} 0 \tag{7}$$
So what's going on here? How can I resolve this contradiction?
Edit 1: A possible resolution.
I've just realized an issue in the algebraic steps. In (2), as per @Maxim's comment, I shifted $z\rightarrow z+(1-u_1-u_2)x$ in order to remove the term linear in $z$. This was fine, since $x$ was real-valued.
However when reversing these steps later, going from (3) to (4), the equivalent shift would be $z\rightarrow z+i(1-u_2-u_2)x$. This is not valid since for arbitrary $x$ the shift may hop over one of the poles! The resulting shifted-contour would not be homotopic to the original one! For it to be valid we would need $x$ to be imaginary.
If we have analytic continuation in mind, we could simply assume at this step that $x$ is actually purely imaginary. But then our original shift $z\rightarrow z+(1-u_1-u_2)x$ would not be valid for the same reason!
I found the solution to this. The integrals along the real and imaginary $z$ axes are of course not equal, the difference being given by the finite pole contribution as mentioned in (6) and (7) of my OP. The logic given in my edit, about why it is not possible to reverse the Feynman parameterization (i.e. going from (3) to (4)), is correct.
The broader issue motivating this post is that in physics, we would naively expect (1) and (4) to be the same. One arrives at integrals similar to (1) when calculating Feynman diagrams in Minkowski signature, whereas one would naively expect the same calculation in Euclidean signature to give rise to the "Wick-rotated" form of the same integral (e.g. (4) compared to (1)). For most 'normal' quantum field theories we expect the Minkowski and Euclidean predictions to coincide. Therefore the fact that these two integrals don't agree seems to be a big issue. However this issue can be resolved by going back to the drawing board, and seeing how these kinds of integrals really arise in the first place.
All this is explained in my PSE answer here.