Using the Sum of Two Squares theorem and Dirichlet Theorem to solve $x^2 + y^2 = k$ for $x,y,k\in \mathbb{Z}/p\mathbb{Z}^*$. Some days ago I asked this question and managed to come up with an answer that determines the existence of a solution to the equation
$$x^2 + y^2 = k$$
where $x,y$ are non-zero elements of $\mathbb{Z}/p\mathbb{Z}$. In fact, such solution exists for all $p>5$ and $k\neq 0$.

What I ask here is to find what the following alternative argument needs for it to be a complete one.
The Sum of Two Squares Theorem tells us that all the primes of the form $4r+1$ can be expressed as the sum of two perfect squares.
So if we find a prime $q$ of the form $4r+1$ such that $q\equiv k \pmod p$ then there are integers $X$ and $Y$ such that $X^2 + Y^2 = q$ and if $x,y$ are the remainders of $X$ and $Y$ modulo $p$, respectively, then $x^2 + y^2 \equiv k \pmod p$.
Now this would be a solution to our original problem as long as $p$ doesn't divide $X$ nor $Y$.
Such prime is always guaranteed to exist due to Dirichlet's Theorem. Indeed, we need to find a prime $q$ such that
$$q\equiv 1 \pmod 4 \\ q \equiv k \pmod p.$$
For $p\neq 2$ then $4$ and $p$ are coprime and therefore there are integers $m,n$ such that
$$4m + pn = 1$$
due to Bézout's Identity. Note that this implies that $(m,p) = (n,4) = 1$. So any prime of the form
$$pn + 4mk + 4pt$$
for some natural number $t$ would solve the above congruences. Since $(4p, pn+4mk) = 1$ then Dirichlet's Theorem assures us that there are infinitely many primes of that form.
Now, if $k$ happens to be a quadratic non-residue then we can take any prime $q$ of that form, find the corresponding $X$ and $Y$ and project them onto $\mathbb{Z}/p\mathbb{Z}$ for a nice solution because we are guaranteed that $x$ and $y$ are non-zero.
However, we need some additional argument to prove that some primes will yield a solution where $p$ doesn't divide $X$ nor $Y$. This argument would have to fail for $p = 5$ since no such solutions exist for that case, but should succeed for any prime $p > 5$.
 A: I'm going to try and answer in a way that simplifies the steps in the result as much as possible. But I'll start with one statement as given: that the Legendre symbol is completely multiplicative, i.e. a product of a residue and a nonresidue is a nonresidue, etc. 
Fact: Let $p$ be any prime $>5$. Then, there exists at least one pair of consecutive residues $\operatorname{mod} p$.
The proof below was borrowed from the answer at the bottom here, by Jack D'Aurizio. His second answer on that page, where he shows to derive an explicit formula for the number of consecutive residue pairs may also provide some insight.
Proof: Since $p$ is a prime $>5$, $2$ and $5$ are nonzero elements $\operatorname{mod} p$. If either $2$ or $5$ is a residue $\operatorname{mod} p$, we are done, because then one of $(1,2)=(1^2,2)$ and $(4,5)=(2^2,5)$ will be our desired pair of consecutive residues. If both $2$ and $5$ are nonresidues, then, since the Legendre symbol is multiplicative, $10=2\cdot 5$ will be a residue $\operatorname{mod} p$. But in that case, $(9,10)=(3^2,10)$ will be a consecutive residue pair. So we are done.
Ok, so we have the existence of consecutive residues $\operatorname{mod} p$. How does that help with the problem at hand? 
So, say we're given some prime $p$ that's $>5$. We know that there'll exist some residues of the form $w$ and $w+1$. Since $w$ is a residue, $w\equiv z^2\operatorname{mod} p$ for some nonzero $z$. But this means $w+1\equiv z^2+1^2\operatorname{mod} p$ and $w+1$ is a residue. So at the very least one residue, $w+1$, can be expressed, via congruence, as the sum of two nonzero squares. 
How can we relate this to the expression of our arbitrary residue, $k$, as a sum of two squares? Notice that, since $k$ and $w+1$ are both residues, $k$ and $(w+1)^{-1}$ are a pair of residues ($(w+1)^{-1}$ denotes the multiplicative inverse of $w+1\operatorname{mod} p$). Then (by the property of the Legendre symbol) $k(w+1)^{-1}$ is a residue $\operatorname{mod} p$, i.e. it is congruent to some $j^2$ (for some nonzero $j$). This gives us $k\equiv j^2(w+1)\operatorname{mod} p$.
Then, by the facts already shown, we have $k\equiv (jz)^2+j^2 \operatorname{mod} p$ making it possible to represent $k$ as the sum of $2$ nonzero squares $\operatorname{mod} p$.  
At this point, we are done. If you observe carefully what we have done above, you'll notice that, fixing $j$, we can convert any expression for $k$ as a sum of two squares to an expression of $w+1$ as a sum of two squares and vice versa. That $w+1$ is part of a consecutive pair of residues is irrelevant to this fact. It is also clear, by the Legendre property, that the same would hold if we replace both $k$ and $w+1$ by a pair of nonresidues. Lastly, it is independent of the exact value of $p$ (whether or not it be $>$ or $\leq 5$). This brings us to a bonus fact: 
Bonus Fact: If $p$ is any prime and $x$ and $y$ are any pair of residues or nonresidues $\operatorname{mod} p$, the number of pairs $(a,b)\in \mathbb Z_p\times \mathbb Z_p$ for which $a^2+b^2\equiv x$ is equal to the number of pairs $(c,d)\in \mathbb Z_p\times \mathbb Z_p$ for which $c^2+d^2\equiv y$.  
This can be used in some counting arguments regarding the number of representations of numbers as sums of two squares $\operatorname{mod}$ a prime.  
EDIT: As for the use of Dirichlet's theorem in this problem. Its (direct) use is only in proving the existence of a sum-of-squares decomposition for nonresidues (as such a decomposition trivially exists for any residue).  For nonresidues, the accepted answer under this post does the job quite nicely.
A: One shouldn't work too hard here. First, if there is a square root of $-1$ in the finite field, $x^2+y^2$ factors and there are non-trivial solutions. If there is no square root of $-1$ in the finite field, then there is one in the unique quadratic extension. By Lagrange's theorem, the Galois norm from $\mathbb F_{p^2}$ to $\mathbb F_p$ is surjective, so, again, the equation has solutions...
