What are the triangle numbers the are squares of other triangle numbers? I have found $1^2=1$ and $6^2=36$, but other than these examples I can't find any other triangle numbers that are squares of other triangle numbers, and I used a program to check this idea into the thousands.
Finding triangle numbers that are squares of other triangle numbers corresponds to finding integers $n$ and $k$ such that $n(n+1)/2=[k(k+1)/2]^2$, or such that
$$2n(n+1)=k^2(k+1)^2 .$$
I believe the only positive integer solutions to this equation are $(1,1)$ and $(3,8)$, but I don't know how to prove it. Are there any others?
The only triangular numbers $T_n = \frac{1}{2} n (n + 1)$ that are squares of triangular numbers are the two you've found, namely $T_1 = 1$ and $T_8 = 36 = 6^2 = T_3^2$.
As you've pointed out, this is equivalent to finding the positive integer solutions to the Diophantine equation $$2 n (n + 1) = k^2 (k + 1)^2 .$$ (Incidentally, if $(k, n)$ is a solution, so are $(k, - n - 1), (- k - 1, n), (- k - 1, - n - 1)$, so one can derive all of the solutions from the positive ones and the cases where $n = 0$ or $k = 0$.)
I wasn't able to find an elementary solution, and I would be grateful to see one. On the other hand, this solution ties the question with a beautiful topic mathematics and exploits some modern machinery (the main reference for the technique here was published in 1996).
The Maple routine
algcurves:-Weierstrassform()gives that the coordinate change $$k = -\frac{18 (y - 4)}{s(x)}, \quad n = \frac{-216 (3 x - 2) (y - 4)}{s(x)^2} + \frac{6 (3 x + 10)}{s(x)} ,$$ where $s(x) := 9 x^2 - 12 x - 68$, is a birational equivalence between the rational quartic curve defined by the Diophantine equation and the elliptic curve $E$ defined by $$y^2 = x^3 - \frac{28}{3} x + \frac{160}{27} .$$ This curve has rank $1$ with free generator $G := (-\frac{1}{3}, 3)$, and its torsion subgroup $\operatorname{Tor}(E(\Bbb Q)) \cong \Bbb Z_2 \times \Bbb Z_2$ with generators $Q_1 := (\frac{8}{3}, 0)$, $Q_2 := (\frac{2}{3}, 0)$ ($Q_1 + Q_2 = (-\frac{10}{3}, 0)$). (Here $+$ refers to the group operation on $E$.) The upshot of this is that every rational point on $E$ has the form $m G + T$ for some integer $m$ and $T \in \operatorname{Tor}(E(\Bbb Q))$, and via the inverse of the first coordinate transformation above these points can be mapped to (all of) the rational solutions to the Diophantine equation; of course, the integer solutions we seek are a subset of these.There is a powerful machinery for finding all of the integer points on such curves, and it amounts to finding an upper bound $M$ for the coefficient $|m|$ for elements $m G + T$ that correspond to integer solutions. With such a bound in hand, it's then a matter of checking whether each of the finitely many points $m G + T$, $|m| \leq M$ corresponds to an integer solution.
Since our original equation is quartic, our situation is covered by the article of N. Tzanakis cited below. In particular, using the coordinate transformation $k = U, n = (V - 1) / 2$, our equation assumes the form $$V^2 = 2 U^4 + 4 U^3 + 2 U^2 + 1 ,$$ which is the form required by the method in that article. (NB a priori this transformation may introduce new integral points, but it will not transform integer solutions to noninteger solutions.)
Following the method detailed in the article, which critically involves computing so-called elliptic logarithms and making reasonable choices for estimates at intermediate steps initially gives an upper bound $M = 8 \cdot 10^{24}$, which gives far too many possibilities to check with even a fast computer. Applying the integral version of the so-called LLL reduction drastically improves the bound in a single step, to $M \leq 17$, and applying it again improves it to $M \leq 12$ (additional iterations don't seem to improve this bound). At this point, there are only $4 \cdot [2(12) + 1] = 100$ points to check, and a CAS can carry out this procedure quickly. Doing so gives exactly 12 integer solutions $(k, n)$ of our original equation, namely, \begin{gather} \mathcal{O} \leftrightarrow ( 0, 0), \quad - G + Q_2 \leftrightarrow \color{#bf0000}{( 1, 1)}, \quad G + Q_1 + Q_2 \leftrightarrow \color{#bf0000}{( 3, 8)}, \end{gather} and the 9 solutions obtained by applying the mentioned symmetries to these three.
"Making reasonable choices" sweeps under the rug a good deal of computational work, and it's difficult to summarize meaningfully the intermediate steps of the computation. (But I can attempt a gross outline if there's any demand for it.) Perhaps it suffices to say (1) using a CAS is virtually necessary here, and (2) the open-source software SAGE has excellent support for elliptic curve computations and was indispensable in carrying out this computation.
To give a flavor for some of the data used in the calculate of the bound, some invariants of the elliptic curve $E$ used in the mention computations: The affine transformation $x = X - \frac{1}{3}, y = Y$ brings $E$ to the minimal form $Y^2 = X^3 - X^2 - 9 X + 9$. Its discriminant is $\Delta = 2^{12} \cdot 3^2$, its conductor is $f = 192 = 2^6 \cdot 3$ (which lets us quickly identify $E$ as elliptic curve $192\textrm{a}2$ in John Cremona's tables of curves with small conductor), and its $j$-invariant is $j = \frac{2^6 \cdot 7^3}{3^2}$.
Edit After some further searching, it seems that this result appears first in a (French-language) 1946 paper of Ljunggren (I couldn't find any copy online), which gives rather difficult proof by considering the field extension $\Bbb Q[\sqrt[4]{2}]$ of $\Bbb Q$. In 1965 Cassels (I can't find an ungated copy online) gave a simpler proof related to realization of an elliptic curve as a pair of quadrics in $3$-space.