distance from the centre of a $n$-cube as $n \rightarrow \infty$

[Note: Added values for five-dimensional hypercube.]

Consider that the space is a hypercube, so each coordinate is independently distributed. The square of that coordinate's difference from $\frac12$ has the pdf

$$ f(x) = \begin{cases} \frac{1}{\sqrt{x}} & 0 \leq x \leq \frac14 \\ 0 & \text{elsewhere} \end{cases} $$

This distribution has a mean of $\frac{1}{12}$ and a variance of $\frac{1}{180}$. As $n$ increases without bound, the squared distance of the point from the hypercube's center is the sum of $n$ independent and identically distributed (i.i.d.) variables with that same distribution, and is thus asymptotically normally distributed (by the central limit theorem) with mean $\frac{n}{12}$ and variance $\frac{n}{180}$. For instance, for $n = 180$, we would have a mean squared distance of $15$ and a variance of $1$. That variance is small enough already that you could just take the square root of the mean squared distance and probably get a very good approximation of the mean distance.

By that logic, the mean distance would be asymptotically $\sqrt{\frac{n}{12}}$, approached from below, since the square root of a nearly normal distribution with a positive mean would be skewed that way.

For $n = 1, 2, 3, 4, 5$, this expression yields approximate mean distances of $0.289, 0.408, 0.500, 0.577, 0.645$, which compares reasonably well with the more accurate values given in the OP ($0.250, 0.383, 0.480, 0.561, 0.631$). These latter values appear to be approaching the asymptotic expression from below, as expected, but are already not too far off.

ETA ($2019$-$02$-$13$, five-dimensional case added $2020$-$04$-$24$): A second-order analysis yields $\sqrt{\frac{5n-1}{60}}$, for which the values for $n = 1, 2, 3, 4, 5$ are $0.258, 0.387, 0.483, 0.563, 0.632$, respectively, showing even closer agreement.

This does not solve the question. This answer contains the details of an attempt that ultimately failed, but can perhaps still be of some help.

We can try to bound it below by the sphere of radius $1/2$ and above by the sphere of radius $\sqrt{n}/2$. The calculations are somewhat simpler in hyperspherical coordinates.

You will have that the bounds $B_n(R)$ are given by

$$ \int_{\phi_{n-1}=0}^{2\pi} \int_{\phi_{n-2}=0}^{\pi} \dots \int_{\phi_1=0}^{\pi} \int_{r=0}^R r^n\,(\sin \phi_1)^{n-2}(\sin \phi_2)^{n-3}\dots\,(\sin \phi_{n-2}) \,dr\,d\phi_{1}\,d\phi_{2}\dots d\phi_{n-1}$$

This is really just a product of integrals. Notice that $\rho(n,r)=\int_{r=0}^R\, r^n\, dr=\frac1{n+1}R^{n+1}$ and let

$$I(n)=\int_{\phi=0}^\pi\,(\sin\phi)^n\,d\phi$$

Then

$$B_n(R)=\frac{2\pi R^{n+1}}{n+1}\cdot \prod_{k=1}^{n-2}I(k)$$

Now, via integration by parts, one can find the following recurrence relation between the $I(k)$, which holds for $k\geq 2$:

$$I(k)=\frac{k-1}{k}\cdot I(k-2)$$

Humm, looks interesting. How do some products $P(n)=\prod_{k=1}^{n-2}I(k)$ end up, in light of this? We have:

\begin{align} P(3)&=I(1)\\ \\ P(4) &=I(1)\cdot I(2)\\ &=\frac12 \cdot I(0)\cdot I(1)\\ \\ P(5) &=\frac12 \cdot I(0)\cdot I(1) \cdot I(3)\\ &=\left[\frac12 \cdot \frac23\right]\cdot I(0)\cdot I(1)^2\\ \\ P(6) &=\left[\frac12 \cdot \frac23\right]\cdot I(0)\cdot I(1)^2 \cdot I(4)\\ &=\left[\frac12 \cdot \frac23\cdot \frac34\right]\cdot I(0)\cdot I(1)^2 \cdot I(2)\\ &=\left[\left(\frac12\right)^2 \cdot \frac23\cdot \frac34\right]\cdot I(0)^2\cdot I(1)^2\\ \\ P(7) &=\left[\left(\frac12\right)^2 \cdot \frac23\cdot \frac34\right]\cdot I(0)^2\cdot I(1)^2 \cdot I(5)\\ &=\left[\left(\frac12\right)^2 \cdot \left(\frac23\right)^2\cdot \frac34\cdot\frac45\right]\cdot I(0)^2\cdot I(1)^3\\ \\ P(8) &=\left[\left(\frac12\right)^2 \cdot \left(\frac23\right)^2\cdot \frac34\cdot\frac45\right]\cdot I(0)^2\cdot I(1)^3 \cdot I(6)\\ &=\left[\left(\frac12\right)^3 \cdot \left(\frac23\right)^2\cdot \left(\frac34\right)^2\cdot\frac45\cdot\frac56\right]\cdot I(0)^3\cdot I(1)^3\\ \end{align}

I hope the pattern shows how induction goes. Noting that $I(0)=\pi$ and $I(1)=2$ we have

$$P(n)=c_n\cdot \pi^{\lfloor\frac{n}2\rfloor-1}\cdot 2^{\lfloor\frac{n-1}2\rfloor}$$

where $c_n$ is the constant given by the product of fractions. The product is reminiscent of a telescoping product, and will feature all fractions of the form $\frac{k}{k+1}$ from $\frac12$ to $\frac{n-3}{n-2}$. The last two fractions will be raised to the power $1$, the preceding two fractions by the power $2$, and so on and so forth. With cancellations, we can recognize that

$$c_n=\frac1{(n-2)!!}$$

where $k!!$ denotes the double factorial of $k$.

Putting it all together we get:

$$B_n(R)=\frac{2^{\lfloor\frac{n+1}2\rfloor}\cdot\pi^{\lfloor\frac{n}2\rfloor}\cdot R^{n+1}}{(n+1)\cdot(n-2)!!} $$

And then $B_n(1/2)\leq A_n\leq B_n(\sqrt{n}/2)$. Would be interesting to get some asymptotic estimates here. It does not look immediately clear to me, and I'm not sure this solve the problem.

---

It appears one may express the double odd factorial in terms of the gamma function via

$$(2k-1)!!=\frac{2^k\cdot\Gamma\left(k+\frac12\right)}{\sqrt{\pi}}$$

Hence, for odd $n=2k+1$ we have

$$(n-2)!!=\frac{2^{\frac{n-1}2}\cdot\Gamma\left(\frac{n}2\right)}{\sqrt{\pi}}$$

It follows that

\begin{align} B_{2k+1}(R) &=\frac{2\cdot\pi^{\left(k+\frac12\right)}\cdot R^{2k+2}}{(2k+2)}\\ &=\frac{2\cdot{\left(\sqrt{\pi}\right)}^{2k+1}\cdot R^{2k+2}}{(2k+2)}\\ &=\frac{1}{(k+1)\,\sqrt{\pi}}\cdot{\left(R\,\sqrt{\pi}\right)}^{2k+2}\\ &=\frac{1}{(k+1)\,\sqrt{\pi}}\cdot{\left(R^2\,\pi\right)}^{k+1} \end{align}

It's clear that as $k\to \infty$, $B_{2k+1}(R)\to 0$ if $R^2\,\pi\leq1$ and $\to\infty$ if $R^2\,\pi>1$. In particular, for $R=1/2$ the limit of the lower bound is $0$, so it does not solve our problem.