Cardinality of hypersphere over $\mathbb{F}_p$ is always divisible by $p$? I recently read an interesting post here (the second answer to this) about proving quadratic reciprocity by considering the set
$$S_p^{n-1} = \{(x_1,\dots,x_n)\in\mathbb{F}_p^n : x_1^2+\dots +x_n^2 = 1\}$$
Which is just the equation for an $(n-1)-$sphere with entries from a finite field.
Then, consider
$$N_p(n):= |S_p^{n-1}|$$
which just gives how many possible solutions to $x_1^2+\dots +x_n^2 = 1$ there are in $\mathbb{F}_p$.
It turns out that
$$N_p(2) \equiv 1 \mod p \iff p \equiv 3 \mod 4$$
and
$$N_p(2) \equiv -1 \mod p \iff p \equiv 1 \mod 4$$
for all $p$, which is interesting.
After reading the proof in the link above, I see why this is the case.
What I do not understand is that for $n>2$, I wrote some code and it seems
$$N_p(n) \equiv 0 \mod p$$
for all $p$.
So my question is

Why is it the case that $p$ divides $N_p(n)$ for all $n>2$? What does having a certain amount of solutions on a hypersphere over $\mathbb{F}_p$ have to do with being divisible by $p$?

 A: Here's a method which is overkill, but quite fun. Since you tagged your question as 'algebraic geometry', it doesn't seem entirely uncalled for.
Let

*

*$X_n=V(x_1^2+\cdots+x_n^2-1)\subseteq \mathbb{A}^n_{\mathbb{F}_p}$

*$\overline{X}_n=V(x_1^2+\cdots+x_n^2-x_{n+1}^2)\subseteq \mathbb{P}^n_{\mathbb{F}_p}$.

Let us note that $\overline{X}_n-X_n$ agrees set theoretically with the closed subscheme

*

*$Z_n=V(x_1^2+\cdots+x_n^2)\subseteq \mathbb{P}^{n-1}_{\mathbb{F}_p}$.

In particular, we see that
$$\# X_n(\mathbb{F}_p)=\# \overline{X}_n(\mathbb{F}_p)-\# Z_n(\mathbb{F}_p).$$
Thus, it suffices to show that
$$\# \overline{X}_n(\mathbb{F}_p)=\# Z_n(\mathbb{F}_p)=1\mod p.$$
Let us make the observation that $\overline{X}_n$ and $Z_n$ are degree $2$ hypersurfaces in $\mathbb{P}^n_{\mathbb{F}_p}$ and $\mathbb{P}^{n-1}_{\mathbb{F}_p}$ respectively. We then start with the following lemma.

Lemma: Let $k$ be any field and let $S\subseteq \mathbb{P}^m_k$ be a degree $2$ hypersuface. If $m\geqslant 2$ then $$H^i(S,\mathcal{O}_S)=\begin{cases}k & \mbox{if}\quad i=0\\ 0 & \mbox{if}\quad \text{otherwise.}\end{cases}$$

Proof: This is a simple calculation. We have the usual exact sequence
$$0\to \mathcal{I}\to \mathcal{O}\to i_\ast\mathcal{O}_S\to 0,$$
where $\mathcal{I}$ is the ideal sheaf of $S$, $\mathcal{O}$ is short for $\mathcal{O}_{\mathbb{P}^m_k}$, and $i\colon S\hookrightarrow \mathbb{P}^m_k$ is the inclusion. Since $S$ has degree $2$ we have that $\mathcal{I}\cong \mathcal{O}(-2)$. The long exact sequence in cohomology gives
$$0\to H^0(\mathcal{O}(-2))\to H^0(\mathcal{O})\to H^0(i_\ast\mathcal{O}_S)\to \cdots\to H^j(\mathcal{O}(-2))\to H^j(\mathcal{O})\to H^j(i_\ast\mathcal{O}_S)\to\cdots,$$
where for shorthand I'm writing $H^i(-):=H^i(\mathbb{P}^m_k,-)$. Now, the standard computations of cohomology of $\mathbb{P}^m_k$ then imply that since $m\geqslant 2$ that $H^j(\mathcal{O})=H^j(\mathcal{O}(-2))=0$ for all $j$, except that $H^0(\mathcal{O})=k$. From this one easily deduces the desired claim since $H^j(S,\mathcal{O}_S)\cong H^j(i_\ast\mathcal{O}_S)$. $\blacksquare$
Using this lemma, and Fulton's trace formula (see [Mustata, Theorem 1.1]) we easily deduce the following which, as previously commented, implies your claim.

Proposition: Let $S$ be a degree $2$ hypersurface in $\mathbb{P}^m_{\mathbb{F}_p}$. If $m\geqslant 2$ then
$$\# S(\mathbb{F}_p)=1\mod p.$$

Proof: Fulton's trace formula implies that
$$\# S(\mathbb{F}_p)=\sum_{i=0}^{\dim(S)}(-1)^i\mathrm{tr}\left(\mathrm{Frob}_p|H^i(S,\mathcal{O}_S)\right)\mod p.$$
By the Lemma we know that the only non-vanishing term on the right is $H^0(S,\mathcal{O}_S)$ which is $\mathbb{F}_p$. Since $\mathrm{Frob}_p$ is the identity on $\mathbb{F}_p$ we deduce that that the right-hand side is $1$, from whence the claim follows. $\blacksquare$

Using the same methods, you can easily deduce the following.

Theorem: Let $d\geqslant 1$ be an integer and $q$ a power of $p$. Then, for any polynomial
$f(x_1,\ldots,x_n)\in\mathbb{F}_q[x_1,\ldots,x_n]$ such that $n>d$,
and for any $c\in \mathbb{F}_q^\times$ one has that
$$\#\left\{(x_1,\ldots,x_n)\in\mathbb{F}_q^n:
 f(x_1,\ldots,x_n)=c\right\}=0\mod p.$$

While not often thought about geometrically, this is all a special case of the classical theorem of Chevalley--Warning:

Theorem (Chevalley--Warning): Let $q$ be a power of $p$, and $f_1,\ldots,f_r$ elements of $\mathbb{F}_q[x_1,\ldots,x_n]$ of total
degrees $d_1,\ldots,d_m$. If $n>d_1+\cdots+d_m$, then
$$\#\left\{(x_1,\ldots,x_n)\in
\mathbb{F}_q^n:f_j(x_1,\ldots,x_n)=0\text{ for
> }j=1,\ldots,r\right\}=0\mod p.$$

You can see [CGS, Theorem 1.1] for an elementary proof.
But, what's cool is that the Fulton trace formula can give results in different directions. For instance, we have the following trivial corollary to Fulton's trace formula (which we used a special case of above in the Proposition).

Corollary(to Fulton's trace formula): Let $q$ be a power of $p$. Suppose that $X$ is a (geometrically connected) proper variety
over $\mathbb{F}_q$ such that $H^i(X,\mathcal{O}_X)$ is $0$ for $i>0$.
Then,
$$\# X(\mathbb{F}_q)=1\mod p.$$

This, in some sense, is more interesting to me than the Chevalley--Warning theorem, even though the latter more directly answers your question, because it really heavily implies the beautiful realization that the geometry of equations affects their number of solutions. This realization being the lynchpin for many of the most deep results we currently know about numbers of solutions to equations over finite fields (e.g. the Weil conjectures).
References:
[CGS] http://alpha.math.uga.edu/~pete/Chevalley_Warning_on_the_Boundary.pdf
[Mustata] http://www.math.lsa.umich.edu/~mmustata/lecture6.pdf
A: To answer your second question of why counting solutions to such an equation mod $p$ might involve congruences mod $p$, we can consider the geometry of the situation. In particular, the Weil conjectures, which motivated massive developments in algebraic geometry in the $20$th century, are precisely about counting points of varieties (things cut out by polynomials) over finite fields. They tell you that the point counts over different finite fields inform you about the intrinsic geometry of your variety, which even translates to geometry about the complex solutions of the same equations.
This is to say, counting points over finite fields can be incredibly deep stuff. For us however, our variety is the sphere, and we want to count the points over $\mathbb{F}_p$.
From this perspective, we should view the $n-1$ sphere as the space of cosets of the orthogonal group, mod the smaller orthogonal group. So $SO(n)$ is $n\times n$ orthogonal matrices of determinant one, and we want the cosets of $SO(n-1)$ embedded as $1$ cross the last $n-1\times n-1$ coordinates being an special orthogonal matrix. Note that it makes sense for this over any field, we just look at the solutions in different fields. Then we can say $S^{n-1}\cong SO(n)/SO(n-1)$, with its $k$ valued points to be the coset space of the $k$ valued points of these orthogonal groups (think about the action of $SO(n)$ on $S^{n-1}$ and the orbit stabiliser theorem). For instance, when $k=\mathbb{R}$ is the reals, then this is the usual sphere.
So then, by the orbit stabiliser theorem, we get that our point count for the sphere is the quotient of the orders of these two orthogonal groups. Here we run into something slightly strange, in that there are two different orthogonal groups over finite fields, so the answer will depend on the $4$ parity of the prime we pick, as we see in the quadratic reciprocity proof. But ignoring that, the point is the higher orthogonal groups have known orders, and the power of $p$ dividing the number of points of $O(2n+1)$ is $q^{n^2}$, and dividing either version of $O(2n)$ is $q^{n(n-1)}$, so we see as soon as $k$ is larger than than $2$, these powers of $p$ don't cancel out, so we are left with powers of $p$ dividing the number of solutions of the sphere.
To justify why $SO(n)$ acts transitively on the sphere, we need some theory of quadratic forms (there might be an easier proof), so this isn't a complete proof of the assertion, but this tells you why we expect such numbers to be divisible by $p$ for large $n$. As far as sources for this stuff, the wikipedia page for the orthogonal group gives the numbers we are using/detail on quadratic forms, and the wikipedia for Weil conjectures gives some nice examples of point counts, and what they indicate geometrically.
Also, its worth pointing our that the method Pig used works more generally for Fermat hypersurfaces, whereas this special case uses the exponent $2$ in a crucial way, and these computations for these hypersurfaces are some of the original motivating examples of Andre Weil that he put forward in the article where he made these beautiful conjectures.
A: I will illustrate this for $n = 3$ for simplicity, but the proof directly generalizes. Let's assume $p > 2$ as well (I didn't check $p = 2$)
First, note that number of solutions to $x^2 = a \mod p$, denoted $N(x^2=a)$, is $1 + \chi(a)$, where $\chi(a) = \left(\frac{a}{p}\right)$ is the Legendre symbol. This allows us to rewrite $N_p(3)$ as a character sum by
$$
\begin{align*}
N_p(3) &= \sum_{a + b+ c = 1} N(x^2 = a) N(x^2 = b) N(x^2 = c) \\
&= \sum_{a+b+c = 1} (1 + \chi(a)) (1 + \chi(b)) (1 + \chi(c))
\end{align*}
$$
Expanding the product, there are three types of terms

*

*"Main term": $\sum_{a+b+c = 1} 1 \equiv 0 \mod p$.

*

*This counts number of solutions to $a+b+c = 1 \mod p$. You can pick anything you like for $a,b$, which then fixes $c = 1 - a - b$, so there's $p^2$ solutions, which is $\equiv 0 \mod p$.



*"Off-diagonal" term: the ones not involving $\chi(a)\chi(b)\chi(c)$. This turns out to be 0. There are two types in this case:

*

*Sums of the form $\sum_{a+b+c = 1} \chi(a)$. Since $\sum_a \chi(a) = 0$, we can sum over $a$ first and get
$$\sum_{a+b+c = 1} \chi(a) = \sum_{a} \chi (a) \sum_{b+c = 1-a} 1 = p^2 \sum_a \chi(a) = 0$$

*Sums of the form $\sum_{a+b+c = 1} \chi(a)\chi(b)$. Similarly,
$$\sum_{a+b+c = 1} \chi(a)\chi(b) = \sum_{a,b} \chi (a) \chi(b) \sum_{c = 1-a-b} 1 = \sum_{a, b} \chi(a)\chi(b) = \left(\sum_a \chi(a)\right)^2 = 0$$



*"Diagonal term": $D_3 = \sum_{a+b+c = 1} \chi(a)\chi(b)\chi(c)$.

*

*This is the most involved case. We will look at the sum
$$D_2 = \sum_{x+y = 1} \chi(x) \chi(y) = -\chi(-1)$$
and
$$E_2 = \sum_{x+y = 0} \chi(x)\chi(y) = (p-1)\chi(-1)$$
When handling $D_3$, sum over $a$ first, so that
$$D_3 = \sum_a \chi(a)\sum_{b+c = 1 - a} \chi(b) \chi(c)$$

*

*When $1 - a = 0$, inner sum is $E_2$.

*When $1 - a \neq 0$, we can write $b = b'(1-a), c = c'(1-a)$, then
$$\chi(b)\chi(c) = \chi(b')\chi(c') \chi(1-a)^2 = \chi(b')\chi(c')$$
since $\chi^2 = 1$ (Legendre symbol), so the inner sum is really $D_2$ in this case.

Hence,
$$D_3 = \chi(1) E_2 + \sum_{a \neq 1} \chi(a) D_2 = E_2 - D_2 = p\chi(-1) \equiv 0 \mod p$$
as well.
