It is easy to prove that for prime $p$ every element of $\mathbb{Z} / p \mathbb{Z}$ can be written as the sum of two squares. An elementary proof is given here: Sum of two squares modulo p
How can we show that, provided further $p \geq 7$, any nonzero element of $\mathbb{Z}/ p \mathbb{Z}$ is the sum of two nonzero squares? I don't see how we could extend the counting argument used in the linked post to this case. Thanks
Adapting Mikhail Ivanov's argument to a slightly different but AFAICT equivalent question to fit here. Some of the elements appeared also in my answer to that question.
Every non-zero element of $\Bbb{Z}/p\Bbb{Z}$ is either a square or a non-square If $a=b^2$ is a non-zero square, then, as $p>5$ we have $$ a=b^2=(3b/5)^2+(4b/5)^2 $$ as a sum of two non-zero squares.
On the other hand, if $a$ is a non-zero non-square then $ab^2$ is a non-zero non-square for any $b\neq0$. Furthermore, we get all the non-squares in this way. If $a=x^2+y^2$ with $xy\neq0$, then $ab^2=(bx)^2+(by)^2$, so it suffices to show that we can write at least one of the non-squares as a sum of two non-zero squares. Let's pretend for one time's sake that $\Bbb{Z}/p\Bbb{Z}$ has an "order", so $0<1<2<\ldots<p-1$. Let $a$ be the smallest non-square in this order. Clearly $a>1$. It follows that $a-1=b^2$ is a square, where $b\neq0$. This implies that $a=1+b^2$ is a sum of two non-zero squares. Therefore so are all the other non-squares.
Let's prove that $1$ is the sum of two nonzero squares.
Think about the geometry of the plane. Any non-vertical line through $(1,0)$ will meet the unit circle in another point. This argument works over $k=\Bbb Z/p\Bbb Z$ too. The "non-vertical" lines through $(1,0)$ are $y=m(x-1)$, and each meets the unit circle in the further point $((m^2-1)/(m^2+1),2m/(m^2+1))$, as long as $m^2\ne-1$. These coordinates are nonzero iff $m\notin\{0,1,-1\}$. There are at most five "bad" values of $m$, so if $p>5$, then $1$ is the sum of two non-zero squares in $k$.
By scaling, each non-zero square is the sum of two non-zero squares in $k$. Of course, if $b$ is not a square, it is still the sum of two squares, and they are both nonzero.
This is not a complete answer, since a crucial part is not proved.
Conjecture: Given an odd prime $p$, the cardinal of the set $$A_p=\{(x,y):0\le x\le(p-1)/2,\, 0\le y\le(p-1)/2,\, x^2+y^2\equiv 1\!\!\!\!\!\mod p\}$$ is $$\begin{cases}(p+3)/4\text{ if } p\equiv 1\pmod 4\\(p+5)/4\text{ if }p\equiv 3\pmod 4\end{cases}$$
This conjecture is based in evidence for odd primes less than $2000$, that I have tried using a program.
Now, assuming this conjecture, let's prove the statement. Take $a\in\Bbb Z_p^\times$ and suppose that $a$ is not the sum of two nonzero squares. Since $a$ is the sum of two squares, we see that $a$ must be a square itself, say $a=b^2$.
The conjecture implies that the number of elements in $A_p$ is at least $3$ for $p\ge 7$. Then, $A_p$ contains an element other than $(1,0)$ and $(0,1)$. This shows that $1$ is the sum of two nonzero squares. That is, $$x^2+y^2\equiv 1\pmod p$$ And, now, multiplying by $b^2$ (which is not $0$), we get $$(bx)^2+(by)^2\equiv a\pmod p$$
Lemma:
Let $7 \leq p$ be a prime number;
then $1$ can be written as the sum of two non-zero squares.
Proof:
See the answer of Lord Shark the Unknown to this question.
Let $a \in \left(\dfrac{\mathbb{Z}}{\mathbb{pZ}}\right)^*$ be arbitrary; let's define the sets $A_a$ and $B$ as follows:
$$ B:= \Bigg{\{} x^2 : x \in \dfrac{\mathbb{Z}}{\mathbb{pZ}} \Bigg{\}} ; \\ A_a:= a-B= \Bigg{\{} a-x^2 : x \in \dfrac{\mathbb{Z}}{\mathbb{pZ}} \Bigg{\}} ; $$
notice that, each of the above sets has exactly $\dfrac{p+1}{2}$ elements.
Claim: $B \cap A_a \neq \phi$,
in other words there is $t \in \dfrac{\mathbb{Z}}{\mathbb{pZ}} $
such that $t^2 \in B$ and $t^2 \in A_a$.
Prood: Suppose on contrary $B \cap A_a = \phi$;
so we can conclude that:
$\Big{|}B \cup A_a \Big{|} = \Big{|} B \Big{|} + \Big{|} A_a \Big{|}$;
on the other hand we know
$B \cup A_a \subseteq \dfrac{\mathbb{Z}}{\mathbb{pZ}}$,
so we can conclude that:
$$ p+1 = \dfrac{p+1}{2} + \dfrac{p+1}{2} = \Big{|} B \Big{|} + \Big{|} A_a \Big{|} = \Big{|} B \cup A_a \Big{|} \leq \Big{|} \dfrac{\mathbb{Z}}{\mathbb{pZ}} \Big{|} = p ; $$
which is an obvious contradiction.
Remark :
Let $7 \leq p$ be a prime number,
then every nonzero element
$a \in \left(\dfrac{\mathbb{Z}}{\mathbb{pZ}}\right)^*$
can be written as the sum of two non-zero squares.
Proof :
By the above claim,
there is $x \in \dfrac{\mathbb{Z}}{\mathbb{pZ}} $
such that $x^2 \in A_a$,
so there is $y \in \dfrac{\mathbb{Z}}{\mathbb{pZ}} $ such that $x^2=a-y^2$;
which proves the remark.
Proof of the original problem:
By the previous remark every $a \in \left(\dfrac{\mathbb{Z}}{\mathbb{pZ}}\right)^*$ can be written in the form $x^2+y^2$; notice that $(x,y)\neq (0,0)$, so we have the following cases:
$x \neq 0$ and $y \neq 0$; in this case there is nothing left to prove. $\checkmark$
$x \neq 0$ and $y = 0$ . By the lemma, $1$ can be written in the form $a^2+b^2$; where both of $c, d$ are non-zero, only notice that: $$ a \overset{p}{\equiv} 1 \cdot a \overset{p}{\equiv} (c^2+d^2) (x^2+0^2) \overset{p}{\equiv} (cx)^2+(dx)^2 . $$
$x = 0$ and $y \neq 0$ . By the lemma, $1$ can be written in the form $a^2+b^2$; where both of $c, d$ are non-zero, only notice that: $$ a \overset{p}{\equiv} 1 \cdot a \overset{p}{\equiv} (c^2+d^2) (0^2+y^2) \overset{p}{\equiv} (cy)^2+(dy)^2 . $$
Since $\mathbf F_p :=\mathbf Z /p\mathbf Z$ is a field, it seems adequate to heartily use the theory of finite fields, more specifically the central result that an extension $\mathbf F_q /\mathbf F_p$, with $q=p^n$, is Galois with cyclic group of order $n$ generated by the Frobenius automorphism $\phi: x \to x^p$. It follows easily that the norm map $N:\mathbf F_q^{*} \to \mathbf F_p^{*}$ is surjective. For later use, I include a quick proof : by Hilbert's 90, $x\in KerN$ iff $x$ is of the form $\phi(z)/z = z^{p-1}$; in particular, since $\mathbf F_q^{*}$ is cyclic, $KerN$ has order $(q-1)/(p-1)$, hence $ImN$ has order $p-1$, i.e. $N$ is surjective.
Apply this for $n=2$ : if $p\neq 2$, then $\mathbf F_{p^2} = \mathbf F_p (\sqrt d)$ where $d$ is the class of a non quadratic residue mod $p$, hence every $x\in \mathbf F_p^{*}$ can be written $x=a^2 -db^2$ = a sum of two squares. But one of these squares could still be zero, so it remains to prove that every element of $(\mathbf F_p^{*})^2$ is a sum of two nonzero squares, or equivalently that $1=x^2+y^2, xy\neq0$. The latter property has been shown in the answers given by @Jyrki Lahtonen and @Lord Shark the Unknown. Although it's not necessary to give all the solutions $(x, y)$, let's provide a parametrization analogous to that of the pythagorean triples ((3, 4, 5) is the "first" such triple). Put $1=a^2-db^2$= a norm from $\mathbf F_p (\sqrt d)$ with $d$ as above. The same argument using Hilbert's 90 will give $a+b\sqrt d= \phi(z)/z$, and writing $z=e+f\sqrt d$, we'll get $a=(e^2+df^2)/(e^2-df^2)$ and $y=-2ef/(e^2-df^2)$ . For example, if $d=-1$, we recover $(\frac 35)^2 + (\frac 45)^2 =1$ by taking $e=2, f=-1$.
