Number of solutions to Markov Diophantine equation mod $p$ I am interested in the number of solutions $(x,y,z)\pmod{p}$ (with $p\ne2,3$) to
$$x^2+y^2+z^2-3xyz\equiv0\pmod{p}$$
See Markov numbers. The solutions are ordered triples and do not include $(0,0,0)$. Let $S(p)$ denote the number of such solution triples. Here are a few values of $S$:
$$S(5)=40,\quad S(7)=28,\quad S(11)=88,\quad S(13)=208$$
It is easy to show that $4\mid S(p)$. Can we prove that $p\mid S(p)$?
Edit: W-t-p has provided an explicit form for $S(p)$ which proves the divisibility claim. The following problem remains:
For $d$ not necessarily prime, if $S(d)$ includes $(0,0,0)$ as a solution, then by the Chinese Remainder Theorem, $S(d)$ is multiplicative: if $p,q$ are coprime, then
$$S(pq)=S(p)S(q)$$
So: is there a closed form for $S(p^r)$? Here are a few values of $S$ including $(0,0,0)$:
$$S(2)=5,S(4)=24,S(8)=96,S(16)=320$$
$$S(3)=9,S(9)=99,S(27)=891$$
$$S(5)=41,S(25)=1280$$
 A: $\newcommand{\L}[1]{\left(\frac{#1}p\right)}$
There is, in fact, a simple and clean explicit formula, which agrees with the numerical data and which readily gives your divisibility claim:
  $$ S(p) = p^2+3\L{-1}p. $$
Here the proof goes.
For any fixed $x,y\in\mathbb F_p$, we can consider $x^2+y^2+z^2=3xyz$ as a quadratic equation in the variable $z$ with the discriminant $9x^2y^2-4x^2-4y^2$, which therefore has $\L{9x^2y^2-4x^2-4y^2}+1$ solutions. It follows that the total number of triples $(x,y,z)$, including the zero triple, satisfying the equation, is
\begin{align*}
  S(p)+1 &= \sum_{x,y} \left(\L{9x^2y^2-4x^2-4y^2} + 1\right) \\
         &= p^2 + 2\sum_x \L{-4x^2} + \sum_{x,y\ne 0} \L{9x^2y^2-4x^2-4y^2} \\
    &= p^2 + 2\L{-1}(p-1) + \sigma,
\end{align*}
where
  $$ \sigma = \sum_{x,y\ne 0} \L{9x^2y^2-4x^2-4y^2} 
     = \sum_{a,b\ne 0} \L{1-a^2-b^2} $$
(for the last equality, we passed to the new variables $a=2/(3x)$ and $b=2/(3y)$). The sum in the right-hand side consists of summand of the form $\L{1-u-v}$, each summand appearing four times if both $u$ and $v$ are quadratic residues, twice if one of $u$ and $v$ is equal to $0$ while another one is a quadratic residue, once if $u=v=0$, and does not appear at all if at least one of $u$ and $v$ is a quadratic non-residue; that is, $\L{1-u-v}$ appears exactly $\left(\L{u}+1\right)\left(\L{v}+1\right)$ times. As a result,
  $$ \sigma = \sum_{u,v\ne 0} \left(\L{u}+1\right) \left(\L{v}+1\right) \L{1-u-v} = \sigma_1 + \sigma_2 + \sigma_3, $$
with
\begin{align*} 
  \sigma_1 &= \sum_{u,v\ne 0} \L{1-u-v}, \\
  \sigma_2 &= 2  \sum_{u,v\ne 0} \L{u}\L{1-u-v}, \\
  \sigma_3 &= \sum_{u,v\ne 0} \L{u}\L{v}\L{1-u-v}.   
\end{align*}
Furthermore,
\begin{align*} 
  \sigma_1 &= \sum_{u\ne 0} \left(- \L{1-u}\right) = 1, \\
  \sigma_2 &= 2\sum_{u\ne 0} \L{u} \left( -\L{1-u} \right) \\
       &= - 2\L{-1} \sum_{u\ne 0} \L{u(u-1)} \\
       &= - 2\L{-1} \sum_{w\ne 0} \L{1-w} \\
       &= 2\L{-1}, 
\end{align*}
and, finally,
\begin{align*}
  \sigma_3 &= \sum_{u,v\ne 0} \L{uv(u+v+1)}  \\
       &= \sum_{u,w\ne 0} \L{w(u+wu+1)} \\
       &= \sum_{w\ne 0} \L{w} \sum_{u\ne 0} \L{(w+1)u+1} \\
       &= \L{-1}(p-1) - \sum_{w\ne 0,-1} \L{w} \\
       &= \L{-1} \left(p-1+1\right) \\
       &= \L{-1} p.
\end{align*}
Assembling all pieces together, 
  $$ \sigma = 1 + 2\L{-1} + \L{-1}p  = \L{-1}(p+2)+1, $$
and then
  $$ S(p)+1 = p^2 + 2\L{-1}(p-1) + \L{-1}(p+2)+1 = p^2 + 3\L{-1}p + 1. $$

Concerning the value of $S(p^n)$ for $n>1$, it seems that Hensel's Lemma does the job; alternatively, just count the number of solutions $\pmod{p^n}$ lying above a given solution $\pmod{p^{n-1}}$. For instance, the trivial solution $(0,0,0)$ modulo $p$ gives raise to $p^3$ solutions modulo $p^2$ (which are all the triples $(x,y,z)$ in $\mathbb Z/p^2\mathbb Z$ with $x,y$, and $z$ all divisible by $p$). Furthermore, the number of solutions arising from a non-trivial solution $(x,y,z)$ modulo $p$ is the number of triples $(a,b,c)\in (\mathbb Z/p\mathbb Z)^3$ such that 
  $$ (x+pa)^2+(y+pb)^2+(z+pc)^2\equiv 3(x+pa)(y+pb)(z+pc)\pmod{p^2}, $$
and it is not difficult to see that there are $p^2$ such triples. This gives
  $$ S(p^2)+1 = p^3 + p^2S(p) = p^3 \left( p+1+3\L{-1} \right),\quad p>3. $$
In particular, $S(25)=1124$, which agrees with the computation I have made.
A: Above equation shown below:
$$x^2+y^2+z^2=3xyz \qquad \text{(1)}$$
Equation $(1)$ has parametric solution shown below:
\begin{align}
x &= w(5k^2-4k+1) \\
y &= w(13k^2-10k+2) \\
z &= 1 \end{align}
where, $w=1/(k^2-3k+1)$.
For, $k=3$ we get $(x,y,z)=(34,89,1)$
For, $k=4$ we get $(x,y,z)=(13,34,1)$
