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)$.
Best Answer
Assume $a=1$ for simplicity. Let $M$ be the minimum distance among these $n$ points. Note $M$ is always at most $\frac1{n-1}$, which occurs when the points are evenly spaced. To calculate $EM$, we first calculate $P(M>m)$, for $0\le m \le 1/(n-1)$.
$P(M>m)$ is, by symmetry, equal to $n!$ times $P(M>m\text{ and } U_1<U_2<\dots<U_n)$, and the probability of the latter is the volume of the below subset $S$ of the unit hypercube $Q=[0,1]^n$: $$ S = \{(x_1,\dots,x_n)\in Q:x_i+m\le x_{i+1}\text{ for }i=1,\dots,n-1\} $$ Now, consider the transformation $f:S\to [0,1]^n$, given by $$ (x_1,\dots,x_n)\mapsto (x_1,x_2-m,x_3-2m,\dots,x_n-(n-1)m) $$ A little thought shows that $f$ is a volume-preserving bijection from $S$ to the region $S'$ of the smaller hypercube $Q'=[0,1-(n-1)m]^n$ given by $$ S'=\{(y_1,\dots,y_n)\in Q':y_i\le y_{i+1}\text{ for }i=2,\dots,n-1\} $$ The volume of $S'$ is, by symmetry of permuting the $y_i$, equal to $1/n!$ times the volume of $Q'$, which is just $(1-(n-1)m)^n$. Putting this all together, $$ P(M>m)=n!\cdot \text{Vol}(S)=n!\cdot \frac1{n!}(1-(n-1)m)^n=(1-(n-1)m)^n $$ and therefore $$ \begin{align} EM =\int_0^{1/(n-1)}P(M>m)\,dm &=\int_0^{1/(n-1)}(1-(n-1)m)^n\,dm\\ &=\frac{1}{n-1}\frac{(1-(n-1)m)^{n+1}}{n+1}\Big|_{0}^{1/(n-1)}\\ &=\boxed{\frac{1}{n^2-1}} \end{align} $$ To get the answer for the interval $[0,a]$, simply multiply this result by $a$.