It is true that, for $n$ large enough, $p_n(x)$ will have distinct real roots, which I will prove by induction on $m$.
Before proceding with the proof, let me first state what it shows about the behaviour of $p_n$ as $n$ increases. Assuming the leading coefficient is positive, $p_n^{(m-2)}$ are quadratics with a unique minimum. As $n$ increases, the minimum decreases monotonically until it becomes negative, at which point $p_n^{(m-2)}$ has distinct real roots. Then, $p_n^{(m-3)}$ is a cubic with a local maximum and a local minimum. As $n$ increases further, the local maximum increases until it is positive and the local minimum decreases until it is negative, at which point $p_n^{(m-3)}$ has distinct real roots. More generally, for each $k=0,1,2,\ldots,m-2$, once $p_n^{(k+1)}$ has distinct real roots then $p_n^{(k)}$ has $m-k-1$ local extrema. Further increasing $n$, the local maxima increase monotonically until they are positive and the local minima decrease until they are negative, at which point $p_n^{(k)}$ has distinct real roots. By induction then, $p_n$ eventually has distinct real roots.
Now, to continue with the proof.
For $m=1$, $p_n(x)$ is of first order so always has a real root.
For $m > 1$, we now suppose that the statement holds for polynomials of degree $m-1$. Then, setting $q_n(x)=p_n^\prime(x)$, we have $q_n(x)=q_{n-1}(x)+q_{n-1}^\prime(x)$. By the induction hypothesis, $q_n(x)$ has $m-1$ distinct real roots for all large enough $n$. I'll denote these by $a_{n,1} > a_{n,2} > \cdots > a_{n,m-1}$. These are the turning points of $p_n(x)$. I'll show that the for large $n$, the signs of $p_n(a_{n,1}),p_n(a_{n,2}),p_n(a_{n,3}),\ldots$ are alternating, from which it follows that $p_n$ has $m$ distinct roots.
By scaling, we assume wlog that $p_n(x)$ has leading coefficient $1$ (i.e., it is monic). Then, $q_n$ is positive on $(a_{n,1},\infty)$, negative on $(a_{n,2},a_{n,1})$, positive on $(a_{n,3},a_{n,2}))$, etc.
Note that $a_{n,i}$ is a local maximum of $(-1)^ip_n(x)$ (for each $i=1,2,\ldots,m-1$). This means that
$$
(-1)^ip_n(a_{n,i}+h)=(-1)^ip_n(a_{n,i})-ch^{2r}+O(h^{2r+1})
$$
for some positive $c$ and positive integer $r$. So,
$$
(-1)^ip_{n+1}^\prime(a_{n,i}+h)=-c(2r)(2r-1)h^{2(r-1)}+O(h^{2r-1}).
$$
This shows that $(-1)^ip_{n+1}(x)$ is decreasing in a neighbourhood of $a_{n,i}$.
Now, for $i=1,2,\ldots,m-2$, we see that $(-1)^ip_{n+1}(x)$ is increasing at $a_{n,i+1}$ and decreasing at $a_{n,i}$. Hence, it has a local maximum in $(a_{n,i+1},a_{n,i})$, which we will denote by $b_{i}$, so $a_{n,i+1} < b_i < a_{n,i}$ and $(-1)^ip_{n+1}(b_i) > (-1)^ip_{n+1}(a_{n-i})$. Similarly, $(-1)^{m-1}p_{n+1}(x)$ is decreasing at $a_{n,m-1}$ and (as $p_{n+1}$ is monic of degree $m$), it tends to $-\infty$ as $x\to-\infty$. Hence, it has a local maximum, $b_{m-1}$ in the range $(-\infty,a_{n,m-1})$. So, $b_{m-1} < a_{n,m-1}$ and $(-1)^ip_{n+1}(b_{m-1}) > (-1)^ip_{n+1}(a_{n,m-1})$. As $b_1 > b_2 > \cdots > b_{m-1}$ are roots of $p_{n+1}^\prime(x)$, we have $a_{n+1,i}=b_i$.
So far, we have shown that
$$
a_{n,1} > a_{n+1,1} > a_{n,2} > a_{n+1,2} > a_{n,3} > \cdots > a_{n,m-1} > a_{n+1,m-1}
$$
and, for each $i$, $(-1)^ip(a_{n,i})$ is strictly increasing in $n$. If, for each fixed $i$, it can be shown that there is an $\epsilon > 0$ with $(-1)^ip_{n+1}(a_{n+1,i})\ge(-1)^ip_n(a_{n,i})+\epsilon$ for all large $n$ then, for $n$ large enough, $(-1)^ip_n(a_{n,i})$ will be positive. So, the signs of $p_n(a_{n,1}),p_n(a_{n,2}),\ldots$ will be alternating in sign and we are done. Let us proceed by contradiction and assume to the contrary that $(-1)^ip_{n+1}(a_{n+1,i})< (-1)^ip_n(a_{n,i})+\epsilon$ (for small enough $\epsilon$ this will give a contradiction regardless of $n$).
If $i < m-1$ then $(-1)^ip_{n+1}(x)\le(-1)^ip_{n+1}(a_{n+1,i})$ on $(a_{n,i+1},a_{n,i})$. Setting $a_{n,m}=-\infty$ then this also holds for $i=m-1$. By the assumption, $(-1)^ip_{n+1}(x)\le(-1)^ip_{n}(a_{n,i})+\epsilon$ in this range so, setting $f(x)=(-1)^i(p_n(a_{n,i})-p_{n}(x))\ge0$,
$$
f(x)+f^\prime(x)=(-1)^i(p_n(a_{n,i})-p_{n+1}(x))\ge-\epsilon.
$$
Therefore,
$$
f(x)=-\int_x^{a_{n,i}}\frac{d}{dy}(e^{y-x}f(y))dy
=-\int_x^{a_{n,i}}e^{y-x}(f(y)+f^\prime(y))dy\le\epsilon\int_x^{a_{n,i}}e^{y-x}dy=\epsilon(e^{a_{n,i}-x}-1).
$$
This shows that in the range $(a_{n,i+1},a_{n,i})$,
\begin{align}
\lvert p_{n}(x)-p_n(a_{n,i})\rvert\le\epsilon(e^{a_{n,i}-x}-1).&&{\rm(1)}
\end{align}
So, $p_n$ is almost constant here for $\epsilon$ small. I'll need the following lemma, which I'll prove in a moment.
Lemma: Let $p$ be a monic polynominal of degree $m\ge1$. Then, for each $a < b$ we have
$$
sup_{x\in(a,b)}\lvert p(x)-p(b)\rvert\ge L(b-a)^m.
$$
where $L$ is a constant depending only on $m$.
We can apply this to (1). For the case $i=m-1$ and any $K>0$, apply to the range $[a_{n,m-1}-K,a_{n,m-1}]$ to get
$$
\epsilon(e^K-1)\ge L K^{m}
$$
which gives a contradiction if $\epsilon$ is small enough.
On the other hand, for $i< m-1$, $(-1)^i(p_n(a_{n,i})-p_n(a_{n,i+1})$ is positive and increasing in $n$ for all large enough $n$, so is greater than some fixed $\delta$. Applying (1),
$$
\delta\le (-1)^i(p_{n,i}(a_{n,i})-p_{n,i}(a_{n,i+1}))\le\epsilon(e^{a_{n,i}-a_{n,i+1}}-1).
$$
Assuming that $\epsilon \le 1$, we choose $K > 0$ such that $\delta\ge e^K-1$. Then, $K\le a_{n,i}-a_{n,i+1}$ and we can apply the lemma to the range $[a_{n,i}-K,a_{n,i}]$ to get
$$
\epsilon(e^K-1)\ge LK^m
$$
again, giving a contradiction for small $\epsilon$. QED
I'll now prove the lemma. For any $\epsilon\in\mathbb{R}$ define the linear operator $T_\epsilon$ of the functions $f\colon\mathbb{R}\to\mathbb{R}$ by $T_\epsilon f(x)=f(x+\epsilon)$. If $f$ is a polynomial of degree $m$ with leading coefficient $c$, then $T_\epsilon f-f$ is a polynomial of degree $m-1$ with leading coefficient $mc\epsilon$. Hence, $(T_\epsilon-1)^mf(x)=m!c\epsilon^m$. On the other hand $\lVert T_\epsilon\rVert=1$ and $(T_\epsilon-1)^mf(x)$ only depends on the values of $f(y)$ with $y$ in the range $[x,x+m\epsilon]$. So, for a monic polynomial $f$ of degree $m$,
$$
m!\epsilon^m=(T_\epsilon-1)^mf(x)\le 2^m\sup_{y\in[x,x+m\epsilon]}\lvert f(y)\rvert
$$
If $p$ is a monic polynomial of degree $m$ and $a < b$, then taking $f(y)=p(y)-p(b)$, $x=a$ and $\epsilon=(b-a)/m$ gives the lemma with $L=m!(2m)^{-m}$.
Best Answer
The characteristic polynomial is crucial to understanding how the function behaves because of how the difference equation works. This answer assumes you are very comfortable with the basics of linear difference equations and that have a good theoretical grasp of how functions work.
We're going to work on getting a closed form for $p_n(x)$. The key realization here is that the value $p_n(1)$, for example, only depends on $p_1(1),p_2(1),p_3(1),\dots,p_{n-1}(1)$. In general, for a fixed $t$, $p_n(t)$ only depends on the values of the other functions at $t$. Therefore, if we fix $t$, we can explore the behavior of the sequence $p_1(t),p_2(t),p_3(t),\dots$ without worrying about how the function behaves in other places.
I'll take your equation, $p_n(x)=(x+2)a_1p_{n-1}(x)+[(d-a_1)x-a_1]p_{n-2}(x)$ as given. If we fix $x$, then we can write $p_n(x)-(x+2)a_1p_{n-1}(x)-[(d-a_1)x-a_1]p_{n-2}(x)=0$, which is a linear difference equation with constant coefficients (remember $x$ is fixed). This means that we can solve it the standard way; assume that $p_n=\lambda^n$ is a solution for some constant $\lambda$ (or, in the world of the function $p_n(x)$, $\lambda$ depends on $x$). This gives us the characteristic function $\lambda^2-(x+2)a_1\lambda-((d-a_1)x-a_1)=0$. This will, of course, give rise to two (either both real or complex conjugate) solutions $\lambda_{1,2}(x)$, and then for suitable constants $c_1,c_2$ we have that $p_n(x)=c_1(\lambda_1(x))^n+c_2(\lambda_2(x))^n$. (This assumes $\lambda_1\neq\lambda_2$: the case where the two are equal is unimportant to the overall explanation and adds extra complexity, so I won't be covering it.)
So far, we've used the same ideas that exist in any beginner's course on difference equations. The only level of abstraction to get your head around is that we are using functions of $x$ instead of sequences. This is important to understand in order to answer the question you posed.
Let's go back to that characteristic polynomial, $\lambda^2-(x+2)a_1\lambda-((d-a_1)x-a_1)=0$. Clearly, a choice of $x,a_1,d$ fixes constant values of $\lambda_{1,2}$. The discriminant of that characteristic polynomial is, as you said, $\Delta=a_1^2x^2+4[a_1(a_1-1)+d]x$, a function of $x$. But what does the discriminant of a quadratic mean? You'll remember that if the discriminant is positive, then there are two distinct real solutions for the quadratic. If the discriminant is negative, then there are two non-real, complex conjugate solutions of the quadratic.
So, if we are seeking to evaluate $p_n(x_1)$, for some fixed real $x_1$, then that choice of $x_1$ will make $\Delta>0$ or $\Delta<0$. (Again, $\Delta=0$ gives the double root, which you should explore in your own time.) If we have $\Delta(x_1)>0$, then the corresponding values $\lambda_{1,2}(x_1)$ will be real and distinct. In that case, we may be able to find a root to $p_n(x_1)=c_1\lambda_1(x_1)^n+c_2\lambda_2(x_1)^n$. But if $\Delta(x)<0$, then the corresponding values $\lambda_{1,2}(x_1)$ will be nonreal and complex conjugates. Let's explore this case in more detail:
Again, we've already chosen $x_1$, so we'll just write $p_n,\lambda_1,\lambda_2$ to save space and keep me sane. Now, we know that $\lambda_{1,2}$ are complex conjugates, and so $\lambda_1^n$ and $\lambda_2^n$ are also complex conjugates. So we can write: $\lambda_1^n=r+bi$ and $\lambda_2^n=r-bi$. So if want $p_n=c_1\lambda_1+c_2\lambda_2=0$, then $r(c_1+c_2)+b(c_1-c_2)i=0$. So either $\lambda_1=\lambda_2=0$ (which is false), or $c_1=c_2=0$ (which is false), or $\lambda_{1,2}$ are pure imaginary and $c_1=c_2$ (which turns out to never be true). So clearly, we can't have $p_n(x_1)=0$.
Let's summarize: when we choose $x_1$ so that $\Delta(x_1)<0$, then we must have complex conjugate solutions $\lambda_{1,2}$ and therefore $p_n(x_1)\neq0$. But if $\Delta(x_1)>0$, then $\lambda_{1,2}$ are real and we might have a root $p_n(x_1)=0$. So, every root $x_0$ of $p_n$ satisfies $\Delta(x_0)>0$. (Technically $\geq 0$, but we're ignoring $\Delta=0$.) But when is $\Delta(x_0)>0$? Why, whenever $x_0$ is in between the two roots $x_1,x_2$ of $\Delta$! So we know that if $p_n(x_0)=0$, then:
$$-4-\frac{4(d-a_1)}{a_1^2}<x_0<0$$
And then we add in the trivial solution $x_0=0$ to get the final inequality.
I hope this helped! Please let me know if you have any questions :)