I was fiddling around with this problem for 3x3 magic squares after seeing another Numberphile video and I got to a point where I'm not sure where the error in proving no such magic squares exists is, so I would appreciate someone pointing it out.
So, consider a generic magic square of the following form:
$\begin{pmatrix} a^2 & b^2 & c^2 \\ d^2 & e^2 & f^2 \\ g^2 & h^2 & i^2 \end{pmatrix}$
and a sum $s$. Using some simple arithmetic and the basic properties of magic squares (the rows, columns and diagonals sum up to $s$), we can see that any magic square of squares can be written as follows:
$\begin{pmatrix} 2 t^2 - \frac{b^2+d^2}{2} & b^2 & t^2 - \frac{b^2-d^2}{2} \\ d^2 & t^2 & 2t^2-d^2 \\ t^2 + \frac{b^2-d^2}{2} & 2 t^2-b^2 & \frac{b^2+d^2}{2} \end{pmatrix}$
where $t^2=\frac{s}{3}$. Now, we can see that the elements in the second diagonal form an arithmetic sequence, and since we know these should be square numbers, we know we can parameterize them using the following replacements (Arithmetic sequence of squares):
$t^2 =(m^2+n^2)^2$
$\frac{b^2-d^2}{2} = 4 m n (m^2-n^2)$ - the congruum
Now, replacing the values of $t$ and $d$ in the matrix and simplifying the terms, we get the elements on the main diagonal to be:
$a^2=(m^2+n^2)^2+\left((m^2+n^2)^2+4 m n (m^2-n^2) - b^2\right)$
$e^2=(m^2+n^2)^2$
$i^2=(m^2+n^2)^2-\left((m^2+n^2)^2+4 m n (m^2-n^2) - b^2\right)$
Again, these are all squares and again form an arithmetic sequence so we can parameterize them in the same way, but since the middle element is already in parameterized form, we know that the congruum of this sequence must be of the form $4 m n (m^2-n^2)$. And now we have two possibilities:
- $(m^2+n^2)^2+4 m n (m^2-n^2) - b^2=4 m n (m^2-n^2)$
which gives $b^2=(m^2+n^2)^2$, but doing the replacement makes all the entries of the magic square equal, or
- $(m^2+n^2)^2+4 m n (m^2-n^2) - b^2=-4 m n (m^2-n^2)$
which when used to express $b^2$ and replace it in the magic square gives the second row to have all the same elements $(m^2+n^2)^2$, all the corner elements to be equal to $(m^2-2mn-n^2)^2$, and $b^2$ and $h^2$ to be equal to $(m^2-n^2)^2-8 m n (m^2-n^2)$.
So, the conclusion is that the only magic squares of squares are the ones with at most three distinct values.