Here is my trial, which is partially successful but still not fully answering to your question.
Using some coordinate change, I derived that
$$ I_{n} := \int_{0}^{\infty} \mathrm{erfc}^{n}(x) \, dx = \frac{1}{\sqrt{n}} \left( \frac{2}{\sqrt{\pi}} \right)^{n-1} \int_{T^{n-1}} \int_{0}^{\infty} s^{n-1}e^{-(|x|^{2}-1)s^{2}} \mathrm{erfc}(s) \, ds d\sigma_{x}, \tag{*} $$
where $T^{n-1} = \{ x \in [0, \infty)^{n} : x_{1} + \cdots + x_{n} = \sqrt{n} \}$ is an $(n-1)$-simplex and $d\sigma_{x}$ is the surface measure on $T^{n-1}$.
Case $n = 1$. Since $T^{0} = \{1\}$ and $d\sigma_{x} = \delta(x-1)$ is a unit mass located at $x = 1$, the equation $\text{(*)}$ gives no new information.
Case $n = 2.$ It is easy to find that
$$ \int_{0}^{\infty} s e^{-(|x|^{2}-1)s^{2}} \mathrm{erfc}(s) \, ds = \frac{1}{2|x|(|x|+1)}$$
and by suitable parametrization of the line segment $T^{1}$, we get
\begin{align*}
I_{2}
&= \sqrt{\frac{2}{\pi}} \int_{T^{1}} \frac{1}{2|x|(|x|+1)} \, d\sigma_{x} \\
&= \sqrt{\frac{2}{\pi}} \int_{0}^{\sqrt{2}} \frac{\sqrt{2} \, dt}{2\sqrt{t^{2} + (\sqrt{2}-t)^{2}} ( \sqrt{t^{2} + (\sqrt{2}-t)^{2}} + 1)} \\
&= \sqrt{\frac{2}{\pi}} \int_{0}^{\frac{\pi}{4}} \frac{d\theta}{1+\cos\theta}
= \sqrt{\frac{2}{\pi}} \tan\left(\frac{\pi}{8}\right)
= \frac{2 - \sqrt{2}}{\sqrt{\pi}}.
\end{align*}
We can also write
$$ I_{2} = \boxed{ \displaystyle \frac{2}{\sqrt{\pi}} - \frac{2\sqrt{2}}{\pi^{3/2}} \arctan(\infty) }. $$
Case $n = 3$. We have
$$ \int_{0}^{\infty} s^{2}e^{-(|x|^{2}-1)s^{2}} \mathrm{erfc}(s) \, ds = \frac{1}{2\sqrt{\pi}(|x|^{2}-1)} \left( \frac{\arctan\sqrt{|x|^{2}-1}}{\sqrt{|x|^{2}-1}} - \frac{1}{|x|^{2}} \right). $$
Now we know that $T^{2}$ is a regular triangle with side length $\sqrt{6}$. Thus by introducing polar coordinate change, we get
\begin{align*}
I_{3}
&= \frac{2}{\pi^{3/2}\sqrt{3}} \int_{T^{2}} \frac{1}{|x|^{2}-1} \left( \frac{\arctan\sqrt{|x|^{2}-1}}{\sqrt{|x|^{2}-1}} - \frac{1}{|x|^{2}} \right) \, d\sigma_{x} \\
&= \frac{4\sqrt{3}}{\pi^{3/2}} \int_{0}^{\frac{\pi}{3}} \int_{0}^{\frac{\sec\theta}{\sqrt{2}}} \frac{1}{r} \left( \frac{\arctan r}{r} - \frac{1}{r^{2} + 1} \right) \, drd\theta \\
&= \frac{4\sqrt{3}}{\pi^{3/2}} \int_{0}^{\frac{\pi}{3}} \left( 1- \frac{\arctan\left(\sec\theta / \sqrt{2}\right)}{\sec\theta / \sqrt{2}} \right) \, d\theta \\
&= \frac{4\sqrt{3}}{\pi^{3/2}} \int_{\frac{1}{\sqrt{2}}}^{\sqrt{2}} \left( 1 - \frac{\arctan u}{u} \right) \frac{du}{u\sqrt{2u^{2} - 1}}.
\end{align*}
I haven't tried the last integral, but Mathematica suggests that
$$ \int_{\frac{1}{\sqrt{2}}}^{\sqrt{2}} \left( 1 - \frac{\arctan u}{u} \right) \frac{du}{u\sqrt{2u^{2} - 1}} = \frac{\sqrt{3}\pi}{4} -\sqrt{\frac{3}{2}} \arctan\left( \sqrt{2} \right).$$
This gives
$$ I_{3} = \boxed{ \displaystyle \frac{3}{\sqrt{\pi}} - \frac{6\sqrt{2}}{\pi^{3/2}} \arctan\left( \sqrt{2} \right) }, $$
which can be numerically checked.
I am trying to generalize this calculation for higher dimensions, but it seems somewhat daunting. For example, if $n = 4$ we have to evaluate
$$ \int_{0}^{\infty} \mathrm{erfc}^{4}(x) \, dx = \frac{1}{\pi^{3/2}} \int_{T^{3}} \frac{2|x|+1}{(|x|+1)^{2}|x|^{3}} \, d\sigma_{x}. $$
Nevertheless, using the formula
$$ \int_{T^{n-1}} f(|x|) \, d\sigma_{x} = n! \int_{0}^{\sqrt{\frac{1}{n-1}}} \int_{0}^{\sqrt{\frac{n}{n-2}}s_{1}} \cdots \int_{0}^{\sqrt{\frac{3}{1}}s_{n-2}} f\left( \sqrt{1+s_{1}^{2} + \cdots + s_{n-1}^{2}} \right) \, ds_{n-1} \cdots ds_{1} $$
together with aid of Mathematica, I was able to evaluate $I_{4}$ and it was
$$ I_{4} = \boxed{ \displaystyle \frac{4}{\sqrt{\pi}} - \frac{24\sqrt{2}}{\pi^{3/2}} \arctan \left( \frac{1}{\sqrt{8}} \right) }. $$
These series of observations lead us to believe that $I_{n}$ is of the form
$$ I_{n} = \frac{n}{\sqrt{\pi}} - \frac{n!\sqrt{2}}{\pi^{3/2}} \arctan \alpha_{n} $$
for some reasonably nice $\alpha_{n}$ (with $\alpha_{1} = 0$, $\alpha_{2} = \infty$, $\alpha_{3} = \sqrt{2}$ and $\alpha_{4} = \frac{1}{\sqrt{8}}$), but inverse symbolic calculations show that this seems no longer the case for $n \geq 5$.
Indeed, for $n = 5$ we have
$$ I_{5} = \frac{5}{\sqrt{\pi}} - \frac{80\sqrt{2}}{\pi^{3/2}} \arctan \sqrt{\frac{2}{3}} + \frac{240\sqrt{2}}{\pi^{5/2}} A\left( \sqrt{\frac{5}{3}}, \sqrt{\frac{3}{2}}, \sqrt{\frac{4}{5}} \right), $$
where $A(p, q, r)$ is the generalized Ahmed's integral. I have no idea whether it will reduce to a familiar closed form or not.
Best Answer
Notice that for $a>0$, we have $$F(ax) = e^{-a^{2}x^{2}}\int_{0}^{ax} e^{y^{2}} \mathrm dy = e^{-a^{2}x^{2}} \int_{0}^{a} u e^{x^{2}u^{2}} \, \mathrm du . \tag{1}$$
Then using $(1)$, we get
$$ \begin{align} I &= \int_{0}^{\infty} F(x) F(x \sqrt{2}) \, \frac{e^{-x^{2}}}{x^{2}} \, \mathrm dx \\&= \int_{0}^{\infty} \int_{0}^{\sqrt{2}} \int_{0}^{1} x e^{-x^{2}} e^{x^{2} y^{2}} x e^{-2x^{2}} e^{x^{2}z^{2}} \, \frac{e^{-x^{2}}}{x^{2}} \, \mathrm dy \, \mathrm dz \, \mathrm dx \\ &= \int_{0}^{\sqrt{2}} \int_{0}^{1} \int_{0}^{\infty} e^{-(4-y^{2}-z^{2})x^{2}} \, \mathrm dx \, \mathrm dy \, \mathrm dz \\ &= \frac{\sqrt{\pi}}{2} \int_{0}^{\sqrt{2}} \int_{0}^{1} \frac{1}{\sqrt{4-y^{2}-z^{2}}} \, \mathrm dy \, \mathrm dz \\&= \frac{\sqrt{\pi}}{2} \int_{0}^{\sqrt{2}} \int_{0}^{\arcsin ( \frac{1}{\sqrt{4-z^{2}}})} \, \mathrm d \theta \, \mathrm d z \tag{2} \\ &= \frac{\sqrt{\pi}}{2} \int_{0}^{\sqrt{2}} \arcsin \left( \frac{1}{\sqrt{4-z^{2}}} \right) \, \mathrm dz \\ &= \frac{\sqrt{\pi}}{2} \left( \frac{ \sqrt{2} \pi}{4} - \int_{0}^{\sqrt{2}} \frac{z^{2}}{\sqrt{3-z^{2}} (4-z^{2})} \, \mathrm dz \right) \tag{3} \\ &= \frac{\sqrt{\pi}}{2} \left( \frac{\sqrt{2} \pi}{4} - \int_{1 / \sqrt{2}}^{\infty} \frac{1}{\sqrt{3u^{2}-1} (4u^{2}-1)} \frac{\mathrm du}{u}\right) \tag{4} \\ &= \frac{\sqrt{\pi}}{2} \left( \frac{\sqrt{2} \pi}{4} - 3 \int_{1 /\sqrt{2}}^{\infty} \frac{1}{ (4w^{2}+1)(w^{2}+1)} \, \mathrm dw\right) \tag{5}\\ &=\frac{\sqrt{\pi}}{2} \left( \frac{\sqrt{2} \pi}{2} - 4 \int_{1/ \sqrt{2}}^{\infty} \frac{1}{4w^{2}+1} + \int_{1/ \sqrt{2}}^{\infty} \frac{1}{w^{2}+1} \, \mathrm dw\right) \\ &= \frac{\sqrt{\pi}}{2} \left[ \frac{\sqrt{2} \pi}{4} - \pi +2 \arctan \left( \sqrt{2} \right) +\frac{\pi}{2} - \arctan \left( \frac{1}{\sqrt{2}} \right) \right] \\ &= \frac{\sqrt{\pi}}{2} \left( \frac{\sqrt{2} \pi}{4} - \pi + 3 \arctan{\sqrt{2}} \right). \end{align}$$
$(2)$ Let $y=\sqrt{4-z^{2}}\sin \theta$.
$(3)$ Integrate by parts.
$(4)$ Let $z = \frac{1}{u}$.
$(5)$ Let $w^{2}=3u^2-1$.
EDIT:
Using the same approach, I get
$$ \int_{0}^{\infty} F(ax) F(bx) \, \frac{e^{-p^{2}x^{2}}}{x^{2}} \, \mathrm dx $$
$$ = \frac{\sqrt{\pi}}{2} \left[b \arcsin \left( \frac{a}{\sqrt{a^{2}+p^{2}}} \right) - \sqrt{a^{2}+b^{2}+p^{2}} \arctan \left(\frac{ab}{p \sqrt{a^{2}+b^{2}+p^{2}}} \right) + a \arctan \left( \frac{b}{p} \right)\right]$$
where $a, b,$ and $p$ are all positive parameters.