Let $x_1, \ldots x_n$ be iid realizations of a $p$-dimensional random column vector $X= (X_1, \ldots, X_p)$ such that $X \sim N_p ( \mu ,\Sigma)$. We can show that $$ \hat{\mu} = \frac{1}{n} \sum_{i=1}^{n} x_i$$ $$\hat{\large \Sigma} = \frac{1}{n} \sum_{i=1}^{n} (x_i - \hat{\mu}) ( x_i - \hat{\mu})^T$$ Denote $x_i = (x_{i1}, \ldots , x_{ip})$, where $x_{ik} \sim X_k$ for $1 \le k \le p$ and $1 \le i \le n$. Then the $ij$-th element is $${\large\hat{\Sigma} _{ij}}= \frac{1}{n} \sum_{l=1}^{n} (x_{li} - \hat{\mu}_i)(x_{lj} - \hat{\mu}_j) \quad \quad (*) $$ Compute $E(\hat{\Sigma}_{ij})$.
So I attempt to compute the expected value of the $ij$th term using properties of covariance. By using the definitions and linear properties $$ Cov(X,Y) = E[(X-E(X))(Y-E(Y))] $$ $$ Cov(X,Y+Z)=Cov(X,Y)+Cov(X,Z)$$ $$ Cov(aX,Y) = aCov(X,Y)$$ and the fact that $$ E(x_{li} - \hat{\mu}_i )= E(x_{lj} - \hat{\mu}_j) = 0 $$ we obtain \begin{align*} E[(x_{li} - \hat{\mu}_i)(x_{lj} - \hat{\mu}_j)] & = Cov[(x_{li} - \hat{\mu}_i, x_{lj} - \hat{\mu}_j)] \\ & = Cov(x_{li}, x_{lj}) - \frac{1}{n}*n Cov(x_{li}, x_{lj}) - \frac{1}{n}*n Cov(x_{li}, x_{lj}) + \frac{1}{n^2} * n^2 Cov_(x_{li}, x_{lj}) \\ &=0 \end{align*} for $1 \le l \le n$.
Hence, $E( \hat{\Sigma}_{ij}) = 0$ from $(*)$?
This did not fit my calculations when I did without using covariance. Or even from internet research, saying that there is a factor. It should have equal to $(n-1)\Sigma_{ij}$. What went wrong?