How do we choose the standard probability current?

Well, here is a property it satisfies. I'm sure the Bohm-Aharonov advection mavens have neat geometric names about such things, but I'm drawing a blank for the moment. Irrotational flow?

Nevertheless, peremptorily normalizing, for simplicity, $$ \mathbf{J}\equiv \psi^* \nabla \psi - \psi \nabla \psi^*= - \nabla \rho + 2 \psi^* \nabla \psi = \nabla \rho - 2 \psi \nabla \psi ^*, $$ you see that $$ \rho ~\nabla \times \mathbf{J}=2\rho (\nabla \psi^* \times \nabla \psi) = \nabla \rho \times \mathbf{J}~, $$ so $$ \nabla \times \mathbf{J} -\frac{\nabla \rho}{\rho} \times \mathbf{J}=0=\rho~~\nabla\times (\mathbf{J}/\rho)~. $$ At a minimum, this might exclude your 2nd and 3rd options, $\mathbf{J}_2,~\mathbf{J}_3$.

Current over density might serve to define some sort of effective velocity of probability flow (Dirac, Landau-Lifschitz), so then with vanishing vorticity.

$\mathbf{J} /\rho$ being irrotational, it may be thought of as a potential flow (a gradient of the phase of the wavefunction), which goes under the name of the Madelung quantum Euler equations' formulation.

Your first option $\mathbf{J}_1$ is dimensionally inconsistent, anyway, so it needs a dimensionfull constant in front of the extra curl term. But there is no good reason the power of the density needs to be one, as you have chosen. If you considered $\nabla \times (\rho^n \mathbf{J})$ instead, then you see the irrotational condition automatically projects out the n = -1 case, but not all other ones.

• (One may always reinstate the normalization $-i\hbar/2m$ in the current, but why?) For broader contexts, see WP.

A simplest exploration is afforded by a stationary system, such as an atom, where the current is divergenceless, $ \nabla \cdot \mathbf{J}=0$, as the only complex part of the wavefunction is the azimuthal exponential, so only the φ component of it survives in polar coordinates: uniform rotational flow for a stationary system!

Firstly, you can derive the continuity equation from Schroedinger's equation ($\hbar = 1$) \begin{equation} -\frac{1}{2m} \nabla^2 \psi + V(x) \psi = i \partial_t \psi \end{equation} by multiplying it by $\overline{\psi}$ and its complex conjugate by $\psi$ and summing them. With a little manipulation, you get:

\begin{equation} \partial_t (\psi \overline{\psi}) -\frac{i}{2m} \nabla \cdot \left( \overline{\psi} \nabla \psi - \psi \nabla \overline{\psi} \right) = 0. \end{equation}

Since you want local conservation of probability density, you define the current $\vec{j}$ as

\begin{equation} \vec{j}_0 = \frac{-i}{2m} \left( \overline{\psi} \nabla \psi - \psi \nabla \overline{\psi} \right) \end{equation}

So there's nothing "special" about it, it just imposes local conservation of probability density, and there is no physical reason to add another term to it, since the current arises naturally from Schroedinger's equation.

Now, the above current has to be redefined depending on the system you are working on. Supposing there is an external magnetic field ($\vec{B} = \nabla \times \vec{A}$), deriving the current in the same way as before you will get another term in $\vec{j}$, proportional to $ \overline{\psi} \vec{A}\psi$.

Now, as said, there is no physical reason to add another term, but there is a special case. In a system where you have an external magnetic field AND you consider its interaction with spin, you define your current with an extra term, proportional to $\nabla \times (\overline{\psi} \vec{S} \psi)$, even though it does not naturally arise try to derive $\vec{j}$ as before (note that $\nabla \cdot \nabla \times \vec{f} = 0$). This is physically justified. In non-relativistic quantum mechanics spin does not arise naturally, but this term exists in the conserved current of the Dirac equation, and survives in the non-relativistic limit. So this is imposed and you get \begin{equation} \vec{j}= \vec{j}_0 - \frac{q}{m} \vec{A} \overline{\psi} \psi\ + \frac{\mu_s}{s} \nabla \times (\overline{\psi} \vec{S} \psi) \end{equation}

P.S.1: Also, remember that Quantum Mechanics is invariant under a global phase transformation $e^{i\alpha}$, and $\vec{j}_0$ is the current associated to this kind of symmetry.

P.S.2: Maybe there are some constants missing, if there are, I apologize.