What's a natural candidate for an analytic function that interpolates the tower function? I know that there are analytic functions whose composition with itself is the exponential function, the so-called functional square root of the exponential function, with the additional property that it is real on the real line. 
Is a 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) = e^{f(n)}, f(1) = 1$. 
 A: 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_{k-1}$ 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_{k-1}$, 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.
A: 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-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=1$$
$$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 Donald 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 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:
     (source)
A: 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 well-behaving 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.
A: Probably the only "natural" candidate so far is due to Kneser, 1949.  It probably has better properties than have been proved so far.
A: 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 and Henryk Trappmann
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)^{m-k}\exp_a^{[k]}(1)$$
$$f(x)=\lim_{n\to\infty}\binom xn\sum_{k=0}^n\frac{x-n}{x-k}\binom nk(-1)^{n-k}\exp_a^{[k]}(1)$$
$$f(x)=\lim_{n\to\infty}\frac{\sum_{k=0}^{2n} \frac{(-1)^k \exp_a^{[k]}(1)}{(x-k)k!(2n-k)!}}{\sum_{k=0}^{2n} \frac{(-1)^k }{(x-k) k!(2n-k)!}}$$
$$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 n-th iteration and $W(x)$ is the Lambert's function.
There is also an 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)}$$
