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,
Let $c \in {\mathbb R}$ and $d \in {\mathbb R}$. Without loss of generality we can consider a slightly different function:
\begin{equation}
f(c,d):= \int\limits_{-\infty}^\infty \exp(-t^2) \mbox{erfc}(t-c) \mbox{erfc}(t-d) dt
\end{equation}
Now we differentiate with respect to $c$. We have:
\begin{eqnarray}
\partial_c f(c,d) &=& \lim\limits_{\epsilon \rightarrow 0} \frac{\sqrt{2}}{\sqrt{\pi}} \exp(-\frac{c^2}{2})
\left( \sqrt{\pi} + 2 \sqrt{\pi} T(\epsilon,\frac{1}{\sqrt{2}},\frac{2d-c}{\sqrt{2}}) - 2 \sqrt{\pi} T(\epsilon,\frac{1}{\sqrt{2}},\frac{-2d+c}{\sqrt{2}})\right)\\
&=&\sqrt{2} \exp(-\frac{c^2}{2}) \left( 1-erf(\frac{c-2 d}{\sqrt{6}})\right)
\end{eqnarray}
where $T(h,a,b)$ is the generalized Owen's T function Generalized Owen's T function .
Now, all we need to do is to integrate. Since $f(-\infty,d)=0$ we integrate over $c$ from minus infinity to $c$. We have:
\begin{eqnarray}
f(c,d) &=& \sqrt{2} \left( \sqrt{\frac{\pi}{2}}(1+erf(\frac{c}{\sqrt{2}}) + 2erf(\frac{d}{\sqrt{2}})) + \sqrt{2 \pi} 2 T(c,\frac{1}{\sqrt{3}},-\frac{2 d}{\sqrt{3}})\right) \\
&=& \frac{1}{3 \sqrt{\pi}} \left( 12 \pi T\left(c,\frac{c-2 d}{\sqrt{3} c}\right)+12 \pi T\left(d,\frac{d-2 c}{\sqrt{3} d}\right)-6 \arctan\left(\frac{c-2 d}{\sqrt{3} c}\right)-6 \arctan\left(\frac{d-2 c}{\sqrt{3} d}\right)+3
\pi \text{erf}\left(\frac{c}{\sqrt{2}}\right)+3 \pi \text{erf}\left(\frac{d}{\sqrt{2}}\right)+4 \pi\right)
\end{eqnarray}
where $T(h,a)$ is Owen's T function https://en.wikipedia.org/wiki/Owen%27s_T_function .
In[1060]:= {c, d} = RandomReal[{-2, 2}, 2, WorkingPrecision -> 50];
NIntegrate[Exp[-t^2] Erfc[t - c] Erfc[t - d], {t, -Infinity, Infinity}]
(4 \[Pi] - 6 ArcTan[(c - 2 d)/(Sqrt[3] c)] -
6 ArcTan[(-2 c + d)/(Sqrt[3] d)] + 3 \[Pi] Erf[c/Sqrt[2]] +
3 \[Pi] Erf[d/Sqrt[2]] + 12 \[Pi] OwenT[c, (c - 2 d)/(Sqrt[3] c)] +
12 \[Pi] OwenT[d, (-2 c + d)/(Sqrt[3] d)])/(3 Sqrt[\[Pi]])
Out[1061]= 0.483318
Out[1062]= 0.48331807775609703646923225386370751256977344715
Update: Now let us take four real numbers $a_1 \in {\mathbb R}$, $a_2 \in {\mathbb R}$, $c\in {\mathbb R}$ and $d \in {\mathbb R}$ and consider a more general integral:
\begin{equation}
f^{(a_1,a_2)}(c,d):= \int\limits_{-\infty}^\infty \exp(-t^2) \mbox{erfc}(a_1 t-c) \mbox{erfc}(a_2 t-d) dt
\end{equation}
Then by doing the same calculations as above we easily arrive at the following formula:
\begin{eqnarray}
&&f^{(a_1,a_2)}(c,d)= \frac{1}{\sqrt{\pi}} \left(\right.\\
&&
4 \pi T\left(\frac{\sqrt{2} c}{\sqrt{a_1^2+1}},\frac{a_1 a_2 c-\left(a_1^2+1\right) d}{c \sqrt{a_1^2+a_2^2+1}}\right)+4 \pi T\left(\frac{\sqrt{2}
\sqrt{a_1^2+1} d}{\sqrt{a_1^2 a_2^2+a_1^2+a_2^2+1}},\frac{a_1 a_2 d-\left(a_2^2+1\right) c}{d \sqrt{a_1^2+a_2^2+1}}\right) +\\
&&-2 \arctan\left(\frac{a_1
a_2 c-\left(a_1^2+1\right) d}{c \sqrt{a_1^2+a_2^2+1}}\right)-2 \arctan\left(\frac{a_1 a_2 d-\left(a_2^2+1\right) c}{d \sqrt{a_1^2+a_2^2+1}}\right)+\\
&&\pi
\text{erf}\left(\frac{\sqrt{a_1^2+1} d}{\sqrt{a_1^2 a_2^2+a_1^2+a_2^2+1}}\right)+\pi
\text{erf}\left(\frac{c}{\sqrt{a_1^2+1}}\right)+\\
&&\pi +2 \arctan\left(\frac{a_1 a_2}{\sqrt{a_1^2+a_2^2+1}}\right)\\
&&\left.\right)
\end{eqnarray}
For[count = 1, count <= 200, count++,
{a1, a2, c, d} = RandomReal[{-5, 5}, 4, WorkingPrecision -> 50];
I1 = NIntegrate[
Exp[-t^2] Erfc[a1 t - c] Erfc[a2 t - d], {t, -Infinity, Infinity}];
I2 = 1/Sqrt[\[Pi]] (\[Pi] +
2 ArcTan[(a1 a2)/Sqrt[(1 + a1^2 + a2^2)]] -
2 ArcTan[ (a1 a2 c - (1 + a1^2) d)/(
Sqrt[(1 + a1^2 + a2^2)] c)] -
2 ArcTan[ (a1 a2 d - (1 + a2^2) c)/(
Sqrt[(1 + a1^2 + a2^2)] d)] + \[Pi] Erf[c/Sqrt[
1 + a1^2]] + \[Pi] Erf[(Sqrt[(1 + a1^2)] d)/ Sqrt[
1 + a1^2 + a2^2 + a1^2 a2^2]] +
4 \[Pi] OwenT[(Sqrt[2] c)/Sqrt[
1 + a1^2], (a1 a2 c - (1 + a1^2) d)/(
Sqrt[(1 + a1^2 + a2^2)] c)] +
4 \[Pi] OwenT[(Sqrt[2] Sqrt[(1 + a1^2)] d)/ Sqrt[
1 + a1^2 + a2^2 + a1^2 a2^2], (a1 a2 d - (1 + a2^2) c)/(
Sqrt[(1 + a1^2 + a2^2)] d)]);
If[Abs[I2/I1 - 1] > 10^(-3),
Print["results do not match", {a1, a2, c, d, {I1, I2}}]; Break[]];
If[Mod[count, 10] == 0, PrintTemporary[count]];
];
Best Answer
Consider the related integral $$ G(a)=\int\limits_{x_0}^\infty dx \ e^{-ax^2}\operatorname{erf}(\alpha x+\gamma) $$
So that your integral is $\left(-\frac{d}{da}\right)^m G(1)$. It turns out $G(a)$ evaluates to a sum of trig and OwenT functions, as done here and used here to evaluate a similar integral.