QM Continuity Equation: Many-Body Version for Density Operator?

$\def\rr{{\bf r}} \def\ii{{\rm i}}$ There is indeed a continuity equation for the particle density $\rho(\rr)=\Psi^\dagger(\rr)\Psi(\rr),$ where the field operator $\Psi^\dagger(\rr)$ creates a particle at position $\rr$. To derive it, you need only the canonical commutation relations for the field \begin{align} [\Psi(\rr),\Psi^\dagger(\rr')]& = \delta(\rr-\rr'),\\ [\Psi(\rr),\Psi(\rr')]&=0 \end{align} together with the correct form of the Hamiltonian, which for free particles reads as $$ H = -\frac{\hbar^2}{2m}\int{\rm d} \rr\; \Psi^\dagger(\rr)\nabla^2\Psi(\rr). $$ Note that here the derivative $\nabla$ is not an operator on the space of quantum states. It acts only on operator-valued (generalised) functions like $\Psi(\rr)$, which themselves act on the space of states. Therefore, the commutator between the field operator and the derivative makes little sense and has no relevance for the problem.

The derivation of the continuity equation proceeds as follows, using the Heisenberg equation of motion for $\rho(\rr)$, \begin{align} \partial_t \rho(\rr') & =\frac{\ii}{ \hbar} [H, \Psi^\dagger(\rr')\Psi(\rr') ]\\ & = \frac{\hbar}{2m\ii} \int{\rm d}\rr\;[\Psi^\dagger(\rr)\nabla^2\Psi(\rr),\Psi^\dagger(\rr')\Psi(\rr')] \\ & =\frac{\hbar}{2m\ii} \int{\rm d}\rr\;\left \lbrace \Psi^\dagger(\rr') [\Psi^\dagger(\rr)\nabla^2\Psi(\rr),\Psi(\rr')] + [\Psi^\dagger(\rr)\nabla^2\Psi(\rr),\Psi^\dagger(\rr')] \Psi(\rr')\right\rbrace\end{align} Now, for simplicity, let's examine one of the terms above. The key is to treat each of $\Psi^\dagger(\rr)$ and $\nabla^2 \Psi(\rr)$ as separate operators, neither of which commute with $\Psi(\rr)$ in general. Nevertheless, you can use integration by parts* to shift the $\nabla^2$ around for convenience, leading to \begin{align} & \int{\rm d}\rr\; \Psi^\dagger(\rr') [\Psi^\dagger(\rr)\nabla^2\Psi(\rr),\Psi(\rr')] \\ & = \int{\rm d}\rr\; \left\lbrace\Psi^\dagger(\rr') \Psi^\dagger(\rr) [\nabla^2\Psi(\rr),\Psi(\rr')] + \Psi^\dagger(\rr') [\Psi^\dagger(\rr),\Psi(\rr')]\nabla^2\Psi(\rr) \right\rbrace \\ & = \int{\rm d}\rr\; \left\lbrace \Psi^\dagger(\rr')\nabla^2\Psi^\dagger(\rr) [\Psi(\rr),\Psi(\rr')] + \Psi^\dagger(\rr') [\Psi^\dagger(\rr),\Psi(\rr')]\nabla^2\Psi(\rr) \right\rbrace\\ & = \int{\rm d}\rr\; \left\lbrace 0 - \Psi^\dagger(\rr')\nabla^2\Psi(\rr) \delta(\rr'-\rr) \right\rbrace \\ & = -\Psi^\dagger(\rr')\nabla^2 \Psi(\rr'). \end{align} Putting it together, you should obtain $$ \partial_t\rho = -\frac{\hbar}{2m\ii}\left(\Psi^\dagger\nabla^2 \Psi - \nabla^2\Psi^\dagger\Psi\right),$$ from which one identifies the particle current operator $$ \mathbf{J} = \frac{\hbar}{2m\ii}\Psi^\dagger\nabla \Psi + {{\rm h.c.}},$$ defined such that $\partial_t\rho + \nabla\cdot{{\bf J}} = 0$.

$\,$

*One also assumes as usual that the fields vanish at infinity, or more strictly, that the Hilbert space only contains states which are annihilated by the field operators at infinity.

The many-body probability continuity equation follows from the Schrödinger equation in the same way as in the one-particle case. Assume n particles with coordinates $x_{1,i},x_{2,i},x_{3,i}$ in 3-D space with the Hamiltonian operator $$H=\sum_{i~=~1}^n\sum_{j~=~0}^3 \left[p_{i,j}^2/2m+W_i(x_{j},t)\right]+V(x_{1,1},x_{2,1},x_{3,1},x_{2,1},\ldots,t)$$ The Schrödinger equation for the multi-particle wave function $\psi(x_{1,1},x_{2,1},x_{3,1},x_{2,1},\ldots,t)$ is $$i\hbar ~\frac {\partial \psi}{\partial t}=H\psi$$ Then the time change of probability density is given by $$\frac{\partial (\psi\psi^*)}{\partial t}=\psi^*\frac{\partial \psi}{\partial t}+\psi\frac{\partial \psi^*}{\partial t}=\frac {i\hbar }{2m}~[\psi^*\Delta \psi-\psi \Delta \psi^*]=\frac {i\hbar }{2m}~ \nabla[\psi^*\nabla \psi-\psi \nabla \psi^*]$$ where $\Delta$ and $\nabla$ are the Laplace and Gradient operator, respectively, in 3N-dimensional configuration space. Thus the probability current density in 3N-d configuration space is $$\vec J=-~\frac {i\hbar }{2m}~[\psi^*\nabla\psi-\psi \nabla\psi^*]$$ and the probability continuity equation in 3N-d configuration space can be written as $$\frac{\partial (\psi\psi^*)}{\partial t}=-~\textrm{div}~{\vec J}$$ where $\textrm{div}$ is the divergence operator in 3N-dimensional configuration space.