Precise asymptotics for moments of order statistics of normal distribution Let $X_1, \cdots, X_n \sim N(0,1)$ be i.i.d. normal random variates.  I am interested in understanding the first two moments of the quasi-range $X_{(n)}-X_{(n-1)}$ (i.e., the maximum value minus the second-largest value) as $n \to \infty$.
I am trying to use the quantile-function method of section 4.6 of David and Nagaraja's Order Statistics, but am having trouble with convergence of some approximations.  The method presented there (originally due to Pearson, later expanded by David and Johnson) is to Taylor-expand the relation $X_{(r)} = Q(U_{(r)})$ about $\mathbb{E}[U_{(r)}]$, where $Q(u) = F^{-1}(u)$ is the quantile function (the inverse of the cdf) and $U_{(r)}$ are the rather-more-well-understood order statistics of the uniform $U[0,1]$ distribution.
The Taylor expansion gives us (Equation 4.6.2), writing $p_r = r/(n+1)$ and $q_r = 1-p_r$, the series $$X_{(r)} = Q(p_r) + \left(U_{(r)}-p_r\right)Q'(p_r) + \frac{1}{2}\left(U_{(r)}-p_r\right)^2Q''(p_r) + \frac{1}{6}\left(U_{(r)}-p_r\right)^2Q'''(p_r) + \cdots.$$
From here, we can use formulas for uniform order statistics to get the first two moments of $X_{(r)}$, as well as covariances $\mathrm{Cov}[X_{(r)},X_{(s)}]$; these are Equations 4.6.3-4.6.5 of David and Nagaraja.  For reference, I give the equation for the expected value: $$\mathbb{E}[X_{(r)}] = Q(p_r) + \frac{p_rq_r}{2(n+2)}Q''(p_r) + \frac{p_rq_r}{(n+2)^2}\left[\frac{1}{3}(q_r-p_r)Q'''(p_r)+\frac{1}{8}p_rq_rQ''''(p_r)\right]+O\left(\frac{1}{(n+2)^3}\right).$$
We can use a neat trick (cf. Example 4.6 of ibid.) to calculate the derivatives of the quantile function for a normal distribution; this boils down to the fact that the pdf $f(x) = (2\pi)^{-1/2}\exp(-x^2/2)$ of a standard normal distribution satisfies $f'(x) = -x\cdot f(x)$.  This gives us the general relation $$Q^{(k)}(u) = \frac{P_k(Q(u))}{f(Q(u))^k},$$ where $P(x)$ is a family of polynomials satisfying $P_1(x) = 1$ and $P_{k+1}(x) = kxP_k(x) + P'_k(x)$.  Importantly, to leading order $P_{k+1}(x) = k!x^k + O(x^{k-1})$.
We can then use the standard approximations $$Q(1-u) \approx \sqrt{\log\left(\frac{1}{2\pi u^2}\right)}, \quad f(Q(1-u)) \approx u\sqrt{\log\left(\frac{1}{2\pi u^2}\right)}, $$ which are valid as $u \to 0$.
So if we want to compute the moments of $X_{(n)}$, say, we would take $p_r = n/(n+1)$ and $q_r = 1/(n+1)$.  The above terms of the expected-value equation then become, writing for convenience $s = \sqrt{\log\left(\frac{(n+1)^2}{2\pi}\right)}$, $$\mathbb{E}[X_{(n)}] \approx s + \frac{1}{s}\left(\frac{n}{2(n+2)} + \frac{2n(n-1)}{3(n+2)^2} + \frac{3n^2}{4(n+2)^2} + \cdots\right) + O\left(\frac{1}{s^2}\right).$$
This is fine for estimating $\mathbb{E}[X_{(n)}] \approx s$, but if we want to look at the expected difference $\mathbb{E}[X_{(n)}-X_{(n-1)}]$, we need to understand the $1/s$ term, because the difference between the $s$-terms in the $X_{(n)}$ and $X_{(n-1)}$ term is on the order of $1/s$.
But it's not clear that the coefficient on the $1/s$ term in the series for $\mathbb{E}[X_{(n)}]$ even converges!  Indeed, David and Nagaraja warn,

From a practical point of view, the most important feature of the expansion is that convergence may be slow or even nonexistent if $r/n$ is too close to 0 or 1.

For the moments of $X_{(n-1)}$, I think that we get coefficients that are geometrically decreasing, so things are fine.  But the same convergence problem plagues our series for $\mathrm{Var}[X_{(n)}]$ and possibly $\mathrm{Cov}[X_{(n)},X_{(n-1)}]$, though I haven't checked the latter series thoroughly.
Question: is there another way to get asymptotic expressions for the first two moments of the normal order statistics $X_{(n)}$ and $X_{(n-1)}$ (as well as their covariance) that has better convergence properties?
Alternate suggestions for understanding the first two moments (or the entire distribution) of the quasi-range $X_{(n)}-X_{(n-1)}$ would also be helpful.  It might be possible to use ad-hoc methods (e.g., Fisher–Tippett–Gnedenko) for $X_{(n)}$ and then the series above for $X_{(n-1)}$, but that would rely on the covariance series converging quickly, which I am not yet confident in.
 A: I call $f$ the density and $F$ the cumulative distribution function of $\mathcal{N}(0,1)$.
Since $X_{(n)} \ge X_{(n-1)}$,
$$X_{(n)}-X_{(n-1)} = \int_\mathbb{R} 1_{[X_{(n-1)} \le x < X_{(n)}]}dx.$$
Taking expectations and applying Fubini theorem for positive functions yields
$$E[X_{(n)}-X_{(n-1)}] = \int_\mathbb{R} P[X_{(n-1)} \le x < X_{(n)}]dx.$$
$$E[X_{(n)}-X_{(n-1)}] = \int_\mathbb{R} nF(x)^{n-1}(1-F(x))dx.$$
More generally, for every integer $d \ge 1$,
$$(X_{(n)}-X_{(n-1)})^d = \int_{\mathbb{R}^d} 1_{[X_{(n-1)} \le x_1 < X_{(n)}]} \times \cdots \times 1_{[X_{(n-1)} \le x_d < X_{(n)}]}dx_1 \cdots dx_d.$$
$$(X_{(n)}-X_{(n-1)})^d = \int_{\mathbb{R}^d} 1_{[X_{(n-1)} \le \min(x_1,\ldots,x_d) \le \max(x_1,\ldots,x_d) < X_{(n)}]} dx_1 \cdots dx_d.$$
Taking expectations and applying Fubini theorem for positive functions yields
$$E[(X_{(n)}-X_{(n-1)})^d] = \int_{\mathbb{R}^d} P[X_{(n-1)} \le \min(x_1,\ldots,x_d) \le \max(x_1,\ldots,x_d) < X_{(n)}] dx_1 \cdots dx_d.$$
$$E[(X_{(n)}-X_{(n-1)})^d] = \int_{\mathbb{R}^d} nF(\min(x_1,\ldots,x_d))^{n-1} \big( 1 - F (\max(x_1,\ldots,x_d)) \big) dx_1 \cdots dx_d.$$
If we want the asymptotic behaviour, we may observe that the function $b$ defined by
$$b(x) := F(x)(1-F(x))/f(x)$$ is in $\mathcal{C}_0(\mathbb{R})$ and equivalent to $1/|x|$ as $|x| \to +\infty$. One can write
$$E[X_{(n)}-X_{(n-1)}] 
= \int_\mathbb{R} nF(x)^{n-2}b(x)f(x)dx 
= \int_0^1 ny^{n-2}b(F^{-1}(y))dy.$$
We get $n/(n-1)$ times the mean value of $b \circ F^{-1}$ for the distribution Bêta$(1,n-1)$, which converges to $\delta_1$ as $n \to +\infty$.
Using that $b(F^{-1}(y)) \to 0$ slowly as $y \to 1$, we will get that $E[X_{(n)}-X_{(n-1)}] \to 0$ slowly.
