Rigorous description of the asymptotics for $\operatorname{erfi}(z)$ as $|z|\to\infty$ According to Wolfram functions, the imaginary error function admits an asymptotic expansion for $|z|\to\infty$ of the form
$$
\tag{1}
\operatorname{erfi}(z)\sim\frac{z}{\sqrt{-z^2}}+\frac{e^{z^2}}{\sqrt\pi z}{_2F_0}\left({1,\frac{1}{2}\atop -};\frac{1}{z^2}\right).
$$
Here, the hypergeometric function converges nowhere and should be understood as a formal description of a divergent asymptotic series.
I am looking for a rigorous derivation of this asymptotic expansion and, in particular, and explanation of where the $z/\sqrt{-z^2}$ term comes from.
Searching this site I came across this answer stating
$$
\operatorname{erfi}(z)\sim i+\frac{e^{z^2}}{\sqrt\pi z}{_2F_0}\left({1,\frac{1}{2}\atop -};\frac{1}{z^2}\right),
$$
which disagrees with $(1)$ since $z/\sqrt{-z^2}=\pm i$ with the sign depending on the argument of $z$. In either case, I also find is curious that the asymptotic expansion has a complex component even if $z\in\Bbb R$.
Further research indicates that the discrepancy here is due to the Stokes phenomenon, which causes the asymptotic behavior of functions to differ depending on the region of the complex plane $z$ resides; however, I am still unsure as to how this phenomenon specifically applies to $\operatorname{erfi}(z)$. Can someone please explain?
 A: From $\operatorname{erf}(z)=1-\operatorname{erfc}(z)$ and $\S7.12(\text{i})$ in the DLMF, it follows that
$$
\operatorname{erf}(z) \sim 1 - \frac{{e^{ - z^2 } }}{{\sqrt \pi  }}\sum\limits_{m = 0}^\infty  {( - 1)^m \frac{{\left( {\frac{1}{2}} \right)_m }}{{z^{2m + 1} }}} 
$$
if $z\to \infty$ and $\Re z>0$, and
$$
\operatorname{erf}(z) \sim  - 1 - \frac{{e^{ - z^2 } }}{{\sqrt \pi  }}\sum\limits_{m = 0}^\infty  {( - 1)^m \frac{{\left( {\frac{1}{2}} \right)_m }}{{z^{2m + 1} }}} 
$$
if $z\to \infty$ and $\Re z<0$. The positive and negative imaginary axes are Stokes lines. On the Stokes lines, the correct expansion is the average of the expansions on the two sides, i.e.,
$$
\operatorname{erf}(z) \sim  - \frac{{e^{ - z^2 } }}{{\sqrt \pi  }}\sum\limits_{m = 0}^\infty  {( - 1)^m \frac{{\left( {\frac{1}{2}} \right)_m }}{{z^{2m + 1} }}} 
$$
if $z\to \infty$ and $\Re z=0$. Since $\operatorname{erfi}(z)=-i\operatorname{erf}(iz)$, we can infer that
$$
\operatorname{erfi}(z) \sim \operatorname{sgn}(\Im z)i + \frac{{e^{z^2 } }}{{\sqrt \pi  }}\sum\limits_{m = 0}^\infty  {\frac{{\left( {\frac{1}{2}} \right)_m }}{{z^{2m + 1} }}} 
$$
as $z\to \infty$, with the usual convention that $\operatorname{sgn}(0)=0$.
