[proofplan] We reduce the quadratic form to the Euclidean norm of a standard Gaussian vector. Since $\Sigma$ is symmetric positive definite, its inverse square root whitens $X$, producing a random vector $Z = \Sigma^{-1/2}X$ with distribution $\mathcal{N}_p(0,I_p)$. The quadratic form then becomes $Z^\top Z$, and the density of $\mathcal{N}_p(0,I_p)$ factors into the product of $p$ one-dimensional standard normal densities. Hence the components of $Z$ are independent standard normal random variables, so their squared sum has the $\chi_p^2$ distribution by definition. [/proofplan] [step:Whiten the Gaussian vector with the inverse square root of $\Sigma$] Since $\Sigma$ is symmetric positive definite, there is a symmetric positive definite matrix $\Sigma^{1/2} \in \mathbb{R}^{p \times p}$ such that \begin{align*} \Sigma^{1/2}\Sigma^{1/2}=\Sigma. \end{align*} Its inverse is denoted by $\Sigma^{-1/2}:=(\Sigma^{1/2})^{-1}$. Define the random vector \begin{align*} Z: \Omega &\to \mathbb{R}^p \\ \omega &\mapsto \Sigma^{-1/2}X(\omega). \end{align*} For each $t \in \mathbb{R}^p$, the moment generating function of $Z$ satisfies \begin{align*} \mathbb{E}\left[e^{t^\top Z}\right] &= \mathbb{E}\left[e^{t^\top \Sigma^{-1/2}X}\right] \\ &= \mathbb{E}\left[e^{((\Sigma^{-1/2})^\top t)^\top X}\right] \\ &= \exp\left(\frac{1}{2} t^\top \Sigma^{-1/2}\Sigma(\Sigma^{-1/2})^\top t\right). \end{align*} Because $\Sigma^{-1/2}$ is symmetric and $\Sigma^{-1/2}\Sigma\Sigma^{-1/2}=I_p$, this becomes \begin{align*} \mathbb{E}\left[e^{t^\top Z}\right] = \exp\left(\frac{1}{2}t^\top t\right). \end{align*} This is the moment generating function of $\mathcal{N}_p(0,I_p)$, so \begin{align*} Z \sim \mathcal{N}_p(0,I_p). \end{align*} [guided] The goal is to remove the covariance matrix from the Gaussian vector. Since $\Sigma$ is symmetric positive definite, it has a symmetric positive definite square root $\Sigma^{1/2}$ with $\Sigma^{1/2}\Sigma^{1/2}=\Sigma$. Its inverse $\Sigma^{-1/2}:=(\Sigma^{1/2})^{-1}$ is also symmetric. Define \begin{align*} Z: \Omega &\to \mathbb{R}^p \\ \omega &\mapsto \Sigma^{-1/2}X(\omega). \end{align*} We compute the distribution of $Z$ through its moment generating function. For $t \in \mathbb{R}^p$, \begin{align*} \mathbb{E}\left[e^{t^\top Z}\right] &= \mathbb{E}\left[e^{t^\top \Sigma^{-1/2}X}\right] \\ &= \mathbb{E}\left[e^{((\Sigma^{-1/2})^\top t)^\top X}\right]. \end{align*} Since $X \sim \mathcal{N}_p(0,\Sigma)$, its moment generating function is \begin{align*} \mathbb{E}\left[e^{s^\top X}\right] = \exp\left(\frac{1}{2}s^\top \Sigma s\right) \end{align*} for every $s \in \mathbb{R}^p$. We apply this with $s=(\Sigma^{-1/2})^\top t$. Then \begin{align*} \mathbb{E}\left[e^{t^\top Z}\right] &= \exp\left(\frac{1}{2}t^\top \Sigma^{-1/2}\Sigma(\Sigma^{-1/2})^\top t\right). \end{align*} Because $\Sigma^{-1/2}$ is symmetric and is the inverse of $\Sigma^{1/2}$, we have \begin{align*} \Sigma^{-1/2}\Sigma(\Sigma^{-1/2})^\top = \Sigma^{-1/2}\Sigma\Sigma^{-1/2} = I_p. \end{align*} Therefore \begin{align*} \mathbb{E}\left[e^{t^\top Z}\right] = \exp\left(\frac{1}{2}t^\top t\right), \end{align*} which is exactly the moment generating function of $\mathcal{N}_p(0,I_p)$. Hence \begin{align*} Z \sim \mathcal{N}_p(0,I_p). \end{align*} [/guided] [/step] [step:Rewrite the Mahalanobis form as the squared Euclidean norm of $Z$] For each $\omega \in \Omega$, the definition $Z(\omega)=\Sigma^{-1/2}X(\omega)$ gives \begin{align*} X(\omega)=\Sigma^{1/2}Z(\omega). \end{align*} Using symmetry of $\Sigma^{1/2}$ and the identity $\Sigma^{-1}=\Sigma^{-1/2}\Sigma^{-1/2}$, we compute \begin{align*} X(\omega)^\top\Sigma^{-1}X(\omega) &= Z(\omega)^\top \Sigma^{1/2}\Sigma^{-1}\Sigma^{1/2}Z(\omega) \\ &= Z(\omega)^\top I_p Z(\omega) \\ &= \sum_{j=1}^{p} Z_j(\omega)^2, \end{align*} where $Z_j: \Omega \to \mathbb{R}$ denotes the $j$-th coordinate random variable of $Z$. [/step] [step:Identify the components of $Z$ as independent standard normal random variables] Since $Z \sim \mathcal{N}_p(0,I_p)$, its density is \begin{align*} f_Z: \mathbb{R}^p &\to [0,\infty) \\ z &\mapsto (2\pi)^{-p/2}\exp\left(-\frac{1}{2}z^\top z\right). \end{align*} Writing $z=(z_1,\dots,z_p)$, this factors as \begin{align*} f_Z(z) &= \prod_{j=1}^{p}\left((2\pi)^{-1/2}\exp\left(-\frac{1}{2}z_j^2\right)\right). \end{align*} The factor corresponding to coordinate $j$ is the density of a one-dimensional standard normal random variable. Hence $Z_1,\dots,Z_p$ are independent and each satisfies $Z_j \sim \mathcal{N}(0,1)$. [/step] [step:Conclude that the quadratic form is chi-squared] By definition, if $Y_1,\dots,Y_p$ are independent random variables with $Y_j \sim \mathcal{N}(0,1)$ for each $j$, then \begin{align*} \sum_{j=1}^{p}Y_j^2 \sim \chi_p^2. \end{align*} Applying this definition to $Y_j=Z_j$ gives \begin{align*} \sum_{j=1}^{p}Z_j^2 \sim \chi_p^2. \end{align*} From the identity established above, \begin{align*} X^\top\Sigma^{-1}X=\sum_{j=1}^{p}Z_j^2, \end{align*} so \begin{align*} X^\top\Sigma^{-1}X \sim \chi_p^2. \end{align*} This proves the theorem. [/step]