I know that there are analytic functions whose composition with itself is the exponential function, so called functional square root of the exponential functions, with the additional property that it is real on the real line. Is similar property possible for a holomorphic function that interpolates the tower function? Tower function on the positive integers is defined recursively by f(n+1) = exp(f(n), f(1) = 1.

First, it is obvious that you cannot have an entire function that tends to infinity at a tower rate, since an entire function that tends to infinity has to be a polynomial. More generally, the best you can hope for is to have an essential singularity at infinity such that the function converges incredibly rapidly to infinity as you approach along the real line. But this means that the coefficients converge to zero faster than exponentially. Let us try to achieve this in a minimal way. We'll choose a very very rapidly increasing sequence $n_1,n_2,\dots$ of integers and we'll choose our coefficients $a_n$ to equal $k^{n}$ when n is between $n_{k1}$ and $n_k$. Now let's estimate the value of the function $\sum a_nz^n$ when z=k. Because $n_k$ is hugely bigger than $n_{k1}$, the dominant part of this sum up to $n_k$ will be approximately $n_k$. As for the rest of the sum, it is at most $\sum_{n>n_k}(k/k+1)^n$, which is bounded above by about k, not that we really care too much (but we need it to be finite). So it looks to me as though you can get a holomorphic function to grow arbitrarily quickly to infinity along the real line. Having done that, one can surely smooth off the above construction to get the growth rate to be whatever one wants. However, the resulting function is likely to be rather artificial and perhaps not what you are hoping for. 


The function you want grows too fast to be interpolated by usual method, but there exists an iterative solution with Cauchy integrals by Dmitry Kouznetsov If you relax the condition so to find a solution for $f(x+1)=a^{f(x)}$ such that $$a \le e^{1/e} $$ then there are multiple expressions for your function: $$f(x)=\sum_{m=0}^{\infty} \binom xm \sum_{k=0}^m \binom mk (1)^{mk}\exp_a^{[k]}(1)$$ $$f(x)=\lim_{n\to\infty}\binom xn\sum_{k=0}^n\frac{xn}{xk}\binom nk(1)^{nk}\exp_a^{[k]}(1)$$ $$f(x)=\lim_{n\to\infty}\frac{\sum_{k=0}^{2n} \frac{(1)^k \exp_a^{[k]}(1)}{(xk)k!(2nk)!}}{\sum_{k=0}^{2n} \frac{(1)^k }{(xk) k!(2nk)!}}$$ $$f(x)=\lim_{n\to\infty} \log_a^{[n]}\left(\left(1\left(\ln \left(\frac{W(\ln a)}{\ln a}\right)\right)^x\right)\frac{W(\ln a)}{\ln a}+\ln \left(\frac{W(\ln a)}{\ln a}\right)\exp_a^{[n]}(1)\right)$$ Always here the number in square brackets designates nth iteration and $W(x)$ is the Lambert's function. There is also expression for inverse function: $$ f^{[1]}(x)=\lim_{n\to\infty} \frac{\ln \left(\frac{\frac{W(\ln a )}{\ln a}+\exp_a^{[n]}(x)}{\frac{W(\ln a)}{\ln a}+\exp_a^{[n]}(1)}\right)}{\ln \ln \left(\frac{W(\ln a)}{\ln a}\right)}$$ 


Probably the only "natural" candidate so far is due to Kneser, 1949. It probably has better properties than have been proved so far. 


The question is often phrased, can tetration or iterated exponentiation be naturally extended to the real and complex numbers. Using the notation $^{1}a=a, ^{2}a=a^a, ^{3}a=a^{a^a}$, how do you compute a number like $^{.5}2$, and what are the properties of $^{x}e$ ? The Derivatives of Iterated Functions Consider the smooth function $f(z): \mathbb{C} \rightarrow \mathbb{C}$ and its iterates $f^{\;\:t}(z), t \in \mathbb{N} $. The standard convention of using a coordinate translation to set a fixed point at zero is invoked, $f(0)\equiv 0$, giving $f(z)=\sum_{n=1}^{\infty} \frac{f_n}{n!} z^n$ for $0\leq z< R$ for some positive $R$. Note that $f(z)$ is the exponential generating function of the sequence $f_0, f_1, \ldots ,f_\infty$, where $f_0=0$ and $f_1$ will be written as $\lambda$. The expression $f_j^k$ denotes $(D^j f(z))^k _{z=0}$ . Note: The symbol $t$ for time assumes $t \in \mathbb{N}$, that time is discrete. This allows the variable $n$ to be used solely in the context of differentiation. Beginning with the second derivative each component will be expressed in a general form using summations and referred to here as Schroeder summations. The First Derivative The first derivative of a function at its fixed point $Df(0)=f_1$ is often represented by $\lambda$ and referred to as the multiplier or the Lyapunov characteristic number; its logarithm is known as the Lyapunov exponent. Let $g(z)=f^{t1}(z)$, then $ Df(g(z)) = f'(g(z))g'(z)$ $ = f'(f^{t1}(z))Df^{t1}(z) $ $ = \prod^{t1}_{k_1=0}f'(f^{tk_11}(z))$ $ Df^t(0) = f'(0)^t $ $ = f_1^t = \lambda^t $ The Second Derivative $D^2f(g(z)) = f''(g(z))g'(z)^2+f'(g(z))g''(z)$ $= f''(f^{t1}(z))(Df^{t1}(z))^2+f'(f^{t1}(z))D^2f^{t1}(z) $ Setting $g(z) = f^{t1}(z)$ results in $ D^2f^t(0) = f_2 \lambda^{2t2}+\lambda D^2f^{t1}(0)$. When $\lambda \neq 0$, a recurrence equation is formed that is solved as a summation. $ D^2f^t(0) = f_2\lambda^{2t2}+\lambda D^2f^{t1}(0)$ $ = \lambda^0f_2 \lambda^{2t2}$ $ +\lambda^1f_2 \lambda^{2t4}$ $+\cdots$ $+\lambda^{t2}f_2 \lambda^2$ $+\lambda^{t1}f_2 \lambda^0$ $ = f_2\sum_{k_1=0}^{t1}\lambda^{2tk_12} $ The Third Derivative Continuing on with the third derivative, $ D^3f(g(z)) = f'''(g(z))g'(z)^3+3f''(g(z))g'(z)g''(z)+f'(g(z))g'''(z)$ $ = f'''(f^{t1}(z))(Df^{t1}(z))^3 $ $ +3f''(f^{t1}(z))Df^{t1}(z)D^2f^{t1}(z)$ $ +f'(f^{t1}(z))D^3f^{t1}(z)$ $ D^3f^t(0) = f_3\lambda^{3t3}+3 f_2^2\sum_{k_1=0}^{t1}\lambda^{3tk_15} +\lambda D^3f^{t1}(0) $ $ = f_3\sum_{k_1=0}^{t1}\lambda^{3t2k_13} +3f_2^2 \sum_{k_1=0}^{t1} \sum_{k_2=0}^{tk_12} \lambda^{3t2k_1k_25} $ Note that the index $k_1$ from the second derivative is renamed $k_2$ in the final summation of the third derivative. A certain amount of renumbering is unavoidable in order to use a simple index scheme. Iterated Functions Putting the pieces together and setting the fixed point at $f_0$ gives, $f^t(z) = \sum_{j=0}^\infty D^j f^t(f_0) (zf_0)^j $ $ = f_0+\lambda^t (zf_0)+( f_2\sum_{k_1=0}^{t1}\lambda^{2tk_12}) (zf_0)^2$ $+ (f_3\sum_{k_1=0}^{t1}\lambda^{3t2k_13} +3f_2^2 \sum_{k_1=0}^{t1} \sum_{k_2=0}^{tk_12} \lambda^{3t2k_1k_25}) (zf_0)^3+ \ldots $ So far we have covered a decent amount of algebra, but still $t \in \mathbb{N}$. The equation $f^t(z)$ , $t \in \mathbb{N}$ is important because it is convergent when $f(z)$ is convergent. Hyperbolic Fixed Points When $\lambda$ is neither zero nor a root of unity $\lambda^t \neq 1, t \in \mathbb{N}$, then the nested summations simplify to $f^t(z)=f_0 + \lambda ^t (zf_0)+\frac{\lambda ^{1+t} \left(1+\lambda ^t\right) f_2}{2 (1+\lambda )} (zf_0)^2 $ $ + \frac{1}{6} \left(\frac{3 \lambda ^{2+t} \left(1+\lambda ^t\right) \left(\lambda +\lambda ^t\right) f_2^2}{(1+\lambda )^2 (1+\lambda )}+\frac{\lambda ^{1+t} \left(1+\lambda ^{2 t}\right) f_3}{1+\lambda ^2}\right) (zf_0)^3+\ldots $ Hyperbolic Tetration Let $a_0$ be a limit point for $f(z)=a^z$, so that $a^{a_0}=a_0$. Also $a_1=\lambda$. This results in a definition for tetration of complex points for all except the set of points with rationally neutral fixed points. For the real numbers $a=e^{e^{1}}\approx 1.44467, a=e^{e}\approx 0.065988 $ have rationally neutral fixed points while $a=1$ is a superattractor. All other real values of $a$ are defined by hyperbolic tetration. $ {}^t a = a_o + \lambda ^t\left(1a_o\right)+\frac{\lambda ^{1+t} \left(1+\lambda ^t\right) \text{Log}\left(a_o\right){}^2}{2 (1+\lambda )}\left(1a_o\right){}^2 $ $ + \frac{1}{6}\text{ }\left(\frac{3 \lambda ^{2+t} \left(1+\lambda ^t\right) \left(\lambda +\lambda ^t\right)\text{ }\text{Log}\left(a_o\right){}^4}{(1+\lambda )^2 (1+\lambda )}+\frac{\lambda ^{1+t} \left(1+\lambda ^t\right) \left(1+\lambda ^t\right)\text{ }\text{Log}\left(a_o\right){}^3}{(1+\lambda ) (1+\lambda )}\right)\left(1a_o\right){}^3+\ldots $ Summary One issue that some researchers have with this approach is that it results in $^x e: \mathbb{R} \rightarrow \mathbb{C} $. Because this derivation is based on the Taylor series of $f^n(z)$, if $f(z)$ is convergent then $f^n(z)$ is convergent where $n \in \mathbb{N}$. 


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), 340353.): $$f(x)=r+\sum_{n=1}^{\infty} \frac{\left(\ln a \right)^{n1}\left(\ln \left(a^r \right)\right)^{nx}\left(1r\right)^n B_n^{x1}}{n!}$$ Where $B_n^x$ are the Bell numbers of xth 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=1$$ $$A_{n+1}^x=\sum_{k=0}^{x1} A_n^x\star A_n^k$$ And then $$B_n^x=A_{n1}^{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(nk)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}^{x1} A_n^x\star A_n^k$ is a polynomial of x and k of degree n1 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: 


Dear fpqc, Check out Mathematics of Paul Erdos Vol 1 (combinatorics and algorithms), page 76. There it says a chemist friend of Erdos gave a clever construction of this functional square root: Let $f(x) = \lim_n a^n (\log_n x a)$, where $a$ is a root of $\exp(x) = x$. I don't know (or forgot my complex analysis) how to show $f$ is holomorphic in a small disk around a so I can use some help here. But then we observe that $f(\exp(x)) = af(x)$ so $\exp(x) = f^{1}(af(x))$ and to get the functional square root we replace $a$ by $\sqrt{a}$. In fact this construction allows you to get arbitrary root of the exponential function and any other wellbehaving functions. There is some German mathmematician in the 50s who constructed such a meromorphic functional square root with the property that it's real on the real line. Note this is not an answer to my original question. 

