Too Long For Comment:
Define the function $\mathcal{I}:\mathbb{N}\times\mathbb{R}_{>0}\rightarrow\mathbb{R}$ via the (convergent) improper integral
$$\mathcal{I}{\left(n,a\right)}:=\int_{0}^{\infty}\mathrm{d}x\,\exp{\left(-a^{2}x^{2}\right)}\left[\operatorname{erf}{\left(x\right)}\right]^{n},$$
where the error function is typically defined via the integral representation
$$\operatorname{erf}{\left(x\right)}:=\frac{2}{\sqrt{\pi}}\int_{0}^{x}\mathrm{d}t\,\exp{\left(-t^{2}\right)};~~~\small{x\in\mathbb{R}},$$
or equivalently
$$\operatorname{erf}{\left(x\right)}=\frac{2x}{\sqrt{\pi}}\int_{0}^{1}\mathrm{d}t\,\exp{\left(-x^{2}t^{2}\right)};~~~\small{x\in\mathbb{R}}.$$
Then by induction, powers of the error function can then be expressed as multiple integrals in the following way:
$$\left[\operatorname{erf}{\left(x\right)}\right]^{n}=\frac{2^{n}x^{n}}{\pi^{n/2}}\int_{[0,1]^{n}}\mathrm{d}t_{1}\dots\mathrm{d}t_{n}\,\prod_{k=1}^{n}\exp{\left(-x^{2}t_{k}^{2}\right)};~~~\small{n\in\mathbb{N}\land x\in\mathbb{R}}.$$
Recalling the definition of the gamma function as
$$\Gamma{\left(z\right)}:=\int_{0}^{\infty}\mathrm{d}t\,t^{z-1}\exp{\left(-t\right)};~~~\small{z>0},$$
we can show that for any $n\in\mathbb{N}\land x\in\mathbb{R}$,
$$\begin{align}
\mathcal{I}{\left(n,a\right)}
&=\int_{0}^{\infty}\mathrm{d}x\,\exp{\left(-a^{2}x^{2}\right)}\left[\operatorname{erf}{\left(x\right)}\right]^{n}\\
&=\int_{0}^{\infty}\mathrm{d}x\,\exp{\left(-a^{2}x^{2}\right)}\frac{2^{n}x^{n}}{\pi^{n/2}}\int_{[0,1]^{n}}\mathrm{d}t_{1}\dots\mathrm{d}t_{n}\,\prod_{k=1}^{n}\exp{\left(-x^{2}t_{k}^{2}\right)}\\
&=\frac{2^{n}}{\pi^{n/2}}\int_{0}^{\infty}\mathrm{d}x\int_{[0,1]^{n}}\mathrm{d}t_{1}\dots\mathrm{d}t_{n}\,x^{n}\exp{\left(-a^{2}x^{2}\right)}\prod_{k=1}^{n}\exp{\left(-x^{2}t_{k}^{2}\right)}\\
&=\frac{2^{n}}{\pi^{n/2}}\int_{[0,1]^{n}}\mathrm{d}t_{1}\dots\mathrm{d}t_{n}\int_{0}^{\infty}\mathrm{d}x\,x^{n}\exp{\left(-a^{2}x^{2}\right)}\prod_{k=1}^{n}\exp{\left(-x^{2}t_{k}^{2}\right)}\\
&=\frac{2^{n}}{\pi^{n/2}}\int_{[0,1]^{n}}\mathrm{d}t_{1}\dots\mathrm{d}t_{n}\int_{0}^{\infty}\mathrm{d}x\,x^{n}e^{-\left(a^{2}+\sum_{k=1}^{n}t_{k}^{2}\right)x^{2}}\\
&=\frac{2^{n-1}}{\pi^{n/2}}\int_{[0,1]^{n}}\mathrm{d}t_{1}\dots\mathrm{d}t_{n}\int_{0}^{\infty}\mathrm{d}u\,u^{(n-1)/2}e^{-\left(a^{2}+\sum_{k=1}^{n}t_{k}^{2}\right)u};~~~\small{\left[x=\sqrt{u}\right]}\\
&=\frac{2^{n-1}}{\pi^{n/2}}\int_{[0,1]^{n}}\mathrm{d}t_{1}\dots\mathrm{d}t_{n}\,\frac{1}{\left(a^{2}+\sum_{k=1}^{n}t_{k}^{2}\right)^{(n+1)/2}}\int_{0}^{\infty}\mathrm{d}t\,t^{(n-1)/2}e^{-t};~~~\small{\left[u=\frac{t}{a^{2}+\sum_{k=1}^{n}t_{k}^{2}}\right]}\\
&=\frac{2^{n-1}}{\pi^{n/2}}\int_{[0,1]^{n}}\mathrm{d}t_{1}\dots\mathrm{d}t_{n}\,\frac{1}{\left(a^{2}+\sum_{k=1}^{n}t_{k}^{2}\right)^{(n+1)/2}}\,\Gamma{\left(\frac{n+1}{2}\right)}\\
&=\frac{2^{n-1}\Gamma{\left(\frac{n+1}{2}\right)}}{\pi^{n/2}}\int_{[0,1]^{n}}\mathrm{d}t_{1}\dots\mathrm{d}t_{n}\,\frac{1}{\left(a^{2}+\sum_{k=1}^{n}t_{k}^{2}\right)^{(n+1)/2}}.\\
\end{align}$$