Suppose I have an N particle wavefunction including spin. How can I coarse grain the spin degrees of freedom? I think simply disregarding the spin parts do not work.
N-particle wave functions for, say, spin-½ particles are constructed from single-particle spin-orbitals, each such orbital consisting of the product of a single-particle orbital φi(r) and a single-particle spin function η(σ); it is denoted by ψi(x), where x = r,σ, in which σ = ±, where + corresponds to spin-up and - to spin-down (thus η(σ) is an eigenstate of, say, the z-component of the single-particle spin operator, where z refers to the arbitrarily chosen spin quantisation axis -- in the presence of an external uniform magnetic field H, the direction of this field defines the z direction). One has ψi(x) = φi(r) η(σ). With xj = rj,σj, it is often convenient to write ψi(xj) = φi(rj) α(j) for σ = + and ψi(xj) = φi(rj) β(j) for σ = -.
An N-particle wave function is therefore denoted as Ψ(x1,x2,...xN), consisting of an anti-symmetrised linear combination of N spin-orbitals of the kind introduced above. Clearly, Ψ(x1,x2,...xN) depends on N continuous variables {r1, r2, ..., rN}, and N discrete variables {σ1, σ2, ..., σΝ}. Course graining of spin at this level, that is at the level of the total wave function, is meaningless, since the symmetric / anti-symmetric nature of Ψ(x1,x2,...xN) under the operation of an even / odd N-permutation of {x1,x2, ...,xN} crucially depends on the discrete values assigned to the elements of the set {σ1, σ2, ..., σΝ}. Naturally for any given spin configuration {σ1, σ2, ..., σΝ}, the distribution Ψ*(x1,x2,...xN) Ψ(x1,x2,...xN) may be coarse grained.
For details, I refer you to the classic text Methods of Molecular Quantum Mechanics, 2nd edition, by R McWeeny (Academic Press, New York, 2001). In particular, consider the section on Young tableau. In this book density matrices of different order are defined (the largest order being N, for which the corresponding matrix is equal to Ψ*(x1,x2,...xN) Ψ(x1,x2,...xN), and the smallest non-trivial order 1, for which the corresponding density matrix coincides with the Slater-Fock density matrix, often denoted by ρ(x1,x2)), in particular spin density matrices, each of which may be coarse grained for a given configuration of the the relevant spin indices.
One last remark, the set of above-mentioned anti-symmetrised N-particle wave functions are directly appropriate for dealing with the cases where spin is a good quantum number. Nonetheless, when this is not the case (for instance in the case of spin-orbit coupling, where the operator to be simultaneously diagonalised with the relevant N-particle Hamiltonian is the total J, where J = L + S) these wave functions can be used as basis functions in for instance perturbative calculations (think of the LS- and JJ-coupling schemes, which can be applied for arbitrary values of N). □