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.
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.
Objective:
The objective is to compute the quantity: $$ \langle n|(a+a^{\dagger})^k|m\rangle $$ the Fock states, $|m\rangle$, normalised such that $\langle n|m\rangle=\delta_{n,m}$, and $[a,a^\dagger]=1$.
Preliminary Result:
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}
Main Calculation:
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, I find it convenient to work in terms of cycle index polynomials, $Z_k(a_s)\equiv Z_k(a_1,a_2,\dots,a_k)$, defined by: $$ Z_k(a_1,a_2,\dots,a_k) = \frac{1}{k!}B_k(0!a_1,1!a_2,\dots,(k-1)!a_k). $$ All the properties of cycle index polynomials we need are given in Appendix B of the long paper. In particular, we need the following result: \begin{equation} \begin{aligned} \partial_z^aZ_k(z+w,1,0,\dots,0)&=Z_{k-a}(z+w,1,0,\dots,0), \end{aligned} \end{equation} which follows immediately from the derivative relation in the long paper, and also: $$ Z_{2p}(0,1,0,\dots,0)=\frac{2^{-p}}{p!},\quad Z_{2p+1}(0,1,0,\dots,0) = 0,\qquad(\star\star\star) $$ for some $p=0,1,2,\dots$, which follow immediately, e.g., from the contour integral representation in equation (B.677a) in the long paper.
Plugging these into the above, some elementary manipulations then lead to: \begin{equation} \begin{aligned} \langle n|(a+a^{\dagger})^k|m\rangle &=\frac{k!}{\sqrt{n!m!}}\partial_z^n\partial_w^m\,Z_k(w+z,1,0,\dots,0) \,e^{wz}\big|_{z=w=0}\\ &=\frac{k!}{\sqrt{n!m!}}\partial_z^n\sum_{b=0}^m\binom{m}{b}\big(\partial_w^b Z_k(w+z,1,0,\dots,0)\big)\big(\partial_w^{m-b}e^{wz}\big)\big|_{z=w=0}\\ &=\frac{k!}{\sqrt{n!m!}}\partial_z^n\sum_{b=0}^m\binom{m}{b} Z_{k-b}(z,1,0,\dots,0)z^{m-b}\big|_{z=0}\\ \end{aligned} \end{equation} where we just used the usual Leibniz rule for differentiating a product $m$ times. Similarly, carrying out the remaining derivatives, \begin{equation} \begin{aligned} \langle n|(a+a^{\dagger})^k|m\rangle &=\frac{k!}{\sqrt{n!m!}}\sum_{a=0}^n\binom{n}{a}\sum_{b=0}^m\binom{m}{b} \big(\partial_z^aZ_{k-b}(z,1,0,\dots,0)\big)\big(\partial_z^{n-a}z^{m-b}\big)\big|_{z=0}\\ &=\frac{k!}{\sqrt{n!m!}}\sum_{a=0}^n\binom{n}{a}\sum_{b=0}^m\binom{m}{b} Z_{k-a-b}(0,1,0,\dots,0)(m-b)!\delta_{n-a,m-b}\\ &=\sum_{a=0}^n\frac{k!\sqrt{n!m!}}{(n-a)!a!(m-n+a)!} Z_{k-m+n-2a}(0,1,0,\dots,0)\\ \end{aligned} \end{equation} where in the last line we used the the Kronecker delta to carry out the sum over $b$. From $(\star\star\star)$ and the above relation, it is clear that unless a positive integer, $r$, can be found such that, $$ k=2r+m-n, $$ the quantity $\langle n|(a+a^\dagger)^k|m\rangle$ will vanish. We can then carry out the sum over $a$, taking into account $(\star\star\star)$, to arrive at: $$ \begin{aligned} \langle n|(a+a^{\dagger})^k|m\rangle &=\sum_{a=0}^n\frac{(2r+m-n)!\sqrt{n!m!}}{(n-a)!a!(m-n+a)!} \frac{2^{-r+a}}{(r-a)!}\\ \end{aligned} $$ We can now sum over $a$ (I used Mathematica for this sum), we arrive at the final answer: $$ \boxed{ \langle n|(a+a^{\dagger})^{Q+m-n}|m\rangle = \left\{ \begin{aligned} &\frac{\sqrt{m!n!}}{2^{Q/2}(m-n)!n!(Q/2)!}\,\,\,\,F_{\!\!\!\!\!\!2\,\,\,\,\,1}\big(-n,-\tfrac{Q}{2};1+m-n;2\big),&\textrm{if}\,\, Q\in 2\mathbf{Z}^+\\ &0&\textrm{if}\,\, Q\in2\mathbf{Z}^++1\\ \end{aligned} \right. } $$ where $m\geq n$.
Consistency Check:
As a consistency check, notice that if $m=n+1$ and $Q=0$ this reduces to the expected answer: $$ \langle n|(a+a^{\dagger})|n+1\rangle = \sqrt{n+1}, $$ which can alternatively be computed using the commutator, $[a,a^\dagger]=1$, and the Fock representation of the state, $|n\rangle = \frac{1}{\sqrt{n!}}(a^\dagger)^n|0\rangle$.
References:
D. Skliros & D. Luest, Handle Operators in String Theory, Appendix B
Acknowledgement:
Thanks to @Codename47 for pointing out that there was an error in a previous result.
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.
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}
