QM Continuity Equation: many-electron in the magnetic field version? 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.
 A: 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).
