Using polar coordinates r and f, the region of integration is given by
{ 0<= r <=2/Cos[f], 0<= f <= 2 \[Pi] }
We can then proceed as follows.
First integration with PrincipalValue
g = 8 Integrate[r/(1 - r^2), {r, 0, 2/Cos[f]}, Assumptions -> 0 < f < \[Pi]/4, PrincipalValue -> True]
(* Out[451]= -4 Log[-1 + 4 Sec[f]^2] *)
Second integration
h = Integrate[g, {f, 0, \[Pi]/4}];
FullSimplify[h]
(* Out[453]= 4 Catalan + 1/2 \[Pi] Log[97 - 56 Sqrt[3]] -
2 I (PolyLog[2, 7 I - 4 I Sqrt[3]] - PolyLog[2, I (-7 + 4 Sqrt[3])]) *)
% // N
(* Out[454]= -4.32381 + 0. I *)
This is a real quantity, and it is in agreement with the real part of one of the values mentioned in the OP.
EDIT #1
Your edit provides a mathematically more natural problem since it avoids the hard ("man-made") cutoff.
Your function to be integrated of the whole space is
f = Exp[-(p1px^2 + p1py^2 + p1pz^2 + p1x^2 + p1y^2 + p1z^2) ] 1/(-2 + p1px^2 + p1py^2 + p1pz^2 + (Sqrt[3]/2 + p1x)^2 + p1y^2 + (1/2 + p1z)^2);
Simplifying the name of your variables a bit gives
f1 = f /. {p1px -> x, p1py -> y, p1pz -> z, p1x -> u, p1y -> v, p1z -> w}
(*
Out[6]= E^(-u^2 - v^2 - w^2 - x^2 - y^2 - z^2)/(-2 + (Sqrt[3]/2 + u)^2 + v^2 + (1/2 + w)^2 + x^2 + y^2 + z^2)
*)
Now let us produce the denominator using
Integrate[Exp[-a h], {a, 0, \[Infinity]}, Assumptions -> h > 0]
(*
Out[97]= 1/h
*)
i.e.
Integrate[Exp[-a h], {a, 0, \[Infinity]}, Assumptions -> h > 0] /.
h -> p + (Sqrt[3]/2 + u)^2 + v^2 + (1/2 + w)^2 + x^2 + y^2 + z^2
(*
Out[99]= 1/(p + (Sqrt[3]/2 + u)^2 + v^2 + (1/2 + w)^2 + x^2 + y^2 + z^2)
*)
Here we have replaced the number (-2) temporarily by the parameter p
.
Now, under the a
-intergal we have all terms in the exponent:
f2 = Exp[-u^2 - v^2 - w^2 - x^2 - y^2 - z^2 -
a (p + (Sqrt[3]/2 + u)^2 + v^2 + (1/2 + w)^2 + x^2 + y^2 + z^2)];
and the integral over the whole space is easily done
f3 = Integrate[
f2, {x, -\[Infinity], \[Infinity]}, {y, -\[Infinity], \[Infinity]}, {z, -\
\[Infinity], \[Infinity]}, {u, -\[Infinity], \[Infinity]}, {v, -\[Infinity], \
\[Infinity]}, {w, -\[Infinity], \[Infinity]}, Assumptions -> a > 0]
(*
Out[104]= (E^(-a (1/(1 + a) + p)) \[Pi]^3)/(1 + a)^3
*)
The expression we are looking for is therefore given by
g[p_] = Integrate[f3, {a, 0, \[Infinity]}, Assumptions -> p > 0]
(*
Out[107]= Integrate[(E^(-a (1/(1 + a) + p)) \[Pi]^3)/(1 + a)^3, {a, 0, \[Infinity]}, Assumptions -> p > 0]
*)
In order for the integral to converge we have tried to get an expression with p>0
hoping to then be able to continue it analytically.
But the integral is not done analytically by Mathematica. Hence we try a numeric approach
gn[p_] := NIntegrate[(
E^(-a (1/(1 + a) + p)) \[Pi]^3)/(1 + a)^3, {a, 0, \[Infinity]}]
For example
gn[1]
(* Out[112]= 7.54936 *)
Now we would need information about a kind of mapping from this divergent integral if p<0
and the procedure of taking the PrincipalValue.
I'll leave the results of my preliminary studies of this question for a later edit.
Sumarizing
We have reduced the multidimensional integral to a one dimensional integral.
It remains to be studied how taking the principal value in the original integral maps to a procedure to give a meaning the divergent one-dimensional integral.
Remark: The reduction to a one dimensional integral described here can of course also be done for the "cut-off" version of the problem.
EDIT #2
The analytic continuation to p < 0
can be done as follows.
After the transformation a -> 1/z-1
the integral g
can be written as
g1[p_] = \[Pi]^3 E^(p - 1) Integrate[z E^z Exp[-p/z], {z, 0, 1}];
Expanding E^z
this becomes
g2[p_] = \[Pi]^3 E^(p - 1)
Sum[1/k! Integrate[z^(k + 1) Exp[-p/z], {z, 0, 1}], {k, 0, \[Infinity]}];
But
Integrate[z^(k + 1) Exp[-p/z], {z, 0, 1}, Assumptions -> p > 0]
(*
Out[212]= ExpIntegralE[3 + k, p]
*)
so that
g3[p_] = \[Pi]^3 E^(p - 1)
Sum[1/k! ExpIntegralE[3 + k, p], {k, 0, \[Infinity]}];
But now, this expression can be continued analytically to p < 0
.
And we obtain numerically to a very good approximation
Table[{k, \[Pi]^3 (E^(-1 + p) ExpIntegralE[3 + k, p])/k! /. p -> -2.}, {k, 0,
15}]
(*
Out[213]= {{0, 1.81404 - 9.69943 I}, {1, 5.01155 - 6.46628 I}, {2,
2.67871 - 1.61657 I}, {3, 0.73738 - 0.215543 I}, {4,
0.140661 - 0.0179619 I}, {5, 0.021617 - 0.00102639 I}, {6,
0.00288102 - 0.0000427664 I}, {7, 0.000342928 - 1.35766*10^-6 I}, {8,
0.0000368633 - 3.39416*10^-8 I}, {9, 3.6023*10^-6 - 6.85689*10^-10 I}, {10,
3.21984*10^-7 - 1.14282*10^-11 I}, {11,
2.64847*10^-8 - 1.59834*10^-13 I}, {12,
2.01624*10^-9 - 1.90279*10^-15 I}, {13,
1.42798*10^-10 - 1.95158*10^-17 I}, {14,
9.4526*10^-12 - 1.74248*10^-19 I}, {15, 5.87243*10^-13 - 1.36665*10^-21 I}}
*)
Summing up
Plus @@ %
(*
Out[214]= {120, 10.4072 - 18.0169 I}
*)
Now we make the bold assumptions that taking the principal value is equivalent to taking the real part here.
Then the original integral becomes numerically
gg = 10.407221182257118 ...
Of course we need to justify the bold assumption. I'll try to do this later
EDIT #3
Gaussian decay instead of cut off
Let us first go back to the original toy example, but now with Gaussian decay instead of hard cut off.
The integral is now
f2g = Integrate[Exp[-x^2 - y^2]/(
1 - x^2 - y^2), {x, -\[Infinity], \[Infinity]}, {y, -\[Infinity], \
\[Infinity]}]
During evaluation of In[286]:= Integrate::idiv: Integral of
E^(-x^2-y^2)/(-1+x^2+y^2) does not converge on
{-[Infinity],[Infinity]}. >>
Using PrincipalValue gives a finite result
f2gp = Integrate[Exp[-x^2 - y^2]/(
1 - x^2 - y^2), {x, -\[Infinity], \[Infinity]}, {y, -\[Infinity], \
\[Infinity]}, PrincipalValue -> True]
(*
Out[287]= (\[Pi] ExpIntegralEi[1])/E
*)
% // N
(* Out[288]= 2.19024 *)
We can also transform to polar coordinates.
The angular part is now trivial and gives a factor 2 [Pi]. The integral is therefore
2 \[Pi] Integrate[(Exp[-r^2] r)/(1 - r^2), {r, 0, \[Infinity]},
PrincipalValue -> True]
(*
Out[289]= (\[Pi] ExpIntegralEi[1])/E
*)
in agreement with the integral in Cartesian coordinates.
Notice that, in contrast to the cut off version, the values is positive.
Now let's do the same integral with replacement of the denominator (as in EDIT #1)
Integrate[Exp[-x^2 - y^2 -
a (p - x^2 -
y^2)], {x, -\[Infinity], \[Infinity]}, {y, -\[Infinity], \[Infinity]}]
(*
Out[322]= ConditionalExpression[-((E^(-a p) \[Pi])/(-1 + a)), Re[a] < 1]
*)
Integrate[%[[1]], {a, 0, \[Infinity]}, Assumptions -> p > 0,
PrincipalValue -> True]
(*
Out[323]= E^-p \[Pi] ExpIntegralEi[p]
*)
% /. p -> 1.
(* Out[324]= 2.19024 *)
Higher dimensions with standard integrand
3 dimensions:
f3 = Integrate[Exp[-x^2 - y^2 - z^2]/(
p - x^2 - y^2 -
z^2), {x, -\[Infinity], \[Infinity]}, {y, -\[Infinity], \[Infinity]}, {z, \
-\[Infinity], \[Infinity]}, PrincipalValue -> True, Assumptions -> p > 0]
(*
Out[310]= 2 \[Pi]^(3/2) (-1 + 2 Sqrt[p] DawsonF[Sqrt[p]])
*)
4 dimensions
f4 = Integrate[Exp[-x^2 - y^2 - z^2 - t^2]/(
p - x^2 - y^2 - z^2 -
t^2), {x, -\[Infinity], \[Infinity]}, {y, -\[Infinity], \[Infinity]}, {z, \
-\[Infinity], \[Infinity]}, {t, -\[Infinity], \[Infinity]},
PrincipalValue -> True, Assumptions -> p > 0]
(*
Out[336]= \[Pi]^2 (-1 + E^-p p ExpIntegralEi[p])
*)
Arbitrary dimension
ff[n_] := Integrate[Exp[-Sum[x[i]^2, {i, 1, n}]]/(p - Sum[x[i]^2, {i, 1, n}]),
Sequence @@ Table[{x[i], 0, \[Infinity]}, {i, 1, n}],
PrincipalValue -> True, Assumptions -> p > 0]
t = Table[ff[n], {n, 1, 6}]
(*
Out[354]= {
(Sqrt[\[Pi]] DawsonF[Sqrt[p]])/Sqrt[p],
1/4 E^-p \[Pi] ExpIntegralEi[p],
1/4 \[Pi]^(3/2) (-1 + 2 Sqrt[p] DawsonF[Sqrt[p]]),
1/16 \[Pi]^2 (-1 + E^-p p ExpIntegralEi[p]),
1/48 \[Pi]^(5/2) (-1 - 2 p + 4 p^(3/2) DawsonF[Sqrt[p]]),
1/128 \[Pi]^3 (-1 - p + E^-p p^2 ExpIntegralEi[p])}
*)
Plot[{t[[1]], t[[3]], t[[5]]}, {p, -2, 3}]
(* 150902_plot_135.jpg *)

Plot[{t[[2]], t[[4]], t[[6]]}, {p, -2, 3}]
(* 150902_plot_246.jpg *)

EDIT #4
We can get a confirmation of the my bold hypothesis by considering another type of regularization of the Gaussian toy example.
$$\text{toyg}=\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }\frac{e^{-x^2-y^2}}{-x^2-y^2+1}dydx$$
This method considers a general power k of the denominator and applies analytic continution of the result to k->-1
Integrate[Exp[-x^2 - y^2] (1 - x^2 - y^2)^
k, {x, -\[Infinity], \[Infinity]}, {y, -\[Infinity], \[Infinity]}]
(*
Out[136]= ConditionalExpression[
E^(-1 - I k \[Pi]) \[Pi] ((-1 + E^(2 I k \[Pi])) Gamma[1 + k] +
E Subfactorial[k]), Re[k] > -1]
*)
f = First[%]
(*
Out[137]= E^(-1 - I k \[Pi]) \[Pi] ((-1 + E^(2 I k \[Pi])) Gamma[1 + k] +
E Subfactorial[k])
*)
Hence toyg should be
tyogk = Limit[f, k -> -1]
(*
Out[138]= (\[Pi] (-2 I \[Pi] - Gamma[0, -1]))/E
*)
% // N
(*
Out[139]= 2.19024 - 3.63082 I
*)
We have obtained a complex quantity.
Comparing this with the PV form
toygpv = Integrate[Exp[-x^2 - y^2]/(1 - x^2 -
y^2), {x, -\[Infinity], \[Infinity]}, {y, -\[Infinity], \[Infinity]},
PrincipalValue -> True]
(*
Out[134]= (\[Pi] ExpIntegralEi[1])/E
*)
% // N
(* Out[135]= 2.19024 *)
shows that the real parts conincide.