Let me interchange $x$ and $y$ for my own convenience. We want to solve $$x^2 - 5y^2 = \pm 4$$ over the integers.
Solving these equations corresponds to finding the elements of norm $\pm 4$ in the quadratic integer ring $\mathbf{Z}[\sqrt{5}]$, where the norm is the function given by $$N(x+\sqrt{5}y) = (x+\sqrt{5}y)(x-\sqrt{5}y) = x^2 - 5y^2.$$
Finding these elements is an exercise in algebraic number theory.
The real quadratic number field $\mathbf{Q}(\sqrt{5})$ has $\mathbf{Z}[\omega]$ with $\omega = (1+\sqrt{5})/2$ as its ring of integers, and $\mathbf{Z}[\sqrt{5}]$ is a subring of this. The field norm on $\mathbf{Q}(\sqrt{5})$ agrees with the norm given above for elements of $\mathbf{Z}[\sqrt{5}]$.
Lemma I.7.2 in Neukirch's Algebraic Number Theory yields that up to multiplication by units in $\mathbf{Z}[\omega]$, there are only finitely many elements of a given norm in $\mathbf{Z}[\omega]$. Since $\mathbf{Z}[\sqrt{5}] \subset \mathbf{Z}[\omega]$ and the norms agree, up to multiplication by units in $\mathbf{Z}[\omega]$ there are only finitely many elements of norm $4$ in $\mathbf{Z}[\sqrt{5}]$.
By Dirichlet's unit theorem the group of units of $\mathbf{Z}[\omega]$ has rank $1$.
A generator of this group, or a fundamental unit of $\mathbf{Q}(\sqrt{5})$, is given by $$\varepsilon = \frac{1+\sqrt{5}}{2},$$ which has norm $-1$.
Since the norm of an element $\alpha$ is the same as the norm of the principal ideal $(\alpha)$, it is useful to determine the number of ideals of norm $4$ in $\mathbf{Z}[\omega]$.
By this answer to an other question this number is $$\sum_{m|4} \chi(m) = \chi(1) + \chi(2) + \chi(4) = \left(\frac{1}{5}\right) + \left(\frac{2}{5}\right) + \left(\frac{4}{5}\right) = 1 - 1 + 1 = 1.$$
Hence if $\alpha, \beta$ are two elements of norm $4$, then $(\alpha) = (\beta)$, so $\beta = u\alpha$ for a unit $u$. That is, up to multiplication by units in $\mathbf{Z}[\omega]$ there is only one element $\alpha$ of norm $4$.
Take $\alpha = 2$; then all the elements of norm $4$ in $\mathbf{Z}[\omega]$ are given by $2\varepsilon^n$, for integer $n$.
But since $2\mathbf{Z}[\omega] \subset \mathbf{Z}[\sqrt{5}]$, all of these elements in fact belong to $\mathbf{Z}[\sqrt{5}]$.
Hence all the solutions to the original equation are the $(x_n, y_n)$ given by $2\varepsilon^n = x_n + \sqrt{5}y_n$.
From the identity $\varphi^n = \frac{L_n + \sqrt{5}F_n}{2}$ of real numbers for nonnegative $n$ mentioned at the end of this section of the Wikipedia article on Lucas numbers it follows that $$2\varepsilon^n = L_n + \sqrt{5}F_n$$ for nonnegative $n$.
For negative $n$ you get extra solutions like $(1,-1)$ and $(-3,1)$, but you could have predicted those from the beginning: if $(x,y)$ is a solution, then so are $(-x,y)$, $(x,-y)$ and $(-x,-y)$.
I should mention that with SAGE you can do calculations in $\mathbf{Q}(\sqrt{5})$,
K.<s> = QuadraticField(5)
eps = (1+s)/2 # = K.units()[0]
for n in range(0,15):
print 2*eps^n
and also with Fibonacci and Lucas numbers:
for n in range(0,15):
print (fibonacci(n), lucas_number2(n,1,-1))
These two pieces of code give the same output (up to formatting).
Edit (01/11/14): A more elementary way to see that there is only one ideal of norm 4 in $\mathbf{Z}[\omega]$ is as follows:
The quadratic field $\mathbf{Q}(\sqrt{5})$ has discriminant $5$ and has no complex embeddings; hence by this inequality we have $N(I) \geq N(x)/\sqrt{5}$ for any ideal $I$ and element $x \in I$.
Since $\mathbf{Z}[\omega]$ is a Dedekind domain we have unique factorization of ideals into primes.
For a prime $\mathfrak{p} \subset \mathbf{Z}[\omega]$ lying over $p$ we get $N(\mathfrak{p}) \geq p^2/\sqrt{5}$.
Since $p^2/\sqrt{5} > 4$ for $p > 2$, the primes of norm at most $4$ must lie over $2$. The minimal polynomial $X^2 - X - 1$ of $\omega$ is irreducible mod $2$, so $2$ is inert in $\mathbf{Z}[\omega]$ by the Kummer-Dedekind theorem. That is, $(2)$ is the only prime with norm at most $4$, and its norm is exactly $4$. By unique factorization into primes and multiplicativity of the norm, $(2)$ is the only ideal of norm $4$ in $\mathbf{Z}[\omega]$.
Let me advertise the following version of the Lifting The Exponent Lemma. As setup, let $K$ be any characteristic zero field with a nonarchimedean $p$-adic valuation $v_p: K \rightarrow \mathbb R \cup \{\infty\}$ with $v_p(p)=1$ (i.e. the restriction of $v_p$ to $\mathbb Q$ is the standard $p$-adic valuation). Then for any two elements $x,y \in K$, and all $n \in \mathbb N$, the following holds:
LTE (Lifting The Exponent) Lemma
If $v_p (x-y) > v_p(x) + \dfrac{1}{p-1}$ then $v_p(x^n -y^n) = (n-1)\cdot v_p(x) + v_p(x-y) +v_p(n)$.
[Special Case $v_p(x) = 0$: If $v_p (x-y) > \dfrac{1}{p-1}$ then $v_p(x^n -y^n) = v_p(x-y) +v_p(n)$.]
A proof goes like this: Via multiplying through with $y^{-n}$ (note that the condition implies $v_p(x) = v_p(y)$), WLOG we are in the special case with $y=1$ i.e. $x = 1+r$ for some $r \in K$ with $v_p(r) > \dfrac{1}{p-1}$.
Then $x^n-y^n = (1+r)^n – 1 = \sum_{k=1}^n \binom{n}{k} r^k$ and it suffices to consider the two cases
Case A. $n=p$, where $v_p (\binom{p}{k} r^k) \begin{cases} =v_p(r) + 1 \qquad \quad \text{ for } k=1 \\ \ge kv_p(r) +1 \qquad \text{ for } 2 \le k \le p-1 \\ = pv_p(r) \qquad \qquad\text{ for }k=p \end{cases}$
By the assumption $v_p(r) > \dfrac{1}{p-1}$ the first one is smaller than all the others, hence $=v_p((1+r)^n – 1)$.
Case B. $\gcd(p,n) =1$, where $v_p (\binom{n}{k} r^k) \begin{cases} = v_p(r) \text{ for } k=1 \\ \ge kv_p(r) \text{ for } 2 \le k \le n \end{cases}$
and the assumption ($v_p(r) >0$ would suffice) again guarantees the first one is smaller than all the others, and hence $=v_p((1+r)^n – 1)$. $\square$
This implies both your assertions about the $P_n = \dfrac{\phi^n - \psi^n}{2\sqrt{2}}$, as follows.
First note that $v_p(\phi) = v_p(\psi) = 0$ for all $p$ (e.g. since $\phi, \psi$ are roots of $X^2+2X-1$), so we are in the "Special Case".
In case $p$ odd we have $v_p(2\sqrt{2}) = 0$, hence $v_p(P_n) = v_p(\phi^n-\psi^n)$ for all $n$. Now your assertion (II) is the lemma applied to $x=\phi^n, y=\psi^n$ (and $n=m$, sorry), noting that the hypothesis of that assertion is $v_p(P_n) \ge 1 >\dfrac{1}{p-1}$.
In case $p=2$ we have $v_2(2\sqrt{2}) = \frac32$, hence $0=v_2(P_1) = v_2(\phi-\psi) – v_2(2\sqrt{2})$ i.e. $v_2(\phi-\psi) = \frac32 > \dfrac{1}{2-1}$. So $$v_p(P_n) = v_p(\phi^n - \psi^n) - \frac32 \stackrel{LTE}= \frac32 + v_p(n) -\frac32 = v_p(n)$$ for all $n$.
Now, your proposed proof works just as well, and is more conceptual by using the $p$-adic logarithm. In fact, I recommended that line of thought here. However, you have to be careful about two things (which, furthermore, might interfere with each other):
Normalizing the extended valuation $v_{\mathfrak p}$, which is $e\cdot v_p$ for the ramification index of the field: That might cause some case distinctions which ultimately don't matter but make a proof cumbersome.
The two facts $$\log(x^n) = n \log(x)$$ and $$v_p(\log(x)) = v_p(x-1)$$ are not as trivial as they are beautiful. In fact, the first one does hold as long as $v_p(x) > 0$, but the second one is valid in general only if $v_p(x-1) > \dfrac{1}{p-1}$ (as it should, because otherwise why would that assumption be needed in our version of LTE). Cf. Corollary 8.6 and Theorem 8.7 in K. Conrad's superb summary of $p$-adic series. This ultimately relies on the same computations as the straightforward algebra above.
In fact, the number $\dfrac{1}{p-1}$ (which annoyingly is $<1$ for odd $p$ but $=1$ for $p=2$) occurs again and again in this theory, as it is the break point where the two operations "raising to the $p$-th power" and "multiplying with $p$" switch places (as to which one has more effect on the valuation).
Best Answer
If $D$ is an integer, then $5p_n^2-4$ has to be a perfect square. It is known that if either $5N^2-4$ or $5N^2+4$ is a perfect square, then $N$ is a Fibonacci number. Since one can prove by induction that $5F_n^2+4(-1)^n=L_n^2$, one can say that $p_n$ is an odd-indexed Fibonacci number $F_{2m+1}$, and that $\sqrt{5p_n^2-4}$ is an odd-indexed Lucas number $L_{2m+1}$.
If we choose $+$, then we have $$D=\frac{3p_n+\sqrt{5p_n^2-4}}{2}\cdot \frac{\Delta}{p_n}=\bigg(\frac 32+\frac{\sqrt{5p_n^2-4}}{2p_n}\bigg)\Delta\gt \Delta$$ which contradicts $D\le\Delta$.
So, we have to choose $-$, and we finally get $$\begin{align}D&=\frac{3p_n-\sqrt{5p_n^2-4}}{2}\cdot \frac{\Delta}{p_n} \\\\&=\frac{3F_{2m+1}-\sqrt{5F_{2m+1}^2-4}}{2}\times (p_1\cdots p_{n-1}) \\\\&=\frac{3F_{2m+1}-L_{2m+1}}{2}\times (p_1\cdots p_{n-1}) \\\\&\color{red}=\frac{3F_{2m+1}-(2F_{2m+2}-F_{2m+1})}{2}\times (p_1\cdots p_{n-1}) \\\\&=(2F_{2m+1}-F_{2m+2})\times (p_1\cdots p_{n-1}) \\\\&=\bigg(2F_{2m+1}-(F_{2m+1}+F_{2m})\bigg)\times (p_1\cdots p_{n-1}) \\\\&=(F_{2m+1}-F_{2m})\times (p_1\cdots p_{n-1}) \\\\&=F_{2m-1}\times (p_1\cdots p_{n-1})\end{align}$$ where the fact that $L_n=2F_{n+1}-F_n$ (which can be proved by induction) was used in the equality in red.