Conceptual explanation for extra/missing $p$ solutions to $x^2+y^2=a \pmod p$ at $a=0$

Question

Throughout, $p$ will denote a prime integer, and $k$ an arbitrary integer.

I have worked through V. Lebesgue's proof of quadratic reciprocity outlined by Keith Conrad in this MO thread, and I feel that I have gained a conceptual understanding of it at all points except for the following "numerical coincidence": it is not difficult to see that the mod-$p$ norm map $\mathbb Z/p\mathbb Z[i] \to \mathbb Z/p\mathbb Z$ restricts to a (multiplicative) group homomorphism on units $(\mathbb Z/p\mathbb Z[i])^\times \to (\mathbb Z/p\mathbb Z)^\times$ and restricts to the $0$-map outside the units; in particular we have that the fibers of the group homomorphism are all the same size ($\implies p-1 \,\big|\, |(\mathbb Z/p\mathbb Z[i])^\times|$), namely

• if $p=4k+1$, then the fibers have size $p-1$ (and the number of non-units in $\mathbb Z/p\mathbb Z[i] \equiv$ solutions to $x^2+y^2=0\pmod p$ equals $2p-1$)
• and if $p=4k+3$, then the fibers have size $p+1$ (and the number of non-units in $\mathbb Z/p\mathbb Z[i] \equiv$ solutions to $x^2+y^2=0\pmod p$ equals $1$).

(This can be proved by using knowledge of how $p$ factors in $\mathbb Z[i]$ and some basic abstract algebra: in the $p=4k+1$ case, $(p) = (\pi_1)(\bar{\pi_1})$ and so by the Chinese remainder theorem $\mathbb Z[i]/(p)$ --- order $p^2$ --- is a product of 2 isomorphic rings, which can only be $\mathbb Z/p\mathbb Z \times \mathbb Z/p\mathbb Z \implies (p-1)^2$ total units; and the $p=4k+3$ case $\mathbb Z[i]/(p)$ is a/the field of order $p^2 \implies p^2-1$ total units.) The "numerical coincidence" is that in both cases, the countfor the "anomaly" case $a=0$ deviates from the count for the "typical" case $a\neq 0$ by $p$.

Goal: I'm looking for a conceptual explanation for this phenomenon ("extra/missing $p$ solutions to the equation $(*): x^2+y^2=a \pmod p$ at $a=0$"), which I'll break down into 3 pieces:

• why the deviation is of the same magnitude in both cases,
• why the $p=4k+1$ deviation is positive and the $p=4k+3$ deviation is negative,
• and why the deviation magnitude is exactly $p$.

I have a somewhat satisfactory half answer for the second item. By the (easy half of the) Fermat Christmas theorem and some knowledge about factorization in the Gaussian integers $\mathbb Z[i]$ (explained conceptually/intuitively/visually/beautifully by Mathologer and 3b1b), we know that for the $p=4k+3$ case, there are no solutions to $x^2+y^2=a$ for any $a=$ a multiple of $p$ between $[p,\frac 12 p^2]$, and so the "anomaly" case $a=0$ of $(*)$ only has 1 solution $(x,y)=(0,0)$. Compared to the average (over $a=0,\ldots, p-1$) number of solutions $\frac{p^2}{p}=p$, this is certainly an anomaly in the negative direction.

Whereas the Fermat Christmas theorem tells us that if $p=4k+1$, we do have solutions to $x^2+y^2=p$. Furthermore, the solutions to $x^2+y^2=a$ for $a=$ a multiple of $p$ in $[p, \frac 12 p^2]$ must come in groups of 8 $(\pm x,\pm y), (\pm y,\pm x)$ (i.e. "diagonal" points $(x,x)$ can't possibly be solutions because $2x^2$ for $x\in \{0, \ldots, \frac{p-1}2\}$ can't possibly equal something times a multiple of $p$, by comparing factorizations into primes in $\mathbb Z$; and the 8 points are seen to be distinct by picturing $(x,y)$ as in the lower left square $\{0,\ldots, \frac{p-1}2\}^2$ under the diagonal $y=x$, and reflecting across the lines $x,y= \frac{p-1}2 + \frac 12$ and $y=x$ to get ther other 7), so we have that $|(\mathbb Z/p\mathbb Z[i]\setminus \{0\}) \smallsetminus (\mathbb Z/p\mathbb Z[i])^\times| = p^2-1-m(p-1)$ for some integer $m \in [0,p+1)$ must be a multiple of $8$. So if $p=8k+5$, $m$ can't possibly be $p$, so it must be $\leq p-1$, so that gives us an anomoly in the positive direction.

---

Maybe it's not so clear why I think my above attempted explanation is more "conceptual". One reason is that the Fermat Christmas theorem dichotomy between $p=4k+1$ vs. $p=4k+3$ ultimately leads to the similar dichotomy in quadratic reciprocity, so I'm happy to use the Fermat result a lot (and I think it makes it intuitive that the $p=4k+3$ case has very few solutions for $(*), a=0$). Also, in the $p=4k+1$ case I exploited some symmetry and made a connection to something people immediately notice upon seeing plots of these solutions, namely the 8-fold symmetry (and I think visual explanations are very valuable/"conceptual").

---

(Crossposted from MSE, but maybe experts here will know some "deeper" facts/theory that explain/are related to this phenomenon --- maybe some algebraic geometry, or something something Weil conjectures, etc.)

Answer 1

I'll begin with a disclaimer that this response doesn't actually resolve the numerical coincidence; I don't provide any direct argument that points gained must equal points lost. But I do show a direct conceptual link between the solutions to $x^2+y^2=a$ and the solutions to $x^2+y^2=0$, which might still be helpful.

Before I get to that I'll offer a direct explanation for the number of solutions to $x^2+y^2=0$ (as opposed to the indirect method of counting and excluding all the solutions to $x^2+y^2\neq 0$). Because the left-hand side can be factored, the equation defines a union of two lines, $$x+iy=0\qquad \text{and} \qquad x-iy=0,$$ where $i$ represents a square root of $-1$ in the algebraic closure of $\mathbb{F}_p$. If $i\in\mathbb{F}_p$ (i.e. if $p=4k+1$) then the number of solutions is $2p-1$: $p$ from each line, minus one for the double-counted intersection point. Otherwise, there is just $(0,0)$ because if $y\in\mathbb{F}_p$ is nonzero then $x=\pm iy\notin \mathbb{F}_p$. This already gives an explanation for why the deviation should be approximately $p$ in each direction: because you're going from one curve (with around $p$ points) to two lines (each with around $p$ points, unless their slopes aren't in $\mathbb{F}_p$).

---

To see how this picture ties in with the solutions of $x^2+y^2=a$, it helps to move to projective coordinates. Let $a\neq 0$ and consider the equation $$X^2+Y^2=aZ^2.$$ Given any point $(X,Y,Z)\in \mathbb{F}_p^3\setminus\{(0,0,0)\}$ satisfying this equation, we can scale it by any $\lambda\in\mathbb{F}_p^\times$ and get another solution to the equation; thus the nonzero solutions come in equivalence classes with $p-1$ points each. There is almost a bijection between solutions $(x,y)$ to $x^2+y^2=a$ (the inhomogeneous equation) and equivalence classes of solutions $(X,Y,Z)$ to $X^2+Y^2=aZ^2$ (the homogeneous equation). Each inhomoegenous solution $(x,y)$ yields a solution $(x,y,1)$ to the homogeneous equation, and each homogeneous solution $(X,Y,Z)$ with $Z$ nonzero determines a solution $(\frac{X}{Z},\frac{Y}{Z})$ of the inhomogeneous equation, where equivalent $(X,Y,Z)$ give the same $(x,y)$.

A key observation is that regardless for any odd prime $p$ and any nonzero $a\in\mathbb{F}_p$, the equation $X^2+Y^2=aZ^2$ always has exactly $p+1$ equivalence classes of solutions (for now I'll just say "proof by stereographic projection," but can provide a more complete explanation if needed). So the difference in the number of solutions to $x^2+y^2=a$ depending on $p\mod 4$ comes from considering the solutions with $Z=0$.

• If $\mathbb{F}_p$ has no square root of $-1$, then $X^2+Y^2=0$ has no nonzero solutions. Therefore all $p+1$ of the equivalence classes have $Z\neq 0$ and so give inhomogeneous solutions.
• If $\mathbb{F}_p$ has a square root $i$ of $-1$, then there are two equivalence classes of solutions with $Z=0$: $(\lambda, i\lambda,0)$ and $(\lambda, -i\lambda, 0)$. Therefore only $p-1$ of the equivalence classes give inhomogeneous solutions, and the remaining two equivalence classes (with $p-1$ points each) provide nontrivial solutions to $X^2+Y^2=0$.

In other words, think of $p=4k+3$ as the default. When we switch to $p=4k-1$, we get an "extra" $2(p-1)$ solutions to $x^2+y^2=0$ (from $1$ to $2p-1$), and we "lose" two solutions of $x^2+y^2=a$ (from $p+1$ to $p-1$). This is no coincidence: the $2(p-1)$ extra solutions to $x^2+y^2=0$ are exactly the two equivalence classes of homogeneous solutions that fail to give inhomogeneous solutions.

In light of this, the numerical coincidence follows from the "conversion rate" between equivalence classes and individual points. If you add the two point difference in $x^2+y^2=a$ to the $2(p-1)$ point difference in $x^2+y^2=0$, you get $2p$, and so $$(p+1)-1=p\qquad\text{implies }\qquad ((p+1)-2)-(1+2(p-1))=-p.$$