[proofplan] Represent the Wishart matrix as a sum of independent rank-one matrices $X_rX_r^\top$. The expectation follows by taking expectations entrywise and using the covariance matrix of a centred Gaussian vector. For the covariance, expand each entry $w_{ij}$ as a sum over observations; independence removes all cross-observation covariance terms, and the remaining single-observation fourth moment is evaluated by Isserlis' formula. The variance statement is the special case $(k,l)=(i,j)$ of the covariance formula. [/proofplan]

[step:Expand the Wishart entries as sums of Gaussian products] For each $r \in \{1,\dots,n\}$, the outer product $X_rX_r^\top$ is the $p \times p$ random matrix whose $(i,j)$ entry is $X_{ri}X_{rj}$. Hence, for every $i,j \in \{1,\dots,p\}$, \begin{align*} w_{ij} = \sum_{r=1}^n X_{ri}X_{rj}. \end{align*} Since $X_r \sim \mathcal{N}_p(0,\Sigma)$, its coordinate means and covariances satisfy \begin{align*} \mathbb{E}[X_{ri}] &= 0,\\ \mathbb{E}[X_{ri}X_{rj}] &= \sigma_{ij}. \end{align*} [/step]

[step:Compute the expectation entry by entry] Fix $i,j \in \{1,\dots,p\}$. By linearity of expectation applied to the finite sum defining $w_{ij}$, \begin{align*} \mathbb{E}[w_{ij}] &= \mathbb{E}\left[\sum_{r=1}^n X_{ri}X_{rj}\right] \\ &= \sum_{r=1}^n \mathbb{E}[X_{ri}X_{rj}] \\ &= \sum_{r=1}^n \sigma_{ij} \\ &= n\sigma_{ij}. \end{align*} Since this holds for every matrix entry, the matrix identity is \begin{align*} \mathbb{E}[W] = n\Sigma. \end{align*} [/step]

[step:Reduce the covariance to one Gaussian observation]Fix $i,j,k,l \in \{1,\dots,p\}$. Define real-valued random variables \begin{align*} A_r &: \Omega \to \mathbb{R}, & A_r(\omega) &= X_{ri}(\omega)X_{rj}(\omega),\\ B_r &: \Omega \to \mathbb{R}, & B_r(\omega) &= X_{rk}(\omega)X_{rl}(\omega), \end{align*} for $r \in \{1,\dots,n\}$. Then \begin{align*} w_{ij} = \sum_{r=1}^n A_r, \qquad w_{kl} = \sum_{r=1}^n B_r. \end{align*} By bilinearity of covariance over finite sums, \begin{align*} \operatorname{Cov}(w_{ij},w_{kl}) = \sum_{r=1}^n \sum_{s=1}^n \operatorname{Cov}(A_r,B_s). \end{align*} If $r \neq s$, then $X_r$ and $X_s$ are independent, so $A_r$ and $B_s$ are independent. Therefore \begin{align*} \operatorname{Cov}(A_r,B_s)=0 \qquad \text{for } r \neq s. \end{align*} Thus only the diagonal terms remain: \begin{align*} \operatorname{Cov}(w_{ij},w_{kl}) = \sum_{r=1}^n \operatorname{Cov}(A_r,B_r). \end{align*} Because the random vectors $X_1,\dots,X_n$ have the same distribution, the summands are equal, and hence \begin{align*} \operatorname{Cov}(w_{ij},w_{kl}) = n\,\operatorname{Cov}(X_{1i}X_{1j},X_{1k}X_{1l}). \end{align*}[/step]

[guided]The entry $w_{ij}$ is a sum over the $n$ independent observations, and the same is true of $w_{kl}$. To make the covariance expansion explicit, define \begin{align*} A_r &: \Omega \to \mathbb{R}, & A_r(\omega) &= X_{ri}(\omega)X_{rj}(\omega),\\ B_r &: \Omega \to \mathbb{R}, & B_r(\omega) &= X_{rk}(\omega)X_{rl}(\omega). \end{align*} Then \begin{align*} w_{ij} = \sum_{r=1}^n A_r, \qquad w_{kl} = \sum_{r=1}^n B_r. \end{align*} Using bilinearity of covariance for finite sums gives \begin{align*} \operatorname{Cov}(w_{ij},w_{kl}) = \sum_{r=1}^n \sum_{s=1}^n \operatorname{Cov}(A_r,B_s). \end{align*} The purpose of introducing $A_r$ and $B_s$ is to separate terms coming from the same observation from terms coming from different observations. If $r \neq s$, then $A_r$ is a measurable function of $X_r$ and $B_s$ is a measurable function of $X_s$. Since $X_r$ and $X_s$ are independent, $A_r$ and $B_s$ are independent, and hence \begin{align*} \operatorname{Cov}(A_r,B_s) = \mathbb{E}[A_rB_s]-\mathbb{E}[A_r]\mathbb{E}[B_s] = 0. \end{align*} Therefore all off-diagonal terms in the double sum vanish: \begin{align*} \operatorname{Cov}(w_{ij},w_{kl}) = \sum_{r=1}^n \operatorname{Cov}(A_r,B_r). \end{align*} Since $X_1,\dots,X_n$ are identically distributed, each covariance in the last sum has the same value. Thus \begin{align*} \operatorname{Cov}(w_{ij},w_{kl}) = n\,\operatorname{Cov}(X_{1i}X_{1j},X_{1k}X_{1l}). \end{align*}[/guided]

[step:Evaluate the single-observation fourth moment by Isserlis' formula]Let $Y : \Omega \to \mathbb{R}^p$ be the random vector $Y := X_1$, and write $Y = (Y_1,\dots,Y_p)$. Then $Y \sim \mathcal{N}_p(0,\Sigma)$, so \begin{align*} \mathbb{E}[Y_a] = 0, \qquad \mathbb{E}[Y_aY_b] = \sigma_{ab} \end{align*} for all $a,b \in \{1,\dots,p\}$. By Isserlis' formula for centred multivariate Gaussian fourth moments (citing a result not yet in the wiki: Isserlis' formula), applied to the four coordinates $Y_i,Y_j,Y_k,Y_l$, \begin{align*} \mathbb{E}[Y_iY_jY_kY_l] = \mathbb{E}[Y_iY_j]\mathbb{E}[Y_kY_l] + \mathbb{E}[Y_iY_k]\mathbb{E}[Y_jY_l] + \mathbb{E}[Y_iY_l]\mathbb{E}[Y_jY_k]. \end{align*} Substituting the covariance entries of $Y$ gives \begin{align*} \mathbb{E}[Y_iY_jY_kY_l] = \sigma_{ij}\sigma_{kl} + \sigma_{ik}\sigma_{jl} + \sigma_{il}\sigma_{jk}. \end{align*} Therefore \begin{align*} \operatorname{Cov}(Y_iY_j,Y_kY_l) &= \mathbb{E}[Y_iY_jY_kY_l] - \mathbb{E}[Y_iY_j]\mathbb{E}[Y_kY_l] \\ &= \left(\sigma_{ij}\sigma_{kl} + \sigma_{ik}\sigma_{jl} + \sigma_{il}\sigma_{jk}\right) - \sigma_{ij}\sigma_{kl} \\ &= \sigma_{ik}\sigma_{jl} + \sigma_{il}\sigma_{jk}. \end{align*}[/step]

[guided]We now need the covariance of two quadratic coordinate products from one centred Gaussian vector. Set $Y := X_1$, so $Y : \Omega \to \mathbb{R}^p$ is a random vector with distribution $\mathcal{N}_p(0,\Sigma)$, and write $Y=(Y_1,\dots,Y_p)$. The covariance matrix condition means \begin{align*} \mathbb{E}[Y_a] = 0, \qquad \mathbb{E}[Y_aY_b] = \sigma_{ab} \end{align*} for every $a,b \in \{1,\dots,p\}$. The relevant fourth moment is computed by Isserlis' formula for centred Gaussian random variables (citing a result not yet in the wiki: Isserlis' formula). Applied to the four centred Gaussian coordinates $Y_i,Y_j,Y_k,Y_l$, it gives \begin{align*} \mathbb{E}[Y_iY_jY_kY_l] = \mathbb{E}[Y_iY_j]\mathbb{E}[Y_kY_l] + \mathbb{E}[Y_iY_k]\mathbb{E}[Y_jY_l] + \mathbb{E}[Y_iY_l]\mathbb{E}[Y_jY_k]. \end{align*} Each expectation on the right is an entry of the covariance matrix $\Sigma$, so \begin{align*} \mathbb{E}[Y_iY_jY_kY_l] = \sigma_{ij}\sigma_{kl} + \sigma_{ik}\sigma_{jl} + \sigma_{il}\sigma_{jk}. \end{align*} Now subtract the product of the means of the two quadratic factors: \begin{align*} \operatorname{Cov}(Y_iY_j,Y_kY_l) &= \mathbb{E}[Y_iY_jY_kY_l] - \mathbb{E}[Y_iY_j]\mathbb{E}[Y_kY_l] \\ &= \left(\sigma_{ij}\sigma_{kl} + \sigma_{ik}\sigma_{jl} + \sigma_{il}\sigma_{jk}\right) - \sigma_{ij}\sigma_{kl} \\ &= \sigma_{ik}\sigma_{jl} + \sigma_{il}\sigma_{jk}. \end{align*} The cancellation of $\sigma_{ij}\sigma_{kl}$ is exactly the subtraction built into covariance.[/guided]

[step:Conclude the covariance and variance formulas] Combining the reduction to one observation with the single-observation covariance calculation gives \begin{align*} \operatorname{Cov}(w_{ij},w_{kl}) &= n\,\operatorname{Cov}(X_{1i}X_{1j},X_{1k}X_{1l}) \\ &= n(\sigma_{ik}\sigma_{jl}+\sigma_{il}\sigma_{jk}). \end{align*} This proves the asserted entrywise covariance formula. Taking $k=i$ and $l=j$ gives \begin{align*} \operatorname{Var}(w_{ij}) &= \operatorname{Cov}(w_{ij},w_{ij}) \\ &= n(\sigma_{ii}\sigma_{jj}+\sigma_{ij}\sigma_{ji}). \end{align*} Since $\Sigma$ is symmetric, $\sigma_{ji}=\sigma_{ij}$, and hence \begin{align*} \operatorname{Var}(w_{ij}) = n(\sigma_{ii}\sigma_{jj}+\sigma_{ij}^2). \end{align*} Together with $\mathbb{E}[W]=n\Sigma$, this completes the proof. [/step]