Let $f(n)$ be the largest prime factor of $n$. The image of function $g(n)=\sqrt{f(n^2+1)}$ is like this:
Question: If we want to draw a horizontal line which bisects the points from $n=1$ to $n=x,$ this line is probably in what position when $x\to \infty$?
The first thing we can find on the image is that there are many approximate straight lines.
Then what's the slope of these straight lines?
I find that the slope of the most inclined straight line is $1,$ hence the equation is $y=x.$ The second line is $y=\dfrac{x}{\sqrt2},$ the third is $y=\dfrac{x}{\sqrt5},$
then $y=\dfrac{x}{\sqrt{10}},$ then $y=\dfrac{x}{\sqrt{13}}.$
$\dfrac{1}{\sqrt{m}}$ can be a slope if and only if $a^2+1\equiv 0\pmod m$ has integer solution.
In fact, the point $(n,\sqrt{f(n^2+1)})$ on $y=x$ means $g(n)\approx n$ hence $n^2+1$ is prime. $y=\dfrac{x}{\sqrt2}$ means $g(n)\approx \dfrac{n}{\sqrt2}$ hence $2\mid n^2+1$ and $\dfrac{n^2+1}{2}$ is prime $\dots$
Whether there are infinity many primes of the form $n^2+1$ is still an open problem. So we can not prove these conclusions, but we can estimate it.
Best Answer
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:
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.