Summary: This turned out to be a rather trivial one indeed. As Marek mentioned in the comment, the continuity equation is trivial. And it indeed turns out be so. Godfrey Miller elaborates on this, showing that the continuity equation merely shows that the transition amplitude is constant with time. All this is only confusing as long as you stick to the belief that the continuity equation is some exclusive indication of probability density.
$\rho = \psi_1^* \psi_2$ is the probability density for transition from state 2 to 1. I am having trouble interpreting the probability though $$dP_{1\leftarrow 2}=\psi_1^* \psi_2 d^3 r $$ does this mean the probability in a region of volume $d^3 r$ to go from state 2 to state 1? That does not make sense at all!
EDIT: I more or less agree with Marek's answer. However, that makes it seem this is not the probability density. How do we reconcile that with the fact that it satisfies the continuity equation. If we have two states:
$$ -i \hbar \partial_t \psi_{1,2}= -\frac{\hbar^2}{2m} \nabla^2 \psi_{1,2} + V\psi_{1,2}$$$$ i \hbar \partial_t \psi_{1,2}= -\frac{\hbar^2}{2m} \nabla^2 \psi_{1,2} + V\psi_{1,2}$$
After making standard manipulations (taking conjugate, multiplying by conjugate, subtracting) we would get $$\partial_t ( \psi_1^*\psi_2) +\frac{\hbar}{2mi} \nabla\cdot \left( \psi_1^* \nabla \psi_2 - (\nabla \psi_1^* )\psi_2\right)=0$$ This seems consistent with the continuity equation: $$\partial_t \rho +\nabla.\mathbb{j} =0$$ However, we seem to be changing its relation to the probability by making the probability $P=\int \int\rho^*(r,r') \rho(r,r') d^3r d^3r' $ instead of the regular $P=\int \rho(r) d^3r$ For a kind of intuitive "dimensional consistency" can we instead say $P=\int \int \; \sqrt{\rho^* \rho} \; d^3r d^3r'$ $$****$$ Full steps: $\psi_1$ and $\psi_2$ are two states of the same system, i.e both satisfy the SE of the same potential: Taking the SE for $\psi_2$ and multiplying by $\psi_1^*$ : $$ \psi_1^*\times\left( i \hbar \partial_t \psi_2 = -\frac{\hbar^2}{2m} \nabla^2 \psi_2 + V\psi_2\right)$$ Taking the SE for $\psi_1$ complex conjugating the entire equation and multiplying by $\psi_2$ : $$\left( -i \hbar \partial_t \psi_1^*= -\frac{\hbar^2}{2m} \nabla^2 \psi_1^* + V\psi_1^*\right)\times \psi_2$$
($V$ is real, no decay processes etc) On subtracting the second equation from the first
$$\partial_t (\psi_1^*\psi_2 )= -\frac{\hbar}{2mi} \left( \psi_1^* \nabla^2 \psi_2 - (\nabla^2 \psi_1^* )\psi_2\right)$$ which is the same as $$ \partial_t (\psi_1^*\psi_2 )= -\frac{\hbar}{2mi} \nabla\cdot \left( \psi_1^* \nabla \psi_2 - (\nabla \psi_1^* )\psi_2\right)$$