Why do we integrate square of wave function over N-1 electron coordinates?

Quote form F.Jensen's Introduction to computational chemistry:

"The electron density is the square of the wavefunction, integrated over N-1 electron coordinates..."

Why do we integrate over N-1 instead of N electorn coordinates? What would happen if we integrated over N coordinates? I suppose we would get a norm if we integrated over N coordinates, or 1 if wavefunction is normed. But why N-1?

• Well I am not an expert in this but I suppose that if you integrate over all coordinates you get a number. The electron density $\rho$ is a function of $(r,\theta,\phi)$, so if you exclude the integration over one set of coordinates then you get a function of one set of coordinates (which is what you want), just like how if you integrate $\int_0^1\int_0^1 (3x + 4y)\,\mathrm{d}x\mathrm{d}y$ you get a number, but if you integrate $\int_0^1 (3x + 4y)\,\mathrm{d}x$ then you get a function of $y$. – orthocresol Sep 19 '16 at 12:23

$\newcommand{\el}{_\mathrm{e}} \newcommand{\dif}{\mathrm{d}} \newcommand{\braket}{\langle{#1}\vert{#2}\rangle}$Consider an $n$-electron system in an arbitrary electronic state represented by the electronic wave function $\psi\el(\vec{q}_{1}, \vec{q}_{1}, \dotsc, \vec{q}_{n})$, where $\vec{q}_{i}$ stands for joint spin-spatial coordinate of $i$-th electron, $\vec{q}_{i} = \{ m_{si}, \vec{r}_{i} \}$.
Integrating the square of the absolute value of the electronic wave function over the spatial coordinates of all the electrons except for electron-one and summing up over the spin coordinates of all of them without any exception gives the probability density of finding electron-one with arbitrary spin at $\vec{r}_{1}$ while having the remaining $n - 1$ electrons with arbitrary spins at arbitrary positions, \begin{equation*} \rho_{\mathrm{e-1}}(\vec{r}_{1}) = \sum_{m_{s1}} \sum_{m_{s2}} \dotsi \sum_{m_{sn}} \iint \dotsi \int \left| ψ\el(\vec{q}_{1}, \vec{q}_{2}, \dotsc, \vec{q}_{n}) \right|^{2} \dif \vec{r}_{2} \dif \vec{r}_{3} \dotsc \dif \vec{r}_{n} \, . \end{equation*} And yes, if we "integrate" over all coordinates in this expression we'll just get an inner product of the wave function with itself, $$\braket{\psi\el}{\psi\el} = \sum_{m_{s1}} \sum_{m_{s2}} \dotsi \sum_{m_{sn}} \iint \dotsi \int \left| ψ\el(\vec{q}_{1}, \vec{q}_{2}, \dotsc, \vec{q}_{n}) \right|^{2} \dif \vec{r}_{1} \dif \vec{r}_{2} \dif \vec{r}_{3} \dotsc \dif \vec{r}_{n} \, ,$$ which is usually equals to 1 since the wave function is usually normalized to unity.
Electrons are, of course, indistinguishable, so that their labeling is arbitrary and can be changed any time. Thus, the expression above should work equally well for calculating the probability density of finding electron-two, electron-three, etc., or, in other words, it gives the probability density of finding one particular electron, no matter which, with arbitrary spin at $\vec{r}_{1}$ while having the remaining $n - 1$ electrons with arbitrary spins at arbitrary positions.
The electron density, defined as the probability density of finding any of the $n$ electrons with arbitrary spin at $\vec{r}_{1}$, is evidently $n$ times as large as the probability of finding one particular electron, \begin{equation} \label{eq:electron_density} \rho(\vec{r}_{1}) = n \sum_{m_{s1}} \sum_{m_{s2}} \dotsi \sum_{m_{sn}} \iint \dotsi \int \left| ψ\el(\vec{q}_{1}, \vec{q}_{2}, \dotsc, \vec{q}_{n}) \right|^{2} \dif \vec{r}_{2} \vec{r}_{3} \dotsb \dif \vec{r}_{n} \, . \end{equation} Note here that, for a wave function that is normalized to unity, if we "integrate" over all coordinates in this expression we'll just get the number of electrons $n$, which is one of the well-known properties of the electron density, \begin{equation} \int \rho(\vec{r}_{1}) \dif \vec{r}_{1} = n \, . \end{equation}