Computing the product $(\frac{d}{dx}+x)^n(-\frac{d}{dx}+x)^n$ I want to compute the product
$$
(\frac{d}{dx}+x)^n(-\frac{d}{dx}+x)^n,
$$
for a natural number $n$. For $n$ equal to 0 or 1, the computation is very simple but for such a low number as 2 the brute force calculation begins to be rather cumbersome and I cannot see any pattern emerging. I tried to find some connection with the Rodrigues' formula for the Hermite polynomials but I could not. 
These operators come up in the algebraic approach to the quantum harmonic oscillator.
Explicit Example
To avoid any misunderstanding, I am going to show explicitly the computation for the case $n=1$:
$$
(\frac{d}{dx}+x)(-\frac{d}{dx}+x)=-\frac{d^2}{dx^2}+1+x\frac{d}{dx}-x\frac{d}{dx}+x^2=-\frac{d^2}{dx^2}+x^2+1.
$$
One can think of a function $f$ the operators are acting on. For example,
$$
(\frac{d}{dx}\circ x) f= (\frac{d}{dx}x)f+x\frac{d}{dx}f=(1+\frac{d}{dx})f,  
$$
then
$$
\frac{d}{dx}\circ x=1+x\frac{d}{dx}
$$
 A: For convenience let us rewrite $x,\partial_x$ with $a,b$ so $[a,b]=x\partial_x-\partial_x x=-1$ (in the operator sense on the Schwartz space).

Lemma. For all $n\in\mathbb N$
$$
[(a+b)^n,a-b]=2n(a+b)^{n-1}
$$

Proof. As darij pointed out, one has $[a+b,a-b]=2$ (i.e. the case $n=1$). The trick then is
$$
\begin{split}
[(a+b)^{n+1},a-b]&=(a+b)[(a+b)^n,a-b]+[a+b,a-b](a+b)^n\\
&=(a+b)2n(a+b)^{n-1}+2(a+b)^{n}=2(n+1)(a+b)^{(n+1)-1}
\end{split}
$$
which concludes the proof via induction. $\square$

Proposition. For all $n\in\mathbb N_0$
$$
(a+b)^n(a-b)^n=\prod_{j=1}^n (a^2-b^2+(2j-1))
$$

Proof. ($n=0$ is obvious). Note that $(a-b)(a+b)=a^2-b^2-1$. Using the previous lemma
$$
\begin{split}
(a+b)^{n+1}(a-b)^{n+1}&=[(a+b)^{n+1},a-b](a-b)^n+(a-b)(a+b)(a+b)^n(a-b)^n\\
&=2(n+1)(a+b)^n(a-b)^n +(a^2-b^2-1)(a+b)^n(a-b)^n\\
&=\big( a^2-b^2+2n+1)(a+b)^n(a-b)^n=\prod_{j=1}^{n+1} (a^2-b^2+(2j-1))
\end{split}
$$
which again concludes the proof via induction. $\square$
This result reproduces the cases (aside from $n=0$, obvious)


*

*$(a+b)(a-b)=a^2-b^2+1$

*Making use of $[a^2,b^2]=-4ab-2$ (similar techniques) one gets
$$
\begin{split}
(a+b)^2(a-b)^2=(a^2-b^2+1)(a^2-b^2+3)&=a^4-a^2b^2-b^2a^2+4a^2-4b^2+b^4+3\\
&=a^4-2a^2b^2-4ab+4a^2-4b^2+b^4+1
\end{split}
$$
etc... I feel like this formula is the best thing one can hope for in terms of structure.
Edit: Thanks darij for the +200 rep!
A: Let $X, Y$ be operators which are central, meaning that $[X,Y] = XY-YX = c$ for some scalar $c$. Then


*

*As a form of binomial theorem, we have
$$ (X+Y)^n
= \sum_{\substack{a,b,m \geq 0 \\ a+b+2m=n}} \frac{n!}{a!b!m!2^m} c^m Y^b X^a
= \sum_{a=0}^{n} \binom{n}{a} P_{n-a}(c,Y)X^a, $$
where $P_n$ is defined by the following sum
$$P_n(c, x) = \sum_{m=0}^{\lfloor n/2\rfloor} \frac{n!}{(n-2m)!m!2^m} c^m x^{n-2m} = \left( c\frac{d}{dx} + x \right)^n \mathbf{1}. $$
In particular, if $c = -1$ then $P_n(-1, x) = \operatorname{He}_n(x)$, where $\operatorname{He}_n$ is the probabilists' Hermite polynomial. Similarly, if $c = 1$, then $P_n(1, x) = i^{-n} \operatorname{He}_n(ix)$.

*Under the same condition, for any polynomials $f, g$ we have
$$ f(X)g(Y) = \sum_{m \geq 0} \frac{c^m}{m!} g^{(m)}(Y)f^{(m)}(X). $$
In our case, $[\frac{d}{dx}, x] = 1$, and so, we can use both formulas to give a complicated, but still explicit expression for the product of $(\frac{d}{dx}+x)^n$ and $(-\frac{d}{dx}+x)^n$. Combining altogether,

\begin{align*}
&\left( \frac{d}{dx} + x \right)^n \left( -\frac{d}{dx} + x \right)^n \\
&= \sum_{p \geq 0} \frac{1}{p!} \Bigg( \sum_{\substack{a_i, b_i, m_i \geq 0 \\ a_i+b_i+2m_i = n-p}} \frac{(-1)^{a_2+m_2} (n!)^2}{a_1!a_2!b_1!b_2!m_1!m_2!2^{m_1+m_2}} x^{b_1+b_2} \left( \frac{d}{dx} \right)^{a_1+a_2} \Bigg).
\end{align*}

The following is a sample Mathematica code, comparing this formula with the actual answer for the case $n = 3$.
Coef[n_, c_, l_] := n!/(l[[1]]! l[[2]]! l[[3]]! 2^l[[3]]) c^l[[3]];
T[n_] := FrobeniusSolve[{1, 1, 2}, n];
(* Compute (x+d/dx) (x-d/dx)^n f(x) *)
n = 3;
Nest[Expand[x # + D[#, x]] &, Nest[Expand[x # - D[#, x]] &, f[x], n], n]
Sum[1/p! Sum[ Sum[ Coef[n, 1, l1] Coef[n, -1, l2]
 (-1)^l2[[2]] x^(l1[[1]] + l2[[1]]) D[f[x], 
 {x, l1[[2]] + l2[[2]]}], {l1, T[n - p]}], {l2, T[n - p]}], {p, 0, n}] // Expand
Clear[n];


