<s>[INST] Pythagorean triples where the sum of the two cubes is also a square Are there any Primitive Pythagorean triple solutions $(a,b,c)$ where the sum of the two cubes is also a square? In other words are there coprime $a,b>0 \in \mathbb{N} \;, (a,b)=1$  where $a^2+b^2=c^2$ and $a^3+b^3=d^2$ for some $c,d \in \mathbb{N}$
 [/INST] Disclaimer: These are some unfinished thoughts I will leave here to work on later, or for others to continue.

Given that $a$ and $b$ are coprime, it follows that $\gcd(a+b,a^2-ab+b^2)$ divides $3$ because
$$\gcd(a+b,a^2-ab+b^2)=\gcd(a+b,3b^2)=\gcd(a+b,3).$$
Suppose towards a contradiction that the gcd equals $3$: Then the factorization
$$d^2=a^3+b^3=(a+b)(a^2-ab+b^2),$$
shows that there exist $e,f\in\Bbb{Z}$ such that
$$a+b=3e^2\qquad\text{ and }\qquad a^2-ab+b^2=3f^2,$$
from which it quickly follows that
$$9e^4=(a+b)^2=a^2+2ab+b^2=3c^2-6f^2,$$
and reducing mod $8$ yields a contradiction, so the gcd is $1$. Hence there exist $e,f\in\Bbb{Z}$ such that
$$a+b=e^2\qquad\text{ and }\qquad a^2-ab+b^2=f^2,$$
and in the same way as before we find that
$$e^4=(a+b)^2=a^2+2ab+b^2=3c^2-2f^2.$$
Luckily $\Bbb{Z}[\sqrt{6}]$ is a UFD, and we have
$$N((3c-2f)+(c-f)\sqrt{6}):=
\left((3c-2f)+(c-f)\sqrt{6}\right)\left((3c-2f)-(c-f)\sqrt{6}\right)
=3c^2-2f^2=e^4.$$
The gcd of two conjugate factors divides $2(3c-2f)$ and $2(c-f)$, and because $c$ and $f$ are coprime it follows that the gcd divides $2$. Because their product $e^4=(a+b)^2$ is odd, the two conjugate factors are in fact coprime. This means there exists some $x+y\sqrt{6}\in\Bbb{Z}[\sqrt{6}]$ such that
$$(3c-2f)+(c-f)\sqrt{6}=(x+y\sqrt{6})^4.\tag{1}$$
This immediately tells us that
$$a+b=e^2=\sqrt{N((3c-2f)+(c-f)\sqrt{6})}=(x+y\sqrt{6})^2(x-y\sqrt{6})^2=(x^2-6y^2)^2.\tag{2}$$
Furthermore, expanding equation $(1)$ yields the two equations
$$3c-2f=x^4+36x^2y^2+36y^4\qquad\text{ and }\qquad c-f=4x^3y+24xy^3.$$
Because $c-f>0$, without loss of generality $x,y>0$. The above tells us that
\begin{eqnarray*}
c&=&x^4-\ 8x^3y+36x^2y^2-48xy^3+36y^4,\\
f&=&x^4- 12x^3y+36x^2y^2-72xy^3+36y^4,
\end{eqnarray*}
and hence that
$$ab=c^2-f^2=(c-f)(c+f)=8xy(x^2+6y^2)(x^4-10x^3y+36x^2y^2-60xy^3+36y^4).\tag{3}$$
This means $a$ and $b$ are the roots of the quadratic polynomial
$$Z^2-(x^2-6y^2)^2Z+8xy(x^2+6y^2)(x^4-10x^3y+36x^2y^2-60xy^3+36y^4).$$
This polynomial has integer roots if and only if its discriminant $\Delta$ is a square, where
$$\Delta=(x^2-6y^2)^4-32xy(x^2+6y^2)(x^4-10x^3y+36x^2y^2-60xy^3+36y^4),$$
which leaves me with the question of when this homogeneous octic polynomial takes on square values.
</s>