Ultimately my goal with this post is to find the orders of magnitude of $E[(\beta_n-\beta)^k]$ for $k\ge 3$ as $n \to \infty$ so that I can be sure that is safe to drop higher order terms in a Taylor expansion based approximation. As a first step towards $E[(\beta_n-\beta)^k]$ I consider $E[(\beta_n-\beta)^k|\mathbf{X}]$ in what follows.
My real interest is in the case of multiple regression (including an intercept term) but I think it is easier to focus on a simpler case first. Suppose we have the simple intercept-free linear regression model $$ Y_i = \beta X_i + \varepsilon_i, \quad \quad i=1,\dots,n, $$ where $E[\varepsilon_i|\mathbf{X}] = 0$, $\text{Var}[\varepsilon_i^2|\mathbf{X}] = \sigma^2$, and $E[\varepsilon_i\varepsilon_j|\mathbf{X}] = 0$ for $i\neq j$.
It is known that $\hat \beta = \beta + A\varepsilon$ where $$ A = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'. $$ In this simple setting the data matrix $\mathbf{X}$ is actually a vector: $$ \mathbf{X} = [x_1,\dots,x_n]^T. $$ So $A\varepsilon$ is simply a scalar in our simple model.
Now consider evaluating $E[(\hat \beta - \beta)^3|\mathbf{X}]$ or $E[(\hat \beta - \beta)^4|\mathbf{X}]$ using the same process (see p. 61, Econometric Analysis 8 ed. by Green or p.29, Econometrics by Hayashi) that is used for finding an expression for $E[(\hat \beta - \beta)^2|\mathbf{X}]$. That is \begin{align} E[(\hat \beta - \beta)^2|\mathbf{X}] &= E[(A\varepsilon)^2|\mathbf{X}] \\ & = A E[\varepsilon \varepsilon^T|\mathbf{X}]] A^T \\ & = A \sigma^2 A^T \quad \quad \quad \quad \quad \quad \quad \quad \text{(using assumption $E[\varepsilon \varepsilon^T|\mathbf{X}]] = \sigma^2$ )}\\ & = \sigma^2 (\mathbf{X}^T\mathbf{X})^{-1}. \end{align}
Now lets apply this procedure to the higher moments: \begin{align} E[(\hat \beta - \beta)^3|\mathbf{X}] &= E[(A\varepsilon)^3|\mathbf{X}] \\ & = A E[\varepsilon (A \varepsilon) \varepsilon^T|\mathbf{X}]] A^T \\ \end{align}
and
\begin{align} E[(\hat \beta - \beta)^4|\mathbf{X}] &= E[(A\varepsilon)^4|\mathbf{X}] \\ & = A E[\varepsilon (A \varepsilon)^2 \varepsilon^T|\mathbf{X}]] A^T \\ \end{align}
In both cases we get stuck with $A$'s inside the expectation. So is there any way to overcome this and obtain explicit expressions for $E[(\hat \beta - \beta)^3|\mathbf{X}]$ and $E[(\hat \beta - \beta)^4|\mathbf{X}]$?} This issue is overcome in the attempt below.
Further attempt: We have $$ A = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}' = ([x_1, \dots, x_n] [x_1, \dots, x_n]^T)^{-1} [x_1, \dots, x_n] = \frac{1}{||\mathbf{X}||^2}[x_1, \dots, x_n]. $$ Then $$ Ae = \frac{1}{||\mathbf{X}||^2}[x_1, \dots, x_n] [\varepsilon_1, \dots, \varepsilon_n]^T = \frac{1}{||\mathbf{X}||^2} \sum_{i=1}^n x_i e_i. $$ Next we have $$ \begin{align} E[(\hat \beta - \beta)^3|\mathbf{X}] &= E[(A\varepsilon)^3|\mathbf{X}] \\ &= \sum_{i=1}^n E\bigg[\frac{(x_i\varepsilon_i)^3}{||\mathbf{X}||^6}\bigg|\mathbf{X}\bigg] \\ \end{align} $$ where all the other terms have been dropped by assuming $E[e_i e_j^2|\mathbf{X}] = 0$ for $i \neq j$.
Then \begin{align} E[(\hat \beta - \beta)^3|\mathbf{X}] &= \frac{1}{||\mathbf{X}||^6} \sum_{i=1}^n x_i^3 E[\varepsilon_i^3|\mathbf{X}] \\ &= \mu_3 \frac{\sum_{i=1}^n x_i^3}{||\mathbf{X}||^6}. \end{align}
Similarly for $k=4$, if we drop terms by setting conditions on the expectation of the various error terms we get \begin{align} E[(\hat \beta - \beta)^4|\mathbf{X}] &= \mu_4 \frac{\sum_{i=1}^n x_i^4}{||\mathbf{X}||^8}. \end{align}
So assuming these expressions are correct, the pattern for the moments is clear. I now need to get the order of magnitude with respect to $n$ of $$ \begin{align} E[(\hat \beta - \beta)^k] &= E[E[(\hat \beta - \beta)^k|\mathbf{X}]] \\ &= \mu_k E\bigg[\frac{\sum_{i=1}^n x_i^k}{||\mathbf{X}||^{2k}}\bigg] \\ &= \dots \end{align} $$
I am not sure how to evaluate this expectation and find its order of magnitude with respect to $n$ as $n\to \infty$.