Why does the mix of Exclusions and SymbolicProcessing->0 sets my integral to $0$?

Question

I am trying to compute a $5$-dimensional integral, which may or may not be zero. I do believe it is finite though, at least for the values $\epsilon$, $\sigma$ used here. The integral is given in the script below and the integrand depends on the variables $r, \theta_1, \theta_2, \phi, \tau_4$. I set "SymbolicProcessing" -> 0 and Exclusions->r==0, since this point is an indeterminate one for some values of $\epsilon$ and $\sigma$. However, in the case I present here there should be nothing special about this point.

I have observed the following behavior: having both SymbolicProcessing and Exclusions as above results in the integral being evaluated to $0$ in about $5$ seconds. There is no error estimate from the IntegrationMonitor, i.e. the variable errors remains empty! Something is fishy about that. When I comment out the SymbolicProcessing but leave Exclusions, I get $-1.3945964 \cdot 10^7$ in $\sim 100$ seconds. When I comment out Exclusions and keep SymbolicProcessing, I obtain $-9.7620595 \cdot 10^6$ in $\sim 90$ seconds. When I comment both out, the result appears to be $-9.7620595 \cdot 10^6$ in $\sim 80$ seconds.

Surely the result with both options cannot be trusted. If I put Exclusions->r==13 instead of Exclusions->r==0, I also get $0$ in $5$ seconds and no error estimate. The other results seem in phase with one another. So why does the integral behave that way when both options are activated?

My code:

Clear[ϵ, σ, r, θ1, θ2, ϕ, x, y, z, \
τ, τ4, d, x15, x24, x25, x45, R, S, a, f, Φ, \
Y245, integrand, errors]

ϵ = 2;
σ = -2;

x = r*Cos[θ1];
y = r*Sin[θ1]*Sin[θ2]*Cos[ϕ];
z = r*Sin[θ1]*Sin[θ2]*Sin[ϕ];
τ = r*Sin[θ1]*Cos[θ2];

d = Sqrt[x^2 + y^2 + z^2] // FullSimplify;
x15 = (1 - x)^2 + y^2 + z^2 + τ^2 // FullSimplify;
x24 = ϵ^2 + σ^2 + τ4^2 // FullSimplify;
x25 = (ϵ - x)^2 + (σ - y)^2 + z^2 + τ^2 // 
   FullSimplify;
x45 = x^2 + y^2 + z^2 + (τ4 - τ)^2 // FullSimplify;
R = x24/x25;
S = x45/x25;
a = 1/4 Sqrt[4*R*S - (1 - R - S)^2];
f = I Sqrt[-((1 - R - S - 4*I*a)/(1 - R - S + 4*I*a))];
Φ = 
  1/a Im[PolyLog[2, f*Sqrt[R/S]] + 
     Log[Sqrt[R/S]]*Log[1 - f*Sqrt[R/S]]];
Y245 = 1/x25 Φ;
integrand = 1/(d*x15^2) (τ^2/x15 - 1) Y245;

NIntegrate[
  integrand, {r, 0, ∞}, {θ1, 0, π}, {θ2, 
   0, π}, {ϕ, 0, 
   2 π}, {τ4, -∞, ∞}, 
  Exclusions -> {r == 0}, 
  Method -> {"GlobalAdaptive", "SymbolicProcessing" -> 0, 
    "SingularityHandler" -> None}, PrecisionGoal -> 2, 
  AccuracyGoal -> 2, WorkingPrecision -> 20, 
  IntegrationMonitor :> ((errors = Through[#1@"Error"]) &)] // Timing
Length@errors
Total@errors

Answer 1

There is a missing Jacobian $r^3 \sin\theta _1^2\sin \theta _2$ in integrand, so after inserting it in a code we have

NIntegrate[
  integrand r^3 Sin[\[Theta]1]^2 Sin[\[Theta]2], {r, 
   0, \[Infinity]}, {\[Theta]1, 0, \[Pi]}, {\[Theta]2, 
   0, \[Pi]}, {\[Phi], 0, 
   2 \[Pi]}, {\[Tau]4, -\[Infinity], \[Infinity]}, 
  Method -> "AdaptiveMonteCarlo"] // AbsoluteTiming

Out[]= {48.4042, -402.2}