Fix a positive integer $a>1$ and let $f\in\mathbb{Z}[x]$ be a polynomial with positive leading coefficient. We define a set $S$ of positive integers,
$$
S=\{n\in\mathbb{Z}^+:n\mid a^{f(n)}-1\}.
$$
Question is, can we compute the density of $S$, that is,
$$
\lim\_{n\to\infty}\frac{|S\cap \{1,2,\dots,n\}|}{n}?
$$
I believe that the answer should be $0$, yet I don't have a proof.



Density of a set of numbers dividing a fixed number with polynomial exponent
Fix $\varepsilon>0$ and choose large prime $q>5\varepsilon^{-1} \cdot \deg(f)$ such that $f$ is non-trivial modulo $q$ (then $f$ has at most $\deg(f)$ roots modulo $q$) and additionally such that $a$ is not a $q$-th perfect power. Say that a prime $p>q$ is $q$-appropriate if $a$ is not $q$-th power modulo $p$. By Chebotarev density theorem the $q$-appropriate primes have positive density. Thus we may choose $q$-appropriate primes $p\_1,\dots,p\_m$ such that $\prod\_{i=1}^m (1-1/p\_i)<\varepsilon/5$. Note that if $n$ divides $a^{f(n)}-1$, then either $q$ divides $f(n)$ (density of such $n$ does not exceed $\deg(f)/q<\varepsilon/5$), or neither $p\_i$ divides $n$. [Indeed, if $p\_i|n|a^{f(n)}-1$, but $q$ does not divide $f(n)$, we may find $s$ such that $sf(n)+1$ is divisible by $q$, that yields $a\equiv a^{1+sf(n)} \pmod {p\_i}$, so $a$ is a $q$-th power modulo $p\_i$, a contradiction.] The density of such $n$'s is exactly $\prod\_{i=1}^m (1-1/p\_i)<\varepsilon/5$. So $S$ is covered by finitely many arithmetic progressions with sum of densities less than $\varepsilon$.


