Integral with exp and erf

Let's try a close variant : differentiate relatively to $d$ instead of $a$ at the start (to avoid the additional factor $(x-d)$ in the integral and complications in the final integration) :

$$I\left(d\right)=\int_{-\infty}^{\infty}e^{-b^{2}\left(x-c\right)^{2}}\ \mathrm{erf}\left(a\left(x-d\right)\right)\,\mathrm{d}x$$

$$\frac{\mathrm{d}I\left(d\right)}{\mathrm{d}d}=\frac {-2a}{\sqrt{\pi}}\int_{-\infty}^{\infty}e^{-b^{2}(x-c)^{2}-a^{2}(x-d)^{2}}\,\mathrm{d}x$$

$$\frac{\mathrm{d}I\left(d\right)}{\mathrm{d}d}=\frac {-2a}{\sqrt{\pi}}\int_{-\infty}^{\infty}e^{-b^{2}(x-c)^{2}-a^{2}(x-d)^{2}}\,\mathrm{d}x$$

$$\frac{\mathrm{d}I\left(d\right)}{\mathrm{d}d}=\frac {-2a}{\sqrt{\pi}}\int_{-\infty}^{\infty}e^{-\left(a^2+b^2\right)\left(x-\frac{b^2c+a^2d}{a^2+b^2}\right)^2+\frac{\left(b^2c+a^2d\right)^2}{a^2+b^2}-\left(b^2c^2+a^2d^2\right)}\,\mathrm{d}x$$

$$\frac{\mathrm{d}I\left(d\right)}{\mathrm{d}d}=\frac {-2a}{\sqrt{\pi}}e^{\frac{\left(b^2c+a^2d\right)^2}{a^2+b^2}-\left(b^2c^2+a^2d^2\right)}\int_{-\infty}^{\infty}e^{-\left(a^2+b^2\right)y^2}\,\mathrm{d}y$$

$$\frac{\mathrm{d}I\left(d\right)}{\mathrm{d}d}=\frac {-2a}{\sqrt{\pi}}e^{\frac{-a^2b^2\left(c-d\right)^2}{a^2+b^2}}\frac{\sqrt{\pi}}{\sqrt{a^2+b^2}}=-2a\frac{e^{\frac{-a^2b^2\left(c-d\right)^2}{a^2+b^2}}}{\sqrt{a^2+b^2}}$$

At this point we have to integrate again relatively to $d$ to get (up to a function $C$ independent of $d$) :

$$I\left(d\right)=\frac {\sqrt{\pi}}b\operatorname{erf}\left(\frac{ab(c-d)}{\sqrt{a^2+b^2}}\right)+C(a,b,c)$$

(Alpha integration check : note that the denominator is $b$ and not $\sqrt{b}$ nor my earlier $ab$ as I checked numerically!)

After that you'll just have to prove that $C(a,b,c)\equiv 0$

Note that the integration relatively to $d$ seems more straightforward than relatively to $a$ in your case (I'm not saying it can't be done your way!).

Hoping it clarified things a little even if it didn't answer your question,