Closed form for $\int_0^1...\int_0^1\frac{1}{\left(1+\sqrt{1+x_1^2+...+x_n^2}\right)^{n+1}}\;dx_1...dx_n$ I am wondering whether there is a closed form for the following integral for $n\in\mathbb{N}$:

$$\gamma(n)=\int_0^1...\int_0^1\frac{1}{\left(1+\sqrt{1+x_1^2+...+x_n^2}\right)^{n+1}}\;dx_1...dx_n\tag{*}$$

Particular values which I am aware of include:
$$\gamma(1)=\frac{4\sqrt{2}-5}{3}$$
$$\gamma(2)=\frac{5}{4}-\frac{9\sqrt{3}}{8}+\frac{\pi}{4}$$
Both of these values were obtained by evaluating different integrals to the above which solved the same problem (see below), and I am not sure how to attack $(*)$ directly; I can see no good way of solving this integral. I am wondering especially about the value of $\gamma(3)$; does it also have a simple closed form? Is there a closed form for all $n$?

Background: I posted a question here about calculating the proportion $p(n)$ of an $n$-cube closer to the centre than to the outside, which seems to me like an interesting problem. The $n=2$ case is simple to solve in terms of the following integral:
$$p(2)=2\int_{0}^{\sqrt{2}-1}\frac{1-x^2}{2}-x\;dx=\frac{4\sqrt{2}-5}{3}$$
I was able to write $p(3)$ as follows:
$$p(3)=6\int_{0}^{\frac{\sqrt{3}-1}{2}}\int_{z}^{\sqrt{2-z^2}-1}\frac{1-x^2-z^2}{2}-x\;dx\;dz$$
and I managed to evaluate this to $\frac{5}{4}-\frac{9\sqrt{3}}{8}+\frac{\pi}{4}$, but I was not able to use my method to solve the problem for higher dimensions. In the comments, however, achille hui made a proposition that we have $p(n)=\gamma(n-1)$ for all $n$ and although I still do not perfectly understand his reasoning, the claim does check out numerically for the two values I know already. Furthermore, the new integral is in a nice simple symmetric form (unlike the methods I had been using which required a case-by-case analysis for every dimension, with ugly bounds on the integrals), which makes me hope for a solution method. However, I really cannot see how to go about it. Thus I ask, is there a method for computing the integral $(*)$?
 A: Just for reference, here is the derivation of the formula $p(n) = \gamma(n-1)$ proposed by @achille hui.


*

*Consider the cube $\mathcal{C} = [-1, 1]^n$ and define $\mathcal{D} = \{x \in \mathcal{C} : |x| < \operatorname{dist}(x, \partial \mathcal{C}) \}$. Then for each point $x \in \partial\mathcal{C}$, we can find the unique point $q(x) \in \partial\mathcal{D}$ that lies on the line segment joining $0$ and $x$.
$\hspace{7.5em}$ 
Using this, define $r(x)$ as the ratio
$$ r(x) = \frac{\text{[length of the line segment from $0$ to $q(x)$]}}{\text{[length of line segment from $0$ to $x$]}} = \frac{|q(x)|}{|x|}. $$

*Let me first give a geometric argument. Let $d\mathcal{S}$ be an infinitesimally small portion of the surface $\partial\mathcal{C}$ near $x$ and consider the cone $d\mathcal{V}$ with the base $d\mathcal{S}$ and the vertex $0$.
$\hspace{7.5em}$ 
Then roughly speaking, the portion of $d\mathcal{V}$ that intersects $\mathcal{D}$ is similar to the cone $d\mathcal{V}$ with ratio $r(x)$. Thus their volume roughly satisfies
$$ |\mathcal{D} \cap d\mathcal{V}|
\approx r(x)^n |d\mathcal{V}|
= r(x)^n \cdot \frac{1}{n}|d\mathcal{S}|. $$
From this, we have
$$ p(n)
= \frac{|\mathcal{D}|}{|\mathcal{C}|}
= \frac{\int |\mathcal{D} \cap d\mathcal{V}|}{\int |d\mathcal{V}|}
= \frac{\int_{\partial \mathcal{C}} r(x)^n \frac{1}{n} \, dA}{\int_{\partial \mathcal{C}} \frac{1}{n} \, dA}. \tag{*} $$

*Here is a rigorous justification of the argument above. Write
\begin{align*}
|\mathcal{D}|
 = \int_{\mathcal{C}} \mathbf{1}_{\mathcal{D}}(x) \, dx
&= \int_{0}^{1} \int_{\partial[-s,s]^n} \mathbf{1}_{\mathcal{D}}(y) \, dyds \\
&= \int_{0}^{1} \int_{\partial\mathcal{C}} \mathbf{1}_{\mathcal{D}}(s\omega) \, s^{n-1}d\omega ds \\
&= \int_{\partial\mathcal{C}} \int_{0}^{1} s^{n-1} \mathbf{1}_{\mathcal{D}}(s\omega) \, ds d\omega.
\end{align*}
The note that $s\omega \in \mathcal{D}$ if and only if $s|\omega| < |q(\omega)|$, or equivalently, $s < r(\omega)$. From this, we have
$$ |\mathcal{D}|
= \int_{\partial\mathcal{C}} \int_{0}^{r(\omega)} s^{n-1} \, ds d\omega
= \int_{\partial\mathcal{C}} \frac{1}{n} r(\omega)^n \, d\omega. $$
Replacing $\mathbf{1}_{\mathcal{D}}$ by $\mathbf{1}$ from the above computation, we also get
$$ |\mathcal{C}|
= \int_{\partial\mathcal{C}} \frac{1}{n} \, d\omega. $$
This proves $\text{(*)}$.

*Now we give an explicit formula for $\text{(*)}$. Consider the top face of $\partial\mathcal{C}$. This face can be written as $[-1,1]^{n-1}\times\{1\}$. Then for each point $x = (x',1) \in [-1, 1]^{n-1}\times\{1\}$, it is not hard to verify that
$$ r(x) = \frac{1}{1 + \sqrt{|x'|^2+1}}. $$
(Using the rotational symmetry, it boils down to proving this when the dimension is $n = 2$.) Plugging this back and exploiting the symmetry, we finally have
$$ p(n) = \int_{[0,1]^{n-1}} r(x',1)^n \, dx' = \int_{[0,1]^{n-1}} \frac{1}{\left( 1 + \sqrt{|x'|^2+1} \right)^n} \, dx' = \gamma(n-1). $$
A: Not a full solution, and too long for a comment, but perhaps a plausible strategy forward. The final formula is written in terms of a finite sum of double integrals, which may (or may not) be cracked.
First, rationalize by multiplying up and down by $\left(1-\sqrt{1+x_1^2+...+x_n^2}\right)^{n+1}$
$$
\gamma(n)=(-1)^{n+1}\int_0^1...\int_0^1\frac{\left(1-\sqrt{1+x_1^2+...+x_n^2}\right)^{n+1}}{\left(x_1^2+...+x_n^2\right)^{n+1}}\;dx_1...dx_n\ .
$$
Next, use the identity
$$
\frac{1}{X^{n+1}}=\int_0^\infty d\xi\frac{\xi ^n \exp (-\xi  X)}{\Gamma (n+1)}
$$
to rewrite your integral as
$$
\gamma(n)=\frac{(-1)^{n+1}}{\Gamma(n+1)}\int_0^\infty d\xi\ \xi^n
\int_0^1...\int_0^1\left(1-\sqrt{1+x_1^2+...+x_n^2}\right)^{n+1} e^{-\xi \left(x_1^2+...+x_n^2\right)}\;dx_1...dx_n\ 
$$
$$
=\frac{(-1)^{n+1}}{\Gamma(n+1)}\int_{1}^{n+1} dr\ (1-\sqrt{r})^{n+1}\int_0^\infty d\xi\ \xi^n
\int_0^1...\int_0^1\delta(r-(1+x_1^2+...+x_n^2)) e^{-\xi \left(x_1^2+...+x_n^2\right)}\;dx_1...dx_n
$$
$$
=\frac{1}{\Gamma(n+1)}\int_{1}^{n+1} dr\ (\sqrt{r}-1)^{n+1}\int_0^\infty d\xi\ \xi^n F_n(r,\xi)\ ,
$$
where 
$$
F_n(r,\xi)=\int_0^1...\int_0^1\delta(r-(1+x_1^2+...+x_n^2)) e^{-\xi \left(x_1^2+...+x_n^2\right)}\;dx_1...dx_n\ .
$$
Therefore
$$
F_n(r,\xi)=e^{-\xi (r-1)}\int_0^1...\int_0^1\delta(r-(1+x_1^2+...+x_n^2)) \;dx_1...dx_n
$$
$$
=e^{-\xi (r-1)}\int_{-\infty}^\infty\frac{dk}{2\pi}e^{\mathrm{i}k(r-1)}\left[\int_0^1 dx\ e^{-\mathrm{i}kx^2}\right]^n\ ,
$$
using the integral representation of Dirac delta.
Writing $\exp(-\mathrm{i}kx^2)=\cos(kx^2)-\mathrm{i}\sin(kx^2)$ and using 
$$
\int_0^1 dx\cos(k x^2)=\frac{\sqrt{\frac{\pi }{2}} C\left(\sqrt{|k|} \sqrt{\frac{2}{\pi }}\right)}{\sqrt{|k|}}
$$
and
$$
\int_0^1 dx\sin(k x^2)=\mathrm{sign}(k)\frac{\sqrt{\frac{\pi }{2}} S\left(\sqrt{|k|} \sqrt{\frac{2}{\pi }}\right)}{\sqrt{|k|}}\ ,
$$
in terms of Fresnel integrals [Mathematica notation], and using the binomial theorem we get
$$
F_n(r,\xi)=e^{-\xi (r-1)}\sum_{\ell=0}^n{n\choose \ell}(-\mathrm{i})^{n-\ell}
\int_{-\infty}^\infty\frac{dk}{2\pi}e^{\mathrm{i}k(r-1)}\left(\frac{\pi}{2|k|}\right)^{n/2}C^\ell\left(\sqrt{|k|} \sqrt{\frac{2}{\pi }}\right)S^{n-\ell}\left(\sqrt{|k|} \sqrt{\frac{2}{\pi }}\right)(\mathrm{sign}(k))^{n-\ell}\ .
$$
Putting everything together
$$
\gamma(n)=\frac{1}{\Gamma(n+1)}\sum_{\ell=0}^n{n\choose \ell}(-\mathrm{i})^{n-\ell}
\int_{-\infty}^\infty\frac{dk}{2\pi}\left(\frac{\pi}{2|k|}\right)^{n/2}C^\ell\left(\sqrt{|k|} \sqrt{\frac{2}{\pi }}\right)S^{n-\ell}\left(\sqrt{|k|} \sqrt{\frac{2}{\pi }}\right)(\mathrm{sign}(k))^{n-\ell}\times\int_1^{n+1}dr (\sqrt{r}-1)^{n+1}
e^{\mathrm{i}k(r-1)}\int_0^\infty d\xi\ \xi^n e^{-\xi (r-1)}\ .
$$
Performing the $\xi$-integral first
$$
\gamma(n)=\frac{1}{\Gamma(n+1)}\sum_{\ell=0}^n{n\choose \ell}(-\mathrm{i})^{n-\ell}
\int_{-\infty}^\infty\frac{dk}{2\pi}\left(\frac{\pi}{2|k|}\right)^{n/2}C^\ell\left(\sqrt{|k|} \sqrt{\frac{2}{\pi }}\right)S^{n-\ell}\left(\sqrt{|k|} \sqrt{\frac{2}{\pi }}\right)(\mathrm{sign}(k))^{n-\ell}\times\int_1^{n+1}dr (\sqrt{r}-1)^{n+1}
e^{\mathrm{i}k(r-1)}\frac{1}{(r-1)^{n+1}}\ .
$$
