The probability density function of the distance between two points chosen randomly on the unit square is given by:
$ P(\ell) = \begin{cases} 2\ell\left(\ell^2 - 4\ell + \pi\right) & 0 \leq \ell \leq 1 \\ 2\ell\left[4\sqrt{\ell^2 - 1} - \left(\ell^2 + 2 - \pi\right) - 4\arctan\sqrt{\ell^2 - 1} \right] & 1 \leq \ell \leq \sqrt{2} \end{cases} $
While this formula checks out empirically, unfortunately, only a personal communication is provided as reference for it.
I would be grateful if someone on here could help me understand how it was derived.
I have the beginnings of what may be a solution:
Let $U = (U_1,U_2)$ and $W=(W_1,W_2)$ be two independent rv's with uniform distributions on the unit square, so $U_i,W_i \sim f(x) := I_{[0,1]}(x)$. Let $D := ((U_1-W_1)^2 + (U_2-W_2)^2)^{1/2}$ be the distance between the two random points. First suppose $d \in [1,\sqrt{2}]$. Write \begin{align} P[D \leq d \mid W] &= \int_{(u_1,u_2) \in B_d(W_1,W_2)} du_1 du_2 \\ &= \int_0^1 \int_0^{W_1 + \sqrt{d^2 - (u_2-W_2)^2}} du_1 du_2 \\ &= \int_0^1 (W_1 + \sqrt{d^2 - (u_2-W_2)^2}) du_2 \\ &= W_1 + \int_0^1 \sqrt{d^2 - (u_2-W_2)^2}) du_2. \end{align} For the last integral on the right, Wolfram alpha gives the anti-derivative $$ \frac{1}{2}\left((u_2-W_2)\sqrt{d^2-(W_2-u_2)^2}+d^2 \arctan\left(\frac{u_2-W_2}{\sqrt{d^2-(W_2-u_2)^2}}\right)\right), $$ and so if $g_1(w) = \sqrt{d^2 - (w-1)^2}$ and $g_0(w) = \sqrt{d^2 - w^2}$ one finds $$ P[D \leq d \mid W] = W_1 + \frac{1}{2}\left( (1-W_2)g_1(W_2)+d^2 \arctan\left(\frac{1-W_2}{g_1(W_2)}\right)\right) -\frac{1}{2}\left( -W_2g_0(W_2)+d^2 \arctan\left(\frac{-W_2}{g_0(W_2)}\right) \right). $$ Then the distribution function for $D$ is \begin{align} P[D \leq d] &= {\bf E}[P[D \leq d \mid W]]. \end{align} Now, \begin{align} {\bf E}\left[ (1-W_2)g_1(W_2) + W_2 g_0(W_2) \right] &= \int_0^1 ((1-w)\sqrt{d^2-(1-w)^2} + w\sqrt{d^2-w^2} )dw \\ &= 2 \int_0^1 w\sqrt{d^2-w^2} dw \\ &= \left. -\frac{2}{3}(d^2-w^2)^{3/2} \right|^1_0 \\ &= \frac{2}{3}(d^3 - (d^2 - 1)^{3/2}), \end{align} and, since \begin{align} {\bf E}\left[\arctan\left(\frac{1-W_2}{g_1(W_2)}\right) \right] &= \int_0^1 \arctan\left(\frac{1-w}{\sqrt{d^2-(1-w)^2}}\right) dw \\ &= \int_0^1 \arctan\left(\frac{y}{\sqrt{d^2-y^2}}\right) dy \\ &= {\bf E}\left[-\arctan\left(\frac{-W_2}{g_0(W_2)}\right) \right], \end{align} one has \begin{align} {\bf E}\left[\arctan\left(\frac{1-W_2}{g_1(W_2)}\right) - \arctan\left(\frac{-W_2}{g_0(W_2)}\right) \right] &= 2\int_0^1 \arctan\left(\frac{y}{\sqrt{d^2-y^2}}\right) dy \\ &= \left. 2\left(\sqrt{d^2 - w^2} + w\arctan\left(\frac{w}{\sqrt{d^2-w^2}}\right) \right) \right|^1_0 \\ &= 2\left(\sqrt{d^2 - 1} + \arctan\left(\frac{1}{\sqrt{d^2-1}}\right) - d \right). \end{align} Therefore \begin{align} P[D \leq d] &= \frac{1}{2} + \frac{1}{3}(d^3 - (d^2 - 1)^{3/2}) + d^2\left(\sqrt{d^2 - 1} + \arctan\left(\frac{1}{\sqrt{d^2-1}}\right) - d \right) \\ &= \frac{1}{2} + \frac{-2d^3 - (d^2 - 1)^{3/2}+3d^2\sqrt{d^2-1}}{3} + d^2 \arctan\left( \frac{1}{\sqrt{d^2-1}}\right) \\ &= \frac{1}{2} + \frac{-2d^3 +(2d^2 +1)\sqrt{d^2-1}}{3} +d^2 \arctan\left( \frac{1}{\sqrt{d^2-1}}\right). \end{align} Differentiating with respect to $d$ should yield the density for $D$: \begin{align} f_D(d) &= -2d^2 + \frac{d(2d^2-1)}{\sqrt{d^2-1}}+2d\arctan\left(\frac{1}{\sqrt{d^2-1}} \right) - \frac{d}{\sqrt{d^2-1}} \\ &= -2d^2 + \frac{2d(d^2-1)}{\sqrt{d^2-1}}+2d\arctan\left(\frac{1}{\sqrt{d^2-1}} \right) \\ &= 2d\left( -d + \sqrt{d^2-1}+\arctan\left(\frac{1}{\sqrt{d^2-1}} \right)\right), \end{align} $1\leq d \leq \sqrt{2}$. Now use the identity $\arctan(1/x) = -\frac{1}{2}\pi - \arctan(x)$ to get
\begin{align} f_D(d) &= 2d\left(\sqrt{d^2-1} -d -\frac{1}{2}\pi-\arctan(\sqrt{d^2-1}) \right), \quad 1 \leq d \leq \sqrt{2}. \end{align} As you can see, this gives a slightly different density in the region $[1,\sqrt{2}]$, but there may be some other trig identity that makes them equal (I really doubt it). Is this helpful? Sorry I couldn't get to the density you have.