Say we are in a box, so that ##\{\hat{C}(\vec{p},\sigma),\hat{C}^{\dagger}(\vec{p}',\sigma')\}=\delta_{\vec{p},\vec{p}'} \delta_{\sigma,\sigma'}## with the ##\delta##'s being Kronecker ##\delta##'s with values ##1## and ##0## rather than Dirac-##\delta## distributions. Then you can just do the calculation: By definition
$$\hat{N}=\hat{C}^{\dagger}(\vec{p},\sigma) \hat{C}(\vec{p},\sigma)$$
leading to
$$\hat{N}^2=\hat{C}^{\dagger} \hat{C} \hat{C}^{\dagger} \hat{C}=\hat{C}^{\dagger} [\{\hat{C},\hat{C}^{\dagger} \}-\hat{C^{\dagger} \hat{C}}] \hat{C}=\hat{C}^{\dagger} \hat{1} \hat{C}=\hat{N},$$
because ##\{C,C \}=2 \hat{C}^2=0##.
Another argument is that ##\hat{N}## has the eigenvalues 0 and 1, i.e.,
$$\hat{N}^2=\sum_{n=0}^1 |n \rangle \langle n|n^2 =\sum_{n=0}^1 |n \rangle \langle n|n=\hat{N},$$
because both ##0^2=0## and ##1^2=1##, i.e., for both eigenvalues ##n^2=n##.