Uniqueness of the probability function for the Schrödinger equation

Question

David Bohm in Section (4.5) of his wonderful monograph Quantum Theory after defining the usual density probability function $P(x,t)=\psi^{*} \psi$ for the Schrödinger equation for the free particle in one dimension: \begin{equation} i \hbar \frac{\partial \psi}{\partial t}= - \frac{\hbar^2}{2m} \frac{\partial^2 \psi}{\partial x^2}, \end{equation} states that $P(x,t)$ is the unique function of $\psi(x,t)$ and the partial derivatives of $\psi$ with respect to $x$ all computed in $(x,t)$ which satisfies the following properties:

1. P is never negative;

2. the probability is large when $|\psi|$ is large and small when $|\psi|$ is small;

3. the significance of $P$ does not depend in a critical way on any quantity which is known on general physical grounds to be irrelevant: in particular this implies (since we are dealing with a nonrelativistic theory) that $P$ must not depend on where the zero of energy is chosen;

4. $\int P(x,t) dx$ is conserved over time, so that by eventually normalizing $P$ we can choose $\int P(x,t) dx=1$ for all $t$.

Bohm gives no mathematical argument at all and actually the statement seems completely unjustified to me.

Does someone know some reason why it should be true?

NOTE (1). Since the Schrödinger equation is a first-order equation, the time evolution of $\psi$ is fixed given the initial state $\psi(x,0)$: this is the reason why we require that the probability $P(x,t)$ depends on the state at time $t$, that is $\psi(x,t)$, and its spatial derivatives. To be explicit, the requirement that $P(x,t)$ is a function only of $\psi(x,t)$ and the partial derivatives of $\psi$ with respect to $x$ all computed in $(x,t)$ means that there exists a function $p$ such that $P(x,t)=p\left(\psi(x,t),\frac{ \partial \psi}{\partial x}(x,t),...,\frac{\partial^m \psi}{\partial x^m}(x,t)\right)$.

NOTE (2). The one given above is not a mathematically rigorous formulation of the problem, but the one originally given by Bohm. So we can feel free to attach a rigorous mathematical meaning to the different properties. In particular, as for property (iv) we can formulate it in a different (and not equivalent) mathematical form, by requiring that a local conservation law holds, in the sense that there exists a function $\mathbf{j}$, such that, if we put $\mathbf{J}(x,t)=\mathbf{j}\left(\psi(x,t),\frac{ \partial \psi}{\partial x}(x,t),...,\frac{\partial^m \psi}{\partial x^m}(x,t)\right)$, we get \begin{equation} \frac{\partial P}{\partial t} + \nabla \cdot \mathbf{J} = 0. \end{equation}

NOTE (3). Similar questions are raised in the posts Nonexistence of a Probability for Real Wave Equations and Nonexistence of a Probability for the Klein-Gordon Equation. Presumably Bohm had in mind the same kind of mathematical argument to tackle these three problems, so the real question is: what mathematical tool did he envisage to use? Maybe some concept from classical field theory or the theory of partial differential equations?

Answer 1

Bohm's assumptions are not mathematically precise, so you must attach a mathematical interpretation to them (especially statements $2$ and $3$). Since you have not done so yourself, I will try to interpret them in a way that I find reasonable.

Definition for $P$: We will require that the probability density $P_\psi(x,t)$ of any smooth function $\psi$ be a local function of its partial derivatives at $(x,t)$.

More formally, let $\psi(x,t)$ be a smooth function and let $j^\infty_{x,t}\psi$ denote the infinite jet prolongation of $\psi$ at $(x,t)$, i.e., its formal Taylor series expansion about $(x,t)$. Then we can write $$P_\psi(x,t) = p(j_{x,t}^\infty\psi),$$ for some function $p$ defined on the jet bundle. In terms of regularity, we will require $p$ to be continuous.

This is essentially the definition you proposed for $P$. In fact, this is already problematic since we will necessarily have to work with wavefunctions which are not smooth. Therefore it is already problematic to require that $P$ depend on higher derivatives since they are not guaranteed to even exist. Nevertheless, we will allow arbitrary dependence on higher partials and apply the consistency conditions for $P$ evaluated on smooth functions. We can then recover $P$ uniquely for arbitrary $L^2$ functions by continuity.

Assumption 1: The function $p$ is a real-valued and non-negative.

Assumption 2: The probability density $P_\psi(x,t)$ is a non-decreasing function of $|\psi(x,t)|$, i.e., $$P_{\psi_1}(x,t) \ge P_{\psi_2}(x,t) \iff |\psi_1(x,t)| \ge |\psi_2(x,t)|.$$

Assumption 3: The function $p$ is invariant under global phase, i.e., $$p(e^{i\theta} j^\infty_{x,t}\psi) = p(j^\infty_{x,t}\psi).$$

Assumption 4: If $\psi(x,t)$ is a normalized function, then $P_\psi(x,t)$ is likewise a normalized function.

Let us go through the assumptions one by one. Assumption $1$ is relatively straightforward as probability densities must be real and non-negative.

Property $2$ is, in my view, the most difficult to properly interpret. The way I have interpreted the property in Assumption $2$ is to say that the magnitude of the probability density at a point is directly reflective of the magnitude of the wavefunction at that point. This is what I feel to be the most direct transcription of Bohm's second property.

This assumption is in fact extremely strong, and it necessarily implies that $p$ is independent of all derivatives of $\psi$. This is essentially because the value of a smooth function and all of its derivatives can be independent prescribed at any point. This was already pointed out by @Kostas.

Lemma: Suppose that $p(j^\infty_{x,t}\psi)$ is a continuous non-decreasing function of $|\psi(x,t)|$. Then $p$ is independent of all derivatives of $\psi$, i.e., $$p(j^\infty_{x,t}\psi) = p(j^0_{x,t}\psi) = p(\psi(x,t)).$$

Proof: By Borel's theorem, given any complex sequence $(a_{n,m})_{n,m=0}^\infty,$ and any point $(x,t)$, there exists a smooth function $\psi$ such that $$\frac{\partial^{n+m}}{\partial x^n \partial t^m}\psi(x,t) = a_{n,m}.$$ Therefore we are able vary the individual entries of the Taylor series completely independently.

Suppose that $p$ is not constant on some partial $\partial_i \psi$. Then by Borel's theorem we can find smooth functions $\psi_1$ and $\psi_2$ such that all Taylor coefficients of $\psi_1$ and $\psi_2$ agree at $(x,t)$ except for $\partial_i$. Then $$p(j_{x,t}^\infty\psi_1) \neq p(j_{x,t}^\infty\psi_2),$$ and without loss of generality we may assume that $$p(j_{x,t}^\infty\psi_1) > p(j_{x,t}^\infty\psi_2).$$ Next, we can find some other smooth function $\psi_3$ which agrees with $\psi_2$ for all Taylor coefficients at $(x,t)$ except for the constant term, $\psi_3(x,t)\neq \psi_2(x,t)$. By continuity, we can choose $\psi_3(x,t)$ slightly larger than $\psi_2(x,t)$ but still such that $$p(j_{x,t}^\infty\psi_1) > p(j_{x,t}^\infty\psi_3).$$ Therefore we have $$|\psi_1(x,t)| = |\psi_2(x,t)| < |\psi_3(x,t)|,\ \ \ \ \text{and}\ \ \ \ \ p(j_{x,t}^\infty\psi_1) > p(j_{x,t}^\infty\psi_3),$$ in contradiction to the monotonicity assumption. $\square$

Therefore we will assume that $P_\psi(x,t) = p(\psi(x,t))$ from now on. Note that at this point, we no longer need the assumption that $\psi$ is smooth.

Assumption $3$ is also essentially just a direct transcription of Property 3. If we change the energy, we effectively change the wavefunction by a global phase factor $e^{iEt}$. Since we can always make the wavefunction real and positive at any given point $(x,t)$ by an appropriate choice of a global phase factor, it follows that $p$ is independent of the phase of $\psi(x,t)$, i.e., $$p(\psi(x,t)) = p(|\psi(x,t)|).$$

Note that this assumption is actually completely unnecessary. We could've deduced the above equation by a slightly modification of our lemma using Assumption $2$. I will keep it just for the sake of completeness.

Finally, we come to Assumption $4$. Bohm's statement for his Property $4$ is that the probabilities should be normalized at all times, namely for all $t$ we should have $$1 = \int_\mathbb{R} P(x,t)\ dx.$$

This has certain ambiguities however. Which time-evolution should we use? Naively, any self-adjoint operator $H$ with a spectrum which is bounded below (so that there is a lowest energy level) should be able to serve as a valid Hamiltonian. If we require that the assignment $\psi \mapsto P_\psi$ be be universally valid, i.e., Hamiltonian independent, then we must require that $P(x,t)$ be normalized with respect to the unitary evolution generated by any Hamiltonian.

It can be shown that given any unitary $U$, there exists some admissible (self-adjoint, bounded below) Hamiltonian $H$ such that $U = e^{iH}$. In fact, we do not even need to consider the set of all admissible Hamiltonians, but rather just the set of all bounded Hamiltonians due to the following theorem.

Theorem: Let $\mathcal{H}$ be a Hilbert space and let $U$ be any unitary operator on $\mathcal{H}$. Then there exists a bounded self-adjoint operator $A$ (with norm at most $\pi$) such that $$U = e^{iA}.$$

Proof: This is a simple consequence of the Borel functional calculus for bounded operators applied to the principal branch of the logarithm. See here for a complete proof. $\square$

Now, let $\psi_1(x)$ be some normalized wavefunction. Let us assume without loss of generality that $P$ is normalized so that $$1 = \int_\mathbb{R} P_{\psi_1}(x)\ dx.$$ Let $\psi_2(x)$ be some other arbitrary normalized wavefunction. Let $U$ be any unitary such that $\psi_2 = U\psi_1$. Then there exists some bounded Hamiltonian such that the time-evolution brings the initial state $\psi(x,t=0) = \psi_1(x)$ to $\psi(x,t=1) = \psi_2(x)$. This means that we must have $$1 = \int_\mathbb{R} P_{\psi_1}(x)\ dx = \int_{\mathbb{R}} P_{\psi_2}(x)\ dx.$$ Since $\psi_2$ was an arbitrary normalized function, it follows that $P_\psi$ is normalized for all normalized $\psi$. We take this as our Assumption $4$.

Physically, this assumption is essentially saying that we should be able to vary the potential of the Hamiltonian so as to drive any normalized wavefunction arbitrarily close to any other normalized wavefunction. Since $P$ is conserved under this evolution, it must be normalized given any normalized wavefunction.

Note that this implies that we must have $p(0) = 0$. Otherwise $p(0) > 0$ will give a divergent integral for any normalized compactly supported function $\psi$.

Now let $y>0$. Define $\psi_y(x,t)$ to be equal to $1/\sqrt{y}$ for $x\in (0,y)$ and zero elsewhere. Then we have $$\int_\mathbb{R} |\psi_{y}(x,t)|^2\ dx = 1 = \int_\mathbb{R} p(|\psi_y(x,t)|)\ dx = \int_0^y p(1/\sqrt{y})\ dx = yp(1/\sqrt{y}).$$ Therefore we must have $$p(|\psi(x,t)|) = |\psi(x,t)|^2.$$

This is the desired claim. Of course, you might disagree with how I've interpreted some of Bohm's statements. But as you've said yourself in the question, some rigorous definitions must be assigned to these physical properties. These are simply what I felt to be the most faithful.

Answer 2

(i) Clearly $P$ is non-negative since $a^*a=|a|^2$ is non-negative for all complex numbers $a$.

(ii) The probability to find the particle in an infinitesimal interval between $x$ and $x+dx$ is given by $Pdx=|\psi(x,t)|^2dx$ where this probability clearly is large for large $|\psi|$ and vice versa.

(iii)Indeed it can be proven that shifting the zero of energy by a constant(i.e., $V(x) \to V(x)+V_o$, where $V_o$ is constant) changes the overall phase of the wavefuntion so that $\psi(x,t) \to \psi(x,t)\exp(-iV_ot/\hbar)$, but this does not affect $P$, since $P_{\text{new}}=\psi_{\text{new}}^*\psi_{\text{new}}=[\psi^*\exp(+iV_ot/\hbar)][\psi\exp(-iV_ot/\hbar)]=\psi^*\psi=P$.

(iv) There're a lot of ways to prove this; One way is to show $\dfrac{d}{dt}{\displaystyle \int_{-\infty}^{+\infty}\psi^*\psi dx }=0$, You can go about doing this by pulling the derivative inside the integral and applying the product rule on $\dfrac{d[\psi^*\psi]}{dt}$ and using Schrodinger equation and integration by parts to prove the desired result.

Answer 3

OP is asking for a mathematical proof of uniqueness of some or all of the properties of P stated by Bohm. It is easy to see that the i), ii), and iii) can be violated separately and even together. Here is a trivial example: $P = (Ψ^* Ψ + |\partial Ψ|^2)^2$ It doesnt satisfy iv) of course. Prooving that no such expression satisifes iv) is hard. One way to do that is to assume P is a polynomial function of Ψ and its derivatives (it is in the example) and proove that all coefficients of the polynomial vanish except for the coefficient of the term with $Ψ^* Ψ$. It is a silly business, but it may help you if you are trying to understand why the Schroedinger equation is the only possible equation!

Answer 4

After a long search in the literature, I must admit that the issue of establishing the uniqueness of the probability density for the Schrödinger equation as discussed by Bohm has attracted very little interest. Actually the only work which explicitly deals with this issue is Uniqueness of Conserved Currents in Quantum Mechanics, even though the uniqueness problem had been previously discussed in the context of de Broglie-Bohm pilot wave theory (see the papers there quoted).

Holland uses relativistic considerations to show the uniqueness of the conserved current for the Klein-Gordon equation, from which he deduces in the non-relativistic limit the analogous uniqueness result for the Schrödinger equation.

Anyway, as the same author suggested to me in a private communication, a more direct proof could be given if we start directly from the Schrödinger equation in the case of a potential $V$: \begin{equation} i \hbar \frac{\partial \psi}{\partial t}(\mathbf{x},t) = -\frac{\hbar^2}{2m} \Delta \psi(\mathbf{x},t) + V(\mathbf{x}) \psi(\mathbf{x},t), \end{equation} and use the method he introduces in his paper to derive all the conserved currents $(P,\mathbf{J})$ which are functions only of $\psi$, the first-order partial derivatives of $\psi$ and eventually $V$, by studying their transformation properties under the Galilei transformations and under possible changes of $V$. Let us note that we have to be careful in applying this procedure since $\psi$ is not invariant under Galilei transformations, but it changes appropriately: see e.g. Commins, Quantum Mechanics or Galilei Invariance of the Schrodinger Equation.

This will give a general set of conserved currents, and we could investigate whether the additional assumption that $P$ must depend only on $\psi$ and not on its first-order partial derivatives (which is a particular consequence of property (ii) required by Bohm) eventually implies a unique expression for $P$ or not. In the last case, some other condition (maybe also derived by property (ii) which is quite vague) should be added in order to obtain the uniqueness.

Anyway, this approach would in any case suffer of a lack of generality, since you must assume from the very beginning that $P$ and $\mathbf{J}$ do not depend on partial derivatives of $\psi$ of order greater than one, an assumption which seems completely absent in Bohm's discussion, even though a very plausible assumption on physical grounds (see the remark made on this point by Holland in his work). For this reason, I am pretty sure that this approach was not the one Bohm had in mind when he made his uniqueness statement, even though I have no idea about the kind of argument he could have envisaged.

Answer 5

After some time, I have tackled again the issue of proving the uniqueness of the probability function for the Schrödinger equation along the lines suggested to me by Peter Holland (see my previous answer), and I discovered that it is actually simpler than I thought.

Consider the three-dimensional Schrödinger equation in the case of a potential $V$: \begin{equation} i \hbar \frac{\partial \psi}{\partial t}(\mathbf{x},t) = -\frac{\hbar^2}{2m} \Delta \psi(\mathbf{x},t) + V(\mathbf{x}) \psi(\mathbf{x},t). \end{equation} Since the Schrödinger equation is a first-order equation, taking as initial data $\mathbf{x} \mapsto \psi(\mathbf{x},0)$, the probability density must be a function of the local value of $\psi$ and of some of its spatial derivatives, in other terms we shall assume that the probability density $P$ is given by $P_\psi(\mathbf{x},t)=p(D_{\mathbf{x}}^{n} \psi(\mathbf{x},t))$, where $p$ is a smooth function, $n$ is a non negative integer, and $D_{\mathbf{x}}^{n}f$ denotes the set of all partial derivatives of the function $f$ with respect to $\mathbf{x}=(x_1,x_2,x_3)$ from order zero (that is the function itself) to order $n$. We shall make the following assumptions, which are the faithful mathematical formulations of Properties (1), (2), and (3) that EuYu gave in his answer:

Assumption 1: The function $p$ is a real-valued and non-negative.

Assumption 2: The probability density $P_\psi(\mathbf{x},t)$ is a non-decreasing function of $|\psi(\mathbf{x},t)|$, i.e., $$P_{\psi_1}(\mathbf{x},t) \ge P_{\psi_2}(\mathbf{x},t) \iff |\psi_1(\mathbf{x},t)| \ge |\psi_2(\mathbf{x},t)|.$$

Assumption 3: The function $p$ is invariant under global phase, i.e., for any $\theta \in \mathbb{R}$ we have $$p(e^{i\theta}D_{\mathbf{x}}^{n} \psi(\mathbf{x},t))= p(D_{\mathbf{x}}^{n} \psi(\mathbf{x},t)).$$

Let us remark that, as noted by EuYu in his answer, Assumption (3) follows from Property (3), since if we shift the potential $V(\mathbf{x})$ by adding a constant $V_0$, the new solution $\tilde{\psi}$ corresponding to the same initial data $\mathbf{x} \mapsto \psi(\mathbf{x},0)$ is $\tilde{\psi}(\mathbf{x},t)=\psi(\mathbf{x},t) \exp(-iV_0 t/ \hbar)$.

Now, the argument given by EuYu in his answer shows that these three assumptions imply that $p$ must depend only on$|\psi(\mathbf{x},t)|$, so that we shall write $P_\psi(\mathbf{x},t)=p(\overline{\psi}(\mathbf{x},t) \psi(\mathbf{x},t))$.

As for Property (4), we shall interpret it as suggested in the NOTE (2) to the post, by requiring that a continuity equation holds for some current density $\mathbf{J}=(J_1,J_2,J_3)$ which is a function of the local value of $\psi$ and of its first-order spatial derivatives. In other terms, if $U$ and $V$ are the real and imaginary part of $\psi$ respectively, we shall assume that the probability current density is given by $\mathbf{J}_\psi (\mathbf{x},t) = j(D_{\mathbf{x}}^1 U(\mathbf{x},t), D_{\mathbf{x}}^1 V(\mathbf{x},t))$, where $j(X_0,\dots,X_3,Y_0,\dots,Y_3)$ is a smooth function (here $X_0$ and $Y_0$ are the variables whose place are occupied respectively by $U$ and $V$, $X_i$ is the variable whose place is occupied by $\partial_{x_i} U$, $i=1,2,3$, and analogously $Y_i$ is the variable whose place is occupied by $\partial_{x_i} V$, $i=1,2,3$).

Assumption 4: The continuity equation

\begin{equation} \frac{\partial P}{\partial t} + \nabla \cdot \mathbf{J} = 0 \end{equation}

must hold for every possible choice of the potential $V$ and every possible choice of the initial data $\mathbf{x} \mapsto \psi(\mathbf{x},0)$.

The continuity equation satisfied by $(p,\mathbf{J})$ gives the constraint (here and in the following $\dot{p}$ denotes simply the derivative of the function $y \mapsto p(y)$): \begin{equation} (\partial_t \overline{\psi} \psi + \overline{\psi} \partial_t{\psi}) \dot{p}+ \sum_{l=1}^{3} \frac{\partial j_l}{\partial X_0} \partial_{x_l} U + \sum_{l=1}^{3} \frac{\partial j_l}{\partial Y_0} \partial_{x_l} V + \sum_{l=1}^{3} \sum_{k=1}^{3} \frac{\partial j_l}{\partial X_{k}} \partial_{x_l x_k} U + \sum_{l=1}^{3} \sum_{k=1}^{3} \frac{\partial j_l}{\partial Y_{k}} \partial_{x_l x_k} V = 0, \end{equation} which can be rewritten by using the Schrödinger equation as: \begin{equation} (I) \quad \frac{i \hbar}{2m} (U \Delta V - V \Delta U)\dot{p} + \sum_{l=1}^{3} \frac{\partial j_l}{\partial X_0} \partial_{x_l} U + \sum_{l=1}^{3} \frac{\partial j_l}{\partial Y_0} \partial_{x_l} V + \sum_{l=1}^{3} \sum_{k=1}^{3} \frac{\partial j_l}{\partial X_{k}} \partial_{x_l x_k} U + \sum_{l=1}^{3} \sum_{k=1}^{3} \frac{\partial j_l}{\partial Y_{k}} \partial_{x_l x_k} V = 0, \end{equation}

Now, we shall make use of the already recalled Borel's Lemma (which is a particular case of the remarkable Whitney's Extension Theorem) which says that for every sequence $(a_{n_1,\dots,n_k})_{n_1,\dots,n_k=0}^{\infty}$ in $\mathbb{R}$ and any point $\mathbf{x}_0 \in\mathbb{R}^k$, there exists a smooth function $f:\mathbb{R}^k \rightarrow \mathbb{R}$ such that \begin{equation} \partial_{x_1}^{n_1} \dots \partial_{x_k}^{n_k} f(\mathbf{x}_0)= a_{n_1,\dots,n_k}, \end{equation} in other terms the value of the function and of its partial derivatives in a given point can be arbitrarily prescribed. From this result and the arbitrariness of the initial data $\mathbf{x} \mapsto \psi(\mathbf{x},0)$, we see that the coefficients of the partial derivatives $\partial_{x_i}^2 U$ and $\partial_{x_i}^2 V$ in Equation (I) must vanish, so that we get: \begin{equation} \frac{\partial j_l}{\partial X_l} = -\frac{\hbar}{m} \dot{p}(X_0^2 + Y_0^2) Y_0 \quad (l=1,2,3), \end{equation} and \begin{equation} \frac{\partial j_l}{\partial Y_l} = \frac{\hbar}{m} \dot{p}(X_0^2 + Y_0^2) X_0 \quad (l=1,2,3), \end{equation}

and again by Borel's Lemma these relations imply that for $l=1,2,3$ we have for some functions $f_l$: \begin{equation} (II) \quad j_l= \frac{\hbar}{m}\dot{p}(X_0^2 + Y_0^2)(X_0 Y_l - Y_0 X_l) + f_l(X_0,X_j,X_k,Y_0,Y_j,Y_k), \end{equation} where $(l, j, k)$ is a permutation of ${1,2,3}$. By using this result, by equating to zero the coefficients in (I) of the partial derivatives $\partial_{x_l x_k} U$ and $\partial_{x_l x_k} V$, with $j \neq k$, we get the conditions: \begin{equation} (F1) \quad \frac{\partial f_l}{\partial X_k} + \frac{\partial f_k}{\partial X_l} = 0, \end{equation} \begin{equation} (F2) \quad \frac{\partial f_l}{\partial Y_k} + \frac{\partial f_k}{\partial Y_l} = 0. \end{equation} Now physics enters. First of all, let us note that by the same reason given above to justify Assumption 3, we must have that the probability current density must be invariant under global phase, that is $\mathbf{j}(e^{i \theta} D_{\mathbf{x}}^1 \psi(\mathbf{x},t))=\mathbf{j}(D_{\mathbf{x}}^1 \psi(\mathbf{x},t))$ for any $\theta \in \mathbb{R}$. Now, if we introduce the complex variable $Z_0=X_0 + i Y_0$ and the complex vector $\mathbf{Z}=(X_1+iY_1,X_2+iY_2,X_3+iY_3)$, and the vector functions $\mathbf{j}(Z_0,Z)=(j_1(Z_0,\mathbf{Z}),j_2(Z_0,\mathbf{Z}),j_3(Z_0,\mathbf{Z}))$ and $\mathbf{f}(Z_0,\mathbf{Z})=(f_1(Z_0,\mathbf{Z}),f_2(Z_0,\mathbf{Z}),f_3(Z_0,\mathbf{Z}))$, we can write equation (II) as \begin{equation} (III) \quad \mathbf{j}(Z_0,\mathbf{Z})=-\frac{i \hbar}{2m}\dot{p} (\overline{Z_0}Z_0)(\overline{Z_0} \mathbf{Z}- Z_0 \overline{\mathbf{Z}}) + \mathbf{f}(Z_0,\mathbf{Z}). \end{equation} If use equation (III), we see that the invariance of $\mathbf{j}$ under global phase and Borel's Lemma imply that for any $\theta \in \mathbb{R}$: \begin{equation} (F3) \quad \mathbf{f}(e^{i \theta} Z_0,e^{i \theta} \mathbf{Z})=\mathbf{f}(Z_0,\mathbf{Z}). \end{equation} Now, let us consider two reference frames $\Sigma$ and $\Sigma'$, with parallel axes and $\Sigma'$ moving at a constant velocity $\mathbf{V}$ withe respect to $\Sigma$, so that the coordinates $(\mathbf{x}',t')$ of an event in $\Sigma'$ are linked to the coordinates $(\mathbf{x},t)$ of the same event in $\Sigma$ by the Galilei transformations $\mathbf{x}'=\mathbf{x}-\mathbf{V}t$, $t'=t$. The wave functions $\psi(\mathbf{x},t)$ and $\psi'(\mathbf{x}',t')$ in the two reference frameworks are linked by \begin{equation} (W) \quad \psi'(\mathbf{x}',t')=\exp \left( -\frac{i}{\hbar} ( \mathbf{V} \cdot \mathbf{x}'+\frac{m V^2 t'}{2} ) \right) \psi(\mathbf{x}' + \mathbf{V} t',t'), \end{equation} as it is proved e.g. in Commins, Quantum Mechanics, $\S 4.8$ (the proof given there makes use of the fact that the probability density is given by $|\psi|^2$, but actually all what we need for the argument to work is that the probability density is a function of $|\psi|$, so that the proof is perfectly suitable in our present context; other suitable proofs of Equation (W) can be found in Esposito, Marmo and Sudarshan, From Classical to Quantum Mechanics, $\S 4.3$).

Clearly the probability density in the two frameworks is the same, so that \begin{equation} (D) \quad P'(\mathbf{x}',t')=P(\mathbf{x}'+\mathbf{V}t',t'). \end{equation}

As for the probability current density, since the probability mass $P d \mathbf{x}$ in the elementary volume $d \mathbf{x}$ near the point $\mathbf{x}$ moves with velocity $\mathbf{v}=\frac{1}{P}\mathbf{J}$, we have that the probability current density $\mathbf{J}=P\mathbf{v}$ transforms as $\mathbf{J}'=P'(\mathbf{v}+\mathbf{V})=P\mathbf{v}+P\mathbf{V}=\mathbf{J}+P\mathbf{V}$, that is \begin{equation} (C) \quad \mathbf{J}'(\mathbf{x}',t')=\mathbf{J}(\mathbf{x}'+\mathbf{V}t',t')+P(\mathbf{x}'+\mathbf{V}t',t')\mathbf{V}. \end{equation} Now, from equation (III) we have in the reference frame $\Sigma$ \begin{equation} (IV) \quad \mathbf{J}=-\frac{i\hbar}{2m}\dot{p}(\overline{\psi}\psi)(\overline{\psi} \nabla \psi - \psi \nabla \overline{\psi})+f(\psi,\nabla \psi), \end{equation}

and analogously in the reference frame $\Sigma'$

\begin{equation} (V) \quad \mathbf{J}'=-\frac{i\hbar}{2m}\dot{p}(\overline{\psi'}\psi')(\overline{\psi'} \nabla \psi' - \psi' \nabla \overline{\psi'})+f(\psi',\nabla \psi'). \end{equation} If insert equation (IV) and (V) in (C) and then we replace everywhere $\psi'$ by the expression given by equation (W) and we use property (F3), after some computations we get: \begin{equation} (VI) \quad \dot{p}(\overline{\psi}{\psi})\overline{\psi}\psi \mathbf{V}+\mathbf{f}\left(\psi,\nabla \psi -\frac{im\psi}{\hbar}\mathbf{V}\right)=p(\overline{\psi}{\psi}) + \mathbf{f}(\psi,\nabla \psi) \end{equation} where $\psi$ is computed in the point $(\mathbf{x}'+\mathbf{V}t',t')$.

Now, by choicing $\mathbf{V}$ in the direction of the first axis, that is $\mathbf{V}=(V_1,0,0)$, since $f_1$ does not depend on $X_1$ and $Y_1$, the first component of Equation (VI) gives $\dot{p}(\overline{\psi}{\psi})\overline{\psi}\psi=p(\overline{\psi}{\psi})$, which by the arbitrariness of $\psi$ gives the differential equation $y\dot{p}(y)=p(y)$, which has solution $p(y)= \alpha y$ for some $\alpha \in \mathbb{R}$. If we consider, as usual, normalized wave functions, that is functions $\psi$ such that $\int |\psi|^2 d \mathbf{x} = 1$, from the condition that the overall probability is one, we get $\alpha = 1$. This concludes the proof of the uniqueness of the probability density function.

Just for sake of completeness, let us finish the discussion of the conditions on the probability current density. From Equation (VI), and by using again Borel's Lemma, we now get that for any $\mathbf{v} \in \mathbb{R}^3$: \begin{equation} (F4) \quad \mathbf{f}(Z_0,\mathbf{Z}+i \mathbf{v} Z_0)=\mathbf{f}(Z_0,\mathbf{Z}). \end{equation} By using the fact that $\dot{p}=1$ we get from Equation (III) \begin{equation} (VII) \quad \mathbf{j}(Z_0,\mathbf{Z})=-\frac{i \hbar}{2m}(\overline{Z_0} \mathbf{Z}- Z_0 \overline{\mathbf{Z}}) + \mathbf{f}(Z_0,\mathbf{Z}). \end{equation} Now, what survives of Equation (I) is the condition \begin{equation} \sum_{l=1}^{3} \frac{\partial j_l}{\partial X_0} \partial_{x_l} U + \sum_{l=1}^{3} \frac{\partial j_l}{\partial Y_0} \partial_{x_l} V = 0, \end{equation} which, by using (VII) and the arbitrariness guaranteed by Borel's Lemma, gives \begin{equation} (F5) \quad \sum_{l=1}^{3} \left( \frac{\partial f_l}{\partial X_0} X_l + \frac{\partial f_l}{\partial Y_0} Y_l \right) = 0. \end{equation} As we see, we do not get a unique expression for the probability current density function, as we should expect, since by adding any field with zero divergence to $\mathbf{J}$ we still get a field satisfying the continuity equation. Indeed, equations (F1)-(F5), together with the fact that each $f_l$ does not depend on $X_l$ and $Y_l$, express the conditions that the function $\mathbf{f}(Z_0,\mathbf{Z})$ must satisfy in order for the function $\mathbf{j}$ given by Equation (VII) to be a "valid" probability current density. The standard expression of the probability current density corresponds to $\mathbf{f}=\mathbf{0}$, which gives: \begin{equation} \mathbf{J}=-\frac{i \hbar}{2m}(\overline{\psi} \nabla \psi- \psi \nabla \overline{\psi}). \end{equation} Another possible choice of the function $\mathbf{f}$ is \begin{equation} \mathbf{f}(Z_0,\mathbf{Z})=\frac{1}{m}(\overline{Z_0} \mathbf{Z}+Z_0 \overline{\mathbf{Z}}) \times \mathbf{s}, \end{equation} where $\mathbf{s} \in \mathbb{R}^3$ is a constant vector. This gives rise to the probability current density \begin{equation} \mathbf{J}=-\frac{i \hbar}{2m}(\overline{\psi} \nabla \psi- \psi \nabla \overline{\psi}) + \frac{1}{m}(\overline{\psi} \nabla \psi + \psi \nabla \overline{\psi}) \times \mathbf{s}, \end{equation} one obtains as the non-relativistic limit of the Dirac current density when the magnetic field is negligible and the system is in a spin eigenstate with spin vector $\mathbf{s}$: see the already quoted work by Holland Uniqueness of Conserved Currents in Quantum Mechanics.

Even though the proof given here has not the (extraordinary!) elegance of the one given by EuYu in his answer, I think there is some worth in it, not only because it is elementary, since it does not make use of functional analysis, but also because it is very much in the spirit and in the style that Bohm adopts in his wonderful book Quantum Theory. His aim in that work is to build quantum mechanics step after step, without framing it in an axiomatic way from the very beginning (as almost all modern texts do), but highlighting the crucial conceptual bricks which constitute the bulk of the theory and that should be always be present in the mind of the physicist who does limit himself to apply a formidable set of mathematical formalisms to the physical system under study, but who also asks himself the conceptual meaning of the theory and its possible developments. Thank you, Master Bohm for this incredible gift to all of us!