Joint distribution of the sum and product of two i.i.d. centered normal random variables 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?

 A: 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}}.
$$
