I am looking at the following integral:
$$I= \int_{-\infty}^\infty d\tau_3 \int_{-\infty}^\infty d\tau_4 \frac{1}{(1+\tau_3^2)^2} \Phi \left(\frac{1+\tau_4^2}{1+\tau_3^2},\frac{(\tau_3-\tau_4)^2}{1+\tau_3^2} \right) \tag{1}$$
where $\Phi(r,s)$ is a complicated dimensionless function given in the code at the end of the question. Importantly, the integrand is finite except when $\tau_3 = \tau_4$ (see left plot below). By defining $(\tau_3 - \tau_4)^2 > \epsilon^2$ (point-splitting regularization), we can extract the divergence of $(1)$ to get:
$$\left. I \right|_\text{div} = - \frac{\pi^2}{2} \log \epsilon^2 \tag{2}$$
Now I would like to reproduce that result numerically, i.e. I redefine the integration limits as follows:
$$\left. I \right|_\text{reg} = \left(\int_{-\infty}^\infty d\tau_3 \int_{-\infty}^{\tau_3-\epsilon/2} d\tau_4 + \int_{-\infty}^\infty d\tau_3 \int_{\tau_3+\epsilon/2}^\infty d\tau_4 \right) \frac{1}{1+\tau_3^2} \Phi \left(\frac{1+\tau_4^2}{1+\tau_3^2},\frac{(\tau_3-\tau_4)^2}{1+\tau_3^2} \right) \tag{3}$$
My idea is the following: use NIntegrate to collect numerical data for different values of $\epsilon$ near $0$, then fit the data to a model $I(\epsilon) = a \cdot \log \epsilon^2 + b$. I have imagined this method myself, and thus can not back it up with references. Neither can I say that this should work for sure. However comparing the 3D plot of the integrand of $(1)$ and a plot for my fit function ($a=-\pi^2/2, b=0$) gives me hope:
The problem comes by the practical part. I thought LocalAdaptive would be a good method, but honestly that is just because of the name, I have no clue how it is evaluating the integral really. Anyhow I get the following data for $\epsilon = 0.001, 0.002, ..., 0.010$:
Not only I did not manage to reproduce the $\log$ behavior, but also when I increase WorkingPrecision, I find that NIntegrate fails to converge, which may suggest that the values I find are too inaccurate near the singularity. And if I include $\epsilon=0$ in the dataset, the integral converges, which obviously should not be the case, or at least it should have a large value compared to the rest of the dataset (I get $32.4685$, in the linear continuity of the plot).
So all in all the question is: how should I perform this numerical integration to obtain reliable data, and is this approach worth pursuing in the first place?
Here is my code so far:
x1 = 1;
R[\[Tau]3_, \[Tau]4_] := (x1^2 + \[Tau]4^2)/(x1^2 + \[Tau]3^2);
S[\[Tau]3_, \[Tau]4_] := (\[Tau]3 - \[Tau]4)^2/(x1^2 + \[Tau]3^2);
a[\[Tau]3_, \[Tau]4_] := 1/4 Sqrt[4*R[\[Tau]3, \[Tau]4]*S[\[Tau]3, \[Tau]4] - (1 - R[\[Tau]3, \[Tau]4] - S[\[Tau]3, \[Tau]4])^2];
F[\[Tau]3_, \[Tau]4_] := I Sqrt[-((1 - R[\[Tau]3, \[Tau]4] - S[\[Tau]3, \[Tau]4] - 4 I*a[\[Tau]3, \[Tau]4])/(1 - R[\[Tau]3, \[Tau]4] - S[\[Tau]3, \[Tau]4] + 4 I*a[\[Tau]3, \[Tau]4]))];
Phi[\[Tau]3_, \[Tau]4_] := 1/a[\[Tau]3, \[Tau]4] Im[PolyLog[2, F[\[Tau]3, \[Tau]4] Sqrt[R[\[Tau]3, \[Tau]4]/S[\[Tau]3, \[Tau]4]]] + Log[Sqrt[R[\[Tau]3, \[Tau]4]/S[\[Tau]3, \[Tau]4]]]*Log[1 - F[\[Tau]3, \[Tau]4] Sqrt[R[\[Tau]3, \[Tau]4]/S[\[Tau]3, \[Tau]4]]]];
Integrand[\[Tau]3_, \[Tau]4_] := 1/(x1^2 + \[Tau]3^2)^2 Phi[\[Tau]3, \[Tau]4];
data = Table[{\[Epsilon], NIntegrate[Integrand[\[Tau]3, \[Tau]4], {\[Tau]3, -\[Infinity], \\[Infinity]}, {\[Tau]4, -\[Infinity], \[Tau]3 - \[Epsilon]/2},Method -> "LocalAdaptive"] + NIntegrate[Integrand[\[Tau]3, \[Tau]4], {\[Tau]3, -\[Infinity], \\[Infinity]}, {\[Tau]4, \[Tau]3 + \[Epsilon]/2, \[Infinity]}, Method -> "LocalAdaptive"]}, {\[Epsilon], 0.001, 0.01, 0.001}];
ListPlot[data]
