In 1-particle non-relativistic QM we have the continuity condition as a per definitionem property for the 1-electron probability current density for an electron in the magnetic field in a stationary state and ignoring spin

$$ \nabla\cdot{\bf j}=0$$

with the definition of

$$ {\bf j}=-\frac{i\hbar}{2m}\left(\psi^*\nabla\psi-\psi\nabla\psi^*\right) +\frac{e}{m}{\bf A}\psi^*\psi $$

My question is what happens in the case of $n$ interacting (non-separable) electrons?

The definition of the $n$-electron current is then as the diagonal of the real part of the kinematic (effective one-particle) momentum density

$$ {\bf j}=diag\{\Re[\hat\pi \rho(\vec{r},\vec{r}')]\}$$

with the $n$-particle kinematic (mechanical) momentum operator.

$$ \hat\pi = \sum_{j=1}^n \big(\hat p_j + e{\bf A}(\vec{r}_j)\big)$$ and the one-particle reduced density matrix (1P-RMD) $$ \rho(\vec{r},\vec{r}') = \int \Psi^*(\vec{r}',\vec{r}_2,\cdots,\vec{r}_2)\Psi(\vec{r},\vec{r}_2,\cdots,\vec{r}_n) d\vec{r}_2 \dots d\vec{r}_n$$

So the question is

How to derive to derive the continuity condition $$ \nabla\cdot{\bf j} = 0$$ for the n-electron effective one-particle probability current density $${\bf j}=diag\{\Re[\hat\pi \rho(\vec{r},\vec{r}')]\}?$$

This question is inspired by this very recent publication where the problem is treated by explicitly assuming a product type expansion for the $n$-electron wave function.

Then there is also this very closely related PSE question on the QM n-particle continuity condition but here without the external magnetic field. However, the highest ranked answer looks very impressive to me and is based on second quantization. So already an extension of this answer to the magnetic field case would be of interest.

I would be mostly interested if it possible to derive the result on basis of a 1P-RDM basis like for example using the von-Neumann equation for a 1P-RDM dynamics.

Starting an attempt to answer:

In analogy to the derivation of the hydrodynamical formulation of the Schrödinger mechanics we start with a "polar form" of the 1P-RDM $$ \rho(\vec{r},\vec{r}';t)=R(\vec{r},\vec{r}';t)\exp{{i\frac{S(\vec r,\vec r';t)}{\hbar}}} $$ Now the task is to insert that into the von Neumann equation, that looks for the static case like $$ -\frac{i}{\hbar}[\hat H,\hat \rho] = 0$$ (where the prefactor of course can be omitted). And about here I come to my limits. I have no clear idea how to insert the 1P-RDM into the commutator ... any help appreciated.

  • 1
    $\begingroup$ Are you silently assuming stationary states? Otherwise, why is your hoped-for continuity equation $\nabla \cdot j = 0$ and not $\nabla \cdot j = - \partial_t \rho$? $\endgroup$
    – ACuriousMind
    Commented Sep 22, 2019 at 9:57
  • $\begingroup$ @ACuriousMind "eigenstate of the full Hamiltonian (=static problem)" was my clumsy way to express that. Maybe I reword it to "for stationary states". $\endgroup$ Commented Sep 22, 2019 at 18:46
  • $\begingroup$ BBGKY hierarchy might obstruct from having an exact closed equation for 1P-RDM. $\endgroup$
    – Sunyam
    Commented Sep 24, 2019 at 21:39

1 Answer 1


Switching to a different formulation can make life easier. To make life easier, I'll do two things:

  • I'll use non-relativistic quantum field theory (NRQFT), also called second quantization. This allows deriving the general continuity equation as a simple operator equation, without any need for multi-electron wavefunctions or reduced density matrices. The operator equation is derived as a direct consequence of the gauge invariance of the Hamiltonian.

  • When it's time to introduce a state, I'll use the full state. This is easier than trying to use an explicit reduced density matrix, and for the time-independent continuity equation, the 1P-RDM result falls out automatically anyway. That's because conceptually, when we work with a reduced density matrix, all we're really doing is agreeing to consider only a limited set of operators — such as the set of "one-particle" operators, which includes the electric current operator.

After deriving the desired result, I'll briefly explain how the derivation relates to the von Neumann equation.

The model, part 1: The Hilbert space and Hamiltonian

Let $\varphi(\mathbf{x})$ be a field operator parameterized by $\mathbf{x}$, which denotes the 3-tuple of coordinates for a point in space. Let $\varphi^\dagger(\mathbf{x})$ denote the adjoint of $\varphi(\mathbf{x})$. The OP specifies electrons, which are fermions, so these operators satisfy the anti-commutation relations $$ \big\{\varphi(\mathbf{x}),\,\varphi^\dagger(\mathbf{y})\big\} = \delta^3(\mathbf{x} - \mathbf{y}) \tag{1} $$ and $$ \big\{\varphi(\mathbf{x}),\,\varphi(\mathbf{y})\big\} = 0 \tag{2} $$ with $\{A,B\} := AB+BA$. This is QFT's version of saying that the wavefunction must be completely antisymmetric. As suggested in the OP,

  • I'll pretend that electrons have spin $0$, so the field operators don't need a spin index.

  • I'll treat the magnetic field as a time-independent external field, expressed in terms of a gauge field $\mathbf{A}(\mathbf{x})$.

For any given choice of the gauge field, the Hamiltonian is $$ H[\mathbf{A}] = H_2[\mathbf{A}] + H_4[\mathbf{A}] \tag{3} $$ \begin{align} H_2[\mathbf{A}] &= \int d^3x\ \frac{\big(\mathbf{D}\varphi(\mathbf{x})\big)^\dagger \big(\mathbf{D}\varphi(\mathbf{x})\big)}{2m} \tag{4} \\ H_4[\mathbf{A}] &= \int d^3x\,d^3y\ \big(\varphi(\mathbf{x})\varphi(\mathbf{y})\big)^\dagger v(\mathbf{x}-\mathbf{y}) \varphi(\mathbf{x})\varphi(\mathbf{y}) \tag{5} \end{align} where $$ \mathbf{D} = \nabla + ie\mathbf{A}(\mathbf{x}) \tag{6} $$ and where $v(\mathbf{x})$ is a real-valued, decreasing positive function of $|\mathbf{x}|$ that describes the repulsive Coulomb interaction between electrons. For any given gauge field $\mathbf{A}$, the Hamiltonian is manifestly non-negative, and a state $|0\rangle$ that satisfies $$ \varphi(\mathbf{x})|0\rangle = 0 \tag{7} $$ for all $\mathbf{x}$ is clearly a state with the lowest possible energy (zero). We can interpret this as the vacuum state, which has no electrons. Applying the operator $$ \varphi^\dagger(f) := \int d^3x\ f(\mathbf{x})\,\varphi^\dagger(\mathbf{x}) \tag{8} $$ to any state adds an electron with wavefunction $f$. The Pauli exclusion principle is enforced by the anti-commutation relations, which imply that applying $\varphi^\dagger(f)$ twice to any state gives zero. In words, no two electrons can have the same wavefunction.

The model, part 2: Gauge invariance and observables

We're really defining a whole family of models, one for each choice of the gauge field $\mathbf{A}$. In this context, an "observable" is really a whole family of operators, one operator $F[\mathbf{A}]$ for each choice of the gauge field $\mathbf{A}$, such that the following condition holds for all functions $\theta(\mathbf{x})$: $$ U[\theta]\,F[\mathbf{A}+\nabla\theta]\,U^{-1}[\theta] = F[\mathbf{A}] \tag{9} $$ where the unitary operator $U[\theta]$ is defined by the conditions \begin{gather} U[\theta]\,\varphi(\mathbf{x})\,U^{-1}[\theta] = \exp\big(ie\,\theta(\mathbf{x})\big)\,\varphi(\mathbf{x}). \tag{10} \end{gather} In words, equation (9) says that observables are gauge invariant. An example is $$ F[\mathbf{A}] = \varphi^\dagger(\mathbf{x}) \left(ie\int_\mathbf{y}^\mathbf{x} \mathbf{A}\cdot d\mathbf{z} \right) \varphi(\mathbf{y}) \tag{11} $$ for a fixed pair of points $\mathbf{x}$ and $\mathbf{y}$. This satisfies the condition (9), so it qualifies as an observable.

Deriving the general (operator) continuity equation

The Hamiltonian $H[\mathbf{A}]$, regarded as a functional of the gauge field, is gauge invariant in the sense defined above. The continuity equation (as an operator equation) can be derived from the gauge invariance of $H[\mathbf{A}]$. To do this, start with the identity $$ U[\theta]\,H[\mathbf{A}+\nabla\theta]\,U^{-1}[\theta] = H[\mathbf{A}], \tag{12} $$ which expresses the gauge invariance of $H[\mathbf{A}]$. Vary both sides with respect to the arbitrary function $\theta$ and then set $\theta=0$ to get $$ \left(\frac{\delta}{\delta\theta} H[\mathbf{A}+\nabla\theta]\right)_{\theta=0} + \left(\frac{\delta}{\delta\theta} U[\theta]H[\mathbf{A}]U^{-1}[\theta]\right)_{\theta=0} = 0. \tag{13} $$ Rewrite the first term using this identity: $$ \left(\frac{\delta}{\delta\theta} H[\mathbf{A}+\nabla\theta]\right)_{\theta=0} \propto \nabla\cdot\mathbf{J} \tag{14} $$ with $$ \mathbf{J}(\mathbf{x}) := \frac{\delta}{\delta \mathbf{A}(\mathbf{x})}H[\mathbf{A}]. \tag{15} $$ Equation (15) defines the current-density operator, and by working out the right-hand side, we see that it resembles the expression shown in the OP — except that here it is an operator, involving the electron field operators in place of a single-particle wavefunction. For the second term in (13), use the fact that any operator $X$ satisfies identity $$ \left(\frac{\delta}{\delta\theta(\mathbf{x})} U[\theta]XU^{-1}[\theta] \right)_{\theta=0} \propto [X,J_0(\mathbf{x})] \tag{16} $$ where $J_0$ is the charge-density operator $$ J_0(\mathbf{x}) := e\varphi^\dagger(\mathbf{x})\varphi(\mathbf{x}). \tag{17} $$ Altogether, this shows that equation (12), which expresses the gauge invariance of the Hamiltonian, implies the continuity equation $$ \nabla\cdot\mathbf{J}(\mathbf{x}) \propto [H,J_0(\mathbf{x})] \tag{18} $$ for the current- and charge-density operators.

Deriving $\nabla\cdot\mathbf{j}=0$

The continuity equation derived above is an operator equation. It doesn't refer to any state. To answer the question in the OP, we need to introduce a multi-electron state and then consider the "one-particle reduced density matrix." Let $|\psi\rangle$ be an arbitrary $N$-electron state-vector: $$ |\psi\rangle = \int d^3x_1\cdots d^3x_N\ \Psi(\mathbf{x}_1,...,\mathbf{x}_N) \varphi^\dagger(\mathbf{x}_1)\cdots \varphi^\dagger(\mathbf{x}_N)|0\rangle \tag{19} $$ where $\Psi$ is the $N$-electron wavefunction. We might as well use the abbreviation $$ \psi(X) := \frac{\langle\psi|X|\psi\rangle}{\langle\psi|\psi\rangle} \tag{20} $$ for any operator $X$, because all predictions that we make with quantum theory ultimately boil down to expressions of this form. I'll follow mathematicians and refer to $\psi(\cdots)$ as the state. It takes any operator as input and returns its expecation value as output. Inserting the operator-valued continuity equation into $\psi(\cdots)$ gives $$ \nabla\cdot\psi(\mathbf{J})\propto\psi\big([H,J_0]\big). \tag{21} $$ If we further assume that $\psi(\cdots)$ is a stationary state — that is, if we assume that $|\psi\rangle$ is an eigenstate of $H$ — then the right-hand side is zero, which leaves $$ \nabla\cdot\mathbf{j} = 0 \hskip1cm\text{with}\hskip1cm\mathbf{j} := \psi(\mathbf{J}). \tag{22} $$ This is the desired result.

Relation to the 1P-RDM

Suppose that $X$ is any operator of the form $$ X = \int d^3x\,d^3y\ \varphi^\dagger(\mathbf{x})M(\mathbf{x}-\mathbf{y}) \varphi(\mathbf{y}) \tag{23} $$ where $M$ is an ordinary function or differential operator. The (components of the) current operator $\mathbf{J}$ is one example. For any such operator $X$, the quantity $\psi(X)$ is equal to the trace of $X$ multiplied by the "one-particle reduced density matrix" defined in the OP. This is clear by inspection of equations (19), (20), and (23), taking the anti-commutation relations (1)-(2) into account.

We don't need to make the reduced density matrix explicit, because whenever we talk about a reduced density matrix, all we're really doing is agreeing to consider only a limited set of operators as inputs to the state $\psi(\cdots)$. In the "one-particle" case, we're agreeing to consider only operators of the form (23), which are sometimes called "one-particle" operators because they annihilate-and-then-create just one particle. The operator $\mathbf{J}$ defined above has this form, so the result (22) is automatically a 1P-RDM equation.

Relation to the von Neumann equation

To relate the preceding derivation to the von Neumann equation, start with $$ \frac{d}{dt}\psi(\cdots) \propto\psi\big([H,\cdots]\big), \tag{24} $$ which may be interpreted as the time-evolution equation for the density matrix. (Remember that working with a reduced density matrix really just means that we've agreed to consider only a limited set of operators.) When we write it as in (24), we don't need to worry about how to define an effective one-particle Hamiltonian, which would be awkward because of the two-particle Coulomb interaction. All we need to do is insert the charge-density operator $J_0$ into (24) and use the identity (21) to get $$ \frac{d}{dt}\psi(J_0) \propto \nabla\cdot\mathbf{j}. \tag{25} $$ For a stationary state, the left-hand side is zero, which leaves the desired result (22).

  • 1
    $\begingroup$ Thank you very much for this detailed answer. This was exactly what I wanted to know. As this result seems to me of some importance to my field of study and research I might want to cite this in talks or papers. Have you got any objections or requirements on the citation? $\endgroup$ Commented Sep 26, 2019 at 8:41
  • 1
    $\begingroup$ @Rudi_Birnbaum That's right. The derivation of equation (18) still holds, and the new term in the Hamiltonian doesn't contribute to (15). The new term will make finding an explicit stationary state of the form (19) even more difficult, but given such a state, the rest of the equations hold. $\endgroup$ Commented Sep 29, 2019 at 13:50
  • 1
    $\begingroup$ @Rudi_Birnbaum Right. Something like $\int \varphi^\dagger(\mathbf{x})V(\mathbf{x})\varphi(\mathbf{x})$ with a prescribed function $V$. $\endgroup$ Commented Oct 4, 2019 at 16:26
  • 1
    $\begingroup$ @Rudi_Birnbaum Both terms on the left side of equation (13) come from the left side of (12). Applying $\delta/\delta\theta$ to the right side of (12) gives zero. $\endgroup$ Commented Oct 8, 2019 at 22:25
  • 1
    $\begingroup$ @Rudi_Birnbaum The derivative of $H[A+\nabla\theta]$ with respect to $A$ is $\propto J$. The derivative with respect to $\theta$ would be the same except for the $\nabla$. Inside the integral, we get $\delta\nabla\theta(x)/\delta\theta(y)=\nabla\delta(x-y)$, and then integrating-by-parts gives $\nabla$ applied to what we would have got from differentiating $H$ with respect to $A$, except for an overall sign. That is, $\propto \nabla\cdot J$. $\endgroup$ Commented Oct 10, 2019 at 0:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.