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.)