# Crank-Nicolson integrator: oscillations with complex matrix

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.

• This is likely due to L-stability. People think of the Crank-Nicholson scheme as "good for stiff" because it's A-stable, but that's a very weak form of stability. In reality, the production-quality ODE solvers for stiff equations all have L-stability, which is the dampening of oscillatory modes (at infinity). This might be an example of where L-stability is required instead of just A-stability. – Chris Rackauckas Nov 24 '17 at 14:14
• Thanks alot for your hint! I will have a look at that. Could you maybe recommend literature for this topic? – Lukk Dec 6 '17 at 9:06
• Hairer's second book on ODE methods (on stiff solvers) has a very good treatment of this topic. – Chris Rackauckas Dec 6 '17 at 15:13