Say we have a mixed-effects model with a single grouping factor, indexed with $i$:

$$ y_i = X_i\beta + Z_ib_i + \epsilon_i \\ \epsilon_i \sim \mathcal N(0, \sigma^2) \\ b_i \sim \mathcal N(0, \Sigma) $$

In this setting, we get a simple closed-form likelihood for each group, because each $y_i$ is multivariate normal:

$$ y_i \sim \mathcal N(X_i\beta,\ Z_i \Sigma Z^T_i + \sigma^2I) $$

In the above the $y_i$s are one-dimensional vectors -- i.e., our model has a univariate response.

Question: What is the likelihood when we instead have a multivariate response?

For example, if we have a bivariate response:

$$ y_i = (Y^T_{i1},Y^T_{i2})^T \\ y_{i1} = X_{i1}B_1 + Z_{i1}b_{i1} + \epsilon_{i1} \\ y_{i2} = X_{i2}B_2 + Z_{i2}b_{i2} + \epsilon_{i2} \\ \epsilon_i \sim \mathcal N(0, \begin{bmatrix} \sigma_1^2 & \rho\sigma_1\sigma_2 \\ \rho\sigma_1\sigma_2 & \sigma_2^2 \end{bmatrix}) \\ b_i \sim \mathcal N(0, \begin{bmatrix} \Sigma_{11} & \Sigma_{12} \\ \Sigma_{12} & \Sigma_{22} \end{bmatrix}) \\ $$

Then I believe $y_i$ should still be multivariate normal, but I don't quite understand how to calculate the covariance:

$$ y_i \sim \mathcal N(X_i\beta,?) $$