Ladder operator identity for $\langle n | (a+a^\dagger)^k | m \rangle$ I would like to know if there is a convenient identity (and what it is) for
$$\langle n | (a+a^\dagger)^k | m \rangle$$
where $| n \rangle, \, | m \rangle$ are energy eigenstates of a simple harmonic oscillator hamiltonian and $a, \, a^\dagger$ are annihilation and creation operators respectively. $k$ is a natural number. I've done problems for $k = 1, 2, 3$ but it's not clear to me how to generalize.
 A: Yes, provided you were willing to perform integrals involving Hermite polynomials in coordinate eigenstates connected to your Fock space.
First, recall 
$$
a+a^\dagger= \sqrt{2}~ \hat{x}.
$$
Then,
$$
\psi_m(x)\equiv \left\langle x \mid m \right\rangle = {1 \over \sqrt{2^m m!}}~ \pi^{-1/4} \exp(-x^2 / 2) H_m(x),
$$ 
so inserting a complete set of coordinate eigenstates,
$$
\langle n | (a+a^\dagger)^k | m \rangle= 2^{k/2}\langle n | \hat{x}^k | m \rangle =  2^{k/2}\int dx ~\langle n | x\rangle x^k \langle x| m\rangle \\
=  2^{k/2} \frac{1}{\sqrt \pi}\frac{1}{\sqrt{2^{m+n} n! m!}}\int dx ~e^{-x^2} x^k   H_n(x) H_m (x).  
$$
You now have to use the correspondingly messy moment identities for Hermite polynomials, but you may check your low-index  specific results for the first few ones, this way. 
A: Consider a coherent state $\exp(it(a^{\dagger}+a))| 0 \rangle$, where one can show that $$a\exp(it(a^{\dagger}+a))| 0 \rangle=it\exp(it(a^{\dagger}+a))| 0 \rangle. $$Then it's easy to work out a partial answer $$\langle n|\exp(it(a^{\dagger}+a)|0\rangle=(it)^ne^{-\frac{1}{2}t^2},$$which leads to $$\langle n|(a^{\dagger}+a)^k|0\rangle=(-i)^k\partial_t^k((it)^ne^{-\frac{1}{2}t^2})_{t=0}. $$
Presumably a similar method can yield for the general case but it is certainly more tedious. 
A: 
13 Sep 2021 edit: the original derivation was incorrect (because instead of $B_k(w+z,1,0,\dots,0)$ in $(\star)$ the quantity $B_k(w+z,0,0,\dots,0)=(w+z)^k$ appeared, which is incorrect). I have corrected the derivation. I have not double checked the updated result, but I will do so when I get a chance.

The objective is to compute the quantity:
$$
\langle n|(a+a^{\dagger})^k|m\rangle
$$
It is convenient to use coherent state techniques, so in particular we primarily need the relation:
\begin{equation}
\boxed{
\begin{aligned}
I(z,w)&\equiv\langle \bar{z} | (a+a^\dagger)^k |w\rangle\\
& = B_k(w+z,1,0,\dots,0) \,e^{wz}
\end{aligned}
}\qquad (\star)
\end{equation}
where $\langle \bar{z}|=\langle 0|e^{za}$ and  $|w\rangle=e^{wa^\dagger}|0\rangle$ are coherent states (with $a|w\rangle = w|w\rangle$ and $\langle\bar{z}|a^\dagger=z\langle\bar{z}|$, and $\langle \bar{z}|w\rangle = e^{wz}$), whereas $B_k(w+z,1,0,\dots,0)$ is a complete Bell polynomial. As usual, we normalise $[a,a^\dagger]=1$.

Derivation of ($\star$):
Since I received some questions about how to derive $(\star)$ let me point out that I made use of the complete Bell polynomial identity, $B_k(a_1,\dots,a_n) = \partial_t^k \exp(\sum_{s=1}^\infty \frac{1}{s!}a_st^s)|_{t=0}$, and the Baker-Campbell-Hausdorff formula, which reduces to $e^{X+Y}=e^X e^Y e^{-\frac{1}{2}[X,Y]}$, when $[X,Y]$ commutes with $X$ and $Y$; from the latter (since the right-hand side must be symmetric in $X,Y$ as is the left-hand side) it also follows that $e^X e^Y = e^Ye^Xe^{[X,Y]}$. In further detail, making use of these results,
\begin{equation}
\begin{aligned}
I(z,w) &\equiv\langle \bar{z} | (a+a^\dagger)^k |w\rangle\\
&=\partial_t^k\langle \bar{z} | e^{(a+a^\dagger)t} |w\rangle\big|_{t=0}\\
&=\partial_t^k\langle \bar{z} | e^{a^\dagger t} e^{at}e^{\frac{1}{2}[a,a^\dagger]t^2} |w\rangle\big|_{t=0}\\
&=\partial_t^ke^{(z+w)t+\frac{1}{2}t^2}\big|_{t=0} \langle \bar{z} |w\rangle\\
&=B_k(w+z,1,0,\dots,0) \,e^{wz}
\end{aligned}
\end{equation}

Returning to ($\star$), we can extract the quantity of interest from it by differentiating it wrt $z$ and $w$ ($n$ and $m$ times respectively) and then setting $z=w=0$; after including relevant normalisations (assuming $\langle n|m\rangle=\delta_{n,m}$):
$$
\boxed{
\begin{aligned}
\langle n|(a+a^{\dagger})^k|m\rangle
&=\frac{1}{\sqrt{n!m!}}\partial_z^n\partial_w^mI(z,w)\big|_{z,w=0}\\
&=\frac{1}{\sqrt{n!m!}}\partial_z^n\partial_w^m\,B_k(w+z,1,0,\dots,0) \,e^{wz}\big|_{z=w=0}\\
\end{aligned}
}\qquad (\star\star)
$$
To evaluate the derivatives notice primarily that,
\begin{equation}
\begin{aligned}
\partial_z^nI(z,w)\big|_{z=0} &= \partial_z^nB_k(w+z,1,0,\dots,0) \,e^{wz}\big|_{z=0}\qquad (\textrm{general Leibniz rule}\,\downarrow)\\
&=\sum_{a=0}^n\binom{n}{a}\partial_z^aB_k(w+z,1,0,\dots,0) \,\partial_z^{n-a}e^{wz}\big|_{z=0}\qquad (\textrm{compl. Bell pol. property}\,\downarrow)\\
&=\sum_{a=0}^n\binom{n}{a}B_{k-a}(w,1,0,\dots,0) \,w^{n-a}\\
\end{aligned}
\end{equation}
The last equality follows immediately from the series representation of complete Bell polynomials.
Proceeding similarly for the remaining, $\partial_w^m$, derivatives one finds,
\begin{equation}
\begin{aligned}
\partial_z^n\partial_w^m\,I(z,w)\big|_{z,w=0} 
&=\partial_w^m\sum_{a=0}^n\binom{n}{a}B_{k-a}(w,1,0,\dots,0) \,w^{n-a}\big|_{w=0} \\
&=\sum_{a=0}^n\binom{n}{a}\sum_{b=0}^m\binom{m}{b}\partial_w^bB_{k-a}(w,1,0,\dots,0) \,\partial_w^{m-b}w^{n-a}\big|_{w=0}\\
%&=\sum_{a=0}^n\binom{n}{a}\sum_{b=0}^m\binom{m}{b}B_{k-a-b}(w,1,0,\dots,0) \,\frac{(n-a)!}{(n-a-m+b)!}\big|_{w=0}\\
%&=\sum_{a=0}^n\sum_{b=0}^m\frac{n!m!}{a!(n-a)!b!(m-b)!}B_{k-a-b}(0,1,0,\dots,0) \,\frac{(n-a)!}{(n-a-m+b)!}\delta_{n-a,m-b}\\
%&=\sum_{a=0}^n\sum_{b=0}^m\frac{n!m!}{a!b!(m-b)!(n-a-m+b)!}B_{k-a-b}(0,1,0,\dots,0)\delta_{n-a,m-b} \\
%&=\sum_{a=0}^n\sum_{b=0}^m\frac{n!m!}{a!(m-n+a)!(n-a)!}B_{k-m+n-2a}(0,1,0,\dots,0)\delta_{n-a,m-b} \\
&=\sum_{a=0}^n\frac{n!m!}{a!(m-n+a)!(n-a)!}B_{k-m+n-2a}(0,1,0,\dots,0) \qquad (\dagger)
\end{aligned}
\end{equation}
Using the defining series for complete Bell polynomials one can in turn show that for $p\geq0$,
$$
B_{2p}(0,1,0,\dots,0) = \frac{(2p)!}{p!},\quad B_{2p+1}(0,1,0,\dots,0) = 0,
$$
and therefore ($\dagger$) and by extension ($\star\star$) vanishes unless a positive integer $r$ can be found such that,
$$
k=2r+m-n.
$$
Given any number eigenstates labelled by $m$ and $n$, we can regard this as a condition on $k$. (E.g., if $m=n$ then only even $k$ contributes.)
\begin{equation}
\begin{aligned}
\partial_z^n\partial_w^m\,I(z,w)\big|_{z,w=0} 
&=\sum_{a=0}^n\frac{n!m!}{a!(m-n+a)!(n-a)!}B_{2(r-a)}(0,1,0,\dots,0)\\
&=\sum_{a=0}^n\frac{n!m!}{a!(m-n+a)!(n-a)!}\frac{(2r-2a)!}{(r-a!)}
\end{aligned}
\end{equation}
The final sum over $a$ can be carried out; I used Mathematica. The result may be written in terms of a generalised hypergeometric function and Gamma functions,
$$
\partial_z^n\partial_w^m\,I(z,w)\big|_{z,w=0} = \frac{m!}{(m-n)!}\frac{2^{2r}\Gamma(r+\frac{1}{2})}{\Gamma(\frac{1}{2})} \phantom{i}_1F_{\,2}\big(-n; 1+m-n,\tfrac{1}{2}-r; \tfrac{1}{4}\big)
$$
Substituting this into ($\star\star$), and taking into account that for $k\neq 2r+m-n$ the result vanishes, we can summarise the above as follows:
\begin{equation}
\boxed{
\begin{aligned}
&\langle n|(a+a^{\dagger})^{C+m-n}|m\rangle=\\
&\quad=
\left\{
\begin{aligned}
&\frac{1}{\sqrt{n!m!}}\frac{m!}{(m-n)!}\frac{2^{C}\Gamma(\frac{C+1}{2})}{\Gamma(\frac{1}{2})} \phantom{i}_1F_{\,2}\big(-n; 1+m-n,\tfrac{1-C}{2}; \tfrac{1}{4}\big)\qquad \textrm{if $C\in2\mathbf{Z}^+$}\\
&=0\qquad \textrm{if $C\in 2\mathbf{Z}^++1$}
\end{aligned}
\right.
\end{aligned}
}
\end{equation}
A: Here's a solution inspired by @Wakabaloola.  It doesn't use Bell numbers and is possibly a little more intuitive but I was unable to obtain a sensible expression for the final sum.
We are looking for an elegant computation of $\langle m\vert \hat x^k \vert n\rangle$, where
$\vert n\rangle $ is a harmonic oscillator ket.
First:
\begin{align}
\vert s\rangle&=e^{s\hat a^\dagger}\vert 0\rangle\, ,\\
    \hat a\vert s\rangle &=s\vert s\rangle 
\end{align}
Next, we start with
\begin{align}
    I(z,w)= \langle 0\vert e^{w\hat a}(\hat a+\hat a^\dagger)^k e^{z\hat a^\dagger}\vert 0\rangle 
\end{align}
and compute
\begin{align}
     e^{w\hat a}(\hat a+\hat a^\dagger)e^{-w\hat a}= \hat a+\hat a^\dagger + w\mathbb{I}
\end{align}
using the usual BCH formula.  Moreover:
\begin{align}
     e^{w\hat a}(\hat a+\hat a^\dagger)^2e^{-w\hat a}&= 
     e^{w\hat a}(\hat a+\hat a^\dagger)e^{-w\hat a}e^{w\hat a}(\hat a+\hat a^\dagger)e^{-w\hat a}\, ,\\
     &=\left(\hat a+\hat a^\dagger + w\mathbb{I}\right)^2
\end{align}
and so by induction
\begin{align}
    e^{w\hat a}(\hat a+\hat a^\dagger)^k&= \left(\hat a+\hat a^\dagger + w\mathbb{I}\right)^k e^{w\hat a}\, .
    \end{align}
Continuing:
\begin{align}
I(z,w)   &= \langle 0\vert\left(\hat a+\hat a^\dagger + w\mathbb{I}\right)^k e^{w\hat a}\vert z\rangle\, ,\\
   &=\langle 0\vert\left(\hat a+\hat a^\dagger + w\mathbb{I}\right)^k \vert z\rangle e^{wz}\, ,\\
   &=\langle 0\vert\left(\hat a+\hat a^\dagger + w\mathbb{I}\right)^k e^{z\hat a^\dagger}\vert 0\rangle
\end{align}
since $\hat a\vert z\rangle=z\vert z\rangle.$
Next, we pass $e^{z\hat a^\dagger}$ to the left of $\left(\hat a+\hat a^\dagger + w\mathbb{I}\right)^k$
using the same trick as before:
\begin{align}
    \left(\hat a+\hat a^\dagger + w\mathbb{I}\right)e^{z\hat a^\dagger}
    &= e^{z\hat a^\dagger} e^{-z\hat a^\dagger}\left(\hat a+\hat a^\dagger + w\mathbb{I}\right)e^{z\hat a^\dagger}\, ,\\
    &=  e^{z\hat a^\dagger}\left(\hat a+\hat a^\dagger + (z+w)\mathbb{I}\right)\, ,\\
   \left(\hat a+\hat a^\dagger + w\mathbb{I}\right)^ke^{z\hat a^\dagger}  
   &=  e^{z\hat a^\dagger}\left(\hat a+\hat a^\dagger + (z+w)\mathbb{I}\right)^k\, ,\\
   \langle 0\vert\left(\hat a+\hat a^\dagger + w\mathbb{I}\right)^k e^{z\hat a^\dagger}\vert 0\rangle&= 
   \langle 0\vert e^{z\hat a^\dagger}\left(\hat a+\hat a^\dagger + (z+w)\mathbb{I}\right)^k\vert 0\rangle\, ,\\
   &=\langle 0\vert\left(\hat a+\hat a^\dagger + (z+w)\mathbb{I}\right)^k\vert 0\rangle
\end{align}
since $\langle 0\vert\hat a^\dagger=0$.  Thus, we now have
\begin{align}
    I(z,w)  &=  e^{wz}
      \langle 0\vert(z+\left(\frac{2m\omega}{\hbar}\right)^{1/2} \hat x + w)^k \vert 0\rangle\, ,\\
      &= e^{wz}
      \sum_{p=0}^k   \left(\frac{2m\omega}{\hbar}\right)^{p/2} 
      \langle 0\vert \hat x^p \vert 0\rangle (z+w)^{k-p}{k \choose p}\, , \\
      &= e^{wz}\sum_{p=0}^k   \left(\frac{2m\omega}{\hbar}\right)^{p/2} 
      \langle 0\vert \hat x^p \vert 0\rangle {k \choose p} \sum_q {{k-p}\choose{q}} z^{q}w^{k-p-q}
      \label{eq:laststep}
\end{align}
where the binomial expansion is justified since $\hat x$ commutes with $(z+w)\mathbb{I}$.
We're getting close.
We want something proportional to $\langle n\vert (\hat a+\hat a^\dagger)^k\vert m\rangle$.  We can produce
$\langle {n}\vert $ and $\vert {m}\rangle $ by taking $n$ partial derivatives of $I(z,w)$ w/r to $w$ and $m$ partial
derivatives of $I(z,w)$ w/r to $z$, respectively:
\begin{align}
    \frac{\partial^{n+m}}{\partial w^n\partial z^m}I(z,w)
    &=\langle 0\vert\hat a^n e^{w\hat a }(\hat a+\hat a^\dagger)^k (\hat a^\dagger)^me^{z\hat a^\dagger}\vert 0\rangle\, ,\\
    &= \sqrt{n!m!} \langle {n}\vert e^{w\hat a }(\hat a+\hat a^\dagger)^k e^{z\hat a^\dagger}\vert {m}\rangle \, ,\\
   \sqrt{n!m!} \langle {n}\vert e^{w\hat a }(\hat a+\hat a^\dagger)^k e^{z\hat a^\dagger}\vert {m}\rangle \Bigl\vert_{w=z=0}&= 
    \sqrt{n!m!} \langle {n}\vert (\hat a+\hat a^\dagger)^k \vert{m}\rangle \, .
    \end{align}
Expanding $e^{wz}$, we can now write
\begin{align} 
 \langle {n}\vert (\hat a+\hat a^\dagger)^k \vert {m}\rangle      
      &=\frac{1}{\sqrt{n!m!}}
 \sum_{p=0}^k \left(\frac{2m\omega}{\hbar}\right)^{p/2} \langle 0\vert x^p\vert 0\rangle\, \,
 \nonumber \\
 &\qquad\qquad\times \partial_z^n\partial_w^m \sum_{q=0}^{k-p} 
 \sum_r \frac{1}{r!}
       z^{q+r} w^{k-p-q+r}{k-p \choose q}{k \choose p}\big|_{s=w=0}
\end{align}
The derivatives are $0$ unless $q+r=n$ and $k-p-q+r=m$, or
\begin{align}
    q=\frac{k-m+n-p}{2}\, ,\qquad r=\frac{m+n+p-k}{2}\, .
\end{align}
When this is the case we have
\begin{align}
    &\partial_z^n\partial_w^m\, \sum_{q=0}^{k-p} 
 \sum_r \frac{1}{r!}
       z^{q+r} w^{k-p-q+r}{k-p \choose q}{k \choose p}\big|_{s=w=0}\, ,\\
       %%%%%%%
       &\qquad \qquad =\partial_z^n\partial_w^m\,  \frac{z^{n}}{(\frac{m+n+p-k}{2})!}
       w^{m}{k-p \choose \frac{k-m+n-p}{2} }{k \choose p}\big|_{s=w=0}\, ,\\
       &\qquad \qquad  = \frac{n!m!}{(\frac{m+n+p-k}{2})!}
      {k-p \choose \frac{k-m+n-p}{2} }{k \choose p}
\end{align}
so that
\begin{align}
   \langle n|(\hat a+\hat a^{\dagger})^k|m\rangle& 
  = \frac{n!m!}{\sqrt{n!m!}}
  \sum_{p=0}^k \left(\frac{2m\omega}{\hbar}\right)^{p/2} \langle 0 \vert x^p \vert 0\rangle \frac{1}{(\frac{m+n+p-k}{2})!}
      {k-p \choose \frac{k-m+n-p}{2} }{k \choose p}
\end{align}
There remains to evaluate $\langle 0\vert x^p \vert 0\rangle$. First note
that $\langle 0\vert\hat x^p\vert 0\rangle\ne 0$ only when $p$ is even by parity.  Next
\begin{align}
      \langle 0\vert\hat x^p\vert 0\rangle=\sqrt{\frac{m\omega}{\hbar\pi}}\int dx e^{-m\omega x^2/\hbar}  x^p\, ,
  \end{align}
so let
\begin{align}
  \xi=\sqrt{\frac{m\omega}{\hbar}}x\, ,\qquad dx= \sqrt{\frac{\hbar}{m\omega}}
  d\xi
  \end{align}
and then
\begin{align}
      \left(\frac{2m\omega}{\hbar}\right)^{p/2} \langle 0\vert \hat x^p\vert 0\rangle&= \frac{1}{\sqrt{\pi}} 
      \left(\frac{2m\omega}{\hbar}\right)^{p/2}
      \left(\frac{\hbar}{m\omega}\right)^{p/2}
      \int d\xi e^{-\xi^2} \xi^p = 
     \frac{2^{p/2}}{\sqrt{\pi}} \int d\xi e^{-\xi^2} \xi^p \, ,\\
     &= 2^{p/2} \frac{(p-1)!!}{2^{p/2}} = (p-1)!!
  \end{align}
Putting all this together:
\begin{align}
   \langle {n}\vert\hat x^k\vert{m}\rangle &=\left(\frac{\hbar}{2m\omega}\right)^{k/2} \langle {n}\vert (\hat a+\hat a^\dagger)^k\vert{m}\rangle \, ,\\
&=\left(\frac{\hbar}{m\omega}\right)^{k/2}\frac{\sqrt{n!m!}}{2^{k/2}}
\sum_{p=0,2,4\ldots}^k 
\frac{(p-1)!!}{(\frac{m+n+p-k}{2})!}
  {k-p \choose \frac{k-m+n-p}{2} }{k \choose p}  
\end{align}
