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]];
];
Let us start with a warm-up exercise. Introduce the functions
$$g_{\pm}(x)=\operatorname{Ai}(x)\pm i\operatorname{Bi}(x).$$
Computing the Wronskian of these two solutions of the Airy equation, one can check that
$$\frac{1}{\operatorname{Ai}^2(x)+\operatorname{Bi}^2(x)}=\frac{\pi}{2i}\left[\frac{g_+'(x)}{g_+(x)}-\frac{g'_-(x)}{g_-(x)}\right]$$
This gives the integral $\mathcal{K}(0)$ as
$$\mathcal{K}(0)=\pi\left[\arg g_+(\infty)-\arg g_+(0)\right]=\pi\left[\pi-\frac{\pi}{3}\right]=\frac{\pi^2}{6}.$$
To compute the integral $\mathcal{K}(3n)$, we will need to develop a more sophisticated approach. First note that (see here)
$$g_{\pm}(x)=-2e^{\mp 2\pi i/3}\operatorname{Ai}\left(e^{\mp2\pi i/3}x\right).$$
Therefore
\begin{align}\mathcal{K}(3n)&=\frac{\pi}{2i}\int_0^{\infty}x^{3n}\left[\frac{g_+'(x)}{g_+(x)}-\frac{g'_-(x)}{g_-(x)}\right]dx=\\
&=\frac{\pi}{2i}\lim_{R\rightarrow\infty}\int_{S_R}z^{3n}\frac{\operatorname{Ai}'(z)}{\operatorname{Ai}(z)}dz,
\end{align}
where the contour $S_R$ in the complex $z$-plane
is composed of two segments: one going from $Re^{2\pi i/3}$ to $ 0$ and another one going from $0$ to $ Re^{-2\pi i/3}$.
It is a well-known fact that the Airy function $\operatorname{Ai}(z)$ has zeros (i.e. our integrand has poles) on the negative real axis only. Therefore by residue theorem our integral is equal to
$$\mathcal{K}(3n)=-\frac{\pi}{2i}\lim_{R\rightarrow \infty}\int_{C_R}z^{3n}\left[\ln\operatorname{Ai}(z)\right]'dz,\tag{1}$$
where $C_R$ is the arc of the circle of radius $R$ centered at the origin going counterclockwise from $Re^{-2\pi i/3}$ to $Re^{2\pi i/3}$.
The limit (1), on the other hand, can be computed using the asymptotics of the Airy function as $z\rightarrow\infty$ for $|\arg z|<\pi$:
\begin{align}
\ln\operatorname{Ai}(z)\sim -\frac23 z^{3/2}-\ln2\sqrt{\pi}-\frac14\ln z+
\ln\sum_{k=0}^{\infty}\frac{(-1)^k\left(\frac16\right)_k\left(\frac56\right)_k}{k!}\left(\frac43 z^{3/2}\right)^{-k}.\tag{2}
\end{align}
Note that if we introduce instead of $z$ the variable $s=\frac43z^{3/2}$, then the integration will be done over the circle of radius $\Lambda=\frac43 R^{3/2}$, i.e. a closed contour in the complex $s$-plane. The corresponding integral can therefore be computed by residues by picking the necessary term in the large $s$ expansion of $\ln \operatorname{Ai}(z)$.
More precisely, we have the following formula:
\begin{align}
\mathcal{K}(3n)=-\frac{\pi}{2i}\lim_{\Lambda\rightarrow \infty}\oint_{|s|=\Lambda}
\left(\frac{3s}{4}\right)^{2n}d\left[-\frac16\ln s+\ln\sum_{k=0}^{\infty}\frac{(-1)^k\left(\frac16\right)_k\left(\frac56\right)_k}{k!}s^{-k}\right]\tag{3}
\end{align}
To compute the residue, it suffices to expand the logarithm-sum up to order $2n$ in $s^{-1}$. Note that the Pochhammer symbol coefficients are in fact some rational numbers.
In the simplest case $n=0$, the residue is determined by the term $-\frac16\ln s$ and we readily reproduce the previous result
$$\mathcal{K}(0)=-\frac{\pi}{2i}\cdot 2\pi i\cdot\left(-\frac16\right)=\frac{\pi^2}{6}.$$
The general formula for arbitrary $n$ would look a bit complicated (but straightforward to obtain) due to the need to expand the logarithm of the sum.
Example. The calculation of the corresponding values $M(n)=\mathcal{K}(3n)$ in Mathematica can be done using the command
\begin{align}\mathtt{\text{ M[n_] := -Pi^2 SeriesCoefficient[
Series[(3 s/4)^(2 n) D[-Log[s]/6 +}} \\ \mathtt{\text{ Log[Sum[(-1)^k
Pochhammer[1/6, k] Pochhammer[5/6, k]/(k! s^k), }} \\
\mathtt{\text{{k, 0, 2 n}]], s], {s, Infinity, 1}], 1]}} \end{align}
This yields, for instance,
$$M(0)=\frac{\pi^2}{6},\quad M(1)=\frac{5\pi^2}{32},\quad M(2)=\frac{565\pi^2}{512},$$ $$\ldots, M(10)=\frac{2\,660\,774\,144\,147\,177\,521\,025\,228\,125\,\pi^2}{2\,199\,023\,255\,552},\ldots$$
and so on.
Added: The large $s$ expansion (2) can also be found directly by using Airy equation. Moreover, by transforming it into an equation for $\ln \operatorname{Ai}(z)$, one can avoid reexpanding the logarithm of the sum. The price to pay will be that the expansion coefficients will be determined by a nonlinear recurrence relation instead of explicit formulas.
Best Answer
It suffices to differentiate, e.g.
$$\begin{gathered}\left(\mathrm{Ai}(x)\mathrm{Hi}'(x)-\mathrm{Ai}'(x)\mathrm{Hi}(x)\right)'= \mathrm{Ai}(x)\mathrm{Hi}''(x)-\mathrm{Ai}''(x)\mathrm{Hi}(x)=\\= \mathrm{Ai}(x)\left(x\mathrm{Hi}(x)+\tfrac1\pi\right)-x\mathrm{Ai}(x)\mathrm{Hi}(x)=\tfrac1\pi\mathrm{Ai}(x)\end{gathered}.$$ We see that the derivatives of the left and right side of the first identity coincide, which implies that the corresponding functions can only differ by a constant. To show that this constant is zero, it suffices to check that the right side vanishes as $x\to\infty$.
P.S. The proof of the second identity is exactly the same. Note that you have copied incorrectly its right side from DLMF.