How to get more terms with the Series[] expansion of InverseErf[x] around x=1?

Question

The inverse error function, which is given by InverseErf[x], is quite important in statistics as it gives the confidence levels around a 1D Gaussian, for example for the $1$ and $2\sigma$ ($\sim 68.3\%$ and $\sim95.5\%$ respectively):

Sqrt[2] InverseErf[0.6827]
=1

Sqrt[2] InverseErf[0.9545]
=2

So, while playing with it and wanting to find some asymptotics close to probability 1 (i.e. many $\sigma$'s), I did a series expansion around $x=1$ (note in hindsight I'm doing the expansion of InverseErf[x]^2 to avoid a square-root when comparing with the result of a paper, see below) and no matter what order I asked for

Series[InverseErf[x]^2, {x, 1, -30}, Assumptions -> x < 1]

or

Series[InverseErf[x]^2, {x, 1, 30}, Assumptions -> x < 1]

Mathematica always gives back the result:

$$ \frac{1}{2} \left(-2 \log (1-x)-\log (-2 \log (1-x)+\log (2)-\log (\pi ))+\log \left(\frac{2}{\pi }\right)\right)+O\left((x-1)^2\right).~~~~~~~~~~~~(1)$$

However, the inverse error function has a well defined (albeit asymptotic) expansion around $x=1$, see for example the seminal paper by Blair et al on approximations of the error function, where in Eq.(2) they give the expression: $$ (\textrm{InverseErf}[x])^2 \sim \eta-\frac12 \ln \eta +\eta^{-1} \left(\frac14 \ln \eta-\frac12\right)+\eta^{-2} \left(\frac1{16} \ln^2 \eta-\frac38 \ln \eta +\frac78\right)+...,~~~~~~~~~~~~(2)$$ where $\eta=-\ln [\pi^{1/2}*(1-x)]$ and in this notation $\ln(x)=\log(x)=\log_e(x)$ (ie the log base e). Also, $\ln^2 \eta=[\ln(\eta)]^2$ (thanks to @user293787).

Using the same variable $\eta$, Mathematica's series expansion, given originally by Eq. (1), can be written as

$$ (\textrm{InverseErf}[x])^2 \sim \eta -\frac{1}{2} \ln \left(\eta +\frac{\ln (2)}{2}\right).~~~~~~~~~~~~(3)$$

So, my question is: why does Mathematica only stop at such low order in the expansion in terms of $x$ and how to coax it to give more terms?

Note 1: In fact the expression from the paper is much more accurate and by keeping up to the $\eta^{-1}$ term (inclusive) one can get an accuracy of $0.00013\%$ at $x \rightarrow 1 - 10^{-30}$ ($\eta\sim68.5051878$) with respect to the exact result, compared to the Series[] expansion which is "only" accurate to $\sim 0.016 \%$.

Note 2: I can in fact just code the expression from the paper and use that if I want or even use the exact arbitrary precision result from InverseErf[x], but the question is how to get Series to give more terms, mainly out of curiosity.

Note 3: The result from the paper also gives a complex $i \pi$ term, but I only consider the real part in the comparison.

Note 4: This question, only discusses the limit around $x=0$, which is already well-known via a MacLaurin series.

Note 5: Tested mainly on v11.2 on Windows but the same issue seems to persist in v13 as well.

Answer 1

This is not strictly speaking an answer, but I thought I provide code that can be used to generate the kind of expansion that OP mentions. Perhaps it is useful for other people here.

I use $y$ as an abbreviation for InverseErf[x]^2. As a function of $\eta$ it satisfies the differential equation $$ y'(\eta) = \sqrt{y} e^{y-\eta} $$ I will use this to determine the coefficients $a_{km}$ in the following ansatz for $\eta \to \infty$. This ansatz is motivated by the expansion that OP has given:

n = 6; (* order of the expansion *)
yansatz = Function[{η}, η - 1/2 * Log[η] +
                   Sum[1/η^k * Sum[a[k,m] * Log[η]^m, {m, 0, k}], {k, 1, n}]];

To plug this into the differential equation and solve for the coefficients, one can use the following code, where I am now using SolveAlways as per @J.M.'s suggestion (thanks!).

eq = {y'[η] - Sqrt[y[η]] * Exp[y[η] - η]};
eq2 = Normal[Series[eq /. {y -> yansatz}, {η, Infinity, n}]] /.
      {Log[η] -> logeta} /. {η -> 1/inveta};
yansatz[η] /. First[SolveAlways[eq2 == 0, {logeta, inveta}]]

For example, with n = 6 this gives

The first few terms are consistent with the expansion given by OP.