My original problem was to find the "coefficient" functions $\varphi_{k,n}(x)$ in
$$ (\partial_x + f(x))^np(x) = \sum_{k=0}^n\varphi_{k,n}(x)\partial_x^kp(x). $$
(i.e. find the coefficients $\varphi_{k,n}$ in the formal power series expansion in $\partial_x$ of the operator $(\partial_x + f(x))^n$).
Without much difficulty one can find $$ \varphi_{k,n}(x) = \binom{n}{k}[(\partial_x + f(x))^{n-k}1] $$
This is not too helpful, however, since now I need $(\partial_x + f(x))^{\mu}1$, for which an expansion in terms of $\partial_x^k\,f(x)$ is impossible due to cross terms. I would love to have a most succinct way of describing $(\partial_x + f(x))^{\mu}1$. At worst I predict that it will end up looking vaguely similar in form (and complexity) to Faà di Bruno's formula for the coefficients in a series of composed functions. At best we may be able to have a closed form for the real numbers that appear in the expression - they seem eerily similar to Stirling numbers of the second kind. Any representation of this, $(\partial_x + f(x))^{\mu}1$, for arbitrary integral $\mu >0$ and analytic $f$, though, I would be very grateful for. The first few iterations $\mu = 0$ to $\mu = 5$ are given below. Thank you so much for any help whatsoever! $$ 1\\ f\\ f'+f^2\\ f''+3 f f'+f^3\\ f'''+4 f f''+6 f^2 f'+3 f'^2+f^4\\ f'''' + 5 f''' f + 10 f^2 f'' + 10 f^3 f' + 15 ff'^2 + 10 f' f'' + f^5 $$
PS. In case anyone is interested in context, I need this because I'd like to know, as an 'end-goal', $q(x)$ in $e^{g(x)\partial_x + f(x)}\psi(x) = \psi(q(x))$ given $g$ and $f$. I have already done this for $f = 0$, but more general differential operators involving mixed powers of the derivative are much harder, it seems. Such an expansion would be helpful in proving supporting theorems describing a more general approach to a more general problem that I play with in my spare time.
An integrating factor seems to help: $$ (\partial_{x}+f)g = g'+fg = e^{-\int fdx}\partial_{x}(e^{\int fdx}g)=(e^{-\int fdx}\partial_{x}e^{\int fdx})g, $$ where $e^{-\int fdx}$ is treated as a multiplication operator. Then $$ (\partial_x+f)^{n} = e^{-\int fdx}\partial_x^n e^{\int f dx}. $$ This helps with your stated end goal. For example, this exponential series looks like a power series, which makes sense for holomorphic functions $f$ and $g$: $$ e^{t(\partial_{x}+f)}g=e^{-\int_{0}^{x} f(r)dr}\left[\sum_{n=0}^{\infty}\frac{t^{n}}{n!}\partial_{x}^{n}\right]e^{\int_{0}^{x} f(r)dr}g(x) = e^{-\int_{0}^{x} f(r)dr}(e^{\int_{0}^{x+t}f(r)dr}g(x+t))=e^{\int_{x}^{x+t}f(r)dr}g(x+t). $$ You're looking at the case where $t=1$. As for your immediate problem, $$ (\partial_x+f)^{n}g=e^{-\int fdx}\frac{d^{n}}{dx^{n}}(e^{\int fdx}g) = e^{-\int fdx}\sum_{k=0}^{n}{{n}\choose{k}}\left[\frac{d^{n-k}}{x^{n-k}}(e^{\int fdx})\right]\partial_{x}^{k}g. $$ But it doesn't really help because the coefficients you want are as expected: $$ \varphi_{k,n} = {{n}\choose{k}}e^{-\int fdx}\frac{d^{n-k}}{dx^{n-k}}(e^{\int fdx})={{n}\choose{k}}(\partial_{x}+f)^{n-k}1. $$ It's explicit, but maybe not so helpful.
These things are related to $C^0$-semigroup theory where one forms $e^{tA}$. Equivalently, you are solving a partial differential equation in $t$ and $x$ because $$ y(x,t)=e^{tA}u(x) $$ solves $$\frac{\partial y}{\partial t} = Ay,\;\;\; y(x,0)=u(x). $$ In your case, you want $y(x,1)$ where $y(x,t)$ is the solution of $$ \frac{\partial y}{\partial t} = g\frac{\partial y}{\partial x}+fy,\;\;\; y(x,0)=\psi(x). $$ This is a first order x-t equation.