Joint distribution of the sum and product of two i.i.d. centered normal random variables

Question

If $X$ and $Y$ are independent normally distributed random variables

$$X,Y\sim\mathcal{N}(0,\sigma^2)$$

How are the sum and product, $X+Y$ and $XY$, co-distributed? You can write the moment generating function

$$V=X+Y\qquad Z=XY$$ \begin{align} M_{V,Z}(s,t)&=\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }\exp\left[s (v+z)+t (vz)\right]\psi (v,0,\sigma^2) \psi (z,0,\sigma^2)\mathrm{d}v\mathrm{d}z\\ &=\frac{\exp\left[\dfrac{s^2 \sigma ^2}{1-\sigma ^2 t}\right]}{\sqrt{1-\sigma ^4 t^2}} \end{align}

But I do not know what probability density function this belongs to, besides that it looks like the MGF for a modified Bessel function of the second kind and normal distribution. Attached is a histogram generated from $10^7$ such pairs of random numbers.

EDIT: Turns out Mathematica can handle this thing, and spits this monster out $$f(v,z;\sigma^2)=\frac{\sqrt[4]{-1} \left[\sqrt{-i \left(v^2-4 z\right)}-i \sqrt{i \left(v^2-4 z\right)}\right] }{\left| v^2-4 z\right| }\frac{\exp\left(-\dfrac{v^2-2 z}{2 \sigma ^2}\right)}{2 \pi \sigma ^2}$$ Does anyone know if this related to some well-known result I could read around?

Answer 1

The usual method works. Assume without loss of generality that $\sigma^2=1$, then, for every measurable bounded function $u$, $$ E(u(V,Z))=E(u(X+Y,XY)), $$ hence $$ E(u(V,Z))=\iint u(x+y,xy)\,\mathrm e^{-(x^2+y^2)/2}\,\frac{\mathrm dx\mathrm dy}{2\pi}, $$ and the task is to transform the RHS into $$ \iint u(v,z)\,g(v,z)\,\mathrm dv\mathrm dz, $$ since then the function $g$ shall be the density of $(V,Z)$.

The change of variable $(x,y)\to(v,z)=(x+y,xy)$ is two-to-one, it has Jacobian $\mathrm dv\mathrm dz=|x-y|\,\mathrm dx\mathrm dy$, and it is inverted through the formulas $x,y=\frac12(v\pm\sqrt{v^2-4z})$ for every $(v,z)$ such that $4z\leqslant v^2$. Finally, $|x-y|=\sqrt{v^2-4z}$ and $x^2+y^2=v^2-2z$ hence $$ E(u(V,Z))=\iint u(v,z)\,2\,\mathrm e^{-(v^2-2z)/2}\,\mathbf 1_{v^2\gt4z}\,\frac{\mathrm dv\mathrm dz}{2\pi\sqrt{v^2-4z}}, $$ and the density of $(V,Z)$ is $$ g(v,z)=\frac1\pi\,\mathrm e^{-v^2/2}\,\mathrm e^{z}\,\frac{\mathbf 1_{v^2\gt4z}}{\sqrt{v^2-4z}}. $$ If $\sigma^2\ne1$, $(V,Z)$ is transformed into $(\sigma V,\sigma^2Z)$ hence the density becomes $$ g(v,z)=\frac1\pi\,\frac1{\sigma^2}\,\mathrm e^{-v^2/(2\sigma^2)}\,\mathrm e^{z/\sigma^2}\,\frac{\mathbf 1_{v^2\gt4z}}{\sqrt{v^2-4z}}. $$ In other words, in the formula you suggest at the end of your post, the funny looking factor $$ \frac{\sqrt[4]{-1} \left[\sqrt{-i \left(v^2-4 z\right)}-i \sqrt{i \left(v^2-4 z\right)}\right] }{\left| v^2-4 z\right| }$$ should read $$ \frac{2\,\mathbf 1_{v^2\gt4z}}{\sqrt{v^2-4z}}. $$