Average minimum distance between $n$ points generate i.i.d. uniformly in the ball Let $U \in \mathbb{R}^3$ be distributed uniformly in the Ball in $\mathbb{R}^3$ centered at zero.   That is $U \sim  f_U(u)= \frac{1}{ \frac{4}{3} \pi R^3}$  for all $\|u\|\le R$  where $R$ is the radius of the ball. 
Now suppose we generate $n$ points i.i.d. according distribution of   $U$. 
My questions is: Can we compute the expected minimum distance between the generated points,
that is
\begin{align}
E\left[ \min_{i,j\in \{1,2,,,n\}} \| U_i-U_j\| \right],
\end{align} 
where  $\| U_i-U_j\|$ is Euclidean distance. 
This question is related to a number of other questions. 
For example,  Average distance between two random points in a square
Average minimum distance between $n$ points generate i.i.d. with uniform dist.
I feel that this question should have been addressed before but not sure where to look.
There is a conjecture that the minimum distance behaves as $\frac{1}{n^{\frac{2}{3}}}$ but I am not sure how to show this? 
Update See a recently add proof of this statement for the case when 'border' effects are negligible. That is the answer is asymptotic.  The question know is how to take into account the border effects? 
Thank you very much. 
 A: Here's an asymptotic bound - hopefully it's tight.
We throw $N$ little balls of diameter $D$ randomly (uniformly) inside the big sphere of radius $R$. 

UPDATE: The new version (lower half) is better than what follows.

Neglecting border effects (reasonable if $N$ is large) the probability that the ball $i$ is "free" (no other overlaps with it) is 
$$P(F_i)=\left(1-v(D)\right)^{N-1}  \tag{1}$$
where $v(D) \triangleq D^3/R^3$
The probability that all balls are free can be bounded as
$$P(\cap F_i)=1- P(\cup F_i^c)\ge 1 - N(1-P(F_i)) \triangleq g(D,N) \tag{2}   $$
For large $N$ 
$$g(D,N) \approx 1 - N^2 \, v(D) \tag{3}  $$
in the range where this is positive, ie. $0\le D \le D_1 \triangleq  R/N^{2/3} $
Now, let $t$ be the minimun distance between the sphere centers. Then
$$P(t \ge D) = P(\cap F_i) \ge g(D,N) \tag{4}$$
And then
$$E(t) = \int_0^{\infty} P(t \ge D) dD  \ge\int_0^{D_1} g(D,N) \, dD \approx  D_1- N^2 \frac{ D_1^4}{4 R^3} = \frac{3 }{ 4 } \frac{R}{N^{2/3}} $$

(Simulation data suggests that the order is right, and so is the bound, but the real coefficient is around $1.12$ - perhaps $9/8$)

Update: (Improved version)
A better approach can be obtained by considering instead of $F_i$ (free ball) the event $S_j\equiv$ "separated pair" (pair of balls are separated, they do not overlap) where $j$ indexes the $M=N(N-1)/2 \approx N^2/2$ pairs.
By the same reasoning:
$$P(S_j)=1-v(D) =1 - \frac{D^3}{R^3} \tag{5}$$
$$P(\cap S_i) \ge \max(1 - M(1-P(S_i)),0)= \max\left(1 - M \frac{D^3}{R^3},0\right) \triangleq h(D,M)  \tag{6}   $$
The range where $h(D,M)$  is positive, ie. $0\le D \le D_2 \triangleq  R/M^{1/3} $
Now, let $t$ be the minimun distance between the sphere centers. Then
$$P(t \ge D) = P(\cap S_i) \ge h(D,M) \tag{7}$$
And then
$$E(t) = \int_0^{\infty} P(t \ge D) dD  \ge\int_0^{D_2} h(D,M) \, dD =\\
= \frac{3}{4} \frac{R}{M^{1/3}}
 \approx 0.945 \frac{R}{\sqrt[3]{N(N-1)}} \approx 0.945 \frac{R}{N^{2/3}} \tag{8}$$

Update 2 : A simple heuristic which seems to produce the correct coefficient:
Following the approach above, we could assume that $S_i$ are asympotically independent, and then:
$$P(\cap S_i) \approx \left(1-\frac{D^3}{R^3}\right)^M \tag{9}$$
Then
$$E(t) \approx \int_0^{R}\left(1-\frac{D^3}{R^3}\right)^M  dD =\\= R \, \Gamma(4/3) \frac{\Gamma(M+1)}{\Gamma(M+4/3)} \approx  R \, \Gamma(4/3) M^{-1/3} \approx 1.12508368 \frac{R}{N^{2/3}} \tag{10}$$

Update 3 : Regarding corrections for border effects. 
(Lets assume $R=1$ to save notation, it's just a scale factor)
If we wished to include border effects we should replace $(5)$ (computing the balls intersection as here) by
$$1-D^3+\frac{9}{16}D^4 -\frac{1}{32}D^6 \hspace{1cm} 0\le D\le 2$$
The integral gets more complicated, but the (first order) asymptotic result is not altered:
Lemma: For any positive differentiable function $g(x)$ which , in $[0,+\infty)$, has global maximum at $g(0)=1$, and which has zero first and second derivates  $g(x)=1-a x^3 + O(x^4)$ we have (variation of Laplace method, see eg here sec 2.1.3)
$$ \int_0^\infty g(x)^M dx = \frac{\Gamma(1/3)}{3 a^{1/3}} M^{-1/3}+ o(M^{-1/3})$$
which again leads as to $(10)$.
