Extended Kalman filter - why not use a 2nd order Taylor expansion? I'm implementing my first EKF and using a Taylor expansion to linearize the model. I'm somewhat new to this. Based on a paper I've read, I'm attempting to linearize the state transition matrix using a Taylor Expansion.

AFAIU, given a state vector X and state transition matrix F, we can write a linearized state transition matrix Φ (probably not the name) in the following manner:

$$
\textbf{F} = \frac{\partial \dot{\textbf{X}}}{\partial \textbf{X}}
$$
$$
\Phi = \textbf{I} + \textbf{F}\Delta t + \textbf{F}^2 \frac{\Delta t^2}{2!} + H.O.T
$$

Such that the state and covariance estimates are:

$$
\hat{\textbf{X}}_k = \Phi \hat{\textbf{X}}_{k-1}
$$

$$
\hat{\textbf{P}}_k = \Phi \hat{\textbf{P}}_{k-1} \Phi^T + \textbf{Q}
$$

Using a 2nd order expansion instead of a 1st, I've found that I get a much more accurate solution. It also allows me to use a larger timestep (minutes instead of seconds). I worry though that in doing so, I'm not actually linearizing, since the 2nd order TSE corresponds to a quadratic approximation.

Assuming the above steps are correct, is there anything wrong with taking the 2nd order TSE?
There is nothing wrong with taking a higher order approximation than first order.
An issue in doing so is that you no longer have Gaussian statistics, so some of the standard Kalman filter equations (which are derived using Gaussianity) no longer hold strictly true.
The usual approach to dealing with (otherwise ignored) higher order terms in the standard EKF is to increase either the process noise or the measurement noise variances from what the model tells you they should be. This is like saying,

I know what the noise power is, but I'm going to say the second and higher order terms in my Taylor series expansion add to it so I'm going to increase the noise variance I use in my filter to cater for it.