The integral $$\int_{0}^{z} e^{-a^{2} x^{2}} {\rm erf}(bx)\, dx$$ is related to the convolution of two half-normal distributions. This can be inferred from this question on MSE. The following expression for the definite version of this integral is known. See, for instance, "A Table of Integrals of the Error Functions" by Edward W. Ng and Murray Geller, equation 4.3.2 (link here). It is as follows: $$ \int_{0}^{\infty} e^{-a^{2} x^{2} } {\rm erf} (bx) \,dx = \frac{\sqrt{\pi}}{2a} - \frac{1}{a \sqrt{\pi}} \tan^{-1} \big{(} \frac{a}{b} \big{)} .$$

However, I don't know about a similar expression for the indefinite integral. Do you know about a paper or other source that adresses the indefinite integral? Or do you know how to compute it yourself?

The reason I ask about this integral is because I'd like to compute the convolution of two half-normal distributions. Please let me know if you have a reference in which the convolution of two half-normal distributions with unequal variances is calculated.