Rationalizing the current operator In eq. 1.191 of Mahan Many-Particle Physics, the current operator is defined as
$$j(r)=\frac{1}{2}\sum_ie_i[v_i\delta(r-r_i)+\delta(r-r_i)v_i].$$
I'm trying to make some sense of this definition. I understand that this can be derived from a continuity equation for quantum probability density, however, I don't have an intuition for what this form means, i.e. what are the delta functions doing.
To make my point more transparent, let's consider the space of one-particle wave functions. Then,
$$j(r)=-\frac{i\hbar}{2}e[\partial_{r_0}\delta(r-r_0)+\delta(r-r_0)\partial_{r_0}].$$
To show the depth of my confusion, let's consider the following manipulations. If we apply this operator to an arbitrary wave function, we have
$$j(r)\psi(r_0)=-\frac{i\hbar}{2}e[\partial_{r_0}\delta(r-r_0)+\delta(r-r_0)\partial_{r_0}]\psi(r_0)$$
$$j(r)\psi(r_0)=-\frac{i\hbar}{2}e[(\partial_{r_0}\delta(r-r_0))+2\delta(r-r_0)\partial_{r_0}]\psi(r_0)$$
$$j(r)\psi(r_0)=-\frac{i\hbar}{2}e[-\delta'(r-r_0)+2\delta(r-r_0)\partial_{r_0}]\psi(r_0)$$
In this form it's no more clear what the role of the delta function is, and I am extra confused by the appearance of the derivative of a delta function.
Can anyone offer an explanation of how to intuitively understand the probability current expression?
 A: The important point is that we need to distinguish positions of particles $\mathbf r_k$ from the position in space $\mathbf x$ where the current density is sought. The delta function that appears in these calculations is
$$
\delta^{(3)}(\mathbf x-\mathbf r_k)
$$
because the charge or current density at $\mathbf x$ of point particle located at $\mathbf r_k$ (electrons are assumed to be point particles in non-relativistic qm) is given by delta distribution (all charge is at a single point $\mathbf r_k$).
In a classical EM theory that would be it for the definition, current density of a point charged particle is simply a singular delta distribution.
However, in quantum theory, we do not have values of particle actual coordinates, instead, there is this function $\psi(\mathbf r_1,...\mathbf r_N,t)$ for all possible coordinates which is used to assign probabilities to all imaginable configurations $\{\mathbf r_k\}$ of $N$ particles, say in an atom, a molecule or a crystal. It is also used to calculate expected averages of various quantities that in classical physics are functions of positions and momenta of the particles. This is the case with current density.
The current density operator is supposed to give us expected net current density at some point of physical space $\mathbf x$, based on the psi function. The supposed usage is
$$
\mathbf j_{exp}(\mathbf x,t) = \langle \psi |\hat{\mathbf j}(\mathbf x,\mathbf r_1,...\mathbf r_N,t)|\psi\rangle = \int \psi^* \hat{\mathbf j} \psi \, d^3\mathbf r_1 ... d^3\mathbf r_N
$$
where the current density operator is
$$
\hat{\mathbf j}  = \sum_{k=1}^N \frac{q_k }{2m}\left(\boldsymbol \pi_k \delta(\mathbf x - \mathbf r_k)+\delta(\mathbf x - \mathbf r_k)\boldsymbol \pi_k\right)
$$
where $\boldsymbol \pi_k$ is operator of kinetic momentum of $k-$th particle.
Because of the integration over all possible positions of the particles, and because of the smooth nature of $\psi$ function (in most cases), this shoul generate a smooth function of $\mathbf j_{exp}(\mathbf x)$ with not trace of a singular delta-like distribution.
So the psi function allows us to perfrom a sort of smearing effect that removes the point concentration of the deltas in the current operator and this produces a smooth function of position $\mathbf x$.
If I simplify, the delta function for particle (say) $k=3$ does single thing: it fixes the $k-$th argument of $\psi$ to be equal to $\mathbf x$ and then the averaging obeys this restriction. All terms of the result correspond to a situation where at least one argument in $\psi$ is fixed to $\mathbf x$ and then weight of contribution of this case is calculated from the psi function, allowing the other particle coordinates to have all possible values.
A: The form of the current operator can be well understood if you are on top of the bra-ket notation. We can start from the single particle Schrodinger equation acting on the vector $\mid\psi\rangle$ with x-representation, one can obtain the current by performing some calculations
$$j(x)=\frac{e}{m}\text{Re}\left[\psi^*(x)\hat{p}\psi(x) \right]=\frac{e}{2m}\left(\langle\psi|x\rangle\langle x|\hat{p}|\psi\rangle+\langle \psi |\hat{p}^{\dagger}|x\rangle\langle x|\psi\rangle\right)$$
where $j(x)$ is the current acting on the state $|\psi\rangle$, name $j(x)=\langle\psi| \hat{j}(x)|\psi\rangle$, here $\hat{j}(x)$ is the single particle current operator expressing as
$$\hat{j}(x)=\frac{e}{2m}\left(|x\rangle\langle x|\hat{p}+\hat{p}^{\dagger}|x\rangle\langle x|\right)$$

*

*How to get the $\delta$-function type expression?

Let's consider the many body system, the current operator acting on the i-th particle state $\psi(x_i)$ is given as $\langle x_i|\hat{j}(x)|\psi\rangle$, which is given as
$$\begin{align}
\langle x_i|\hat{j}(x)|\psi\rangle&=\frac{e}{2m}\left(\langle x_i|x\rangle\langle x|\hat{p}|\psi\rangle+\langle x_i|\hat{p}^{\dagger}|x\rangle\langle x|\psi\rangle\right)\\&=\frac{e}{2m}\left(\delta(x_i-x)\frac{\hbar}{i}\frac{\partial\psi(x)}{\partial x}+\frac{\hbar}{i}\frac{\partial\delta(x_i-x)}{\partial(x_i-x)}\psi(x)\right)\\&=J(x)\psi(x)
\end{align}$$
where $J(x)$ is the current operator
$$J(x)=\frac{e\hbar}{2mi}\left(\delta(x_i-x)\frac{\partial}{\partial x}+\frac{\partial}{\partial x_i}\delta(x_i-x)\right)$$
ADDENDUM
Now considering the system with applied external electric field, by employing the form of the current operator, one can proof that the Hamiltonian can be expressed as $H=H_0+H_{ext}$. where $H_0$ is the unperturbed Hamiltonian, $H_{ext}$ is the Hamiltonian with external field, which is given as
$$\hat{H}_{ext}=-\int_x \hat{j}(x)\cdot \hat{A}(x)$$
where $\hat{A}(x)$ is the operator of the external field, as a function of position $\vec{x}=x$, （for simplicity consider the 1-D single particle case).
$$\hat{H}=\frac{e}{2m}\left(\hat{p}-e\hat{A}(x)\right)^2+\hat{V}(x)$$
is the total Hamiltonian (choose $c=1$), which gives the current operator with the external field
$$\begin{align}
\hat{j}(x)&=\frac{1}{2m}\left(|x\rangle\langle x|\hat{\pi}+\hat{\pi}^{\dagger}|x\rangle\langle x|\right)\\&=\hat{j}_0(x)+\hat{j}_A(x)
\end{align}$$
where $\hat{\pi}=\hat{p}-e\hat{A}(x)$. Here we just consider the first order of $A$,
$$\hat{H}_{ext}=-\frac{e}{2m}(\hat{p}\cdot\hat{A}(x)+\hat{A}(x)\cdot\hat{p})$$
Remember that all the operators we care about in this case behave like $n\times n$ matrices ($n\rightarrow +\infty $), to avoid confusions, one needs to calculate the term like $\langle x_i|-\int_x\hat{j}(x)\cdot\hat{A}(x)|\psi\rangle$ and proof that equals to $\langle x_i|H_{ext}|\psi\rangle$, which is straight by using the expansion of $\hat{j}_0(x)$:
$$\begin{align}-\int_x\langle x_i|\hat{A}(x)\cdot\hat{j}_0(x)|\psi\rangle&=\int_x\int_y\langle x_i|\hat{A}(x)|y\rangle\langle y|\hat{j}_0(x)|\psi\rangle\\&=-\frac{e}{2m}\int_x\int_y\left(A(y)\delta(x_i-y)\delta(y-x)\frac{1}{i}\frac{\partial\psi(x)}{\partial x}\right)\\&-\frac{e}{2m}\int_x\int_y\left(A(y)\delta(x_i-y)\frac{1}{i}\frac{\partial \delta(y-x)}{\partial(y-x)}\psi(x)\right)\\&=-\frac{e}{2m}\left(A(x)\frac{1}{i}\frac{\partial\psi(x_i)}{\partial x_i}+\frac{1}{i}\frac{\partial A(x_i)\psi(x_i)}{\partial x_i}\right)\\&=\langle x_i|\hat{H}_{ext}|\psi\rangle
\end{align}$$
