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?
There is an error in your calculations involving $\hat\mu_j$ and $\hat\mu_i$ (i.e., the last three terms in your expansion). Using a non-clashing index of summation, we have $$ \hat\mu_j=\frac1n\sum_kx_{kj}, $$ so that $${\rm Cov}(x_{li},\hat\mu_j)=\frac1n\sum_k {\rm Cov}(x_{li},x_{kj}).\tag1$$ The terms in the sum (1) are zero when $k\ne l$ (by independence), and are equal to ${\rm Cov}(x_{li},x_{lj})=:\Sigma_{i,j}$ only when $k=l$. So $${\rm Cov}(x_{li},\hat\mu_j)=\frac1n\Sigma_{ij}.$$ By similar reasoning we obtain $${\rm Cov}(x_{lj},\hat\mu_i)=\frac1n\Sigma_{ij}$$ and $${\rm Cov}(\hat\mu_i,\hat\mu_j)=\frac1n\Sigma_{ij}.$$ Putting this all together you should obtain $$E(\hat\Sigma_{ij})=\frac{n-1}n\Sigma_{ij}.$$