The time-dependent Schroedinger equation (TDSE) depends linearly on the system's initial state $\vert \psi(0) \rangle$, such that the solution can be generally written as
$$ \vert \psi(t) \rangle = \exp \left( \frac{i}{\hbar} \int_0^t \text{d}t' \widehat{H(t')} \right) \vert \psi(0) \rangle \quad , \quad (1) $$
where the [possibly time-dependent] Hamiltonian operator is given by $\widehat{H(t)}$.
A common technique to numerically solve the above equation is to use the Zassenhaus formula to approximately split the exponential for an infinitesimal time $\Delta t$. The error one gets depends on how this splitting is performed: simple schemes, such as the famous leapfrog method, have an [energy] error of $\mathcal{O}(\Delta t^2)$. Yoshida, Strang, Feit and others showed that smart splitting schemes can produce methods of arbitrary accuracy, albeit paying the price in computational cost. These splitting methods are known as "spectral methods". After choosing a splitting scheme the final wavefunction is obtained recursively and, for separable Hamiltonians expressible as a sum of kinetic plus potential energies, solved by moving back and forth in momentum and position representations by the repeated use of Fast Fourier Transforms (FFTs).
A much more general scheme is the Crank-Nicholson method, which in principle can be used to solve equations that are a lot more complicated than the TDSE and is a stable, also $\mathcal{O}(\Delta t^2)$ [in energy] method. In general, the solution to these equations cannot be written in the simple form (1), and spectral methods such as the ones described above are not available. In these cases it makes sense to think about Crank-Nicolson, but for the TDSE there is a plethora of super sophisticated/efficient spectral schemes available.
Why, then, do I keep stumbling into papers that employ Cranck-Nicholson to solve the TDSE?
