"Nearly" Fermat triples: case cubic

Any such matrix $M$ would give rise to an automorphism of the cubic surfaces $$a^3 + b^3 = c^3 \pm d^3 \quad \subset \mathbb{P}^3.$$ These are both just different ways of writing the Fermat cubic surface $$x_0^3 + x_1^3 + x_2^3 + x_3^3 = 0 \quad \subset \mathbb{P}^3.$$ It is well-known that the automorphism group of the Fermat cubic surface over $\mathbb{C}$ is generated by the "obvious automorphisms", namely by permuting coordinates and multiplying coordinates by a third root of unity (See Table 9.6 of "Dolgachev - Classical algebraic geometry" for this claim).

Hence, as you are interested with matrices with integer entries, we see that the only such $M$ are the 6 permutation matrices which permute $a,b,c$, with also possibly multiplying some variable by $-1$ to fix the signs.

For example \begin{bmatrix}0&1&0 \\ 1&0&0 \\ 0&0&1\end{bmatrix} is an example of such a matrix. The other matrices can be written down analogously.

This writeup addresses the orginal question, as stated, just marginally. However, I hope to answer in more detail Amdeberhan's "other thoughts when he did ask the question, as he wrote.

Actually, for every rational integer solution $(a,b,c)$ to $a^3+b^3-c^3=1$ there is a matrix $M$ of infinite order, depending on $(a,b,c)$, such that for all integer $n$ we have that $(a',b',c')=M^n\cdot (a,b,c)^T$ is another triple of integers satisfying the equation.

To see this, we can observe that there are parametrizations $x_i=f_i(p,q)$ of solutions to $x_1^3+x_2^3+x_3^3+x_4^3=0$, such that:

• $f_1$, $f_2$, $f_3$, $f_4$ are quadratic polynomials in $p,q$;
• $f_4(p,q)=p^2-Dq^2$ for some squarefree $D$;
• $f_1(1,0)=-a$, $f_2(1,0)=-b$, $f_3(1,0)=c$.

Then we solve the quadratic Pell equation $f_4(p,q)=1$. In other words, we find a primitive solution $(p_1,q_1)$ and we generate the sequence $(p_n,q_n)$ by $p_n+\sqrt D q_n = (p_1+\sqrt D q_1)^n$.

We plug in the parametrization and we find $(a_n,b_n,c_n)$. Now, the key observation is that the transition from $n$ to $n+1$ is a linear transformation of the triple, so we get our matrix that does the job.

Below, I give a numerical example. I use the parametrization (5) given in this MSE post, with $a=-10, b=12, c=-9, d=1$, after the linear change of coordinates $p=x+22y$, $q=2y$. Unfortunately a solution to the Pell equation with $D=85$ is quite big: $p_1=285769$, $q_1=30996$. But it does give a matrix of integers $$M=\begin{pmatrix}-331671644135 & -512714878704& -314190334080\\ 266762362608 & 412374813481 & 252702205056\\ 259637126112 & 401360260896 & 245952516097\end{pmatrix} $$ that predictably works: $M\cdot (-10,12,-9)= (-8149096378, 6554290188, 6379224759)$, and $(-8149096378)^3+6554290188^3=(-6379224759)^3+1$.

To help comparing with Stadnicki's analysis, I report that $M$ has characteristic polynomial $(x - 1) \cdot (x^2 - 326655685442 x + 1)$ and eigenvectors $v_{max}=(-\frac 1 8 (1+ \sqrt{85}),-\frac 1 8 (1- \sqrt{85}), 1)$, $v_1=(-16,-16,43)$, $v_{min}=(-\frac 1 8 (1- \sqrt{85},-\frac 1 8 (1+ \sqrt{85}),1)$.