Set of linear equations in quadratic integers: is a solution possible? ("Magic square of squares")

Question

This problem is for me more a "recreational" one, but might even have no solution or better: is part of an open problem. The set of equations comes from the problem of "magic 3x3 square of squares" where the $9$ unknowns are written as $$M = \small \begin{bmatrix} a^2 & b^2 & c^2 \\ d^2 & e^2 & f^2 \\ g^2 & h^2 & i^2 \\ \end{bmatrix} $$ where row, column and diagonal sums shall equal the same value (say $3 e^2$). By the rules of magic squares all values $a²,b²,...,i²>0$ shall be pairwise distinct.

---

I managed to get the circular sets of linear equations in quadratics $$\tag 1 \small \begin{bmatrix} -1 & 1 & 1 \\ 4 & -1 & -2 \\ 2 & 0 & -1 \\ \end{bmatrix} \cdot \small \begin{bmatrix} e^2 \\ h^2 \\ i^2 \\ \end{bmatrix} = \small \begin{bmatrix} c^2 \\ f^2 \\ a^2 \\ \end{bmatrix} $$

$$\tag 2 \small \begin{bmatrix} 2 & 3 & -4 \\ 1 & 2 & -2 \\ 2 & 1 & -2 \\ \end{bmatrix} \cdot \small \begin{bmatrix} c^2 \\ f^2 \\ a^2 \\ \end{bmatrix} = \small \begin{bmatrix} b^2 \\ d^2 \\ g^2 \\ \end{bmatrix} $$

$$\tag 3 \small \begin{bmatrix} -1/2 & 1 & 1/2 \\ -2 & 2 & 1 \\ 1/2 & 0 & 1/2 \\ \end{bmatrix} \cdot \small \begin{bmatrix} b^2 \\ d^2 \\ g^2 \\ \end{bmatrix} = \small \begin{bmatrix} e^2 \\ h^2 \\ i^2 \\ \end{bmatrix} $$ $ \qquad \qquad $ (Matrix (2) has been corrected against first version due to the comment)

I think that this is an interesting diophantine problem in itself: to handle connected equations in $4$ squares. By L. Euler we know already parametrizations for possible solutions for $4$-squares-equations, but I'm not familiar enough with that to try to apply this here an to get -possibly- contradictions or an infinite ascent (or a solution). Likewise we have at least two connected Pell-equations (where some zero matrix-coefficients occur). I know already that all squares must be congruent $1$ to modulus $24$ (the unsquared values $a,b,c...,i= \pm1 \pmod 6$) and we have at most three free parameters (I start with $e^2,h^2,i^2$ as independents).

For instance, we can extract a set of connected equalities, each only involving three variables, which are very similar to Pell-equations.

$$\small \begin{array} {} 2i^2 &= b^2+d^2 \\ 2g^2 &= b^2+f^2 \\ 2e^2 &= a^2+i^2 &= b^2+h^2 &= c^2+g^2 &= d^2+f^2 \\ 2c^2 &= d^2+h^2 \\ 2a^2 &= f^2+h^2 \\ \end{array}$$

Q1: How could I possibly improve that ansatz?

or even:

Q2: Is a nonzero solution (besides a trivial solution with all $a=b=...=i=1$ which is excluded by definition of the problem) possible?

---

Looking at what I've observed so far

1) Hmm. I do not really look farther with this, but when I write the system of equations in a joint structure of a matrix-diagonalization, I get first this: $$ M= \small \begin{bmatrix} 0 & 0 & 0 & 0 & 0 & 0 & -1 & 1 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 4 & -1 & -2 \\ 0 & 0 & 0 & 0 & 0 & 0 & 2 & 0 & -1 \\ 2 & 3 & -4 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 2 & -2 & 0 & 0 & 0 & 0 & 0 & 0 \\ 2 & 1 & -2 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & -1/2 & 1 & 1/2 & 0 & 0 & 0 \\ 0 & 0 & 0 & -2 & 2 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1/2 & 0 & 1/2 & 0 & 0 & 0 \end{bmatrix}\cdot \small \begin{bmatrix} c^2\\f^2\\a^2\\b^2\\d^2\\g^2\\e^2\\h^2\\i^2 \end{bmatrix} = \small \begin{bmatrix} c^2\\f^2\\a^2\\b^2\\d^2\\g^2\\e^2\\h^2\\i^2 \end{bmatrix} $$ By this $M$ the vector of the squares of the variables $c^2,f^2,...i^2$ is an eigenvector to the eigenvalue $1$; the eigenvalues are the three cube-roots of complex unit and each occurs threefold. We have only three real eigenvectors and they just encode the composition of $c^2,f^2,...i^2$ by the ("independents") $e^2,h^2,i^2$.
...

2) The third power $M^3$ is just the identity matrix, so no new insights come from here.

3) ... here my "expertise" in linear algebra is at its limit at the moment and possibly a helpful hint could emerge from that observations/speculations so far.

Answer 1

Generating the hyperplanes for the eigenvectors and searching through them would be indistinguishable from simply trying solutions to the original equations.

This problem can be formulated in many branches of mathematics. For example, those sequences of squares "which are very similar to Pell-equations" can be written as $$b^2-i^2=i^2-d^2\quad=\quad g^2-e^2=e^2-c^2\quad = \quad f^2-a^2=a^2-h^2$$ $$b^2-g^2=g^2-f^2\quad=\quad i^2-e^2=e^2-a^2\quad = \quad d^2-c^2=c^2-h^2$$ $$b^2-e^2=e^2-h^2\quad\text{and}\quad d^2-e^2=e^2-f^2$$ That is to say, there are 8 arithmetic sequences of 3 square numbers. In fact, if we rearrange the square to be $$M'=\begin{bmatrix}b^2&i^2&d^2\\g^2&e^2&c^2\\f^2&a^2&h^2\end{bmatrix}$$ each row, column, and the two main diagonals form an arithmetic sequence. Additionally, the column sequences share the same common difference and the row sequences share the same common difference. If we let $a,b,...,h,i$ be complex numbers, then the matrix $M'$ forms a $3\times 3$ grid in the complex plane and can be written as $$\begin{bmatrix}b^2&i^2&d^2\\g^2&e^2&c^2\\f^2&a^2&h^2\end{bmatrix}= A\begin{bmatrix}1&1&1\\1&1&1\\1&1&1\end{bmatrix} +B\begin{bmatrix}1&1&1\\0&0&0\\-1&-1&-1\end{bmatrix} +C\begin{bmatrix}1&0&-1\\1&0&-1\\1&0&-1\end{bmatrix} $$ The matrices you found can be understood as symmetries of this "grid".

The common differences are called $\textit{congrua}$ and have been fully parametrized for a long time. $$\text{if }r^2-s^2=s^2-t^2\quad \text{then there exists integers } m,n,k \text{ such that}$$ $$r=k(m^2+2mn-n^2)$$ $$s=k(m^2+n^2)$$ $$t=k(m^2-2mn-n^2)$$ $k$ is pretty much a scaling factor. $m$ and $n$ are "doing the work" so to speak. Interestingly, if we pick the right complex number: $\omega =\sqrt{k}(m+ni)$ then we get the parametrization: $$r=\text{Re}[\omega^2(1-i)]$$ $$s=|\omega|^2$$ $$t=\text{Re}[\omega^2(1+i)]$$ (We should have a function for going between $\omega$ and $r,s,$ and $t$ so let $g(\omega)=(r,s,t)$.) The common difference can even be expressed with $\omega$: $$r^2-s^2=s^2-t^2=\text{Im}[\omega^4]$$ The last interesting property to be mentioned is that $rt=\text{Re}[\omega^4]$.

Returning to the magic square of squares, there's complex numbers like $\omega$ for each of those 8 sequences. Picking 4, let $$g(\alpha)=(a,e,i),\quad g(\beta)=(b,e,h), \quad g(\gamma)=(c,e,g), \quad\text{and}\quad g(\delta)=(d,e,f)$$ These 4 were picked because they each have the same middle element $e$, which means $$|\alpha|=|\beta|=|\gamma|=|\delta|=\sqrt{e}$$ And if we care enough to work out the implications of the other linear equations, it follows that $$\text{Im}[\alpha^4]+\text{Im}[\beta^4]+\text{Im}[\gamma^4]=\text{Im}[\alpha^4]+\text{Im}[\delta^4]-\text{Im}[\gamma^4]=0$$ or equivalently, that $$\alpha^4+\beta^4+\gamma^4,\quad \alpha^4+\delta^4-\gamma^4\quad\in \mathbb{R}$$ If you want a reason that no one has found a solution by brute force, notice that the 4th powers of complex numbers $\omega^4$ get spread out very quickly in the complex plane as you move away from the origin.

The last point to be made (and this is the farthest I've gotten going down this route of complex number parametrization) is that if we restrict $\alpha, \beta,$ and $\gamma$ to have integer parts, then the fact that $|\alpha|=|\beta|=|\gamma|$ tells us that they have the same prime factorization up to conjugation and rotation in the Gaussian integers. So without loss of generality, let $w,x,y,z\in \mathbb{Z}[i]$ such that $$\alpha=w\bar{x}yz \quad \beta=wx\bar{y}z \quad \gamma=wxy\bar{z}$$ (We are categorizing how the primes are conjugated in $\alpha, \beta,$ and $\gamma$. $w$ is the product of all the primes with the same conjugation in all 3. $x$ is the product of all the primes with the same conjugation in $\beta$ and $\gamma$ but not in $\alpha$. Etc.)

Rotation need not be worried about since raising a Gaussian integer to the 4th power brings the 4 roots of unity to $1$. (I.e. $1^4=i^4=(-1)^4=(-i)^4=1$)

The earlier equation - with the correct manipulations - can now be transformed as follows $$\alpha^4+\beta^4+\gamma^4\in \mathbb{R}$$ $$\Updownarrow$$ $$\text{Im}[(w\overline{xyz})^4]=4\text{Im}[w^4]\text{Re}[x^4]\text{Re}[y^4]\text{Re}[z^4]+4\text{Re}[w^4]\text{Im}[x^4]\text{Im}[y^4]\text{Im}[z^4]$$ This has some interesting consequences, such as if we restrict to $w=1$, then the equation becomes $$\text{Im}[(xyz)^4]=-4\text{Im}[x^4]\text{Im}[y^4]\text{Im}[z^4]$$ The RHS is divisible by $4\cdot 24^4$ which gives us information about $xyz$ on the LHS. (I used this to do a computer search on on $xyz=m+ni$ corresponding to $a,b,c,d,e,f,g,h,i<10^{14}$)

Some other interesting areas for tackling this problem are elliptic curves and modular arithmetic. Unfortunately, there isn't much literature published explicitly on the magic square of squares. I've had the chance to talk to two professional number theorists about the problem and they didn't have anything to say - at the time - beyond what has been written here more or less.