[proofplan]
We use the [one-sample Hotelling confidence ellipsoid](/theorems/4026) for $\mu$ and show that every linear functional $m \mapsto a^\top m$ is bounded on that ellipsoid by an explicit support-function calculation. On the Hotelling event, the displacement $\mu-\bar X$ has controlled length in the $S^{-1}$-inner product. A finite-dimensional Cauchy-Schwarz argument in this inner product gives the bound for each fixed $a$, and because the same Hotelling event is used for all $a \in \mathbb{R}^p$, the resulting intervals are simultaneous.
[/proofplan]
[step:Define the Hotelling ellipsoid event with its exact coverage]
Define the critical constant $c_\alpha \in (0,\infty)$ by
\begin{align*}
c_\alpha
&:= \frac{p(n-1)}{n-p}F_{p,n-p;1-\alpha}.
\end{align*}
Let $E \in \mathcal{F}$ be the event
\begin{align*}
E
:=
\left\{
S \text{ is positive definite and }
n(\bar X-\mu)^\top S^{-1}(\bar X-\mu) \le c_\alpha
\right\}.
\end{align*}
By the one-sample Hotelling confidence ellipsoid theorem (citing a result not yet in the wiki: One-Sample Hotelling Confidence Ellipsoid Theorem), applied to the independent sample $X_1,\dots,X_n$ from $\mathcal{N}_p(\mu,\Sigma)$ with $n>p$ and positive definite $\Sigma$, we have
\begin{align*}
\mathbb{P}(E)=1-\alpha.
\end{align*}
[guided]
The role of Hotelling's theorem is to produce one event, not one event for each direction $a$. Define
\begin{align*}
c_\alpha
&:= \frac{p(n-1)}{n-p}F_{p,n-p;1-\alpha}.
\end{align*}
The Hotelling confidence ellipsoid event is
\begin{align*}
E
:=
\left\{
S \text{ is positive definite and }
n(\bar X-\mu)^\top S^{-1}(\bar X-\mu) \le c_\alpha
\right\}.
\end{align*}
The hypotheses needed for the one-sample Hotelling confidence ellipsoid theorem are exactly the hypotheses in the statement: $X_1,\dots,X_n$ are independent and identically distributed as $\mathcal{N}_p(\mu,\Sigma)$, the covariance matrix $\Sigma$ is positive definite, and $n>p$. Therefore that theorem gives
\begin{align*}
\mathbb{P}(E)=1-\alpha.
\end{align*}
Everything that follows is deterministic on this single event $E$. This is what makes the final confidence intervals simultaneous in $a$.
[/guided]
[/step]
[step:Control every linear contrast on the ellipsoid]
Fix an outcome $\omega \in E$. Define $h \in \mathbb{R}^p$ by
\begin{align*}
h := \mu-\bar X(\omega).
\end{align*}
Since $S(\omega)$ is positive definite, define the bilinear map
\begin{align*}
(\cdot,\cdot)_{S(\omega)} : \mathbb{R}^p \times \mathbb{R}^p &\to \mathbb{R} \\
(u,v) &\mapsto u^\top S(\omega)^{-1}v .
\end{align*}
This is an inner product on $\mathbb{R}^p$.
We record the needed Cauchy-Schwarz estimate in this inner product. For $u,v \in \mathbb{R}^p$, the non-negativity of
\begin{align*}
(u-tu,u-tu)_{S(\omega)}
\end{align*}
is not the useful quadratic; instead take $u-tv$ with $t \in \mathbb{R}$. Then
\begin{align*}
0
&\le (u-tv,u-tv)_{S(\omega)} \\
&= (u,u)_{S(\omega)} -2t(u,v)_{S(\omega)} + t^2(v,v)_{S(\omega)}.
\end{align*}
If $v \ne 0$, this quadratic in $t$ has non-positive discriminant, so
\begin{align*}
|(u,v)_{S(\omega)}|^2
\le
(u,u)_{S(\omega)}(v,v)_{S(\omega)}.
\end{align*}
The case $v=0$ is immediate.
Now fix $a \in \mathbb{R}^p$. Since $S(\omega)$ is symmetric, we have
\begin{align*}
a^\top h
&= h^\top a \\
&= h^\top S(\omega)^{-1}S(\omega)a \\
&= (h,S(\omega)a)_{S(\omega)}.
\end{align*}
Applying the estimate above with $u=h$ and $v=S(\omega)a$ gives
\begin{align*}
|a^\top h|^2
&\le
(h,h)_{S(\omega)}(S(\omega)a,S(\omega)a)_{S(\omega)} \\
&=
h^\top S(\omega)^{-1}h \, a^\top S(\omega)a.
\end{align*}
Because $\omega \in E$,
\begin{align*}
h^\top S(\omega)^{-1}h
\le
\frac{c_\alpha}{n}.
\end{align*}
Therefore, for every $a \in \mathbb{R}^p$,
\begin{align*}
|a^\top(\mu-\bar X(\omega))|
\le
\sqrt{\frac{c_\alpha}{n}\,a^\top S(\omega)a}.
\end{align*}
[guided]
Fix one outcome $\omega \in E$. From this point to the end of the step, the argument is deterministic linear algebra. Define the displacement vector $h \in \mathbb{R}^p$ by
\begin{align*}
h := \mu-\bar X(\omega).
\end{align*}
Since $\omega \in E$, the matrix $S(\omega)$ is positive definite, so $S(\omega)^{-1}$ exists and is positive definite. Define
\begin{align*}
(\cdot,\cdot)_{S(\omega)} : \mathbb{R}^p \times \mathbb{R}^p &\to \mathbb{R} \\
(u,v) &\mapsto u^\top S(\omega)^{-1}v .
\end{align*}
This is an inner product because $S(\omega)^{-1}$ is symmetric positive definite.
We need to bound $a^\top h$, where $a \in \mathbb{R}^p$ is arbitrary. The right inner product is the one induced by $S(\omega)^{-1}$ because the Hotelling ellipsoid controls $h^\top S(\omega)^{-1}h$. For completeness, we prove the Cauchy-Schwarz estimate in this inner product. For $u,v \in \mathbb{R}^p$ and $t \in \mathbb{R}$,
\begin{align*}
0
&\le (u-tv,u-tv)_{S(\omega)} \\
&= (u,u)_{S(\omega)} -2t(u,v)_{S(\omega)} + t^2(v,v)_{S(\omega)}.
\end{align*}
If $v \ne 0$, this quadratic polynomial in $t$ is non-negative for every real $t$, so its discriminant is non-positive. Hence
\begin{align*}
4(u,v)_{S(\omega)}^2
-
4(u,u)_{S(\omega)}(v,v)_{S(\omega)}
\le 0,
\end{align*}
which gives
\begin{align*}
|(u,v)_{S(\omega)}|^2
\le
(u,u)_{S(\omega)}(v,v)_{S(\omega)}.
\end{align*}
If $v=0$, the same inequality holds because both sides are zero.
Now take an arbitrary $a \in \mathbb{R}^p$. The expression $a^\top h$ can be rewritten as an inner product with first entry $h$:
\begin{align*}
a^\top h
&= h^\top a \\
&= h^\top S(\omega)^{-1}S(\omega)a \\
&= (h,S(\omega)a)_{S(\omega)}.
\end{align*}
Applying the just-proved Cauchy-Schwarz estimate with $u=h$ and $v=S(\omega)a$ gives
\begin{align*}
|a^\top h|^2
&\le
(h,h)_{S(\omega)}(S(\omega)a,S(\omega)a)_{S(\omega)} \\
&=
h^\top S(\omega)^{-1}h \, a^\top S(\omega)S(\omega)^{-1}S(\omega)a \\
&=
h^\top S(\omega)^{-1}h \, a^\top S(\omega)a.
\end{align*}
Because $\omega \in E$, the Hotelling inequality gives
\begin{align*}
n h^\top S(\omega)^{-1}h \le c_\alpha,
\end{align*}
and hence
\begin{align*}
h^\top S(\omega)^{-1}h \le \frac{c_\alpha}{n}.
\end{align*}
Substituting this into the previous estimate yields
\begin{align*}
|a^\top(\mu-\bar X(\omega))|
\le
\sqrt{\frac{c_\alpha}{n}\,a^\top S(\omega)a}.
\end{align*}
The direction $a$ was arbitrary, and the outcome $\omega$ was fixed only under the condition $\omega \in E$, so this bound holds for every $a \in \mathbb{R}^p$ on the same event $E$.
[/guided]
[/step]
[step:Translate the contrast bound into the simultaneous interval]
For every $\omega \in E$ and every $a \in \mathbb{R}^p$, the previous step gives
\begin{align*}
-\sqrt{\frac{c_\alpha}{n}\,a^\top S(\omega)a}
\le
a^\top\mu-a^\top\bar X(\omega)
\le
\sqrt{\frac{c_\alpha}{n}\,a^\top S(\omega)a}.
\end{align*}
Adding $a^\top\bar X(\omega)$ throughout gives
\begin{align*}
a^\top\bar X(\omega)
-
\sqrt{\frac{c_\alpha}{n}\,a^\top S(\omega)a}
\le
a^\top\mu
\le
a^\top\bar X(\omega)
+
\sqrt{\frac{c_\alpha}{n}\,a^\top S(\omega)a}.
\end{align*}
Since
\begin{align*}
\frac{c_\alpha}{n}
=
\frac{p(n-1)}{n(n-p)}F_{p,n-p;1-\alpha},
\end{align*}
this is exactly the stated interval. Because $\mathbb{P}(E)=1-\alpha$ and the assertion holds for every $a \in \mathbb{R}^p$ on $E$, the intervals have simultaneous coverage probability $1-\alpha$.
[/step]