Let's say I have the covariance matrix for a 4-mode system, $\sigma$, which is an $8 \times 8 $ matrix (the first moments are zero). I want to calculate the expectation value for an operator, for example, $\hat O = \hat a^{\dagger} \hat a \otimes \hat 1 \otimes \hat c^{\dagger} \hat c \otimes \hat 1 $, $\hat a$ and $\hat c$ being the annihilation operators for the first and the third mode.
I am new to QuTiP, so I don't know how to do this directly or even indirectly. Can I at least obtain the density matrix and calculate the expectation value? Or is there another way to do it? I know there is a command to calculate the Wigner function, but how do I go further than that?