Let $Y$ be a $n\times p$ random matrix with expected value $\mathbb{E}(Y)=M$ and variance $\mathbb{V}(vec(Y))=V \otimes U$, with $V$ and $U$ positive semidefinite $p\times p$ and $n\times n$ matrices, respectively. I found in an article that $$ \mathbb{E} \left(Y \otimes Y\right)= \operatorname{vec}(U) \operatorname{vec}(V)' + M \otimes M $$ but the authors do not provide a proof and I cannot understand how this relationship can be shown. How can I calculate this exectation?
I found a similar question here but it involves random vectors instead of matrices.
Let $M\equiv\mathsf{E}[Y\otimes Y]$. Then, the $((i-1)p+j)$-th column of $M$ can be written as \begin{align} [M]_{(i-1)p+j}&=\operatorname{vec}(\mathsf{E} [Y]_j[Y]_i^{\top})=\operatorname{vec}(\operatorname{Cov}([Y]_j,[Y]_i)+\mathsf{E}[Y]_j\mathsf{E}[Y]_i^{\top}) \\ &=\operatorname{vec}(V_{ij}U)+\operatorname{vec}(\mathsf{E}[Y]_j\mathsf{E}[Y]_i^{\top}) \\ &=[\operatorname{vec}(U)\operatorname{vec}(V)^{\top}]_{(i-1)p+j}+[M\otimes M]_{(i-1)p+j}. \end{align}