[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*}[/step]

[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: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]