Do I catch all solutions for the generalized Pell-equation $a^2+b^2 = 2 c^2$ by this matrix-method?

As one question in a somewhat bigger analysis I want to characterize the set of solutions of the generalized Pell-equation in the title: $$a^2+b^2 = 2 c^2 \tag 1$$ I'm not much fluent with the Pell-solving technology, and from some initial brute-force solutions I constructed the following scheme, which seems to me as if it could produce all solutions by changing a vector of parameters.
My question is, whether the set of solutions that I find is exhaustive or not - tests in small numbers seem to suggest that.

I introduce two matrices; $$T= \begin{bmatrix} 1&0&0\\0&3&2\\0&4&3 \end{bmatrix} \tag 2$$ $$C= \begin{bmatrix} 0&1&0\\1&0&0\\0&0&1 \end{bmatrix} \tag 3$$ The starting solution is the vector $A_0=[a,b,c]=[1,1,1]$.

With this I create further solutions by $$A_k=A_0 \cdot T^k \qquad k \in \mathbb Z \tag {4.1}$$ which gives a list with all solution having $a=1$. But from each of those solutions one can derive a new initial $A_0$ from where then again a complete list of solutions can be drawn: $$A_{k_1,k_2,k_3}=A_0 \cdot T^{k_1} \cdot C \cdot T^{k_2}\cdot C \cdot X^{k_3} \tag {4.2}$$ and so on, for arbitrarily many indexes/exponents $k_j$.
Using Pari/GP I implemented a function, which takes the vector of arbitrary length containing the indexes/exponents as argument in the following form $$ A_x = \operatorname{Pell2Sol}([k_1,k_2,...,k_n]) \tag 5$$ I don't know, whether the resulting structure can be called a tree or whether it has another name.

Q: As I said above, it seems to me that this indexing of solutions is exhaustive (assuming the number of indexes $n$ is taken from $1$ to $\infty$). Does someone know whether this is true?

_{Context: if this is a valid method, I'll try to extend/adapt it to some similar equations which occur in a combined set of equations in quadratics.}

\\ Update: initial vector A0 is configurable to have a gcd-common factor
\\ vector E takes the indexes
\\ env=0 says: give only that solution, 
\\ env=1 says: printout the set of neigboured solutions towards index +- env.
\\ A0 can be configured by a common factor in the function call
{Pell2Sol(E,env=0,gcd=1)=my(A0=gcd*[1,1,1], T=[1,0,0;0,3,2;0,4,3],C=[0,1,0;1,0,0;0,0,1]);
     if(#E==0,return(A0));
     A0*=prod(k=1,#E,C*T^E[k]); 
     if(!env,return(A0));
     return(Mat(vectorv(2*env+1,r,A0 * T^(r-1-env))));}

Giving Pell2Sol([1,0],3) shows the following sequence of solutions:

  a    b    c
 --------------
  7  -601  425
  7  -103   73
  7   -17   13
  7     1    5
  7    23   17
  7   137   97
  7   799  565

Here the negative signs should not be significant, since we use the squares anyway. But the signs are significant when computing the lists of solutions using the powers of the matrix T and the matrix C. I didn't supress the signs here to show that internal effect.

and Pell2Sol([1,5,1,0],3) gives

     a         b        c
  76793  -5375863  3801697
  76793   -920801   653365
  76793   -148943   118493
  76793     27143    57593
  76793    311801   227065
  76793   1843663  1304797
  76793  10750177  7601717

Perhaps a better picture than the two above is the following, (where I also draw the motivating connection to the problem of finding magic squares-of-squares) (Screenshot):

I'll assume that you are only considering positive integer solutions $(a,b,c)$, as your method only generates positive integer solutions. You already note that your method preserves $\gcd(a,b,c)$, and so the seed value $(1,1,1)$ can only generated primitive triples. The question then remains whether your method generates all primitive triples.

First let's characterize all primitive positive integral solutions to your equation. If $a$, $b$ and $c$ are positive integers with $\gcd(a,b,c)=1$ such that $$a^2+b^2=2c^2,$$ then as noted in the comments, we must have $a\equiv b\pmod{2}$ and so $x=\frac{a-b}{2}$ and $y=\frac{a+b}{2}$ are integers. They are easily verified be coprime, and if $(a,b,c)\neq(1,1,1)$ then we may interchange $a$ and $b$ if necessary so that $x$ and $y$ are both positive. It is easily verified that $$x^2+y^2=c^2,$$ meaning that $(x,y,c)$ is a primitive Pythagorean triple. Then it is well known that there exist coprime positive integers $m$ and $n$, not both odd and with $m>n$, such that $$x=m^2-n^2,\qquad y=2mn,\qquad c=m^2+n^2,$$ possibly with the roles of $x$ and $y$ reversed. Conversely, any such choice of $m$ and $n$ yields a primitive Pythagorean triple, and by extension any triple $$a=m^2+2mn-n^2,\qquad b=n^2+2mn-m^2,\qquad c=m^2+n^2,$$ is a solution to your original equation, and these are all primitive solutions.

Next let's consider the solutions generated from the seed $(1,1,1)$ by your method. Because $C^2=I$, the set of all generated solutions has the structure of a binary tree, by generating the two solutions $AT$ and $ATC$ from any previous solution $A$.

From every level of the binary tree to the next, we see that $c\mapsto\ 2b+3c$. This means that for every solution $(a,b,c)$ at the $k$-th level we have $c\geq3^k$, where level $k=0$ corresponds to the seed $A=(1,1,1)$. Because these are levels of a binary tree, this shows that your method generates at most $2^{k+1}$ solutions with $c<3^k$.

To conclude that your method does not generate all primitive solutions, it now suffices to show that the number of pairs of coprime positive integers $(m,n)$ with $m>n$ that satisfy $$m^2+n^2<3^k,$$ is greater than $2^{k+1}$. For this fact I refer to this paper, that shows as an example application of a more general theorem, that the number of such solutions is asymptotically proportional to $\frac{3^k}{2\pi}$.

Partial answer conjecture: No, this method does not catch up all solutions.

• By samples of observations: it seems, it only catches only all solutions with no common divisor in $[a,b,c]$. If we insert in the function's initialization part instead of $A_0=[1,1,1]$ a multiple like $A_0=17 \cdot [1,1,1]$ , then we get the analogue tree, but with every number triple having that same common divisor.

Remark: the existence of these solutions with common gcd are sometimes nontrivial - in the problem of "magic-squares-of-squares" ("msqosq") we have triples made from some entries with a common divisor and can only be found when $A_0=gcd \cdot [1,1,1]$ is given as starting point. For instance in the only known msqosq with $7$ quadratic entries ($2$ are not squares) there is a triple $[h,b,e]=[-289,527,425]$ which can be found by the proposed function, by changing initialization to $A_0=17 \cdot[1,1,1]$ and then perform Pell2Sol([-2,-1,1]).

I have no idea so far, whether this gcd-aspect is the only part that must be introduced into the question's indexing method.

Update
Triggered by the doubts as user Servaes presented in his/her answer, I tried to find hints, as when the Pythagorean-triple generation rules would construct one triple, which is not/cannot be in the Pell-solution- scheme as described in my question.

The first list is using the first generating-matrix ("WP1") as shown in wikipedia, starting at $A_{py,0}=[1 ,0, 1] $ $$ \begin{array} {l} \text{Pyth WP1 } \to \quad \text{Pell} \\ \small \begin{array} {rrr|rrr} x & y & c & a & b & c \\ \hline \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 1 & 0 & 1 & 1 & 1 & 1 \\ 3 & 4 & 5 & 7 & -1 & 5 \\ 5 & 12 & 13 & 17 & -7 & 13 \\ 7 & 24 & 25 & 31 & -17 & 25 \\ 9 & 40 & 41 & 49 & -31 & 41 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ \end{array} \end{array} $$

The second list is using the second generating-matrix ("WP2") as shown in wikipedia, again starting at $A_{py,0}=[1 ,0, 1] $

$$ \begin{array} {l} \text{Pyth WP2 } \to \quad \text{Pell} \\ \small \begin{array} {rrr|rrr} x & y & c & a & b & c \\ \hline \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 1 & 0 & 1 & 1 & 1 & 1 \\ 3 & 4 & 5 & 7 & -1 & 5 \\ 21 & 20 & 29 & 41 & 1 & 29 \\ 119 & 120 & 169 & 239 & -1 & 169 \\ 697 & 696 & 985 & 1393 & 1 & 985 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ \end{array} \end{array} $$

The third list is using the third generating-matrix ("WP3") as shown in wikipedia, now starting at $A_{py,0}=[3 ,4, 5] $. (A starting triple $[1,0,1]$ runs into a list of zeros)

$$ \begin{array} {l} \text{Pyth WP3 } \to \quad \text{Pell} \\ \small \begin{array} {rrr|rrr} x & y & c & a & b & c \\ \hline \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 3 & 4 & 5 & 7 & -1 & 5 \\ 15 & 8 & 17 & 23 & 7 & 17 \\ 35 & 12 & 37 & 47 & 23 & 37 \\ 63 & 16 & 65 & 79 & 47 & 65 \\ 99 & 20 & 101 & 119 & 79 & 101 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ \end{array} \end{array} $$

All three Pell-lists look very innocent, and trying some examples from this showed, that all generated Pell-triples could be generated by appropriate keys for the function $Pell2Sol()$