Given two normally distributed variables $X,Y$ in $\mathcal{N}\left(0,1\right)$, it is stated in Wikipedia that the product distribution of $Z=XY$ is a modified Bessel function of the second kind
$$f_{Z}\left(z\right)=\dfrac{K_{0}\left(\left|z\right|\right)}{\pi}$$
Now given two variables $Z_{1}, Z_{2}$ with this Bessel distribution, I was wondering if there is a close-form expression for the distribution of their sum $Z=Z_{1}+Z_{2}$. This boils down to evaluating the convolution
$$f_{Z}\left(z\right)=\dfrac{1}{\pi^{2}}\displaystyle\int_{-\infty}^{\infty}{\rm d}q\;K_{0}\left(\left|q\right|\right)K_{0}\left(\left|z-q\right|\right)$$
This problem originates in atomic physics in the study of absorption lines. The transverse components of the wave-vector of a finite Gaussian beam are distributed normally. Also the velocities of thermal atoms follow a normal distribution. I am looking for the distribution of the Doppler shift $\boldsymbol{k}\cdot\boldsymbol{v}=k_{x}v_{x}+k_{y}v_{y}$.