Calculating the value of a triple integral

At first we study the denominater of the integrand. \begin{gather*} |(x,y,z-1)|^3|(x,y,z+1)|^3=(x^2+y^2+z^2+1-2z)^{3/2}(x^2+y^2+z^2+1+2z)^{3/2} =\\[2ex]((x^2+y^2+z^2+1)^2-4z^2)^{3/2}. \end{gather*} If we change to spherical coordinates we get

(see https://en.wikipedia.org/wiki/Spherical_coordinate_system) \begin{gather*} I = \int_{\mathbb R^3}\dfrac{x^2+y^2+1-z^2}{((x^2+y^2+z^2+1)^2-4z^2)^{3/2}}\, dxdydz = \\[2ex] 2\pi\int_{0}^{\infty}\int_{0}^{\pi}\dfrac{r^2(\sin^2\theta-\cos^2\theta)+1}{((r^2+1)^2-4r^2\cos^2\theta)^{3/2}}r^{2}\sin\theta\, d\theta dr. \end{gather*} If we substitute $\cos\theta =w $ we get an even integrand wrt $w$. Thus \begin{gather*} I=4\pi\int_{0}^{\infty}\int_{0}^{1}\dfrac{r^{2}(r^2+1-2r^2w^2)}{((r^2+1)^2-4r^2w^2)^{3/2}}\, dw dr. \end{gather*} If we proceed straight forward we will be disturbed by $|r^2-1|$. To avoid that we split $I=I_1 +I_2$ where \begin{equation*} I_1 =4\pi\int_{0}^{1}\int_{0}^{1}\dfrac{r^{2}(r^2+1-2r^2w^2)}{((r^2+1)^2-4r^2w^2)^{3/2}}\, dw dr \end{equation*} and \begin{equation*} I_2 =4\pi\int_{1}^{\infty}\int_{0}^{1}\dfrac{r^{2}(r^2+1-2r^2w^2)}{((r^2+1)^2-4r^2w^2)^{3/2}}\, dw dr. \end{equation*} In $I_2$ we change $r$ to $1/r$. We get \begin{equation*} I_2 =4\pi\int_{0}^{1}\int_{0}^{1}\dfrac{r^2+1-2w^2}{((r^2+1)^2-4r^2w^2)^{3/2}}\, dw dr. \end{equation*} Then we add $I_1$ and $I_2$ and get \begin{equation*} I= 4\pi\int_{0}^{1}\int_{0}^{1}\dfrac{(r^2+1)^{2}-2(r^4+1)w^2}{((r^2+1)^2-4r^2w^2)^{3/2}}\, dw dr. \end{equation*} Finally we make the substitution $w=\dfrac{r^2+1}{2r}\sin v.$ \begin{gather*} I=2\pi\int_{0}^{1}\int_{0}^{\arcsin\frac{2r}{r^2+1}}\dfrac{1-\dfrac{r^4+1}{2r^2}\sin^2 v}{r(1-\sin^2v)^{3/2}}\cos v\, dvdr =\\[2ex] 2\pi\int_{0}^{1}\int_{0}^{\arcsin\frac{2r}{r^2+1}}\left(\dfrac{r^4+1}{2r^3}-\dfrac{(1-r^2)^2}{2r^3\cos^2v}\right)\, dvdr =\\[2ex] \pi\int_{0}^{1}\left[\dfrac{r^4+1}{r^3}v-\dfrac{(1-r^2)^2}{r^3}\tan v\right]_{0}^{\arcsin\frac{2r}{r^2+1}}\, dr =\\[2ex] \pi\int_{0}^{1}\left(\dfrac{r^4+1}{r^3}\arcsin\frac{2r}{r^2+1}-\dfrac{(1-r^2)^2}{r^3}\dfrac{2r}{1-r^2}\right)\, dr = [\mbox{ integration by parts }]\\[2ex] \pi\left[\dfrac{r^4-1}{2r^2}\arcsin\dfrac{2r}{r^2+1}+r+\dfrac{1}{r}\right]_{0} ^{1} =2\pi. \end{gather*}

Probably the fastest method: Do not rewrite the denominator and integrate with respect to $D$ directly: \begin{align} I &= \pi \int \limits_{-\infty}^{\infty} \int \limits_0^\infty \frac{D + 1 - z^2}{[(D + (1-z)^2)(D + (1 + z)^2)]^{3/2}} \, \mathrm{d} D \, \mathrm{d} z \\ &= \pi \int \limits_{-\infty}^{\infty} \left[\frac{D - (1 - z^2)}{2\sqrt{(D + (1-z)^2)(D + (1 + z)^2)}}\right]_{D=0}^{D=\infty} \, \mathrm{d} z \\ &= \pi \int \limits_{-\infty}^{\infty} \frac{1 + \operatorname{sgn}(1-z^2)}{2} \, \mathrm{d} z = \pi \int \limits_{-1}^1 \, \mathrm{d} z = 2 \pi \, . \end{align} However, it seems difficult to come up with this antiderivative without a CAS or a lucky guess.

---

Alternatively, we can use Feynman parameters to write \begin{align} I &= \pi \int \limits_{-\infty}^{\infty} \int \limits_0^\infty (D + 1 -z^2)\left[\frac{8}{\pi} \int \limits_0^1 \frac{\sqrt{u(1-u)}}{[(D+(1-z)^2)u + (D+(1+z)^2)(1-u)]^3} \, \mathrm{d} u\right] \, \mathrm{d} D \, \mathrm{d} z \\ &= 8 \int \limits_{-\infty}^{\infty} \int \limits_0^\infty \int \limits_0^1 \frac{(D + 1 -z^2)\sqrt{u(1-u)}}{[D + (1+z)^2 - 4 z u]^3} \, \mathrm{d} u \, \mathrm{d} D \, \mathrm{d} z \, . \end{align} Fubini's theorem allows us to integrate with respect to $D$ first, which is now straightforward (partial fractions). We find \begin{align} I &= 8 \int \limits_0^1 \sqrt{u(1-u)} \int \limits_{-\infty}^\infty \frac{1 + z - 2 z u}{[(1+z)^2 - 4 z u]^2} \, \mathrm{d} z \, \mathrm{d} u \\ &= 8 \int \limits_0^1 \sqrt{u(1-u)} \int \limits_{-\infty}^\infty \frac{1 + (1-2u) z}{[(z + 1-2u)^2 + 4 u (1-u)]^2} \, \mathrm{d} z \, \mathrm{d} u \\ &= 8 \int \limits_0^1 \sqrt{u(1-u)} \int \limits_{-\infty}^\infty \frac{4 u (1-u) \color{red}{+ (1-2u) \zeta}}{[\zeta^2 + 4 u (1-u)]^2} \, \mathrm{d} \zeta \, \mathrm{d} u \end{align} after introducing $\zeta = z + 1 - 2u$ . The integral of the red term vanishes by symmetry, so the substitution $\zeta = 2 \sqrt{u(1-u)} \tau$ yields $$ I = 4 \int \limits_0^1 \int \limits_{-\infty}^\infty \frac{\mathrm{d} \tau}{(1+\tau^2)^2} \, \mathrm{d} u = 4 \int \limits_{-\infty}^\infty \frac{\mathrm{d} \tau}{(1+\tau^2)^2} \, .$$ The final integral is a typical exercise in residue calculus, but integration by parts also does the trick: $$ I = 4 \int \limits_{-\infty}^\infty \left[\frac{1}{1+\tau^2} - \tau \frac{\tau}{(1+\tau^2)^2}\right] \, \mathrm{d} \tau = 4 \left[\pi - \frac{1}{2} \int \limits_{-\infty}^\infty \frac{\mathrm{d}\tau}{1+\tau^2}\right] = 2 \pi \, .$$