Counting points on the Klein quartic
$\def\FF{\mathbb{F}}\def\PP{\mathbb{P}}\def\ZZ{\mathbb{Z}}$ Throughout this answer, write $\zeta$ for a primitive $7$-th root of unity. (What field $\zeta$ is in will vary from one part of the answer to another.)
A great reference for the material I'm discussing here is Elkies' notes on the number theoretic properties of the Klein quartic.
Short explanation for $(\ast)$ Set $X = \{ (u:v:w) : u+v+w=0 \} \subset \PP^2$. There is a map $\phi: K \to X$ given by $\phi(x:y:z) = (x^3 y: y^3 z : x z^3)$. This is a $7$ to $1$ covering, branched over $(1:-1:0)$, $(0:1:-1)$ and $(-1:0:1)$. Formally, this map seems to be undefined at $(1:0:0)$, $(0:1:0)$ and $(0:0:1)$, but one can check that $\phi$ extends to these points and maps them to the branch points. I'll write $X^{\circ}$ and $K^{\circ}$ for $X$ and $K$ with these special points removed. See also Elkies Section 2.1.
For $(u:1-u:1) \in X^{\circ}$, the preimage of $(u:1-u:1)$ is $$(x:y:1) = \left(\sqrt[7]{\frac{u^3}{1-u}}:\ \frac{u}{x^2}: 1 \right).$$
Since $\phi$ is defined over $\FF_q$ for any $q$, we have $\phi(K(\FF_q)) \subset X(\FF_q)$. Now assume $q \not \equiv 1 \bmod 7$. Then every element of $\FF_q$ has a unique $7$-th root in $\FF_q$, so the displayed equation above shows that $\phi$ is bijective from $K^0(\FF_q)$ to $X^0(\FF_q)$. The three missing points in both $K$ and $X$ are also defined over $\FF_q$, so $\phi$ is a bijection $X(\FF_q) \to K(\FF_q)$, and $K(\FF_q)$ clearly has $q+1$ points.
Moduli interpretation The map $\phi$ is the quotient of $K$ by the automorphism $(x:y:z) \mapsto (\zeta x: \zeta^2 y: \zeta^4 z)$. This order $7$ automorphism must correspond to the element $\left( \begin{smallmatrix} 1 & 1 \\ 0 & 1 \end{smallmatrix} \right)$ in $PSL_2(\FF_7) \cong \mathrm{Aut}(K)$. So $X$ is the quotient of $K$ by $\left( \begin{smallmatrix} 1 & \ast \\ 0 & 1 \end{smallmatrix} \right)$ and we see that $X = X_1(7)$. So the question requiring modular interpretation is why $X(7)(\FF_q) \to X_1(7)(\FF_q)$ is bijective for $q \not \equiv 1 \bmod 7$. (Actually, there is a slight error here! See below.)
I think I can answer this, but my answer requires you are comfortable with group schemes. Throughout this answer, I will ignore the cusps of the modular curves.
We write $Z/7$ for the cyclic group scheme of order $7$ and $\mu_7$ for the group scheme of $7$-th roots of unity. So $\#(\ZZ/7)(\FF_q) =7$ but $\# \mu_7(\FF_q) = GCD(q-1,7)$. Then (see Elkies, section 4.1), $X(7)(\FF_q)$ parametrizes elliptic curves $E$ defined over $\FF_q$ with a chosen isomorphism $E[7] \cong \ZZ/7 \times \mu_7$, compatible with the Weil pairing. Meanwhile, $X_1(7)(\FF_q)$ parametrizes $E$ defined over $\FF_q$ with a chosen embedding of $\ZZ/7$ into $E[7]$. See also this answer by Pete Clark.
We must show the following:
Let $q \not \equiv 1 \bmod 7$. Given an elliptic curve $E$ defined over $\FF_q$ and an embedding $\ZZ/7 \to E[7]$, there is exactly one way to extend this embedding to an isomorphism $\ZZ/7 \times \mu_7 \to E[7]$, compatible with the Weil pairing.
Proof sketch: I'll assume $q \neq 7$. Consider the action $F$ of Frobenius on $E[7]$. This is given by a $2 \times 2$ matrix $\Phi$ with determinant $q$. Since $\ZZ/7$ embeds in $E[7]$, one of the eigenvalues of $\Phi$ is $1$. Since $q \not \equiv 1 \bmod 7$, we get that $\Phi$ is diagonalizable with eigenvalues $1$ and $q$. The eigenspace for $q$ is isomorphic to $\mu_7$. So $E[7] \cong \ZZ/7 \times \mu_7$. Since $q \not \equiv 1 \bmod 7$, there are no nontrivial maps between $\ZZ/7$ and $\mu_7$, so the isomorphism is unique up to twisting by an automorphism on each factor. The isomorphism on the $\ZZ/7$ factor is given to us, and the Weil pairing condition fixes the isomorphism on the $\mu_7$ factor. $\square$
Error corrected From section 4.2 of Elkies, I learned that this is not quite right! The map $K \to X$ corresponds not to forgetting the $\mu_7$ factor of $E[7]$, but to forgetting the $\ZZ/7$ factor! So $X$ is actually the moduli space of curves $E$ with an injection $\mu_7 \to E[7]$. It happens that both $X$ and $X_1(7)$ are isomorphic to $\PP^1$, but on the level of moduli, we should be taking a curve $E$ with an injection $\mu_7 \to E[7]$ and lifting it to an isomorphism $\ZZ/7 \times \mu_7 \cong E[7]$. Fortunately, the same argument works for this.
Ah! So, it suffices to understand what the $b_q$ are doing. The link above notes that the elliptic curve $E/\mathbf{Q}$ in question has $\mathrm{End}(E) = \mathbf{Z}[\frac{1+\sqrt{-7}}{2}]$, that is, it is a CM curve. My above observation about quadratic progressions is then explained by this paper (page 1, no less!)
In this case, though, it seems possible to get an even more explicit description...
$\def\FF{\mathbb{F}}$I'm switching to another answer to discuss relations with $X_0(49)$ and complex multiplication. Again, Elkies notes are a good reference and $\zeta$ is a primitive $7$-th root of unity.
Let $\Delta \subset PSL_2(\FF_7)$ be the group of diagonal matrices. Let $Y = \Delta \backslash K$. Then $Y$ is the modular curve corresponding to the arithmetic subgroup $\left( \begin{smallmatrix} a & 7b \\ 7c & d \end{smallmatrix} \right)$. Conjugating by $\left( \begin{smallmatrix} 7 & 0 \\ 0 & 1 \end{smallmatrix} \right)$ turns this group into $\left( \begin{smallmatrix} a & b \\ 49c & d \end{smallmatrix} \right)$ which is $\Gamma_0(49)$, so $Y \cong X_0(49)$. On the level of moduli, when our base field contains a $7$-th root of unity, $K$ parametrizes elliptic curves with a basis for $E[7]$ (up to a condition on the determinant) and $Y$ parametrizes curves with a chosen splitting of $E[7]$ into two cyclic subgroups. See the penultimate paragraph of Elkies Section 4.2 for why $Y \cong X_0(49)$ from a moduli perspective.
We can also use this perspective to explain why $Y$ has CM by $(1+\sqrt{-7})/2$. The curve $Y$ has genus $1$; make it into an elliptic curve by choosing the cusp at $\infty$ as the origin. Let $y \in Y$ and let $z \in K$ be a preimage of $y$. Write $\pi$ for the map $K \to Y$. Define $$\alpha(y) = \pi \left( \begin{pmatrix} 1 & 1 \\ 0 & 1 \end{pmatrix} z \right)+ \pi \left( \begin{pmatrix} 1 & 2 \\ 0 & 1 \end{pmatrix} z \right)+ \pi \left( \begin{pmatrix} 1 & 4 \\ 0 & 1 \end{pmatrix} z \right)$$ and $$\beta(y) = \pi \left( \begin{pmatrix} 1 & 3 \\ 0 & 1 \end{pmatrix} z \right)+ \pi \left( \begin{pmatrix} 1 & 5 \\ 0 & 1 \end{pmatrix} z \right)+ \pi \left( \begin{pmatrix} 1 & 6 \\ 0 & 1 \end{pmatrix} z \right)$$ where the sums are in the group law on $Y$.
Note that replacing $z$ by a different preimage $\left( \begin{smallmatrix} a & 0 \\ 0 & a^{-1} \end{smallmatrix} \right) z$ would just turn $\alpha(y)$ into $$\pi \left( \begin{pmatrix} a & a^{-1} \\ 0 & a^{-1} \end{pmatrix} z \right)+ \pi \left( \begin{pmatrix} a & 2 a^{-1} \\ 0 & a^{-1} \end{pmatrix} z \right)+ \pi \left( \begin{pmatrix} a & 4 a^{-1} \\ 0 & a^{-1} \end{pmatrix} z \right)=$$ $$\pi \left( \begin{pmatrix} 1 & a^{-2} \\ 0 & 1 \end{pmatrix} z \right)+ \pi \left( \begin{pmatrix} 1 & 2 a^{-2} \\ 0 & 1 \end{pmatrix} z \right)+ \pi \left( \begin{pmatrix} 1 & 4 a^{-2} \\ 0 & 1 \end{pmatrix} z \right) = \alpha(y).$$ To get to the second line, we used that $\pi$ is the quotient by $\Delta$, so $\pi\left( \left( \begin{smallmatrix} a & 0 \\ 0 & a^{-1} \end{smallmatrix} \right) w \right) = \pi(w)$. So $\alpha$ and $\beta$ are well defined endomorphisms of $Y$.
We claim that $\alpha+\beta+1=0$ and $\alpha \beta = \beta \alpha = 2$. So $\alpha$ and $\beta$ generate a quadratic ring, with $\alpha$ and $\beta$ corresponding to $\frac{-1 \pm \sqrt{-7}}{2}$.
We verify that $\alpha+\beta+1=0$: Writing this out, we must show that $$\sum_{b=0}^6 \pi\left( \begin{pmatrix} 1 & b \\ 0 & 1 \end{pmatrix} z \right) =0.$$
Write $\psi$ for the map $Y \to X_0(7)$ which forgets one of the two summands into which $Y$ decomposes $E[7]$. Then the above sum is $\sum_{\psi(y')= \psi(y) } y'$. But $X_0(7)$ has genus $0$, so the class of $\sum_{\psi(y') = w} y'$ is independent of the choice of point $w \in X_0(7)$, and taking $y$ to be the cusp shows that this constant value is $0$.
We sketch the verification of $\alpha \beta = \beta \alpha = 2$. Multiplying out both sides, we must show $$3 y + \sum_{b=1}^6 \pi\left( \begin{pmatrix} 1 & b \\ 0 & 1 \end{pmatrix} z \right) = 2y.$$
Cancelling $2y$ from both sides, this turns into our previous verification.
Of course, all of this is just using that $\mathbb{Z}[\zeta]$ acts on the Jacobian of $K$ and that $\zeta+\zeta^2+\zeta^4 = \frac{-1+\sqrt{-7}}{2}$, but I think it is fun to see it in terms of matrices.