# Obtaining the total charge density from a multi-species many-body wavefunction

One can exactly solve the two-body wavefunction describing the interaction of an electron and proton through the following Hamiltonian

$$H=-\frac{\hbar^2}{2m_p}\nabla^2_{p} + -\frac{\hbar^2}{2m_e}\nabla^2_{e} -\frac{e^2}{\vert \mathbf{r}_p-\mathbf{r}_e\vert}$$

The resulting two-particle wavefunction is then a function of both the electron and proton coordinates $\psi(\mathbf{r}_p,\mathbf{r}_e)$. The strange thing about this wavefunction is that it describes the probability amplitude of both the electron and proton, which have different charge. However, textbooks often ignore this fact, and state the the total charge density is given by $\vert \psi(\mathbf{r}_p,\mathbf{r}_e)\vert^2$. This cannot possibly be correct, because if we integrate the total charge density over all space we should get identically zero for the electron-proton system.

My question is: what is the correct method for obtaining the charge density of a many-body wavefunction which has different charge species contained in it? If it were purely a wavefunction of electrons, then you could look at the one-body density matrix which contains all the information you need:

$$n(\mathbf{r}_1) = -Ne\int d\mathbf{r}_2...d\mathbf{r}_k \vert \psi(\mathbf{r}_1,\mathbf{r}_2,...,\mathbf{r}_k)\vert^2$$

But what is the correct procedure for multi-species wavefunctions (i.e. wavefunctions of particles with different charge/quantum numbers)? Naively I would assume the charge density is given by the following equation, but I am not sure if this is rigorously true.

$$\rho(\mathbf{r}) = +e\left(\int d\mathbf{r}_e \vert \psi(\mathbf{r}_p,\mathbf{r}_e)\vert^2\right)\rvert_{\mathbf{r}_p=\mathbf{r}} -e\left(\int d\mathbf{r}_p \vert \psi(\mathbf{r}_p,\mathbf{r}_e)\vert^2\right) \rvert_{\mathbf{r}_e=\mathbf{r}}$$

Your final expression, $$\rho(\mathbf{r}) = +e\left(\int d\mathbf{r}_e \vert \psi(\mathbf{r}_p,\mathbf{r}_e)\vert^2\right)\rvert_{\mathbf{r}_p=\mathbf{r}} -e\left(\int d\mathbf{r}_p \vert \psi(\mathbf{r}_p,\mathbf{r}_e)\vert^2\right) \rvert_{\mathbf{r}_e=\mathbf{r}},$$ is indeed correct, but it's more helpful to rephrase it in the form $$\rho(\mathbf{r}) = \sum_{j=p,e} \int d\mathbf{r}_e d\mathbf{r}_p \: q_j \delta (\mathbf r-\mathbf r_j) \: \vert \psi(\mathbf{r}_p,\mathbf{r}_e)\vert^2,$$ which makes it much more obvious how to connect it to a full formal operator expectation value, namely $$\rho(\mathbf r) = \langle \psi | \hat \rho(\mathbf r) |\psi\rangle$$ where $$\hat \rho(\mathbf r) = \sum_{j=p,e} q_j \delta (\mathbf r-\hat {\mathbf r}_j).$$ This operator version is then the obvious quantum-mechanical version of the classical charge density $\rho(\mathbf r) = \sum_{j=p,e} q_j \delta (\mathbf r-{\mathbf r}_j)$ of two point charges $q_e$ and $q_p$ at positions $\mathbf r_e$ and $\mathbf r_p$, respectively.

And, of course, it generalizes transparently to the case of $N$ particles, both distinguishable and indistinguishable, and it reduces to the (correct) formula you give for indistinguishable particles when you do have the symmetry in place.

Textbooks often ... state the the total charge density is given by $|\psi({\bf r}_p,{\bf r}_e)|^2$. This cannot possibly be correct, because if we integrate the total charge density over all space we should get identically zero for the electron-proton system.

Any textbook making this claim (do you have an example?) is seriously confused. First of all, $|\psi({\bf r}_p,{\bf r}_e)|^2$ isn't a charge density at all, but rather a probability density. Second of all, even if you multiply by $e$, then $e |\psi({\bf r}_p,{\bf r}_e)|^2$ still isn't a spatial charge density; it's a charge density over configuration space, which is conceptually very different. For example, spatial charge densities have dimension charge/volume, while for an $n$-particle system, configuration-space charge densities have dimension charge/(volume)$^n$.

Probably the best way to physically interpret the quantity $e |\psi({\bf r}_p,{\bf r}_e)|^2$ is to hold ${\bf r}_p$ fixed and think of it as a function $e |\psi_{{\bf r}_p}({\bf r}_e)|^2$ of a single (vector) argument. In this case, the function can be roughly thought of as the conditional spatial charge density of the electron alone given the location of the proton (up to a sign that depends on your conventions).

(You could of course instead hold ${\bf r}_e$ fixed and consider $e |\psi_{{\bf r}_e}({\bf p}_e)|^2$ as the proton's charge density conditioned on the electron's location, but this is much less useful because the proton is so much heavier that in practice its wavefunction is much more localized than the electron's, so in order to understand atomic (as opposed to nuclear) physics you can treat it as a classical point particle.)

The interpretation of a "conditional spatial density" is a bit subtle; note that it doesn't have the units of a true spatial density, and in order to convert it into a true spatial density you need to integrate it with respect to the parameter ${\bf r}_p$ over some possible volume $V_p$ where the proton could be located. Integrating ${\bf r}_p$ over all space (the second term in your expression) gives the marginal charge density of the election alone, with no restrictions on the proton's location.

Your full expression for $\rho({\bf r})$ (with the obvious generalization for more species) is indeed the correct expectation value for the total charge density operator at point ${\bf r}$, which is the most natural way to assign a "spatial charge density" to the system. But note that even in the single-particle case, naively considering $e |\psi({\bf r})|^2$ as a charge density is a useful heuristic but can be a bit dangerous; a wavefunction is fundamentally quantum mechanical, and an electron maintains a particle-like nature in certain senses. (E.g. if you precisely measure its location then you'll always observe it to be localized; you'll never measure it to simultaneously be on opposite sides of the nucleus, as with a classical charged fluid.)

• I should be more precise, textbooks will separate the wavefunction into relative and c.o.m. parts, with the relative part $\vert\psi(r_e-r_p)\vert^2$ given the role of the charge density. It is inherently assumed that the proton has a (classically) fixed position so r_p=const. That's practically true, but if one were to look at the case of an electron and positron (or muon), where that assumption breaks down, then the well known hydrogen orbitals do not describe the charge density nor position of the electron. – KF Gauss May 16 '18 at 1:41
• If I recall correctly, the charge density can be written as an operator that obeys all the usual Heisenberg time evolution, so it really isn't that dangerous right? – KF Gauss May 16 '18 at 1:45
• Actually after thinking about it, would the formula be derivable from $\rho= \Sigma Z_i \psi^{\dagger}_i(r) \psi_i(r')$ and the canonical commutation relations for fermions $[\psi^{\dagger}_i(r),\psi_i(r')]=\delta_{r,r'}$ – KF Gauss May 16 '18 at 1:56
• @user157879 Writing the wavefunction as $\psi({\bf r}_c - {\bf r}_p)$ with ${\bf r}_p$ a fixed constant completely changes its interpretation, and even its units, so that's a very different question. In that case it's just a one-particle wavefunction and the effective charge density for the electron is indeed just $e |\psi({\bf r}_c - {\bf r}_p)|^2$ with no integration needed. – tparker May 16 '18 at 1:57
• sure but while the relative coordinate wavefunction is still valid for an electron-positron system and is a traditional hydrogen orbital, there the interpretation would no longer be the effective charge density of electron right? – KF Gauss May 16 '18 at 2:00

The charge density is the probability of finding a particle at position $\vec r$, regardless of where the other particles are. For the wave function above this is $\int{ d\vec r' (e |\psi(\vec r,\vec r') |^2 - e |\psi(\vec r',\vec r )|^2})$. For an $n$-particle wave function this is $\sum_{i=1}^n{ q_i \int{ d\vec r_1 .. d\vec r_{i-1} d\vec r_{i+1} .. d\vec r_n |\psi(\vec r_1,..,\vec r_{i-1},\vec r,\vec r_{i+1}..,\vec r_n |^2} }$. The $q_i$ have to be factored in manually in the Schrödinger - also in the Dirac - picture.

Triggered by tparkers answer: $|\psi(\vec r_p,\vec r_e) |^2$ gives the joint probability of finding the proton at $\vec r_p$ and the electron at $\vec r_e$.