Finding the Expected Average Distance from the Maximum given a distribution For a given sample set $S$ with $N$ individual samples $x_i$, I can easily find the average distance from the maximum by doing something like this:
$\sigma_{max_N}:=\sqrt{\frac{1}{N}\sum\limits_{i=1}^N {\left(x_i - \max{\left(S\right)}\right)}^2}$
Trying to do so for larger and larger sample sizes, I (unsurprisingly) end up with different limit values given different distributions.
For instance, for a uniform distribution on a unit interval, I get a value around 0.577 with digits after that varying wildly.
A normal distribution with the same mean $\left(\frac{1}{2}\right)$ and variance $\left(\frac{1}{12}\right)$ gives me something around 1.5, though the exact value shifts around a lot more.
For an exponential distribution with parameter $\lambda=2\sqrt{3}$, I end up at around 4.5.
Is there a way to get an expected value of $\sigma_{max_N} \to \sigma_{max}$ for a given distribution as the sample size $N \to \infty$? It doesn't necessarily have to be analytic.
 A: I presume we are talking about a collection of i.i.d random variables, with finite moments. Using a bit more convenient notation, let's consider (we can afterwards apply the continuous mapping theorem for the square root),
$$\sigma^2_{max_n}:=\frac{1}{n}\sum\limits_{i=1}^n {\left(X_i - X_{(n)}\right)}^2$$
$$= \frac{1}{n}\sum\limits_{i=1}^n X_i^2 -2X_{(n)}\bar X+X_{(n)}^2$$
where $X_{(n)}$ is the maximum order statistic, and $\bar X$ is the sample mean. We want the probability limit of this expression, which will be, if it exists, the sum of the probability limits of the three components.  
We have easily that 
$$\frac{1}{n}\sum\limits_{i=1}^n X_i^2 \xrightarrow{p} E(X_1^2)$$
but the main question is 

Does the maximum order statistic converges in probability
  somewhere?

(while the usual issue is whether extreme order statistics, properly transformed, converge in distribution - to some variant of the Extreme Value Distribution).
Applying the definition of convergence in probability, and searching for the value of $k$, we would want that $\exists \;k:$
$$\begin{align}P[|X_{(n)}-k|< \epsilon]=1 &\Rightarrow P[k-\epsilon<X_{(n)}<k+\epsilon]=1\\ &\Rightarrow  F_{Y_{(n)}}(k+\epsilon)- F_{Y_{(n)}}(k-\epsilon)=1 \\&\Rightarrow [F_X(k+\epsilon)]^n-[F_X(k-\epsilon)]^n =1 \end{align}$$
since $F_{Y_{(n)}}(y) = [F_X(y)]^n$.
RANDOM VARIABLES with support bounded from above 
We see what happens here: for random variables with bounded support from above, denote this upper bound simply $\theta$, if we set $k=\theta$ then, by the properties of the cumulative distribution function, the equality above will hold, and so 
$$X_{(n)} \xrightarrow{p} \theta \Rightarrow X_{(n)}^2 \xrightarrow{p} \theta^2$$
by the continuous mapping theorem. The sample mean will converge to the expected value so in all, for this class of random variables,
$$\sigma^2_{max_n} \xrightarrow{p} E(X_1^2) -2\theta E(X_1) + \theta^2$$
Verification: For a sample of i.i.d uniforms $U(0,1)$ we have $E(X_1^2) =1/3$, $E(X_1) =1/2$, $\theta=1$.
So (using again the continuous mapping theorem)
$$\sigma_{max_n}  \xrightarrow{p} \sqrt{\frac 13 -2\cdot1\cdot \frac 12 +1^2} = \sqrt{1/3} \approx 0.577$$
as the OP found.  
RANDOM VARIABLES with support unbounded above
Intuitively, here the maximum order statistic will tend to $+\infty$ as $n\rightarrow \infty$. So we have to accept approximations, since our sample will be finite after all. Treating the continuous case the OP treated, for a normal $N(0.5, 1/12)$ a virtual "zero-tail probability" value is   $\theta=2$. Moreover, 
$$E(X_1^2)= \operatorname{Var}(X) + \mu^2 = 1/12 + (1/2)^2 = 1/3$$. So here
$$\sigma_{max_n}  \approx \sqrt{\frac 13 -2\cdot2\cdot \frac 12 +2^2} =\sqrt{2.33}\approx 1.52$$
This matches the OP simulations, but of course it is very sensitive to the exact tail probability chosen. Why not choose $\theta=3$?... This is not a robust approximation. Moreover we are silent about the size of $n$ -although in finite samples, we will never observe an actual infinity, with really large $n$, well, there is a whole world in the interval $(2,\infty)$. So again, the above value may match the specific batch of simulations, but see @whuber's comments below for a way out (and also, the OP appears to have found a way...). In any case, this valuet does not represent an upper bound on the $\sigma_{max_n}$.
A: This alternative is easier to analyze:
$$
\newcommand{\E}{\mathbb{E}}
d(S) := \frac{1}{N} \sum_{i=1}^N (\max(S) - x_i)
$$
Then we have
$$\begin{align*}
\E d(S)
&= \frac{1}{N} \sum_{i=1}^N (\E \max(S) - \E x_i)
= \E \max(S) - \E x
\end{align*}$$
$\E \max(S)$ might or might not be easy to find for a given distribution.


*

*For uniform distribution on $[a, b]$, it of course converges quickly to $b$; you can find an explicit expression without too much difficulty.

*For normals, you can get a good approximation.

*For exponentials, it's the $n$th harmonic number over the rate: basic argument, more rigorous.

If you want to stick to the root-mean-square type distance, I'll call it $d_2$ because I think $\sigma_{\max_N}$ is misleading:
$$\begin{align*}
d_2^2
&= \frac{1}{N} \sum_{i=1}^N (x_i - \max(S))^2
\\&= \frac{1}{N} \sum_{i=1}^N \left( x_i^2 + \max(S)^2 - 2 x_i \max(S) \right)
\\&= \frac{1}{N} \sum_{i=1}^N x_i^2 + \max(S)^2 - 2 \frac{1}{N} \sum_{i=1}^N x_i \max(S)
\end{align*}$$
Then
$$\begin{align*}
\E d_2^2
&= \E x^2 + \E( \max(S)^2 ) - 2 \frac{1}{N} \sum_{i=1}^N \E \left( x_i \max(S) \right)
\end{align*}$$
Now,
$$\begin{align*}
\E x_i \max(S)
&= \E\left[ x_i \max(S) \mid x_i = \max(S) \right] \Pr(x_i = \max(S))
\\&\qquad   + \E\left[ x_i \max(S) \mid x_i \ne \max(S) \right] \Pr(x_i \ne \max(S))
\\&= \E\left[ \max(S)^2 \right] \frac{1}{N}
   + \E\left[ x_i \max(S) \mid x_i \ne \max(S) \right] \frac{N-1}{N}
\\&\approx \frac{1}{N} \E\left[ \max(S)^2 \right]
   + \frac{N-1}{N} \E\left[ x_i \right] \E\left[ \max(S) \right]
\end{align*}$$
when $N$ is large.
Then
$$\begin{align*}
\E d_2^2
&\approx \E x^2 + \E( \max(S)^2 ) - 2 \frac{1}{N} \E\left[ \max(S)^2 \right] - 2 \frac{N-1}{N} \E\left[ x_i \right] \E\left[ \max(S) \right]
\\&= \E x^2 + \frac{N - 2}{N} \E\left[ \max(S)^2 \right] - 2 \frac{N-1}{N} \E\left[ x_i \right] \E\left[ \max(S) \right]
.
\end{align*}$$
Via Jensen's inequality, $\E d_2 = \E \sqrt{d_2^2} \le \sqrt{ \E d_2^2}$.
So you can at least get an approximate upper bound from the first and second moments of $x$ and $\max S$.
