When calculating the vector-matrix-vector product $$\mathbb{E}\left[x_1^TWx_2\right],$$ where $x_1$ and $x_2$ are $n \times 1$ vectors of random variables and $W$ is a $n \times n$ constant matrix, one can use the trace to "factor out" $W$ in the following way $$\begin{align}\mathbb{E}\left[x_1^TWx_2\right] &= \text{Tr}\left(\mathbb{E}\left[x_1^TWx_2\right]\right)\\ &= \mathbb{E}\left[\text{Tr}\left(x_1^TWx_2\right)\right]\\ &= \mathbb{E}\left[\text{Tr}\left(Wx_2x_1^T\right)\right]\\ &= \text{Tr}\left(\mathbb{E}\left[Wx_2x_1^T\right]\right)\\ &= \text{Tr}\left(W\mathbb{E}\left[x_2x_1^T\right]\right). \end{align}$$
However, if we want to calculate $$\mathbb{E}\left(X_1^TWX_2\right)$$ where now $X_1$ and $X_2$ are random matrices of shape $n \times m$, it is not possible to apply the trace operator to factor out $W$ anymore, since the expected value is now a matrix and not a scalar.
How can we calculate this expected value?