# Efficient geometric integrators for nonadiabatic quantum dynamics. I. The adiabatic representation

## TL;DR

This paper introduces implicit geometric integrators for nonadiabatic quantum dynamics that are energy-conserving, unitary, and more accurate than traditional methods, demonstrated on a two-surface NaI model.

## Contribution

The paper develops and analyzes new implicit geometric integrators applicable to non-separable Hamiltonians in the adiabatic representation, enhancing accuracy and efficiency.

## Key findings

- Integrators conserve energy exactly and are symplectic and unitary.
- Numerical tests show significant speedup over traditional methods.
- Higher-order integrators achieve arbitrary accuracy with controlled computational cost.

## Abstract

Geometric integrators of the Schr\"{o}dinger equation conserve exactly many invariants of the exact solution. Among these integrators, the split-operator algorithm is explicit and easy to implement, but, unfortunately, is restricted to systems whose Hamiltonian is separable into a kinetic and potential terms. Here, we describe several implicit geometric integrators applicable to both separable and non-separable Hamiltonians, and, in particular, to the nonadiabatic molecular Hamiltonian in the adiabatic representation. These integrators combine the dynamic Fourier method with recursive symmetric composition of the trapezoidal rule or implicit midpoint method, which results in an arbitrary order of accuracy in the time step. Moreover, these integrators are exactly unitary, symplectic, symmetric, time-reversible, and stable, and, in contrast to the split-operator algorithm, conserve energy exactly, regardless of the accuracy of the solution. The order of convergence and conservation of geometric properties are proven analytically and demonstrated numerically on a two-surface NaI model in the adiabatic representation. Although each step of the higher order integrators is more costly, these algorithms become the most efficient ones if higher accuracy is desired; a thousand-fold speedup compared to the second-order trapezoidal rule (the Crank-Nicolson method) was observed for wavefunction convergence error of $10^{-10}$. In a companion paper [J. Roulet, S. Choi, and J. Van\'{\i}\v{c}ek (2019)], we discuss analogous, arbitrary-order compositions of the split-operator algorithm and apply both types of geometric integrators to a higher-dimensional system in the diabatic representation.