As noted by @Zacky the first result can be adapted from this answer from @DrZafarAhmedDSc. Using the series expansion
\begin{equation}
\operatorname{erf}(z)= e^{-z^2} \sum_{n=0}^{\infty} \frac{z^{2n+1}}{\Gamma(n+3/2)}.
\end{equation}
we have
\begin{align}
\mathcal{I}_1&=\int_{0}^{\frac{\pi}{2}}\operatorname{erf}(\sqrt{a}\cos(x))\operatorname{erf}(\sqrt{a}\sin(x))\sin(2x)dx\\
&=2e^{-a}\int_{0}^{\frac{\pi}{2}}\sum_{n=0}^\infty\sum_{m=0}^\infty a^{n+m+1}\frac{\cos^{2n+2}x}{\Gamma(n+3/2)}\frac{\sin^{2m+2}x}{\Gamma(m+3/2)}\,dx\\
&=e^{-a}\sum_{n=0}^\infty\sum_{m=0}^\infty a^{n+m+1}\frac{B(n+3/2,m+3/2)}{\Gamma(n+3/2)\Gamma(m+3/2)}\\
&=e^{-a}\sum_{n=0}^\infty\sum_{m=0}^\infty\frac{a^{n+m+1}}{(n+m+2)!}
\end{align}
where the Beta integral was used. Now the double summation can be rearranged as
\begin{align}
\mathcal{I}_1&=a^{-1}e^{-a}\sum_{m=0}^\infty\sum_{p=m}^\infty\frac{a^{p+2}}{(p+2)!}\\
&=a^{-1}e^{-a}\sum_{p=0}^\infty\sum_{m=0}^{p}\frac{a^{p+2}}{(p+2)!}\\
&=a^{-1}e^{-a}\sum_{p=0}^\infty\frac{(p+1)a^{p+1}}{(p+2)!}\\
&=a^{-1}e^{-a}\left( a\sum_{p=0}^\infty\frac{a^{p+2}}{(p+1)!} -\sum_{p=0}^\infty\frac{a^{p+2}}{(p+2)!}\right)\\
&=a^{-1}e^{-a}\left[a\left( e^a-1 \right)-\left( e^a-1-a \right)\right]\\
&=\frac{1}{a}\left( e^{-a}-1+a\right)
\end{align}
as proposed.
For the second formula, we use the integral representation of $\operatorname{erf}^2$ here:
\begin{equation}
\int_{0}^{1}\frac{e^{-\alpha t^{2}}}{t^{2}+1}\mathrm{d}t=\frac{\pi}{4}e^{a}\left(1-(%
\operatorname{erf}\sqrt{\alpha})^{2}\right)
\end{equation}
(This expression can be derived from the integral definition of $\operatorname{erf}$, by interpreting the product of integrals as a double integral in the unit square and by expressing it in polar coordinates). Then,
\begin{equation}
\operatorname{erf}^2(\sqrt{a}\cos(x))=1-\frac{4}{\pi}\int_0^1e^{-a(1+t^2)\cos^2x}\frac{dt}{1+t^2}
\end{equation}
which can be plug into the integral
\begin{align}
\mathcal{I}_2&=\int_{0}^{\frac{\pi}{2}}\operatorname{erf}^2(\sqrt{a}\cos(x))\cos^2(x)\,dx\\
&=\int_{0}^{\frac{\pi}{2}}\left[1-\frac{4}{\pi}\int_0^1e^{-a(1+t^2)\cos^2x}\frac{dt}{1+t^2}\right]\cos^2x\,dx\\
&=\frac{\pi}{4}-\frac{2}{\pi}\int_0^1\frac{e^{-\frac{a}{2}(1+t^2)}dt}{1+t^2}\int_{0}^{\frac{\pi}{2}}e^{-\frac{a}{2}(1+t^2)\cos2x }\left( 1+\cos 2x \right)\,dx\\
&=\frac{\pi}{4}-\frac{1}{\pi}\int_0^1\frac{e^{-\frac{a}{2}(1+t^2)}dt}{1+t^2}\int_{0}^{\pi}e^{\frac{a}{2}(1+t^2)\cos y }\left( 1-\cos y \right)\,dy\\
&=\frac{\pi}{4}-\int_0^1\frac{e^{-\frac{a}{2}(1+t^2)}}{1+t^2}\left[I_0\left(\frac{a}{2}(1+t^2) \right)-I_1\left(\frac{a}{2}(1+t^2) \right)\right]\,dt
\end{align}
where the classical integral representation for the modified Bessel functions (DLMF) was used. This result corresponds to te proposed formula.
