Q:

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:

An easy derivation is to consider:
\begin{equation}
\begin{aligned}
I(z,w)&=\langle \bar{z} | (a+a^\dagger)^k |w\rangle\\
&=(w+z)^ke^{zw}
\end{aligned}
\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}|$. 
Differentiating this wrt $z$ and $w$ ($n$ and $m$ times respectively) and then setting $z=w=0$ yields the quantity of interest after including relevant normalisations (assuming $\langle n|m\rangle=\delta_{n,m}$):
$$
\boxed{\langle n|(a+a^{\dagger})^k|m\rangle=\frac{1}{\sqrt{n!m!}}\partial_z^n\partial_w^m\,(w+z)^ke^{zw}\big|_{z=w=0}}
$$
Using the binomial theorem, $(z+w)^k=\sum_{b=0}^k\binom{k}{b}z^{k-b}w^b$, and also $e^{zw}=\sum_{a=0}^{\infty}\frac{1}{a!}(zw)^a$,
\begin{equation}
\begin{aligned}
\langle n|(a+a^{\dagger})^k|m\rangle 
&= \frac{1}{\sqrt{n!m!}}\sum_{b=0}^{k}\binom{k}{b}\sum_{a=0}^{\infty}\frac{1}{a!}\big(\partial_z^n\,z^{k+a-b}\big)\big(\partial_w^m\,w^{a+b}\big)\Big|_{z=w=0}\\
&= \frac{1}{\sqrt{n!m!}}\sum_{b=0}^{k}\binom{k}{b}\sum_{a=0}^{\infty}\frac{1}{a!}\big(n!\delta_{k+a-b,n}\big)\big(m!\delta_{a+b,m}\big)
\end{aligned}
\end{equation}
Examining the range of the two sums we learn that $\langle n|(a+a^{\dagger})^k|m\rangle =0$ unless integers $N,M$ can be found such that:
$$
n=k+N-M,\qquad m=N+M,\qquad N\in\{0,1,2,\dots\},\quad M\in\{0,1,2,\dots,k\}
$$
The non-vanishing expectation values are therefore:
\begin{equation}
\begin{aligned}
\boxed{\langle k+N-M|(a+a^{\dagger})^k|N+M\rangle 
= \frac{k!\sqrt{(N+M)!(k+N-M)!}}{N!(k-M)!M!}\,\,}
\end{aligned}
\end{equation}
As a consistency check note that for $k=0$ the right-hand side equals $\delta_{M,0}$, and so $\langle N|N\rangle=1$ as expected. 

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. 

