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

Question

I am trying to brush up my rusty intuition on second quantization and many-particle systems and i came across the following problem:

In 1-particle QM we have the continuity equation $$ \frac{\partial}{\partial t}\left(\psi\psi^*\right)=\frac{i\hbar}{2m}\left(\psi^*\triangle\psi-\psi\triangle\psi^*\right) $$ Now, in many-particle physics (free particles!) i also expect the spatial density operator (or: number operator in spatial basis) to somehow evolve or "diffuse", if i start with spatially non-homogeneous initial conditions. Therefore i started to wonder, what the evolution equation for this operator actually looks like. My naive expectation is a direct analogy to the wavefunction: $$ \frac{\partial}{\partial t}\left(\hat{\psi}\hat{\psi}^{\dagger}\right)=\frac{i\hbar}{2m}\left(\hat{\psi}^{\dagger}\triangle\hat{\psi}-\hat{\psi}\triangle\hat{\psi}^{\dagger}\right) $$ If I try to actually calculate it, I get: $$ i\hbar\frac{\partial}{\partial t}\left(\hat{\psi}\hat{\psi}^{\dagger}\right)=\left[\hat{\psi}\hat{\psi}^{\dagger},\hat{H}\right]=\left[\hat{\psi}\hat{\psi}^{\dagger},\frac{\hbar^{2}}{2m}\triangle\hat{\psi}\hat{\psi}^{\dagger}\right] $$ or $$ \frac{\partial}{\partial t}\left(\hat{\psi}\hat{\psi}^{\dagger}\right)=\frac{-i\hbar}{2m}\left[\hat{\psi}\hat{\psi}^{\dagger},\triangle\hat{\psi}\hat{\psi}^{\dagger}\right]=\frac{-i\hbar}{2m}\left[\hat{\psi}\hat{\psi}^{\dagger},\triangle\right]\hat{\psi}\hat{\psi}^{\dagger}-\frac{i\hbar}{2m}\triangle\left[\hat{\psi}\hat{\psi}^{\dagger},\hat{\psi}\hat{\psi}^{\dagger}\right] $$ and finally $$ \frac{\partial}{\partial t}\left(\hat{\psi}\hat{\psi}^{\dagger}\right)=\frac{-i\hbar}{2m}\left[\hat{\psi}\hat{\psi}^{\dagger},\triangle\right]\hat{\psi}\hat{\psi}^{\dagger}. $$ Now the question actually is: What is the commutator of the number operator and the Laplacian? (Why I cannot answer this myself: I have no intuition about how the Laplacian acts on a many-particle state.)

Answer 1

$\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.

Answer 2

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.