Suppose I have a $p \times n$ (with $p \geq n$) matrix $\bf U$ such that ${\bf U}' {\bf U} = {\bf I}_{n}$ and that $\bf U$ is uniformly distributed on the Stiefel manifold $V_{n,p}$. I would like to know the expected values $\Bbb{E} [ U_{ij} U_{k\ell} ]$ because I need to compute $\Bbb{E} [ \operatorname{tr}( {\bf U} {\bf A} {\bf U} {\bf B}) ]$ and $\Bbb E [ \operatorname{tr}( {\bf U} {\bf A} {\bf U}' {\bf B}) ]$.
I know that $\Bbb{E} [{\bf U}]=0$, $\Bbb E[{\bf U} {\bf U}']=\frac{n}{p} {\bf I}_{p}$ and the moments $\Bbb{E} [ P_{ij} P_{k \ell}]$ with ${\bf P} = {\bf U} {\bf U}'$ are given in Mardia & Khatri JMA 1977 but I need $\Bbb{E} [ U_{ij} U_{k\ell} ]$. Is the latter $\frac{1}{p}\delta_{ik}\delta_{j\ell}$?
I think I have a partial answer: ${\bf U}$ has a spherical distribution and hence has the same distribution as ${\bf \Gamma U \Omega}$ for any orthogonal matrices ${\bf \Gamma} \in O(p)$ and ${\bf \Omega} \in O(n)$, where $O(p)= \{ {\bf{Q}} (p \times p); \; {\bf{Q}}'{\bf{Q}}={\bf{I}}_{p} \}$ is the orthogonal group. Let $\operatorname{vec}({\bf U})$ be the operator that stacks the columns of ${\bf U}$. Then the covariance matrix ${\bf C} = \Bbb{E} [\operatorname{vec}({\bf U}) \operatorname{vec}({\bf U})']$ satisfies ${\bf C} = \Bbb{E} [ \operatorname{vec}({\bf \Gamma U \Omega}) \operatorname{vec}({\bf \Gamma U \Omega})'] = ({\bf \Omega}' \otimes {\bf \Gamma}) {\bf C} ({\bf \Omega} \otimes {\bf \Gamma}')$. Therefore ${\bf C}$ commutes with any matrix of the form ${\bf \Omega}' \otimes {\bf \Gamma}$ with ${\bf \Gamma} \in O(p)$ and ${\bf \Omega} \in O(n)$. Note that ${\bf \Omega}' \otimes {\bf \Gamma} \in O(pn)$. However it does not mean that ${\bf C}$ commutes with every matrix in the orthogonal group $O(pn)$, in which case I think I could assert that ${\bf C}$ is a scalar times the identity matrix. So I would need to prove that if ${\bf C}$ commutes with every ${\bf \Omega}' \otimes {\bf \Gamma}$ then necessarily ${\bf C} = \alpha {\bf I}_{pn}$. If I can do that I am done since $\operatorname{tr}({\bf C}) = \Bbb{E} [\operatorname{tr}(\operatorname{vec}({\bf U}) \operatorname{vec}({\bf U})')] = \Bbb{E} [\operatorname{vec}({\bf U})' \operatorname{vec}({\bf U})] = \Bbb{E} [\operatorname{tr}({\bf U}'{\bf U})] = \Bbb{E} [\operatorname{tr}({\bf I}_{n})]=n$ and hence ${\bf C} = \frac{1}{p} {\bf I}_{pn}$ which implies that $\Bbb{E} [U_{ij} U_{k\ell} ]=\frac{1}{p}\delta_{ik}\delta_{j\ell}$.