Let $p$ be a prime. How many solutions $(x, y)$ are there to the equation $x^2 + 3y^2 \equiv 1 \pmod{p}$? Here $x, y \in \{0, 1, \ldots p-1\}$. This paper (https://arxiv.org/abs/1404.4214) seems like it may give some ideas, but focuses on equations of the form $\sum_{i = 1}^n x_i^2 \equiv 1 \pmod{p}$, i.e. it does not allow for the possibility of coefficients on the variables.

The count is $2$ and $6$ for $p=2$ and $p=3$ respectively, and otherwise $p-1$ or $p+1$ according as $p$ is $1$ or $-1 \bmod 3$.

More generally, it is well-known that a smooth plane conic
over ${\bf Z} / p {\bf Z}$ has $p+1$ points in the *projective* plane,
so we need only subtract the number of points at infinity,
which here is the number of square roots of $-3 \bmod p$.
(The primes $p=2,3$ are special because $x^2+3y^2-1$ factors as
$(x+y-1)^2 \bmod 2$ and $(x+1)(x-1) \bmod 3$.)

For general $x^2 - D y^2 = 1 \bmod p$ (with $p \not\mid 2D$) the count is $p - (D/p)$ where $(D/p)$ is the Legendre symbol. For $D=-3$ we can also choose between $p-1$ and $p+1$ by observing that the solutions come in triples $\{ (x,y), \frac12(-x-3y,x-y), \frac12(-x+3y,-x-y) \}$ so the count must be a multiple of $3$. This trick is possible here because $(-1+\sqrt{-3})/2$ is a cube root of unity; likewise for $D = -1$ there's a fourth root of unity $\sqrt{-1}$, and the solutions come in quadruples $\{(\pm x, \pm y), (\pm y, \mp x)\}$ so for any odd prime $p$ the number of solutions of $x^2 + y^2 \equiv 1 \bmod p$ is whichever of $p \pm 1$ is a multiple of $4$.

You may simply rationally parametrize this curve: if $y=0$ we have 2 solutions, otherwise we may denote $x=1+ty$, and $x^2+3y^2=1$ reads as $2t+(3+t^2)y=0$. If $-3$ is a quadratic non-residue modulo $p\ne 3$, this has unique solution $y=-2t/(3+t^2)$ for each $t\ne 0$. If $-3$ is a quadratic residue, there are two non-zero values of $t$ for which there is no solution.

Here is an alternate approach, more algebraic and less geometric. As in Noam's answer we'll consider the more general equation $x^2 - Dy^2 \equiv 1 \bmod p$. Consider the $\mathbb{F}_p$-algebra $A = \mathbb{F}_p[\sqrt{D}] \cong \mathbb{F}_p[\alpha]/(\alpha^2 - D)$, which is not necessarily a field. It has elements of the form $x + y \sqrt{D}$ and a multiplicative norm map

$$N : A \ni x + y \sqrt{D} \mapsto N(x + y \sqrt{D}) = x^2 - Dy^2 \in \mathbb{F}_p.$$

We want to compute $|N^{-1}(1)|$, which can be done as follows. Things split into a few cases, the last few of which are degenerate. In most of the cases $N$ will turn out to be surjective which gives $|N^{-1}(1)| = \frac{|A^{\times}|}{|\mathbb{F}_p^{\times}|}$ and we will compute this.

**Case 1:** $\gcd(D, p) = 1$ and $\left( \frac{D}{p} \right) = -1$ (which implies $p$ odd). Then $A$ is the finite field $\mathbb{F}_{p^2}$, whose multiplicative group is cyclic of order $p^2 - 1$, and the norm map $N : \mathbb{F}_{p^2} \to \mathbb{F}_p$ is given by $N(\alpha) = \alpha \overline{\alpha} = \alpha^{p+1}$. Hence if $\alpha$ is a generator of the multiplicative group of $\mathbb{F}_{p^2}$ then $\alpha^{p+1}$ has order $p-1$ and so is a generator of the multiplicative group of $\mathbb{F}_p$. It follows that

$$\left| \{ (x, y) \in \mathbb{F}_p : x^2 - Dy^2 = 1 \} \right| = \frac{p^2 - 1}{p - 1} = p + 1.$$

**Case 2:** $\gcd(D, p) = 1$ and $\left( \frac{D}{p} \right) = 1$ and $p$ is odd. Then $\alpha^2 - D = (\alpha - \sqrt{D})(\alpha + \sqrt{D})$ splits into a polynomial with two distinct roots which gives $A \cong \mathbb{F}_p^2$. Under this isomorphism the norm map can be identified with the multiplication map $\mathbb{F}_p^2 \ni (x, y) \mapsto xy \in \mathbb{F}_p$ which is clearly surjective. It follows that

$$\left| \{ (x, y) \in \mathbb{F}_p : x^2 - Dy^2 = 1 \} \right| = \frac{(p - 1)^2}{p - 1} = p - 1.$$

**Case 3a:** $p \mid D$ and $p$ is odd. Then $\alpha^2 - D \equiv \alpha^2 \bmod p$ is a polynomial with repeated roots so $A \cong \mathbb{F}_p[\alpha]/\alpha^2$ and the norm map is $N(x + y \alpha) = x^2$. Here the norm map is not surjective: its image on units is the subgroup of squares which has order $\frac{p-1}{2}$. This gives

$$\left| \{ (x, y) \in \mathbb{F}_p : x^2 = 1 \} \right| = \frac{p(p - 1)}{ \frac{p-1}{2} } = 2p$$

which of course can be obtained in a more elementary way; I'm including this case for the sake of completeness.

**Case 3b:** $p = 2$. Here $\alpha^2 - D$ is either $\alpha^2 \bmod 2$ or $(\alpha - 1)^2 \bmod 2$ so $A \cong \mathbb{F}_2[\varepsilon]/\varepsilon^2$. The norm map is either $N(x + y \sqrt{D}) = x^2 = x$ if $D$ is even or $N(x + y \sqrt{D}) = x^2 - y^2 = x + y$ if $D$ is odd; either way it is clearly surjective and we get

$$\left| \{ (x, y) \in \mathbb{F}_2 : x^2 - Dy^2 = 1 \} \right| = 2$$

which is even easier to obtain but again is being included for the sake of completeness.

For any for any nonzero residues $a,b,c\not\equiv 0\pmod{p}$, the number of solutions of the congruence $ax^2+by^2\equiv c\pmod{p}$ equals $p-\left(\frac{-ab}{p}\right)$. For the proof, see my post here. Note also that this proof generalizes easily to any number of variables.

We give yet another approach to proving the result, using nothing more advanced than properties of the Legendre symbol.

Proposition.Suppose $p$ is an odd prime and $D$ is not a multiple of $p$.Then the number of solutions of $x^2 - Dy^2 \equiv 1 \bmod p$ is$p - (D/p)$, i.e. $p-1$ if $D$ is a square mod $p$ and $p+1$ if not.

*Proof*:
It is clear that the answer depends only on the Legendre symbol $(D/p)$,
because for any nonzero $c \bmod p$ the invertible change of variables
$(x,y) \to (x,cy)$ shows that the count is the same for $D$ and $c^2 \! D$.
Let $N_+$, then, be the number of solutions for $(D/p) = +1$,
and let $N_-$ be the number of solutions for $(D/p) = -1$.

Next observe that $N_+ + N_- = 2p$, because each $x$ that contributes $2$, $1$, or $0$ solutions to $N_+$ contributes $0$, $1$, or $2$ solutions respectively to $N_-$. So it is enough to prove the Proposition for one of $N_+$ and $N_-$.

But $N_+$ is the number of solutions mod $p$ of
$1 \equiv x^2 - y^2 = (x+y) (x-y)$. Write $(r,s) = (x+y, x-y)$,
so $x,y \equiv (r\pm s)/2$. Now the equation
$rs \equiv 1$ has $p-1$ solutions mod $p$:
each nonzero $r \bmod p$ determines uniquely the value of
$s \equiv r^{-1} \bmod p$. Therefore the equation
$x^2 - y^2 \equiv 1 \bmod p$ also has $p-1$ solutions,
namely $\bigl(\frac12(r+r^{-1}), \frac12(r-r^{-1})\bigr)$.
This proves that $N_+ = p-1$,
whence it follows that $N_- = 2p - N_+ = p+1$. **QED**

We give an elementary approach to the more general equation $x^2+Dy^2=1 \mod p$.

Firstly, if $p=2$ we have $2$ solutions as either $x^2=1$ or $x^2+y^2=1$ and

if $p|D$ we have $p$ values for $y$ and $x=\pm1$ giving $2p$ solutions.

The remaining case is $(D,p)=1$ with $p$ odd.

The equation $x^2+Dy^2=z^2$ can be solved by re-writing as $z^2-x^2=Dy^2$ or $AB=Dy^2$ where $z-x=A$, $z+x=B$. An easy calculation shows this has a total of $p^2$ solutions. (Consistent with the Chevalley-Warning theorem which says the number of solutions must be divisible by $p$)

However we need to remove the $z=0$ solutions. We have $2p-1$ to remove if -D is a quadratic residue and 1 if -D is a non-residue giving $(p-1)^2$ and $p^2-1$ solutions respectively.

Since $z\neq0$ has been excluded we can divide by $p-1$ to get the number of solutions to the original equation yielding $p-1$ if $-D$ is a residue and $p+1$ otherwise.