In the electron density, what exactly do you mean by $\sigma$ mathematically? Suppose there is a many-body eigenvalue problem
$$
H(r_1,\cdots,r_n) = - \frac{1}{2} \nabla^2|_{r_1\cdots,r_n} + V(r_1,\cdots,r_n).
$$
This potential is symmetric and the hamiltonian does not depend on spin, but the solution should be for spin 1/2 particles. I am struggling to interpret the symbol $\sigma_i$ mathematically. Perhaps you can help?
It seems many authors take the «non-degenerate ground-state electron density» to be the following symbol
$$
n(r) = \sum_{\sigma_1} \sum_{\sigma_2} \cdots  \sum_{\sigma_n} \int |\psi(r_1,\sigma_1 , r_2,\sigma_2,\cdots, r_n,\sigma_n)|^2 dr_2 \cdots  dr_n.
$$
In that context, I will present 3 scenarios of how $\psi(r_1,\sigma_1 , r_2,\sigma_2,\cdots, r_n,\sigma_n)$ might be related  to the hamiltonian. Under each scenario, I have entered a few questions. It would be helpful if you could indicate which scenario is correct, and answer the corresponding questions under the correct scenario.
Scenario (1): For every $r_1,\sigma_1,\cdots,r_n,\sigma_n$,
$$
H\psi(r_1,\sigma_1 , r_2,\sigma_2,\cdots, r_n,\sigma_n) = E_{\sigma_1,\cdots, \sigma_n} \psi(r_1,\sigma_1 , r_2,\sigma_2,\cdots, r_n,\sigma_n),
$$
and $E_{\sigma_1,\cdots, \sigma_n}$ is different for every combination of $\sigma_1,\cdots,\sigma_n$.
(a) Why do all $\psi(r_1,\sigma_1 , r_2,\sigma_2,\cdots, r_n,\sigma_n)$ appear in the ground state density, when the ground state can only be one of the $\psi(r_1,\sigma_1 , r_2,\sigma_2,\cdots, r_n,\sigma_n)$ (by virtue of nondegeneracy)? Do we simply define that the $2^n$ smallest eigenvalues correspond to the ground state? (b) Does $\psi(r_1,\sigma_1 , r_2,\sigma_2,\cdots, r_n,\sigma_n)$ necessarily factor into a spatial part and a spin part for each $\sigma_1,\cdots,\sigma_n$? (c) Is there reason to expect that $\psi(r_1,\sigma_1 , r_2,\sigma_2,\cdots, r_n,\sigma_n)$ antisymmetric for every $\sigma_1,\cdots,\sigma_n$?
Scenario (2) For every $r_1,\sigma_1,\cdots,r_n,\sigma_n$,
$$
H\psi(r_1,\sigma_1 , r_2,\sigma_2,\cdots, r_n,\sigma_n) = E_{\sigma_1,\cdots, \sigma_n} \psi(r_1,\sigma_1 , r_2,\sigma_2,\cdots, r_n,\sigma_n),
$$
but $E_{\sigma_1,\cdots, \sigma_n}$ is the same for every $\sigma_1,\cdots, \sigma_n$. (a) Is there any reason to expect that the eigenspace of $E_{\sigma_1,\cdots, \sigma_n}$ is at least $2^n$ or eactly $2^n$? (b) Do physicist say that the ground state is degenerate if and only if the dimension of the eigenspace is larger than $2^n$? (c) Is there reason to expect that $\psi(r_1,\sigma_1 , r_2,\sigma_2,\cdots, r_n,\sigma_n)$ antisymmetric for every $\sigma_1,\cdots,\sigma_n$?
Scenario (3): There a function $\psi$ which is the ground state that is a formal linear combination of terms
$$
\psi(r_1,\cdots,r_n) = \sum_{\sigma_1} \sum_{\sigma_2} \cdots  \sum_{\sigma_n} \psi(r_1,\sigma_1 , r_2,\sigma_2,\cdots, r_n,\sigma_n),
$$
such that $H\psi = E\psi$, but the individual terms $\psi(r_1,\sigma_1 , r_2,\sigma_2,\cdots, r_n,\sigma_n)$ no not necessarily solve this eigenvalue problem for any $r_1,\sigma_1,\cdots,r_n,\sigma_n$. If this interpretation is correct, I have some questions again: (a) In this case the functions $\psi(r_1,\sigma_1 , r_2,\sigma_2,\cdots, r_n,\sigma_n)$ do not correspond to different eigenvalues for each $\sigma_1,\sigma_2,\cdots,\sigma_n$, so I cannot see a reason why they would be orthogonal. Are they antisymmetric or even normalised for some reason? (b) Is the function $\psi(r_1,\cdots,r_n)$ antisymmetric or normalised in this case? (c) What 'gained' by writing the ground state as as sum of these $2^n$ functions? (d) How can I arrive at the expression for the electron density in this case?
Hoping for some answers!
 A: Perhaps a very concrete example will be helpful:
Consider the ground state of a hydrogen atom, that is, the electron in a $1s$ orbital. The spatial wavefunction is:
$$\psi(r) = \phi_{1s}(r)$$
But this isn't the electron's wavefunction because the electron also has a spin. Let's suppose the electron is spin up. Then (just like how a particle of definite position has a position wavefunction that is a Dirac delta) the spin wavefunction is:
$$\chi(\sigma) = \delta(\sigma,+1/2)$$
where $\delta(\sigma_1,\sigma_2)$ is a Kronecker delta (we're also specifically working in a basis where spins take values $\pm 1/2$).
There are two degenerate ground states, one of which is:
$$\psi_+(r,\sigma) = \phi_{1s}(r) \delta(\sigma, +1/2)$$
and the other has a minus sign in the spin delta function. These are both individually eigenvalues of $H$ and we can write $H\psi_{\pm}=E^{1s}_{\pm}\psi_{\pm}$. Now in the time-reversal invariant case these are degenerate but we could add a magnetic field to split them via the Zeeman effect. Okay, so with this setup let's quickly go through the OP's scenarios (out of order):


*This scenario is that $E^{1s}_+ = E^{1s}_-$ but we've said we can split them with a Zeeman term.


*This doesn't really make sense. The Hamiltonian $H$ contains a Zeeman term $H_{Zeeman} \sim \sigma\cdot B$. The reduced wavefunction you suggest no longer has spin space information, it's not in the vector space on which the Hamiltonian acts.
This leaves us with:

*

*This is accurate, but you seem confused about what exactly you're saying when you compute the electron density. In particular, you don't use both the wavefunction $\psi_+$ and $\psi_-$. Suppose the ground state has the spin aligned with our $B$ field so is $\psi_+$. Then $\sum_\sigma |\psi_+|^2 = \sum_\sigma \delta^2(\sigma,1/2)|\phi_{1s}|^2$. A Kronecker delta is either $0$ or $1$ so $\delta^2 = \delta$ and $\sum_\sigma \delta(\sigma,\sigma')=1$ so the net effect of all of this is just to recover $|\psi_{1s}|^2$ which is what we want. Note that we didn't need to consider $\psi_-$ anywhere for the ground state density calculation.

