Cutoff-Scheme Renormalization and Order of Integration in QFT The following is the result of Fubini's Theorem, describing when you can replace a double integral with an iterated integral safely:
For a set $X \times Y \subset \mathbb{R}^2$, if $\iint |f(x,y)| d(x,y)$ is a finite number the we have the equality:
$$
\iint_{X\times Y} f(x,y) d(x,y) = \int_X \left[ \int_Y f(x,y) dy \right] dx = \int_Y \left[ \int_X f(x,y) dx \right] dy
$$
So this tells us: We can safely switch the order of integration if integrating the absolute value of the integrand gives a finite number.
If the condition is not met, then the iterated integrals may converge to different values. For example (as given in the wiki) we have the following:
$$
\int_{0}^{1} \left[ \int_{0}^{1} \frac{x^2-y^2}{(x^2+y^2)^2} dx \right] dy = - \frac{\pi}{4}\\
\int_{0}^{1} \left[ \int_{0}^{1} \frac{x^2-y^2}{(x^2+y^2)^2} dy \right] dx = \frac{\pi}{4}
$$
Renormalization: When we encounter a divergent Feynman integral, a renormalization method is to introduce a cutoff so as to characterize 'how' the integral diverges (ie, like a logarithm, or power-law etc.). For example; $\int_{\mathbb{R}^4} \frac{d^4p}{p^2 - m^2}$ is UV divergent, and diverges like $\propto \Lambda^2 - m^2 \log\left( \frac{\Lambda^2}{m^2} \right)$ (for $\Lambda \gg m$ a UV-cutoff).
My Question: What if we have a Feynman integral that diverges differently depending on the order in which you integrate it?
As an example, I have the following contrived integral:
$$
\mathcal{I} := \int_0^\infty dt \int_0^\infty dx\ \frac{\ln(x)^2}{x(1-x^2)(\ln(x)^2+t^2)}
$$
This integral diverges at a few places, but let's just consider how it diverges near $(t,x) = (0,0)$.
CASE I: Integrate $t$ away first, and we are left with the integral $\int_0^\infty \frac{dx}{x(1-x^2)|\ln(x)|}$. This integral diverges at $x=0$ in the following way $\propto \ln\left(\ln(\tfrac{1}{x})\right)$ (for $x \approx 0$) (you can check the asymptotics on this).
CASE II If you integrate $x$ away first (define a new variable $q = \ln(x)$ and manipulate it to get $\int_0^\infty \frac{q^2}{(q^2+t^2)^2} = \frac{\pi}{4|t|}$), we are left with the integral $\int_0^\infty \frac{\pi}{4t}$, which diverges like $\propto \log(t)$ for $t \approx 0$.
So we see that the integral diverges differently depending on the way in which you integrate this.
What if one were to encounter such a situation physically? Is it not allowed for some reason? How do we renormalize things if the way our integral diverges can vary? Furthermore, what if the variable upon which we need to place our cutoff is different depending on the order of integration? This seems very concerning!
Is this situation too contrived, and would not happen in an actual physically relevant QFT problem?
 A: You are asking if perturbation theory can generate conditionally convergent integrals and the answer is yes. The comment about Fubini being safe after the integral is made finite is a bit too fast. Yes if you avoid the regions where $\iint |f(x,y)| \, dx dy$ diverges then $\iint f(x,y) \, dx dy$ can be done in either order. But this is trading an ambiguity in the ordering for an ambiguity in the regulator. In other words, the common statement that we can regulate integrals in any convenient way is false. What we can do is regulate divergences in any convenient way. The same integral might have many divergences with the same physical origin and these need to be regulated consistently.
The example integral I'm thinking of is
$$
I = \int \langle \epsilon(0) \mathcal{O}(x) \mathcal{O}(y) \epsilon(\infty) \rangle \, d^2x d^2y.
$$
This comes up in two loop conformal perturbation theory where we would interpret $\mathcal{O}$ as a classically marginal operator. There are singularities as the operators approach eachother. But if we keep them a distance $a$ apart, the coefficient of the $\log(a)$ divergence of this integral has a nice physical meaning: it is the anomalous dimension of $\epsilon$. (Note that I am using a short distance cutoff here but this is equivalent to a large momentum cutoff if we Fourier transform.) Using the fact that $\mathcal{O}$ is marginal, we can rescale by $|y|$ and substitute $z = x / |y|$ to get
\begin{align}
I &= \int |y|^{-4} \langle \epsilon(0) \mathcal{O}(x / |y|) \mathcal{O}(y / |y|) \epsilon(\infty) \rangle \, d^2x d^2y \\
&= 2\pi \int \frac{d|y|}{|y|} \int |y|^{-2} \langle \epsilon(0) \mathcal{O}(x / |y|) \mathcal{O}(\hat{e}) \epsilon(\infty) \rangle \, d^2y \\
&\sim 2\pi \log(a) \int \langle \epsilon(0) \mathcal{O}(z) \mathcal{O}(\hat{e}) \epsilon(\infty) \rangle \, d^2z.
\end{align}
This integral came up in a study of the long-range Ising model where the 4pt function is
$$
\langle \epsilon(0) \mathcal{O}(z) \mathcal{O}(\hat{e}) \epsilon(\infty) \rangle = \frac{|1 + z|^2}{4|z| |1 - z|^4}.
$$
This has a power law divergence as $z \to 1$ and we don't care about those (remember we're only trying to get the coefficient of the log) but here comes the subtlety. In order to add a counterterm which renormalizes this away, we need to first isolate the divergence with a regulator. And this is already constrained by what we did before. To get the $\log(a)$ above, we pulled out a $2\pi$ with polar co-ordinates which means we had a circular cutoff. Doing this for $y \to 0$ means we also have to do it for $z \to 1$. So now we finally get to
$$
J = \int \frac{1}{|1 - z|^4} \left ( \frac{|1 + z|^2}{4|z|} - 1 \right ) d^2z
$$
which includes the counterterm. If one were to set $z = re^{i\theta}$ and split it up as
$$
\int d^2z \to \int_0^1 \int_0^{2\pi} rdrd\theta + \int_1^\infty \int_0^{2\pi} rdrd\theta
$$
it is relatively straightforward to confirm that the result is $\frac{\pi}{4}$. But   ignoring the measure zero domain $\{ r = 1 \}$ like this means we are approaching the $z = 1$ singularity at different rates depending on the angle thus violating the circular cutoff prescription. So we should be skeptical of this $\frac{\pi}{4}$. Instead, try substituting different polar co-ordinates defined by $1 - z = \rho e^{i\phi}$ and then integrating numerically. The answer might surprise you!
Oh and in case this seems esoteric, results for this anomalous dimension (which favor the $1 - z = \rho e^{i\phi}$ ordering) have been found with Monte Carlo.
A: Usually, when calculating Feynman's integrals, we perform a Wick's rotation so that:
$$\int f(p^2)d^4p\rightarrow\int f(p_E^2)d^4p_E$$
where: $p^2 = p_0^2-\vec{p}^2$, $p_E^2 = p_0^2+\vec{p}^2$, then do the integration in the 4-dimensional Euclidean space by using the 4-dimensional solid angle: $d^4p_E = p_E^3d\Omega_4$, where $p_E$ is the norm in 4-dimensional space having the meaning of energy.
The function in the Feynman integral with $p_E$ is symmetric in term of $p_0$ and $\vec{p}$, so, i think it is not necessary to concern about the Fubini's theorem here.
