(Too long for a comment)

I don't have good knowledge about literature on this problem, but here are some of my thoughts:

1) For the "typical" distances, the set of $n$ uniformly chosen points is not so different from $n$ points forming a lattice in the unit square. So we should expect that the profile of $t\mapsto d_{k(t)}$ for $k(t) = \bigl\lfloor\binom{n}{2}t\bigr\rfloor$ converges, in an appropriate sense, to the inverse CDF of the distance between two randomly chosen points in the unit square. Numerical simulations indeed support this picture:

I suspect that a kind of coarse-graining argument might help verify this, although this approach would be too crude for obtaining good concentration results.

2) The small values of $d_k$ will be achieved only for occasional pairs that are close to each other, and we would expect that such pairs occur almost independently of each other. More precisely, fix $\ell \gg 1$ and let $m \sim n/\ell $ as $n\to\infty$. Also, let $\mathsf{B}_{i,j} = [\frac{i-1}{m}, \frac{i}{m}]\times[\frac{j-1}{m}, \frac{j}{m}]$ be subsquares of $[0, 1]^2$. Then the event $A_{i,j}$ that there are at least two points lying in $\mathsf{B}_{i,j}$ is approximately $n^2/2m^4$, so there are typically $n^2/2m^2 \sim \ell^2/2$ number of subsquares containing two or more points.

This suggests that $n d_k$'s for "small" $k$'s will behave like the distance between the origin and the $k$-th closest point in the Poisson point process on $\mathbb{R}^2$ with intensity $\frac{1}{2}$, or equivalently, the arrivals in the inhomogeneous Poisson process on $[0, \infty)$ with intensity measure $\lambda(\mathrm{d}t) = \mathrm{d}(\pi t^2/2)$. A numerical simulation also confirms this heuristics:

(The blue dots are simulated values of $nd_k$'s for $k=1,2,\ldots,200$ with $n=1000$, and the orange line is the graph of the function $k\mapsto\sqrt{2k/\pi}$.)

(Too long for a comment)

I don't have good knowledge about literature on this problem, but here are some of my thoughts:

1) For the "typical" distances, the set of $n$ uniformly chosen points is not so different from $n$ points forming a lattice in the unit square. So we should expect that the profile of $t\mapsto d_{k(t)}$ for $k(t) = \bigl\lfloor\binom{n}{2}t\bigr\rfloor$ converges, in an appropriate sense, to the inverse CDF of the distance between two randomly chosen points in the unit square. Numerical simulations indeed support this picture:

I suspect that a kind of coarse-graining argument might help verify this, although this approach would be too crude for obtaining good concentration results.

2) The small values of $d_k$ will be achieved only for occasional pairs that are within $\mathcal{O}(1/n)$ distance to each other, and we would expect that such pairs occur almost independently of each other. More precisely, fix $\ell \gg 1$ and let $m \sim n/\ell $ as $n\to\infty$. Also, let $\mathsf{B}_{i,j} = [\frac{i-1}{m}, \frac{i}{m}]\times[\frac{j-1}{m}, \frac{j}{m}]$ be subsquares of $[0, 1]^2$. Then the event $A_{i,j}$ that there are at least two points lying in $\mathsf{B}_{i,j}$ is approximately $n^2/2m^4$, so there are typically $n^2/2m^2 \sim \ell^2/2$ number of subsquares containing two or more points.

This suggests that $n d_k$'s for "small" $k$'s will behave like the distance between the origin and the $k$-th closest point in the Poisson point process on $\mathbb{R}^2$ with intensity $\frac{1}{2}$, or equivalently, the arrivals in the inhomogeneous Poisson process on $[0, \infty)$ with intensity measure $\lambda(\mathrm{d}t) = \mathrm{d}(\pi t^2/2)$. A numerical simulation also confirms this heuristics:

(The blue dots are simulated values of $nd_k$'s for $k=1,2,\ldots,200$ with $n=1000$, and the orange line is the graph of the function $k\mapsto\sqrt{2k/\pi}$.)