I know that if I have a single fermion then the expectation value of the spin along the $\hat{n}$ direction can be computed by $\sigma \cdot \hat{n}$, where $\sigma_i$ are the pauli matrices.
Now if I have a multiparticle system, is there a generalization of the pauli matrices (lets call them $\rho_i$) such that the expectation value of total spin along the $\hat{n}$ direction can be computed by $\rho \cdot \hat{n}$?
Assuming you want an $N$ particle system of spin 1/2 fermions, then the corresponding operator is $\sum_{j=1}^N (\vec \sigma_j \cdot \hat n)$, where $\vec \sigma_j = \hat x \sigma_{jx}+ \hat y \sigma_{jy}+\hat z\sigma_{jz} $ are the Pauli operators for particle $j$. You need $2^N$ basis states, and the usual Pauli matrices will tell you what the corresponding matrix elements are. As an example $\langle \uparrow \downarrow | \sigma_{2x} |\uparrow \uparrow\rangle = 1$, Since the usual Pauli operator $\sigma_x$ flips the single spin, $\sigma_{2x}$ flips the second spin.
Typically, it is easier to use the operators to directly calculate the needed matrix elements rather than to write out $2^N \times 2^N$ matrices.
There are two quantum numbers you will have to consider:
intrinsic spin
azimuthal
If you use the z axis, then the x and y axis' quantum numbers will be undefined as per the Pauli exclusion principle.
For two fermions there are four spin states:
αα, ββ, αβ, βα
The first two are symmetric, the last two are not symmetric or antisymmetric when electrons are exchanged so linear combination has to be made.
$\chi_{-}:=\frac{1}{\sqrt{2}}\left[\alpha(1)\beta(2)-\alpha(2)\beta(1)\right]$
$\chi_{\alpha}:=\alpha(1)\alpha(2),\quad\chi_{\beta}:=\beta(1)\beta(2),\quad\chi_{+}:=\frac{1}{\sqrt{2}}\left[\alpha(1)\beta(2)+\alpha(2)\beta(1)\right]$
These all have the same energy level.
