How to prove each element of the following sequence is a perfect square?

As Umberto P. commented, it is sufficient to prove by induction that $$a_{n+2}=b_n^2\tag1$$ where $b_{n+2}=4b_{n+1}-b_{n}$ with $b_0=1, b_1=5$.

The base case : $a_2=1=b_0^2$ and $a_3=25=b_1^2$

Suppose that $a_{n+2}=b_n^2$ and $a_{n+1}=b_{n-1}^2$.

Then, using that $$b_n^2-4b_nb_{n-1}+b_{n-1}^2=6\tag2$$ we get $$\begin{align}a_{n+3}&=14a_{n+2}-a_{n+1}+12\\&=14b_{n}^2-b_{n-1}^2+12\\&=16b_n^2-8b_nb_{n-1}+b_{n-1}^2-2b_n^2+8b_nb_{n-1}-2b_{n-1}^2+12\\&=(4b_n-b_{n-1})^2-2(b_n^2-4b_nb_{n-1}+b_{n-1}^2-6)\\&=b_{n+1}^2\qquad\blacksquare\end{align}$$

Let us prove $(2)$ by induction.

The base case : $b_1^2-4b_1b_0+b_0^2=25-20+1=6$

Suppose that $b_n^2-4b_nb_{n-1}+b_{n-1}^2=6$. Then,

$$\begin{align}b_{n+1}^2-4b_{n+1}b_{n}+b_{n}^2&=(4b_{n}-b_{n-1})^2-4(4b_{n}-b_{n-1})b_n+b_n^2\\&=16b_{n}^2-8b_{n}b_{n-1}+b_{n-1}^2-16b_n^2+4b_nb_{n-1}+b_n^2\\&=b_n^2-4b_nb_{n-1}+b_{n-1}^2\\&=6\qquad\blacksquare\end{align}$$

I don't know a simpler way than breaking down the generating function: let $b_n=a_n+1$; then $b_{n+2}$ $= a_{n+2}+1$ $= 14a_{n+1}-a_n+12+1$ $=14(b_{n+1}-1)-(b_n-1)+13$ $= 14b_{n+1}-b_n$, so the recurrence relation for $b_n$ is homogeneous. The characteristic equation is $x^2-14x+1=0$ with roots $7\pm4\sqrt{3}$, which are themselves the squares of $2\pm\sqrt{3}$; note that $(2+\sqrt{3})(2-\sqrt{3})=1$. This isn't coincidence: if you look at the equation $c_{n+2}=4c_{n+1}-c_n$ (whose characteristic equation $x^2-4x+1=0$ has $2\pm\sqrt{3}$ for its roots) with initial conditions $c_1=-1$, $c_2=1$ you'll find that you have $a_n=c_n^2$.

This works because we have $c_n=x(2+\sqrt{3})^n+y(2-\sqrt{3})^n$, with $x=\frac12(3\sqrt{3}-5)$ and $y=-\frac12(3\sqrt{3}+5)$and so when we square $c_n$ we get a result of the form $x^2(2+\sqrt{3})^{2n}+y^2(2-\sqrt{3})^{2n}+xy\left((2+\sqrt{3})(2-\sqrt{3})\right)^n$ — but $xy$ is $-1$ and the last term in parenthesesis just $1$ for all $n$, so $c_n^2+1$ is of the form $sz^n+t\bar{z}^n$ (where $\bar{z}$ denotes the conjugate with respect to $\sqrt{3}$), and so it satisfies a second-degree homogeneous linear equation with integer coefficients — specifically, the equation for $b_n$.

The key here is that the (conjugate) roots of the equation have norm $1$ in $\mathbb{Q}(\sqrt{3})$, so the cross term will always be $1$; there's no exponential growth there, so we can offset the square by a constant factor to get another homogeneous linear equation. The same trick would work for e.g. an equation like $x^2-66x+1=0$ (i.e., $c_{n+2}=66c_{n+1}-c_n$) with conjugate roots $z,\bar{z}=33\pm8\sqrt{17}$ that satisfy $z\bar{z}=1$; with the right initial conditions you'll find exactly the same phenomenon that the solutions of $d_{n+2}=4354d_{n+1}-d_n$ (whose roots are $z^2$ and $\bar{z}^2$) satisfy $d_n=c_n^2-1$.

The scenic route.

This begins with Umberto's comment that $ x_{n+2} - 4 x_{n+1} + x_n = 0, $ so I knew I wanted a 2 by 2 matrix of trace $4$ and determinant $1,$ with preference for one with the two diagonal terms equal, because those serve as automorphism matrices for Pell type quadratic forms. I thought of $$ \left( \begin{array}{rr} 2 & 3 \\ 1 & 2 \end{array} \right) $$ and began thinking of $x^2 - 3 y^2.$ Looking at the first few items, $1,25,361,$ in each case we get triple a square just by adding $2,$ as in $3,27,363.$ So, $x^2 + 2 = 3 y^2,$ or $x^2 - 3 y^2 = -2.$

Take $x^2 - 3 y^2$ as a quadratic form, and solve for $x^2 - 3 y^2 = -2.$ The "automorphism matrix" says that $$ \left( \begin{array}{r} x_{n+1} \\ y_{n+1} \end{array} \right) = \left( \begin{array}{rr} 2 & 3 \\ 1 & 2 \end{array} \right) \left( \begin{array}{r} x_{n} \\ y_{n} \end{array} \right), $$ or, from Cayley Hamilton, $$ x_{n+2} - 4 x_{n+1} + x_n = 0, $$ $$ y_{n+2} - 4 y_{n+1} + y_n = 0. $$

If we simply multiply out the entries we get $$ \left( \begin{array}{c} x_{n+1}^2 \\ x_{n+1} y_{n+1} \\ y_{n+1}^2 \end{array} \right) = \left( \begin{array}{rrr} 4 & 12 & 9 \\ 2 & 7 & 6 \\ 1 & 4 & 4 \end{array} \right) \left( \begin{array}{c} x_{n}^2 \\ x_{n} y_{n} \\ y_{n}^2 \end{array} \right) $$ This 3 by 3 matrix is sort of famous, as we have taken a 2 by 2 matrix in the modular group and produced a 3 by 3 image in $SL_3 \mathbb Z;$ the map is an isomorphism. You can find this construction on page 58 of Zagier or page 23 of Magnus. Actually, Magnus and Zagier use transposes of each other's matrices, can't have everything. With the letters permuted, page 301 of Cassels.

The coefficient matrix has determinant $1.$ In fact, the characteristic polynomial is $$ \lambda^3 - 15 \lambda^2 + 15 \lambda - 1 = (\lambda - 1)(\lambda^2 - 14 \lambda + 1). $$

The Cayley-Hamilton Theorem says that $$ x_{n+3}^2 - 15 x_{n+2}^2 + 15 x_{n+1}^2 - x_{n}^2 = 0. $$

NOW, suppose I have any sequence $z_n$ with $$ z_{n+3} - 15 z_{n+2} + 15 z_{n+1} - z_{n} = 0, $$ and I define $$ w_n = z_{n+2} - 14 z_{n+1} + z_n. $$ THEN $$ w_{n+1} = z_{n+3} - 14 z_{n+2} + z_{n+1}, $$ AND $$ w_{n+1} - w_n = z_{n+3} - 15 z_{n+2} + 15 z_{n+1} - z_{n} = 0,$$ or $w_n$ is constant! Taking $z_n = x_n^2,$ we get $$ x_{n+2}^2 - 14 x_{n+1}^2 + x_{n}^2 = \mbox{CONSTANT}. $$ Check with the triple $1,5,19$ beginning $x_n,$ we find $$ x_{n+2}^2 - 14 x_{n+1}^2 + x_{n}^2 = 12. $$

It is not difficult to confirm that, given an element of the modular group (our "automorphism matrix") $$ \left( \begin{array}{rr} \alpha & \beta \\ \gamma & \delta \end{array} \right), $$ so that $$ \alpha \delta - \beta \gamma = 1, $$ then construct $$ M = \left( \begin{array}{ccc} \alpha^2 & 2 \alpha \beta & \beta^2 \\ \alpha \gamma & (\alpha \delta + \beta \gamma) & \beta \delta \\ \gamma^2 & 2 \gamma \delta & \delta^2 \end{array} \right), $$ that the characteristic polynomial of $M$ is $$ \lambda^3 - \left( (\alpha + \delta)^2 - 1 \right) \lambda^2 + \left( (\alpha + \delta)^2 - 1 \right) \lambda - 1 = (\lambda - 1) \left( \lambda^2 - \left( (\alpha + \delta)^2 - 2 \right) \lambda + 1 \right) $$ Part of this is the claim that $1$ is an eigenvalue of $M.$ Here is an eigenvector with eigenvalue $1:$ $$ X = \left( \begin{array}{c} 2 \beta \\ \delta - \alpha \\ - 2 \gamma \end{array} \right), $$ with $MX=X,$ and $X$ is nonzero unless $W = \pm I,$ in which case $M=I$ and all eigenvalues become $1.$

Here is page 301 in Cassels. He switches $\beta, \gamma$ in the 3 by 3 matrix, compared with what I typed above.

Here is a start:

jagy@phobeusjunior:~$ ./Pell_Target_Fundamental
  Automorphism matrix:  
    2   3
    1   2
  Automorphism backwards:  
    2   -3
    -1   2

  2^2 - 3 1^2 = 1

 x^2 - 3 y^2 = -2

Tue Aug 23 11:43:01 PDT 2016

x:  1  y:  1 ratio: 1  SEED   KEEP +- 
x:  5  y:  3 ratio: 1.66667
x:  19  y:  11 ratio: 1.72727
x:  71  y:  41 ratio: 1.73171
x:  265  y:  153 ratio: 1.73203
x:  989  y:  571 ratio: 1.73205
x:  3691  y:  2131 ratio: 1.73205
x:  13775  y:  7953 ratio: 1.73205
x:  51409  y:  29681 ratio: 1.73205
x:  191861  y:  110771 ratio: 1.73205
x:  716035  y:  413403 ratio: 1.73205
x:  2672279  y:  1542841 ratio: 1.73205
x:  9973081  y:  5757961 ratio: 1.73205

Tue Aug 23 11:43:31 PDT 2016

 x^2 - 3 y^2 = -2

jagy@phobeusjunior:~$

Calculations for showing that the mapping from 2 by 2 matrices to 3 by 3, as in Magnus page 23, is an isomorphism.

parisize = 4000000, primelimit = 500509
? m = [ a,b; c,d]
%1 = 
[a b]

[c d]

? m3 = [ a^2, 2 * a * b, b^2; a * c, a * d + b * c, b * d; c^2, 2 * c * d, d^2]
%2 = 
[a^2 2*b*a b^2]

[c*a d*a + c*b d*b]

[c^2 2*d*c d^2]

? matdet(m3) - ( matdet(m) )^3
%24 = 0

? n = [ e,f; g,h]
%3 = 
[e f]

[g h]

? n3 = [ e^2, 2 * e * f, f^2; e * g, e * h + f * g, f * h; g^2, 2 * g * h, h^2]
%4 = 
[e^2 2*f*e f^2]

[g*e h*e + g*f h*f]

[g^2 2*h*g h^2]

? m * n
%5 = 
[e*a + g*b f*a + h*b]

[e*c + g*d f*c + h*d]

? 
? 
? m3 * n3
%6 = 
[e^2*a^2 + 2*g*e*b*a + g^2*b^2 2*f*e*a^2 + (2*h*e + 2*g*f)*b*a + 2*h*g*b^2 f^2*a^2 + 2*h*f*b*a + h^2*b^2]

[(e^2*c + g*e*d)*a + (g*e*c + g^2*d)*b (2*f*e*c + (h*e + g*f)*d)*a + ((h*e + g*f)*c + 2*h*g*d)*b (f^2*c + h*f*d)*a + (h*f*c + h^2*d)*b]

[e^2*c^2 + 2*g*e*d*c + g^2*d^2 2*f*e*c^2 + (2*h*e + 2*g*f)*d*c + 2*h*g*d^2 f^2*c^2 + 2*h*f*d*c + h^2*d^2]

? t = m * n
%7 = 
[e*a + g*b f*a + h*b]

[e*c + g*d f*c + h*d]

? t[1,1]
%8 = e*a + g*b
? t[1,2]
%9 = f*a + h*b
? t[2,1]
%10 = e*c + g*d
? t[2,2]
%11 = f*c + h*d
? 
? u = m3 * n3
%12 = 
[e^2*a^2 + 2*g*e*b*a + g^2*b^2 2*f*e*a^2 + (2*h*e + 2*g*f)*b*a + 2*h*g*b^2 f^2*a^2 + 2*h*f*b*a + h^2*b^2]

[(e^2*c + g*e*d)*a + (g*e*c + g^2*d)*b (2*f*e*c + (h*e + g*f)*d)*a + ((h*e + g*f)*c + 2*h*g*d)*b (f^2*c + h*f*d)*a + (h*f*c + h^2*d)*b]

[e^2*c^2 + 2*g*e*d*c + g^2*d^2 2*f*e*c^2 + (2*h*e + 2*g*f)*d*c + 2*h*g*d^2 f^2*c^2 + 2*h*f*d*c + h^2*d^2]

? u[1,1]
%13 = e^2*a^2 + 2*g*e*b*a + g^2*b^2
? u[1,2]
%14 = 2*f*e*a^2 + (2*h*e + 2*g*f)*b*a + 2*h*g*b^2
? u[1,2] - 2 * t[1,1] * t[1,2]
%15 = 0
? u[1,1] - t[1,1] * t[1,1]
%16 = 0
? u[1,3] - t[1,2] * t[1,2]
%17 = 0
? u[2,1] - t[1,1] * t[2,1]
%18 = 0
? u[2,2] - t[1,1] * t[2,2] - t[1,2] * t[2,1]
%19 = 0
? u[2,3] - t[1,2] * t[2,2]
%20 = 0
? u[3,1] - t[2,1] * t[2,1]
%21 = 0
? u[3,2] - 2 * t[2,1] * t[2,2]
%22 = 0
? u[3,3] -  t[2,2] * t[2,2]
%23 = 0
?