Evaluate $\iiint_{[0,1]^3}\frac{dx\,dy\,dz}{(1+x^2+y^2+z^2)^2}$

As in the title I have to evaluate this triple integral:

$$\iiint_{[0,1]^3}\frac{dx\,dy\,dz}{(1+x^2+y^2+z^2)^2}$$

I've been trying to solve this since a week ago.

The first thing I've done was understand the meaning of the integral. I think this integral represents the mass (as an example) of a unitary cube which contains materials of different density.

The values of the materials density are, point for point, the inverse of the square of spheres centered in the origin plus one.

The max value of the density is $1$ in the origin of the cube and the min value is $\frac{1}{16}$ on the opposite vertex.

I suppose that the value of the integral is $\frac{\pi^2}{32}$

I've tried to use simple substitutions without any results, so I tried to change the coordinates with spherical and cylindrical systems. The spherical coordinates give me an incredibly long sum of integrals and I doubt that they're all integrable as elementary functions.

The cylindrical gives me the following result

$$\frac{\pi^2}{16}-\int_0^\frac{\sqrt2}{2}{\frac{\arctan{\sqrt{\frac{u^2-1}{u^2-2}}}}{\sqrt{2-u^2}}du},$$

which I'm not able to solve.

My instinct tells me that there is a trick in some steps where I can observe that a difficult integral actually is exactly half of another one simpler but I can't figure out where.

I'll appreciate any kind of suggestions.

Solution 1:

Here's an intuitive solution involving surfaces of hypercubes.

Let $$v := \iiint_{[0,1]^3}\frac{dxdydz}{(1+x^2+y^2+z^2)^2}$$

be the value you're looking for.

Suppose a point-shaped light source shines uniformly from the origin of 4D space with intensity of $1$ watt per cubic radian. Let $p$ be a parametrization of a cubic piece $P$ of hyperplane, specifically $p(x,y,z) = (1,x,y,z)$ for $x, y, z \in [0,1]$.

It is clear that no part of this piece of hyperplane is in shadow. If we compute the irradiance $I(x,y,z)$ of arriving light in watts per cubic unit, we find

$$I(x,y,z) = \frac{1}{(1+x^2+y^2+z^2)^2}.$$

If we integrate this irradiance over the surface of $P$, we find that the total power arriving at $P$ is $v$ watts per cubic unit.

Now consider the surface of a hypercube with corners at $(\pm 1, \pm 1, \pm 1, \pm 1)$. We know that $P$ is an octant of one of the eight cells of this surface, and every other octant receives the same irradiance, so we know that the total power arriving at the hypercube is $64v$. However, the hypercube captures every bit of light emitted by our light source, so we know that $64v$ equals $2\pi^2$, the total wattage of light emitted by our light source.

Therefore $v = \frac{2\pi^2}{64} = \frac{\pi^2}{32}$.

Solution 2:

Let us denote the integral by $I$, $$I=\int_{0}^{1} \int_{0}^{1} \int_{0}^{1} \frac{dx dy dz}{(1+x^2+y^2+z^2)^2}=\int_{0}^{1} \int_{0}^{1} \int_{0}^{1} \int_{0}^{\infty} t \exp[-(1+x^2+y^2+z^2)t] dx dy dz dt$$ Due to the symmetry of $I$ in $x,y,z$; we can write $$I=\int_{0}^{\infty} dt ~t e^{-t} \left ( \int_{0}^{1} e^{-tx^2} dx \right)^3=\frac{\pi \sqrt{\pi}}{8} \int_{0}^{\infty} t e^{-t} \mbox{erf}~^3(\sqrt{t}) dt= \frac{\pi^2}{8} \int_{0}^{1} v^3 dv=\frac{\pi^2}{32}. $$ Here $\mbox{erf}(x)=\frac{2}{\sqrt{\pi}} \int_{0} ^{x} e^{-t^2} dt$ and we have used $\mbox{erf}(\sqrt{t})=v$.

I also post the screen shot of from Mathematica,

Solution 3:

Since for any $a>0$ we have $\frac{1}{a^2}=\int_{0}^{+\infty} w e^{-aw}\,dw$, by Fubini's theorem the original integral can be written as

$$ \int_{0}^{+\infty} w e^{-w}\left(\frac{\sqrt{\pi}\,\text{Erf}(\sqrt{w})}{2\sqrt{w}}\right)^3\,dw\stackrel{w\mapsto w^2}{=}\frac{\pi\sqrt{\pi}}{4}\int_{0}^{+\infty}e^{-w^2}\text{Erf}^3(w)\,dw $$ and the RHS is clearly $$ \frac{\pi\sqrt{\pi}}{4}\left[\frac{\sqrt{\pi}}{8}\,\text{Erf}^4(w)\right]_{0}^{+\infty} =\color{red}{\frac{\pi^2}{32}}.$$

Interesting side-note: the same approach in dimension $2$ gives a relation between $\iint_{(0,1)^2}\frac{dx\,dy}{(1+x^2+y^2)^2}=\frac{1}{\sqrt{2}}\arctan\frac{1}{\sqrt{2}}$ and $\int_{0}^{+\infty}\left(1-\text{Erf}^3(w)\right)\,dw$.
In particular it proves that $$ \int_{0}^{+\infty}\left(1-\text{Erf}^3(x)\right)\,dx = \frac{6\sqrt{2}}{\pi\sqrt{\pi}}\,\arctan\frac{1}{\sqrt{2}} $$ which is not recognized by Mathematica (or at least by my version).

Solution 4:

Another solution:

Let us use $x=r \cos \phi,~ y= r \sin \phi$ and write $$I=\int_{0}^{1}\int_{0}^{1} \int_{0}^{1} \frac{dx dy dz}{(1+x^2+y^2+z^2)^2}= 2\int_{0}^{1} \int_{0}^{\pi/4} \int_{0}^{\mbox{sec}\phi} \frac{r dr d\phi dz}{(1+r^2+z^2)^2}~~~(1).$$ $$ \Rightarrow I =\int_{0}^{1} \int_{0}^{\pi/4} \left(\frac{1}{1+z^2}-\frac{1}{(1+z^2+\mbox{sec}^2\phi)}\right) d\phi dx=\int_{0}^{1} \int_{0}^{1} \frac{1}{(1+z^2)(2+z^2+t^2)} dz~ dt,~~~~(2)$$ $t=\mbox{sec}\phi$ taken here. Next, interchange $t$ and $z$ to get $$I=\int_{0}^{1} \int_{0}^{1} \frac{dz dt}{(1+t^2)(2+z^2+t^2)}~~~~(3)$$ Adding (2) and (3), we get $$I=\frac{1}{2} \int_{0}^{1} \int_{0}^{1} \frac{1}{(2+z^2+t^2)} \left( \frac{1}{1+z^2}+\frac{1}{1+t^2} \right)dx~dt=\frac{1}{2} \int_{0}^{1} \int_{0}^{1} \frac{ dz~ dt}{(1+z^2)(1+t^2)}= \frac{1}{2} \left ( \int_{0}^{1}\frac {dz}{1+z^2} \right)^2=\frac{\pi^2}{32}.$$