In a reply to Markiyan Hirnyk's recent question on this topic, Alec Mihailovs posted solutions 
for  n=1..7 for the mean distance between random points in an n-cube.

I also posted some code to simulate the value for arbitrary n.

In this post I'm going to show the results of simulations for n=1..10, along with histograms 
of the results (with superimposed normal pdf curve) and normal probability plots, since the 
simulations rely on the central limit theorem. % errors, for n=1..7 where the result is known, 
and 95% confidence intervals for those where the result is not known (by me), are also shown.

ncubem:=proc(n::posint,N::posint,m)
  # procedure simulates the mean value between two random points A and B 
  # in a n-cube of face length m, N times.
  local A,B,s,i,RV,dim:
  RV := Statistics:-RandomVariable(Uniform(0,m)):
  s:=0:
  for dim from 1 to n do
    A[dim]:=Statistics:-Sample(RV,N);
    B[dim]:=Statistics:-Sample(RV,N);
  end do:
  for i from 1 to N do
    s:=s+sqrt(sum((A[dimm][i]-B[dimm][i])^2,dimm=1..n))
  end do:
  s/N;
end proc:

unitncubesim:=proc(n::posint,N::posint,monte::posint,exact)
# procedure samples the mean value of a unit n-cube N times
# exact is an expression for the known exact value. If this 
# is provided, the statistical error is computet. If it is not
# then a 95% confidence interval for the mean is computed.
local sample,i,mu,sigma,qtl,ciL,ciU,p1,p2,p3,XX,A,pr1;
sample:=seq(ncubem(n,monte,1),i=0..N):
mu:=Statistics:-Mean([sample]);
sigma:=Statistics:-StandardDeviation([sample]):
printf("simulated value for %d-cube: %0.6f",n,mu);

if type(exact,algebraic) then
  printf("; known value: %0.6f,    %0.4f%% error",evalf(exact),100*evalf(abs(mu-exact)/exact,6));
else
  qtl:=Statistics:-Quantile(Normal(0,1),evalf(1-((1-(95/100))/2))):
  ciL:=evalf(mu-(qtl*sigma/sqrt(N)),5):ciU:=evalf(mu+(qtl*sigma/sqrt(N)),5):
  printf(" with 95%% confidence interval [%0.6f,%0.6f]",ciL,ciU);
end if;

p1:=Statistics[Histogram]([sample],bincount=15):

XX:=Statistics:-RandomVariable(Normal(mu,sigma)):
p2:=plot(Statistics:-PDF(XX,x),x=mu-(4*sigma)..mu+(4*sigma)):
p1:=plots:-display(p1,p2);
p3:=Statistics:-ProbabilityPlot([sample],Normal(mu,sigma),thickness=2,color=BLACK):
print(plots:-display(Array([p1,p3])));
end proc:




f:=n->evalf(2^n*Int(sqrt(add((x||i)^2,i=1..n))*mul(1-x||i,i=1..n),
    [seq(x||i=0..1,i=1..n)])):
#HT Alec Mihailovs for this solution for n=1..7

F:=seq(f(n),n=1..7):

Digits:=15:

t:=time():
unitncubesim(1,2000,500,F[1]);

    simulated value for 1-cube: 0.333516; known value: 0.333333,    0.0549% error

unitncubesim(2,1500,500,F[2]);

    simulated value for 2-cube: 0.521228; known value: 0.521405,    0.0339% error

unitncubesim(3,1000,500,F[3]);

    simulated value for 3-cube: 0.661723; known value: 0.661707,    0.0024% error

unitncubesim(4,900,400,F[4]);

    simulated value for 4-cube: 0.777830; known value: 0.777666,    0.0211% error

unitncubesim(5,800,300,F[5]);

    simulated value for 5-cube: 0.878373; known value: 0.878531,    0.0180% error

unitncubesim(6,700,300,F[6]);

    simulated value for 6-cube: 0.969059; known value: 0.968942,    0.0121% error

unitncubesim(7,650,300,F[7]);

    simulated value for 7-cube: 1.052150; known value: 1.051584,    0.0542% error

unitncubesim(8,600,300,"");

    simulated value for 8-cube: 1.128114 with 95% confidence interval [1.12700,1.12920]

unitncubesim(9,550,300,"");

    simulated value for 9-cube: 1.200337 with 95% confidence interval [1.199100,1.201500]

unitncubesim(10,500,300,"");

    simulated value for 10-cube: 1.267209 with 95% confidence interval [1.266000,1.268400]

time()-t:
%/60;

                  8.0751