Matrix element of powers of position operator for quantum harmonic oscillator A similar question has been asked here before, but that did not contain the particular solution I am after and is now closed. I was wondering if there is a compact analytical formula for matrix elements of the form
$$
\langle m|\hat{x}^k|n \rangle,
$$
where $|m \rangle$ and $|n \rangle$ are the standard energy eigenstates of a quantum harmonic oscillator. I am aware that elements for a given particular $k$ can be found by using the expression $\hat{x} = \sqrt{\hbar / (2 m \omega)} (\hat{a} + \hat{a}^{\dagger})$ for the position operator in terms of ladder operators and the action of these ladder operators on the number states. However, is there a general closed-form expression for arbitrary $k$? Since, for example, there is an analytical formula for $\langle m|e^{i \hat{x}}|n \rangle$, I thought there might also be one for the matrix element that I am interested in.
More generally, is there an analytical expression for elements of the form
$$
\langle m|(\alpha \hat{a}^{\dagger} + \alpha^* \hat{a})^k|n \rangle,
$$
where $\alpha$ is a complex number? I have searched for references in textbooks or journal articles containing such expressions, but could not find any.
 A: Just to make you aware of the "bridge" expressions, the real Hermite functions (not Hermite polynomials on which they are based), $\psi_n(x)=\langle x|n\rangle$. You then have
$$\langle m| \hat x ^k|n\rangle= \langle m|\int\!\! dx~ x^k ~|x\rangle \langle x |n\rangle = \int\!\! dx~ x^k ~ \psi_m(x) \psi_n(x) . 
$$
There is a vast literature on such integrals, possibly on this site or its math sibling but it's not clear to me you prefer integrals to the recursions of @Emilio 's answer.
A further trick you could empty in your quest for the more general expression is to consider the exponential generating function of your expression,
$$
\langle m|e^{t(\alpha \hat{a}^{\dagger} + \alpha^* \hat{a})}|n \rangle,
$$
whose k t-derivative at t=0 yields your expression. Now,
$$
e^{t(\alpha \hat{a}^{\dagger} + \alpha^* \hat{a})}= e^{t^2 \alpha \alpha^* /2}  e^{t\alpha \hat a^\dagger}  e^{t\alpha^* \hat a} ,
$$
or the analogous expression for $\hat x$ and $\hat p$...
A: The answer is basically that there are formulas for these matrix elements, but they are too complicated to be useful. However, the expectation $\langle n|x^{2p}|n\rangle$ has a reasonably simple formula.
Generating Function
First of all we can calculate the exponential generating function referred to in the question and in Cosmas Zachos's answer.
The position operator is $ x = \sqrt{\hbar/2m\omega}(a^\dagger + a) $, so setting $ \lambda = t\sqrt{\hbar/2m\omega} $ gives $e^{tx} = e^{\lambda(a^\dagger + a)}$. We can use the Baker–Campbell–Hausdorff formula to normal order this operator,
\begin{align}
e^{\lambda a^\dagger}e^{\lambda a} = e^{\lambda(a^\dagger + a) + \frac{1}{2}[\lambda a^\dagger, \lambda a]} &\implies e^{\lambda(a^\dagger + a)} = e^{\frac{\lambda^2}{2}}e^{\lambda a^\dagger}e^{\lambda a}\\
&\implies \langle m|e^{tx}|n\rangle = e^{\frac{\lambda^2}{2}}\color{blue}{\langle m|e^{\lambda a^\dagger}}\color{red}{e^{\lambda a}|n\rangle}.
\end{align}
The annihilation operators satisfy
$$ a^k|n\rangle = \left(\sqrt{n-k+1}\ldots\sqrt{n}\right)\,|n-k\rangle = \left(\frac{n!}{(n-k)!}\right)^{1/2}\,|n-k\rangle, $$
so
\begin{align}
\color{blue}{\langle m|e^{\lambda a^\dagger}} &\color{blue}{= \sum_{l=0}^m \frac{\lambda^l}{l!}\left(\frac{m!}{(m-l)!}\right)^{1/2}\langle m-l|}, & \color{red}{e^{\lambda a}|n\rangle} &\color{red}{= \sum_{k=0}^n \frac{\lambda^k}{k!} \left(\frac{n!}{(n-k)!} \right)^{1/2} |n-k\rangle}.
\end{align}
Then
\begin{align}
\langle m|e^{tx}|n\rangle &= e^{\frac{\lambda^2}{2}}\color{blue}{\sum_{l=0}^m \frac{\lambda^l}{l!}\left(\frac{m!}{(m-l)!}\right)^{1/2}}\color{red}{\sum_{k=0}^n \frac{\lambda^k}{k!}\left(\frac{n!}{(n-k)!}\right)^{1/2}}\underbrace{\color{blue}{\langle m-l}|\color{red}{n-k\rangle}}_{=\delta_{l,k+m-n}}\\
&= e^{\frac{\lambda^2}{2}}\sum_{k=0}^n \frac{(m!\,n!)^{1/2}\;\lambda^{2k+m-n}}{(k+m-n)!\,(n-k)!\,k!}\tag{$\dagger$}\label{$\dagger$}\\
&= \left(\frac{n!}{m!}\right)^{1/2}\,\lambda^{m-n}e^{\frac{\lambda^2}{2}}\sum_{k=0}^n {m\choose n-k} \frac{\left(\lambda^2\right)^k}{k!}.
\end{align}
Using the definition of a generalized Laguerre polynomial,
$$ L_n^{(\alpha)}(x) = \sum_{k=0}^n (-1)^k { n+\alpha \choose n-k } \frac{x^k}{k!}, $$
we find that the generating function is
$$ \langle m|e^{tx}|n\rangle = \left(\frac{n!}{m!}\right)^{1/2}\,\lambda^{m-n}e^{\frac{\lambda^2}{2}}\,L_n^{(m-n)}\left(-\lambda^2\right).$$
If $n = m$, then this reduces to $ \langle n|e^{tx}|n\rangle = e^{\frac{\lambda^2}{2}}L_n\left(-\lambda^2\right)$.
Matrix Elements
The matrix elements $\langle m|x^p|n\rangle$ can be calculated by expanding $\langle m|e^{tx}|n\rangle$ and equating coefficients of $t^p$. Writing $ [f(t)]_{t^p}$ for the coefficient of $t^p$ in the expansion of $f(t)$, ($\ref{$\dagger$}$) gives
\begin{align}
\frac{1}{p!}\langle m|x^p|n\rangle &= \left[\color{orange}{e^{\frac{\lambda^2}{2}}}\sum_{k=0}^n \frac{(m!\,n!)^{1/2}\;\lambda^{2k+m-n}}{(k+m-n)!\,(n-k)!\,k!}\right]_{t^p}\\
&= \left( \frac{\hbar}{2m\omega} \right)^{p/2}\left[\color{orange}{\sum_{l=0}^\infty \frac{1}{l!}\left( \frac{\lambda^2}{2} \right)^l} \sum_{k=0}^n \frac{(m!\,n!)^{1/2}\;\lambda^{2k+m-n}}{(k+m-n)!\,(n-k)!\,k!}\right]_{\lambda^p}\\
&= \left( \frac{\hbar}{2m\omega} \right)^{p/2} \sum_{k=0}^N \frac{(m!\,n!)^{1/2}}{2^{(p+n-m)/2-k}\left(\frac{p+n-m}{2}-k\right)!\,(k+m-n)!\,(n-k)!\,k!}\tag{$\star$}\label{$\star$}
\end{align}
where $ N = \text{min}(n,\frac{p+n-m}{2})$. Mathematica can evaluate this sum for us to give the following expression involving the hypergeometric function:
$$ \langle m|x^p|n\rangle =  \left( \frac{\hbar}{m\omega} \right)^{p/2}\cdot \frac{2^{(m-n)/2}p!}{2^p\left( \frac{p+n-m}{2} \right)!} \cdot {m \choose n} \cdot \left( \frac{n!}{m!} \right)^{1/2}{}_2F_1\left( -n,\tfrac{m-n-p}{2};1+m-n;2 \right), $$
(yuck!) which pretty much rules out the possibility that there is a nice closed-form expression.
For the most interesting case $n=m$, we get a fairly nice result (replacing $p\rightarrow 2p$ because expectations of odd moments of $x$ vanish):
$$ \boxed{\langle n| x^{2p} |n\rangle = \left(\frac{\hbar}{m\omega}\right)^p\cdot\frac{(2p)!}{4^p\,p!}\cdot{}_2F_1(-n,-p;1;2).} $$
The equivalent form
$$ \boxed{\langle n| x^{2p} |n\rangle = \left(\frac{\hbar}{m\omega}\right)^p\cdot\frac{(2p)!}{4^p\,p!}\cdot\sum_{k=0}^\text{min(n,p)} {n \choose k}{p \choose k} 2^k,}$$
which follows from ($\ref{$\star$}$), is probably more useful (e.g. for finding asymptotics). It's also possible to use this form to work out the expectation of $x^{2p}$ in the canonical ensemble, which has a neat answer:
$$ \langle x^{2p} \rangle_\beta = \frac{1}{Z}\sum_{n=0}^\infty \langle n|x^{2p}|n\rangle e^{-\beta E_n} = \left( \frac{\hbar}{m\omega} \right)^p\cdot \frac{(2p)!}{4^p\,p!}\cdot \coth^p \left( \frac{\beta \hbar\omega}{2} \right) $$
where $\beta = 1/(k_BT)$ and $Z = \sum e^{-\beta E_n}$ is the partition function.
A: This is fully determined by the action of the creation and annihilation operators on the number basis,
\begin{align}
\hat a|n⟩ & = \sqrt{n}|n-1⟩ \\
\hat a^\dagger|n⟩ & = \sqrt{n+1}|n+1⟩.
\end{align}
Thus, for elements of linear combinations of the two, you have
\begin{align}
\langle m|\alpha \hat{a}^{\dagger} + \alpha^* \hat{a}|n \rangle
& =
\alpha\langle m| \hat{a}^{\dagger} |n \rangle
+  \alpha^*\langle m| \hat{a}|n \rangle
\\ & =
\alpha\sqrt{n+1}\delta_{m,n+1}
+ \alpha^* \sqrt{n}\delta_{m,n-1}.
\end{align}
For powers of $\hat x$, simply "rinse and repeat". This will probably get quite tedious if you want to calculate this for a high power of $\hat x$, or if you want a systematic formula that covers all of them, but if that's the case then I would suggest you to consider carefully whether you actually need it.
A: No and yes.  There is a “shortcut” of sorts using the Bargmann map
where $\hat a^\dagger \mapsto \xi$ and $\hat a\mapsto \partial_\xi$, which preserves the commutation relations.  (Note in above link $\hat a^\dagger$ is denoted by $\hat a^*$). Here $\xi$ is just a dummy variable (they use $z$ in the linked page) with no physical meaning.
Thus
\begin{align}
\langle n\vert \mapsto \frac{\partial_\xi^n}{\sqrt{n!}}\, ,\qquad 
\vert n\rangle \mapsto \frac{\xi^n}{\sqrt{n!}}
\end{align}
and thus your expression can be evaluated as
\begin{align}
\langle m\vert (\alpha \hat a^\dagger + \alpha^* \hat a)^k \vert n\rangle \mapsto
\frac{\partial^m_\xi}{\sqrt{m!}}\left(\alpha\, \xi +\alpha^* \partial_\xi\right)^k \frac{\xi^n}{\sqrt{n!}}\bigl\vert_{\xi=0}\, .
\end{align}
You still have to work through $k$ power of $\alpha \xi +\alpha^* \partial_\xi$ but if you have access to a symbolic calculator (like Mathematica) then this is actually quite quick.
Even expanding the power by hand is doable if tedious.
