How does one compute the convolution of two half-normal distributions?

The half-normal distribution is a special case of the folded normal distribution. If $X$ follows an ordinary normal distribution $N(0, \sigma^{2})$, then $Y = |X|$ follows a half-normal distribution.

The probability density function associated with a half-normal distribution $Y$ is is given by: $$f_{Y} (y, \sigma) = \frac{\sqrt{2}}{\sigma \sqrt{\pi}} e^{- \frac{y^{2}}{2 \sigma^{2}}} \quad .$$

Let $H(0, \sigma^{2})$ denote the half-normal distribution. Let $Y_{1} \sim H(0, \sigma_{1}^{2})$ and $Y_{2} \sim H(0, \sigma_{2}^{2})$. I'd like to figure out what the probability density function is that is associated with the sum of these two random variables: $Z = Y_{1} + Y_{2}$. In other words, I want to compute the convolution $h(z) = f_{Z} (z) = (f_{Y_{1}} * f_{Y_{2}}) (z)$.

Here's what I tried so far:

$\begin{equation} \begin{split} h(z) & = (f_{Y_{1}} * f_{Y_{2}}) (z) \\ &= \int_{- \infty}^{\infty} f(t) g(z-t) dt \\ & = \int_{- \infty}^{\infty} \Big{(} \frac{1}{\sigma_{1}} \sqrt{\frac{2}{\pi}} e^{-\frac{t^{2}}{2 \sigma_{1}^{2}}} \Big{)} \Big{(} \frac{1}{\sigma_{2}} \sqrt{ \frac{2}{\pi} } e^{- \frac{(z-t)^{2}}{2 \sigma_{2}^{2}}} \Big{)} dt \\ &= \int_{- \infty}^{\infty} \frac{1}{\sigma_{1} \sigma_{2}} \cdot \frac{2}{\pi} e^{- \frac{t^{2}}{2 \sigma_{1}^{2}}} e^{- \frac{z^{2}-2t + t^{2}}{2 \sigma_{2}^{2}}} dt \\ &= \int_{- \infty}^{\infty} \frac{1}{\sigma_{1} \sigma_{2}} \cdot \frac{2}{\pi} e^{-\frac{1}{2} \big{(}\frac{\sigma_{1}^{2}z^{2} + (\sigma_{1}^{2} + \sigma_{2}^{2})t^{2} - 2 \sigma \sigma_{1}^{2} t}{\sigma_{1}^{2} \sigma_{2}^{2}} \big{)}} dt \\ &= \frac{1}{\sigma_{1} \sigma_{2}} \cdot \frac{2}{\pi} e^{-\frac{1}{2} \big{(} \frac{\sigma_{1}^{2} z^{2}}{\sigma_{1}^{2} \sigma_{2}^{2}} \big{)}} \int_{\ -\infty}^{\infty} e^{- \frac{1}{2} \big{(} \frac{(\sigma_{1}^{2}+\sigma_{2}^{2})t^{2} - 2 \sigma_{1}^{2}t}{\sigma_{1}^{2} \sigma_{2}^{2}} \big{)}} dt \qquad (1) \\ &= \frac{1}{\sigma_{1} \sigma_{2}} \cdot \frac{2}{\pi} e^{-\frac{1}{2} \big{(} \frac{\sigma_{1}^{2} z^{2}}{\sigma_{1}^{2} \sigma_{2}^{2}} \big{)}} \Bigg{[} \frac{\sqrt{\frac{\pi}{2}} \sigma_{1} \sigma_{2} e^{\frac{\sigma_{1}^{2}}{2\sigma_{2}^{2}(\sigma_{1}^{2}+\sigma_{2}^{2})}}erf\big{(} \frac{\sigma_{1}^{2}(t-1) + \sigma_{2}^{2}t}{\sqrt{2}\sigma_{1} \sigma_{2}} \big{)}} \qquad {\sqrt{\sigma_{1}^{2} + \sigma_{2}^{2}}} \Bigg{]}^{t = \infty}_{t=-\infty} \qquad (2) \\ &= \frac{2 \sqrt{\frac{2}{\pi}}}{\sqrt{\sigma_{1}^{2}+\sigma_{2}^{2}}} e^{- \frac{1}{2} \big{(} \frac{z^{2}}{\sigma_{2}^{2}} - \frac{\sigma_{1}^{2}}{\sigma_{2}^{2} ( \sigma_{1}^{2} + \sigma_{2}^{2})} \big{)}} \\ &= \frac{2 \sqrt{\frac{2}{\pi}}}{\sqrt{\sigma_{1}^{2}+\sigma_{2}^{2}}} e^{- \frac{1}{2} \big{(} \frac{z^{2}(\sigma_{1}^{2} + \sigma_{2}^{2})}{\sigma_{2}^{2}(\sigma_{1}^{2} + \sigma_{2}^{2})} + \frac{\sigma_{1}^{2}}{\sigma_{2}^{2}(\sigma_{1}^{2}+\sigma_{2}^{2})} \big{)}} \qquad (3) \end{split} \end{equation} .$

In $(2)$, $erf(\cdot)$ is the so-called error function. I have to admit I did not calculate the integral in $(1)$ (which yielded equality $(2)$) myself; I did it with Wolframalpha (see the calculation here, with $\sigma_{1}$ and $\sigma_{2}$ changed to $a$ and $b$ respectively for convenience).

After having arrived at equality $(3)$, I don't know how to proceed anymore. If $\frac{\sigma_{1}^{2}}{\sigma_{2}^{2}(\sigma_{1}^{2} + \sigma_{2}^{2})}$ in the exponent would have been multiplied by $z^{2}$, then it seems I would get close to showing that $Z \sim H(0, \sigma_{1}^{2} + \sigma_{2}^{2})$, but the whole calculation does not support the validity of this extra multiplication (and there's still the extra factor $2$ to be dealt with).

I have a number of questions:

1. Did I make a mistake somewhere in my calcuations? Or should I perhaps have done something a bit different in one or more steps?
2. Did Wolframalpha perhaps make a mistake in calculating equation $(1)$ ?

3. Is there perhaps an (other) way to calculate the convolution of two half-normal distributions? Perhaps by means of Laplace/Fourier transforms and the convolution theorem? Can you show it over here directly or can you point me towards an external reference?

The convolution of two half-normal densities is not a half-normal density. Intuitively, this is because there isn't enough "mass" (in fact, there is no mass) on the left-hand side of zero to contribute to the convolution. In particular, the mass of the convolution around zero "slumps". An explicit computation yields a term with error functions that accounts for this "slump" about zero (plot it and see!). $$(f_1 * f_2)(y) = 2\sqrt{\frac{\sigma^2 + \tau^2}{2\pi}} \exp\Big(\frac{-y^2}{2(\sigma^2 + \tau^2)}\Big) \Big( \mathrm{erf}\big(\frac{y\tau}{\sqrt{2\sigma^4 + 2\sigma^2\tau^2}}\big) + \mathrm{erf}\big(\frac{y\sigma}{2\sigma^2\tau^2 + 2 \tau^4}\big)\Big)$$ Indeed, from the note following (23) of On the folded normal distribution:
• I used $$\sigma$$ and $$\tau$$ in place of $$\sigma_1$$ and $$\sigma_2$$ for notational brevity;
• and used $$y$$ as the convolution variable in place of $$z$$).