What is the product of two multivariate normal random variables with different length? Suppose $X$ is a $p$-dimensional multivariate normal random variable with conditional distribution
$$f_X (x\mid Z=z) = \mathcal{N}(\mu + Az, I\sigma^2),$$
Further suppose that
$$Z\sim \mathcal{N}(0, I),$$
where $Z$ is a $q$-dimensional multivariate normal random variable. (So in fact, $A$ is a $p\times q$ matrix.)
What is the joint distribution of $X$ and $Z$? In other words, what is the pdf
$$f_{X_1,\cdots,X_p,Z_1,\cdots,Z_q}(x_1,\cdots,x_p,z_1,\cdots,z_q)?$$
We can assume all the independence that we want (i.e. all $x_i$'s and $z_j$'s are independent)
I spent quite some time googling but all the examples I came across seem to implicitly assume that $X$ and $Y$ have the same dimension. For example, in this Matrix Cookbook page 41, it talks about "product of Gaussian density" but assumes that two input vectors need to have same dimensions.
 A: $$
X \mid Z=z \sim \mathcal{N}(\mu + Az, \sigma^2 I)
$$
$$
(X- AZ) \mid (Z=z) \sim\mathcal{N}(\mu,\sigma^2 I)
$$
Since the conditional distribution of $X-AZ$ given $Z$ does not depend on $Z,$ we must conclude that (1) $X-AZ$ and $Z$ are independent, and (2) the marginal (or "unconditional") distribution of $X-AZ$ is that same distribution.
Thus we have
$$
X-AZ\sim\mathcal N(\mu,\sigma^2 I)
$$
and
$$
(X-Az) \text{ and } Z \text{ are independent.}
$$
Thus
$$
\begin{bmatrix} X-AZ \\ Z \end{bmatrix} \sim \mathcal N \left( \begin{bmatrix} \mu \\ 0 \end{bmatrix}, \begin{bmatrix} \sigma^2 I & 0 \\ 0 & I \end{bmatrix} \right).
$$
$$
\begin{bmatrix} X \\ Z \end{bmatrix} = \begin{bmatrix} I & A \\ 0 & I \end{bmatrix} \begin{bmatrix} X-AZ \\ Z \end{bmatrix}
$$
Therefore $\begin{bmatrix} X \\ Z \end{bmatrix}$ is multivariate normal and
$$
\operatorname{E}\begin{bmatrix} X \\ Z \end{bmatrix} = \begin{bmatrix} I & A \\ 0 & I \end{bmatrix} \operatorname{E} \begin{bmatrix} X-AZ \\ Z \end{bmatrix} = \begin{bmatrix} I & A \\ 0 & I \end{bmatrix} \begin{bmatrix} \mu \\ 0 \end{bmatrix} = \begin{bmatrix} \mu \\ 0 \end{bmatrix},
$$
and
$$
\operatorname{var} \begin{bmatrix} X \\ Z \end{bmatrix} = \begin{bmatrix} I & A \\ 0 & I \end{bmatrix} \begin{bmatrix} \sigma^2 I & 0 \\ 0 & I \end{bmatrix} \begin{bmatrix} I & 0 \\ A^T & I \end{bmatrix} = \begin{bmatrix} \sigma^2 I + A A^T & A \\ A^T & I \end{bmatrix}.  
$$
