can we operator split within VODE?

Right now we solve the system:

$X^{n+1}_k = \dot{\omega}_k (\rho, T^{n+1}, X_k^{n+1})$
$e^{n+1} = \dot{S} (\rho, T^{n+1}, X_k^{n+1})$

but what if we treated the energy/temperature update explicitly. We are substepping in VODE, so we could still be accurate / stable (?) if we did:

$X^{n+1}_k = \dot{\omega}_k (\rho, T^{n}, X_k^{n+1})$
$e^{n+1} = \dot{S} (\rho, T^{n}, X_k^{n+1})$

the advantage of this is that the Jacobian then does not need any temperature derivatives, significantly simplifying it and making issues like #729 a non-issue.

Start by locating the VODE implementation and its current coupled energy/species update and Jacobian path. Analyze whether the proposed explicit temperature treatment preserves the required accuracy and stability, and document the resulting design and validation criteria for the operator split.