Given a natural number $n$ and an element $k \in \mathbb{Z}_n$, how many solutions are there in $\mathbb{Z}_n$ to the equation $x^2+y^2 =k$? That is, I'm wondering whether there is a mod-$n$ version of the sum of squares function discussed in this post.
As an example, to write 1 as a sum of squares modulo 7, one has permutations of $0^2+1^2$ but also $2^2+2^2$ and $2^2+5^2$ and $5^2+5^2$.
A formula for the quantity you are considering is given is computed by elementary means in https://arxiv.org/abs/1404.4214 .
In fact, they consider the same counting problem for general quadratic congruences of the form
$$ a_1 x_1^2 + \dots + a_k x_k^2 \equiv b \mod n .$$
Here is a quick way to compute it in the case of an odd sphere (so even number of variables) and $p$ is an odd prime. Note that for a general odd number its easy to reduce to this case by the Chinese reminder theorem and an easy lifting from $p$ to $p^m$. The equation $x^2+y^2=a$ is equivalent to the equation $N_{k[i]/k}(\alpha)=a$ where $k=\mathbb{F}_p$ is the field with $p$-elements. If $p\equiv 3 mod(4)$ then this equation has exactly $p+1$ solutions if $a \ne 0$ and one solution if $a=0$ because the Norm map $Nm:k[i]^\times\to k^\times$ is a surjective homomorphism.
If we denote by $f_\ell(a)$ the function that count the solution with $\ell$ variables, then byt clearly we have
$$f_{2\ell}(a)=\sum_{a_1+...+a_\ell = a}f_2(a_1)f_2(a_2)...f_2(a_\ell)$$, because we can count conditionally on the value of $x_{2i-1}^2+x_{2i}^2$ for every $1\le i\le \ell$ and sum the results. In other words, $f_{2\ell}$ is the $\ell$-fold convolution $f_2*f_2*...*f_2$. Let $I$ denote the constant function with value 1 on $k$, and $\delta_0$ the function which is $1$ at $0$ and $0$ otherwise. The calculation we did for $f_2$ gives exactly that $f_2=(p+1)I -p\delta_0$. Using the linearity of the convolution and the convolutions $I*I=pI$, $\delta_0*f=f$ for every $f$ we get that $$((p+1)I - p\delta_0)^{*\ell}=\sum_{k=0}^{\ell-1} {\ell \choose k} (-1)^{\ell - k}(p+1)^kp^{k-1}p^{\ell-k}I+(-1)^\ell p^\ell \delta_0=$$
$$=p^{\ell-1} \sum_{k=1}^{\ell} {\ell \choose k} (-1)^{\ell - k}(p+1)^kI+(-1)^\ell p^\ell\delta_0=p^{\ell-1}(p^\ell-(-1)^\ell)I+(-1)^\ell p^\ell \delta_0$$
and this (if I don't have any computation mistake!) gives you the answer in the case $p\equiv3 mod 4$. A similar argument works for $p\equiv 1 mod 4$ with different $f_2$ (here the algebra $k[i]$ splits into $k\oplus k$)
For even spheres, its slightly more complicated since eventually you have to convolve with the function $f_1$ which is the sum of $I$ and a quadratic character, but this is not a big difference when you notice that if $\chi$ is such a character then $\chi*I=0$ and $(I+\chi)*(I+\chi)=I+2\chi+\chi*\chi=f_2$ so you have all the convolution formulas you need.
