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