We want to find the numbers that divide natural numbers in the form of $n^2+1$ and solve for their natural density.
Using Wolfram Mathematica, I found divisors from $n=0$ to $1000000$ and eliminated repeated divisors. Here is the list
$$\left\{1,2,5,10,13,17,25,26,29,34,37,41,50,53,58,61,65,73,74,82,85,89,97,101,106,109,113,122,125,130,137,145,146,149,157,169,170,173,178,181,185,193,194,197,202,205,218,221,226,229,233,241,250,257,265,269,274,277,281….. \right\}$$
After looking at the list extensively my guess is the density is zero.
Is there a mathematical way of finding this without computer programming?
Best Answer
As I understand it, you claim that every divisor of a number of the form $n^2+1$ belongs to one of some set of residue classes (mod 72).
This appears to be true. For instance if any multiple of 3 were to divide $n^2+1$ then we would have $n^2+1 \equiv 0 \pmod{3}$ and in turn $n^2 \equiv -1 \pmod{3}$. But -1 is not a quadratic residue (mod 3), so we can rule out classes like 72m, 72m+3, etc. as potential divisors of $n^2+1$.
I didn't verify your list, but clearly some such list exists.
If you wanted to go from here to calculating the asymptotic density, you'd need to prove something stronger, say that every (perhaps sufficiently large) number belonging to such a class divides some number of the form $n^2+1$.
UPDATE
In fact let's consider that question. Say you have some number $72m + r$ and claim it divides a number of the form $n^2 + 1$. That is, for some k we have $(72m + r)k = n^2 + 1$. Then $-1 + (72m+r)k$ is a perfect square for some $k>0$. But this is the same as $-1$ being a quadratic residue $\pmod{72m+r}$.
My elementary number theory is rusty, but I think that's equivalent to every odd prime divisor of $72m+r$ being of the form $1 \pmod{4}$. So you could try and see whether this applies to all numbers in such a class.