How to integrate $\exp\left(-ax^2\right)\operatorname{erf}(bx+c)$?

How to integrate $$\int\limits_0^\infty e^{-ax^2} \operatorname{erf}(bx + c) \, dx$$

to get a closed-form solution?

This integration comes from the second integration when I calculate the possibility $P(|A|<|B|),$ where $A \sim N(\mu _A, \sigma _A^2)$ and $B \sim N(\mu _B, \sigma _B^2)$. I have already obtained the exact numerical result, but now I need a closed-form solution.

Thanks in advance.

I am skeptical about a possible closed form.

In the case where $c=0$, there is an explicit result (assuming $a>0$ and $b>0$) $$\int\limits_0^\infty e^{-ax^2} \operatorname{erf}(bx ) \, dx=\frac{\tan ^{-1}\left(\frac{b}{\sqrt{a}}\right)}{\sqrt{\pi a} }$$ Assuming that $c$ is small, $\operatorname{erf}(bx +c)$ could be expanded as series $$\operatorname{erf}(bx +c)=\operatorname{erf}(bx)+\frac{2 e^{-b^2 x^2}}{\sqrt{\pi }}\left(c-b xc^2+\frac{1}{3} \left(2 b^2 x^2-1\right)c^3-\frac{1}{6} b x \left(2 b^2 x^2-3\right)c^4+O\left(c^5\right) \right)$$

This would lead to $$\frac{\tan ^{-1}\left(\frac{b}{\sqrt{a}}\right)}{\sqrt{\pi a}}+k-\frac{b }{\sqrt{\pi }}k^2-\frac{a }{3}k^3+\frac{ b\left(3 a+b^2\right)}{6 \sqrt{\pi }}k^4+O\left(c^5\right) \quad\text{with} \quad k=\frac{c}{\sqrt{a+b^2}}$$ The key problem is to find the general term of the expansion of $\operatorname{erf}(bx +c)$ around $c=0$

In any manner, doing it, we should face integrals $$I_n=\int_0^\infty x^n \,e^{ -\left(a+b^2\right)x^2}\,dx=\frac{1}{2} \Gamma \left(\frac{n+1}{2}\right) \left(a+b^2\right)^{-\frac{n+1}{2}}$$

Edit

We can make things a bit simpler.

Let $bx=y$ and $\alpha=\frac a {b^2}$ to make $$\int\limits_0^\infty e^{-ax^2} \operatorname{erf}(bx+c ) \, dx=\frac 1b\int\limits_0^\infty e^{-\alpha y^2} \operatorname{erf}(y+c ) \, dy$$

$$\operatorname{erf}(y+c )=\operatorname{erf}(y)+\frac{2 e^{-y^2}}{\sqrt{\pi }}c\sum_{n=0}^\infty \frac{2^{\left\lfloor \frac{n+1}{2}\right\rfloor }}{(n+1)!}\,P_{n}(y)\, c^n$$

The polynomials are given below $$\left( \begin{array}{cc} n & P_n(y) \\ 0 & 1 \\ 1 & -y \\ 2 & 2 y^2-1 \\ 3 & 3 y-2 y^3 \\ 4 & 4 y^4-12 y^2+3 \\ 5 & -4 y^5+20 y^3-15 y \\ 6 & 8 y^6-60 y^4+90 y^2-15 \\ 7 & -8 y^7+84 y^5-210 y^3+105 y \\ 8 & 16 y^8-224 y^6+840 y^4-840 y^2+105 \\ 9 & -16 y^9+288 y^7-1512 y^5+2520 y^3-945 y \\ 10 & 32 y^{10}-720 y^8+5040 y^6-12600 y^4+9450 y^2-945 \\11 & -32 y^{11}+880 y^9-7920 y^7+27720 y^5-34650 y^3+10395 y \\ 12 & 64 y^{12}-2112 y^{10}+23760 y^8-110880 y^6+207900 y^4-124740 y^2+10395 \end{array} \right)$$ and $$J_n=\int_0^\infty y^n \,e^{(\alpha +1) \left(-y^2\right)}\,dy=\frac{1}{2} (\alpha +1)^{-\frac{n+1}{2}} \Gamma \left(\frac{n+1}{2}\right)$$

It "just" remains to identify these polynomials.

After integration, this would lead tp $$\frac{\tan ^{-1}\left(\frac{1}{\sqrt{\alpha }}\right)}{\sqrt{\pi } \sqrt{\alpha }}+k-\frac{k^2}{\sqrt{\pi }}-\frac{\alpha k^3}{3}+\frac{(1247400 \alpha +415800) k^4}{2494800 \sqrt{\pi }}+\frac{\alpha ^2 k^5}{10}-$$ $$\frac{\left(415800 \alpha ^2+277200 \alpha +83160\right) k^6}{2494800 \sqrt{\pi }}-\frac{\alpha ^3 k^7}{42}+$$ $$\frac{\left(103950 \alpha ^3+103950 \alpha ^2+62370 \alpha +14850\right) k^8}{2494800 \sqrt{\pi }}+\frac{\alpha ^4 k^9}{216}-$$ $$\frac{\left(20790 \alpha ^4+27720 \alpha ^3+24948 \alpha ^2+11880 \alpha +2310\right) k^{10}}{2494800 \sqrt{\pi }}-\frac{\alpha ^5 k^{11}}{1320}+$$ $$\frac{\left(3465 \alpha ^5+5775 \alpha ^4+6930 \alpha ^3+4950 \alpha ^2+1925 \alpha +315\right) k^{12}}{2494800 \sqrt{\pi }}+O\left(k^{13}\right)$$ where $k=\frac{c}{\sqrt{\alpha +1}}$.