How can I find the inverse function of $f(x)=x^{x^x}$? Has this inverse function ever been defined / studied? are there any asymptotic expansions?
It would be nice if the inverse of $f(x)$ could be expressed in terms of standard mathematical functions, but it would be enough for me to know even just a few properties. The inverse of $x^{x^x}$ can be linked to the inverse of $(x + a)^x$ and to the inverse of $x e^x +a x$, and many other functions, so knowing its properties can have many applications in solving a wide range of equations
Thanks
edit:
The series expansion of $x^{x^x} $ at $x=0$ can be expressed as:
$$x^{x^x}=x \sum_{j=0}^{\infty} \frac{\left( \ln(x^x)\right)^j}{j!}B_j\left( \ln(x)\right)\tag{1}$$
$B_j(x)$ is the Bell Polynomial.
By the General Leibniz rule we have that the nth derivative of $x^{x^x}$ can be shown as
$$\sum_{k=1}^n \binom{n}{k}f^{(n-k)}(x)P_j^{(k)}(\ln(x))$$
$$f(x)=x^{j+1}\Rightarrow f^{(n-k)}(x)=x^{j+1-n+k} \frac{(j+1)!}{(j+1-n+k)!}$$ $$P_j(\ln(x))=S_{j}^{(1)}\ln(x)^{j+1}+S_{j}^{(2)}\ln(x)^{j+2}+\dots+S_{j}^{(j)}\ln(x)^{j+j}$$ $S_j^{(k)}$ is the Stirling number of the second kind
$$P^{(n)}_j(\ln(x))=\sum_{k=1}^j \ \sum_{r=n-k-j}^{n-1}\frac{(j+k)!}{(j+k-n+r)!}s_{n}^{(n-r)}S_{j}^{(k)}\frac{\ln(x)^{j+k-n+r}}{x^n}$$ $s_j^{(k)}$ is the Stirling number of the first kind
$$P_j^{(n)}(\ln(1))=P_j^{(n)}(0)=\sum_{k=1}^j (j+k)!s_{n}^{(j+k)}S_{j}^{(k)}$$
Therefore the n th derivative of $x^{x^x}$ in $x = 1$ can be expressed as: $$D_n(1)=\sum_{j=1}^n\sum_{k=1}^n\sum_{h=1}^n(j+1)\frac{(h+j)!}{(j+1+k-n)!} \binom{n}{k}s_{n}^{(j+k)}S_{j}^{(k)}$$
With $n$ between $2$ and $10$ we have $D_n(1)={2,9,32,180,954,6524,45016,360144,3023640} $ A179230 obtaining an explicit form for the Taylor expansion coefficients for $x\to 1:$
$$ x^{x^x}=f(x)=x+\sum_{n=2}^{\infty}\frac{D_{n}(1)}{n!}(x-1)^n$$
Since the series has no constant terms, it is possible to express the inverse function $f^{-1} (x)$ by Series Reversion
$$f^{-1}(x)=x-(-1+x)^2+ \frac{1}{2}(-1+x)^3+\frac{7}{6}(-1+x)^4-\frac{17}{4} (-1+x)^5+O(x^6) $$
Basic recursive method:
This is not a closed form, but it actually works. It uses a strategy of solving for a variable recursively. See this graph for proof as the recursive function converges to a vertical or horizontal line as the solution is x=recursive(x). It will be easier to calculate without splitting the bases of the logarithm:
$$x^{x^x}=y=y_0\implies x=A_0(y_0)=\log_x\left(\log_x\left(y_0\right)\right)=A_1(y_0)\implies x=\lim_{n\to\infty}A_n(y_0), A_n= \log_{A_n(y_0)}\left(\log_{A_n(y_0)}\left(y_0\right)\right)\implies A_2(y_0)= \log_{\log_x\left(\log_x\left(y_0\right)\right)}\left(\log_{\log_x\left(\log_x\left(y_0\right)\right)}\left(y_0\right)\right), A_3(y_0)= \log_{\log_{\log_x\left(\log_x\left(y_0\right)\right)}\left(\log_{\log_x\left(\log_x\left(y_0\right)\right)}\left(y_0\right)\right)}\left(\log_{\log_{\log_x\left(\log_x\left(y_0\right)\right)}\left(\log_{\log_x\left(\log_x\left(y_0\right)\right)}\left(y_0\right)\right)}\left(y_0\right)\right),…\implies x=A_\infty (x_0)$$
See graph 2 to see how to solve the the value that gets you y. The graph converges to a line as well like in paragraph 1:
$$x^{x^x}=y, B_0(a)=x=a\implies x=B_1(a)=\frac{\ln(\ln(a))-\ln\left(\ln\left(x\right)\right)}{\ln\left(x\right)}\implies x=\lim_{n\to \infty} B_n(a), B_{n+1}(a)= \frac{\ln(\ln(a))-\ln\left(\ln\left(B_n(a)\right)\right)}{\ln\left(B_n(a)\right)}\implies x=B_{\infty}(a)\implies B_2(a)= \frac{\ln(\ln(a))-\ln\left(\ln\left({\frac{\ln(\ln(a))-\ln\left(\ln\left(x\right)\right)}{ln\left(x\right)}}\right)\right)}{\ln\left({\frac{\ln(\ln(a))-\ln\left(\ln\left(x\right)\right)}{\ln\left(x\right)}}\right)}, B_3(a)= \frac{\ln(\ln(a))-\ln\left(\ln\left({ \frac{\ln(\ln(a))-\ln\left(\ln\left({\frac{\ln(\ln(a))-\ln\left(\ln\left(x\right)\right)}{\ln\left(x\right)}}\right)\right)}{\ln\left({\frac{\ln(\ln(a))-\ln\left(\ln\left(x\right)\right)}{\ln\left(x\right)}}\right)}}\right)\right)}{\ln\left({\frac{\ln(\ln(a))-\ln\left(\ln\left({\frac{\ln(\ln(a))-\ln\left(\ln\left(x\right)\right)}{\ln\left(x\right)}}\right)\right)}{\ln\left({\frac{\ln(\ln(a))-\ln\left(\ln\left(x\right)\right)}{\ln\left(x\right)}}\right)}}\right)},…\implies x=B_\infty(a)$$
$\def\srt{\operatorname{srt}}$ Analytic solutions:
With Lagrange reversion, we get:
$$y^{y^y}=z\mathop\iff^{y=e^w} w=\ln(y)e^{-we^w}\implies y=1+\sum_{n=1}^\infty\frac{\ln^n(z)}{n!}\left.\frac{d^{n-1}}{dt^{n-1}}e^{t\left(1-ne^t\right)}\right|_0$$
General Leibniz rule and a Maclaurin series of $e^y$ uses Kronecker $\delta_{m,k}$ and factorial power $a^{(b)}$:
$$\lim_{t\to0}\sum_{m=0}^\infty\sum_{k=0}^{n-1}\frac{(-n)^m}{m!}\binom{n-1}k\left.\frac{d^{n-1}t^m}{dt^{n-1}}\right|_0\left.\frac{d^{n-1}e^{(m+1)t}}{dt^{n-1}}\right|_0=\sum_{m=0}^\infty\sum_{k=0}^{n-1}\binom{n-1}k\frac{(-n)^m m^{(k)}\delta_{m,k}}{(m+1)^{k+1-n}m!}$$
Quantile mechanics (93) to (96) hints at $\delta_{m,k}$ removing the $m$ sum and replacing $m=k$:
$$\sum_{m=0}^\infty\sum_{k=0}^{n-1}\binom{n-1}k\frac{(-n)^m m^{(k)}\delta_{m,k}}{(m+1)^{k+1-n}m!}=\Gamma(n)\sum_{k=0}^{n-1}\frac{(-1)^k(k+1)^{n-k}n^k}{(k+1)!\Gamma(n-k)}$$
Therefore:
$$\boxed{\sqrt[3]z_s=\srt_3(z)=1-\sum_{n=1}^\infty\sum_{k=1}^{n-2}\frac{(-1)^kk^{n-k+1}n^{k-2}\ln^n(z)}{(n-k)!k!}}$$
shown here. Both sums are interchangeable if $k\to\infty$ and the region of convergence is near $|z|=1$. Also, @Dark Malthorp found an integral representation.