Zeta-function of $x_0^3 + x_1^3 + x_2^3 = 0$ over $\mathbb{F}_4$. Consider the projective curve over the finite field $\mathbb{F}_4$
$$C/\mathbb{F}_4\colon x_0^3 + x_1^3 + x_2^3 = 0.$$
Let's calculate the zeta function
$$Z (C/\mathbb{F}_q) = \exp \left(\sum_{n \ge 1} \frac{\# C (\mathbb{F}_{q^n})\,t^n}{n}\right)$$
(where $q = 4$ in this case).
I see two ways of doing that, but neither is elementary.


*

*The curve is smooth of degree $3$, so its genus is $1$, and in this case the zeta function depends only on $\# C (\mathbb{F}_q)$: it is given by
$$Z (C/\mathbb{F}_q) = \frac{1 - at + qt^2}{(1-t)\,(1-qt)}, \quad a = 1 + q - \# C (\mathbb{F}_q)$$
(see e.g. [Silverman, GTM 106, Section V.2]). It is easy to count that $\# C (\mathbb{F}_4) = 9$, and therefore
$$Z (C/\mathbb{F}_4) = \frac{(1 + 2t)^2}{(1-t)\,(1-4t)}.$$

*Alternatively, for any hypersurface of the form $x_0^m + x_1^m + \cdots + x_n^m = 0$ there is a formula in terms of Gauss sums (see e.g. [Ireland-Rosen, GTM 84, Section 11.3]), which in this case boils down to
$$Z (C/\mathbb{F}_4) = \frac{P (t)}{(1-t)\,(1-4t)},$$
where
$$P (t) = \left(1 + \frac{1}{4}\,g (\chi_1)^3\,t\right)\,\left(1 + \frac{1}{4}\,g (\chi_2)^3\,t\right),$$
and $\chi_1,\chi_2$ are the nontrivial characters $\mathbb{F}_4^\times \to \mathbb{C}$. A little calculation gives $P (t) = (1+2t)^2$, as expected.
However, once we have the zeta function, we see that the number of points is given by a rather easy formula:
\begin{align*}
\# C (\mathbb{F}_4) & = 9 = (2+1)^2,\\
\# C (\mathbb{F}_{4^2}) & = 9 = (2^2-1)^2,\\
\# C (\mathbb{F}_{4^3}) & = 81 = (2^3+1)^2,\\
\# C (\mathbb{F}_{4^4}) & = 225 = (2^4-1)^2,\\
\# C (\mathbb{F}_{4^5}) & = 1089 = (2^5+1)^2,\\
\# C (\mathbb{F}_{4^6}) & = 3969 = (2^6-1)^2,\\
 & \cdots
\end{align*}
So my question is:

could anybody give an elementary counting argument explaining this?

 A: Using some facts about elliptic curves it could go as follows.
Let's look at the cubic projective curve given by the equation
$$
X^3+Y^3+Z^3=0.
$$
Putting this into Weierstrass form is easy. All we need to do is to substitute $Z=Y+U$, and (recalling that we are in characteristic two) we get
$$
X^3+U^3+U^2Y+UY^2=0.
$$
If we dehomogenize $x=X/U,y=Y/U$, we get an equation in the Weierstrass form
$$
E:x^3+1=y^2+y.\qquad(*)
$$
The above process maps $\Bbb{F}_{4^m}$-rational points to each other bijectively, so we might as well count the number of points $\#E(\Bbb{F}_{4^m})$.
The curve $E$ has a very special form. If $P=(x_0,y_0)\in E$, then the negative of this point (= the other point on the same vertical line) is $[-1]P=(x_0,y_0+1)$.
Furthermore, implicit differentiation tells us that the tangent of $E$ at $P$ has slope $x_0^2$. The usual process then leads to the very simple point doubling formula
$$
[2](x_0,y_0)=(x_0^4,y_0^4+1).\qquad(**)
$$
In other words doubling can be achieved by applying the square of the Frobenius followed by negation.
If $(x_0,y_0)\in E(\Bbb{F}_{4^m})$ then repeated application of $(**)$ tells us that
$$
[2^m](x_0,y_0)=(x_0^{4^m},y_0^{4^m}+m\cdot1)=
\begin{cases}P,&\ \text{if $m$ is even, and}\\
[-1]P,&\ \text{if $m$ is odd.}
\end{cases}
$$

This means that all the points $P$ of $E(\Bbb{F}_{4^m})$ satisfy the equation $[M]P=0$ for $M=2^m-(-1)^m$.

With a few basic facts about elliptic curves in place this implies that
$$E(\Bbb{F}_{4^m})=E[M],\qquad(**)$$
and thus has $M^2$ points. More precisely, the highlighted result tells that we have an inclusion "$\subseteq$" in $(**)$, so $\#E(\Bbb{F}_{4^m})$ must be a factor of $M^2$. But then the Hasse-Weil bound $|\#E(\Bbb{F}_{4^m})-(4^m+1)|\le2\cdot 2^m$ rules out all proper factors.

I'm afraid this is unsatisfactory as an answer here. Hasse-Weil bound was needed to conclude, and that is at the same depth as the use of Zeta functions. I posted it chiefly, because I discussed an argument related to elliptic curve in the comments.
The calculation I recalled must have been about another elliptic curve defined over $\Bbb{F}_2$. Sorry. Anyway, this does lead to the conclusion, but overall the method is kludgier than those listed in the question.
A: The idea in my latter comment (under main) leads to a "pure" counting argument. Unfortunately only in the case of an odd value of $m$. I abbreviate $q=2^m$.
Because we think projectively we can normalize the first non-zero homogeneous
coordinate of $[x_0:x_1:x_2]$ to be equal to one.
Let us first look at the points with one of the coordinates equal to zero.
The field $E=\Bbb{F}_{q^2}$ contains a primitive third root of unity $\omega$. If we for example assume that $x_0=0, x_1=1$, we see that $x_2$ must be a power of $\omega$. This shows that the projective curve has exactly nine points with the property that one coordinate vanishes, namely the cyclic shifts of $[0:1:1]$, $[0:1:\omega]$ and
$[0:1:\omega^2]$. For the rest of this post we consider the points with all the 
homogeneous coordinates $\neq0$. In other words, we want to count the number of solutions
of
$$
x_1^3+x_2^3=1\tag{1}
$$
with $x_1,x_2\in E^*$.
The multiplicative group $E^*$ is cyclic of order $q^2-1=(q-1)(q+1)$. As $\gcd(q-1,q+1)=1$ we can split it into a direct product
$$
E^*=F^*\times S,
$$
where $F=\Bbb{F}_{q}$ is the unique subfield such that $[E:F]=2$, and $S$ is the
kernel of the relative norm map $N_{E/F}:x\mapsto x^{q+1}$. 
If we assume that $2\nmid m$, then $3\nmid 2^m-1=q-1$ implying that cubing permutes the elements of $F^*$. Another key observation is that if $\eta_1$ and $\eta_2$ are two distinct elements of $S$, then their ratio is not in $F$, and hence $\{\eta_1,\eta_2\}$
is an $F$-basis of $E$. Observe that $3\mid 2^m+1$ implying that $\omega\in S$.
Anyway, we can uniquely write $x_i=y_i\eta_i$ with $y_i\in F^*$, $\eta_i\in S$. Plugging these in $(1)$ gives the equation $y_1^3\eta_1^3+y_2^3\eta_2^3=1$. Because cubing permutes the set of available $y$-components, and we are only in the interested in the number of solutions, we can equally well replace $y_i^3$ with $y_i$. So we want the number of solutions of
$$
y_1\eta_1^3+y_2\eta_2^3=1.\tag{2}
$$
On the other hand, cubing is a 3-to-1 mapping on $S$. This allows us to count the number of solutions of $(2)$:


*

*Assume first that $\eta_1^3\neq1$. This holds for all choices of $\eta_1$ except powers of $\omega$, so there are $q+1-3=q-2$ choices for $\eta_1$. Assume further that $\eta_2^3\notin\{1,\eta_1^3\}$. To a fixed $\eta_1$ this leaves $q+1-6=q-5$ choices for $\eta_2$. Anyway, in this main case $\eta_1^3$ and $\eta_2^3$ for a basis of $E$ over $F$. Therefore $(2)$ has a unique solution $(y_1,y_2)\in F^2$. Here both $y_1,y_2$ are necessarily non-zero. For if, say $y_1=0$, then $(2)$ implies that $y_2\eta_2^3=1$ implying that $\eta_2^3\in F^*$ and thus forcing $\eta_2^3=1$ contradicting the choices we made.

*If $\eta_1^3\neq\eta_2^3$ but exactly one of these cubes happens to equal to $1$, then the argument from the last line of the preceding bullet shows that one of $y_1,y_2$ vanishes, and we get no new solutions to $(2)$.

*On the other hand in the case $\eta_1^3=1=\eta_2^3$ (three choices for both $\eta_1$ and $\eta_2$, so a total of nine combinations), the equation $(2)$ becomes
$$y_1+y_2=1.$$ As we prescribed $y_1$ and $y_2$ to both be non-zero, this equation holds for $q-2$ pairs $(y_1,y_2)$ in the field $F$.


The conclusion is thus that $(2)$ has
$$
(q-2)(q-5)+9(q-2)=q^2+2q-8
$$
solutions. Adding the nine points with a single vanishing homogeneous coordinate we see that
$$
\#C(E)=q^2+2q-8+9=(q+1)^2
$$
verifying the answer gotten by more high-powered methods.

Unfortunately I couldn't make this argument work in the case $2\mid m$ on either of first two attempts. This time cubing is no longer bijective in $F^*$. It is bijective in $S$, but the analogue of $(2)$:
$$
y_1^3\eta_1+y_2^3\eta_2=1\tag{2'}
$$
is not as easy to analyze:


*

*Using linear dependence of $\mathcal{B}=\{\eta_1,\eta_2\}$ when applicable is not straight forward for it is difficult to analyze the occasions where the coordinates of $1$ with respect to the basis $\mathcal{B}$ should be cubes in $F$.

*There is the further observation that in the special case $\eta_1=1=\eta_2$ the equation $(2')$ becomes $(1)$, but with the variables ranging over the smaller field $F$. This is not a problem per se, as we can use induction on $m$ to handle those points. However, it also implies that the number of solutions of the type $\eta_1\neq1, \eta_2\neq1$ will depend on the residue class of $m$ modulo four. I didn't spot a way forward reflecting this.

*An alternative (initially my main line of attack) would be to use $(2')$ to solve one of the $\eta_i$s (the equation is linear!). All the elements of $S$ satisfy the equation $\eta^q=\eta^{-1}$. Together with $y_i^q=y_i$ an application of the Frobenius shows that the solution $\eta_2$ of $(2')$ is an element of $S$ if and only if the solutions (=the choices of $\eta_1$) of the quadratic
$$
\eta^2+\frac{1+y_1^6+y_2^6}{y_1^3}+1=0
$$
are not elements of $F$ (when they will necessarily be in $S$). This leads to a trace condition, and it looks like the number of solutions can only be countes with the aid of a character sum, which is something I wanted to avoid.


Anyway, if I come up with another way of handling $(2')$, I will edit this. Apologies for failing to deliver.
