Closed form for $H(n,k)=\frac{1}{k!}\sum_{j=0}^\infty H(n-k, j)j^k$ where $H(1,k)=\frac{1}{k!}$

Question

Let $H(n,k)$ be defined such that $$H(1,k)=\frac{1}{k!}\text{, and }H(n,k)=\frac{1}{k!}\sum_{j=0}^\infty H(n-1,j)j^k$$

As pointed out in the comments, I should mention that we must define $0^0=1$ as well, to be absolutely formal about $H(n,0)$ (to fit the properties outlined in the motivation section below).

Question:

Does there exist a reasonably "closed-form" for $H(n,k)$, particularly one which admits a natural extension over $\mathbb{R}$?

By this definition, the computational complexity (of the partial sums) grows by the power of $n$ as $n$ increases, so this isn't even numerically that feasible.

Motivation: It is not to hard to show that $$\begin{align*} \sum_{k=0}^\infty H(1,k) x^k&=e^x,\\ \sum_{k=0}^\infty H(2,k) x^k&=e^{e^x},\\ \text{and in general, }\sum_{k=0}^\infty H(n,k) x^k&=e^{⋰^x}\}^{{}_n} \end{align*}$$

These are very much related to bell numbers (in fact $eB_k=H(2,k)k!$). There seems to be hint of some "closed form" for $H(3,k)$ as well, as $H(3,0)=e^e$ and $H(3,1)=e^{e+1}$, however other than those I see no other nice values.

---

My work: Of course, by means of the Gamma function, one can extend $H(n,k)$ over $k\in\mathbb{R}$, however extending over $n\in\mathbb{R}$ is rather difficult, due to the inductive nature of their definition.

However, it is clear from playing around with the partial terms that they almost provide an interpolation for us.

For instance, if we consider the sum $$\frac{1}{k!}\sum_{j=0}^J \frac{j^k}{j!},$$ (which is precisely the partial sums of $H(2,k)$), notice that when $J=1$ this sum is precisely $H(1,k)$. Hence, intuitively as $J$ goes from $1$ to $\infty$, we interpolate (albeit discretely) between $H(1,k)$ and $H(2,k)$. We can convert this into a continuous interpolation $J\mapsto x$ by the taking the limit $$\begin{align*} \lim_{J\to\infty}\frac{1}{k!}\sum_{j=0}^J \frac{j^k}{j!} - \frac{1}{k!}\sum_{j=0}^J \frac{(j+1+x)^k}{(j+1+x)!} &= \frac{1}{k!}\sum_{j=0}^\infty \left(\frac{j^k}{j!} - \frac{(j+1+x)^k}{(j+1+x)!}\right) \end{align*}$$

Thus given some monotone bijection $\phi:[1,2)\mapsto [1,\infty)$, we can guess $$x\in[1,2)\implies H(x,k)\approx\frac{1}{k!}\sum_{j=0}^\infty \left(\frac{j^k}{j!} - \frac{(j+1+\phi(x))^k}{(j+1+\phi(x))!}\right),$$ which is a good start, but seems rather contrived and arbitrary. I'm really just hoping that we get a natural extension over $\mathbb{R}$ for free by finding a more "closed-form" solution for $H$.

Answer 1

Maclaurin series gives: $$H(n,k)=\frac1{k!}\left.\frac{d^k}{dx^k}\underbrace{e^{\dots^ x}}_{n\ “e”\text s}\right|_0$$

Repeated application of a Stirling S2 formula:

$$\frac{d^jf(e^x)}{dx^j}=\sum_{m=1}^j e^{jx}s_j^{(m)}f^{(j)}(e^x)$$

gives:

$$k!H(5,k)=e^{e^{e^e}}\sum_{k_1=1}^k \sum_{k_2=1}^k \sum_{k_3=1}^k \sum_{k_4=1}^k s_k^{(k_1)}s_{k_1}^{(k_2)} s_{k_2}^{(k_3)} s_{k_3}^{(k_4)} e^{k_1}e^{k_2e}e^{k_3 e^e}$$

which is true. So, the closed form is:

$$\bbox[border: 2pt solid gray]{H(n,k)=\frac 1{k!}\underbrace{e^{\dots^e}}_{n-1\ “e”\text s}\sum_{k_{1,\dots,n-1}=1}^ks_k^{(k_1)}\prod_{j=1}^{n-2}s_{k_j}^{(k_{j+1})}e^{k_j\underbrace{e^{\dots^e}}_{j-1\ “e”\text s}}}$$

For $n=1,2$, assume an empty product and for $n=1$, the sum evaluates to $1$. This formula matches the special cases in the OP’s question