The number of solution of $x_1^2 + \cdots + x_k^2 \equiv \lambda \bmod q$

Yes, this is standard. More generally, we have the following

Theorem. Let $p$ be an odd prime, and let $a_1,\dots,a_k\in\mathbb{F}_p^\times$. Then the number of solutions of the equation $a_1x_1^2+\cdots +a_kx_k^2=1$ in $\mathbb{F}_p$ equals \begin{align*} p^{k-1}-\left(\frac{a_1\dots a_k}{p}\right)p^{\frac{k-2}{2}},&\qquad k\equiv 0\pmod{4};\\ p^{k-1}+\left(\frac{a_1\dots a_k}{p}\right)p^{\frac{k-1}{2}},&\qquad k\equiv 1\pmod{4};\\ p^{k-1}-\left(\frac{-a_1\dots a_k}{p}\right)p^{\frac{k-2}{2}},&\qquad k\equiv 2\pmod{4};\\ p^{k-1}+\left(\frac{-a_1\dots a_k}{p}\right)p^{\frac{k-1}{2}},&\qquad k\equiv 3\pmod{4}. \end{align*}

P.S. I am sure there is a reference for this, but I found it easier to copy from my notes.

Added 1. As Gerry Myerson pointed out (in a comment he deleted), the special case $a_1=\dots=a_k=1$ is Proposition 8.6.1 in Ireland-Rosen: A classical introduction to modern number theory (2nd edition).

Added 2. For the sake of completeness, I provide the proof. Let $e_p:\mathbb{F}_p\to\mathbb{C}^\times$ be the standard additive character, $\chi:\mathbb{F}_p^\times\to\{\pm 1\}$ the nontrivial quadratic character, and $$\lambda:=\begin{cases} 1,&\qquad p\equiv 1\pmod{4};\\ i,&\qquad p\equiv 3\pmod{4}.\end{cases}$$ Let us remark that $\lambda^2=\chi(-1)$. If $n$ denotes the number of solutions in the theorem, then \begin{align*} pn&=\sum_{x_1,\dots,x_k\in\mathbb{F}_p}\ \sum_{m\in\mathbb{F}_p}\ e_p(m(a_1x_1^2+\cdots +a_kx_k^2-1))\\ &=\sum_{m\in\mathbb{F}_p}e_p(-m)\sum_{x_1\in\mathbb{F}_p}e_p(ma_1x_1^2)\ \cdots\sum_{x_k\in\mathbb{F}_p}e_p(ma_kx_k^2)\\ &=p^k + \sum_{m\in\mathbb{F}_p^\times}e_p(-m)\bigl\{\chi(ma_1)\lambda\sqrt{p}\bigr\}\cdots \bigl\{\chi(ma_k)\lambda\sqrt{p}\bigr\}\\ &=p^k + \chi(a_1\dots a_k)\lambda^k p^{\frac{k}{2}}\sum_{m\in\mathbb{F}_p^\times}e_p(-m)\chi(m)^k. \end{align*} The inner sum equals $-1$ or $\chi(-1)\lambda\sqrt{p}$ depending on whether $k$ is even or odd, therefore $$n=\begin{cases} p^{k-1}-\chi((-1)^{\frac{k}{2}}a_1\dots a_k)p^{\frac{k-2}{2}},&k\equiv 0\pmod{2};\\ p^{k-1}+\chi((-1)^{\frac{k-1}{2}}a_1\dots a_k)p^{\frac{k-1}{2}},&k\equiv 1\pmod{2}. \end{cases}$$

I just want to share a character-free proof - very accessible and too long for a comment (also - probably not original). Of course, it's just a restatement of a character-full proof.

Let $\vec{v} \in \mathbb{R}^{\mathbb{F}_q}$ be the vector whose $i$'th coordinate is the number of solutions to $x^2 \equiv i \mod q$. Then the distribution of $\sum_{i=1}^{k} x_i^2$ is given by the convolution (coming from the additive structure of $\mathbb{F}_q$) of $\vec{v}$ with itself $k$ times, i.e. $\underbrace{\vec{v} * \cdots * \vec{v}}_{k \text{ times}}$.

The main observation is that the convolution of $\vec{v}$ with itself is almost constant, in the following sense: there are only 2 (possibly equal) values appearing in it, one corresponding to the number of solutions to $x_1^2+x_2^2=0 \mod q$ and the other for all the rest:

• If $q \equiv 1 \mod 4$, there is a root of $-1$ in $\mathbb{F}_q$, denote it by $i$. The equation $x_1^2+x_2^2 =a\mod q$ factors as $(x_1+ix_2)(x_1-ix_2)=a$, which is equivalent (after an $\mathbb{F}_q$-linear transformation) to $y_1 y_2 = a \mod q$. In this form the result is evident.
• If $q \equiv -1 \mod 4$, $i:=\sqrt{-1}$ still exists but in $\mathbb{F}_{q^2}$. The equation still factors as $(x_1 +ix_2)(x_1-ix_2) = a$ in $\mathbb{F}_{q^2}$. Raising to $q$'th power is a (Frobenius) automorphism acting as $(x_1+ix_2)^q=x_1-ix_2$, so we can rewrite this as $z^{q+1} = a$ where $z \in \mathbb{F}_{q^2}$. Because of the cyclic structure of $\mathbb{F}_{q^2}^{\times}$ and because $\mathbb{F}_{q}^{\times}$ is a subgroup of index $q+1$, the number of solutions is $q+1$ when $0 \neq a \in \mathbb{F}_{q}$ (there are $q+1$ roots of unity of order dividing $q+1$ in $\mathbb{F}_{q^2}$).

We can summarize both cases as $(\vec{v}*\vec{v})_a = \begin{cases} q+(\frac{-1}{q})(q-1)& a =0 \\ q-(\frac{-1}{q}) & a\neq 0 \end{cases}$.

The second obserivation is that if $\{ \vec{u_i}\}_{i=1}^{2}$ are two vectors of the form $\vec{v_i} = (a_i, \underbrace{b_i, \cdots, b_i}_{q-1 \text{ times}})$, then their convolution is of the form $(a, \underbrace{b, \cdots, b}_{q-1 \text{ times}})$, where $a,b$ satisfy $a+(q-1)b=(a_1+(q-1)b_1)(a_2+(q-1)b_2)$ and $a-b=(a_1-b_1)(a_2-b_2)$. In other words, the functional $(a, \underbrace{b, \cdots, b}_{q-1 \text{ times}}) \to \begin{cases} a+(q-1)b\\ a-b \end{cases}$ respect convolution and turn it into multiplication.

Now it is only a computational matter: if $k$ is even, we need to convolve $\vec{v} * \vec{v}$ with itself $\frac{k}{2}$ times using the second observation. We get the vector $(a, \underbrace{b, \cdots, b}_{q-1 \text{ times}})$ where $a,b$ are given as a solution of the following linear system: $a+(q-1)b=q^{k/2}, a-b=(q(\frac{-1}{q}))^{k/2}$.

If $k$ is odd, there's an additional complication - we first compute $\underbrace{v * \cdots * v}_{k-1 \text{ times}}=(a, \underbrace{b, \cdots, b}_{q-1 \text{ times}})$, and then we need to convolve it once again with $\vec{v}$ - the $j$'th coordinate will be $av_j + (q-v_j)b= qb + (a-b)(1+(\frac{j}{q}))$.