The integral $\int_{0}^{\pi/2} \sin 2\theta ~ \mbox{erf}(\sin \theta)~ \mbox{erf}(\cos \theta)~d\theta=e^{-1}$

The antiderivative can be found by hand. Integrate by parts $$\newcommand{\erf}{\operatorname{erf}} \mathrm{d}u=\sin 2\theta\,\mathrm{d}\theta,\quad v=\erf(\cos \theta)\erf(\sin \theta) $$ gives $$ -\frac12\cos 2\theta\erf(\sin \theta)\erf(\cos \theta)+ \frac1{\sqrt\pi}\int\cos(2\theta)\left[e^{-\sin^2 \theta}\cos \theta\erf(\cos \theta)-e^{-\cos^2 \theta}\sin \theta\erf(\sin \theta)\right]\,\mathrm{d}\theta $$ Hence we are reduced to computing $$ \int\cos(2\theta)\left[e^{-\sin^2 \theta}\cos \theta\erf(\cos \theta)-e^{-\cos^2 \theta}\sin \theta\erf(\sin \theta)\right]\,\mathrm{d}\theta\\ =e^{-1}\int\cos(2\theta)\left[e^{\cos^2 \theta}\cos \theta\erf(\cos \theta)-e^{\sin^2 \theta}\sin \theta\erf(\sin \theta)\right]\,\mathrm{d}\theta $$ But note that \begin{align*} &\frac{\mathrm{d}}{\mathrm{d}\theta}\left[e^{\cos^2 \theta}\sin \theta\erf(\cos \theta)\right]\\ &=(-2\sin \theta\cos \theta)e^{\cos^2 \theta}\sin \theta\erf(\cos \theta) +e^{\cos^2 \theta}\cos \theta\erf(\cos \theta) +e^{\cos^2 \theta}\sin \theta\frac2{\sqrt\pi}e^{-\cos^2 \theta}(-\sin \theta)\\ &=(-2\sin^2 \theta+1)e^{\cos^2 \theta}\cos \theta\erf(\cos \theta) -\frac{2\sin^2 \theta}{\sqrt\pi}\\ &=\cos (2 \theta)e^{\cos^2 \theta}\cos \theta\erf(\cos \theta) -\frac{2\sin^2 \theta}{\sqrt\pi} \end{align*} and the similar equation with $\cos \theta$ and $\sin \theta$ interchanged. So putting these two together, \begin{align*} &\int\cos(2\theta)\left[e^{\cos^2 \theta}\cos \theta\erf(\cos \theta)-e^{\sin^2 \theta}\sin \theta\erf(\sin \theta)\right]\,\mathrm{d}\theta\\ &=\left( e^{\cos^2 \theta} \erf(\cos \theta) \sin \theta - e^{\sin^2 \theta} \erf(\sin \theta) \cos \theta\right)+\int\frac2{\sqrt\pi}(\sin^2 \theta+\cos^2 \theta)\,\mathrm{d}\theta \end{align*} and of course you know how to do the final integral.

The final upshot is \begin{align*} &\int\sin(2\theta)\erf(\sin\theta)\erf(\cos\theta)\,\mathrm{d}\theta\\ &=\frac2{e\pi}\theta-\frac12\cos 2\theta\erf(\sin \theta)\erf(\cos \theta)\\ &\quad+\frac1{e\sqrt\pi} \left( e^{\cos^2 \theta} \erf(\cos \theta) \sin \theta - e^{\sin^2 \theta} \erf(\sin \theta) \cos \theta\right) \end{align*}

---

For $I$, this approach doesn't work, and I don't think there is a closd form solution to the definite integral. However, we can express the MacLaurin series $$ \tanh z=\sum_{k\geq 0} 2\frac{(-1)^k}{\pi^{2k+2}}\left(1-\frac1{4^{k+1}}\right)\zeta(2k+2)z^{2k+1} $$ so $$ I=4\sum_{k,k'\geq 0}(1-4^{-k-1})(1-4^{-k'-1})\frac{(-1)^{k+k'}\zeta(2k+2)\zeta(2k'+2)\Gamma(k+\frac12)\Gamma(k'+\frac12)}{\Gamma(k+k'+1)\pi^{2(k+k')+4}} $$ which doesn't really simplify the task.

$$J=\int_{0}^{\pi/2} \sin 2\theta \mbox{erf}(\sin \theta) \mbox{erf}(\cos \theta). $$ I get clue from the fact that though $\tanh x$ and erf$(x)$ are similar, but the latter one enjoys beautiful expansions! Expanding $$\mbox{erf}(x)= e^{-x^2} \sum_{n=0}^{\infty} \frac{x^{2n+1}}{\Gamma(n+3/2)}.$$ We re-write $$J=2e^{-1}\sum_{m=0}^{\infty} \sum _{n=0}^{\infty} \int_{0}^{\pi/2} \frac{ \sin^{2m+2} \theta ~\ cos^{2n+2} \theta} {\Gamma(m+3/2)~\Gamma(n+3/2)} ~d \theta.$$
By Beta-integral, we get $$J=e^{-1}\sum_{m=0}^{\infty} \sum_{n=0}^{\infty} \frac{1}{\Gamma(m+n+2)}=e^{-1}\sum_{k=1}^{\infty} \sum_{m=0}^{k-1} \frac{1}{(m+1)!}=e^{-1}\sum_{k=1}^{\infty} \frac{k}{(k+1)!}.$$ $$\Rightarrow e^{-1} \sum_{k=0} \left( \frac{1}{k!}-\frac{1}{(k+1)!} \right)=e^{-1}.$$