# 3.2: Integrating the Schrödinger Equation Directly

$$\newcommand{\vecs}{\overset { \rightharpoonup} {\mathbf{#1}} }$$ $$\newcommand{\vecd}{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}}$$$$\newcommand{\id}{\mathrm{id}}$$ $$\newcommand{\Span}{\mathrm{span}}$$ $$\newcommand{\kernel}{\mathrm{null}\,}$$ $$\newcommand{\range}{\mathrm{range}\,}$$ $$\newcommand{\RealPart}{\mathrm{Re}}$$ $$\newcommand{\ImaginaryPart}{\mathrm{Im}}$$ $$\newcommand{\Argument}{\mathrm{Arg}}$$ $$\newcommand{\norm}{\| #1 \|}$$ $$\newcommand{\inner}{\langle #1, #2 \rangle}$$ $$\newcommand{\Span}{\mathrm{span}}$$ $$\newcommand{\id}{\mathrm{id}}$$ $$\newcommand{\Span}{\mathrm{span}}$$ $$\newcommand{\kernel}{\mathrm{null}\,}$$ $$\newcommand{\range}{\mathrm{range}\,}$$ $$\newcommand{\RealPart}{\mathrm{Re}}$$ $$\newcommand{\ImaginaryPart}{\mathrm{Im}}$$ $$\newcommand{\Argument}{\mathrm{Arg}}$$ $$\newcommand{\norm}{\| #1 \|}$$ $$\newcommand{\inner}{\langle #1, #2 \rangle}$$ $$\newcommand{\Span}{\mathrm{span}}$$

Okay, how do we evaluate the time-propagator and obtain a time-dependent trajectory for a quantum system? Expressions such as the time-ordered exponentials are daunting, and there are no simple ways in which to handle this. One cannot truncate the exponential because usually this is not a rapidly converging series. Also, the solutions oscillate rapidly as a result of the phase acquired at the energy of the states involved, which leads to a formidable integration problem. Rapid oscillations require small time steps, when in fact the time scales. For instance in a molecular dynamics problem, the highest frequency oscillations may be as a result of electronically excited states with periods of less than a femtosecond, and the nuclear dynamics that you hope to describe may occur on many picosecond time scales. Rather than general recipes, there exist an arsenal of different strategies that are suited to particular types of problems. The choice of how to proceed is generally dictated by the details of your problem, and is often an art-form. Considerable effort needs to be made to formulate the problem, particularly choosing an appropriate basis set for your problem. Here it is our goal to gain some insight into the types of strategies available, working mainly with the principles, rather than the specifics of how it’s implemented.

Let’s begin by discussing the most general approach. With adequate computational resources, we can choose the brute force approach of numerical integration. We start by choosing a basis set and defining the initial state $$\psi_0$$. Then, we can numerically evaluate the timedependence of the wavefunction over a time period $$t$$ by discretizing time into $$n$$ small steps of width $$\delta t = t / n$$ over which the change of the system is small. A variety of strategies can be pursed in practice.

One possibility is to expand your wavefunction in the basis set of your choice

$| \psi (t) \rangle = \sum _ {n} c _ {n} (t) | \varphi _ {n} \rangle \label{2.32}$

and solve for the time-dependence of the expansion coefficients. Substituting into the right side of the TDSE,

$i \hbar \frac {\partial} {\partial t} | \psi \rangle = \hat {H} | \psi \rangle \label{2.33}$

and then acting from the left by $$\langle k |$$ on both sides leads to an equation that describes their time dependence:

$i \hbar \frac {\partial c _ {k} (t)} {\partial t} = \sum _ {n} H _ {k n} (t) c _ {n} (t) \label{2.34}$

or in matrix form $$i \hbar \dot {c} = H c$$. This represents a set of coupled first-order differential equations in which amplitude flows between different basis states at rates determined by the matrix elements of the time-dependent Hamiltonian. Such equations are straightforward to integrate numerically. We recognize that we can integrate on a grid if the time step forward ($$\delta t$$) is small enough that the Hamiltonian is essentially constant. Then Equation \ref{2.34} becomes

$i \hbar \delta c _ {k} (t) = \sum _ {n} H _ {k n} (t) c _ {n} (t) \delta t \label{2.35}$

and the system is propagated as

$c _ {k} ( t + \delta t ) = c _ {k} (t) + \delta c _ {k} (t) \label{2.36}$

The downside of such a calculation is the unusually small time-steps and significant computational cost required.

Similarly, we can use a grid with short time steps to simplify our time-propagator as

$\hat {U} ( t + \delta t , t ) = \exp \left[ - \frac {i} {\hbar} \int _ {t}^{t + \delta t} d t^{\prime} \hat {H} \left( t^{\prime} \right) \right] \approx \exp \left[ - \frac {i} {\hbar} \delta t \hat {H} (t) \right] \label{2.37}$

Therefore the time propagator can be written as a product of $$n$$ propagators over these small intervals.

\begin{align} \hat {U} (t) & = \lim _ {\delta t \rightarrow 0} \left[ \hat {U} _ {n} \hat {U} _ {n - 1} \cdots \hat {U} _ {2} \hat {U} _ {1} \right] \label{2.38A} \\[4pt] & = \lim _ {n \rightarrow \infty} \prod _ {j = 0}^{n - 1} \hat {U} _ {j} \label{2.38B} \end{align} Here the time-propagation over the jth small time step is

\left.\begin{aligned} \hat {U} _ {j} & = \exp \left[ - \frac {i} {\hbar} \delta t \hat {H} _ {j} \right] \\[4pt] \hat {H} _ {j} & = \hat {H} ( j \delta t ) \end{aligned} \right. \label{2.39}

Note that the expressions in Equations \ref{2.38A} and \ref{2.38B} are operators time ordered from right to left, which we denote with the “+” subscript. Although Equation \ref{2.38B} is exact in the limit $$\delta t \rightarrow 0$$ (or $$n→∞$$), we can choose a finite number such that $$H(t)$$ does not change much over the time $$\delta t$$. In this limit the time propagator does not change much and can be approximated as an expansion

$\hat {U} _ {j} \approx 1 - \frac {i} {\hbar} \delta t \hat {H} _ {j} \label{2.40}.$

In a general sense this approach is not very practical. The first reason is that the time step is determined by $$\delta \mathrm {t} < \hbar / | H |$$ which is typically very small in comparison to the dynamics of interest. The second complication arises when the potential and kinetic energy operators in the Hamiltonian don’t commute. Taking the Hamiltonian to be $$\hat {H} = \hat {T} + \hat {V}$$

\left.\begin{aligned} e^{- i \hat {H} (t) \delta t / h} & = e^{- i ( \hat {T} (t) + \hat {V} (t) ) \delta t / h} \\[4pt] & \approx e^{- i \hat {T} (t) \delta t / \hbar} e^{- i \hat {V} (t) \delta t / h} \end{aligned} \right. \label{2.41}

The second line makes the Split Operator approximation, what states that the time propagator over a short enough period can be approximated as a product of independent propagators evolving the system over the kinetic and potential energy. The validity of this approximation depends on how well these operators commute and the time step, with the error scaling like $$\frac {1} {2} [ \hat {T} (t) , \hat {V} (t) ] ( \delta t / \hbar )^{2}$$ meaning that we should use a time step, such that $$\delta t < \left\{2 \hbar^{2} / [ \hat {T} (t) , \hat {V} (t) ] \right\}^{1 / 2}$$

This approximation can be improved by symmetrizing the split operator as

$e^{- i \hat {H} (t) \delta t / h} \approx e^{- i \hat {V} (t) \frac {\delta t} {2} / h} e^{- i \hat {T} (t) \delta t / h} e^{- i \hat {V} (t) \frac {\delta t} {2} / h} \label{2.42}$

Here the error scales as $$\frac {1} {12} ( \delta t / \hbar )^{3} \left\{[ \hat {T} , [ \hat {T} , \hat {V} ] ] + \frac {1} {2} [ \hat {V} , [ \hat {V} , \hat {T} ] ] \right\}$$. There is no significant increase in computational effort since half of the operations can be combined as

$e^{- \frac {i \hat {V}} {h} \frac {( j + 1 ) \delta t} {2}} e^{- \frac {i \hat {V}} {\hbar} \frac {j \delta t} {2}} \approx e^{- i \hat {V} j \delta t / \hbar}$

to give $$U (t) \approx e^{- i \hat {V} \frac {n \delta t} {2} / h} \left[ \prod _ {j = 1}^{n} e^{- i \hat {V} j \delta t / h} e^{- i \hat {T} j \delta t / h} \right] e^{- i \hat {V} \frac {\delta t} {2} / h} \label{2.44}$$