Matrix element of powers of position operator for quantum harmonic oscillator

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 & = alphalangle m| hat{a}^{dagger} |n rangle + alpha^*langle m| hat{a}|n rangle & = alphasqrt{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.

No and yes. There is a “shortcut” of sorts using the Bargmann map where $hat a^dagger mapsto xi$ and $hat amapsto 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 nvert mapsto frac{partial_xi^n}{sqrt{n!}}, ,qquad vert nrangle mapsto frac{xi^n}{sqrt{n!}} end{align} and thus your expression can be evaluated as begin{align} langle mvert (alpha hat a^dagger + alpha^* hat a)^k vert nrangle mapsto frac{partial^m_xi}{sqrt{m!}}left(alpha, xi +alpha^* partial_xiright)^k frac{xi^n}{sqrt{n!}}biglvert_{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.

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|nrangle$. You then have $$langle m| hat x ^k|nrangle= langle m|int!! dx~ x^k ~|xrangle langle x |nrangle = 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^{talpha hat a^dagger} e^{talpha^* hat a} , $$ or the analogous expression for $hat x$ and $hat p$...

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}|nrangle$ 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/2momega}(a^dagger + a) $, so setting $ lambda = tsqrt{hbar/2momega} $ 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}|nrangle = e^{frac{lambda^2}{2}}color{blue}{langle m|e^{lambda a^dagger}}color{red}{e^{lambda a}|nrangle}. end{align} The annihilation operators satisfy $$ a^k|nrangle = left(sqrt{n-k+1}ldotssqrt{n}right),|n-krangle = left(frac{n!}{(n-k)!}right)^{1/2},|n-krangle, $$ 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}|nrangle} &color{red}{= sum_{k=0}^n frac{lambda^k}{k!} left(frac{n!}{(n-k)!} right)^{1/2} |n-krangle}. end{align}

Then

begin{align} langle m|e^{tx}|nrangle &= 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-krangle}}_{=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 {mchoose n-k} frac{left(lambda^2right)^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}|nrangle = left(frac{n!}{m!}right)^{1/2},lambda^{m-n}e^{frac{lambda^2}{2}},L_n^{(m-n)}left(-lambda^2right).$$

If $n = m$, then this reduces to $ langle n|e^{tx}|nrangle = e^{frac{lambda^2}{2}}L_nleft(-lambda^2right)$.

Matrix Elements

The matrix elements $langle m|x^p|nrangle$ can be calculated by expanding $langle m|e^{tx}|nrangle$ 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|nrangle &= 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}{2momega} 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}{2momega} right)^{p/2} sum_{k=0}^N frac{(m!,n!)^{1/2}}{2^{(p+n-m)/2-k}left(frac{p+n-m}{2}-kright)!,(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|nrangle = left( frac{hbar}{momega} right)^{p/2}cdot frac{2^{(m-n)/2}p!}{2^pleft( frac{p+n-m}{2} right)!} cdot {m choose n} cdot left( frac{n!}{m!} right)^{1/2}{}_2F_1left( -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 $prightarrow 2p$ because expectations of odd moments of $x$ vanish):

$$ boxed{langle n| x^{2p} |nrangle = left(frac{hbar}{momega}right)^pcdotfrac{(2p)!}{4^p,p!}cdot{}_2F_1(-n,-p;1;2).} $$

The equivalent form

$$ boxed{langle n| x^{2p} |nrangle = left(frac{hbar}{momega}right)^pcdotfrac{(2p)!}{4^p,p!}cdotsum_{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}|nrangle e^{-beta E_n} = left( frac{hbar}{momega} right)^pcdot frac{(2p)!}{4^p,p!}cdot coth^p left( frac{beta hbaromega}{2} right) $$

where $beta = 1/(k_BT)$ and $Z = sum e^{-beta E_n}$ is the partition function.