Result for covariance between elements of the sample covariance matrix Given data matrix $X_{n,p}$ where $n$ is the sample size, $p$ is the dimension, the sample covariance is 
$$\hat{\Sigma}=\frac{1}{n-1}X^\top X=\{\hat{\sigma}_{ij}\}_{1\le i,j\le p}$$
Is there any result on the quantity like $cov(\hat{\sigma}_{i_0,j_0},\hat{\sigma}_{i_0,j_0+1})$? I.e. what is the distribution? Or asymptotic distribution? And for even more general terms like $cov(\hat{\sigma}_{i_1,j_1},\hat{\sigma}_{i_2,j_2})$? We can assume Gaussian distribution first, but not necessarily. 
 A: If you're assuming Normality, you want to have a look at the Wishart distribution.
If $\mathbf{X}_{n,p}$ ($n \ge p$) has rows that are iid multivariate normal with mean $0$ and variance $\Sigma$, then
$$
\mathbf{X}^T\mathbf{X} \sim \text{Wishart}_p(n,\Sigma).
$$
This is for the particular case you have where the mean of each row is the zero vector. Also


*

*$\widehat{\Sigma} \sim \text{Wishart}_p(n,\frac{1}{n-1}\Sigma)$ so you probably want to divide by $n$ rather then $n-1$ because $E\hat{\Sigma} = \frac{n}{n-1}\Sigma$. This is biased.

*The $\text{vec}$ operator the kronecker product and all their properties will come in handy for looking at higher order moments. Closed form expressions exist for all covariances and variances.


If you want to take a closer look at the formula for all the variances and covariances, it helps to write $\mathbf{X} = \mathbf{Z}\mathbf{C}$, where $\mathbf{C}$ is the Cholesky decomposition of $\Sigma = \mathbf{C}^T\mathbf{C}$. The variance is
\begin{align*}
\text{Var}\left( \text{vec} \left[\mathbf{X}^T\mathbf{X} \right]\right) &=
\text{Var}\left( \text{vec} \left[\mathbf{C}^T\mathbf{Z}^T\mathbf{Z} \mathbf{C} \right]\right) \\
&= \text{Var}\left(\left[\mathbf{C}^T\otimes \mathbf{C}^T \right] \text{vec} \left[\mathbf{Z}^T\mathbf{Z} \right]\right)\tag{prop.s of vec/kron} \\
&= \left[\mathbf{C}^T\otimes \mathbf{C}^T \right] \text{Var}(\text{vec}(\mathbf{Z}^T\mathbf{Z})) \left[\mathbf{C}^T\otimes \mathbf{C}^T \right] ^T \tag{properties of var}\\
&= \left[\mathbf{C}^T\otimes \mathbf{C}^T \right] \left[I_p \otimes I_p + M_p \right] \left[\mathbf{C}\otimes \mathbf{C} \right] 
\end{align*}
where $\mathbf{M}_p$ is a $p^2 \times p^2$ matrix that is described further in this document. If you want the variance of $\hat{\Sigma}$, just scale the last expression by $(n-1)^{-2}$.
