A007781(n) = (n+1)^(n+1) - n^n.
A007781(6n+1) is not squarefree for n > 0. a(n) is the largest square divisor of
A007781(6n+1). All terms belong to
A003215 Hex (or centered hexagonal) numbers: 3n(n+1)+1 (crystal ball sequence for hexagonal lattice). It appears that a(n) =
A003215(2n) = 6n(2n+1)+1.
A007781(6n+1)/
A003215(2n)^2 = ((6n+2)^(6n+2)-(6n+1)^(6n+1))/(6n(2n+1)+1)^2 = {44193, 2904899682603, 6378521938392937343349, 128847538453506016002947264859159, 13183819636551142123977274666051092130410345, ...}. Prime terms of a(n) belong to
A002407. Factorizations of the terms of a(n) are {19, 61, 127, 7*31, 331, 7*67, 631, 19*43, 13*79, 13*97, 7*7*31, 1801, 7*7*43, 2437, 2791, 3169, 3571, 7*571, 4447, 7*19*37, 5419, 13*457, 13*499, 7067, 7*1093, 8269, 7*19*67, 61*157, 10267, 79*139, ...}. All prime factors of a(n) are of the form 6k+1.
