Inverse function of $\operatorname{li}(x)$ over $x>\mu$? How can I get the inverse function of $\operatorname{li}(x)$ over $x>\mu$?
Where $$\operatorname{li}(x)=\int_{0}^{x}\frac{ds}{\ln(s)}$$ is the so-called logarithmic integral, and $\operatorname{li}(\mu)=0$.
 A: There are two answers. There is such an inverse function, and it is real analytic. However, as J. M. indicates, there is no evidence, in the usual places, that anyone has found an attractive asymptotic expansion for the inverse of the exponential integral function. There is a very careful treatment of this in PECINA
Just to include one item I like, for $x > 1,$ from 5.1.10 in Abramowitz and Stegun, we have
$$  \mbox{li} \; x = \gamma + \log \log x + \sum_{n=1}^\infty \; \frac{(\log x)^n}{n \, n!}    $$
where $\gamma = 0.5772156649...$ is the Euler-Mascheroni constant.
A: Background Information:
This answer is out of the $x>\mu$ range, but still is worth posting. Here is a “closed form” given that $0\le x\le1$ using the Inverse of the Regularized Gamma function, introduced in $1996$, used to calculate quantiles in Regularized Incomplete Gamma function based cdf distributions with a $0\le \text{cdf}(\text{quantile})\le 1$. Here are the definitions:
$$Q(a,z)=\frac{\int_z^\infty e^t t^{a-1}dt}{\int_0^\infty e^t t^{a-1}dt}=\frac{}x\implies z=Q^{-1}(a,x)$$
$$\Gamma(a,z)= \int_z^\infty e^t t^{a-1}dt $$
Since the Exponential Integral can be written as $$\text{Ei}(x)=-\Gamma(0,-x)$$
but a limit must occur to compensate for the value of $Q(0,x)=\frac{\Gamma(0,x)}{\Gamma(0)}=0$:
$$\text{Ei}(x)=-\Gamma(0,-x) =-\lim_{a\to 0}\Gamma(a)Q(a,-x)=y\implies x=\text{Ei}^{-1}(y)=-\lim_{a\to0}Q^{-1}\left(a,-\frac y{\Gamma(a)}\right)$$
Here is a plot of $$\text{Ei}\left(-\lim_{a\to0}Q^{-1}\left(a,-\frac y{\Gamma(a)}\right)\right):$$

Plot of  $$-\lim_{a\to0}Q^{-1}\left(a,-\frac y{\Gamma(a)}\right): $$

so $$\text{Ei}\left(-\lim\limits_{a\to0}Q^{-1}\left(a,-\frac y{\Gamma(a)}\right)\right)=x\le0$$
Here is a differential equation for the Inverse of the Exponential Integral ignoring constants of integration and using the Product Logarithm:
$$-\frac d{dx} \lim\limits_{a\to0}Q^{-1}\left(a,-\frac y{\Gamma(a)}\right)= -\lim\limits_{a\to0}Q^{-1}\left(a,-\frac y{\Gamma(a)}\right) e^{\lim\limits_{a\to0}Q^{-1}\left(a,-\frac y{\Gamma(a)}\right)}\implies \text W(-y’)+y=0\iff y’’y+y’^2y-y’^2=0\implies y=\text{Ei}^{-1}(x)$$
Using $\text{Ei}(\ln(x))=\text{li}(x)$:
$$x=\text{li}(y)\implies y=\text{li}^{-1}(x)=e^{-\lim_{a\to0}Q^{-1}\left(a,-\frac x{\Gamma(a)}\right)}$$
Here is a plot of $\text {li}\left(e^{-\lim\limits_{a\to0}Q^{-1}\left(a,-\frac x{\Gamma(a)}\right)} \right)=x$:

Plot of $$e^{-\lim\limits_{a\to0}Q^{-1}\left(a,-\frac x{\Gamma(a)}\right)}: $$

again which only works for $0\le y\le1$
Finally, ignoring the constant of integration
$$\frac d{dx} e^{-\lim_{a\to0}Q^{-1}\left(a,-\frac x{\Gamma(a)}\right)} =-\lim\limits_{a\to 0}Q^{-1}\left(a,\frac x{\Gamma(a)}\right)\iff y\implies y’=\ln(y)\implies y’-yy’’=0\implies y=\text{li}^{-1}(x)$$
These inverses work in Wolfram Alpha and other softwares. A domain extension of the inverses may be possible with this sum definition of $Q^{-1}(a,x)$. The inverse Logarithmic Integral and inverse Exponential Integral are quantile functions.
Also see the original strategy:

How to solve $$x=\lim\limits_{t\to0} Q^{-1}(t,t)\implies \text{Ei}(-x)=-1\implies Γ(0,x)=1?$$

Summary:
It turns out that we can simplify. Again, the inverse function does not work in $x>\mu$, but the following still is a “closed form with a limit” using Mathematica functions
$$\boxed {x=\text{Ei}(y<0)\implies y=\text {Ei}^{-1}(x<0)=-\lim_{a\to0}Q^{-1}(a,-ax)}$$
Here is a demo
and
$$\boxed{x=\text{li}(0\le y<1)\implies y=\text {li}^{-1}(x<0)=e^{-\lim\limits_{a\to0}Q^{-1}(a,-ax)}}$$
Here is a demo
where the limit is not $0$ and the series seems to have no effect on the domain of the function. The correct decimal places that the formula gives is based on the precision of $a$. These inverse functions relate to the quantile of a gamma type distribution with a limit.
Series Expansion:
Note the expansion with the Euler Gamma/Mascheroni Constant $\gamma$ with the limit evaluated. Note the branch of the inverse function being used and the above definitions
$$\boxed{\text{Ei}^{-1}(x)=-e^{x-γ}- e^{2(x-γ)}\frac 54 e^{3(x-γ)}-\frac{31}{18} e^{4(x-γ)}- \frac{361}{144}e^{5(x-γ)}-\frac{4537}{1200} e^{6(x-γ)}-…=-\sum_{n=1}^\infty a_n e^{n(x-γ)}}$$
$$\boxed{\text{li}^{-1}(x)=e^{-e^{x-γ}- e^{2(x-γ)}\frac 54 e^{3(x-γ)}-\frac{31}{18} e^{4(x-γ)}- \frac{361}{144}e^{5(x-γ)}-\frac{4537}{1200} e^{6(x-γ)}-…}=e^{-\sum_{n=1}^\infty a_n e^{n(x-γ)}}=\prod_{n=0}^\infty e^{-a_n e^{n(x-γ)}}}$$
where $a_n$ depends on the coefficients of the expansion of $Q^{-1}(a,z)$ and has no closed form yet. We used:$$\lim_{x\to0}x!^\frac1x=e^{-\gamma}$$
Please correct me and give me feedback!
