Your last statement provides an important clue: not only would $D$ be diagonal, it would have to have $p$ units on the diagonal and zeros elsewhere. So there must be something special about $X(X^\prime X)^{-1}X^\prime$. To see what, look at the Singular Value Decomposition of $X$,
$$X = U \Sigma V^\prime$$
where $U$ and $V$ are orthogonal (that is, $U^\prime U$ and $V^\prime V$ are identity matrices) and $\Sigma$ is diagonal. Use this to simplify
$$(X^\prime X)^{-1} = \left( \left( U \Sigma V^\prime \right)^\prime \left( U \Sigma V^\prime \right)\right)^{-1} = \left( V \Sigma^\prime U^\prime U \Sigma V^\prime \right)^{-1} = \left( V \Sigma^\prime \Sigma V^\prime \right)^{-1}$$
and employ that to compute
$$X(X^\prime X)^{-1}X^\prime =\left( U \Sigma V^\prime \right) \left( V \Sigma^\prime \Sigma V^\prime \right)^{-1}\left( U \Sigma V^\prime \right)^\prime = U\left(\Sigma\left(\Sigma^\prime \Sigma\right)^{-1}\Sigma^\prime\right) U^\prime.$$
This exhibits the covariance matrix of $\epsilon$ as being conjugate (via the similarity induced by $U$) to $\Sigma\left(\Sigma^\prime \Sigma\right)\Sigma^\prime$, which (since $\Sigma$ is diagonal) has $\text{rank}(\Sigma)$ ones along the diagonal and zeros elsewhere: in other words, the distribution of $\epsilon$ is that of an orthogonal linear combination of $\text{rank}(\Sigma) = \text{rank}(X)$ independent, identically distributed Normal variates.
Orthogonal transformations such as $U$ preserve sums of squares. Provided $X$ has full rank--which is $\min(p,n)=p$--the distribution of $E$ therefore is that of the sum of $p$ squares of independent standard Normal variables, which by definition is $\chi^2(p)$. More generally, $E \sim \chi^2(\text{rank}(X)).$
This algebraic argument is one way of finding out that ordinary least squares is just Euclidean geometry: this result is a rediscovery of the Pythagorean Theorem.
I found by problem by generating random matrices and checking every equality.
The problem is in the second last equality. Turns out,
$$
L'^{-1}QR(R'R)^{-1}R'Q'L^{-1}\neq L'^{-1}L^{-1}
$$
because
$$
R(R'R)^{-1}R'
=\begin{bmatrix}
I_p&0\\
0&0
\end{bmatrix}
\neq I
$$
where $I_p$ is the identity matrix of dimension $p<N$, where $N$ is the dimension of the full matrix.
Best Answer
In classical statistics the parameter value $\beta$ in a linear regression model is an unknown constant. The value $\hat{\beta}$ is not a parameter - it is an estimator of the parameter, which is a function of the data. The reason this estimator is normally distributed is that it is a linear function of the underlying error vector (as written in the equation you have shown), which is normally distributed under the model assumptions. (Note that even if you relax of the normality assumption, the parameter estimator will still be a summation quantity for which you can invoke the CLT under fairly general conditions; so the distribution of the parameter estimator will converge to the normal distribution under broad conditions even if the model is misspecified.)