Asymptotic inversion of ExpIntegralEi function

Question

I'm looking at the small-x and large-x asymptotic expansions of the inverse of exponential integral $E_1$ (https://dlmf.nist.gov/6.2#E1)

$$\begin{array}{lll} E_1 & = & \int_z^\infty \frac{e^{-t}}{t}\mathrm{d}t\\ E_{1\text{small}}^{ - 1} (x) &\sim& - \log x - \log ( - \log x) - \frac{{\log ( - \log x) - 1}}{{\log x}} + \mathcal{O}\!\left( {\frac{{(\log ( - \log x))^2 }}{{\log ^2 x}}} \right)\\ E_{1\text{large}}^{ - 1} (x) &\sim& e^{ - x - \gamma } + e^{ - 2x - 2\gamma } + \frac{5}{4}e^{ - 3x - 3\gamma } + \cdots \end{array} $$

Gergo Nemes derived them by hand here:

I expressed −𝛾−𝐸1(𝑧) via dlmf.nist.gov/6.6.E2, took the exponential of both sides, expanded the exponential of the power series and finally employed series reversion to solve for 𝑧

I'm interested in improving these expansions by adding terms and generally learning how to automate the process. What is the most elegant way to do this in Mathematica?

E1[z_] := -ExpIntegralEi[-z]; (* https://dlmf.nist.gov/6.2#E1 *)

E1inv[y_] := x /. First@Flatten@Solve[{E1[x] == y, x > 0}, x];
E1invSmall[x_] := -Log[x] - Log[-Log[x]] - (Log[-Log[x]] - 1)/
  Log[x];(* + O(Log[-Log[x]]^2/Log[x]^2)*)
E1invLarge[x_] := 
 Exp[-x - EulerGamma] + Exp[-2 x - 2 EulerGamma] + 
  5/4 Exp[-3 x - 3 EulerGamma];

Plot[{E1inv[y], E1invSmall[y], E1invLarge[y]}, {y, 0, 2}, 
 PlotLabel -> 
  "Approximating \!\(\*SubscriptBox[SuperscriptBox[\(E\), \(-1\)], \
\(1\)]\)", 
 PlotLegends -> {"\!\(\*SubscriptBox[SuperscriptBox[\(E\), \(-1\)], \
\(1\)]\)", 
   "\!\(\*SubscriptBox[SubscriptBox[SuperscriptBox[\(E\), \(-1\)], \
\(1\)], \(small\)]\)", 
   "\!\(\*SubscriptBox[SubscriptBox[SuperscriptBox[\(E\), \(-1\)], \
\(1\)], \(large\)]\)"}]

edit July 21

Large x expansion:

E1[z_] := EulerGamma - ExpIntegralEi[-z];
series1 = Series[E1[z], {z, 0, 4}, Assumptions -> z > 0] // Normal;
result = InverseSeries@Series[Exp[series1], {z, 0, 4}] // Normal;
result /. (z -> Exp[z + EulerGamma]) // TraditionalForm

$$\frac{907}{240} e^{-6 z-6 \gamma }+\frac{361}{144} e^{-5 z-5 \gamma }+\frac{31}{18} e^{-4 z-4 \gamma }+\frac{5}{4} e^{-3 z-3 \gamma }+e^{-2 z-2 \gamma }+e^{-z-\gamma }$$

Answer 1

You got correct hints on Math.SE. Define:

E1[z_] := EulerGamma - ExpIntegralEi[-z];
E2[z_] := -ExpIntegralEi[-z] Exp[z];

They can be inverted like this

Series[E1[z], {z, 0, 4}, Assumptions -> z > 0]
(*-Log[z]+z-z^2/4+z^3/18-z^4/96+O[z]^5*)
InverseSeries[Series[Exp[z - z^2/4 + z^3/18 - z^4/96], {z, 0, 4}]/z]
(*1/z+(1/z)^2+5/(4 z^3)+31/(18 z^4)+361/(144 z^5)+O[1/z]^7*)

Some plot for the confirmation

ParametricPlot[{{E1[z], z}, {z, E1[z]}, 
 {z, Log[Exp[z - z^2/4 + z^3/18 - z^4/96]/z]}, 
 {Log[z],1/z + (1/z)^2 + 5/(4 z^3) + 31/(18 z^4) + 361/(144 z^5)}, 
 {z, (361 E^(-5 z))/144 + (31 E^(-4 z))/18 + (5 E^(-3 z))/4 + E^(-2 z) +
     E^-z}}, 
 {z, 0.01, 5}, AspectRatio -> 1, 
 PlotRange -> {{0, 5}, {0, 5}}]

Now let us focus on the other limit. E2 has a nice series expansion at infinity

s2 = Series[E2[x], {x, Infinity, 1}]

$$ \frac{1}{x}+O\left(\frac{1}{x^2}\right)$$

The dominant decay of E1 at large values of z comes therefore from the exponent---seek the solution in the form

$$z= -\log(y)+u(y).$$

Substitute this in the defining series expansion of E1. We get an equation of the form

$$e^u=\text{rational function of $u$}.$$

For instance, to the leading order we have

$$e^{u_0}=\frac{1}{\log(z)}.$$

Determine $u_0$, substitute back, determine the next term $u_1$, iterate.