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

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.)

• that's.... a weird continuity equation for the single-particle case; isn't it missing a divergence? In any case, since you're happy to just hat the $\psi$s, why not keep the exact same form as the single particle but with hats on? – Emilio Pisanty Nov 25 '16 at 18:23
• The missing divergences were a typo (nabla instead of laplacian) - thanks for pointing it out. So after correcting this typo, the guess for the second quantized version actually is the same as the single-particle one with hats added. The question is just whether it is actually correct and if yes, how to prove it properly... – Wave and Matter Nov 25 '16 at 21:38
• I would normally leave it in the div-of-grad form, but that's a matter of taste. @ping me in a couple of days here if you don't get an answer. – Emilio Pisanty Nov 25 '16 at 21:47

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

• Man, second quantization is bizarre. – Emilio Pisanty Nov 26 '16 at 14:01
• @EmilioPisanty Really? But the result looks exactly like you'd expect, right? PS. I'm pretty sure I made multiple sign errors here which may or may not have cancelled out: I wrote this on a bus and don't have time to check through it right now. – Mark Mitchison Nov 26 '16 at 14:07
• I guess that that just makes it a bit more unsettling. I mean, yes, it looks formally like it's all the same, but what does it really mean? What does $\partial_t\rho + \nabla\cdot{{\bf J}} = 0$ actually tell you if the thing on the left isn't a density and the thing on the right isn't a current? (and, also, more mathematically, what the heck is an 'operator-valued generalized function'?) – Emilio Pisanty Nov 26 '16 at 14:17
• @EmilioPisanty But the thing on the left is the density, and the thing on the right is the current. To me, the many-body version above looks closer to the classical continuity equation (CE) than the probability conservation equation of single-particle QM. It expresses the fact that the Hamiltonian locally conserves particular number. If you take its average you get back exactly the classical CE, but as an operator identity it also applies to the higher moments, so you can't "cheat" particle conservation by quantum fluctuations. If you start in one number sector you always stay there. – Mark Mitchison Nov 26 '16 at 14:32
• I should try to update the answer with some more details about this, but I won't have time until Monday at the earliest. – Mark Mitchison Nov 26 '16 at 14:32

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.

• Okay, so the actual equivalent to the single-particle continuity equation is the same equation as before, but formulated for the multi-particle wave-function. This is already interesting and tells me that the evolution equation for $\hat{\psi}\hat{\psi}^{\dagger}$ is probably not the same, because the quantity i look at already carries less information than the actual probability density $\psi*\psi$ made from full many-particle wavefunctions. Nevertheless what i am really interested in is the time-evolution of the number-operator $\hat{\psi}\hat{\psi}^{\dagger}$. – Wave and Matter Nov 26 '16 at 10:59
• That's not really in the spirit of what was asked, though. Surely you can fold that 3n-dimensional divergence and current into real 3d space, with an actual many-body 3d current vector. – Emilio Pisanty Nov 26 '16 at 11:27