Let $X \sim \mathcal N(0, \Sigma)$ where we take $0$ to be the mean w.l.o.g., and where $\Sigma$ is a symmetric positive-definite covariance matrix of dimension $d\times d$. For $d = 1$, it is known that it holds $P(|X| \leq 1.96\sqrt{\Sigma}) = 0.95$, so that $1.96\sqrt{\Sigma}$ is the reference value for computing a 95% confidence interval. Let now $d > 1$. What should I employ, in this case, for a confidence interval? The square root $\sqrt{\Sigma}$, such that $\sqrt{\Sigma}(\sqrt{\Sigma})^\top = \Sigma$ is not unique in this case. Can the Cholesky factor $L \in \mathbb R^{d \times d}$, such that $LL^\top = \Sigma$, be employed?
The application is the following. Imagine I have a $d$-dimensional random variable $X \sim \mathcal N(m, \Sigma)$, with known $m$ and $\Sigma$. I would like to plot the mean $m$ (simply against $(1,\cdots, d)$) and a "confidence interval" around the mean $m$, i.e., a "shaded" area around $m$, where $X$ belongs with high probability (e.g. 95%). I've seen this done by taking $m \pm 1.96\sqrt{\mathrm{diag}(\Sigma)}$, but this seems to me arbitrary, because one neglects completely correlations. Would taking the diagonal of the Cholesky factor be better? If not, which of the square roots of $\Sigma$ would be the go-to choice?
By eigendecomposition $\Sigma=W\Lambda W^{-1}$ you obtain the principal components (PCs) $V=WX$. Define the standardized PCs $\tilde{V}=V\Lambda^{-\frac{1}{2}}$. These are $\sim \mathcal{N}(0,1)$ and uncorrelated. We want to know for a given $\alpha \in (0,1)$ $$P(\tilde{V} \in B_{r}(0))=\alpha$$ where $B_{r}(0)$ is a $d$-ball of radius $r$ (our unknown). This is equivalent to $$P(\|\tilde{V}\|^2<r^2)=\alpha$$ Recall that the sum of $d$ uncorrelated standard normal rvs is $\sim\chi^2_d$ so we have $$r^2=(\chi_{d}^2)^{-1}(\alpha)$$ The confidence region contour is recovered by stretching the $\bar{B}_{r}(0)$ along the axes using the square root of the eigenvalues and in the direction of the eigenvectors, i.e. for any $(x_1,x_2,...,x_d)\in \bar{B}_{r}(0)$ we apply $$W\begin{bmatrix} x_1\sqrt{\lambda_1}\\ x_2\sqrt{\lambda_2}\\ ...\\ x_d\sqrt{\lambda_d} \end{bmatrix}$$