I finally found what seems to be the answer. QMechanic's answer solidifies this as paradox - the Euclidean and Minkowski Feynman integrals "appear" to give the exact same result as per his answer, which we naively expect, but at the same time there is an unaccounted pole contribution as per my original post which explicitly does not vanish, and this pole contribution precisely represents the difference between the Minkowski and Euclidean integrals. So they're equal... but not equal? What's going on? Only one answer can be correct. Spoiler alert: they are not equal.
This paradox was brought up nearly verbatim in a recent-ish paper$^1$ by Carlson & Freid (2017). [1] The authors found this tension in the context of calculating a certain loop-correction in both Euclidean and Minkowski signatures. There, the unaccounted pole-contribution (which this post is all about) gave rise to a crucial IR-divergence, which suggested that calculating the correlator in Euclidean signature may be invalid, or at least require greater care. What the author didn't seem to point out was that this issue is extremely general and would affect nearly every Euclidean loop calculation.
Just a couple months later, the issue was resolved in a paper by R. Briceño et al (2017). [2] Section 3 of this paper fully explains and resolves this issue. The upshot is the following:
We need to be sure we are calculating a quantity which will actually be the same in both Euclidean and Minkowski signatures.
We need to be careful with how we go about calculating this quantity in both signatures.
I'll now explain the relevant points given in both of these papers.
Section 1: $I_M\neq I_E$
Let's consider a simple theory of two scalars $\phi,\chi$ in $d=2$ dimensions. Let each scalar have non-zero pole-mass $m_\phi,m_\chi$, and let the interaction simply be:
$$S_{\textrm{int}}=\int d^2x\, \frac{g}{2!} \phi^2 (x)\chi(x) \tag{1}$$
Consider the one-loop contribution to the two-point function coming from
a $\phi-\chi$ loop, as illustrated below.
In Minkowski (mostly plus) signature, this (amputated) contribution is given by:
$$I_M(p^2)=g^2 \int \frac{d^2 k}{(2\pi)^2}\frac{1}{\left(k^2+m_{\chi}^2-i\epsilon \right)\left((k-p)^2+m_{\phi}^2-i\epsilon \right)} \tag{2}$$
where
$$\begin{align}
k_\mu &= \left(k_0,k_1\right) & p_\mu &= \left(\sqrt{m_{\phi}^2+p_1^2}, p_1\right) \\
k^2&=-k_0^2+k_1^2 & p^2&=-m_\phi^2
\end{align} \tag{3}$$
We can solve this integral using Feynman parameterization, but for transparency I'll evaluate each integral directly. Let's first fix $k_1$ and do the $k_0$ integral. In the complex $k_0$ plane we have 4 poles, just as in my original post above. They are
$$k_0^{\pm} = \pm\left(\omega_{k}^\chi - i\epsilon\right) \tag{4}$$
$$\tilde{k_0}^{\pm} = \omega_p^\phi \pm\left(\omega_{k-p}^\chi - i\epsilon\right) \tag{5}$$
where I have used the shorthand notation $\omega_q^{\phi,\chi}\equiv \sqrt{m_{\phi,\chi}^2+q_1^2}$. We can close the contour in the UHP and get:
$$I_M(p^2)=\frac{-ig^2}{4\pi}\int_{-\infty}^{\infty}dk_1 \left[\frac{1}{\omega_{k}^{\chi}\left(\left(\omega_{k}^{\chi}+\omega_{p}^{\phi}\right)^2-\left(\omega_{k-p}^{\phi}\right)^2\right)}+\frac{1}{\omega_{k-p}^{\phi}\left(\left(\omega_{p}^{\phi}-\omega_{k-p}^{\phi}\right)^2-\left(\omega_{k}^{\chi}\right)^2\right)}\right] \tag{6}$$
A plot of this function is shown below. As required, the result is only a function of the invariant masses $m_{\phi,\chi}$, i.e. it is independent of $p_1$. (below I set $m_\chi = 1$)
Naively, the same calculation in Euclidean signature would lead to the "Wick-rotated" integral $I_E$, with the $k_0$ contour running over the imaginary axis rather than the real axis.
$$I_E(p^2)\overset{?}{=} i g^2 \int \frac{d^2 k}{(2\pi)^2}\frac{1}{\left(k^2+m_{\chi}^2-i\epsilon \right)\left((k-p)^2+m_{\phi}^2-i\epsilon \right)} \tag{7}$$
where in the above expression all square-momenta are Euclidean, so $p_0=ip_2=i\omega_p^\phi$. However, we must remember that in $I_M$, for certain values of $k_1$ the pole $\tilde{k_0}^-$ will reside in the first quadrant.
$$(k_1-p_1)^2 < p_1^2 \implies \textrm{Re}\left( \tilde{k_0}^-\right) >0 \tag{8}$$
In trying to smoothly deform (rotate CCW) the Minkowskian contour to the "Euclidean" one, we will encircle this pole in the first quadrant. The contribution from this pole will therefore give the difference between $I_M$ and $I_E$. The result is:
$$I_M-I_E\equiv \Delta I =\frac{-ig^2}{4\pi}\int_{-p_1}^{p_1} dk_1 \frac{1}{\omega_{k}^{\phi}\left(\left(\omega_{p}^{\phi}-\omega_{k}^{\phi}\right)^2-\left(\omega_{k+p}^{\chi}\right)^2\right)} \tag{9}$$
A plot of this difference is shown below. Again, $m_\chi=1$.
Note that it is not zero! Not only that, but it depends explicitly on $p_1$. Clearly, $I_E$ does not represent a physically relevant quantity.
But what about the analytical continuation, as in QMechanic's answer? Can we not regard $I_M(p_0)$ as analytic in $p_0$, and then analytically continue to the imaginary axis $p_0\rightarrow i p_2$? Can we not evaluate $I_M(p_0)$ and $I_E(p_0)$ in their respective ranges of validity, and then claim that since they agree over a certain domain in the complex $p_0$ plane, by analytic continuation they must agree everywhere?
The answer to both of these is a resounding no, as evidenced by the previous direct numerical calculations. I have directly calculated $I_M$ and $I_E$ at the physical point $p_0=\omega_{p}^\phi$, and they are finite but not equal.
The expression I calculated for $I_M(p_0)$ is valid for all real $p_0$. We could have equally calculated via Feynman parameterization, and obtained $I_M(p_0)$ for real $p_0$. However if we did the same calculation via Feynman parameterization for $I_E(p_0)$, it would only be valid for imaginary $p_0$, as in QMechanic's answer. The intersection of these two domains of validity is the single point $p_0=ip_2=0$. Analytic continuation from one domain to another would require agreement of the two functions over at least an accumulation point in the complex $p_0$ plane. We do not have that, and therefore we cannot justify analytic continuation.
Section 2: Resolving the paradox
Let's consider the quantum theory in the Heisenberg picture. From a Hamiltonian perspective, assuming a time-independent Hamiltonian for scalar particles, the only difference between the Euclidean and Minkowski theories is through the time evolution of operators.
$$\begin{align}
\mathcal{O}(t)&=e^{iHt}\,\mathcal{O}(0)\,e^{-iHt}, & \mathcal{O}(\tau)&=e^{H\tau}\,\mathcal{O}(0)\,e^{-H\tau}
\end{align} \tag{10}$$
where $t=-i\tau$. The instantaneous Hamiltonian is the same in both signatures, ergo so is the spectrum. If we thus ask for a correlator which has no explicit time dependence, the result should be exactly the same in both theories. For example, in the two-scalar theory from earlier, $\phi,\phi$ overlap:
$$\langle P,\phi|P,\phi\rangle \tag{11}$$
should be the same in both signatures, as there is no time-dependence whatsoever involved. Of course, this matrix element is sort of trivial, since by assumption it gives a delta-function, but this is the matrix element which we use to renormalize the kinetic Lagrangian (by demanding triviality).
In Minkowski signature, we calculate overlap amplitudes such as (11) via the LSZ prescription.
$$\langle P,\phi|P,\phi\rangle = \lim_{P^2\rightarrow -m_{\phi}^2} \frac{P^2+m_{\phi}^2}{iZ_{\phi}}\frac{P^2+m_{\phi}^2}{iZ^*_{\phi}} \langle \phi(P)\phi(-P)\rangle\tag{12}$$
where
$$\langle \phi(P)\phi(-P)\rangle = \int d^4 y \,e^{-iP_\mu y^\mu} \int d^4 x \,e^{iP_\mu x^\mu} \langle 0 |\textrm{T}\left\{ \phi(y)\phi(x)\right\} |0\rangle \tag{13a}$$
$$Z_\phi = \langle 0 |\phi(0)|\vec p, \phi\rangle \tag{13b}$$
Note that it is for (12) that we use standard momentum-space Feynman rules. It is precisely with this that, in Minkowski signature, we get $I_M$ as in (2), after amputating external propagators of course via the LSZ prescription (11). More importantly, if we calculated this quantity in Euclidean signature, we would indeed get $I_E$ as in (7).
However LSZ is not applicable in Euclidean signature - in fact it doesn't even make sense. The desired matrix element can be accessed through a different prescription. Consider the following Euclidean-time dependent propagator.
$$C(\tau',\tau,\vec P)=\langle \phi(\vec p,\tau' ) \phi (\vec p,\tau )\rangle = \int d^3 y\,e^{-i\vec P \cdot \vec y} \,\int d^3 x\,e^{i\vec P \cdot \vec x} \, \langle \phi(y,\tau' ) \phi (x,\tau )\rangle \tag{14}$$
On the one hand, we can use the delta-function identity
$$f(\tau)=\int \frac{dP_4}{2\pi}\,e^{\pm i P_4 \tau}\,\int dx_4\,e^{\mp iP_4 x_4} f(x_4) \tag{15}$$
to rewrite (14) in terms of something resembling $I_E$.
$$C(\tau',\tau,\vec P)=\int \frac{dP_4'}{2\pi}\,e^{iP_4'\tau'} \int \frac{dP_4}{2\pi}\,e^{-iP_4\tau} \langle \phi( p' ) \phi (p)\rangle \tag{16}$$
where
$$ \langle \phi( p' ) \phi (p)\rangle=\int d^4 y\,e^{-iP \cdot y} \,\int d^4 x\,e^{iP \cdot x} \, \langle \phi(y ) \phi (x)\rangle \tag{17}$$
with $P=(\vec P, P_4)$ and $P'=(\vec P, P_4')$ off-shell momenta. The dot-products above represent euclidean contractions $a\cdot b = \sum a_i b_i$. Again, note that if we were to set $P_4=P_4'=i\omega_p^{\phi}$ on-shell in (17), we would essentially calculate $I_E$ as in (7).
On the other hand, we may insert the identity operator twice into (14), and extract the Euclidean time dependence via (10).
$$C(\tau',\tau,\vec P)=ZZ^* e^{-\omega_p (\tau' - \tau)} \langle \vec P,\phi|\vec P,\phi\rangle + O\left(e^{-E'\tau'+E\tau}\right) \tag{18}$$
where $E,E'>\omega_p$ are multi-particle-state energies which are necessarily greater than that of the relevant single-particle state$^2$. To extract the desired matrix element, we pick out the leading exponential dependence. Putting this together with (16) gives us:
$$ \langle \vec P,\phi|\vec P,\phi\rangle = \lim_{\tau'\rightarrow \infty \\ \tau\rightarrow -\infty} \frac{1}{Z_\phi Z_\phi^*} e^{\omega_p (\tau' - \tau)}\int \frac{dP_4'}{2\pi}\,e^{iP_4'\tau'}\,\int \frac{dP_4}{2\pi}\,e^{-iP_4\tau}\,\langle \phi( p' ) \phi (p)\rangle \tag{19}$$
This is the Euclidean prescription, which is to be compared with the Minkowski LSZ prescription (12). For a particular Feynman diagram, the prescription is as follows:
Calculate the full, unamputated diagram with off-shell $P_4$ for the incoming state and $P_4'$ for the outgoing state.
Integrate over $P_4$ and $P_4'$, with the exponential factors as in (12).
Pick out the leading $\tau,\tau'$ dependence, i.e. the term which scales as $e^{-\omega_{\textrm{out}}\tau'+\omega_{\textrm{in}}\tau}$.
If one carries out this procedure for the Euclidean diagram in section 1 which contributed to $\langle \phi(P)\phi(-P)\rangle$, one will find two terms. The first will be $I_E$, and the second will exactly be $\Delta I$ as found in (9). I will not carry out this explicit calculation, but you can see it in section 3 of [2].
And thus, the paradox is resolved. :)
Footnotes.
$^1$ Given the generality of this issue, I wouldn't be surprised if this issue has already appeared countless times in the literature, unbeknownst to each new discoverer.
$^2$ This point requires more care in a theory with massless particles.
