The original split operator paper from Feit&Fleck actually did not do propagation, but tried to calculate eigenstates of the Hamiltonian. The basic idea is very simple: Calculate the wavefunction psi(t) on a dense time grid, and Fourier-transform this series to get some psi(omega), which for a given eigenfrequency/-energy omega gives an approximation to the eigenstate at this energy.
In practice, things are a bit more tricky, with degenerate or almost degenerate states (you may want to orthogonalize w.r.t. some states), and of course you need to have a good guess for the energy to start with, because you cannot store psi(t), but need to Fourier-transform immediately. Maybe you also want to do this iteratively, refining an initial guess.
In any case, this sounds like a pretty pecular but fun scheme to get arbitrary eigenstates with high accuracy. Could be implemented and tested out.