[proofplan] We first whiten and recenter the observations so that the problem becomes the corresponding statement for independent Gaussian vectors with identity covariance and common mean $\delta=\Sigma^{-1/2}(\mu-\mu_0)$. An orthogonal change of coordinates in the sample index separates the transformed sample mean from the transformed residuals, giving an independent noncentral Gaussian vector and a central Wishart matrix. The remaining distributional reduction is a block decomposition of a Wishart matrix: for $Z\sim \mathcal N_p(\theta,I_p)$ independent of $W\sim W_p(I_p,n-1)$, the quadratic form $Z^\top W^{-1}Z$ has the same law as $\chi^2_p(|\theta|^2)/\chi^2_{n-p}$ with independent numerator and denominator. This is exactly the defining representation of the noncentral $F$ distribution after multiplying by $(n-p)/p$. [/proofplan]

[step:Whiten the observations and rewrite Hotelling's statistic]Let $\Sigma^{1/2}$ denote the unique symmetric positive definite square root of $\Sigma$, and let $\Sigma^{-1/2}$ denote its inverse. Define \begin{align*} \delta &= \Sigma^{-1/2}(\mu-\mu_0) \in \mathbb{R}^p, \end{align*} and define transformed random vectors \begin{align*} Y_i : \Omega &\to \mathbb{R}^p \\ \omega &\mapsto \Sigma^{-1/2}(X_i(\omega)-\mu_0) \end{align*} for $i=1,\dots,n$. Since affine transformations of multivariate Gaussian random vectors are Gaussian, the random vectors $Y_1,\dots,Y_n$ are independent with common distribution $\mathcal N_p(\delta,I_p)$. Define \begin{align*} \bar Y &= \frac{1}{n}\sum_{i=1}^n Y_i, \\ S_Y &= \frac{1}{n-1}\sum_{i=1}^n (Y_i-\bar Y)(Y_i-\bar Y)^\top. \end{align*} Then \begin{align*} \bar Y = \Sigma^{-1/2}(\bar X-\mu_0), \end{align*} and \begin{align*} S_Y = \Sigma^{-1/2}S\Sigma^{-1/2}. \end{align*} Since $S_Y^{-1}=\Sigma^{1/2}S^{-1}\Sigma^{1/2}$ whenever $S$ is invertible, Hotelling's statistic may be written as \begin{align*} T^2 &= n(\bar X-\mu_0)^\top S^{-1}(\bar X-\mu_0) \\ &= n\bar Y^\top S_Y^{-1}\bar Y. \end{align*} Thus it suffices to determine the distribution of $n\bar Y^\top S_Y^{-1}\bar Y$ in the standardized model.[/step]

[guided]The first move is to remove the covariance matrix from the problem. Because $\Sigma$ is symmetric positive definite, it has a unique symmetric positive definite square root $\Sigma^{1/2}$, and $\Sigma^{-1/2}$ is well-defined. We introduce the vector \begin{align*} \delta = \Sigma^{-1/2}(\mu-\mu_0), \end{align*} which is the mean shift measured in covariance-standardized coordinates. For each observation, define the whitened random vector \begin{align*} Y_i : \Omega &\to \mathbb{R}^p \\ \omega &\mapsto \Sigma^{-1/2}(X_i(\omega)-\mu_0). \end{align*} The transformation $x\mapsto \Sigma^{-1/2}(x-\mu_0)$ is affine. Therefore each $Y_i$ is Gaussian with mean \begin{align*} \mathbb{E}[Y_i] = \Sigma^{-1/2}(\mu-\mu_0)=\delta \end{align*} and covariance \begin{align*} \Sigma^{-1/2}\Sigma(\Sigma^{-1/2})^\top = I_p, \end{align*} because $\Sigma^{-1/2}$ is symmetric. Independence is preserved because each $Y_i$ is a measurable function of $X_i$ alone. Hence $Y_1,\dots,Y_n$ are independent with common distribution $\mathcal N_p(\delta,I_p)$. Now define \begin{align*} \bar Y &= \frac{1}{n}\sum_{i=1}^n Y_i, \\ S_Y &= \frac{1}{n-1}\sum_{i=1}^n (Y_i-\bar Y)(Y_i-\bar Y)^\top. \end{align*} From the definition of $Y_i$, \begin{align*} \bar Y=\Sigma^{-1/2}(\bar X-\mu_0). \end{align*} Also, \begin{align*} Y_i-\bar Y=\Sigma^{-1/2}(X_i-\bar X), \end{align*} so \begin{align*} S_Y &= \frac{1}{n-1}\sum_{i=1}^n \Sigma^{-1/2}(X_i-\bar X)(X_i-\bar X)^\top \Sigma^{-1/2} \\ &= \Sigma^{-1/2}S\Sigma^{-1/2}. \end{align*} Therefore, on the almost sure event where $S$ is invertible, \begin{align*} S_Y^{-1}=\Sigma^{1/2}S^{-1}\Sigma^{1/2}. \end{align*} Substituting these identities into Hotelling's statistic gives \begin{align*} T^2 &= n(\bar X-\mu_0)^\top S^{-1}(\bar X-\mu_0) \\ &= n\bar Y^\top S_Y^{-1}\bar Y. \end{align*} Thus the original statistic is unchanged by whitening, provided we express it using the whitened sample covariance matrix. The noncentrality parameter also becomes \begin{align*} \lambda &= n(\mu-\mu_0)^\top \Sigma^{-1}(\mu-\mu_0) \\ &= n|\delta|^2. \end{align*}[/guided]

[step:Separate the sample mean from the residual directions]Choose an orthogonal matrix $Q\in \mathbb{R}^{n\times n}$ whose first row is \begin{align*} \left(\frac{1}{\sqrt n},\dots,\frac{1}{\sqrt n}\right). \end{align*} Define random vectors $Z_1,\dots,Z_n : \Omega \to \mathbb{R}^p$ by \begin{align*} Z_a = \sum_{i=1}^n Q_{ai}Y_i \end{align*} for $a=1,\dots,n$. Since $(Y_1,\dots,Y_n)$ is jointly Gaussian with independent components across the sample index and common covariance $I_p$, orthogonality of $Q$ implies that $Z_1,\dots,Z_n$ are independent Gaussian random vectors. Their means are \begin{align*} \mathbb{E}[Z_1] &= \sqrt n\,\delta, \\ \mathbb{E}[Z_a] &= 0 \qquad \text{for } a=2,\dots,n, \end{align*} and each has covariance matrix $I_p$. Hence \begin{align*} Z_1 \sim \mathcal N_p(\sqrt n\,\delta,I_p), \end{align*} while $Z_2,\dots,Z_n$ are independent $\mathcal N_p(0,I_p)$ random vectors, independent of $Z_1$. Moreover, \begin{align*} Z_1=\sqrt n\,\bar Y, \end{align*} and orthogonality of $Q$ gives the residual sum-of-squares identity \begin{align*} \sum_{i=1}^n (Y_i-\bar Y)(Y_i-\bar Y)^\top = \sum_{a=2}^n Z_aZ_a^\top. \end{align*} Define \begin{align*} W = \sum_{a=2}^n Z_aZ_a^\top. \end{align*} Then $W$ is independent of $Z_1$, has the central Wishart distribution $W_p(I_p,n-1)$, and \begin{align*} S_Y = \frac{1}{n-1}W. \end{align*} Consequently, \begin{align*} T^2 = n\bar Y^\top S_Y^{-1}\bar Y = (n-1)Z_1^\top W^{-1}Z_1. \end{align*}[/step]

[guided]The purpose of this step is to isolate the part of the data containing the sample mean from the part containing the sample covariance. We do this by rotating the $n$ observations in the sample-index direction. Choose an orthogonal matrix $Q\in \mathbb{R}^{n\times n}$ whose first row is \begin{align*} \left(\frac{1}{\sqrt n},\dots,\frac{1}{\sqrt n}\right). \end{align*} For each $a=1,\dots,n$, define \begin{align*} Z_a=\sum_{i=1}^n Q_{ai}Y_i. \end{align*} This is a linear transformation of the jointly Gaussian vector $(Y_1,\dots,Y_n)$, so $(Z_1,\dots,Z_n)$ is jointly Gaussian. Orthogonality of $Q$ preserves the covariance in the sample-index direction. More explicitly, for $a,b\in\{1,\dots,n\}$, \begin{align*} \operatorname{Cov}(Z_a,Z_b) &= \sum_{i=1}^n\sum_{j=1}^n Q_{ai}Q_{bj}\operatorname{Cov}(Y_i,Y_j) \\ &= \sum_{i=1}^n Q_{ai}Q_{bi} I_p \\ &= \delta_{ab}I_p, \end{align*} where $\delta_{ab}$ denotes the Kronecker delta. Since the vector is jointly Gaussian, zero covariance between distinct $Z_a$ and $Z_b$ implies independence. The means are computed from $\mathbb{E}[Y_i]=\delta$: \begin{align*} \mathbb{E}[Z_a] &= \sum_{i=1}^n Q_{ai}\delta \\ &= \left(\sum_{i=1}^n Q_{ai}\right)\delta. \end{align*} For $a=1$, the row sum is $\sqrt n$, so $\mathbb{E}[Z_1]=\sqrt n\,\delta$. For $a\ge 2$, the row of $Q$ is orthogonal to the first row, hence its entries sum to $0$, so $\mathbb{E}[Z_a]=0$. Therefore \begin{align*} Z_1\sim \mathcal N_p(\sqrt n\,\delta,I_p), \end{align*} and $Z_2,\dots,Z_n$ are independent central $\mathcal N_p(0,I_p)$ random vectors, independent of $Z_1$. The first transformed vector is the scaled sample mean: \begin{align*} Z_1=\sum_{i=1}^n \frac{1}{\sqrt n}Y_i=\sqrt n\,\bar Y. \end{align*} The remaining transformed vectors span the orthogonal complement of the constant sample direction, so they contain exactly the centered residual information. Orthogonality of $Q$ gives \begin{align*} \sum_{i=1}^n Y_iY_i^\top = \sum_{a=1}^n Z_aZ_a^\top. \end{align*} Since $Z_1Z_1^\top=n\bar Y\bar Y^\top$, subtracting this rank-one mean part gives \begin{align*} \sum_{i=1}^n (Y_i-\bar Y)(Y_i-\bar Y)^\top = \sum_{a=2}^n Z_aZ_a^\top. \end{align*} Define \begin{align*} W=\sum_{a=2}^n Z_aZ_a^\top. \end{align*} Because $Z_2,\dots,Z_n$ are independent $\mathcal N_p(0,I_p)$ random vectors, $W$ has the central Wishart distribution $W_p(I_p,n-1)$. It is independent of $Z_1$ because it is a measurable function of $Z_2,\dots,Z_n$, which are independent of $Z_1$. Finally, \begin{align*} S_Y=\frac{1}{n-1}W, \end{align*} and therefore \begin{align*} T^2 &= n\bar Y^\top S_Y^{-1}\bar Y \\ &= Z_1^\top \left(\frac{1}{n-1}W\right)^{-1}Z_1 \\ &= (n-1)Z_1^\top W^{-1}Z_1. \end{align*} This reduces the theorem to a distributional statement about one noncentral Gaussian vector independent of one central Wishart matrix.[/guided]

[step:Reduce the inverse Wishart quadratic form to independent chi-square variables]Let $m=n-1$ and let $\nu=m-p+1=n-p$. We prove the following distributional fact in the present setting: if $Z:\Omega\to\mathbb{R}^p$ satisfies $Z\sim\mathcal N_p(\theta,I_p)$, if $W\sim W_p(I_p,m)$ is independent of $Z$, and if $m\ge p$, then \begin{align*} Z^\top W^{-1}Z \sim \frac{U}{V}, \end{align*} where $U\sim\chi^2_p(|\theta|^2)$, $V\sim\chi^2_{\nu}$, and $U,V$ are independent. Condition on the value $Z=z\in\mathbb{R}^p$. If $z=0$, then $z^\top W^{-1}z=0$. If $z\neq 0$, choose an orthogonal matrix $R_z\in\mathbb{R}^{p\times p}$ such that $R_zz=|z|e_1$, where $e_1=(1,0,\dots,0)^\top$. Since $W\sim W_p(I_p,m)$ is orthogonally invariant, $R_zWR_z^\top$ has the same distribution as $W$. Therefore the conditional distribution of $z^\top W^{-1}z$ equals the distribution of \begin{align*} |z|^2 e_1^\top W^{-1}e_1. \end{align*} Represent $W$ as \begin{align*} W=\sum_{r=1}^m G_rG_r^\top, \end{align*} where $G_1,\dots,G_m:\Omega\to\mathbb{R}^p$ are independent $\mathcal N_p(0,I_p)$ random vectors. Let $A\in\mathbb{R}^{m\times p}$ be the random matrix whose $r$-th row is $G_r^\top$, so $W=A^\top A$. Write $A=(a,B)$, where $a:\Omega\to\mathbb{R}^m$ is the first column and $B:\Omega\to\mathbb{R}^{m\times(p-1)}$ is the matrix of the remaining columns. Then $a$ is a standard Gaussian vector in $\mathbb{R}^m$, $B$ is independent of $a$, and almost surely $B$ has rank $p-1$. The Schur complement formula for the inverse of $W=A^\top A$ gives \begin{align*} e_1^\top W^{-1}e_1 = \frac{1}{a^\top (I_m-P_B)a}, \end{align*} where $I_m$ is the identity matrix on $\mathbb{R}^m$ and $P_B:\mathbb{R}^m\to\mathbb{R}^m$ is the [orthogonal projection](/theorems/437) onto the column space of $B$. Conditional on $B$, the matrix $I_m-P_B$ is the orthogonal projection onto a subspace of dimension \begin{align*} m-(p-1)=m-p+1=\nu. \end{align*} Since $a\sim\mathcal N_m(0,I_m)$ and $a$ is independent of $B$, the quadratic form \begin{align*} a^\top(I_m-P_B)a \end{align*} has the $\chi^2_{\nu}$ distribution, conditionally on $B$ and hence unconditionally. Denote this random variable by $V$. Its conditional distribution does not depend on $z$, and $W$ is independent of $Z$, so $V$ may be chosen independent of $Z$. Thus, conditionally on $Z=z$, \begin{align*} Z^\top W^{-1}Z \sim \frac{|z|^2}{V}, \end{align*} with $V\sim\chi^2_{\nu}$ independent of $Z$. Since $Z\sim\mathcal N_p(\theta,I_p)$, the squared norm \begin{align*} U=|Z|^2 \end{align*} has the noncentral chi-square distribution $\chi^2_p(|\theta|^2)$. Hence \begin{align*} Z^\top W^{-1}Z \sim \frac{U}{V}, \end{align*} with $U$ and $V$ independent.[/step]

[guided]We now need the exact distribution of the quadratic form $Z^\top W^{-1}Z$, where $Z$ is Gaussian and $W$ is an independent central Wishart matrix. The key point is that, after conditioning on $Z=z$, the Wishart matrix is rotationally symmetric. This lets us rotate $z$ onto the first coordinate axis and compute only the upper-left entry of $W^{-1}$. Let $m=n-1$ and $\nu=m-p+1=n-p$. Consider the general situation where \begin{align*} Z\sim\mathcal N_p(\theta,I_p), \end{align*} where $\theta\in\mathbb{R}^p$, and where $W\sim W_p(I_p,m)$ is independent of $Z$. Since $m=n-1\ge p$, the Wishart matrix $W$ is invertible almost surely. Fix $z\in\mathbb{R}^p$ and condition on the event $Z=z$ in the regular conditional distribution sense. If $z=0$, then $z^\top W^{-1}z=0$, so there is nothing to rotate. If $z\neq 0$, choose an orthogonal matrix $R_z\in\mathbb{R}^{p\times p}$ satisfying \begin{align*} R_zz=|z|e_1, \end{align*} where $e_1=(1,0,\dots,0)^\top$. The central Wishart distribution with scale matrix $I_p$ is invariant under orthogonal conjugation: $R_zWR_z^\top$ has the same distribution as $W$. Therefore \begin{align*} z^\top W^{-1}z &= (R_zz)^\top (R_zWR_z^\top)^{-1}(R_zz) \\ &\sim |z|^2 e_1^\top W^{-1}e_1 \end{align*} conditionally on $Z=z$. It remains to compute $e_1^\top W^{-1}e_1$. Represent the Wishart matrix as \begin{align*} W=\sum_{r=1}^m G_rG_r^\top, \end{align*} where $G_1,\dots,G_m$ are independent $\mathcal N_p(0,I_p)$ random vectors. Equivalently, if $A\in\mathbb{R}^{m\times p}$ is the matrix whose $r$-th row is $G_r^\top$, then \begin{align*} W=A^\top A. \end{align*} Write the columns of $A$ as \begin{align*} A=(a,B), \end{align*} where $a:\Omega\to\mathbb{R}^m$ is the first column and $B:\Omega\to\mathbb{R}^{m\times(p-1)}$ contains the remaining columns. The entries of $A$ are independent standard normal random variables, so $a$ is a standard Gaussian vector in $\mathbb{R}^m$, $B$ is independent of $a$, and $B$ has rank $p-1$ almost surely because $m\ge p$. With respect to this block decomposition, \begin{align*} W=A^\top A = \begin{pmatrix} a^\top a & a^\top B \\ B^\top a & B^\top B \end{pmatrix}. \end{align*} The Schur complement formula says that the upper-left entry of $W^{-1}$ is the reciprocal of the Schur complement of $B^\top B$: \begin{align*} e_1^\top W^{-1}e_1 = \frac{1}{a^\top a-a^\top B(B^\top B)^{-1}B^\top a}. \end{align*} Let $P_B:\mathbb{R}^m\to\mathbb{R}^m$ be the orthogonal projection onto the column space of $B$. Since \begin{align*} P_B=B(B^\top B)^{-1}B^\top \end{align*} almost surely, the denominator becomes \begin{align*} a^\top(I_m-P_B)a. \end{align*} Conditional on $B$, the operator $I_m-P_B$ is the orthogonal projection onto the orthogonal complement of the column space of $B$. Since $B$ has rank $p-1$ almost surely, this orthogonal complement has dimension \begin{align*} m-(p-1)=m-p+1=\nu. \end{align*} Because $a\sim\mathcal N_m(0,I_m)$ and is independent of $B$, the squared Euclidean length of the projection of $a$ onto this $\nu$-dimensional subspace has the $\chi^2_\nu$ distribution. Thus \begin{align*} a^\top(I_m-P_B)a\sim\chi^2_\nu. \end{align*} Call this denominator $V$. The conditional distribution of $V$ does not depend on $z$, and the entire Wishart construction is independent of $Z$. Therefore $V$ is independent of $Z$. We have shown that, conditionally on $Z=z$, \begin{align*} Z^\top W^{-1}Z \sim \frac{|z|^2}{V}, \end{align*} where $V\sim\chi^2_\nu$ is independent of $Z$. Finally, since $Z\sim\mathcal N_p(\theta,I_p)$, its squared norm \begin{align*} U=|Z|^2 \end{align*} has the noncentral chi-square distribution $\chi^2_p(|\theta|^2)$. Hence \begin{align*} Z^\top W^{-1}Z \sim \frac{U}{V}, \end{align*} where $U\sim\chi^2_p(|\theta|^2)$, $V\sim\chi^2_\nu$, and $U,V$ are independent.[/guided]

[step:Identify the noncentral F distribution]Apply the previous step with \begin{align*} Z=Z_1,\qquad \theta=\sqrt n\,\delta,\qquad m=n-1,\qquad \nu=n-p. \end{align*} Then \begin{align*} |\theta|^2 = n|\delta|^2 = n(\mu-\mu_0)^\top\Sigma^{-1}(\mu-\mu_0) = \lambda. \end{align*} Therefore there exist independent random variables \begin{align*} U &\sim \chi^2_p(\lambda), \\ V &\sim \chi^2_{n-p} \end{align*} such that \begin{align*} Z_1^\top W^{-1}Z_1 \sim \frac{U}{V}. \end{align*} Using $T^2=(n-1)Z_1^\top W^{-1}Z_1$, we obtain \begin{align*} \frac{n-p}{p(n-1)}T^2 &= \frac{n-p}{p}Z_1^\top W^{-1}Z_1 \\ &\sim \frac{U/p}{V/(n-p)}. \end{align*} By definition, the distribution of $(U/p)/(V/(n-p))$, with $U\sim\chi^2_p(\lambda)$ independent of $V\sim\chi^2_{n-p}$, is the noncentral $F$ distribution $F_{p,n-p}(\lambda)$. Hence \begin{align*} \frac{n-p}{p(n-1)}T^2 \sim F_{p,n-p}(\lambda). \end{align*} This proves the theorem.[/step]

[guided]We now substitute the objects obtained from the data into the distributional reduction. In our case, \begin{align*} Z=Z_1,\qquad \theta=\sqrt n\,\delta,\qquad m=n-1. \end{align*} The denominator degrees of freedom from the previous step are \begin{align*} \nu=m-p+1=(n-1)-p+1=n-p. \end{align*} The noncentrality parameter of $|Z_1|^2$ is the squared length of its mean: \begin{align*} |\theta|^2 &= |\sqrt n\,\delta|^2 \\ &= n|\delta|^2 \\ &= n(\mu-\mu_0)^\top\Sigma^{-1}(\mu-\mu_0) \\ &= \lambda. \end{align*} Therefore the previous step gives independent random variables \begin{align*} U &\sim \chi^2_p(\lambda), \\ V &\sim \chi^2_{n-p} \end{align*} such that \begin{align*} Z_1^\top W^{-1}Z_1 \sim \frac{U}{V}. \end{align*} From the sample-mean and residual decomposition, we already proved \begin{align*} T^2=(n-1)Z_1^\top W^{-1}Z_1. \end{align*} Multiplying by the normalizing factor in the theorem gives \begin{align*} \frac{n-p}{p(n-1)}T^2 &= \frac{n-p}{p}Z_1^\top W^{-1}Z_1 \\ &\sim \frac{n-p}{p}\frac{U}{V} \\ &= \frac{U/p}{V/(n-p)}. \end{align*} The defining representation of the noncentral $F$ distribution with numerator degrees of freedom $p$, denominator degrees of freedom $n-p$, and noncentrality parameter $\lambda$ is precisely \begin{align*} \frac{U/p}{V/(n-p)}, \end{align*} where $U\sim\chi^2_p(\lambda)$ and $V\sim\chi^2_{n-p}$ are independent. Hence \begin{align*} \frac{n-p}{p(n-1)}T^2\sim F_{p,n-p}(\lambda). \end{align*} This is the desired distributional identity.[/guided]