We have two independent Gaussian random variables with zero mean and variance $\sigma^2$, i.e., $\theta_V \sim \mathcal{N}(0,\sigma^2)$ and $\theta_H \sim \mathcal{N}(0,\sigma^2)$.
Let $X=\theta_V^2+\theta_H^2$ and $Y=2\theta_V^2+2\theta_H^2+\alpha \theta_V$, where $\alpha$ is a real number. How can we derive the joint probability density function (PDF) or cumulative density function (CDF) of $X$ and $Y$?
This answer summarizes my comments: One way is to use the Jacobian matrix method. For a transformation $(V,H)\rightarrow(X,Y)$, if there are $n$ distinct points $(v_i,h_i)$ that map to the same $(x,y)$ value then $$ f_{X,Y}(x,y) = \sum_{i=1}^n f_{V,H}(v_i,h_i)|J(v_i,h_i)|^{-1}$$ where $|J(v_i,h_i)|$ is the determinant of the $2 \times 2$ Jacobian matrix for the point $(v_i,h_i)$: $$ J(v_i,h_i) = \left[\begin{array}{}dx/dv \quad dx/dh\\ dy/dv \quad dy/dh\end{array}\right]\left|_{(v_i,h_i)}\right.$$
In this particular case we have: $$ f_{V,H}(v,h) = \frac{1}{2\pi \sigma^2}e^{-\frac{(v^2+h^2)}{2\sigma^2}} \quad \forall (v,h) \in \mathbb{R}^2$$ Since $X=V^2+H^2$ and $Y=2V^2+2H^2 + \alpha V$, we can compute the Jacobian to get $|J(v,h)|^{-1} = 1/|2\alpha h|$, and for $(x,y)$ values with $n=2$ we get: $$ f_{X,Y}(x,y) = \sum_{i=1}^2 \frac{f_{V,H}(v_i,h_i)}{2|\alpha h_i|}$$ where it suffices to solve for the corresonding $(v_1, h_1)$ and $(v_2, h_2)$ points.