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}}$$