Expectation of maximum of multivariate Gaussian

Given a multivariate Gaussian $$\mathbf{X} \sim \mathcal{N}(\mathbf{\mu},\Sigma)$$, I believe it is a difficult question to find a closed form for $$\mathbb{E}[ \max\{X_1,\ldots,X_d\}].$$

However, the case I have at hand is perhaps combinatorially nicer: my Gaussian vector is $$(X_1,\ldots,X_n)$$ where $$X_j = Z_j - \frac{\sum_{i \neq j} Z_i}{n-1}$$

where $$Z_1,\ldots,Z_n$$ are i.i.d. standard normals.

In the $$n = 2$$ case, the vector $$(X_1,X_2)$$ is distributed like $$\frac{1}{\sqrt{2}}(Z,-Z)$$ where $$Z$$ is a standard normal, and thus the maximum is simply $$\frac{1}{\sqrt{2}} \mathbb{E}[|Z|] = \frac{1}{\sqrt{ \pi}}.$$

Can anything like this be done to compute $$c_n = \mathbb{E}\left[\max\{X_1,\ldots,X_n\} \right]$$

for general $$n$$? Note also that this is the same (up to scaling) of taking $$(X_1,\ldots,X_n)$$ to be i.i.d. normals conditioned on summing to $$0$$. It wouldn't be too surprising to me if this is possible to compute precisely (and someone has already done so!)

[this is cross-posted from math stack exchange]

For $$S:=\sum_1^n Z_j$$, we have $$X_j=Z_j-\frac{S-Z_j}{n-1}=-\frac{S}{n-1}+\frac{n}{n-1}\,Z_j,$$ whence $$\max_1^n X_j=-\frac{S}{n-1}+\frac{n}{n-1}\,\max_1^n Z_j$$ and $$E\max_1^n X_j=\frac{n}{n-1}\,EM_n,\quad M_n:=\max_1^n Z_j.$$ In turn, $$EM_n=\int_0^\infty [P(M_n>x)-P(M_n<-x)]\,dx =\int_0^\infty [1-\Phi(x)^n-\Phi(-x)^n]\,dx,$$ where $$\Phi$$ is the standard normal cdf. Alternatively, we can write $$EM_n=\int_{-\infty}^\infty x\, dP(M_n cf. Robert Israel's answer.