[proofplan] We center the sample by writing $Y_i := X_i-\mu$, so that each $Y_i$ has mean zero and covariance matrix $\Sigma$. Linearity of expectation gives the unbiasedness of $\bar X$, and independence makes all cross-covariance terms vanish in the covariance of $\bar X$. For the sample covariance, we rewrite each residual as $X_i-\bar X = Y_i - n^{-1}\sum_{j=1}^n Y_j$, expand the outer products, and use the same vanishing of cross terms to show that the expected numerator is $(n-1)\Sigma$. [/proofplan] [step:Center the random sample and record the second moment identities] For each $1 \leq i \leq n$, define the centered random vector \begin{align*} Y_i : \Omega &\to \mathbb{R}^p \\ \omega &\mapsto X_i(\omega)-\mu. \end{align*} Since $X_i \sim \mathcal{N}_p(\mu,\Sigma)$, each $Y_i$ is integrable, has finite second moments, and satisfies \begin{align*} \mathbb{E}[Y_i] &= 0, & \mathbb{E}[Y_iY_i^\top] &= \Sigma. \end{align*} Because the random vectors $X_1,\dots,X_n$ are independent, the centered random vectors $Y_1,\dots,Y_n$ are also independent. Hence, for $i \neq j$, \begin{align*} \mathbb{E}[Y_iY_j^\top] = \mathbb{E}[Y_i]\mathbb{E}[Y_j]^\top = 0. \end{align*} [/step] [step:Compute the expectation and covariance matrix of the sample mean] Using the definition \begin{align*} \bar X = \frac{1}{n}\sum_{i=1}^{n}X_i, \end{align*} linearity of expectation gives \begin{align*} \mathbb{E}[\bar X] &= \frac{1}{n}\sum_{i=1}^{n}\mathbb{E}[X_i] \\ &= \frac{1}{n}\sum_{i=1}^{n}\mu \\ &= \mu. \end{align*} Also, \begin{align*} \bar X-\mu = \frac{1}{n}\sum_{i=1}^{n}Y_i. \end{align*} Therefore \begin{align*} \operatorname{Cov}(\bar X) &= \mathbb{E}\big[(\bar X-\mu)(\bar X-\mu)^\top\big] \\ &= \mathbb{E}\left[\left(\frac{1}{n}\sum_{i=1}^{n}Y_i\right)\left(\frac{1}{n}\sum_{j=1}^{n}Y_j\right)^\top\right] \\ &= \frac{1}{n^2}\sum_{i=1}^{n}\sum_{j=1}^{n}\mathbb{E}[Y_iY_j^\top]. \end{align*} The diagonal terms are $\Sigma$, and the off-diagonal terms vanish by independence and centering. Thus \begin{align*} \operatorname{Cov}(\bar X) &= \frac{1}{n^2}\sum_{i=1}^{n}\Sigma \\ &= \frac{1}{n}\Sigma. \end{align*} [guided] The sample mean is an average of the random vectors $X_i$, so its expectation is computed by linearity of expectation. Since each $X_i$ has mean $\mu$, we obtain \begin{align*} \mathbb{E}[\bar X] &= \mathbb{E}\left[\frac{1}{n}\sum_{i=1}^{n}X_i\right] \\ &= \frac{1}{n}\sum_{i=1}^{n}\mathbb{E}[X_i] \\ &= \frac{1}{n}\sum_{i=1}^{n}\mu \\ &= \mu. \end{align*} For the covariance matrix, we must compute the second centered moment of $\bar X$. Since \begin{align*} \bar X-\mu = \frac{1}{n}\sum_{i=1}^{n}(X_i-\mu)=\frac{1}{n}\sum_{i=1}^{n}Y_i, \end{align*} we expand the outer product: \begin{align*} \operatorname{Cov}(\bar X) &= \mathbb{E}\big[(\bar X-\mu)(\bar X-\mu)^\top\big] \\ &= \mathbb{E}\left[\left(\frac{1}{n}\sum_{i=1}^{n}Y_i\right)\left(\frac{1}{n}\sum_{j=1}^{n}Y_j\right)^\top\right] \\ &= \frac{1}{n^2}\sum_{i=1}^{n}\sum_{j=1}^{n}\mathbb{E}[Y_iY_j^\top]. \end{align*} There are two kinds of terms. When $i=j$, the term is \begin{align*} \mathbb{E}[Y_iY_i^\top]=\Sigma. \end{align*} When $i \neq j$, independence gives \begin{align*} \mathbb{E}[Y_iY_j^\top] = \mathbb{E}[Y_i]\mathbb{E}[Y_j]^\top = 0. \end{align*} Thus only the $n$ diagonal terms remain: \begin{align*} \operatorname{Cov}(\bar X) &= \frac{1}{n^2}\sum_{i=1}^{n}\Sigma \\ &= \frac{1}{n}\Sigma. \end{align*} [/guided] [/step] [step:Expand the sample covariance numerator in centered variables] Define the centered sum random vector \begin{align*} A : \Omega &\to \mathbb{R}^p \\ \omega &\mapsto \sum_{i=1}^{n}Y_i(\omega). \end{align*} Then \begin{align*} \bar X-\mu = \frac{1}{n}A, \end{align*} and, for each $1 \leq i \leq n$, \begin{align*} X_i-\bar X = Y_i-\frac{1}{n}A. \end{align*} Therefore the numerator of $S$ satisfies \begin{align*} \sum_{i=1}^{n}(X_i-\bar X)(X_i-\bar X)^\top &= \sum_{i=1}^{n}\left(Y_i-\frac{1}{n}A\right)\left(Y_i-\frac{1}{n}A\right)^\top \\ &= \sum_{i=1}^{n}Y_iY_i^\top -\frac{1}{n}\sum_{i=1}^{n}Y_iA^\top -\frac{1}{n}\sum_{i=1}^{n}AY_i^\top +\frac{1}{n^2}\sum_{i=1}^{n}AA^\top. \end{align*} Since $\sum_{i=1}^{n}Y_i=A$, this becomes \begin{align*} \sum_{i=1}^{n}(X_i-\bar X)(X_i-\bar X)^\top &= \sum_{i=1}^{n}Y_iY_i^\top -\frac{1}{n}AA^\top -\frac{1}{n}AA^\top +\frac{1}{n}AA^\top \\ &= \sum_{i=1}^{n}Y_iY_i^\top-\frac{1}{n}AA^\top. \end{align*} [guided] The residual $X_i-\bar X$ is easier to handle after centering at the true mean $\mu$. Define \begin{align*} A : \Omega &\to \mathbb{R}^p \\ \omega &\mapsto \sum_{i=1}^{n}Y_i(\omega). \end{align*} Then the sample mean satisfies \begin{align*} \bar X-\mu = \frac{1}{n}A, \end{align*} so each residual can be rewritten as \begin{align*} X_i-\bar X &= (X_i-\mu)-(\bar X-\mu) \\ &= Y_i-\frac{1}{n}A. \end{align*} Now expand the outer product exactly: \begin{align*} \sum_{i=1}^{n}(X_i-\bar X)(X_i-\bar X)^\top &= \sum_{i=1}^{n}\left(Y_i-\frac{1}{n}A\right)\left(Y_i-\frac{1}{n}A\right)^\top \\ &= \sum_{i=1}^{n}Y_iY_i^\top -\frac{1}{n}\sum_{i=1}^{n}Y_iA^\top -\frac{1}{n}\sum_{i=1}^{n}AY_i^\top +\frac{1}{n^2}\sum_{i=1}^{n}AA^\top. \end{align*} The sums involving $A$ simplify because $A=\sum_{i=1}^{n}Y_i$. Hence \begin{align*} \sum_{i=1}^{n}Y_iA^\top &= AA^\top, \\ \sum_{i=1}^{n}AY_i^\top &= AA^\top, \\ \sum_{i=1}^{n}AA^\top &= nAA^\top. \end{align*} Substituting these identities gives \begin{align*} \sum_{i=1}^{n}(X_i-\bar X)(X_i-\bar X)^\top &= \sum_{i=1}^{n}Y_iY_i^\top -\frac{1}{n}AA^\top -\frac{1}{n}AA^\top +\frac{1}{n}AA^\top \\ &= \sum_{i=1}^{n}Y_iY_i^\top-\frac{1}{n}AA^\top. \end{align*} This identity is the algebraic reason for the denominator $n-1$ in the unbiased sample covariance matrix. [/guided] [/step] [step:Take expectations to obtain the unbiasedness of the sample covariance matrix] Taking entrywise expectations in the identity from the previous step gives \begin{align*} \mathbb{E}\left[\sum_{i=1}^{n}(X_i-\bar X)(X_i-\bar X)^\top\right] &= \sum_{i=1}^{n}\mathbb{E}[Y_iY_i^\top] -\frac{1}{n}\mathbb{E}[AA^\top]. \end{align*} The first term is \begin{align*} \sum_{i=1}^{n}\mathbb{E}[Y_iY_i^\top]=n\Sigma. \end{align*} For the second term, expand $AA^\top$: \begin{align*} \mathbb{E}[AA^\top] &= \mathbb{E}\left[\left(\sum_{i=1}^{n}Y_i\right)\left(\sum_{j=1}^{n}Y_j\right)^\top\right] \\ &= \sum_{i=1}^{n}\sum_{j=1}^{n}\mathbb{E}[Y_iY_j^\top] \\ &= n\Sigma. \end{align*} Thus \begin{align*} \mathbb{E}\left[\sum_{i=1}^{n}(X_i-\bar X)(X_i-\bar X)^\top\right] &= n\Sigma-\frac{1}{n}n\Sigma \\ &= (n-1)\Sigma. \end{align*} Since $n \geq 2$, division by $n-1$ is valid. By the definition of $S$, \begin{align*} \mathbb{E}[S] &= \frac{1}{n-1}\mathbb{E}\left[\sum_{i=1}^{n}(X_i-\bar X)(X_i-\bar X)^\top\right] \\ &= \frac{1}{n-1}(n-1)\Sigma \\ &= \Sigma. \end{align*} Combining this with the previous computation of $\mathbb{E}[\bar X]$ and $\operatorname{Cov}(\bar X)$ proves all three asserted identities. [/step]