Denote $\square_m=\{\pmb{x}=(x_1,\dots,x_m)\in\mathbb{R}^m: 0\leq x_i\leq1,\,\,\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 $\Vert\cdot\Vert$stands for the Euclidean norm, then is this true? $$\int_{\square_{2n-1}}\frac{d\pmb{x}}{(1+\Vert\pmb{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}\,d\vec{r},\;\;\text{we need}\;\;C_{2n-1}(-2n).\\ &\text{the box integral is}\;\;B_m(s)=\int_{[0,1]^m}|\vec{r}|^s\,d\vec{r},\\ &\text{related to our integral by}\;\;2n C_{2n-1}(-2n)=\lim_{\epsilon\rightarrow 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}\,d\vec{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 $p\geq 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{d\theta_1 d\theta_2\cdots d\theta_k}{(1+1/\cos^2\theta_1+1/\cos^2\theta_2\cdots+1/\cos^2\theta_k)^{1/2}}=\frac{k!^2\pi^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{2\pi}}\,e^{-z^2/2}$, and $\Phi(z):=\int_{-\infty}^z\phi(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_1\Big)^{2n-1} \\ &=\frac{\pi^{n-1/2}}{\Gamma(n)}\,\int_0^\infty u^{-1/2} e^{-u}\,du \Big(\Phi(\sqrt{2u})-\frac12\Big)^{2n-1} \\ &=\frac{2\pi^{n}}{\Gamma(n)}\,\int_0^\infty \Big(\Phi(z)-\frac12\Big)^{2n-1}\,\phi(z)\,dz \\ &=\frac{2\pi^{n}}{\Gamma(n)}\,\int_0^\infty \Big(\Phi(z)-\frac12\Big)^{2n-1}\,d\Phi(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 $a\mapsto \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 $a\mapsto e^{-u\,a}$ for $a>0$, where $\mu$ is any finite, possibly signed ("mixing") measure, will do just as well.
