Suppose $X_1,\dots,X_n$ are independent samples from some distribution with known absolutely continuous CDF $F:\mathbb{R}\rightarrow[0,1]$. Let $X_{(1)},\dots,X_{(n)}$ denote the order statistics, i.e. the ordered sample. Defining the column vector $X_{(\cdot)}=[X_{(1)},\dots,X_{(n)}]'$, we want to numerically calculate $\mathbb{E}[(X_{(\cdot)}-\mathbb{E}X_{(\cdot)})(X_{(\cdot)}-\mathbb{E}X_{(\cdot)})']$.
The most obvious approach uses the fact that for $j,k\in\{1,\dots,n\}$, $j<k$:
$$\mathbb{E}X_{(k)}=\int{x\frac{n! [F(x)]^{k-1}[1-F(x)]^{n-k}}{(k-1)!(n-k)!} dF(x)},$$ $$\mathbb{E}X_{(k)}^2=\int{x^2\frac{n! [F(x)]^{k-1}[1-F(x)]^{n-k}}{(k-1)!(n-k)!} dF(x)},$$ $$\mathbb{E}X_{(j)}X_{(k)}=\int{\int{xy\frac{n! [F(x)]^{j-1}[F(y)-F(x)]^{k-1-j}[1-F(y)]^{n-k}}{(j-1)!(k-j-1)!(n-k)!} 1[x\le y] dF(x) } dF(y)},$$
where $1[\cdot]$ is the indicator function. See e.g. Wikipedia here for an informal proof.
Using these formulae with standard numerical integration methods works well for small $n$. However, for large $n$ the resulting covariance matrix often ends up non-positive definite unless implausibly many integration nodes are used. This is despite the individual elements of the covariance matrix usually being (loosely) close to a covariance matrix generated via a naïve Monte Carlo approach that guarantees positive definiteness (i.e. draw such a sample of length $n$, then sort it, repeat this lots of times, take the covariance).
Is it possible to express these integrals in such a way that the resulting covariance matrix is guaranteed to be positive definite?
E.g. is it the case that:
$$\mathbb{E}[(X_{(\cdot)}-\mathbb{E}X_{(\cdot)})(X_{(\cdot)}-\mathbb{E}X_{(\cdot)})']=\int{\int{g(u,v) dF(u)}dF(v)},$$
for some function $g:\mathbb{R}^2\rightarrow \mathbb{R}^{n\times n}$ where $g(u,v)$ is positive semi-definite for all $u,v\in\mathbb{R}$.