Consider going a small distance $h$ along the diagonal away from a corner. The probability of finding a point in the corner triangle with height $h$ is $n h^2$, and such a point is at least $x = \sqrt{2} - 2h$ away from a point in the same triangle at the opposite corner. If $D$ is the maximum distance, we have for small enough $h$ $$\Pr[D > \sqrt{2} - 2h] \approx n h^2 .$$ Rewriting $h = (\sqrt{2}-x)/2$, we hence get the approximate distribution function $$F(x) = \Pr[D \le x] \approx 1 - \frac{n}{4} (\sqrt{2} - x)^2$$ with the lowest acceptable $x_0$ where $F(x_0) = 0$ or $x_0 = \sqrt{2} - 2/\sqrt{n}$. Integrating then gives the expectation, $$E[D] \approx \int_{x_0}^{\sqrt{2}} \! x \, F'(x) \, dx = \sqrt{2} - \frac{4}{3\sqrt{n}}.$$
The idea here is that there is a density $n_1$ of points for which the maximum distance is a.s. larger than unity. In that case, only points in opposing corners need to be considered. From the approximation above, $x_0 = 1$ implies $n_1 \approx 12 + 8\sqrt{2} \approx 23.3$, so that is where we can expect the result to become reasonable.
Here is a plot of the approximation and simulated distances, which seems to indicate that this is more or less what's actually going on:
Of course, there is then the further approximation of the true distance by the projection onto the diagonal. But the more we move into the corners, the better that approximation becomes.