[proofplan]
The proof converts the continued-fraction periodicity of $\sqrt{d}$ into an algebraic identity. First, Galois's theorem applies to $a_0 + \sqrt{d}$ (whose conjugate lies in $(-1, 0)$) to show that the expansion of $\sqrt{d}$ has the special form $[a_0, \overline{a_1, \ldots, a_n}]$. Evaluating $\sqrt{d}$ as the $(n+1)$-st truncation yields a Möbius identity relating $\sqrt{d}$ to the convergents $p_n/q_n$ and $p_{n-1}/q_{n-1}$. Equating rational and irrational parts in this identity produces the linear system from which the Pell identity $p_{n-1}^2 - d q_{n-1}^2 = (-1)^n$ drops out. For even $n$, $(p_{n-1}, q_{n-1})$ solves Pell's equation; for odd $n$, pass to the doubled period. Iterating the period at integer multiples $kn$ with $kn$ even yields infinitely many solutions.
[/proofplan]
[step:Show $\sqrt{d}$ has a continued fraction of the form $[a_0, \overline{a_1, \ldots, a_n}]$]
Let $a_0 := \lfloor \sqrt{d} \rfloor \geq 1$ (valid since $d \geq 2$; the case $d = 1$ is excluded because $d$ is not a perfect square). Consider
\begin{align*}
\phi &:= a_0 + \sqrt{d}.
\end{align*}
Then $\phi > a_0 + a_0 = 2a_0 \geq 2 > 1$ (using $\sqrt{d} > a_0$, which holds because $d$ is not a perfect square so $\sqrt{d} \notin \mathbb{Z}$, so $\sqrt{d} > a_0$ strictly). The Galois conjugate of $\phi$ is
\begin{align*}
\phi' &= a_0 - \sqrt{d} \in (a_0 - (a_0 + 1), a_0 - a_0) = (-1, 0),
\end{align*}
where the bound $a_0 < \sqrt{d} < a_0 + 1$ comes from $a_0 = \lfloor \sqrt{d} \rfloor$.
Hence $\phi = a_0 + \sqrt{d}$ satisfies the hypotheses of [Galois — Purely Periodic Characterisation](/theorems/1764): $\phi > 1$ and $-1 < \phi' < 0$. Therefore its continued fraction is purely periodic:
\begin{align*}
a_0 + \sqrt{d} &= [\overline{b_0, b_1, \ldots, b_{n-1}}]
\end{align*}
for some period $n \geq 1$ and partial quotients $b_j$. The first partial quotient is $b_0 = \lfloor a_0 + \sqrt{d} \rfloor = 2 a_0$. Since $\sqrt{d} = (a_0 + \sqrt{d}) - a_0$, subtracting $a_0$ from the first partial quotient gives the expansion of $\sqrt{d}$:
\begin{align*}
\sqrt{d} &= [a_0, \overline{b_1, b_2, \ldots, b_{n-1}, 2 a_0}] = [a_0, \overline{a_1, a_2, \ldots, a_n}],
\end{align*}
where we have relabelled $a_j := b_j$ for $1 \leq j \leq n-1$ and $a_n := b_0 = 2 a_0$. Thus $\sqrt{d}$ has eventually periodic continued fraction expansion with preperiod exactly one and period $n$.
[/step]
[step:Express $\sqrt{d}$ via a Möbius identity in the $(n+1)$-st complete quotient]
From the expansion $\sqrt{d} = [a_0, a_1, \ldots, a_n, \theta_{n+1}]$ where $\theta_{n+1}$ is the $(n+1)$-st complete quotient, the [Matrix Identity for Convergents](/theorems/1759) gives
\begin{align*}
\sqrt{d} &= \frac{p_n \theta_{n+1} + p_{n-1}}{q_n \theta_{n+1} + q_{n-1}}.
\end{align*}
By the periodicity established in Step 1, $\theta_{n+1} = \theta_1$ (the expansion repeats with period $n$ after the preperiodic $a_0$), and $\theta_1$ is the first complete quotient:
\begin{align*}
\theta_1 &= \frac{1}{\sqrt{d} - a_0}.
\end{align*}
Rationalising the denominator, $\theta_1 = (\sqrt{d} + a_0)/(d - a_0^2)$; but it is more useful to keep $\theta_1 = 1/(\sqrt{d} - a_0)$ in this form.
[/step]
[step:Expand the Möbius identity and separate rational and irrational parts]
Substitute $\theta_{n+1} = 1/(\sqrt{d} - a_0)$ into the Möbius identity:
\begin{align*}
\sqrt{d} &= \frac{p_n \cdot \frac{1}{\sqrt{d} - a_0} + p_{n-1}}{q_n \cdot \frac{1}{\sqrt{d} - a_0} + q_{n-1}} = \frac{p_n + p_{n-1}(\sqrt{d} - a_0)}{q_n + q_{n-1}(\sqrt{d} - a_0)}.
\end{align*}
Cross-multiplying by the denominator,
\begin{align*}
\sqrt{d}\bigl(q_n + q_{n-1}(\sqrt{d} - a_0)\bigr) &= p_n + p_{n-1}(\sqrt{d} - a_0).
\end{align*}
Expanding the left-hand side,
\begin{align*}
q_n \sqrt{d} + q_{n-1} d - a_0 q_{n-1} \sqrt{d} &= p_n - a_0 p_{n-1} + p_{n-1} \sqrt{d}.
\end{align*}
Group rational and $\sqrt{d}$ terms:
\begin{align*}
\bigl(q_{n-1} d\bigr) + \bigl(q_n - a_0 q_{n-1}\bigr) \sqrt{d} &= \bigl(p_n - a_0 p_{n-1}\bigr) + p_{n-1} \sqrt{d}.
\end{align*}
Since $\{1, \sqrt{d}\}$ is $\mathbb{Q}$-linearly independent (as $d$ is not a perfect square), the rational parts and the coefficients of $\sqrt{d}$ must match separately:
\begin{align*}
q_{n-1} d &= p_n - a_0 p_{n-1}, \\
q_n - a_0 q_{n-1} &= p_{n-1}.
\end{align*}
The first equation rearranges to $p_n = a_0 p_{n-1} + d q_{n-1}$. The second rearranges to $p_{n-1} = q_n - a_0 q_{n-1}$.
[/step]
[step:Derive the Pell identity $p_{n-1}^2 - d q_{n-1}^2 = (-1)^n$]
Compute
\begin{align*}
p_{n-1}^2 - d q_{n-1}^2 &= p_{n-1} (q_n - a_0 q_{n-1}) - q_{n-1} (p_n - a_0 p_{n-1}),
\end{align*}
where we used $p_{n-1} = q_n - a_0 q_{n-1}$ in the first term and $d q_{n-1} = p_n - a_0 p_{n-1}$ in the second term (substituted only after multiplying by $q_{n-1}$). Expanding,
\begin{align*}
p_{n-1}^2 - d q_{n-1}^2 &= p_{n-1} q_n - a_0 p_{n-1} q_{n-1} - q_{n-1} p_n + a_0 p_{n-1} q_{n-1} \\
&= p_{n-1} q_n - p_n q_{n-1}.
\end{align*}
By the [Determinant Identity for Convergents](/theorems/1757), $p_n q_{n-1} - p_{n-1} q_n = (-1)^{n-1}$, hence $p_{n-1} q_n - p_n q_{n-1} = (-1)^n$. Therefore
\begin{align*}
p_{n-1}^2 - d q_{n-1}^2 &= (-1)^n.
\end{align*}
[guided]
We have just completed the critical algebraic computation: the quantity $p_{n-1}^2 - d q_{n-1}^2$ — the Pell form evaluated at the convergent $(p_{n-1}, q_{n-1})$ — is exactly $(-1)^n$.
Let us retrace. From Step 3 we had two integer identities:
\begin{align*}
p_n &= a_0 p_{n-1} + d q_{n-1}, \\
p_{n-1} &= q_n - a_0 q_{n-1}.
\end{align*}
The second says $p_{n-1} + a_0 q_{n-1} = q_n$, and the first says $d q_{n-1} = p_n - a_0 p_{n-1}$. Squaring the first and using the second to eliminate $q_n$ would give a direct formula, but the cleaner route is to observe that $p_{n-1}^2 - d q_{n-1}^2$ factors as $p_{n-1} \cdot p_{n-1} - q_{n-1} \cdot (d q_{n-1})$, and we can substitute both $p_{n-1}$ (using the second identity) and $d q_{n-1}$ (using the first identity):
\begin{align*}
p_{n-1}^2 - d q_{n-1}^2 &= p_{n-1} \underbrace{(q_n - a_0 q_{n-1})}_{= p_{n-1}} - q_{n-1} \underbrace{(p_n - a_0 p_{n-1})}_{= d q_{n-1}}.
\end{align*}
Wait — this substitution is not quite consistent. Let me do it more carefully. We substitute into $p_{n-1}^2 - d q_{n-1}^2$ by replacing **one** factor of $p_{n-1}$ in $p_{n-1}^2$ by $q_n - a_0 q_{n-1}$, and the factor $d q_{n-1}$ in $d q_{n-1}^2$ by $p_n - a_0 p_{n-1}$. That is:
\begin{align*}
p_{n-1}^2 - d q_{n-1}^2 &= p_{n-1} \cdot (q_n - a_0 q_{n-1}) - q_{n-1} \cdot (p_n - a_0 p_{n-1}) \\
&= p_{n-1} q_n - a_0 p_{n-1} q_{n-1} - q_{n-1} p_n + a_0 p_{n-1} q_{n-1} \\
&= p_{n-1} q_n - p_n q_{n-1}.
\end{align*}
The cross-terms $-a_0 p_{n-1} q_{n-1}$ and $+a_0 p_{n-1} q_{n-1}$ cancel. What remains is $p_{n-1} q_n - p_n q_{n-1}$, which is the determinant of the matrix $\begin{pmatrix} p_{n-1} & p_n \\ q_{n-1} & q_n \end{pmatrix}$. By the [Determinant Identity for Convergents](/theorems/1757), $p_n q_{n-1} - p_{n-1} q_n = (-1)^{n-1}$, so $p_{n-1} q_n - p_n q_{n-1} = (-1)^n$. This gives the clean conclusion $p_{n-1}^2 - d q_{n-1}^2 = (-1)^n$.
Why is this called the "Pell identity" for $\sqrt{d}$? Because when $n$ is even, the right-hand side is $+1$, and the pair $(p_{n-1}, q_{n-1})$ is a positive integer solution to $x^2 - d y^2 = 1$ — Pell's equation.
[/guided]
[/step]
[step:Produce a solution for any parity by doubling the period]
Case $n$ even: Directly, $(x, y) = (p_{n-1}, q_{n-1})$ satisfies
\begin{align*}
x^2 - d y^2 &= (-1)^n = 1.
\end{align*}
Both $x = p_{n-1} > 0$ and $y = q_{n-1} > 0$ (convergents of $\sqrt{d} > 0$ are positive from index $\geq 0$). This is a positive integer solution.
Case $n$ odd: $(p_{n-1}, q_{n-1})$ solves the negative Pell equation $x^2 - d y^2 = -1$, not the equation we want. We use the doubled period. Since $\sqrt{d} = [a_0, \overline{a_1, \ldots, a_n}]$ is periodic with period $n$, it is also periodic with period $2n$ (concatenating two copies of the period block):
\begin{align*}
\sqrt{d} &= [a_0, \overline{a_1, \ldots, a_n, a_1, \ldots, a_n}].
\end{align*}
Applying Step 4 with period $2n$ in place of $n$,
\begin{align*}
p_{2n-1}^2 - d q_{2n-1}^2 &= (-1)^{2n} = 1.
\end{align*}
Hence $(x, y) = (p_{2n-1}, q_{2n-1})$ is a positive integer solution. (Alternatively, one may compute $(p_{n-1} + q_{n-1} \sqrt{d})^2 = (p_{n-1}^2 + d q_{n-1}^2) + 2 p_{n-1} q_{n-1} \sqrt{d}$ and verify that the pair $(p_{n-1}^2 + d q_{n-1}^2,\, 2 p_{n-1} q_{n-1})$ satisfies $X^2 - d Y^2 = (p_{n-1}^2 - d q_{n-1}^2)^2 = 1$; this is the norm form of $\mathbb{Z}[\sqrt{d}]$, which will be used below.)
[/step]
[step:Generate infinitely many solutions from the fundamental one]
Let $(x_1, y_1)$ be the positive integer solution constructed above (i.e., $(p_{n-1}, q_{n-1})$ if $n$ even, $(p_{2n-1}, q_{2n-1})$ if $n$ odd), so $x_1^2 - d y_1^2 = 1$ and $x_1, y_1 > 0$. Define $\alpha := x_1 + y_1 \sqrt{d} \in \mathbb{Z}[\sqrt{d}]$. Its norm is
\begin{align*}
N(\alpha) &:= \alpha \cdot \alpha' = (x_1 + y_1 \sqrt{d})(x_1 - y_1 \sqrt{d}) = x_1^2 - d y_1^2 = 1.
\end{align*}
The norm is multiplicative: for $\alpha, \beta \in \mathbb{Z}[\sqrt{d}]$, $N(\alpha \beta) = N(\alpha) N(\beta)$ (a routine consequence of $\sigma$ being a ring homomorphism). Hence for every $k \geq 1$,
\begin{align*}
N(\alpha^k) &= N(\alpha)^k = 1.
\end{align*}
Write $\alpha^k = x_k + y_k \sqrt{d}$ with $x_k, y_k \in \mathbb{Z}$ (this representation exists and is unique because $\mathbb{Z}[\sqrt{d}]$ is a free $\mathbb{Z}$-module on $\{1, \sqrt{d}\}$). Then
\begin{align*}
x_k^2 - d y_k^2 &= N(\alpha^k) = 1,
\end{align*}
so $(x_k, y_k)$ satisfies Pell's equation for every $k \geq 1$.
Finally, since $\alpha = x_1 + y_1 \sqrt{d} > 1$ (because $x_1 \geq 1$ and $y_1 \geq 1$, so $\alpha \geq 1 + \sqrt{d} > 1$), the sequence $(\alpha^k)_{k \geq 1}$ is strictly increasing and unbounded. Hence the pairs $(x_k, y_k)$ are pairwise distinct, and $y_k \to \infty$. We need $x_k, y_k > 0$: this holds because $\alpha > 0$ and $\alpha > 1$ imply $\alpha^k = x_k + y_k \sqrt{d} > 0$ with growing $y_k > 0$, and then $x_k = \alpha^k - y_k \sqrt{d}$; using $\alpha^k > 1$ and a standard induction (or direct computation from the recursion $x_{k+1} = x_1 x_k + d y_1 y_k$, $y_{k+1} = x_1 y_k + y_1 x_k$), $x_k > 0$ as well.
Therefore $\{(x_k, y_k) : k \geq 1\}$ is an infinite family of positive integer solutions to $x^2 - d y^2 = 1$. This completes the proof.
[/step]
