The question would, naturally, be very straight forward if there was a $2xy$ instead of a $4xy$. Then it would simply be a matter of doing:
$$ x^2+2xy+y^2=1\\ (x+y)^2=1\\ \sqrt{(x+y)^2}=\sqrt{1}\\ |x+y|=1 $$
Clearly, there are infinitely many pairs of integers that differ by one so there is an infinite number of integer solutions to $x^2+2xy+y^2=1$. Unfortunately, the same principle does not apply to $x^2+4xy+y^2=1$ where if I attempt to construct a similar proof all I can do is:
$$ x^2+4xy+y^2=1\\ (x+y)^2+2xy=1\\ 2xy=1-(x+y)^2\\ 2xy=(1+x+y)(1-x-y) $$
And, from there, I have no idea how to proceed to complete the proof that there are an infinite number of integer solutions. I'm wondering whether I'm approaching the question entirely in the wrong way or if I am simply missing something. A nudge in the right direction would be much appreciated!
Hint $\ $ Consider $\,f(x) = x^2 + 4y\ x + y^2\!-\!1\,$ as a quadratic in $\,x,$ where $y$ is constant. By Vieta its roots $\,x,x'$ satisfy $\ x+x' = -4y.\,$ Thus if $\,x\,$ is a root then so too is $\,x' = -4y-x.$
This yields a reflection symmetry $\ \, (x,y) \mapsto (-4y-x,\,y)\,$ on the solution space. Composing
this with the reflection symmetry $\,(x,y)\mapsto (y,x)\,$ yields the map $\,(x,y)\mapsto (-4x-y,x)$
which, iterated starting at solution $(1,0),\,$ yields infinitely many solutions
$$ (1,0),\ (-4,1),\ (15,-4),\ (-56,15),\ (209,-56),\ (-780,209),\ (2911,-78),\ \ldots$$
Sequence $\, 0,1,4,15,56,209,\ldots$ satsifies the recursion $\,f_{n+2} = 4 f_{n+1} - f_{n}\,$ as is easily derived.
This can be transformed to the Pell equation $\ X^2\! - 3 Y^2 = 1\ $ and studied using standard results on Pell equations. See also the comments on this sequence at OEIS sequence A001353.