Characterizing integer solutions to $a^2 \mid \bigl((b^2)^2 + (b^2+1)^2\bigr)$ I’m looking to characterize all integers $a$ and $b$ satisfying $$a^2 \mid \bigl((b^2)^2 + (b^2+1)^2\bigr), \qquad a > b \ge 1.$$
Brute force searches have so far turned up the two solutions $(a,b) = (13,7)$ and $(a,b)=(9601,5441)$. I’m guessing there are an infinite number, but I’d like to reduce the numbers I have to check. Obviously, $a$ must be odd and the sum of two squares, so I’m already taking that into account.
Any thoughts about what the next solution might be, or how to attack the characterization of any/all solutions?
If $b$ is odd, then there exists an integer $c$ such that $a=b+2c$. Substituting, we find $a \parallel (32c^4+8c^2+1) = \bigl((4c^2)^2+(4c^2+1)^2\bigr)$, and then $a^2$ divides a related quotient. Based on that, I feel like there might be a Vieta-jumping solution… but I still can’t quite find it.
 A: We know that $(b)^2$ and $(b+1)^2$ are of opposite parity  and this makes the sum of their squares odd so $(a)^2$ must be odd.  Brute force, even limited to $15$ digits, shows solutions not in the OP question.
Here are the first $3$ $\space (a,b)$-values where
$\quad 2 b^4 + 2 b^2 + 1 \le 987654321054321.$
$$
(13 ,7 )\quad 
(533 ,162 )\quad 
(1585 ,914 )\quad 
$$
A sum need not be a perfect square but only contain a perfect square as one of its factors.
Here: $\quad f(7)/13^2=29\\
 f(162)/533^2=4849\\
 f(914)/1585^2=555593$
and here is the program that generated them
  10 print "Enter H1" : input h1 : print
  20 for a1 = 3 to h1
  30    for b1 = 2 to a1-1
  40    if 2*b1^4+2*b1^2+1 <= 987654321054321
  50       d1 = (2*b1^4+2*b1^2+1)/a1^2
  60       if d1 = int(d1)
  70          print "(" a1 "," b1 ")\quad "
  80       endif
  90     endif
  100    next b1
  110 next a1

Note: In an earlier edit of this answer I overlooked rounding errors that made most of my solutions invalid where
$\quad 2 b^4 + 2 b^2 + 1\ge 10^{15},\quad $ I have limited my "solutions" here but, in any language with arbitrary precision, this restriction should not be needed.
A: If $a^2 | 2b^4+2b^2+1$, then for any prime factor $p$ of $a$:
$ p | 2b^4+2b^2+1$
We can find the possible values of $b$ modulo $p$. Values only occur when $p = 1$ mod 4, and there are either 0,2, or 4 possible values,
as shown in the following analysis:
Set $ X = b^2$, then $2X^2+2X+1 = 0$ mod $p$
Multiply by 2 (clearly p is odd) and rearrange to get  $(2X+1)^2 = -1$ mod $p$
Find $w$, a fourth root of unity mod $p$ (this is only possible if $p = 1$ mod 4), then
$2X+1 = \pm w$ mod $p$
rearranging, $b^2 = X = (\pm w-1)(\frac{p+1}{2})$ mod $p$
$b$ has 0, 2 or 4 solutions mod p, depending on whether there are solutions to the above equation.
This suggests an efficient search algorithm:

*

*for primes $p = 1$ mod 4, find and store the 2 or 4 values of $b$ mod $p$


*for $b$ in a range of values, find the primes $p$ that satisfy $ p | 2b^4+2b^2+1$,


*from those, find the ones where $ p^{2n} | 2b^4+2b^2+1$


*using the values $p,b$ found in step 3, construct values $a$ (primes or composite) satisfying $ a^2 | 2b^4+2b^2+1 $
I have implemented this a python 3 program, searching primes up to $10^8$ and values of $b$ up to $10^8$.
The only solutions found have been those already identified by @Tomita
The Python 3 program is on github
https://github.com/armchaircaver/Characterizing-integer-solutions-to-a-squared-divides-polynomial-in-b
Addition 30 Aug 2021:
There are techniques to find a square root modulo a prime power (see for example https://www.johndcook.com/blog/quadratic_congruences/ ), which i have used to find values of b modulo $p^2$ satisfying $p^2 | 2b^4+2b^2+1$ directly, making the search even more efficient, as there are 2 or 4 values of $b$ mod $p^2$. The more efficient version is "stackexchange puzzle sieve p squared.py", and has been added to the same github repository as above.
Using the more efficient version I have searched primes up to $10^8$ and values of $b$ up to $10^9$ with no new solutions found.
A: This problem in a way is related to equation $x^2+(x+1)^2=y^2$ which has infinite solution.There is a proof in this site.I will try to show this relation by an examples:
$119^2+120^2=13^4$
$119=7\times 17$
$(7\times17)^2+(7\times17+1)=13^4$$\space\space\space\space\space\space(1)$
and we see that:
$7^4+(7^2+1)^2=29\times 13$
Relation (1) is symmetric for 7 and 17 so relation $b^4+(b^2+1)^2= km^2$ must also be true for 17, we check this:
$17^4+(17^2+1)^2=86021=509\times 13^2$
that is if $b$ is a primes multiple of a primitive solution, that prime is also a solution. $x$ and $x+1$ make following sets:
$x\in\{3, 119, 4059, 137903\cdot\cdot\cdot\}$
$(x+1)\in\{4, 120, 4060, 137904\cdot\cdot\cdot\}$
