I am working in a multivariate linear regression setting with the following model: $$Y_{n \times q} = X_{n \times p} B_{p\times q} + \mathcal{E}_{n \times q}.$$
Let $\mathcal{E}$ have variance-covariance matrix $\Sigma$.
Each row of $Y$ represents one observation of a $q$-dimensional random variable (so we have $n$ observations total). Suppose that the distribution of each row of $Y$ is independent and identically distributed as $N_q(XB, \Sigma)$ (the $q$ dimensional multivariate normal distribution).
Then the distribution of $Y$ given $X$, $B$, and $\Sigma$ is proportional to the following: $$f(T|B, X, \Sigma)\propto \left(\frac{1}{\sqrt{|{\Sigma}|}}\right)^n\text{exp}\left\{-\frac{1}{2}\sum_{m = 1}^n ({Y}_m - {X}_m{B}){\Sigma}^{-1} ({Y}_m - {X}_m{B})^{T} \right\}$$
I am trying to write the portion in the exponential function in terms of the vec operator and Kronecker products. It seems like it should be pretty simple, but for some reason I cannot see it.
Ultimately, of course, I am trying to reach a connection between the summation in the exponential and the least squares estimate of $B$ given by $$\textbf{vec}(\widehat{B}) = (I_q \otimes (X^TX)^{-1}X^T)\textbf{vec}(Y)$$
Is there a way to write $\sum_{m = 1}^n ({Y}_m - {X}_m{B}){\Sigma}^{-1} ({Y}_m - {X}_m{B})^{T}$ as a combination of Kronecker products and the vec operator?