What is the analytic expression for PDF of joint distribution of two Gaussian random vectors? I know that if $X$ and $Y$ are random variables with respective PDFs,
$$
f_X(x) = \frac{1}{\sqrt{2\pi\sigma_x^2}}\exp\left\{-\frac{\left(x-\mu_x\right)^2}{2\sigma_x^2}\right\} \\
f_Y(y) = \frac{1}{\sqrt{2\pi\sigma_y^2}}\exp\left\{-\frac{\left(y-\mu_y\right)^2}{2\sigma_y^2}\right\}
$$
Then their joint PDF is written as
$$
f_{XY}(x,y) =
      \frac{1}{2 \pi  \sigma_x \sigma_y \sqrt{1-\rho^2}}
      \exp\left(
        -\frac{1}{2(1-\rho^2)}\left[
          \frac{(x-\mu_x)^2}{\sigma_x^2} -
          \frac{2\rho(x-\mu_x)(y-\mu_y)}{\sigma_x \sigma_y} +
          \frac{(y-\mu_y)^2}{\sigma_y^2}
        \right]
      \right)
$$
But when $\mathbf{x}$ and $\mathbf{y}$ are random vectors with PDFs
$$
f_{\mathbf x}(x_1,\ldots,x_k)\, =
\frac{1}{(2\pi)^{k/2}|\boldsymbol\Sigma_x|^{1/2}}
\exp\left(-\frac{1}{2}({\mathbf x}-{\boldsymbol\mu_x})^T{\boldsymbol\Sigma_x}^{-1}({\mathbf x}-{\boldsymbol\mu_x})
\right) \\
f_{\mathbf y}(y_1,\ldots,y_k)\, =
\frac{1}{(2\pi)^{k/2}|\boldsymbol\Sigma_y|^{1/2}}
\exp\left(-\frac{1}{2}({\mathbf y}-{\boldsymbol\mu_y})^T{\boldsymbol\Sigma_y}^{-1}({\mathbf y}-{\boldsymbol\mu_y})
\right)
$$
How do you express their joint PDF?
$$
f_{\mathbf xy}(\mathbf {x,y})\, = \, ?
$$
 A: In the case where you only assume that $\mathbf{X}$ and $\mathbf{Y}$ are marginally Gaussian, you can't say much about the joint density of $(\mathbf{X},\mathbf{Y})$, and you certainly can't conclude that the joint density is a Gaussian density. In the answer below I've added the additional assumption that the joint distribution is indeed Gaussian.
Assume that $\mathbf{X}=(X_1,\ldots,X_k)\sim\mathcal{N}_k(\boldsymbol{\mu}_x,\mathbf{\Sigma}_x)$ and $\mathbf{Y}=(Y_1,\ldots,Y_k)\sim\mathcal{N}_k(\boldsymbol{\mu}_y,\mathbf{\Sigma}_y)$ and that $\mathbf{Z}=(\mathbf{X},\mathbf{Y})$ follows a Gaussian distribution, then the joint density is given by
$$
f_\mathbf{Z}(\mathbf{z})=
\frac{1}{(2\pi)^{2k/2}|\boldsymbol\Sigma_z|^{1/2}}
\exp\left(-\frac{1}{2}({\mathbf z}-{\boldsymbol\mu_z})^T{\boldsymbol\Sigma_z}^{-1}({\mathbf z}-{\boldsymbol\mu_z})
\right),\quad\mathbf{z}\in\mathbb{R}^{2k},
$$
where
$$
\boldsymbol{\mu_z}=(\boldsymbol{\mu}_x,\boldsymbol{\mu}_y)
$$
and
$$
\mathbf{\Sigma}_z=\begin{bmatrix}
\mathbf{\Sigma}_x & \mathbf{\Sigma}_{xy} \\
\mathbf{\Sigma}_{yx} & \mathbf{\Sigma}_y
\end{bmatrix}
$$
written in block form. Here $\mathbf{\Sigma}_{xy}$ is the $k\times k$ matrix whose $(i,j)$th entry is $$
\mathbf{\Sigma}_{xy}^{i,j}=\mathrm{Cov}(X_i,Y_j).
$$
Note that this also holds when $\mathbf{X}$ is $k$-dimensional and $\mathbf{Y}$ is $n$-dimensional for $n\neq k$.
