How to understand vector $J$ as probability current from Schroedinger equation? Multiplying the Schroedinger equation by $\psi^*$:
\begin{align}
\psi^*[i\hbar\frac{\partial\psi}{\partial t}] = \psi^*[-\frac{\hbar^2}{2m}\nabla^2 + V]\psi
\end{align}
Its complex conjugate by $\psi$: 
\begin{align}
\psi[-i\hbar\frac{\partial\psi^{\ast}}{\partial t}]= \psi[-\frac{\hbar^2}{2m}\nabla^2 + V]\psi^*
\end{align}
and subtracting: 
\begin{align}i\hbar[\psi^*\frac{\partial\psi}{\partial t} + \psi\frac{\partial\psi^{\ast}}{\partial t}] = \psi^*(-\frac{\hbar^2}{2m}\nabla^2\psi) + \psi^*V\psi\; -\; \psi(-\frac{\hbar^2}{2m}\nabla^2\psi^*) - \psi V\psi^*\end{align}
The potential terms cancel, and dividing by $\hbar$ and multiplying by $-i$:
\begin{align}[\psi^*\frac{\partial\psi}{\partial t} + \psi\frac{\partial\psi^{\ast}}{\partial t}] = i\frac{\hbar}{2m}[\psi^*\nabla^2\psi\; -\; \psi\nabla^2\psi^*]\end{align}
The left side is the derivative of the product $\psi^*\psi$, what is by definition the density of probability $\rho$:  
\begin{align}\frac{\partial\rho}{\partial t} = \frac{\partial}{\partial t}[\psi^*\psi ] = i\frac{\hbar}{2m}[\psi^*\nabla^2\psi  - \psi\nabla^2\psi^*]\end{align} 
But now things are not so clear: creating the vectors $\nabla\psi$  and $\nabla\psi^*$, the right side can be expressed as the divergence of a vector:
\begin{align}\nabla\cdot[\psi^*\nabla\psi  - \psi\nabla\psi^*] = \nabla\psi^*\cdot\nabla\psi + \psi^*\nabla^2\psi - \nabla\psi\cdot\nabla\psi^* -  \psi\nabla^2\psi^* =  \psi^*\nabla^2\psi -  \psi\nabla^2\psi^*\end{align}
Of course if we call $-i\hbar/(2m)[\psi^*\nabla\psi  - \psi\nabla\psi^*] = \mathbf{J}$, the expression becomes similar to the continuity equation:  
\begin{align}\frac{\partial\rho}{\partial t} = i\frac{\hbar}{2m}\nabla\cdot[\psi^*\nabla\psi  - \psi\nabla\psi^*] => \frac{\partial\rho}{\partial t} + \nabla\cdot\mathbf{J} = 0\end{align}
I understand the idea: if the probability decreases in any region, it must flow outside, so that the total probability is always one. But how such a flow relates to the vector difference in the brackets? 
I mean: the left side is a known subject, the square of the wave function is associated to the probability, but it is not the case of the right side. 
 A: Let $\psi = |\psi|e^{\mathfrak{j} \alpha}$, then 
$\ln \left( \frac{\psi}{\psi^*} \right)=2\mathfrak{j}\alpha $ and $$
\\ 2\mathfrak{j} \nabla \alpha  = \nabla \ln \left( \frac{\psi}{\psi*} \right)=\nabla \ln\psi - \nabla \ln\psi^* \\=\frac{1}{\psi}\nabla \psi - \frac{1}{\psi^*}\nabla \psi^*.$$ 
Now multiply both sides with $\rho=\psi\psi^*$
$$\psi\psi^*2\mathfrak{j} \nabla \alpha = \psi\psi^* \left(\frac{1}{\psi}\nabla \psi - \frac{1}{\psi^*}\nabla \psi^* \right)\\
= \psi^* \nabla \psi -\psi \nabla\psi^*.$$ 
But since $-\mathfrak{j}\hbar/(2m)[\psi^*\nabla\psi  - \psi\nabla\psi^*] = \mathbf{J}$ we also have $$\mathbf{J} = -\mathfrak{j}\frac{\hbar}{2m}\psi\psi^*2\mathfrak{j} \nabla \alpha \\= \rho \frac{\hbar}{m}\nabla\alpha = \rho \mathbf{v}$$
Here I defined a "velocity" $\mathbf{v} = \frac{\hbar}{m} \nabla \alpha$ with which the probability density $\rho$ is "convected". Notice that $\mathbf{v}$ is proportional with the gradient of the equiphase surfaces $\alpha$, in other words the velocity field comprise the orthogonal trajectories, rays, of these surfaces.
