I am doing some physics calculations, and I need to deal with integrals that can be reduced to this general form: $$ \lim_{\epsilon \to 0^{+}} \int_{-a}^{b} \int_{-a}^{b} \frac{f(x,y)}{(x-i\epsilon)(y-i\epsilon)}dxdy, $$ where in my case $f(x,y)$ is symmetrical with respect to interchanging $x$ and $y$, and I also only need the real part of these integrals. I am not sure how the Sokhotski–Plemelj theorem would apply in such cases.
I tried to experiment numerically with the following expression instead: $$ \lim_{\epsilon_1 \to 0^{+}} \lim_{\epsilon_2 \to 0^{+}} \int_{-a}^{b} \int_{-a}^{b} \frac{f(x,y)}{(x-i\epsilon_1)(y-i\epsilon_2)}dxdy, $$ since I was not completely sure if $\epsilon_1$ and $\epsilon_2$ had to be the same in my formulas. In the latter case, I could imagine applying the Sokhotski–Plemelj theorem twice, so I would obtain the following: $$\begin{align} Re \left( \lim_{\epsilon_1 \to 0^{+}} \lim_{\epsilon_2 \to 0^{+}} \int_{-a}^{b} \int_{-a}^{b} \frac{f(x,y)}{(x-i\epsilon_1)(y-i\epsilon_2)}dxdy \right) & = \mathrm{PV}\int_{-a}^{b} \frac{1}{x} \left[ \mathrm{PV}\int_{-a}^{b} \frac{f(x,y)}{y} dy \right] dx \\ & - \pi^2f(0,0), \end{align}$$ but the corresponding numerical results seemed to be wrong from the physics point of view.
It is actually much more likely that I need to work with the first expression, and I suspect there are some subtleties that I am missing. How do I take the limit in the first equation properly when $\epsilon_1$ and $\epsilon_2$ coincide?
Also, if I consider a simple case of $f(x,y)=1$ and $a=b=1$, the answer seems to be $-\pi^2$, which I find rather strange, so I realize that I do not have a good intuition for this kind of things.
Any thoughts/ideas would be greatly appreciated!
Prelude:
The product of two principal value integrals (denoted by $\mathcal{P}$ in the following) is special due to the fact that in a double integral the order of integration matters. In particular, the Poincare-Bertrand theorem states that $$\int\!dx\,\frac{\mathcal{P}}{x-u}\int\!dy\,\frac{\mathcal{P}}{y-x} f(x,y) = \int\!dy\,\int\!dx\,\frac{\mathcal{P}}{x-u}\frac{\mathcal{P}}{y-x} f(x,y) - \pi^2 f(u,u)\,; \tag{1}$$ for a proof see, e.g., here. The Poincare-Bertrand theorem is important, as it shows that the Hilbert transform $$ \mathcal{H}(f)(x) = \frac{1}{\pi}\int\!\frac{\mathcal{P}}{y-x}f(y)$$ is an involution. Indeed $$ \mathcal{H}(\mathcal{H}(f))(u) =\frac{1}{\pi^2} \int\!dx\,\frac{\mathcal{P}}{u-x}\int\!\frac{\mathcal{P}}{y-x}f(y) = f(u)$$ due to (1) and the fact that $$\int\!dx\,\frac{\mathcal{P}}{u-x}\frac{\mathcal{P}}{y-x} = 0\,.$$
Main:
In your specific case, it is not possible to simply use the Sokhotski–Plemelj theorem for a double integral due to the fact that the double principal value integral depends on the order of integration (see above). In particular, the Sokhotski–Plemelj theorem for double integrals that takes care of the Poincare-Bertrand theorem is $$ \lim_{\epsilon_1, \epsilon_2 \to 0^+} \frac{1}{x-i \epsilon_1}\frac{1}{y-i \epsilon_2} = \left[ \frac{\mathcal{P}}{x} +i \pi \delta(x) \right] \left[ \frac{\mathcal{P}}{y} +i \pi \delta(y) \right] + \pi^2 \delta(x)\delta(y)\,. $$ Note that the last term, due to the Poincare-Bertrand theorem, cancels the product of $\delta$-distributions.