[proofplan] We first compute the distribution of the sample mean under $H_0$ using characteristic functions and independence. Then we standardize the sample mean by the symmetric inverse square root of $\Sigma$, obtaining a standard Gaussian vector in $\mathbb{R}^p$. Finally, we identify the statistic $Q$ as the squared Euclidean norm of this standard Gaussian vector, hence as the sum of squares of $p$ independent standard normal random variables. [/proofplan] [step:Compute the Gaussian distribution of the sample mean under $H_0$] For a random vector $Y: \Omega \to \mathbb{R}^p$, define its characteristic function $\phi_Y:\mathbb{R}^p \to \mathbb{C}$ by \begin{align*} \phi_Y(t) := \mathbb{E}\left[e^{i t^\top Y}\right]. \end{align*} Under $H_0$, each $X_k$ has distribution $\mathcal{N}_p(\mu_0,\Sigma)$, so \begin{align*} \phi_{X_k}(t) = \exp\left(i t^\top \mu_0-\frac{1}{2}t^\top \Sigma t\right) \end{align*} for every $t \in \mathbb{R}^p$. Since $X_1,\dots,X_n$ are independent, \begin{align*} \phi_{\bar{X}}(t) &= \mathbb{E}\left[ \exp\left(i t^\top \frac{1}{n}\sum_{k=1}^{n}X_k\right) \right] \\ &= \prod_{k=1}^{n} \mathbb{E}\left[ \exp\left(i \left(\frac{t}{n}\right)^\top X_k\right) \right] \\ &= \prod_{k=1}^{n} \exp\left( i\left(\frac{t}{n}\right)^\top \mu_0 -\frac{1}{2}\left(\frac{t}{n}\right)^\top \Sigma \left(\frac{t}{n}\right) \right) \\ &= \exp\left( i t^\top \mu_0 -\frac{1}{2}t^\top \left(\frac{\Sigma}{n}\right)t \right). \end{align*} This is the characteristic function of $\mathcal{N}_p(\mu_0,\Sigma/n)$, so \begin{align*} \bar{X} \sim \mathcal{N}_p\left(\mu_0,\frac{\Sigma}{n}\right). \end{align*} [guided] The point of this step is to justify the distribution of $\bar{X}$ rather than simply quote it. For a random vector $Y: \Omega \to \mathbb{R}^p$, define its characteristic function $\phi_Y:\mathbb{R}^p \to \mathbb{C}$ by \begin{align*} \phi_Y(t) := \mathbb{E}\left[e^{i t^\top Y}\right]. \end{align*} Under the null hypothesis $H_0:\mu=\mu_0$, each observation satisfies $X_k \sim \mathcal{N}_p(\mu_0,\Sigma)$. Therefore, for every $t \in \mathbb{R}^p$, \begin{align*} \phi_{X_k}(t) = \exp\left(i t^\top \mu_0-\frac{1}{2}t^\top \Sigma t\right). \end{align*} Now compute the characteristic function of the sample mean. Since \begin{align*} \bar{X}=\frac{1}{n}\sum_{k=1}^{n}X_k, \end{align*} we have \begin{align*} \phi_{\bar{X}}(t) &= \mathbb{E}\left[ \exp\left(i t^\top \frac{1}{n}\sum_{k=1}^{n}X_k\right) \right]. \end{align*} The exponent separates into a sum: \begin{align*} i t^\top \frac{1}{n}\sum_{k=1}^{n}X_k = \sum_{k=1}^{n} i\left(\frac{t}{n}\right)^\top X_k. \end{align*} Because the random vectors $X_1,\dots,X_n$ are independent, the expectation of the product of the corresponding exponential functions factors: \begin{align*} \phi_{\bar{X}}(t) &= \prod_{k=1}^{n} \mathbb{E}\left[ \exp\left(i \left(\frac{t}{n}\right)^\top X_k\right) \right]. \end{align*} Using the characteristic function of each $X_k$ at $t/n$ gives \begin{align*} \phi_{\bar{X}}(t) &= \prod_{k=1}^{n} \exp\left( i\left(\frac{t}{n}\right)^\top \mu_0 -\frac{1}{2}\left(\frac{t}{n}\right)^\top \Sigma \left(\frac{t}{n}\right) \right) \\ &= \exp\left( i t^\top \mu_0 -\frac{1}{2}t^\top \left(\frac{\Sigma}{n}\right)t \right). \end{align*} This is exactly the characteristic function of a multivariate Gaussian random vector with mean $\mu_0$ and covariance matrix $\Sigma/n$. Hence \begin{align*} \bar{X} \sim \mathcal{N}_p\left(\mu_0,\frac{\Sigma}{n}\right). \end{align*} [/guided] [/step] [step:Standardize the sample mean using the inverse square root of $\Sigma$] Let $I_p \in \mathbb{R}^{p \times p}$ denote the identity matrix. Since $\Sigma$ is symmetric positive definite, the spectral theorem gives an orthonormal eigenbasis with positive eigenvalues; applying the positive function $\lambda \mapsto \lambda^{-1/2}$ to these eigenvalues produces a symmetric positive definite matrix $\Sigma^{-1/2} \in \mathbb{R}^{p \times p}$ satisfying \begin{align*} \Sigma^{-1/2}\Sigma\Sigma^{-1/2}=I_p. \end{align*} Define the deterministic matrix $A\in\mathbb{R}^{p\times p}$ and deterministic vector $b\in\mathbb{R}^p$ by \begin{align*} A := \sqrt{n}\,\Sigma^{-1/2}, \qquad b := -\sqrt{n}\,\Sigma^{-1/2}\mu_0. \end{align*} The dimensions are compatible because $\bar{X}:\Omega\to\mathbb{R}^p$, $A$ maps $\mathbb{R}^p$ to $\mathbb{R}^p$, and $b\in\mathbb{R}^p$. Define the standardized random vector $Z:\Omega \to \mathbb{R}^p$ by \begin{align*} Z := A\bar{X}+b = \sqrt{n}\,\Sigma^{-1/2}(\bar{X}-\mu_0). \end{align*} Since $\bar{X}\sim \mathcal{N}_p(\mu_0,\Sigma/n)$ and $Z$ is the affine image of $\bar{X}$ under the deterministic compatible-dimensional map $x\mapsto Ax+b$, $Z$ is Gaussian with mean \begin{align*} \sqrt{n}\,\Sigma^{-1/2}(\mu_0-\mu_0)=0 \end{align*} and covariance \begin{align*} n\,\Sigma^{-1/2}\left(\frac{\Sigma}{n}\right)\Sigma^{-1/2} = I_p. \end{align*} Therefore \begin{align*} Z \sim \mathcal{N}_p(0,I_p). \end{align*} [guided] We now remove the mean and covariance from $\bar{X}$. Let $I_p \in \mathbb{R}^{p \times p}$ denote the identity matrix. Since $\Sigma$ is symmetric positive definite, the spectral theorem applies: $\Sigma$ has an orthonormal eigenbasis and all eigenvalues are positive. Applying the positive function $\lambda \mapsto \lambda^{-1/2}$ to these eigenvalues defines a symmetric positive definite matrix $\Sigma^{-1/2}\in\mathbb{R}^{p\times p}$ such that \begin{align*} \Sigma^{-1/2}\Sigma\Sigma^{-1/2}=I_p. \end{align*} Define the deterministic matrix $A\in\mathbb{R}^{p\times p}$ and deterministic vector $b\in\mathbb{R}^p$ by \begin{align*} A := \sqrt{n}\,\Sigma^{-1/2}, \qquad b := -\sqrt{n}\,\Sigma^{-1/2}\mu_0. \end{align*} The matrix-vector products are well-defined: $\bar{X}:\Omega\to\mathbb{R}^p$, $A$ is a $p\times p$ deterministic matrix, and $b$ is a deterministic vector in $\mathbb{R}^p$. Define the random vector $Z:\Omega \to \mathbb{R}^p$ by \begin{align*} Z := A\bar{X}+b = \sqrt{n}\,\Sigma^{-1/2}(\bar{X}-\mu_0). \end{align*} This definition is chosen so that the covariance normalizes to the identity. From the previous step, \begin{align*} \bar{X}\sim \mathcal{N}_p\left(\mu_0,\frac{\Sigma}{n}\right). \end{align*} The affine stability property for Gaussian random vectors applies because the transformation $x\mapsto Ax+b$ is deterministic and maps $\mathbb{R}^p$ to $\mathbb{R}^p$. Hence $Z$ is Gaussian. The mean of $Z$ is \begin{align*} \mathbb{E}[Z] = \sqrt{n}\,\Sigma^{-1/2}(\mu_0-\mu_0) = 0, \end{align*} and its covariance matrix is \begin{align*} \operatorname{Cov}(Z) &= n\,\Sigma^{-1/2}\left(\frac{\Sigma}{n}\right)\Sigma^{-1/2} \\ &= \Sigma^{-1/2}\Sigma\Sigma^{-1/2} \\ &= I_p. \end{align*} Thus $Z$ is a centered Gaussian random vector with identity covariance: \begin{align*} Z\sim \mathcal{N}_p(0,I_p). \end{align*} [/guided] [/step] [step:Identify the quadratic statistic as a chi-squared random variable] Write $Z=(Z_1,\dots,Z_p)$, where each $Z_j:\Omega\to\mathbb{R}$ is the $j$-th coordinate random variable of $Z$. Since $Z\sim \mathcal{N}_p(0,I_p)$, the coordinates $Z_1,\dots,Z_p$ are independent standard normal random variables. Using the definition of $Z$ and the symmetry of $\Sigma^{-1/2}$, \begin{align*} Z^\top Z &= \left(\sqrt{n}\,\Sigma^{-1/2}(\bar{X}-\mu_0)\right)^\top \left(\sqrt{n}\,\Sigma^{-1/2}(\bar{X}-\mu_0)\right) \\ &= n(\bar{X}-\mu_0)^\top \Sigma^{-1/2}\Sigma^{-1/2} (\bar{X}-\mu_0) \\ &= n(\bar{X}-\mu_0)^\top \Sigma^{-1} (\bar{X}-\mu_0) \\ &= Q. \end{align*} Therefore \begin{align*} Q = Z^\top Z = \sum_{j=1}^{p} Z_j^2. \end{align*} By the definition of the chi-squared distribution with $p$ degrees of freedom as the sum of squares of $p$ independent standard normal random variables, \begin{align*} Q \sim \chi^2_p. \end{align*} This proves the claim. [guided] We now translate the standardized Gaussian vector into the stated quadratic statistic. Write $Z=(Z_1,\dots,Z_p)$, where each coordinate map $Z_j:\Omega\to\mathbb{R}$ is defined by projecting $Z$ onto its $j$-th coordinate. Since $Z\sim \mathcal{N}_p(0,I_p)$ and $I_p$ is diagonal with diagonal entries equal to $1$, the defining coordinate description of the standard multivariate normal law gives that $Z_1,\dots,Z_p$ are independent real-valued random variables with $Z_j\sim\mathcal{N}(0,1)$ for every $j\in\{1,\dots,p\}$. It remains to check that the statistic in the theorem is exactly the squared Euclidean norm of $Z$. Using the definition of $Z$ and the symmetry of $\Sigma^{-1/2}$, we compute \begin{align*} Z^\top Z &= \left(\sqrt{n}\,\Sigma^{-1/2}(\bar{X}-\mu_0)\right)^\top \left(\sqrt{n}\,\Sigma^{-1/2}(\bar{X}-\mu_0)\right) \\ &= n(\bar{X}-\mu_0)^\top \Sigma^{-1/2}\Sigma^{-1/2} (\bar{X}-\mu_0) \\ &= n(\bar{X}-\mu_0)^\top \Sigma^{-1} (\bar{X}-\mu_0) \\ &= Q. \end{align*} Because $Z^\top Z$ is the squared Euclidean norm of $Z$, it is the sum of the coordinate squares: \begin{align*} Q=Z^\top Z=\sum_{j=1}^{p} Z_j^2. \end{align*} By the definition of the chi-squared distribution with $p$ degrees of freedom as the sum of squares of $p$ independent standard normal random variables, this gives \begin{align*} Q\sim \chi^2_p. \end{align*} This proves the claim. [/guided] [/step]