Nicer expression for the following differential operator

Stirling numbers $\newcommand\s[2]{\left\{#1\atop #2\right\}} \def\d{\partial_t}$

As @HenningMakholm notes, it is not hard to see that $$D_n = \sum_{k=1}^n a_k^n t^k \d^k,$$ where $a_k^n \in \mathbb{N}$. Consider $D_{n+1} = D_1 D_n$. We find $$D_{n+1} = \sum_{k=1}^{n+1} a_k^{n+1} t^k \d^k$$ where $$\begin{equation*} a_k^{n+1} = a_{k-1}^n + k a_k^n \tag{1} \end{equation*}$$ with the condition $a_0^n = a_{n+1}^n = 0$ for $n\ge 1$. We also have $a_1^1 = 1$, since $D_1 = t \d$.

Using (1) it is possible to show that an equivalent set of boundary conditions is $a_0^0 = 1$ and $a_0^n = a_n^0 = 0$. This implies, as @RahulNarain notes in the comments, that the $a_k^n$ are the Stirling numbers of the second kind, $$\begin{eqnarray*} D_n &=& \sum_{k=1}^n \s{n}{k} t^k \d^k \\ &=& t \d + (2^{n-1}-1)t^2\d^2 + \ldots + \frac{1}{2}n(n-1) t^{n-1}\d^{n-1} + t^n\d^n. \end{eqnarray*}$$

It is not clear what is the ultimate goal of applying this operator to the generating function for the Hermite polynomials. One can show, for example, that $$D_n e^{2xt-t^2} = \sum_{k=1}^n \s{n}{k} \sum_{m=0}^\infty H_{m+k}(x) \frac{t^{m+k}}{m!}.$$

Change of variables

Another way to go is to let $t = e^s$. Then, $D_n = \partial_s^n$. Thus, we must calculate $$\partial_s^n \exp(2x e^s - e^{2s}) = \partial_s^n \sum_{k=0}^\infty H_k(x) \frac{e^{k s}}{k!}.$$ We find $$D_n e^{2xt-t^2} = \sum_{k=1}^\infty H_k(x) \frac{k^n t^k}{k!}.$$

The Weyl "algebra" $\mathbb{Z}[t, \partial_t]$ acts faithfully on the polynomial ring $\mathbb{Z}[t]$ in the obvious way. The latter is graded by degree, $t$ raises degree by $1$, $\partial_t$ lowers degree by $1$, and $t \partial_t$ preserves degree: in fact $$(t \partial_t) t^m = m t^m.$$

Thus $$(t \partial_t)^n t^m = m^n t^m.$$

A basis for the space of degree-preserving elements of the Weyl algebra is given by $t^k \partial_t^k, k \in \mathbb{Z}_{\ge 0}$, and $$(t^k \partial_t^k) t^m = m(m-1)...(m-k+1) t^m = (m)_k t^m$$

where $(m)_k$ denotes the falling factorial. Thus in order to find coefficients $a_{n,k}$ such that $$(t \partial_t)^n = \sum_k a_{n,k} t^k \partial_t^k$$

it is necessary and sufficient to find coefficients $a_{n,k}$ such that $$m^n = \sum_k a_{n,k} (m)_k$$

for all $m$. Since the polynomials $(m)_k$ form a basis of the space of polynomials, the coefficients $a_{n,k}$ exist uniquely, and in fact this is one way to define the Stirling numbers of the second kind.

The combinatorial interpretation is as follows: $m^n$ counts the number of functions $[n] \to [m]$, where $[n] = \{ 1, 2, ... n \}$. The above identity groups these functions together according to the size of their range: there are $(m)_k$ possible ranges of size $k$ and $a_{n,k}$ functions $[n] \to [m]$ having range a fixed subset of $[m]$ of size $k$. (This is the same thing as an equivalence relation on $[n]$ with $k$ equivalence classes by taking preimages.)

Using the Leibniz rule $\partial_t t=1+t\partial_t$ I get $$t\partial_t t^n\partial_t^n = n t^n\partial_t^n + t^{n+1}\partial_t^{n+1}$$ so your $D_n$ is $\sum_{i=1}^n a_i t^i \partial_t ^i$ for some combinatorial coefficients $a_i \in \mathbb N$.

Is that the direction you're looking for?

Here are some coefficients:

n=1   1
n=2   1   1
n=3   1   3   1
n=4   1   7   6    1
n=5   1  15  25   10    1
n=6   1  31  90   65   15   1
n=7   1  63 301  350  140  21   1
n=8   1 127 966 1701 1050 266  28   1

We see that $a_1=1$, $a_2=2^n-1$, and $a_{n-1}$ are the triangular numbers. In between it doesn't look very familiar. (But see Rahul's comment below!)