Vrftsjtryk

Denote $square_m={pmb{x}=(x_1,dots,x_m)inmathbb{R}^m: 0leq x_ileq1,,,forall i}$ be an $m$-dimensional cube.

It is all too familiar that $int_{square_1}frac{dx}{1+x^2}=frac{pi}4$.

QUESTION. If $VertcdotVert$stands for the Euclidean norm, then is this true?
$$int_{square_{2n-1}}frac{dpmb{x}}{(1+Vertpmb{x}Vert^2)^n}
=frac{pi^n}{4^nn!}.$$

Yes, this is true, it follows from formulas in Higher-Dimensional Box Integrals, by Jonathan M. Borwein, O-Yeat Chan, and R. E. Crandall.

begin{align}
&text{define};;C_{m}(s)=int_{[0,1]^m}(1+|vec{r}|^2)^{s/2},dvec{r},;;text{we need};;C_{2n-1}(-2n).\
&text{the box integral is};;B_m(s)=int_{[0,1]^m}|vec{r}|^s,dvec{r},\
&text{related to our integral by};;2n C_{2n-1}(-2n)=lim_{epsilonrightarrow 0}epsilon B_{2n}(-2n+epsilon)equiv {rm Res}_{2n}.
end{align}
The box integral $B_n(s)$ has a pole at $n=-s$, with "residue" ${rm Res}_n$. Equation 3.1 in the cited paper shows that the residue is a piece (an $n$-orthant) of the surface area of the unit $n$- sphere, so it is readily evaluated,
begin{align}
&{rm Res}_{n}=frac{1}{2^{n-1}}frac{pi^{n/2}}{Gamma(n/2)},\
&text{hence};;C_{2n-1}(-2n)=frac{1}{2n}frac{1}{2^{2n-1}}frac{pi^{n}}{Gamma(n)}=frac{1}{4^n}frac{pi^n}{n!};;text{as in the OP}.
end{align}

More generally, we can consider

$$C_{p-1}(-p)=int_{[0,1]^{p-1}}(1+|vec{r}|^2)^{-p/2},dvec{r}=frac{1}{p},{rm Res}_p=frac{1}{p} frac{1}{2^{p-1}}frac{pi^{p/2}}{Gamma(p/2)}. $$
This formula holds for even and odd $pgeq 2$, as surmised by Gro-Tsen in a comment.

There are similarly remarkable hypercube integrals where this came from, for example

$$int_{[0,pi/4]^k}frac{dtheta_1 dtheta_2cdots dtheta_k}{(1+1/cos^2theta_1+1/cos^2theta_2cdots+1/cos^2theta_k)^{1/2}}=frac{k!^2pi^k}{(2k+1)!}.$$

Here is an elementary proof. Using the formula
begin{equation}
frac1{a^n}=frac1{Gamma(n)},int_0^infty u^{n-1} e^{-u,a},du tag{1}
end{equation}
with $a=1+|x|^2$,
denoting your integral by $J_n$, letting $I:=[0,1]$, $phi(z):=frac1{sqrt{2pi}},e^{-z^2/2}$, and $Phi(z):=int_{-infty}^zphi(t),dt$, and making substitutions $x_1=z_1/sqrt{2u}$ and then $u=z^2/2$, we have
begin{align}
J_n&=frac1{Gamma(n)},int_0^infty u^{n-1} e^{-u},du int_{I^{2n-1}}e^{-u,|x|^2},dx \
&=frac1{Gamma(n)},int_0^infty u^{n-1} e^{-u},du Big(int_I e^{-u,x_1^2},dx_1Big)^{2n-1} \
&=frac{pi^{n-1/2}}{Gamma(n)},int_0^infty u^{-1/2} e^{-u},du Big(Phi(sqrt{2u})-frac12Big)^{2n-1} \
&=frac{2pi^{n}}{Gamma(n)},int_0^infty Big(Phi(z)-frac12Big)^{2n-1},phi(z),dz \
&=frac{2pi^{n}}{Gamma(n)},int_0^infty Big(Phi(z)-frac12Big)^{2n-1},dPhi(z) \
&=frac{pi^{n}}{4^n Gamma(n+1)},
end{align}
as desired.

This derivation obviously holds whenever $2n-1$ is a positive integer.

Further comments:

More generally, in view of Bernstein's theorem on completely monotone functions, one can quite similarly express the box integral
begin{equation}
int_{I^n} g(|x|^2),dx
end{equation}
as an ordinary integral over $[0,infty)$ -- for any completely monotone function $g$. Indeed, Bernstein's theorem states that any completely monotone function is a positive mixture of decreasing exponential functions; identity (1) is then such an explicit Bernstein representation of the completely monotone function $amapsto frac1{a^n}$.

Now, furthermore, one does not have to insist that the mixture be positive, since we are dealing here with identities rather than inequalities. So, any function $g$ that can be represented as the mixture
begin{equation}
g(a)=int_0^infty e^{-u,a},mu(du)
end{equation}
of decreasing exponential functions $amapsto e^{-u,a}$
for $a>0$, where $mu$ is any finite, possibly signed ("mixing") measure, will do just as well.