For the QuTip backend it would be helpful to get the density matrix of multiple qubits and not just one

In the QuTip backend, current the the `density_operator` method takes one qubit at a time. This prevents getting the density matrix of a joint system. This can be  easily modified to get the density matrix for a list of qubits which would make it more useful.