IMC 2025 · Day 2 · P10

Throughout the solution, we use the Vinogradov notation $A \ll B$ to mean $A = O(B)$, which in turn means that there exists a constant $C > 0$, independent of the quantities $A$ and $B$, such that $|A| \le C |B|$, on the entirety of the domain where $A$ and $B$ are defined (for us, this will always be the interval $[1, \infty)$.)

We will show that the limit equals 1, corresponding to the trivial solutions $a = b$. Note that $(a^2 + a)(b^2 + b)$ is a perfect square if and only if $a^2 + a = d z_1^2$ and $b^2 + b = d z_2^2$ for some square-free $d$ and $z_1, z_2 \in \mathbb{Z}_{>0}$. From this point on, all sums over $d$ will be over square-free positive integers. Multiplying the equations by 4 and setting $y_i = 2 z_i$, we get $$S_N = \sum_{d \ll N^2} c_d(N)^2 + O(1),$$ where $c_d(N)$ is the number of solutions to $(2k + 1)^2 - d y^2 = 1$ with $1 \le k \le N$ and $1 \le y \le N/2$ with $y$ even. Other than for the purpose of identifying the trivial solutions, we will work with Pell's equation $x^2 - d y^2 = 1$ with $1 \le x, y \ll N$. Split the sum as $$\sum_{\substack{d \ll N^2 \\ c_d(N) \le 1}} c_d(N) + \sum_{\substack{d \ll N^2 \\ c_d(N) > 1}} c_d(N)^2.$$ Note that if $d \gg N$, then the size of the second solution $x_2 = x_1^2 + d y_1^2$ (coming from $x_2 + y_2 \sqrt{d} = (x_1 + y_1 \sqrt{d})^2$, where $x_1 + y_1 \sqrt{d}$ is the fundamental solution) is $\gg d \gg N$. Hence we may assume that $d \ll N$ if $c_d(N) > 1$ (with a suitable choice of hidden constants). Denote the second sum by $E$ (for error, which we will bound momentarily). The first sum is easily manipulated into being asymptotic to $N$ (up to the error $E$), using the fact that fixing $x = 2a + 1 \le 2N + 1$ fixes the square-free $d$ and the square $y^2$, namely $$\sum_{\substack{d \ll N^2 \\ c_d(N) \le 1}} c_d(N) = \sum_{d \ll N^2} \sum_{\substack{1 \le a \le N \\ 1 \le y \le N/2}} \chi_{(2a+1)^2 - d y^2 = 1} + O(E) = \sum_{a=1}^{N} \sum_{d, y} \chi_{(2a+1)^2 - 1 = d y^2} + O(E) = N + O(E).$$ Here $\chi_\cdot$ denotes the characteristic function (taking the value 0 if $\cdot$ is not satisfied, and 1 otherwise).

Now we bound the error sum $E$. Note that solutions to Pell's equation $x^2 - d y^2 = 1$ grow exponentially, hence we have $c_d(N) \ll \log N$. This means we may assume that $N \gg d \gg N^{1-\delta}$ for some small enough fixed $\delta > 0$, since the contribution of $d \ll N^{1-\delta}$ is bounded by $N^{1-\delta} \log N$. By $x^2 - d y^2 = 1$, we have that $d \gg N^{1-\delta}$ implies $y \ll N^{1/2 + \delta}$.

Fixing each $y \ll N^{\delta/2}$ gives $\ll N^{1-\delta}$ choices for $x$ (hence also for $d$). So we may assume $N^{1/2+\delta} \gg y \gg N^{\delta/2}$, since the contribution of $y \ll N^{\delta/2}$ is bounded by $N^{1-\delta/2}$.

By placing $x$ in residue classes modulo $y^2$ and splitting the interval $[1, 2N+1]$ into intervals of length $y^2$, we get that each choice of $N^{\delta/2} \ll y \ll N^{1/2+\delta}$ gives $\ll N g(y) / y^2$ choices for $x$ (hence also for $d$) by $y^2 \mid x^2 + 1$, where $g(y) = |\{ 1 \le x \le y^2 : x^2 + 1 \equiv 0 \pmod{y^2} \}|$.

By elementary number theory, $g$ is multiplicative and $g(p^k) \le 2$ for all prime powers $p^k$. In particular we obtain $g(n) \le \tau(n) \ll n^\varepsilon$ for any $\varepsilon > 0$ (this is not hard to prove directly for $g$, but may be used as a well-known fact for the divisor function $\tau$). Therefore the contribution of such $y$ is $$\ll \sum_{N^{\delta/2} \ll y \ll N^{1/2+\delta}} \frac{N g(y)}{y^2} \ll N^{1 - \delta/2 + \varepsilon},$$ which is acceptable by choosing $\varepsilon > 0$ small enough. We conclude that $S_N = N (1 + o(1))$, as desired.

Remark. There is a secondary infinite family of solutions of “size” $\sqrt{N}$, namely given by $a = 4b(b+1)$. This shows that $$\limsup_{N \to \infty} \frac{S_N - N}{\sqrt{N}} > 0.$$