The only constant polynomial that satisfies the relation is $P\equiv 0$, so take $P$ nonconstant. We can suppose wlog that $P$ is monic( otherwise divide by the dominant coefficient) of degree $k$: $P(x)=x^k+b_{k-1}x^{k-1}+...+b_1x+b_0$. Let $f(n)$ be that the integer such that $P(f(n))=aP(n)$. Note that $P$ is increasing and positive from a certain point on, so from $P(f(n))=aP(n)>P(n)$, we deduce that $f(n)>n$ for $n>M$. Fix a positive integer $n>M$ and let $g(t)=f^t(n) $. We have that $P(g(t))=a^tP(n)$. As $\dfrac{P(x)}{x^k}=1$ when $x\rightarrow \infty$, from $$a^k\cdot \dfrac{P(g(t))}{g(t)^k}= \dfrac{ P(g(t+k))}{g(t+k)^k}\cdot \left ( \dfrac{g(t+k)}{g(t)} \right )^k$$we infer that $\dfrac{g(t+k)}{g(t)}\rightarrow a$ when $t\rightarrow \infty$. Now, we have that on one hand $$P(g(t+k))-P(ag(t))=a^kP(g(t))-P(ag(t))=b_{k-1}(a^k-a^{k-1})g(t)^{k-1}+Q(g(t))$$with $\mathrm{deg}\ Q<k-1$. On the other hand $$P(g(t+k))-P(ag(t))=(g(t+k)-ag(t)) \left ( g(t+k)^{k-1}+g(t+k)^{k-2}ag(t)+...+(ag(t))^{k-1}+R(g(t), g(t+k) )\right )$$with $\mathrm{deg}\ R<k-1$. Therefore, $$g(t+k)-ag(t)=\dfrac{b_{k-1}(a^k-a^{k-1})+\frac{Q(g(t))}{g(t)^{k-1}}}{\left ( \frac{g(t+k)}{g(t)}\right )^{k-1}+\left ( \frac{g(t+k)}{g(t)}\right )^{k-2}a+...+a^{k-1}+\frac{R(g(t), g(t+k) )}{g(t)^{k-1}}}$$so when $t\rightarrow \infty$, we get that $g(t+k)-ag(t)\rightarrow \dfrac{b_{k-1}(a-1)}{k}$. However $(g(t+k)-ag(t))_{t\ge 1}$ is a sequence of integers, so from the convergence we infer that it is eventually constant. Hence, for $t>M_2$, we have that $g(t+k)=ag(t)+b$ with fixed $b$, or equivalently $g(t+k)+\dfrac{b}{a-1}=a\left (g(t)+\dfrac{b}{a-1} \right )$. Noting $h(t)=g(t)+\dfrac{b}{a-1}$, we have that $h(t+k)=ah(t)$. Take $P(g(t))=P \left ( g(t)+\dfrac{b}{a-1}-\dfrac{b}{a-1} \right )=P \left ( h(t) -\dfrac{b}{a-1} \right )=S(h(t))$, where $S$ is a polynomial having the same degree as $P$. We have that $S(ah(t))=S(h(t+k))=P(g(t+k))=a^kP(g(t))=a^kS(h(t))$ for all $t>M_2$, hence $S(ax)=a^kS(x)$. If $S$ has a root $\omega$ different from zero, then $a\omega$ is also a root and $|a\omega|>|\omega|$. Iterating, we get that $S$ has an infinite number of roots, whence $S\equiv 0$, or $P\equiv 0$, contradiction. So $S$ has only $0$ as root and therefore $$P(g(t))=S(h(t))=h(t)^k=\left ( g(t)+\dfrac{b}{a-1} \right )^k$$for all $t>M_2$, which implies that $P(x)=(x+c)^k$. From $a(n+c)^k=aP(n)=P(f(n))=(f(n)+c)^k$ for $n$ big enough, we infer that $f(n)+c=\sqrt[k]{a}n+\sqrt[k]{a}c$. Writing the same relation for $n+1$ and substracting, we get that $\sqrt[k]{a}\in \mathbb{Z}$, hence $a=d^k$. This implies that $c\in \mathbb{Q}$; in turn, for any $c\in \mathbb{Q}$ we can find an $a$: write $c=\dfrac{p}{q}$ with $p,q\in \mathbb{Z},\ q>0$ and take $a=(q+1)^k$ to get $f(n)=(q+1)n+p$. Therefore, the answer is $P(x)=\alpha (x+q)^k$ with $\alpha\in \mathbb{R},\ q\in \mathbb{Q}$.

Clearly the zero polynomial is the only constant solution of the problem, so assume that $P$ is a solution with $P$ nonconstant. Dividing $P$ by its leading coefficient, we may assume that $P$ is monic, say of degree $d.$ By assumption there is a sequence of integers $z_n$ such that $aP(n) = P(z_n)$ for all $n\geq 1.$ Let $A =\sqrt[d] a$ and write $P(x)=x^d+bx^{d-1}+Q(x)$ with $\operatorname {deg}Q \leq d-2$. Clearly $\lim_{n \rightarrow \infty} \vert z_n\vert=\infty$ Write the equation in the form : $$ a\cdot \frac {P(n)}{n^d} = \frac {P(z_n)}{z_n^d} \cdot \left ( \frac {z_n}{n} \right )^d$$ Taking the limit at ${n \rightarrow \infty}$ we have $\lim_{n \rightarrow \infty}\frac {\vert z_n \vert}{n} \rightarrow A$. Consider the sequence $x_n \overset {\text {def}}{:=} \frac {z_n}{An}$. Note that $x_n$ is bounded. Writing the original equation in terms of $x_n$ and taking the limit when ${n \rightarrow \infty}$ (calculations omitted) gives us the following equation : $$\lim_{n \rightarrow \infty} An(x_n^d-1)+b(x_n^{d-1}-A) =0$$ Call the above equation $\spadesuit$. Consider the following cases based on the parity of $d$. Case 1 : $d$ is odd. Note that this implies that $\lim_{n \rightarrow \infty}x_n \rightarrow 1$. Using this in the equation obtained in $\spadesuit$, we get $:=$ $$\lim_{n \rightarrow \infty} An(x_n-1) \rightarrow \frac {b(A-1)}{d} \implies \lim_{n \rightarrow \infty} z_n-An \rightarrow \frac {b(A-1)}{d} \overset {\text{def}}{:=} B$$ This means that for all $n$ large enough; we have $z_n=An+B$ (because $z_n$ is a integer). It's now easy to finish off the problem. The original equation rearranges to give $A^dP(n)=P(An+B)$. Consider $c= \frac {B}{1-A}$ and consider $Q(X)$ such that $Q(X)= P(X+c)$. Note that we have $:=$ $$Q(Ax)= P(Ax+c)= P(A(x+c)+B)= A^dP(x+c)=A^dQ(x)$$ From here it easily follows that $Q(x)=x^d$ (say by comparing coefficients ) and hence we get that $\boxed {P(x) \equiv C(x-c)^d}$ is a solution in this case. ‐ Case 2: $d$ is even. Assume that $d=2k$. Note that we cannot deduce $x_n\rightarrow 1$ for $n \rightarrow \infty$. (However $x_n = O(\frac 1n)$.) Note that this time $\spadesuit$ gives us : $$\lim_{n \rightarrow \infty} Ank(x_n^2-1)+ b(x_n-A) = 0$$ Set $y_n= z_n + \frac {b}{2k}$. Then we have that the above rearranges to $\lim_{n \rightarrow \infty} \frac {y_n^2}{An}-An \rightarrow \frac {bA}{k}$. Note that we also have $\lim_{n \rightarrow \infty} \left \vert \frac {y_n}{An} \right \vert \rightarrow 1$. This means we have $\lim_{n \rightarrow \infty} \left( \vert y_n \vert -An\right) \rightarrow \frac {bA}{2k}$ Now consider some distinct integers $i,j$ with $\vert i - j\vert \leq 2$ and $x,y$ large. This means that $\vert y_i - y_j \vert = \pm (z_i-z_j).$However note that $\vert y_i\vert -Ai - \vert y_j\vert -Aj$ goes to $0$ . Hence since the $y_i$'s are integers and $\vert i - j\vert$ takes only finitely many values, hence $2A$ is a integer. Now we are ready to finish. Consider a increasing sequence $k_n$ and $\varepsilon \in \{1,-1\}$ such that $\varepsilon z_{k_n} >0$ for all large $n$. Note that $\varepsilon z_{k_n}-Ak_n \rightarrow \frac {b(A-\varepsilon)}{2k}$. Since $2A$ is a integer, we deduce that $$\varepsilon z_{k_n}-Ak_n=\frac {b(A-\varepsilon)}{2k} \ \overset {\text{def}}{:=} C$$holds for all $n$ large enough. This way the original equation translates to : $$ A^d P(k_n)= (\varepsilon A)^d = P(\varepsilon Ak_n + \varepsilon C)$$ Proceeding in a similar way as the previous case , we conclude that $\boxed {P=c(x-x_0)^d}$ for some $c\in \mathbb R ,x_0 \in \mathbb Q$. Having exhausted all cases, we are done. $\blacksquare$