I'm working on a Real-Time TDDFT implementation and I am currently comparing different propagation schemes for the propagation of the Kohn-Sham wave function, $$ \phi(t+\Delta t) = \hat{\mathcal{U}}\phi(t). $$ An external field is incorporated by a vector potential and the kinetic energy matrix element is modified such that the hamilton matrix gets a complex part: $$ \mathcal{H}_{ij} = \mathcal{H}_{ij}^{(R)} + i \mathcal{Z}_{ij} $$ where $\mathcal{H}^{(R)}$ denotes all the real parts in the hamilton matrix and $\mathcal{Z}$ is a matrix that incorporates the imaginary part of the terms arising due to the vector potential. I think details are not that important here. Among other things, I use the instantaneous total energy $E_\text{tot}(t)$ to compare different propagators$^*$.
Most implemented integration schemes are exponential-based, i.e. $$ \mathcal{U}=\exp\left(f(\mathcal{H})\right), $$ e.g. Exponential Midpoint with $f(\mathcal{H})=-i\mathcal{S}^{-1}\mathcal{H}(t+\Delta t/2)$ or Enforced Time-Reversal Symmetry with $f(\mathcal{H})=-i\Delta t/2(\mathcal{S}^{-1}\mathcal{H}(t+\Delta t)+\mathcal{S}^{-1}\mathcal{H}(t))$ ($\mathcal{S}$ is the overlap matrix). These methods mostly compute the same results but when I use the Crank-Nicolson propagator, $$ \mathcal{U}=\frac{\mathcal{S}-i\mathcal{H}(t)\Delta t/2}{\mathcal{S}+i\mathcal{H}(t)\Delta t/2}, $$ it only reproduces the results when $\mathcal{Z}=0$, else I get (a while after the pulse ends, constant) amplitude oscillations that even persist after a pulse is switched off. Reducing the step size does not change anything.
So my question is what the reason for this could be. I don't know too much about details of applying these propagation schemes to the given kind of problem.
$^*$ $E_\text{tot}(t) = \sum_i f_i \epsilon_i(t) - F[\rho(t)] + E_{nn}$ where $\epsilon_i$ are the eigenvalues of the time-dependent hamiltonian, $F$ denotes some functionals incorporating XC and electrostatic corrections and $E_{nn}$ is the internuclear repulsion energy.
