Repeated Indefinite Integration of Gaussian Integral

I have an integral that can be solved via recursive integration by parts. In my case, $\mathrm{d}v=e^{-ax^{2}}$. Question: Is there a solution or special function defined as the n-th indefinite integration of the Gaussian function?

\begin{equation}
\mathbf{i}^{n} \exp\{-ax^{2}\}= \mathop{\int\dots\int}\limits_{n \ \text{times}} e^{-ax^{2}} \ \mathrm{d}x^{n} = \ ?
\end{equation}

Anyways, I am aware that there is the solution for this problem http://dlmf.nist.gov/7.18 for the complementary error function:

\begin{equation}
\mathbf{i}^{n}\ \mathrm{erfc}\, z= \int_{z}^{\infty} \mathbf{i}^{n-1} \mathrm{erfc}\, t \ \mathrm{d}t = \frac{2}{\sqrt{\pi}}\int_{z}^{\infty} \frac{(t-z)^{n}}{n!} e^{-t^{2}} \ \mathrm{d}t
\end{equation}

I thought this may be useful. Any help would be great.

Answer

Firstly let us reformulate the problem in a more rigorous way.
Lent $n\ge 1$ be an integer.
We are seeking to solve a following differential equation:
\begin{equation}
\frac{d^n f_n(x)}{d x^n} = \rho(x) := \frac{e^{-x^2/2}}{\sqrt{2 \pi}}
\end{equation}
subject to $f^{(p)}(\infty)=0$ for $p=0,\cdots,n-1$.

Using the definition of the error function https://en.wikipedia.org/wiki/Error_function and integration by parts we easily establish the result for low values of $n$. We have:
\begin{eqnarray}
f_n(x)&=&
\left\{
\begin{array}
-\frac{1}{2} {\rm erfc}\left(\frac{x}{\sqrt{2}}\right) & \mbox{if $n=1$}\\
\rho(x) - \frac{x}{2} {\rm erfc}\left(\frac{x}{\sqrt{2}}\right) & \mbox{if $n=2$}\\
\frac{x}{2} \rho(x) - \frac{1}{4}(1+x^2) {\rm erfc}\left(\frac{x}{\sqrt{2}}\right) & \mbox{if $n=3$}\\
\frac{1}{6}(2+x^2) \rho(x) - \frac{x}{12}(3+x^2) {\rm erfc}\left(\frac{x}{\sqrt{2}}\right) & \mbox{if $n=4$}

\end{array}
\right.
\end{eqnarray}

In[2226]:= rho[x_] := Exp[-x^2/2]/Sqrt[2 Pi];
ll = {-(1/2) Erfc[x/Sqrt[2]], 
   rho[x] - 
    1/2 x Erfc[x/Sqrt[2]], (1/2 x rho[x] - 
     1/4 (1 + x^2) Erfc[x/Sqrt[2]]), (1/6 (2 + x^2) rho[x] - 
     1/12 x (3 + x^2) Erfc[x/Sqrt[2]])};

Table[D[ll[[p]], {x, p}], {p, 1, Length[ll]}] // Simplify


Out[2228]= {E^(-(x^2/2))/Sqrt[2 \[Pi]], E^(-(x^2/2))/Sqrt[
 2 \[Pi]], E^(-(x^2/2))/Sqrt[2 \[Pi]], E^(-(x^2/2))/Sqrt[2 \[Pi]]}

From the above we already see a pattern emerging. Therefore we make a following induction ansatz:
\begin{eqnarray}
f_n(x) = W_1^{(n-1)}(x) \cdot {\rm erfc}\left(\frac{x}{\sqrt{2}}\right) + W_2^{(n-2)}(x) \cdot \rho(x)

\end{eqnarray}
where $W_1^{(n-1)}(x):=\sum\limits_{p=0}^{n-1} a_1^{(n)}(p) x^p$ and $W_2^{(n-2)}(x):=\sum\limits_{p=0}^{n-2} a_2^{(n)}(p) x^p$ are polynomials.

Now we have:
\begin{eqnarray}
&&f_{n+1}(\infty)- f_{n+1}(x) = \\
&&\int\limits_x^\infty \left(W_1^{(n-1)}(\xi) \cdot {\rm erfc}\left(\frac{\xi}{\sqrt{2}}\right) + W_2^{(n-2)}(\xi) \cdot \rho(\xi)\right) d\xi = \\
&&\sum\limits_{p=0}^{n-1} a_1^{(n)}(p) \int\limits_x^\infty \xi^p {\rm erfc}\left(\frac{\xi}{\sqrt{2}}\right) d\xi + \sum\limits_{p=0}^{n-2} a_2^{(n)}(p) \int\limits_x^\infty\xi^p \rho(\xi) d\xi = \\
%%%%
&&\sum\limits_{p=0}^{\lfloor(n-1)/2\rfloor} a_1^{(n)}(2p) \int\limits_x^\infty \xi^{2 p} {\rm erfc}\left(\frac{\xi}{\sqrt{2}}\right) d\xi +\\

&&
\sum\limits_{p=0}^{\lfloor(n-2)/2\rfloor} a_1^{(n)}(2p+1) \int\limits_x^\infty \xi^{2 p+1} {\rm erfc}\left(\frac{\xi}{\sqrt{2}}\right) d\xi +\\
&&\sum\limits_{p=0}^{\lfloor (n-2)/2 \rfloor} a_2^{(n)}(2p) \int\limits_x^\infty\xi^{2 p}\rho(\xi) d\xi + \\
&&
\sum\limits_{p=0}^{\lfloor (n-3)/2 \rfloor} a_2^{(n)}(2p+1) \int\limits_x^\infty\xi^{2 p+1}\rho(\xi) d\xi = \\
&&-\sum\limits_{p=0}^{\lfloor(n-1)/2\rfloor} a_1^{(n)}(2p) \left(\frac{{\rm erfc}\left(\frac{x}{\sqrt{2}}\right) x^{2 p+1}}{2 p+1}-\frac{1}{2p+1} \rho(x)\sum\limits_{q=0}^p (p)_{(q)} 2^{q+1} x^{2p-2q} \right) +\\
&&
-\sum\limits_{p=0}^{\lfloor(n-2)/2\rfloor} a_1^{(n)}(2p+1) \left(\frac{{\rm erfc}\left(\frac{x}{\sqrt{2}}\right) \left(x^{2 p+2}-2^p p! \binom{p+\frac{1}{2}}{p}\right)}{2 (p+1)} - \frac{1}{2p+2} \rho(x) \sum\limits_{q=0}^p (p+1/2)_{(q)} 2^{q+1} x^{2p+1-2 q}\right) +\\
&&-\sum\limits_{p=0}^{\lfloor (n-2)/2 \rfloor} a_2^{(n)}(2p) \left(-2^{p-2} (p-1)! \binom{p-\frac{1}{2}}{p-1} {\rm erfc}\left(\frac{x}{\sqrt{2}}\right)-\rho(x)\sum\limits_{q=0}^{p-1} (p-1/2)_{(q)} 2^q x^{2 p-1-2 q} \right) + \\
&&

-\sum\limits_{p=0}^{\lfloor (n-3)/2 \rfloor} a_2^{(n)}(2p+1) \left(-\rho(x) \sum\limits_{q=0}^p (p)_{(q)} 2^q x^{2 p-2 q} \right)
\end{eqnarray}

Now the only thing that remains is to collect the terms being proportional to the complementary error function and those being proportional to the Gaussian density and then retrieve recurrence relations for the coefficients of our polynomials $W_1^{(n-1)}(x)$ and $W_2^{(n-2)}(x)$. One thing that is already seen from the last line above is that indeed the orders of polynomials are $n$ and $n-1$ for the complementary error function and for the Gaussian density-- as should be.

Update: The polynomials in question read:
\begin{eqnarray}
W_1^{(n-1)}(x)&:=& -\frac{1}{(n-1)!} \sum\limits_{p=0}^{\lfloor (n-1)/2 \rfloor} \binom{n-1}{2 p} (p-1/2)_{(p-1)} 2^{p-2} x^{n-1-2 p}\\
W_2^{(n-2)}(x)&:=& -\frac{1}{(n-1)!} \sum\limits_{p=0}^{\lfloor (n-2)/2 \rfloor} {\mathcal C}^{(n)}_p p! 2^p x^{n-2-2 p}
\end{eqnarray}

where
\begin{eqnarray}
{\mathcal C}^{(n)}_p:=
\sum\limits_{\xi=p+1}^{\lfloor (n-1)/2\rfloor} \binom{n-1}{2 \xi} \binom{\xi-1/2}{p}-
\sum\limits_{\xi=p}^{\lfloor (n-2)/2\rfloor} \binom{n-1}{2\xi+1} \binom{\xi}{p}
\end{eqnarray}

M = 10;
rho[x_] := Exp[-x^2/2]/Sqrt[2 Pi];
CC[n_, p_] := (Sum[

     Binomial[n - 1, 2 xi] Binomial[xi - 1/2, p], {xi, p + 1, 
      Floor[(n - 1)/2]}] - 
    Sum[Binomial[n - 1, 2 xi + 1] Binomial[xi + 0/2, p], {xi, p, 
      Floor[(n - 2)/2]}]);
W2[x_, n_] :=
     (-1)^(1)/(n - 1)! Sum[
    CC[n, p] p! 2^(p) (x)^(n - 2 - 2 p) , {p, 0, Floor[(n - 2)/2]}];
W1[x_, n_] := (-1)^(1)/(n - 1)! Sum[
    Binomial[n - 1, 2 p] If[p == 0, 2, 
      Binomial[p - 1/2, p - 1] (p - 1)!] 2^(p - 2) (x)^(n - 1 - 

        2 p), {p, 0, Floor[(n - 1)/2]}];
ll = Table[W2[x, n] rho[x] + W1[x, n] ( Erfc[x/Sqrt[2]]), {n, 1, M}];
Table[D[ll[[p]], {x, p}], {p, 1, Length[ll]}] // Simplify
Collect[ll, {rho[x], Erfc[x/Sqrt[2]]}, Simplify] // MatrixForm