$$f(x)=r+\sum_{n=1}^{\infty} \frac{\left(\ln a \right)^{n-1}\left(\ln \left(a^r \right)\right)^{nx}\left(1-r\right)^n B_n^x}{n!}$$B_n^{x-1}}{n!}$$Where B_n^x are the Bell numbers of x-th order and r=\frac{W(-\log (a))}{\log (a)} (W(x) is the Lambert function). Here: http://arxiv.org/abs/0812.4047 one can read about Bell numbers of higher orders. The problem is that Bell numbers are only defined for integer order. We can easily generalize that to any real number by induction as follows:$$A_0^x=1A_{n+1}^x=\sum_{k=0}^{x-1} A_n^x\star A_n^k$$And then$$B_n^x=A_{n-1}^{x+1}$$where f(n)\star g(n) is the binomial convolution as described by David Knuth:$$f(n)\star g(n)=\sum_{k=0}^n \binom nkf(n-k)g(k)$$To obtain the value for any real x, we can note that the right part in A_{n+1}^x=\sum_{k=0}^{x-1} A_n^x\star A_n^k is a polynomial of x and k of degree n-1 and integer coefficients and we can take indefinite sum of it symbolically following the rule$$\sum_x cx^n=\frac{B_{c+1}(x)}{c+1}$$Where B_c(x) are the Bernoulli polynomials. Unfortunately this method also works only for a \le e^{1/e} in f(x+1)=a^{f(x)}. Here is the plot of the function, for a=\sqrt{2}, obtained with this method and 5 terms: 6 deleted 4 characters in body In addition to the above formulas, we can also use this very old formula, dating back to 1945 ( J. Ginsburg, Iterated exponentials, Scripta Math. 11 (1945), 340-353.):$$f(x)=r+\sum_{n=1}^{\infty} \frac{\left(\ln a \right)^{n-1}\left(\ln \left(a^r \right)\right)^{nx}\left(1-r\right)^n B_n^x}{n!}$$Where B_n^x are the Bell numbers of x-th order and r=\frac{W(-\log (a))}{\log (a)} (W(x) is the Lambert function). Here: http://arxiv.org/abs/0812.4047 one can read about Bell numbers of higher orders. The problem is that Bell polynomials numbers are only defined for integer order. We can easily generalize that to any real number by induction as follows:$$A_0^x=1A_{n+1}^x=\sum_{k=0}^{x-1} A_n^x\star A_n^k$$And then$$B_n^x=A_{n-1}^{x+1}$$where f(n)\star g(n) is the binomial convolution as described by David Knuth:$$f(n)\star g(n)=\sum_{k=0}^n \binom nkf(n-k)g(k)$$To obtain the value for any real x, we can note that the right part in A_{n+1}^x=\sum_{k=0}^{x-1} A_n^x\star A_n^k is a polynomial of x and k of degree n-1 and integer coefficients and we can take indefinite sum of it symbolically following the rule$$\sum_x cx^n=\frac{B_{c+1}(x)}{c+1}$$Where B_c(x) are the Bernoulli polynomials. Unfortunately this method also works only for a \le e^{1/e} in f(x+1)=a^{f(x)}. Here is the plot of the function, for a=\sqrt{2}, obtained with this method and 5 terms: 5 added 37 characters in body In addition to the above formulas, we can also use this very old formula, dating back to 1945 ( J. Ginsburg, Iterated exponentials, Scripta Math. 11 (1945), 340-353.):$$f(x)=r+\sum_{n=1}^{\infty} \frac{\left(\ln a \right)^{n-1}\left(\ln \left(a^r \right)\right)^{nx}\left(1-r\right)^n B_n^x}{n!}$$Where B_n^x are the Bell numbers of x-th order and r=\frac{W(-\log (a))}{\log (a)} (W(x) is the Lambert function). Here: http://arxiv.org/abs/0812.4047 one can read about Bell numbers of higher orders. The problem is that Bell polynomials are only defined for integer order. We can easily generalize that to any real number by induction as follows:$$B_1^x=1$$A_0^x=1$$ $$B_{n+1}^x=\sum_{k=0}^{x-1} B_n^x\star B_n^k$$$A_{n+1}^x=\sum_{k=0}^{x-1} A_n^x\star A_n^k$$And then$$B_n^x=A_{n-1}^{x+1}$$where f(n)\star g(n) is the binomial convolution as described by David Knuth:$$f(n)\star g(n)=\sum_{k=0}^n \binom nkf(n-k)g(k)$$To obtain the value for any real x, we can note that the right part in B_{n+1}^x=\sum_{k=0}^{x-1} B_n^x\star B_n^k A_{n+1}^x=\sum_{k=0}^{x-1} A_n^x\star A_n^k is a polynomial of x and k of degree n-1 and integer coefficients and we can take indefinite sum of it symbolically following the rule$$\sum_x cx^n=\frac{B_{c+1}(x)}{c+1}$$Where B_c(x) are the Bernoulli polynomials. Unfortunately this method also works only for$a \le e^{1/e}$in$f(x+1)=a^{f(x)}$. Here is the plot of the function, for$a=\sqrt{2}\$, obtained with this method and 5 terms: