The Second Order Differencing scheme (SOD) [3]

The most natural ansatz to approximating the propagation operator

$$\widehat{U}=e^{-\frac{i}{\hbar }\widehat{H}\Delta t}=1-\frac{i}{\hbar }\widehat{H}\Delta t+\cdots$$

would be the straightforward expansion into a Taylor series. This has, however, proven to be numerically unstable because of its asymmetry with respect to time inversion. Therefore the method has to be symmetrised by considering one forward and one backward step to first order in the Taylor expansion:

$$\textstyle \psi (t+\Delta t)-\psi (t-\Delta t) = \left( e^{-\frac{i}{\hbar }\widehat{H}\Delta t}-e^{\frac{i}{\hbar }\widehat{H}\Delta t}\right) \psi (t),$$

$$\rightarrow \, \psi (t+\Delta t) \approx \psi (t-\Delta t)-2\frac{i}{\hbar }\Delta t\widehat{H}\psi (t).$$ (6)

There are two ways of obtaining the second wavefunction $\psi (\Delta t)$ from $\psi (0)$:

1. (SOD) Propagate by a first order scheme for half a time-step and from there propagate with SOD for another half time-step.
2. (SODS) Propagate with SOD half a step forward to get $\psi \left( \frac{\Delta t}{2}\right)$ and half a step backward to get $\psi \left( -\frac{\Delta t}{2}\right)$. The final result is the arithmetic mean value of $\psi \left( t+\frac{\Delta t}{2}\right)$ and $\psi \left( t-\frac{\Delta t}{2}\right)$. This method is more symmetric.
3. $\begin{array}{ccc} \psi (0) & \rightarrow (1st\, order)\rightarrow & \psi \lef... ...)\rightarrow & \psi (0),\, \psi \left( \frac{\Delta t}{2}\right) . \end{array}$
   this initialization stops there because the lower two steps become cyclic because of the time-reversability of the SOD method.

The method is unitary and conserves norm and energy. It can also be very easily and naturally adjusted to involve interactions between potential surfaces a and b due to the electromagnetic field E(t) and the (approximately constant) magnetic dipole moment $\mu$: [7]

$$\psi _{a}(t+\Delta t) = \psi _{a}(t-\Delta t)-2\frac{i}{\hbar }\Delta t\widehat{H}_{a}\psi _{a}(t)-2\frac{i}{\hbar }\Delta t\mu E(t)\psi _{b}(t),$$

$$\psi _{b}(t+\Delta t) = \psi _{b}(t-\Delta t)-2\frac{i}{\hbar }\Delta t\widehat{H}_{b}\psi _{b}(t)-2\frac{i}{\hbar }\Delta t\mu E(t)\psi _{a}(t).$$