How to compute the expected distance to a nearest neighbor in an array of random vectors? Let us have $k$ independent random vectors $x_1, x_2, \dots, x_k$ with uniform distribution over $ \left[0;1 \right]^n $. Then the distance (preferrably Manhattan) between an arbitrary vector $x_a$ and its nearest neighbor among these vectors is a well defined random variable. How can we estimate its momenta? Especially the expected value?
Since I have no idea how to start let us begin with the case of two vectors $a$ and $b$ in two dimensions:
In this case one is the nearest neighbor to the other, so their distance is:
$$
d \left(a, b\right) = \left | a_1 - b_1 \right | + \left | a_2 - b_2 \right | 
$$ 
Where $a_1, a_2, b_1, b_2 \sim  U \left[0;1 \right]$ and are independent.
In this case we can compute the expected value as a double integral:
$$
I := \iint_{\left[0;1 \right]^2}^{} \left | x - y \right | + \left | x - y \right | dxdy
$$
Since the integrand is symmetric along $x = y$:
$$
I = 2 \iint_{\left[0;1 \right]^2}^{} \left | x - y \right | dxdy = 4 \int_{0}^{1} \left ( \int_{0}^{x} x - y  dy \right ) dx = 4 \int_{0}^{1} \frac{x^2}{2} dx = \frac{2}{3}
$$
It would probably not be so hard to extend this result to $n$ dimensions. The real challenge is to add other vectors to the picture.
If we have three vectors $a,b,c$, then both $b$ and $c$ have equal probability of being the closest individual to $a$. We can then compute the expected distance as the weighted sum of conditional expected values.
$$
\frac{1}{2} E \left [ d \left( a, b \right) | b \text{ is closer}  \right] + \frac{1}{2} E \left [ d \left( a, c \right) | c \text{ is closer}  \right]
$$
Both $E \left [ d \left( a, b \right) | b \text{ is closer}  \right]$ and $E \left [ d \left( a, c \right) | c \text{ is closer}  \right]$ could be computed as integrals in 6 dimensions somehow.
 A: Consider the variable $D=|A-B|$ where $A$,$B$, and iid uniform on $[0,1]$
Then, $D$ has a triangular density $f_D(d)=2(1-d) \,[0\le d\le1]$ and characteristic function
$$\psi_D(t)=\frac{2}{t^2}(1+i\, t -e^{i t})= 2\left(\frac{1}{2!}+\frac{i\,t}{3!}+\frac{(i \, t)^2}{4!}+\frac{(i \, t)^3}{5!}+\cdots \right)$$
From this we get the first  moments of $D$ : $E(D)=1/3$, $E(D^2)=1/6$ or $\sigma_D^2=1/18$, etc
Then, in $n$ dimensions, the Manhattan distance between two points is given by the sum of $n$ variables with that density. The resulting density is obtained by the $n$-th convolution of $f_D(d)$, and the CF as $\psi_D^n(t)$. 
Given that, you are insterested in the minimum of $k-1$ realizations of such (independent) variables. That's well known problem.
These steps allow you to obtain the density of the nearest neighbour distance from $x_1$, and hence the moments. Of course, all this is not very straighforward to compute, but I doubt that there is a shorcut (to obtain the first moments only). Perhaps for large $n$ and/or $k$ some asymptotics can be found.
See also this related question
Update: If the number of points $k$ is big (and the dimension $n$ is not, so that that the selected point is, with high probabily, surrounded by other points, and not too near the box edges; specifically, if $k > M^n$ with  $M$ , say, greater than 10 or 30) we can get an approximation by computing the volume of the "diamond" shape that corresponds to a bounded Manhattan distance (cross polytope). So, if $X$ is the Manhattan distance of the nearest neighbour, we'd get
$$P(X \le x) \approx 1 - \left(1 - x^n \frac{2^n}{n!}\right)^{k-1}, \;\; 0\le x\le x_M=\frac{(n!)^{1/n}}{2}$$
And so, for example, we could approximate $$E(X) = \int_0^\infty (1-F_X(x))\, dx \approx \int_0^{x_M}  \left(1 - x^n \frac{2^n}{n!}\right)^{k-1} dx$$
