Expected value of $xp$ in harmonic oscillator I wanna find out the expected value of the $xp$ operator for the $n$-th excited state of the harmonic oscillator, i.e. calculate the value $\langle n|xp|n \rangle$. I express the position and momentum operators in terms of the ladder operators:
$$x = \frac{x_0}{\sqrt{2}} \left( a+a^{\dagger}\right)$$
$$p = \frac{p_0}{i\sqrt{2}} \left( a-a^{\dagger}\right)$$
Where $x_0,p_0$ are the scales, $x_0p_0 = \hbar$. Since $$\langle n|a^m|n \rangle =0$$
For every $m$, and likewise for $a^\dagger$. So the only term that doesn't yield a zero expectation value after multiplying will be the commutator term, yielding $-1$. Which still leaves me with an imaginary answer of $\frac{i \hbar}{2}$. Where's the mistake?
 A: Since you want the expectation value on a state with definite energy, say $|n\rangle$, then only $aa^\dagger$ and $a^\dagger a$ will give you some meaningful contributions (of course you can discard $a^\dagger a$ if $n=0$). So since it is an excited state you're interested in I will assume that $n\geq0$. Your expectation value is then
$$\langle n|xp|n\rangle=i\frac\hbar2\langle n|(aa^\dagger - a^\dagger a)|n\rangle.$$
Now recall that $[a,a^\dagger]=1$, whence
$$\langle n|xp|n\rangle=i\frac\hbar2\langle n|1|n\rangle = i\frac\hbar2.$$
There is nothing wrong with this, since, as someone pointed out $xp$ is not a self-adjoint operator in QM. So $\langle n|\ \cdot\ |n\rangle$, which defines a state on the C*-algebra of the observables is just a positive linear functional from operators to $\mathbb C$, which restricts to $\mathbb R$ on the self-adjoint elements.
A: Is it really a mistake? Let me see using QM. 
For calculating 
$-i \hbar \int dx \ x \ \psi^* \frac {∂ \psi}{∂x},$ 
we can integrate by parts, use the fact that the solutions vanish at $+\infty$ and $-\infty$, and obtain 
$i \hbar \int dx \frac {∂ (x\ \psi ^*)}{∂x}\psi$ 
= $i \hbar \int dx \ \psi ^* \psi$ + $i \hbar \int dx \ x \frac {∂\psi ^*}{∂x}\psi .$ 
So, look what I get: 
$-i \hbar \int dx \ x \ Re[\psi^* \frac {∂ \psi}{∂x}] = \frac {1}{2} i \hbar \int dx \ \psi ^* \psi = \frac {i \hbar }{2} .$
But, aren't the oscillator eigenfunctions real?
A: As you have already worked out, $(xp)^\dagger = p^\dagger x^\dagger = px \neq xp$, so that  the product of two self-adjoint operators is not needfully self-adjoint, and indeed from this equation we see that the product is self adjoint if and only if $p$ and $x$ commute.
To complete your thinking, how indeed do we calculate the product of the two measurements if we impart measurement $p$ first, then $x$? As you have already worked out for yourself, it can't be $\langle n|xp|n \rangle$.
Suppose our beginning quantum state is $\psi$ and impart $p$, getting measurement $\mu(p,\,\alpha)$ if the measurement forces our state into eigenvector $\psi_\alpha$, with probability $|\langle \psi_\alpha|\psi \rangle|^2$. Then, we impart $x$, and the mean value of this measurement is $\langle\psi_\alpha |x|\psi_\alpha\rangle$.
So the mean value of the product of the two measurements is:
$$\langle\psi|\left(\sum\limits_\alpha \langle\psi_\alpha |x|\psi_\alpha\rangle\,\mu(p,\,\alpha)\,|\psi_\alpha\rangle\langle\psi_\alpha|\right)|\psi\rangle$$
which we can write $\langle\psi|P_x|\psi\rangle$ where $P_x$ is the self adjoint observable:
$$P_x=\sum\limits_\alpha \langle\psi_\alpha |x|\psi_\alpha\rangle\,\mu(p,\,\alpha)\,|\psi_\alpha\rangle\langle\psi_\alpha|$$
This observable, interestingly, has the same eigenvectors as $p$ (being a superposition of projectors $|\psi_\alpha\rangle\langle\psi_\alpha|$ onto the eigenstates $\psi_\alpha$ of $p$), but with eigenvalues, i.e. measurements $\langle\psi_\alpha |x|\psi_\alpha\rangle\,\mu(p,\,\alpha)$ instead of $\mu(p,\,\alpha)$ (as is the case with the operator $p$). Naturally in the above $\alpha$ is an index variable ranging over the set of appropriate cardinality and integrals can replace the sums where appropriate.
Now, if we did measurement $x$ first, then $p$, the "observable" that would work out the appropriate mean value of the product of the measurements would be:
$$X_p=\sum\limits_\alpha \langle\tilde{\psi}_\alpha |p|\tilde{\psi}_\alpha\rangle\,\mu(x,\,\alpha)\,|\tilde{\psi}_\alpha\rangle\langle\tilde{\psi}_\alpha|$$
an operator with the same eigenvectors as $x$  (being a superposition of projectors $|\tilde{\psi}_\alpha\rangle\langle\tilde{\psi}_\alpha|$ onto the eigenstates $\tilde{\psi}_\alpha$ of $x$), but with eigenvalues, i.e. measurements $\langle\tilde{\psi}_\alpha |p|\tilde{\psi}_\alpha\rangle\,\mu(x,\,\alpha)$ instead of $\mu(x,\,\alpha)$. So in this case the mean of the product will be different from when we impart the measurements in the opposite order.
In your case, the application of these formulas is going to be tricky and needs to be handled with the appropriate theory of tempered distributions, since $x$ and $p$ do not have $\mathbf{L}^2(\mathbb{R})$ eigenfunctions. Unfortunately, this is not a calculation I have ever done, so I can't write up a full solution yet.
