Does there exist polynomials $P$ s.t. $P(k)\mid k!$ holds for only finitely many $k\in\Bbb N$? 
Question: Does there exist polynomials $P$ with integer coefficients such that  $$P(k)\mid k!\tag{*}$$ holds for only finitely many $k\in\Bbb N$?

I think such polynomial doesn't exist (don't know how to prove it). For some specific polynomials, we can prove there exists infinitely many $k\in\Bbb N$ satisfying $(*)$ by constructing $k=\text{a polynomial}$:


*

*Example: $P(k)=k^n-1$, $n\in \Bbb N$
Choose distinct primes $p_1,p_2,\dots,p_j$ coprime to $n$ such that
$$\ell=\prod_{i =1}^j \frac{p_i-1}{p_i}<\frac 1n.$$
Let $m=p_1p_2\cdots p_j$ and $k=r^{m}$, $r\in\Bbb N$, then
$$P(k)=P(r^m)=r^ {mn}-1=\prod_{d\mid mn 
}\Phi_d(r),$$
where $\Phi_d(r)$ are cyclotomic polynomials. Note that
$$\deg \Phi_d=\phi(d)\le \phi(mn)\le m\phi(n)\le mn \ell<m,$$
Therefore for integer $r>N$, we have $r^m>\Phi_d(r)$ and thus
$$\left(\prod \Phi_d(r)\right)\mid\left(r^m\right)!\quad \iff\quad k^n-1\mid k!.$$

With the aid of Mathematica, "surprisingly" I found that the growth rate of the number of integers smaller than $n$ satisfying $(*)$ approximately equals $C\cdot n$ for some constants $C$:
\begin{align*}
P(k)&=k^2-1,\quad \#(k\le n \text{ satisfying } (*))\sim 0.85 n,\\
P(k)&=k^2+1,\quad \#(k\le n \text{ satisfying }(*))\sim 0.23 n.
\end{align*}
Here are some graphs for different $P$:
$x$-axis: $n$
$y$-axis: $\#(k\le n\text{ satisfying }(*))$

Graph 1: $P(k)=k^2-1$

Graph 2: $P(k)=k^6-6k^4+9k^2-1$

Conjecture: For any polynomial $P$ with integer coefficients, there are infinitely many $k\in \Bbb N$ satisfying $(*)$, and the growth rate 
  $$\#(k\le n\text{ satisfying }(*))\sim C\cdot n$$
  for large $n$ and a constant $C$.

 A: We will be looking at different values of $\mathrm{deg}(P)$ where $P \in \mathbb{Z}(x)$. For $\mathrm{deg}(P)=0$, the question is trivial as $\forall$ $x \geqslant c$ where $P(x)=c$ , we will have $P(x) \mid x!$

Claim: There are infinitely many $x$ such that $P(x) \mid x!$ where $P \in \mathbb{Z}[x]$ and $\deg(P)=1$. 
Let $P(x)=ax+b$ for $a,b \in \mathbb{Z}$. We have: 
$$P(bx)=a(bx)+b = b(ax+1)$$ 
It suffices to show that $P(bx) \mid (bx)!$ infinitely often. 
$$P(bx) \mid (bx)! \iff b(ax+1) \mid (bx)! \iff (ax+1) \mid (bx-1)!$$ 
since $b \mid bx$ and $\gcd(ax+1,x)=1$. Now, since $a$ is finite in magnitude, we choose a prime $p$ such that $p \nmid a$. Next, we set: 
$$x=\frac{p^{k\phi(a)}-1}{a}$$ 
Clearly, $x$ is integral $\forall$ $k \in \mathbb{N}$ since $a \mid p^{\phi(a)}-1 \mid p^{k\phi(a)}-1$. Now, we have: 
$$ax+1=p^{k\phi(a)} \implies \nu_p(ax+1)=k\phi(a)$$
$$\nu_p((bx-1)!)=\sum_{i=1}^\infty \bigg\lfloor\frac{bx-1}{p^i} \bigg\rfloor>\frac{bx-1}{p}-1=\frac{b(p^{k\phi(a)}-1)}{a}-\frac{p+1}{p}$$ 
We see that the power of $p$ that divides $ax+1$ is a polynomial in $k$ while the power dividing $(bx-1)!$ is an exponential. Thus, $\exists$ sufficiently large $M$ s.t. $\forall$ $k \geqslant M$, we have $(ax+1) \mid (bx-1)!$. 
QED

Now, we call a polynomial smooth if there are infinitely many $x$ s.t. all the primes dividing $P(x)$ are less than $x^\epsilon$ $\forall$ $\epsilon>0$. We can see:
$(i)$ If $P(x) \mid x!$ for infinitely many $x \in \mathbb{N}$, then $P$ is specifically smooth for $\epsilon = 1$
Proof : For prime $p$, we have $p \nmid x! \iff p>x$. Thus, we will have $p \nmid P(x)$ infinitely often for any $p>x$ which shows specific smoothness at $\epsilon=1$.
$(ii)$ If $P$ is smooth, then $P(x) \mid x!$ for infinitely many $x \in \mathbb{N}$.
Proof : We can find sufficiently large $M$ s.t. $x^M > P(x)$ $\forall$ $x>1$ where $x$ is a positive integer. Since $P$ is smooth, we can have $\epsilon = \frac{1}{k}$ for sufficiently small $k$. This means that any prime $p$ dividing $P(x)$ satisfies $p<x^\frac{1}{k}$ for infinitely many $x$ which will be the values we will consider. Thus:
$$\log_p(x)<\log_p(x^M)=M\log_p(x) \tag{1}$$
$$\nu_p(x!)=\sum_{i=1}^\infty \big\lfloor\frac{x}{p^i} \big\rfloor > \frac{x}{p} \tag{2}$$
We can see that $\frac{x}{p}>M\log_p(x) \iff \frac{x}{M\log_p(x)}>p$. However, we know that $x^\frac{1}{k}>p$ and $\frac{x}{M\log_p(x)}>x^\frac{1}{k}$ will hold infinitely often by comparison of power of $x$ on both sides. Thus, we are done.

Our problem has been solved for $\mathrm{deg}(P)<2$ above. For $\mathrm{deg}(P)=2$, one can refer the paper which is attached below (also present in first comment of OP's question).
Proof for $\mathrm{deg}(P)=2$ : Smoothness of Quadratic Polynomials 
For $\mathrm{deg}(P)>2$, we are to first try to prove the weaker condition of specific smoothness of higher degree polynomials for $\epsilon=1$. However, this seems to be an open problem, only solvable for specific polynomials (such as the example provided by OP) based on the looks of the above paper and personal search and browsing by myself and @TheSimpliFire. Thus, we finally conclude:

The problem has an affirmative answer for $\mathrm{deg}(P) \leqslant
 2$ and is open for $\mathrm{deg}(P) > 2$ due to no answer for specific
   smoothness of integer-coefficient polynomials for $\epsilon=1$.

