Squarefree parts of integers of the form $xy(x+2y)(y+2x)$ The motivation for this question comes from Theorem 3.3 of the 1995 paper Tilings of Triangles by M. Laczkovich, which states:

Let $x$ and $y$ be non-zero integers such that $x+2y\neq 0\neq y+2x$. Then there is  a positive integer $k$ such that the equilateral triangle can be dissected into $n=|xy(x+2y)(y+2x)k^2|$ congruent triangles.

I am curious about the realizable squarefree parts of such integers, hereafter $s(x,y)$. Since we focus on the squarefree part, it suffices to consider the case where $\gcd(x,y)=1$.
In the range $|x|,|y|<10,000$, the realizable values of $s(x,y)$ are
$$1, 5, 6, 10, 11, 13, 14, 15, 17, 19, 21, 22, 23, 29, 30, 33,\ldots$$
of which the squarefree integers missing are
$$2, 3, 7, 26, 31, 38, 43, 51, 53,\ldots$$
This sequence is not in OEIS, but I am very far from confident that it is complete; the smallest solution for $s(x,y)=19$ is given by $x=578,y=-225$ and for $s(x,y)=37$ it is $x=5929,y=648$. However, if there are any squarefree integers missing from the image of $s$ then no compatible OEIS sequences exist, even allowing for an initial $0$ term.
I would be curious to see a proof that specific values like $7$ are not realizable by this function, even if it does not generalize to a complete characterization. Alternatively, pointers to open problems that render solving this very difficult or conjectures that would imply certain results here would also be welcome.

As an aside, here is a sketch of a rough heuristic argument (not a proof!) that the equation $s(x,y)=k$ should have only finitely many solutions for a fixed $k$ with $x$ and $y$ coprime:
Consider the squarefree parts of $x$, $y$, $2x+y$, and $x+2y$, hereafter $a,b,c,d$ respectively. No prime greater than $3$ can be shared among two of $a,b,c,d$ without violating coprimality, so we have finitely many ways to partition the factors of $k$ among the four values times at most $4^4$ ways to assign them all an additional factor of $1,2,3,$ or $6$. So we are seeking a solution to one of a finite number of equation systems of the form $x=\pm a\cdot p^2, y=\pm b\cdot q^2, x+2y=\pm c\cdot r^2, 2x+y=\pm d\cdot s^2$.
Let's estimate the expected number of solutions to one such equation system with $\frac z2 < |x|,|y|\le z$. We have $O(\sqrt{z})$ choices for each of $x$ and $y$ in this range, hence $O(z)$ choices for the pair. For each such choice of $x$ and $y$, if we model an integer $n$ of having an $O(1/n)$ chance to be square, and assume independence, we'll get an $O(1/z^2)$ chance of success (ie that both $c(x+2y)$ and $d(2x+y)$ ended up being perfect squares) for every such $x,y$ we pick. Adding up, we have (heuristically) an expected $O(1/z)$ total solutions for valid $x,y$ in this range. If we sum this over the ranges $[\frac z2=1,z=2],[\frac z2=2,z=4],[4,8],[8,16],\ldots$ we find that the sum converges. Hence, the expected number of solutions across all the equations, and so the expected number of solutions to $s(x,y)=k$, is finite, though the expected size of a solution if one exists is infinite. This accords fairly well with the data so far, including the isolated very large solutions.
(Of course, there might well be substantial correlations or patterns in the factors of $x+2y$ and $2x+y$ that make this model of the terms as independent random draws extremely
bad! But absent reasons to suspect such patterns in one direction or another, I'd guess this is roughly the behavior of this function.)
Cross-posted from Math StackExchange here.
 A: If a given squarefree integer $s$ is given, determining whether there is a pair of integers $(x,y)$ for which $s(x,y) = s$ is a problem about rational points on elliptic curves. In particular, there is a bijection to solutions to
$$
  xy(x+2y)(y+2x) = sz^{2}
$$
and rational points on the elliptic curve $E_{s} : Y^{2} = X^{3} + 5sX^{2} + 4s^{2} X$, where $X = \frac{2sx}{y}$ and $Y = \frac{2s^{2} z}{y^{2}}$. This elliptic curve always has four points on it (namely $(-4s,0)$, $(-s,0)$, $(0,0)$ and the point at infinity) that form a subgroup isomorphic to $\mathbb{Z}/2\mathbb{Z} \times \mathbb{Z}/2\mathbb{Z}$. If $s = 1$, the torsion subgroup is larger, but if $s > 1$, the existence of an additional rational point (and hence a value of $(x,y)$ with $s(x,y) = s$ with $x + 2y \ne 0$ and $y + 2x \ne 0$) is equivalent to at least one of $E_{s}$ or $E_{-s}$ having positive rank.
These curves each have rank zero for $s \in \{ 2, 3, 7, 26, 31, 38, 51 \}$ and hence it is impossible to solve $s(x,y) = s$ for these values of $s$. For $s = 43$, $E_{-43}$ has positive rank and one rational point has $X = \frac{926171883}{29604481}$. This produces $x = -21538881$ and $y = 59208962$ for which
$$
  xy(x+2y)(y+2x) = -43 \cdot 215287751668080^{2}.
$$
A solution for $s = 53$ can be found with $x = -28322$ and $y = 223925$ which gives
$$
  xy(x+2y)(y+2x) = -53 \cdot 2897871340^{2}.
$$
The most complicated solution I found (within the range of Cremona's database of elliptic curves) was with $s = 139$ where it is necessary to take $x = -4659047399556225$ and $y = 15599782505124098$.
It is worth noting that the heuristic you propose at the end of your question is actually not correct, and if $E_{s}$ or $E_{-s}$ has positive rank, then there will be infinitely many solutions to $s(x,y) = s$, but it's also the case that elliptic curves "just barely" have infinitely many rational points.
Finally, all of these elliptic curves are quadratic twists of $E_{1} : Y^{2} = X^{3} + 5X^{2} + 4X$, and for this reason, there is the possibility of an analogue of Tunnell's theorem that would tell you (conditional on the Birch and Swinnerton-Dyer conjecture) for precisely which values of $s$ there is a solution and for which there isn't.
