Start from $1105 = 5\cdot 13\cdot 17$ whose prime factors all has the form
$4k+1$.
Decompose the prime factors in $\mathbb{Z}$ over gaussian integers $\mathbb{Z}[i]$ as
$$5 = (2+i)(2-i),\quad 13 = (3+2i)(3-2i)\quad\text{ and }\quad17 = (4+i)(4-i)$$
Recombine the factors of $1105^2$ over $\mathbb{Z}[i]$ in different order
and then turn them to sum of two squares. We get:
$$\begin{array}{rc:rr}
1105^2
= & 943^2 + 576^2 & ((2+i)(3+2i)(4+i))^2 = & -943 + 576i\\
= & 817^2 + 744^2 & ((2-i)(3+2i)(4+i))^2 = & 817 + 744i\\
= & 1073^2 + 264^2 & ((2+i)(3-2i)(4+i))^2 = & 1073 + 264i\\
= & 1104^2 + 47^2 & ((2-i)(3-2i)(4+i))^2 = & -47 - 1104i\\
\end{array}
$$
A counter-example for the speculation that an integer can appear in at most two primitive Pythagorean triples.
Here is an approach to generate solutions where all the terms in each individual sum of square are relatively prime. It's not a direct formula, but an algorithm that uses Pythagorean triples as inputs. We can, of course, generate the Pythagorean triples from the well-known formula.
Begin with two primitive Pythagorean triples. Here I will use
$3^2+4^2=5^2$
$5^2+12^2=13^2$
We multiply these together, using the Brahmagupta-Fibonacci identity on the left side:
$(3×5+4×12)^2+(3×12-4×5)%2=5^2×13^2$
$63^2+16^2=65^2$
Now we get a second solution by going back to our original triples and rewriting obe of them in reverse order:
$3^2+4^2=5^2$
$\color{blue}{12^2+5^2=13^2}$
Now with the Brahmagupta-Fibonacci identity our product relation is
$(3×12+4×5)^2+(3×5-4×12)^2=5^2×13^2$
$56^2+33^2=65^2$
Putting these results together we now have
$63^2+16^2=56^2+33^2=65^2.$
To get twice a square we then replace $a^2+b^2$ with its doubled value $(a+b)^2+(a-b)^2$, and so with the $3-4-5$ and $5-12-13$ triples as inputs we end with
$\color{blue}{79^2+47^2=89^2+23^2=2×65^2}.$
$6241+2209=7921+529=8450☆$
These are the relatively prime combinations with $x=65$, but there are more solutions where the component squares have a common factor. Go back to our Pythagorean triples and replace one with its $c^2+0^2=c^2$ form, keeping its hypotenuse the same as in the original triple. If, in the pair we are considering, we replace $3^2+4^2$ with $5^2+0^2$, we get:
$\color{blue}{5^2+0^2=5^2}$
$5^2+12^2=13^2$
Hence
$(5×5)^2+(5×12)^2=5^2×13^2$
$25^2+60^2=65^2$
$\color{blue}{85^2+35^2=2×65^2}$
Similarly, keeping $3^2+4^2$ but replacing $5^2+12^2$ with $13^2+0^2$ will yield
$\color{blue}{91^2+13^2=2×65^2}$
Thus ultimately, we get all the following for $x=65$:
$\color{blue}{79^2+47^2=85^2+35^2=89^2+23^2=91^2+13^2=2×65^2}.$
The reader can verify that all members if this equation equal $8450$.
It might be noted that while this approach is attractive for giving solutions, it does not mean such solutions can be readily pieced together into a $3×3$ magic square of squares. The numbers generated in the above solution are much smaller than those that are now known to be required if such a square were to exist.
Best Answer
Hint: For $a^2+b^2=2c^2$, observe that $a, b$ have the same parity. Therefore there exist integers $u, v$ such that $a = u+v$ and $b = u-v$. Expand...