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?
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} $$
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.