The multivariate normal distribution of a vector $x\in \mathbf{R}^n$ is defined as $$f(x) = \frac{1}{(2\pi)^n |\Sigma|} e^{-1/2 \cdot (x-\mu)^T\Sigma^{-1}(x-\mu)}$$
where $\Sigma\in \mathbf{R}^{n\times n}$ is the covariance matrix where $\Sigma_{ij} = \mathbf{Cov}(x_i,x_j)$ (note the special case when $i=j$ is just $\mathbf{var}(x_i)$), $|\Sigma|$ is the determinant of the matrix, and $\mu\in\mathbf{R}^n$ is the mean.
With all this out of the way, the answer to your second question about how to represent the distribution in $XY$ rather than $UV$ is already answered by the definition of the normal distribution. Specifically,
$$\mu = [\mu_x, \mu_y]^T$$
$$\Sigma = \begin{bmatrix}\sigma_x^2 & \sigma_{xy}^2\\ \sigma_{yx}^2 & \sigma_y^2\end{bmatrix}$$
$$|\Sigma| = \sigma_x^2\sigma_y^2 - \sigma_{xy}^2 \cdot \sigma_{yx}^2$$
Plug in $x = [X,Y]^T$ and that is distribution.
For your first question, notice that we can relate $UV$ and $XY$ by a linear (technically affine) transformation
$$\begin{bmatrix}
U\\
V
\end{bmatrix} =
\begin{bmatrix}
1/\sigma_x & 0\\
0 & 1/\sigma_y
\end{bmatrix}
\begin{bmatrix}
X\\
Y
\end{bmatrix}
+ \begin{bmatrix}
-\mu_x /\sigma_x\\
-\mu_y / \sigma_y
\end{bmatrix}
$$
More succinctly, we have $u = Ax + b$. It just so happens that a linear combination (plus a possible constant) of Gaussian random variables, is in fact Gaussian (this is not obvious). And since we know the Gaussian PDF is defined by the two parameters $\Sigma$ and $\mu$, if we can find those, then we are done.
It should be clear from linearity of expectation that the mean of $u\in \mathbf{R}^2$ is just $A\mu +b$. For the covariance we have the identity $\mathbf{cov}(Az + b) = A\mathbf{cov}(z)A^T$ (this also more or less falls out from linearity of expectation and you can find numerous proofs online)
Doing everything out we find that the distribution of $UV$ is Gaussian with mean vector
$$\mu_{UV} =
\begin{bmatrix}
1/\sigma_x & 0\\
0 & 1/\sigma_y
\end{bmatrix}
\begin{bmatrix}
\mu_x\\
\mu_y
\end{bmatrix} +
\begin{bmatrix}
-\mu_x /\sigma_x\\
-\mu_y / \sigma_y
\end{bmatrix} =
\begin{bmatrix}
0\\
0
\end{bmatrix}
$$
$$\Sigma_{UV} =
\begin{bmatrix}
1/\sigma_x & 0\\
0 & 1/\sigma_y
\end{bmatrix}
\begin{bmatrix}\sigma_x^2 & \sigma_{xy}^2\\
\sigma_{yx}^2 & \sigma_y^2
\end{bmatrix}
\begin{bmatrix}
1/\sigma_x & 0\\
0 & 1/\sigma_y
\end{bmatrix} =
\begin{bmatrix}
1 & \rho \\
\rho & 1
\end{bmatrix}
$$
You can go ahead and compute the inverse and determinant and see that you get the joint PDF you have described