Let us fix real numbers \(s\) and \(t\). The joint cumulative distribution function \(F_{U,V}(s,t)\) satisfies
\begin{eqnarray*}
F_{U,V}(s,t)&=&\frac1{2\pi}\iint_D e^{-\frac{x^2+y^2}2}\,dxdy,
\end{eqnarray*}
where \(D\) is the the region in \(\mathbb R^2\) that satisfies
\begin{eqnarray*}
D=\left\{(x,y): \alpha x+\beta y \leq s,\; \beta x+\gamma y\leq t\right\}.
\end{eqnarray*}
Let us introduce the substitution \(u=\alpha x+ \beta y\) and \(v=\beta x+\gamma y\). Then the bounds for integration are \begin{eqnarray*}
u&\in&(-\infty,s]\\
v&\in&(-\infty,t].\end{eqnarray*}
Let us denote \[A=\left[\begin{array}{cc}\alpha&\beta\\ \beta &\gamma\end{array}\right].\]
Then the relation between pairs \((x,y)\) and \((u,v)\) can be written as \[\left[\begin{array}{c}u\\ v\end{array}\right]=A\cdot \left[\begin{array}{c}x\\ y\end{array}\right].\]
The solution to the system is \[\left[\begin{array}{c}x\\ y\end{array}\right]=A^{-1}\left[\begin{array}{c}u\\ v\end{array}\right],\] where \(A^{-1}\) is the inverse of the matrix \(A\). The inverse \(A^{-1}\) of the matrix \(A\) is
\begin{eqnarray*}A^{-1}&=&\left[\begin{array}{cc} \frac{\gamma}{\alpha\gamma-\beta^2}& \frac{-\beta}{\alpha\gamma-\beta^2} \\
\frac{-\beta}{\alpha\gamma-\beta^2} & \frac{\alpha}{\alpha\gamma-\beta^2}\end{array}\right].\end{eqnarray*}
The variables \(x\) and \(y\) are explicitly expressed in terms of \(u\) and \(v\) as
\begin{eqnarray*}
x&=& \frac{\gamma}{\alpha\gamma-\beta^2}u- \frac{\beta}{\alpha\gamma-\beta^2}v\\
y&=&- \frac{\beta}{\alpha\gamma-\beta^2}u+ \frac{\alpha}{\alpha\gamma-\beta^2}v.
\end{eqnarray*}
The Jacobian is the absolute value of the determinant of the matrix \(\left[\begin{array}{cc} \frac{\gamma}{\alpha\gamma-\beta^2}& \frac{-\beta}{\alpha\gamma-\beta^2} \\
\frac{-\beta}{\alpha\gamma-\beta^2} & \frac{\alpha}{\alpha\gamma-\beta^2}\end{array}\right]\)
\[J=\left|\mbox{det}\left[\begin{array}{cc} \frac{\gamma}{\alpha\gamma-\beta^2}& \frac{-\beta}{\alpha\gamma-\beta^2} \\
\frac{-\beta}{\alpha\gamma-\beta^2} & \frac{\alpha}{\alpha\gamma-\beta^2}\end{array}\right]\right|= \left|\frac{\alpha\gamma-\beta^2}{(\alpha\gamma-\beta^2)^2}\right|=\frac1{|\alpha\gamma-\beta^2|}.\]
The integral becomes
\begin{eqnarray*}
F_{U,V}(s,t)&=&\frac1{2\pi}\cdot \int_{-\infty}^{t}\int_{\infty}^se^{-\frac12\left(\left( \frac{\gamma}{\alpha\gamma-\beta^2}u- \frac{\beta}{\alpha\gamma-\beta^2}v\right)^2+\left(- \frac{\beta}{\alpha\gamma-\beta^2}u+ \frac{\alpha}{\alpha\gamma-\beta^2}v\right)^2\right)}
\cdot\frac1{\alpha\gamma-\beta^2}\,dudv\\
&=&
\frac1{2\pi\cdot\left(\alpha\gamma-\beta^2\right)}\cdot \int_{-\infty}^{t}\int_{\infty}^s
e^{-\frac{\left( \left(\alpha^2+\gamma^2\right)u^2-2\beta(\alpha+\gamma)uv+\left(\alpha^2+\beta^2\right)v^2
\right)}{2\left(\alpha\gamma-\beta^2\right)^2}}
\,dudv.
\end{eqnarray*}
Therefore the density function is
\begin{eqnarray*}f_{U,V}(s,t)&=&\frac1{2\pi\cdot\left(\alpha\gamma-\beta^2\right)}\cdot
e^{-\frac{\left( \left(\alpha^2+\gamma^2\right)s^2-2\beta(\alpha+\gamma)st+\left(\alpha^2+\beta^2\right)t^2
\right)}{2\left(\alpha\gamma-\beta^2\right)^2}}.
\end{eqnarray*}
Observe that since \(\Sigma=A^2\) we have that \(\mbox{det}(\Sigma)=\mbox{det}(A)^2=\left(\alpha\gamma-\beta^2\right)^2\). Therefore
\(\alpha\gamma-\beta^2=\sqrt{\mbox{det}(\Sigma)}\).
It remains to prove that
\begin{eqnarray*}\frac{\left(\alpha^2+\gamma^2\right)s^2-2\beta(\alpha+\gamma)st+\left(\alpha^2+\beta^2\right)t^2}{\left(\alpha\gamma-\beta^2\right)^2}&=&\left[\begin{array}{cc} s &t \end{array}\right] \cdot \Sigma^{-1}\cdot
\left[\begin{array}{c}s \\ t\end{array}\right].
\end{eqnarray*}
This follows by simple algebraic calculations from \begin{eqnarray*}\Sigma^{-1}&=&\left(A^2\right)^{-1}=\left(A^{-1}\right)^2
\end{eqnarray*} and from the representation \(A^{-1}=\left[\begin{array}{cc} \frac{\gamma}{\alpha\gamma-\beta^2}& \frac{-\beta}{\alpha\gamma-\beta^2} \\
\frac{-\beta}{\alpha\gamma-\beta^2} & \frac{\alpha}{\alpha\gamma-\beta^2}\end{array}\right]\).