A Leibniz-like formula for $(f(x) \frac{d}{dx})^n f(x)$? Let $f(x)$ be sufficiently regular (e.g. a smooth function or a formal power series in characteristic 0 etc.). In my research the following recursion made a surprising entrance
$$
f_1(x) = f(x),\ f_{n+1}(x) = f(x) f'_{n}(x)
$$
Thus I would like to understand the sequence
$$
\left(f(x) \frac{d}{dx} \right)^n f(x)
$$
This looks like a classical question that must have been already studied in dynamical systems, or Weyl algebras (say for $f(x) \in R[x]$, $R$ a commutative ring of characteristic $0$, and the derivation $\partial := f(x) \frac{d}{dx} \in A_1(R)$), or generating functions in combinatorics. But I have been unable to pinpoint it. My question is this:

Is there a known formula in the spirit of the general Leibniz formula which expresses
$$
\left(f(x) \frac{d}{dx} \right)^n f(x)
$$
in terms of $f$ and its derivatives $f',f'',\dots, f^{(n)}$?

Any references would be also very much appreciated!
 A: Revamped Feb. 12, 2022:
I posted an answer to this (perennial) question in detail in the old MO-Q "Formula for n-th iteration of dx/dt=B(x)" and pointed out a common conflation of related but distinct number arrays, all related to 'natural growth' of rooted trees (and therefore Lehmer codes, see Adler ref below, p. 12). With decreasing order of refinement, they are
1) OEIS A139002: the Connes-Moscovici weights of the Connes-Moscovici Hopf algebra, enumerating forests of 'naturally-grown' rooted trees (see Hivert et al. in Taylor's answer, eqns. 34 and 36-40)
2) OEIS A139605: coefficients for the Scherk-Comtet partition polynomials, the normal-ordered operator expansion of the diff op $(f(z)\partial_z)^n$ in terms of the derivatives of $f$ and the derivative operator
3) OEIS A145271: the refined Eulerian numbers, coefficients for the expansion of $(f(z)\partial_z)^n \; f(z)$, or, equivalently, $(f(z)\partial_z)^{n+1} \; z$, as a partition polynomial in terms of the derivatives of $f(z)$ (see Hivert et al., eqns. 15-20). This is the expansion the OP, M.G., is addressing.
The paper by Hivert et al. presented by Peter Taylor addresses item 1 and illustrates item 3; however, there is no discussion of a multinomial-type formula for the refined Eulerian numbers--only eqns. 28 and 35 are given, essentially stating that item 3 is a coarsening/reduction of item 1. As far as I can tell, beyond eqn. 40, Hivert et al. deal with trees in relation to generic statistics generated by various codes and don't come back to the illustration of item 3. The Findstat entry created by Hivert that Taylor links to does give examples of the refined Eulerian numbers generated by the Lehmer code algorithm. This is an algorithm, rather than a formula of the multinomial-type (correct me if I'm mistaken)--formulas that are available for item 1 and for specialized Lagrange inversion partition polynomials (see notes III and IV below) generated from item 3, but not for item 3 directly. Giving the algorithm rather than a formula is analogous to saying a number array counts the number of perfect matchings of the vertices of the n-simplices and giving an algorithm that constructs the matchings and then counts them, or saying the number array enumerates a sum over the non-intersecting dissections of the convex polygons and giving a construction algorithm, without giving a numerical formula in terms of the number of vertices of the constructs. A direct formula for the refined Eulerian numbers is still an open question, I believe.
The paper by Domininci that skbmoore (any relation to eecummings?) cites deals with related but distinct expansions, giving as the most general example the associahedra partition polynomials I point out below. Although very important, this is not the expansion the OP is addressing. (Dominici doesn't address the connections to the associahedra nor flow equations nor does he go into any detail on trees or other combinatorial constructs, givibg just a ref. on trees)
Additional refs and further notes:
I) "Set partitions and integrable hierarchies" by Adler presents the first few partition polynomials of the expansion the OP desires on p. 11 and discusses them in the context of 'natural growth' sequences $T_n$. He also alludes to Lehmer codes in the context of $T_n$. Cayley's rooted trees and the related Connes-Moscovici weights are other manifestations of 'natural growth'.
II) OEIS A145271: the partition polynomials for $(g(z)\partial_z)^n g(z)$ with $g(z) = 1/f'(z)$ (so just change notation). I called (unaware of the Hivert et al. paper until now) the coefficients of the polynomials the refined Eulerian numbers since they naturally reduce to the Eulerian numbers A008292, or A123125. My blog post "A Creation Op, Scaled Flows, and Operator Identities" contains detail on related flow functions, p.d.e.s, and more as a prelude to introducing the action of the iterated generalized Lie derivative $q(z)+g(z)\partial_z$. See also the recent MO-Q "How are Sheffer polynomials related to Lie theory?" on connections to the Sheffer polynomials, the core polynomials related to the umbral / finite operator calculus.
III) Particular series reps for $f(z)$ in $g(z) = 1/f'(z)$ of the refined Eulerian partition polynomials lead to the classic Lagrange inversion partition polynomials (LIPs) A134685 for compositional inversion of functions and formal Taylor series / e.g.f.s with $f(z)= a_1 z + a_2 \frac{z^2}{2!}+ a_3 \frac{z^3}{3!}+...$ ; the associahedra LIPs A111785 (renormalized A133437), for formal power series / o.g.f.s with $f(z)= b_1 z + b_2 z^2+ b_3 z^3+...$ ; the LIPs A133932 with the log rep $f(z) = c_1 z + c_2 \frac{z^2}{2}+ c_3 \frac{z^3}{3}+...$ ; and finally, but not least, the noncrossing partition LIPs / free cumulant partition polynomials A134264, so prominent in free probability and related quantum theory and random matrix theory, with $f(z) = z/h(z) = z/( h_1 z + h_2 z^2+ h_3 z^3+...)$. {For more info, see, e.g., my answer to the MO-Q "Important formulas in combinatorics" and the recent posts "Ruling the Inverse Universe, the inviscid Hopf-Burgers evolution equation ..." and "A Taste of Moonshine in Free Moments".}
IV) There are multinomial-type expressions for each numerical coefficient of the normal-ordering expansion of the operator $(g(z)\partial_z)^{n}$ (item 2), which are presented in the MO-Q "Differential operator power coefficients". In addition, direct simple multinomial-type expressions exist for the coefficient of any given partition monomial for all the LIPs listed in III. However, as far as I know, no such multinomial-type formula currently exists for the refined Eulerian numbers even though each full partition polynomial can be calculated independently of the others using the corresponding partition polynomials for the other LIPs in III and the permutahedra partition polynomials of A133314 (or the o.g.f. version A263633) or calculated via a matrix computation presented in the MO-Qs "Сlosed formula for $(g\partial)^n$" and "Expansions of iterated, or nested, derivatives, or vectors--conjectured matrix computation".
V) A noncommutative version of A139605 (link therein) was presented by Kentaro Ihara in "Derivations and automorphisms on non-commutative power series".
VI) Around 1853, when Lie was about ten, Charles Graves published the elegant generalized Taylor shift formula
$$ e^{t \; g(z)\partial_z} H(z) = H[f^{(-1)}(f(z)+t)]$$
and, in the 1850s also, published the operator commutator
$$[L,R] =LR-RL= 1$$ (e.g., $L=\partial_z$ and $R=x$) from which the Graves-Lie-Pincherle commutator $$[h(L),R] = \frac{dh(L)}{dL} = h'(L)$$
can be inferred and the dual
$$[L,h(R)] = \frac{dh(R)}{dR} = h'(R).$$
For Sheffer polynomial sequences, $p_n(z)$, the lowering/annihilation/destruction op $L$ and the raising/creation op $R$ are defined by $L \; p_n(z) = n \; p_{n-1}(z)$ and $R \; p_n(z) = p_{n+1}(z)$.
Obviously, $R= g(z)\partial_z$ is the raising op for the partition polynomials formed from $(g(z)\partial_z)^n g(z)$, and, with $g(z)=1/f'(z)$,
$$[g(z)\partial_z,f(z)] = 1,$$
implying, similar to the commutator identities above, that
$$[(g(z)\partial_z)^n,f(z)] = n \; (g(z)\partial_z)^{n-1}.$$
Then we have the functional identity
$$[(g(z)\partial_z)^n,f(z)]g(z) =(g(z)\partial_z)^n f(z)g(z) - f(z)(g(z)\partial_z)^n g(z) = n \; (g(z)\partial_z)^{n-1}g(z).$$
Recall the commutator acting on a function as in
$$[\partial_z,H(z)]K(z) = \partial_zH(z)K(z) - H(z)\partial_zK(z)$$
is the Newton-Leibniz formula in disguise
$$\partial_zH(z)K(z) =  H(z)\partial_zK(z) + [\partial_z,H(z)]K(z) =  H(z)\partial_zK(z) + K(z)\partial_zH(z).$$
A: Suppose that we don't reorder the parts within terms and that we always post-multiply the $f$. So e.g. $$\begin{eqnarray*}f_1 &=& f \\
f_2 &=& (Df)f \\
f_3 &=& (D^2f)ff + (Df)(Df)f \\
f_4 &=& (D^3f)fff + (D^2f)(Df)ff + (D^2f)f(Df)f + (D^2f)(Df)ff + (Df)(D^2f)ff + (Df)(Df)(Df)f\end{eqnarray*}$$etc. Each term of $f_n$ is an ordered product of $n$ subterms which generates $n$ terms in $f_{n+1}$ by the product rule. We can label each term by the sequence of indices of the subterm whose exponent of $D$ increased; then the terms of $f_n$ are labelled by sequences $(a_1, a_2, \ldots, a_{n-1})$ where $1 \le a_i \le i$ and a sequence produces a term $$\prod_{j=1}^n D^{|\{ i\,:\, a_i = j \}|}f$$
But these sequences could equally be interpreted as inversion codes (nomenclature varies; in Wikipedia they're left inversion count vectors, in Mathworld they're inversion vectors) for permutations.
If we now allow grouping of equal terms we get that the coefficient of $\prod_i (D^i f)^{b_i}$ is the number of permutations on $n-1$ elements whose left inversion count vector contains each $i$ with multiplicity $b_i$.
The description of the FindStat statistic which I linked earlier in comments is

Number of permutations whose sorted list of non zero multiplicities of the Lehmer code is the given partition.

but its reference  Hivert, F., Novelli, J.-C., Thibon, J.-Y. Multivariate generalizations of the Foata-Schützenberger equidistribution mentions some equivalent descriptions. The same coefficients also appear in section 3.1 in the context of a different process relating to a differential equation, and attributed to Cayley. (The publishers are, ridiculously, asking 47€ for access to an 1857 paper, so I haven't been able to follow up the reference).
A: In the paper 'Nested Derivatives:' A simple method for computing series expansions of inverse function' by D. Dominici, with arXiv version
https://arxiv.org/pdf/math/0501052.pdf
has something very close to what the OP wants.  Define
$$ \cal{D}^0[f](x) = 1$$
$$ \cal{D}^n[f](x) = \frac{d}{dx}\big(f(x)\cal{D}^{n-1}[f](x)\big) $$
The only difference is that the OP has an additional multiplication by $f$ at the end of the chain of operations.  That paper has some closed-form formulas for select $f$ and is oriented towards inverse functions.
A: While not pretending to answer the OP, the following is too long to fit in a comment while it might contain elements of interest to the poster.
If $f$ is a convergent object (smooth or analytic), then you encounter these quantities while computing the flow of the vector field $X(x):=f(x)\frac{\partial}{\partial x}$. According to the well-known Lie formula, one has for small $t\in (\mathbb{R},0)$ (or  $t\in (\mathbb{C},0)$ according to your setting) that the solution to the differential equation $$\dot{x}(t)=f(x(t))~~,~x(0)=x_*$$ satisfies $$x(t)=\left(\sum_{n\geq 0}\frac{t^n}{n!}X\cdot^n x\right)|_{x:=x_*}.$$ The time-independent term $X\cdot^n x$ is given by the application of $f(x)\frac{d}{dx}$ $n$ times to the identity function $x\mapsto x$.
More generally, for a function $g$ you have $$g(x(t))=\left(\sum_{n\geq 0}\frac{t^n}{n!}X\cdot^ng(x)\right)|_{x:=x_*}.$$ So the sought quantities $\left(f(x)\frac{d}{dx}\right)^nf(x)$  are basically the coefficients of the series expansion of $f(x(\bullet))$ where $x(\bullet)$ is the flow of the vector field/derivation $X$ with initial value $x(0)=x$ (sorry for the loose notations here). Because for each fixed $t$ the change of variables $x_*\mapsto x(t)$ is a symmetry of $X$, one has finally $$f(x(t))=f(x_*)\times \frac{dx(t)}{dx_*}.$$
In case $f$ is a formal object, the flow $x(t)$ might only be defined for $t$ in a discrete lattice of $\mathbb{R}$ (or $\mathbb{C}$), or just only for $t:=0$.
A: It also related to the so-called elementary differentials which appear in the algebraic setting for Runge-Kutta methods. See for example the related chapter in the book Hairer, Wanner and Lubich.
