The Magnus Expansion
The goal of all real-time electronic dynamics methods is to solve the time-dependent Schrödinger equation (TDSE)
\[i\hbar \frac{\partial}{\partial t} \psi (t) = H(t) \psi(t)\]where \(H(t)\) is the time-dependent Hamiltonian and \(\psi(t)\) is the time-dependent wave function. The goal of the Magnus expansion is to find a general solution for the time-dependent wave function in the case where \(H\) is time-dependent, and, more crucially, when \(H\) does not commute with itself at different times, e.g. when \(\left[H(t_1),H(t_2)\right] \neq 0\). In the following we will follow closely the notation of Blanes, et al..
First, for simplicity we redefine \(\tilde{H}(t) \equiv \frac{-i}{\hbar}H(t)\) and introduce a scalar \(\lambda = 1\) as a bookkeeping device, so that
\[\begin{aligned} \frac{\partial}{\partial t} \psi (t) &= \lambda \frac{-i}{\hbar} H(t) \psi(t) \\ &= \lambda \tilde{H}(t) \psi(t) \label{eq:tdse}\end{aligned}\]At the heart of the Magnus expansion is the idea of solving the TDSE by using the quantum propagator \(U(t,t_0)\) that connects wave functions at different times, e.g. \(\psi(t) = U(t,t_0)\psi(t_0)\) Furthermore, the Magnus expansion assumes that \(U(t,t_0)\) can be represented as an exponential, \(U(t,t_0) = \text{exp} \left( \Omega(t,t_0) \right)\) This yields the modified TDSE
\[\label{eq:modtdse} \frac{\partial}{\partial t} U(t, t_0) = \lambda \tilde{H}(t) U(t,t_0); \quad U(t_0, t_0) = I\]Now, for scalar \(H\) and \(U\), the above has a simple solution, namely
\[U(t,t_0) = \text{exp} \left( \lambda \int\limits^{t}_{t_0} \tilde{H}(t') dt' \right)\]However, if \(H\) and \(U\) are matrices this is not necessarily true. In other words, for a given matrix \(A\) the following expression does not necessarily hold:
\[\frac{\partial}{\partial t } \left( \text{exp}\left(A(t)\right) \right) = \left(\frac{\partial}{\partial t }A(t) \right)\text{exp}\left(A(t)\right) = \text{exp}\left(A(t)\right) \left(\frac{\partial}{\partial t } A(t) \right)\]because the matrix \(A\) and its derivatives do not necessarily commute. Instead, Magnus proved that in general \(\Omega (t, t_0)\) satisfies
\[\frac{\partial}{\partial t} \left( \Omega(t,t_0) \right) = \lambda \tilde{H}(t) + \lambda \sum\limits_{k=1}^{\infty} (-1)^k \frac{B_k}{k!}\overbrace{[\Omega(t,t_0), [ \cdots[\Omega(t,t_0),}^{k-\text{times}} \tilde{H}(t)]] \cdots ]; \quad \Omega(t_0,t_0) = 0\]where \(B_k\) are the Bernoulli numbers. This equation may be solved by integration, and iterative substitution of \(\Omega(t,t_0)\). While it may appear that we are worse off than when we started, collecting like powers of \(\lambda\) (and setting \(\lambda = 1\)) allows us to obtain a power-series expansion for \(\Omega(t,t_0)\),
\[\begin{aligned} \label{eq:magnus} \Omega(t,t_0) &= \int\limits_{t_0}^{t} \tilde{H}_1 dt_1 \\ &+ \frac{1}{2} \int\limits_{t_0}^{t} dt_1 \int\limits_{t_0}^{t_1} dt_2 \left[\tilde{H}_1,\tilde{H}_2\right] \\ &+ \frac{1}{6} \int\limits_{t_0}^{t} dt_1 \int\limits_{t_0}^{t_1} dt_2 \int\limits_{t_0}^{t_2} dt_3 \left( [\tilde{H}_1, [\tilde{H}_2, \tilde{H}_3]] + [\tilde{H}_3,[\tilde{H}_2,\tilde{H}_1]] \right) + \cdots\end{aligned}\]This is the Magnus expansion, and here we have given up to the third-order terms. We have also made the notational simplification that \(\tilde{H}_k = \tilde{H}(t_k)\). This is the basis for nearly all numerical methods to integrate the many-body TDSE in molecular physics. Each subsequent order in the Magnus expansion is a correction that accounts for the proper time-ordering of the Hamiltonian.
The Magnus expansion immediately suggests a route to many numerical integrators. The simplest would be to approximate the first term by
\[\int\limits_{t}^{t + \Delta t} \tilde{H}_1 dt_1 \approx \Delta t \tilde{H}(t)\]leading to a forward-Euler-like time integrator of
\[\psi(t + \Delta t) = \text{exp}\left(\Delta t \tilde{H}(t)\right)\psi(t)\]which we can re-write as
\[\psi(t_{k+1}) = \text{exp}\left(\Delta t \tilde{H}(t_{k})\right)\psi(t_{k})\]where subscript \(k\) gives the node of the time-step stencil. This gives a first-order method with error \(\mathcal{O}(\Delta t)\). A more accurate second-order method can be constructed by approximating the first term in the Magnus expansion by the midpoint rule, leading to an \(\mathcal{O}({\Delta t}^2)\) time integrator
\[\psi(t_{k+1}) = \text{exp}\left(\Delta t \tilde{H}(t_{k+1/2})\right)\psi(t_{k})\]Modifying the stencil to eliminate the need to evaluate the Hamiltonian at fractional time steps (e.g. change time step to \(2 \Delta t\)) leads to the modified midpoint unitary transformation (MMUT) method
\[\psi(t_{k+1}) = \text{exp}\left(2 \Delta t \tilde{H}(t_{k})\right)\psi(t_{k-1})\]which is a leapfrog-type unitary integrator. Note that the midpoint method assumes \(\tilde{H}\) is linear over its time interval, and the higher order terms (containing the commutators) in this approximation go to zero. There are many other types of integrators based off the Magnus expansion that can be found in the literature. The key point for all of these integrators is that they are symplectic, meaning they preserve phase-space relationships. This has the practical effect of conserving energy (within some error bound) in long-time dynamics, whereas non-symplectic methods such as Runge-Kutta will experience energetic “drift” over long times. A final note: in each of these schemes it is necessary to evaluate the exponential of the Hamiltonian. In real-time methods, this requires computing a matrix exponential. This is not a trivial task, and, aside from the construction of the Hamiltonian itself, is often the most expensive step in the numerical solution of the TDSE. However, many elegant solutions to the construction of the matrix exponential can be found in the literature.
References
Blanes, S., Casas, F., Oteo, J.A. and Ros, J., 2010. A pedagogical approach to the Magnus expansion. European Journal of Physics, 31(4), p.907.
Magnus, W., 1954. On the exponential solution of differential equations for a linear operator. Communications on pure and applied mathematics, 7(4), pp.649-673.