I found a nice and useful approximation of squared error function $$ \mathrm{erf}^{2}\!\left(x\right)=1-\exp\!\left(-\frac{\pi^{2}}{8}x^{2}\right)+\varepsilon\!\left(x\right). $$
I checked numerically that maximum error is bounded by $\left|\varepsilon\!\left(x\right)\right| < 61\cdot10^{-4}$ but I was asked if this could be somehow proved analytically. Or at least if order of the error could be determined in such manner.
Error function is defined as $$ \mathrm{erf}\left(x\right)=\frac{2}{\sqrt{\pi}}\int_{0}^{x}\exp\left(-t^{2}\right)\,\mathrm{d}t = 2\Phi\!\left(x\sqrt{2}\right)-1, $$ where $\Phi\!\left(x\right)$ is normal cumulative distribution function.
And more generally: is numerical check not enough in such cases?
Such kind of inequalities can be proved through Fubini's theorem, since:
$$\operatorname{erf}(x)^2 = \frac{4}{\pi}\iint_{[0,x]^2}e^{-(u^2+v^2)}\,du\,dv=\frac{1}{\pi}\iint_{[-x,x]^2}e^{-(u^2+v^2)}\,du\,dv $$ and the last integral, over a square, can be effectively approximated with an integral over a suitable circle, whose explicit form is given by: $$ 1-\exp\left({-\rho^2 x^2}\right). $$ For instance, it is trivial that: $$\operatorname{erf}(x)^2 \geq \frac{1}{\pi}\iint_{u^2+v^2\leq x^2}e^{-(u^2+v^2)}\,du\,dv = 1-e^{-x^2}. $$ as well as: $$\operatorname{erf}(x)^2 \leq \frac{1}{\pi}\iint_{u^2+v^2\leq 2x^2}e^{-(u^2+v^2)}\,du\,dv = 1-e^{-2x^2}. $$ A tight approximation is given by considering the circle that has the same area of the original square of integration, $[-x,x]^2$:
One may even try to find "the best" constant $\rho$ such that: $$\operatorname{erf}(x)^2 \approx 1-e^{-\rho^2 x^2}$$ by solving a least-square minimization problem. Numerically, such optimal value is around: $$ \rho = 1.1131 $$ that is extremely close to $\rho=\frac{\pi}{\sqrt{8}}=1.11072\ldots$ given by your approximation.