A problem about the largest prime factor of $n^2+1$

What you are trying to calculate is the median largest prime factor of $f(n)=n^2+1$, where $1\leq n\leq x$. This will be a function of $x$, which we denote by $M_{n^2+1}(x)$, and this will grow to infinity as $x$ grows to infinity.

To get a handle on things, lets first take a step back, and ask an easier question:

What is the median largest prime factor of the integers $n$ in the range $1\leq n\leq x$?

Let us denote this function by $M(x)$. This question was raised on Math Overflow, and I proved that we have the asymptotic $$M(x)\sim e^{\frac{\gamma-1}{\sqrt{e}}}x^{\frac{1}{\sqrt{e}}},\ \ \ \ \ \ \ \ \ \ (1)$$ where $\gamma$ is the Euler Mascheroni constant. See the corresponding paper on the arXiv for details.

So what about $n^2+1$? This seems like it must be significantly more difficult, as now we are dealing with the prime factors of a polynomial, which usually increases the complexity of a problem significantly. Modifying section $4$ of the paper I mentioned above, we can relate the problem to $$\psi_{f}(x,y)=\left|\{n\leq x:\ P(f(n))\leq y\}\right|$$ when $f(n)=n^2+1$, where $P(n)$ denotes the largest prime factor of $n$. Unfortunately, not too much is known about $\psi_f(x,y)$ when $f$ is a polynomial of degree $2$ or higher. The specific case $f(n)=n^2+1$ has been looked at by Dartyge, and Harman, where they give lower bounds on the order of magnitude for $u=\frac{\log x}{\log y}$ in a fairly small range. To obtain an asymptotic however, we would need something stronger - a precise asymptotic for $\psi_{n^2+1}(x,y)$, where $f(n)=n^2+1$, with error of size $\frac{x}{\log^2 x}$. (That is, we need to know the $\frac{x}{\log x}$ term as well) Under a major assumption, Martin proves that we do in fact have $$\psi_F(x,y)\sim x\rho\left(\frac{\log F(x)}{\log y}\right)$$ when $F$ is an irreducible polynomial, and where $\rho$ is the Dickman de Bruijn rho function. ($u=\frac{\log x}{\log y}$ will lie in a particular bounded range) Using Martin's result, and the methods I mentioned above, we are able to obtain that your quantity is $$M_{n^2+1}(x)\approx x^{2/\sqrt{e}}.$$ Note that this requires assuming a version of Hardy and Littlewood's second conjecture.

Personally, I am not entirely satisfied without an asymptotic. I believe that Martin's work can be modified to show that $$\psi_F(x,y)=x\Lambda\left(F(x),y\right)+O_F\left(\frac{x}{\log^2 x}\right)$$ for $u=\frac{\log x}{\log y}$ in a particular bounded range, and where $\Lambda(x,y)$ is a function appearing in de Bruijn's original paper on the subject. Using the expansion for $\Lambda(x,y)$ appearing in Eric Saias' 1989 paper, we would be able to recover the asymptotic

$$M_{n^2+1}(x) \sim e^{\frac{\gamma-1}{\sqrt{e}}} x^{2/\sqrt{e}}, \ \ \ \ \ \ \ \ \ (2)$$ which is very similar to equation $(1)$ above.

Moreover, I believe that given any irreducible integer polynomial $F$ of degree $d$, we have that the median largest prime factor of $F(n)$ in the interval $[1,x]$ satisfies the asymptotic $$M_F(x)\sim e^{\frac{\gamma-1}{\sqrt{e}}} x^{d/\sqrt{e}}.\ \ \ \ \ \ \ \ \ (3)$$

Under the assumptions in Martin's paper, the above outline should be able to prove that equation $(3)$ holds when the degree is less than or equal to $2$. However, additional issues arise when $d\geq 3$. Specifically, there are problems regarding the range of $u$ which I did not discuss so thoroughly above, and the previous bounds we had on $\psi_F(x,y)$ are not in a large enough range of $u$ for the proof to work.

I hope that answers your question, and that equation $(2)$ is what you are looking for.

Here's some non-optimized code

prim = Table[First[Last[FactorInteger[n^2 + 1]]], {n, 1, 100000}];   
ListPlot[Table[With[{temp = Sort[Take[prim, 2 k]]}, (temp[[k]] + temp[[k + 1]])/2], 
  {k, 2, 20000}]]

That leads to the following picture. I'd say that's too irregular to predict long-term behavior.

EDIT -- I was multiplying by 2, twice, which skewed the values. Here's better code.

prim = Table[First[Last[FactorInteger[n^2 + 1]]], {n, 1, 100000}];  
sortedprimes = {};   
vals = {};   
Monitor[
  Do[   
  addprime = Sort[{prim[[k - 1]], prim[[k]]}];
  sortedprimes = Sort[Join[sortedprimes, addprime]];
  newval = (sortedprimes[[k/2]] + sortedprimes[[k/2 + 1]])/2;
  vals = Append[vals, newval],
  {k, 2, 100000, 2}],
k]
ListPlot[Transpose[{Range[50000] 2, vals}]]

Gives the following picture

We can try fitting the data in various ways. Fitting to $A x + B x/log(x)$ didn't immediately go bad when more terms were added.

ListPlot[{#[[1]], #[[2]] - (25.1 #[[1]] - (185.1 #[[1]])/
   Log[#[[1]]])} & /@ Transpose[{Range[50000] 2, vals}]]

Leading to the picture