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.)