$(X, Y)$ PDF is in the Multivariate Normal form $\implies$ $(X, Y)$ multivariate normal

Definition. $\mathbf{X} = \begin{bmatrix}X \\ Y\end{bmatrix}$ is multivariate normal (equivalently, $X$ and $Y$ are jointly normal) if there is a vector $\boldsymbol\mu \in \mathbb{R}^2$, a $k \times 2$ matrix $\mathbf{A}$ (with $k \leq 2$), and a random vector $\mathbf{z} \in \mathbb{R}^k$ of $\mathcal{N}(0, 1)$ random variables such that $$\mathbf{X} = \boldsymbol\mu+\mathbf{A}^{T}\mathbf{z}$$ ($T$ denoting the transpose operator) with $\mathbf{A}^{T}\mathbf{A} = \boldsymbol\Sigma$. We write $\mathbf{X} \sim \text{MVN}_{2}(\boldsymbol\mu, \boldsymbol\Sigma)$ if this is the case.

Theorem. Suppose $\mathbf{X} \sim \text{MVN}_{2}(\boldsymbol\mu, \boldsymbol\Sigma)$ with $\boldsymbol\Sigma$ non-singular. Then $\mathbf{X}$ has density $$f_{\mathbf{X}}(\mathbf{x}) = \left[\dfrac{1}{(2\pi)^2\det(\boldsymbol\Sigma)} \right]^{1/2}\exp\left[-\dfrac{1}{2}(\mathbf{x}-\boldsymbol\mu)^{T}\boldsymbol\Sigma^{-1}(\mathbf{x}-\boldsymbol\mu) \right]\text{.}\tag{1}$$

Question. Suppose I have two random variables $X, Y$ with support in $\mathbb{R}^2$ such that their (variance-)covariance matrix is non-singular. I am able to write their joint density in the form of $(1)$ above. Is this sufficient to show that $\mathbf{X} = \begin{bmatrix} X \\ Y \end{bmatrix}$ is multivariate normal?

If more context is needed, let me know. I have not been able to find an answer to this in any searching I've done, and apologize that I haven't been able to provide much work for this. I understand how the Theorem above is proven: use transformation techniques, but the problem with this is that you need to know what the distribution of a random variable already is prior to the transformation, so you can't easily "reverse" this proof.

11/26: I've considered pursuing this by computing the MGF of $\mathbf{X}$ and showing that it matches the MGF of a multivariate normal distribution (invoking the theorem that MGFs uniquely determine the distribution), but this is seems to be a lot more work than it is worth.

But you do know the density of $(X, Y)$, which is all you need. I will write down the general proof for the $n$-dimensional normal distribution.

Since $\Sigma$ is a symmetric matrix, we can diagonalize it orthogonally, i.e. $\Sigma = Q^T \Lambda Q$, where $\Lambda = \operatorname{diag}(\lambda_1, \ldots, \lambda_n)$ are the positive eigenvalues of $\Sigma$. This allows us to define a square root of $\Sigma$ as $A := Q^T \Lambda^{1/2} Q$, where $\Lambda^{1/2} := \operatorname{diag}(\sqrt{\lambda_1}, \ldots, \sqrt{\lambda_n})$.

Note that $A$ is a symmetric matrix that satisfies $$A \Sigma^{-1} A = (Q^t \Lambda^{1/2} Q) (Q^T \Lambda^{-1} Q) (Q^T \Lambda^{1/2} Q) = Q^T \Lambda^{1/2} \Lambda^{-1} \Lambda^{1/2} Q = Q^TQ = I_n.$$

Also, $\det A = \sqrt{\lambda_1\cdots \lambda_n} = \sqrt{\det \Sigma}$ holds.

Now consider the random variable $Z = A^{-1}(X - \mu)$. From the density transformation theorem we know that the density of $Z$ is given by

$$(\det A) \cdot \frac{1}{(2\pi)^{n/2} \sqrt{\det \Sigma}} \exp(-\tfrac{1}{2} z^t A^t \Sigma^{-1} A z) = (2\pi)^{-n/2} \exp(-\tfrac{1}{2} z^t z).$$

But this means that $Z$ consists of independent standard normal distributions, which we wanted to prove.