Dgdrxrt

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

$$left1,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?

The set in question, call it $A$, is a proper subset of the set, call it $B$, of numbers which can be written as a sum of two squares (the difference between the two sets is that $B$ is the set of all numbers of the form $2^rPQ^2$ where $r$ is a nonnegative integer, the prime factors of $P$ are all $1bmod4$, and the prime factors of $Q$ are all $3bmod4$; $A$ is the same, but with $rle1$ and $Q=1$).

$B$ has density zero (and so, a fortiori, $A$ has density zero). This is discussed at https://mathoverflow.net/questions/205862/sums-of-two-squares-positive-lower-density and at the links to be found there. It goes back to Landau.

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 pmod3$ and in turn $n^2 equiv -1 pmod3$. 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 $pmod72m+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 pmod4$. So you could try and see whether this applies to all numbers in such a class.