Essentially just this.

Consider pell equatios: $x^2-2y^2=-1$, let $x_n=\frac{\alpha^{2n+1}+\beta^{2n+1}}{2}$, $y_n=\frac{\alpha^{2n+1}-\beta^{2n+1}}{\alpha-\beta}$, where, $\alpha=1+\sqrt{2}$, $\beta=1-\sqrt{2}$, $\alpha \beta=-1$. We claim: if $2k+1\mid 2n+1$, then $y_k \mid y_n$. Indeed, notice: $$(\frac{\alpha^{2k+1}-\beta^{2k+1}}{\alpha-\beta})(\alpha^{2n-2k}+\beta^{2n-2k})=\frac{\alpha^{2n+1}-\beta^{2n+1}+(\alpha^{2n+1-2(2k+1)}-\beta^{2n+1-2(2k+1)}}{\alpha-\beta}).$$It follows that $$gcd(\frac{\alpha^{2k+1}-\beta^{2k+1}}{\alpha-\beta}, \frac{\alpha^{2n+1}-\beta^{2n+1}}{\alpha-\beta})=gcd(\frac{\alpha^{2k+1}-\beta^{2k+1}}{\alpha-\beta},\frac{\alpha^{2n+1-2(2k+1)}-\beta^{2n+1-2(2k+1)}}{\alpha-\beta})=\cdots=\frac{\alpha^{2k+1}-\beta^{2k+1}}{\alpha-\beta}.$$So we have $y_n \mid y_{3n+1}$, and $y_n=O(x_{3n+1}^{\frac{1}{3}})< x_{3n+1} \cdot \pi^{-2019}$, $n\rightarrow +\infty$, and $\frac{y_{3n+1}}{x_n}=O(x_{3n+1}^{\frac{2}{3}})<x_{3n+1} \cdot \pi^{-2019}$, $n\rightarrow +\infty$. We did it.

Let us prove that for every $\varepsilon > 0$ there exist infinitely many integers $n > 0$ such that $P(n^2 + 1) < \varepsilon n$. Obviously, it is enough to prove that there exists one such $n$. As you noticed, if $n = 2m^2$, for some integer $m > 0$, then $n^2 + 1 = (2m^2 - 2m + 1)(2m^2 + 2m + 1)$. Hence, it is enough to prove that we can find $m$ such that $P(2m^2 - 2m + 1) < 2\varepsilon m^2$ and $P(2m^2 + 2m + 1) < 2\varepsilon m^2$. It is well known and not difficult to prove that for every polynomial $f \in \mathbb{Z}[x]$ the set of prime numbers $p$ such that $p$ divides $f(k)$ for some integer $k > 0$ is infinite. Therefore, we can find an integer $\ell > 0$ such that $p_1 := P(2\ell^2 - 2\ell + 1) > 2/\varepsilon$ and $p_2 := P(2\ell^2 + 2\ell + 1) > 2/\varepsilon$. Now letting $m := \ell + t p_1 p_2$, for some integer $t > 0$, we get that $p_1$ divides $2m^2 - 2m + 1$ and $p_2$ divides $2m^2 + 2m + 1$. Hence, $P(2m^2 - 2m + 1) < \max(p_1, (2m^2 - 2m + 1) / p_1)$ and $P(2m^2 + 2m + 1) < \max(p_2, (2m^2 + 2m + 1) / p_2)$. Taking $t$ sufficiently large, we get that $P(2m^2 + 2m + 1) < (2m^2 + 2m + 1) / p_2 < \varepsilon(2m^2 + 2m + 1) / 2 < \varepsilon 4m^2 / 2 < 2\varepsilon m^2$ and, similarly, $P(2m^2 - 2m + 1) < 2\varepsilon m^2$, as required. (given c>0 and nonconstant polynomials f1,f2,…,fk ∈ Z[x], we can prove that there exists an integer n>0 such that each of f1(n),…,fk(n) has a prime factor greater than c. This can be done by induction, since once we know that f1(n) is divisible by a prime p>c, then we have that f1(n+tp) is divisible by p, so we have to chose t such that f2(n+tp) has a prime factor greater than c... and so on )