I am interested in solutions to $$ \frac{d}{dt} \Psi = -iH \Psi $$ for $H$ hermitian and time independent. This boils down to evaluating $$ \Psi(t) = e^{-iHt}\Psi_0 $$ at points of interest $t_n$. I want to quickly compute $\Psi(t)$ up to a large final time $T$ by spacing the $t_n$ randomly and advancing the propagator over the nonuniform grid.

To be specific, fix $\Delta t$ and draw $t_n$ uniformly at random from $\Delta t, 2\Delta t, ..., N\Delta t$. The gaps $t_{n+1}-t_{n}$ will be distributed geometrically (ie we will have large gaps).

Any ideas how I could compute $e^{-iHt}$ accurately and efficiently when $t$ may be large?

Fast approximations that I know of (eg Strang splitting, Trotter product) are very efficient but I have found them to be horribly inaccurate and unstable for large gaps.

Machine precision accuracy is achieved by expanding $e^{-iHt}$ in Chebyshev polynomials, but the computational cost (measured by counting applications of $H$) scales linearly with the timestep size and we get to $T$ no faster.

I am beginning to think that spectrally accurate approximations of $e^{-iHt}$ provably require linearly more applications of $H$ as $T$ increases. Does anyone know if this is proven?

Thanks for solving my thesis problem,

Ryan