Determine all pairs of positive integers $(a,b)$ such that $a^2+b^2+ab$ is a perfect square. Consideration via mod $4$ shows that $a,b=0$ mod $4$ or one of between $a,b$ is $=1$ mod $4$ while the other is $=0 $ mod $4$.
Considering $(a+b)^2,a^2+an+b^2,(a+b-1)^2$ as we can deduce that $a,b>2$.
 A: Characterization of all solutions (non-trivial) to the Diophantine equation $$a^2+ab+b^2=x^2\quad a,b,x\in\mathbb{Z}$$ using Algebraic Number Theory

$\mathrm{WLOG}$ assume that $\gcd(a,b)=1$.We know that $\mathbb{Z}[\omega]$ is a $\mathrm{UFD}$ where $\omega=\frac{1+\sqrt{-3}}{2}$. Now we can factorize $$a^2+ab+b^2=(a-b\omega)(a-b\omega^2)=(a-b\omega)(a+b+b\omega)$$ in $\mathbb{Z}[\omega]$. Let $a^2+ab+b^2=x^2$ for some $x\in\mathbb{N}$. If $\pi\mid(a-b\omega)$ and $\pi\mid(a+b+b\omega)$ then we can show that $\pi\mid1$ and hence $\gcd(a-b\omega,a+b+b\omega)=1$. Then by unique factorization property we conclude $\exists$ $\alpha-\beta\omega,\gamma-\delta\omega\in\mathbb{Z}[\omega]$ such that $a-b\omega=u(\alpha-\beta\omega)^2=u(\alpha^2-2\alpha\beta\omega-\beta^2(1+\omega))=u((\alpha^2-\beta^2)-(2\alpha\beta+\beta^2)\omega)$ and $a+b+b\omega=v(\gamma-\delta\omega)^2$ where $u,v$ are units and hence $u,v\in\{\pm1,\pm\omega,\pm\omega^2\}$. We consider the case $u=1$ only, the other cases gives similar solutions. Comparing we get $a=\alpha^2-\beta^2$ and $b=\beta^2+2\alpha\beta$. We see that in fact $\gamma-\delta\omega=\overline{\alpha-\beta\omega}=\alpha-\beta\omega^2$ and hence $(a^2+ab+b^2)=(\alpha^2+\alpha\beta+\beta^2)^2$. So the set $$\{(\alpha^2-\beta^2,\beta^2+2\alpha\beta,\alpha^2+\alpha\beta+\beta^2):\alpha,\beta\in\mathbb{Z},\gcd(\alpha,\beta)=1\}$$ describes all solutions $(a,b,x)$ to the Diophantine equation $$a^2+ab+b^2=x^2$$ Also we observe that $(a,b,x)$ is a solution if and only if $(a,b,-x),(-a,-b,-x),(-a,-b,x),(b,a,x),(b,a,-x),(-b,-a,-x),(-b,-a,x)$ are solutions. These solutions will be obtained by considering other units except $u=1$. These observations are trivial.

A: Primitive solutions to $a^2+ab+b^2=c^2.$
Let $x=a/c, y=b/c.$
$x^2+xy+y^2=1\tag{1}$
Substitute  $x = t , y=1+kt$  to equation $(1)$, then we get  $$t = -\frac{(1+2k)}{(1+k+k^2)}$$
Thus, we get a parametric solution.  
$(x,y,z) = (-\frac{(1+2k)}{(1+k+k^2)}, -\frac{(k-1)(k+1)}{(1+k+k^2)})$ 
Hence $(a,b,c)=(1+2k, k^2-1, 1+k+k^2)$
$k$ is arbitrary. 
Further, take $k=\frac{\alpha}{\beta}$ then we get
$(a,b,c)=(\beta^2+2\alpha\beta,\alpha^2-\beta^2
,{\beta}^{2}+\alpha\,\beta+{\alpha}^{2})$
A: Just another way to approach the problem.
Since $$a^2+b^2+ab=(a+b)^2-ab $$ is a squared integer, there exists an integer $n$ such that
$$ (a+b)^2-(a+b-n)^2=ab$$
Solving it for $a$, we have
$$a = \frac{n (2 b - n)}{b - 2 n} $$
with the immediate condition  $b \neq 2 n$. However, we also have to set two other conditions. Because $a$ is positive, we must have $$ (2b-n)/(b-2n)>0$$ which for positive  $n$ leads to
$$b<\frac{n}{2} \,\,\,\text{or} \,\,\, b>2n$$
In addition, because $a$ is integer, we must have
$$(b-2n) |n (2 b - n)$$
Since $n(2b-n)$ can be written as $$n(2b-4n+3n)=2n(b-2n) +3n^2$$ it is divisible by $(b-2n)$ if and only if $(b-2n)$ divides $3n^2$. This occurs only when $b$ is of the form $$b=2n+k$$
where $k$ is a positive integer that divides $3n^2$.
Collecting all results, to find the pairs $a,b$ asked in the OP, we can:


*

*set an arbitrary positive integer $n$;

*take all positive integers $k$ that are divisors of $3n^2$;

*assign to $b$ the values $2n+k$;

*calculate for each $b$ the value $a=n(2b-n)/(b-2n)$.


For example, setting $n=1$, we get that the possible values of $k$, i.e. the divisors of $3 \cdot 1^2=3$, are $1$ and $3$. So we can assign to $b$ the values $2\cdot 1 +1=3$ and $2\cdot 1 +3=5$. The corresponding values of $a$ are $(6-1)/(3-2)=5$ and $(10-1)/(5-2)=3$. So the case $n=1$ yields the symmetric pairs $5,3$ and $3,5$. 
As another example, setting $n=6$, we get that the possible values of $k$, i.e. the divisors of $3 \cdot 6^2=108 $, are $1$, $2$, $3$, $4$, $6$, $9$, $12$, $18$, $27$, $36$, $54$ and $108$. Proceeding as above, we can assign to $b$ the values $2\cdot 6 +k$ and calculate the corresponding values of $a$. This yields the $a,b$ pairs $(120,13)$, $(66,14)$, $(48,15)$, $(39,16)$, $(30,18)$, $(24,21)$, plus the other six symmetric pairs. 
You can find here a Wolfram confirmation of the $12$ pairs generated by the case $n=6$.
A: $a^2 +ab+b^2 = c^2$   when  $m>n>0$  are coprime integers and
$$  a=m^2 - n^2 $$
$$ b = 2mn + n^2 $$
$$  c = m^2 + mn + n^2  $$
The phrase Eisenstein Triples is for the very similar $a^2 - ab + b^2 = c^2$
https://en.wikipedia.org/wiki/Eisenstein_triple
https://en.wikipedia.org/wiki/Integer_triangle#Integer_triangles_with_a_120%C2%B0_angle
