Why Runge-Kutta for Quaternion integration in Kalman filter?
I think the confusion comes from the authors not parameterizing things clearly.  Furthermore, by switching to geometric algebra rather than quaternions, some additional confusion can be cleared up.

The main difference between normal vector algebra and geometric algebra is that we can multiply vectors. So of $e_x$, $e_y$, and $e_z$ are our (orthonormal) basis vectors, we also have $e_x e_y$ which is not a vector but a bivector which can be thought of as an oriented plane element.  (By "pure vector quaternion" they mean a bivector; you are correct that that means $0$ $w$ part and it's reasonable to say it has a $0$ "real" part.)  A key property is that the basis vectors anti-commute with each other, i.e. $e_x e_y = -e_y e_x$.  This leads to $(e_x e_y)^2 = -1$.  In general, any unit bivector squares to $-1$, which means we can apply Euler's formula: $$R = e^{\theta B} = \cos\theta + B\sin\theta$$ where $B$ is a unit bivector and we call $R$ a rotor.  A unit quaternion is just a 3D rotor.  (Complex numbers are just the even sub-algebra of the 2D geometric algebra.  Here we are looking at quaternions, the even sub-algebra of 3D GA.)

If we want a rotation by $\theta$ in the plane of a unit bivector $B$ we use the rotor $R = e^{\frac{\theta}{2}B}$.  (The half comes in because we rotate a vector $v$ via $Rv\widetilde{R}$; see the link above.)  Of course, we want to allow both the plane and the angle to vary so we define the (non-unit) bivector $\Theta(t) = \theta(t)B(t)$ so if $$R(t) = e^{\frac{1}{2}\Theta(t)}$$ then $$\dot{R}(t) = \frac{1}{2}\Omega(t)R(t)$$ where $\Omega = \dot{\Theta}$.  This is their $f = q\omega/2$; their $q$ is our $R$ and their $\omega$ is our $\Omega$.  Now if we have $R(t_0)$ and we want $R(t_0 + \Delta{}t)$ we need to integrate $\dot{R}$ from $t_0$ to $t_0 + \Delta{}t$.

The scheme you describe is roughly analogous to doing forward Euler integration, which essentially assumes their $q$ and $\omega$ our $R$ and $\Omega$ is constant over $\Delta{}t$.  RK4 is just a better integration method.  I doubt there is any special reason they chose to use RK4 as opposed to other integration methods.  It's just the default choice typically.