Trying the following code in qutip
rho = basis(6, 5) * basis(6, 5).dag()
x = (create(6) + destroy(6)) / np.sqrt(2)
Now, expect(x, rho) gives the expected answer 0. But when I try expect(x * x, rho) I get 2.4999999999999996, instead of the expected $5.5$.
What is going wrong?
P.S. I am running qutip 4.7.0.
This answer arrives a bit late, but I will nevertheless post an answer in case someone lands on this page.
The problem here is that the physical ladder operators don't know anything about a cut-off in the number of excitations. However, destroy(6) and create(6) do. If we expand the expectation value we see
$\langle 5|(a+a^\dagger)^2|5\rangle = \langle 5|a^2+(a^\dagger)^2+aa^\dagger + a^\dagger a|5\rangle = \langle 5|aa^\dagger + a^\dagger a|5\rangle$
In QuTiP, the operators $a$ and $a^\dagger$ have a finite representation, which means that in your Hilbert space only the states from $|0\rangle$ to $|5\rangle$ exist. However, in order to evaluate the first term we need to go through the $|6\rangle$ state, which is just outside of the excitation cut-off.
There are two possible solutions: either increase the total number of excitations beyond 6, or rewrite the expression in normal order, such that
$\langle 5|aa^\dagger + a^\dagger a|5\rangle = \langle 5|[a,a^\dagger] + 2a^\dagger a|5\rangle$.
The second option is more efficient, especially for very large composite systems.
