Numerical solution of the linear semiclassical Schrödinger equation on the real line

This paper is concerned with the numerical solution of the linear Schrödinger equation in the semiclassical regime, describing the motion of an electron in a quantum system,

$${\mathrm{i}}\varepsilon\frac{\partial u}{\partial t}=-\varepsilon^{2}\Delta u+V(\mbox{\boldmath$x$\unboldmath})u,\qquad t\geq 0,\;\mbox{\boldmath$x$\unboldmath}\in\mbox{\Bbb R}^{d},$$

(1.1)

where the initial condition $u(\mbox{\boldmath$x$\unboldmath},0)=u_{0}(\mbox{\boldmath$x$\unboldmath})\in\mathrm{L}_{2}(\mbox{\Bbb R}^{d})$ for all $\mbox{\boldmath$x$\unboldmath}\in\mbox{\Bbb R}^{d}$ is given. The semiclassical parameter $\varepsilon>0$ is a small number which describes the square root of the ratio between the mass of an electron and the total mass of the system, and $V:\mbox{\Bbb R}^{d}\rightarrow\mbox{\Bbb R}$ is the interaction potential which is assumed to be smooth for the purposes of this paper. Since $|u(\mbox{\boldmath$x$\unboldmath},t)|^{2}$ gives the probability density of the electron residing at $x$ at time $t$, the system is required to satisfy,

$$\int_{\mbox{\sBbb R}^{d}}|u(\mbox{\boldmath$x$\unboldmath},t)|^{2}\,\mathrm{d}\mbox{\boldmath$x$\unboldmath}\equiv 1,$$

(1.2)

and any physically relevant numerical solution must be consistent with this conservation law. To read more, (?) is an excellent, up to date review of both the equation (1.1) and its numerical solution.

Respecting the unitarity property (1.2) underlies the importance of geometric numerical integration methodologies in this context and has been central to modern treatment of the linear Schrödinger equation in the semiclassical, $0<\varepsilon\ll 1$, and the atomistic, $\varepsilon=1$, regimes alike (?, ?, ?, ?, ?). However, all these publications are focussed on a subtly different problem: instead of being defined in $\mbox{\Bbb R}^{d}$, the equation (1.1) is set on a torus, typically $\mbox{\Bbb T}^{d}$, with periodic boundary conditions. This is of crucial importance to splitting techniques, a common denominator to all these methodologies, because the free Schrödinger equation

$${\mathrm{i}}\frac{\partial u}{\partial t}=-\varepsilon\Delta u,$$

(1.3)

given with periodic boundary conditions, can be approximated very rapidly, affordably and precisely by means of the Fast Fourier Transform (FFT).

Our contention is that the periodic setting imposes unwelcome limitations on the solution, which might lead to altogether false outcomes, and this becomes problematic once a solution over a long time interval is sought (e.g. in quantum control). The underlying reason is the tension arising from the nature of the differential equation and of the initial condition, both predicated by quantum-mechanical considerations. The differential equation itself is dispersive: different waves travel at different speeds, dependent on their wavelengths, which can span a very wide range, all the way from ${\cal O}\!\left(1\right)$ to ${\cal O}\!\left(\varepsilon^{-1}\right)$. The initial condition is typically a linear combination of highly localised (and rapidly oscillating) wave packets. Recall that $|u(\mbox{\boldmath$x$\unboldmath},t)|^{2}$ represents the probability of a particle residing at $x$ in time $t$: while it is a central tenet of quantum mechanics that a particle cannot be completely localised, typically $|u(\mbox{\boldmath$x$\unboldmath},t)|^{2}$ is a linear combination of narrowly-concentrated Gaussian-like structures. These Gaussian-like structures travel at different speeds and, provided the equation is solved for sufficiently long time, some of them eventually reach the boundary. At this very moment periodicity becomes a foe because the wave packet reaches the boundary and ‘pops out’ at the other end — this is not physical!

An alternative to periodic boundary conditions is to impose zero Dirichlet or zero Neumann boundary conditions. However, the following argument shows that this approach is also problematic. Consider an initial condition $u_{0}\in\mathrm{H}^{1}_{0}(0,1)$ and potential $V\in\mathrm{H}^{1}(0,1)$. Now consider the following two initial boundary value problems, the first of which has zero Dirichlet boundary conditions, the second of which has periodic boundary conditions:

	$\displaystyle{\mathrm{i}}\varepsilon\frac{\partial v}{\partial t}$	$\displaystyle=$	$\displaystyle-\varepsilon^{2}\frac{\partial^{2}v}{\partial x^{2}}+V(x)v,\qquad x\in[0,1],$		(1.4)
	$\displaystyle v(0,t)$	$\displaystyle=$	$\displaystyle 0,\qquad v(1,t)=0,\qquad t>0,$
	$\displaystyle v(x,0)$	$\displaystyle=$	$\displaystyle u_{0}(x),\qquad x\in[0,1],$

and,

	$\displaystyle{\mathrm{i}}\varepsilon\frac{\partial w}{\partial t}$	$\displaystyle=$	$\displaystyle-\varepsilon^{2}\frac{\partial^{2}w}{\partial x^{2}}+V(|x|)w,\qquad x\in[-1,1],$		(1.5)
	$\displaystyle w(-1,t)$	$\displaystyle=$	$\displaystyle w(1,t),\qquad\partial_{x}w(-1,t)=\partial_{x}w(1,t),\qquad t>0,$
	$\displaystyle w(x,0)$	$\displaystyle=$	$\displaystyle\mathrm{sign}(x)u_{0}(|x|),\qquad x\in[-1,1].$

The relationship between $v(x,t)$ and $w(x,t)$ for $x\in[0,1]$ and $t>0$ is rather simple. Clearly the oddness of $w(x,0)$ is preserved since the second derivative and multiplication by $V(|x|)$ preserve oddness. Combining oddness with periodicity implies that $w(0,t)=0=w(1,t)$ for all time. It therefore follows from uniqueness of solution to (1.4) that $w(x,t)=v(x,t)$ for $x\in[0,1]$ and $t>0$. So now let us return to the notion of a wave packet moving towards the boundary, but this time with zero Dirichlet boundary conditions imposed. The solution to the odd extension implies that this wave packet will be reflected back and its sign reversed — while this physically happens in the case of an infinite potential barrier, it is not the correct behaviour when posed in free space! A similar construction can be made for Neumann boundary conditions.

We hope this has convinced the reader: no matter what we do, and no matter how rapidly and accurately we can solve Schrödinger’s equation posed on a bounded set, the result of truncating the domain from $\mbox{\Bbb R}^{d}$ to such a set destroys the physics of the problem over a large enough time interval. This is the raison d’être for this paper: solve (1.1) without compromising its setting in $\mbox{\Bbb R}^{d}$. Throughout the remainder of the paper we assume that (1.1) is presented in a single space dimension, $d=1$. A generalisation to a modest number of space dimensions can be accomplished with tensor products along the lines of (?), while generalisation to a large number of dimensions would require a raft of additional techniques and is beyond the scope of the current paper.

To achieve this aim, we will extend the framework of symmetric Zassenhaus splittings, which has been developed for (1.1) on the torus 𝕋 (?, ?), to (1.1) posed on the whole real line. This is not a straightforward exercise, because we cannot use special properties of the Fourier basis. In Section 2 we derive these Zassenhaus splittings under more general assumptions, allowing for bases other than the Fourier basis to be used. In Section 3, we discuss the solution of the free Schrödinger equation (1.3), focusing on two bases which are orthonormal on the real line: Hermite functions and Malmquist–Takenaka functions. In Section 5 we demonstrate how these pieces can be put together to construct practical numerical solvers on the real line.

As it was shown in (?) Schrödinger equations in semiclassical regime produce oscillations in space of frequency ${\cal O}\!\left(\varepsilon^{-1}\right)$, which places restrictions on the discretisation in space depending on which basis is used, because we must employ sufficiently fine discretisation to resolve these oscillations. If the spatial variable is discretised using the Fourier basis then this necessitates $\mathcal{O}(\varepsilon^{-1})$ basis elements, which in turn, leads to the conclusion that after discretisation, operators of type $\partial_{x}^{n}$ have a spectral radius which scales like ${\cal O}\!\left(\varepsilon^{-n}\right)$. As we discuss in Section 3, for other bases it is not necessarily the case that $\mathcal{O}(\varepsilon^{-1})$ basis elements can resolve spatial oscillations of frequency ${\cal O}\!\left(\varepsilon^{-1}\right)$ (indeed the Fourier basis is the optimal basis for resolving periodic oscillations). As such, we will not make assumptions about the number of basis elements, but rather, make assumptions directly on the spectral radius of the partial derivative operator (an assumption which holds in both of our examples discussed in Section 3)

Since the potential $V(x)$ can in principle be an unbounded function on the real line, we must be careful that our expansions be treated locally in $x$.

We can now make sense of “h.o.t.” in the sBCH formula by bounding the magnitude element of the Hall basis for the free Lie algebra generated by $A=-\varepsilon^{-1}V$ and $B=\varepsilon\partial_{x}^{2}$ (i.e. $A$, $B$, $[A,B]$, $[[B,A],A]$, $[[B,A],B]$ $\ldots$ (?, ?)).

Before we proceed with the proof of the theorem, note that all elements of Hall basis (?) of commutators constructed of letters $A$ and $B$ live in the set

$$\mbox{\gothic G}=\left\{\sum_{k=0}^{K}y_{k}(x)\partial_{x}^{k}\,:K\in\mbox{\Bbb Z}_{+},y_{0},\ldots,y_{K}\in\mathrm{C}^{\infty}_{\mathrm{loc}}(\mbox{\Bbb R})\right\},$$

by applying the product rule (for differentiation). For example,

	$\displaystyle[B,A]=-[\partial_{x}^{2},V]=-\left(V^{(2)}+2V^{(1)}\partial_{x}+V\partial_{x}^{2}-V\partial_{x}^{2}\right)=-V^{(2)}-2V^{(1)}\partial_{x}$
	$\displaystyle[[B,A],A]=\varepsilon^{-1}[[\partial_{x}^{2},V],V]=\varepsilon^{-1}2(V^{(1)})^{2}$
	$\displaystyle[[B,A],B]=-\varepsilon[[\partial_{x}^{2},V],\partial_{x}^{2}]=\varepsilon\left(V^{(4)}+4V^{(3)}\partial_{x}+4V^{(2)}\partial_{x}^{2}\right),$
	$\displaystyle[[[B,A],A],A]=0,$

where $V^{(k)}=\partial_{x}^{k}V$. We define the height of the commutator $C$ as the largest index of non-zero coefficient $y_{K}(x)$:

$${\rm ht}(C)={\rm ht}\!\left(\sum_{k=0}^{K}y_{k}(x)\partial_{x}^{k}\right)=K,\text{ where }y_{K}(x)\not\equiv 0,$$

One can observe, that ${\rm ht}(A)=0$, ${\rm ht}(B)=2$, ${\rm ht}([B,A])=1$, ${\rm ht}([[B,A],A])=0$ and ${\rm ht}([[B,A],B])=2$.

In the proof we will also refer to the formula elaborated in (?)

	$\displaystyle\left[\sum_{i=0}^{n}f_{i}(x)\partial^{i}_{x},\sum_{j=0}^{m}g_{j}(x)\partial^{j}_{x}\right]=$	$\displaystyle\sum_{i=0}^{n}\sum_{j=0}^{m}\sum_{\ell=0}^{i}\binom{i}{\ell}f_{i}(x)\left(\partial^{i-\ell}_{x}g_{j}(x)\right)\partial^{\ell+j}_{x}$
		$\displaystyle\mbox{}-\sum_{j=0}^{m}\sum_{i=0}^{n}\sum_{\ell=0}^{j}\binom{j}{\ell}g_{j}(x)\left(\partial^{j-\ell}_{x}f_{i}(x)\right)\partial^{\ell+i}_{x}.$		(2.9)

Proof (of Theorem 3) Let us assume, that a certain non-zero commutator $C$ in Hall basis is built of $N_{A}$ letters $A$ and $N_{B}$ letters $B$. We show by induction on $N_{A}+N_{B}$, that

$${\rm ht}(C)\leq N_{B}-N_{A}+1.$$

(2.10)

The cases in which $N_{A}+N_{B}=1$ are obtained explicitly as ${\rm ht}(A)=0$ and ${\rm ht}(B)=2$, thus (2.10) is satisfied for the generators of the free Lie algebra. Now let us assume that a given non-zero commutator $C$ satisfies (2.10), so can be written as

$$C=\sum_{k=0}^{K}y_{k}(x)\partial_{x}^{k},$$

where $0\leq K\leq N_{B}-N_{A}+1$ and $y_{K}\not\equiv 0$. Then by (2.9),

	$\displaystyle[A,C]$	$\displaystyle=\varepsilon^{-1}\sum_{k=0}^{K-1}\left(\sum_{j=k}^{K}\binom{j}{k}y_{j}(x)V^{(j-k)}(x)\right)\partial_{x}^{k},$
	$\displaystyle[B,C]$	$\displaystyle=\varepsilon\sum_{k=0}^{K}y_{k}^{\prime\prime}(x)\partial_{x}^{k}+2y_{k}^{\prime}(x)\partial_{x}^{k+1},$

Therefore, ignoring the cases where these commutators vanish identically, we see that (2.10) is satisfied for $[A,C]$ and $[B,C]$ by the inductive hypothesis. This, in fact, completes the induction step for the entire Hall basis, because any commutator in the Hall basis can be written as a linear combination of words of the form

$$[a_{1},[a_{2},[\dots,[a_{n-1},a_{n}]\dots]]],$$

where $a_{k}\in\{A,B\}$ for all $k$, by the Jacobi identity (this is known as the Dynkin basis).

Next we show that every non-zero commutator $C$ in the Hall basis scales like ${\cal O}\!\left(\varepsilon^{-1}\right)$. Indeed, when $C$ is made up of $N_{A}$ letters $A=-\varepsilon^{-1}V$ and $N_{B}$ letters $B=\varepsilon\partial_{x}^{2}$, the linearity of commutators implies the equality of commutators $C$ and $\varepsilon^{N_{B}-N_{A}}\bar{C}$, where $\bar{C}$ has the same structure as $C$, but with $\bar{A}=-V$ and $\bar{B}=\partial_{x}^{2}$ instead of $A$ and $B$. Obviously ${\rm ht}(C)={\rm ht}(\bar{C})$. Now by Assumption 1, $\partial_{x}$ scales like $\varepsilon^{-1}$ after discretisation and by Assumption 2 we have that all variable coefficients $y_{k}$ lie in $\mathrm{C}^{\infty}_{\mathrm{loc}}(\mbox{\Bbb R})$ (so all derivatives are locally bounded). Therefore $\bar{C}=\mathcal{O}(\varepsilon^{-{\rm ht}(C)})$. Since for non-zero $C$, we have ${\rm ht}(C)\leq N_{B}-N_{A}+1$, we conclude that,

$$C=\varepsilon^{N_{B}-N_{A}}\bar{C}=\mathcal{O}(\varepsilon^{N_{B}-A_{N}-{\rm ht}(\bar{C})})=\mathcal{O}(\varepsilon^{N_{B}-N_{A}-(N_{B}-N_{A}+1)})=\mathcal{O}(\varepsilon^{-1}),$$

which concludes the proof of the theorem.        $\Box$

An immediate consequence of Theorem 3 and (2.8) is that

$${\mathrm{e}}^{\frac{1}{2}\tau A}{\mathrm{e}}^{\tau B}{\mathrm{e}}^{\frac{1}{2}\tau A}={\mathrm{e}}^{\tau(A+B)+{\cal O}\!\left(h^{3}\varepsilon^{-1}\right)}={\mathrm{e}}^{\tau(A+B)}+{\cal O}\!\left(h^{3}\varepsilon^{-1}\right).$$

This means that taking the time step size $h={\cal O}\!\left(\varepsilon\right)$ in the Strang splitting (2.7) yields a local truncation error of $\mathcal{O}(h^{2})$ or equivalently, $\mathcal{O}(\varepsilon^{2})$. However, a time step $h={\cal O}\!\left(\varepsilon\right)$ is overly expensive. If instead, one took a more reasonable $h={\cal O}\!\left(\varepsilon^{1/2}\right)$, then the local truncation error is effectively $\mathcal{O}(h)$ or equivalently, $\mathcal{O}(\varepsilon^{1/2})$. In summary, unless the time step is unacceptably reduced, the effective error of the Strang splitting is larger than that suggested by an analysis which ignores the smallness of $\varepsilon$.

This order reduction for the Strang splitting in the case of Hamiltonians in a semi-classical setting motivates the quest for higher order splittings. A systematic approach is to calculate higher order symmetric Zassenhaus splittings, first proposed in (?). Using this methodology we will derive two splittings for the solution operator $\exp(\tau(A+B))$ where $A=-\varepsilon^{-1}V$, $B=\varepsilon\partial_{x}^{2}$ and $\tau={\mathrm{i}}h$, of order ${\cal O}\!\left(h^{3}\varepsilon^{-1}\right)$ and ${\cal O}\!\left(h^{5}\varepsilon^{-1}\right)$ respectively, in the family of symmetric Zassenhaus splittings.

1. 1.

   To derive the first symmetric Zassenhaus splitting, we apply the sBCH formula in the following way.

   $${\mathrm{e}}^{-\frac{1}{2}\tau A}{\mathrm{e}}^{\tau A+\tau B}{\mathrm{e}}^{-\frac{1}{2}\tau A}={\mathrm{e}}^{\mathrm{sBCH}(-\tau A,\tau A+\tau B)}$$

   (2.11)

   where

   		$\displaystyle\mathrm{sBCH}(-\tau A,\tau A+\tau B)=$		(2.12)
   		$\displaystyle=\tau B+\tau^{3}\frac{1}{24}[[B,A]A]+\tau^{3}\frac{1}{12}[[B,A],B]$
   		$\displaystyle-\tau^{5}\frac{1}{720}[[[[B,A],A],B],B]-\tau^{5}\frac{1}{720}[[[[B,A],B],B],B]$
   		$\displaystyle-\tau^{5}\frac{1}{480}[[[B,A],A],[B,A]]-\tau^{5}\frac{1}{240}[[[B,A],B],[B,A]]+{\cal O}\!\left(h^{7}\varepsilon^{-1}\right).$

   Substituting (2.12) into (2.11) results in the first symmetric Zassenhauss splitting, which coincides with Strang splitting,

   	$\displaystyle{\mathrm{e}}^{\tau(A+B)}$	$\displaystyle={\mathrm{e}}^{\frac{1}{2}\tau A}{\mathrm{e}}^{\mathrm{sBCH}(-\tau A,\tau A+\tau B)}{\mathrm{e}}^{\frac{1}{2}\tau A}$		(2.13)
   		$\displaystyle={\mathrm{e}}^{\frac{1}{2}\tau A}{\mathrm{e}}^{\tau B}{\mathrm{e}}^{\frac{1}{2}\tau A}+{\cal O}\!\left(h^{3}\varepsilon^{-1}\right).$

2. 2.

   To derive the second symmetric Zassenhaus splitting, we split the inner term of (2.13) by the same approach as above, that is

   $\displaystyle{\mathrm{e}}^{-\frac{1}{2}\tau B}{\mathrm{e}}^{\mathrm{sBCH}(-\tau A,\tau A+\tau B)}{\mathrm{e}}^{-\frac{1}{2}\tau B}$

   $\displaystyle={\mathrm{e}}^{\mathrm{sBCH}(-\tau B,\mathrm{sBCH}(-\tau A,\tau A+\tau B))}$

   which leads to,

   $${\mathrm{e}}^{\tau A+\tau B}={\mathrm{e}}^{\frac{1}{2}\tau A}e^{\frac{1}{2}\tau B}{\mathrm{e}}^{\mathrm{sBCH}(-\tau B,\mathrm{sBCH}(-\tau A,\tau A+\tau B))}{\mathrm{e}}^{\frac{1}{2}\tau B}e^{\frac{1}{2}\tau A},$$

   (2.14)

   where

   	$\displaystyle\mathrm{sBCH}(-\tau B,\mathrm{sBCH}(-\tau A,\tau B+\tau A))$
   	$\displaystyle=\frac{1}{24}\tau^{3}[[B,A],A]+\frac{1}{12}\tau^{3}[[B,A],B]$
   	$\displaystyle-\frac{19}{2880}\tau^{5}[[[[B,A],A],B],B]-\frac{17}{1440}\tau^{5}[[[[B,A],B],B],B]$
   	$\displaystyle-\tau^{5}\frac{1}{480}[[[B,A],A],[B,A]]-\tau^{5}\frac{1}{240}[[[B,A],B],[B,A]]+{\cal O}\!\left(h^{7}\varepsilon^{-1}\right).$

   Observe that by Theorem 3, the first two commutators (which involve three letters) scale like ${\cal O}\!\left(h^{3}\varepsilon^{-1}\right)$ and the remainder scales like ${\cal O}\!\left(h^{5}\varepsilon^{-1}\right)$. Therefore, these first two terms are what will appear in this Zassenhaus splitting. However, the commutator,

   $\displaystyle[[B,A],B]=[[\varepsilon\partial_{x}^{2},-\varepsilon^{-1}V],\varepsilon\partial_{x}^{2}]$

   $\displaystyle=\varepsilon\left(V^{(4)}+4V^{(3)}\partial_{x}+4V^{(2)}\partial_{x}^{2}\right),$

   will not be skew-Hermitian after discretisation (which would result in loss of unitarity of the method), and therefore cannot be substituted into (2.14). For this reason, as proposed in (?), we use a substitution rule of the following kind:

   $$y(x)\partial_{x}=-\frac{1}{2}\left[\int_{x_{0}}^{x}y(s){\rm d}s\right]\partial_{x}^{2}-\frac{1}{2}\partial_{x}y(x)+\frac{1}{2}\partial_{x}^{2}\left[\int_{x_{0}}^{x}y(s){\rm d}s\,\cdot\right],$$

   and obtain terms that remain skew-Hermitian after discretisation:

   		$\displaystyle\mathrm{sBCH}(-\tau B,\mathrm{sBCH}(-\tau A,\tau B+\tau A))$
   		$\displaystyle=\tau^{3}\varepsilon^{-1}\frac{1}{12}(V^{(1)})^{2}+\tau^{3}\varepsilon\frac{1}{12}V^{(4)}+\tau^{3}\varepsilon\frac{1}{3}V^{(3)}\partial_{x}+\tau^{3}\varepsilon\frac{1}{3}V^{(2)}\partial_{x}^{2}+{\cal O}\!\left(h^{5}\varepsilon^{-1}\right)$
   	$\displaystyle=$	$\displaystyle\tau^{3}\varepsilon^{-1}\frac{1}{12}(V^{(1)})^{2}+\frac{1}{6}\tau^{3}\varepsilon\underbrace{\left\{V^{(2)}\partial_{x}^{2}+\partial_{x}^{2}\left[V^{(2)}\cdot\right]\right\}}_{{\cal O}\!\left(\varepsilon^{-2}\right)}-\frac{1}{12}\tau^{3}\varepsilon V^{(4)}+{\cal O}\!\left(h^{5}\varepsilon^{-1}\right).$

   In the final form of the splitting (2.14) the small ${\cal O}\!\left(h^{3}\varepsilon\right)$ term involving $V^{(4)}$ can be discarded.

Summing up these two derivations, we have the splittings,

$$u^{k+1}(x)={\mathrm{e}}^{\mathcal{R}_{0}}{\mathrm{e}}^{2\mathcal{R}_{1}}{\mathrm{e}}^{\mathcal{R}_{0}}u^{k}(x)+{\cal O}\!\left(h^{3}\varepsilon^{-1}\right)$$

(2.15)

and

$$u^{k+1}(x)={\mathrm{e}}^{\mathcal{R}_{0}}{\mathrm{e}}^{\mathcal{R}_{1}}{\mathrm{e}}^{2\mathcal{R}_{2}}{\mathrm{e}}^{\mathcal{R}_{1}}{\mathrm{e}}^{\mathcal{R}_{0}}u^{k}(x)+{\cal O}\!\left(h^{5}\varepsilon^{-1}\right),$$

(2.16)

where, letting $\tau={\mathrm{i}}h$,

	$\displaystyle\mathcal{R}_{0}$	$\displaystyle=$	$\displaystyle-\frac{1}{2}\tau\varepsilon^{-1}V,$
	$\displaystyle\mathcal{R}_{1}$	$\displaystyle=$	$\displaystyle\frac{1}{2}\tau\varepsilon\partial_{x}^{2},$
	$\displaystyle\mathcal{R}_{2}$	$\displaystyle=$	$\displaystyle\color[rgb]{0,0,0}\frac{1}{12}\tau^{3}\varepsilon\left\{\partial_{x}^{2}[V^{(2)}\,\cdot\,]+V^{(2)}\partial_{x}^{2}\right\}+\frac{1}{24}\tau^{3}\varepsilon^{-1}(V^{(1)})^{2}.$

Note that $\mathcal{R}_{0}={\cal O}\!\left(h\varepsilon^{-1}\right)$, $\mathcal{R}_{1}={\cal O}\!\left(h\varepsilon^{-1}\right)$, $\mathcal{R}_{2}={\cal O}\!\left(h^{3}\varepsilon^{-1}\right)$.

It is also possible to derive even higher order methods, such as

$$u^{n+1}(x)={\mathrm{e}}^{\mathcal{R}_{0}}{\mathrm{e}}^{\mathcal{R}_{1}}{\mathrm{e}}^{\mathcal{R}_{2}}{\mathrm{e}}^{2\mathcal{R}_{3}}{\mathrm{e}}^{\mathcal{R}_{2}}{\mathrm{e}}^{\mathcal{R}_{1}}{\mathrm{e}}^{\mathcal{R}_{0}}u^{n}(x)+{\cal O}\!\left(h^{7}\varepsilon^{-1}\right),$$

(2.17)

where

	$\displaystyle\mathcal{R}_{3}$	$\displaystyle=$	$\displaystyle-\frac{1}{120}\tau^{5}\varepsilon^{-1}V^{(2)}(V^{(1)})^{2}+\frac{1}{24}\tau^{3}\varepsilon V^{(4)}$
			$\displaystyle+\frac{1}{120}\tau^{5}\varepsilon\left\{\partial_{x}^{2}\left[\left(7(V^{(2)})^{2}+V^{(3)}V^{(1)}\right)\,\cdot\right]+\left(7(V^{(2)})^{2}+V^{(3)}V^{(1)}\right)\partial_{x}^{2}\right\}$
			$\displaystyle+\frac{1}{60}\tau^{5}\varepsilon^{-3}\left\{\partial_{x}^{4}\left[V^{(4)}\,\cdot\,\right]+V^{(4)}\partial_{x}^{4}\right\}.$

Note that $\mathcal{R}_{3}={\cal O}\!\left(h^{5}\varepsilon^{-1}\right)$. We refer the reader to (?) for discussion of deriving such higher order methods via a sequence of skew-Hermitian operators $\mathcal{R}_{0},\mathcal{R}_{1},\ldots$. Our new analysis encapsulated in Theorem 3 shows that each term $\mathcal{R}_{\ell}$ is actually of size ${\cal O}\!\left(h^{2\ell-1}\varepsilon^{-1}\right)$ for $\ell=1,2,\ldots$. In Section 5, we will discuss how to go about computing ${\mathrm{e}}^{\mathcal{R}_{\ell}}$ for each $\ell$.

Solving (1.1) by spectral methods based upon symmetric Zassenhaus splittings (2.15) or (2.16) involves three ingredients: the splitting itself into $\mathcal{R}_{0},\mathcal{R}_{1},\mathcal{R}_{2},\ldots$, the choice of spectral basis, and the means to compute the exponentials ${\mathrm{e}}^{\mathcal{R}_{\ell}}$. The generalisation of each to the new setting requires new ideas and substantial effort. In this subsection we are concerned with the choice of the spectral basis.

We seek a set $\Phi=\{\varphi_{n}\}_{n=0}^{\infty}$ which forms an orthonormal basis of $\mathrm{L}_{2}(\mbox{\Bbb R})$ – this means that any $f\in\mathrm{L}_{2}(\mbox{\Bbb R})$ can be expanded in the form

$$f(x)=\sum_{n=0}^{\infty}\hat{f}_{n}\varphi_{n}(x),\qquad\mbox{where}\qquad\hat{f}_{n}=\int_{-\infty}^{\infty}f(x)\overline{\varphi_{n}(x)}\,\mathrm{d}x,\quad n\in\mbox{\Bbb Z}_{+}.$$

For the time being we require the $\varphi_{n}$s to be real, although this will be lifted as necessary (with suitable changes). In addition we require that $\Phi$ has a tridiagonal differentiation matrix (which, it is easy to prove, must be skew-symmetric),

$$\varphi_{n}^{\prime}=-b_{n-1}\varphi_{n-1}+b_{n}\varphi_{n+1},\qquad n\in\mbox{\Bbb Z}_{+},$$

(3.1)

where $b_{-1}=0$ and $b_{n}>0$, $n\in\mbox{\Bbb Z}_{+}$. This makes both computation and analysis considerably easier.

A comprehensive theory of such orthogonal systems has been developed in (?, ?). The main issue, making (3.1) compatible with orthonormality, can be explicated by considering Fourier transforms of the $\varphi_{n}$s. Specifically, let $w(\xi)\,\mathrm{d}\xi$ be a Borel measure over ℝ and its Radon–Nikodym derivative $w\geq 0$ be absolutely continuous and even. Furthermore assume that all the moments of this measure are finite. Such measure generates a system of orthonormal polynomials $\{p_{n}\}_{n=0}^{\infty}$,

$$\int_{-\infty}^{\infty}p_{n}(\xi)p_{m}(\xi)w(\xi)\,\mathrm{d}\xi=0,\quad m\neq n,\qquad\int_{-\infty}^{\infty}p_{n}^{2}(\xi)w(\xi)\,\mathrm{d}\xi=1.$$

Then the scaled inverse Fourier transform,

$$\varphi_{n}(x)=\frac{(-{\mathrm{i}})^{n}}{\sqrt{2\pi}}\int_{-\infty}^{\infty}p_{n}(\xi)g(\xi){\mathrm{e}}^{{\mathrm{i}}x\xi}\,\mathrm{d}\xi,\qquad n\in\mbox{\Bbb Z}_{+},$$

(3.2)

where $g$ is any function satisfying $|g(\xi)|^{2}=w(\xi)$, forms an orthonormal system on the real line which satisfies (3.1). Note that this system is real-valued if and only if $g$ has even real part and odd imaginary part, for example $g(\xi)=\sqrt{w(\xi)}$. The constants $b_{n}$ in (3.1) are inherited from the recurrence relation for orthonormal polynomials,

$$b_{n}p_{n+1}(\xi)=\xi p_{n}(\xi)-b_{n-1}p_{n-1}(\xi),\qquad n\in\mbox{\Bbb Z}_{+}.$$

The orthonormal system given by (3.2) need not be dense in ℝ – as a matter of fact, it is dense in the Paley–Wiener space $\mathcal{PW}_{\mbox{supp}(w)}(\mbox{\Bbb R})\subseteq\mathrm{L}_{2}(\mbox{\Bbb R})$ which is the space of $L_{2}(\mbox{\Bbb R})$ functions whose Fourier transforms vanish outside of the support of $w$. Therefore, the system is a basis of $\mathrm{L}_{2}(\mbox{\Bbb R})$ if and only if the weight function $w$ is positive on the whole real line.

Complete orthonormal bases can be formed also from polynomials $P=\{p_{n}\}_{n=0}^{\infty}$ orthogonal on the half-line $[0,\infty)$ (?), e.g. the Laguerre polynomials whose orthogonality measure is ${\mathrm{e}}^{-\xi}\,\mathrm{d}\xi$, $\xi\geq 0$: The representation (3.2) survives intact but, to render the system dense in $\mathrm{L}_{2}(\mbox{\Bbb R})$, we need to complement $P$ with orthogonal polynomials with respect to the mirror image of $w$ in the left half-line, $w(-\xi)\,\mathrm{d}\xi$ for $\xi\leq 0$. The new system $\Phi$ is enumerated by $n\in\mbox{\Bbb Z}$ and in place of (3.1) we have

$$\varphi_{n}^{\prime}=-b_{n-1}\varphi_{n-1}+{\mathrm{i}}c_{n}\varphi_{n}+b_{n}\varphi_{n+1},\qquad n\in\mbox{\Bbb Z},$$

with $b_{n}>0$, $n\neq 0$, $b_{0}=0$ and real $c_{n}$ – note that the new differentiation matrix is skew-Hermitian.

Given an orthonormal system $\Phi$ on the real line, we denote by $\psi_{n}$, $n\in\mbox{\Bbb Z}_{+}$, the solution of the free Schrödinger equation (1.3) with the initial condition $\varphi_{n}$ – in other words,

$$\frac{\partial\psi_{n}}{\partial t}=-{\mathrm{i}}\varepsilon\frac{\partial^{2}\psi_{n}}{\partial x^{2}},\qquad\psi_{n}(x,0)=\varphi_{n}(x),\;x\in\mbox{\Bbb R}.$$

(3.3)

We call $\Psi(t)=\{\psi_{n}(\cdot,t)\}_{n=0}^{\infty}$ the free Schrödinger evolution (FSE) of $\Phi$.

The exact solution of (3.3) via the Fourier transform is well known and can be easily verified by direct differentiation:

$$\psi_{n}(x,t)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}\hat{\varphi}_{n}(\eta){\mathrm{e}}^{{\mathrm{i}}\eta^{2}\varepsilon t+{\mathrm{i}}\eta x}\,\mathrm{d}\eta,$$

(3.4)

where

$$\hat{\varphi}_{n}(\eta)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}\varphi_{n}(\xi){\mathrm{e}}^{-{\mathrm{i}}\eta\xi}\,\mathrm{d}\xi$$

is the familiar Fourier transform of $\varphi_{n}$.

On the face of it, our job is done: any mention of the phrase “Fourier transform” elicits from a numerical analyst the instinctive response “Fast Fourier Transform!”. This, however, is somewhat rash. An FFT computes rapidly the discrete Fourier transform which, in turn, is a very precise (at any rate, for very smooth functions) approximation of the Fourier transform of a periodic function in a compact interval, while our setting is the entire real line. One possibility is to clip the real line, approximating it by a sufficiently large interval and disregarding the Gibbs effect at the endpoints. This immediately begs the question “how large” which, while not beyond the ken of numerical reasoning, presents its own challenges. In this paper we adopt an alternative – and arguably more effective – point of view, seeking the exact solution of (3.4) for specific orthonormal systems $\Phi$. While this approach cannot be expected to apply to each and every $\Phi$ consistent with the setting of Subsection 2.1, it does so with the two most interesting orthonormal systems: Hermite functions and Malmquist–Takenaka functions.

Once FSEs $\Psi(t)$ are known, the solution of the free Schrödinger equation (1.3) with the initial condition $u(x,kh)$ proceeds as follows: The function $u(x,kh)$ is expanded in the orthonormal basis $\Phi$,

$$u(x,kh)\approx\sum_{n=0}^{N}\hat{u}_{n}\varphi_{n}(x)$$

(3.5)

for a sufficiently large truncation parameter $N$. Having done so, linearity of (1.3) implies that

$$u(x,(k+1)h)\approx\sum_{n=0}^{N}\hat{u}_{n}\psi_{n}(x,h).$$

(3.6)

We get the coefficients for free because they do not change — it is the basis which changes. The choice of $N$ is governed by approximation properties of the spectral basis, and its ability to approximate spatial oscillations of frequency ${\cal O}\!\left(\varepsilon^{-1}\right)$ as discussed in the introduction.

Indeed, orthonormal systems are not all of equal value: more specifically, they can approximate functions at different speeds. While standard spectral methods on a torus are known to converge (for analytic functions) at an exponential speed, equivalent theory does not exist yet on the real line. Recalling from Section 1 that solutions of (1.1) are typically composed of wave packets, it is instructive to enquire how well different orthonormal systems approximate wave packets. This is investigated in (?) for the two families $\Phi$ described in the sequel: in both cases we can prove exponential convergence to any set error tolerance.

We note for further reference that the computation of (3.6) (once $N$ and $h$ have been appropriately chosen) requires both the knowledge of $\Psi(h)$ and the means to evaluate an expansion as in (3.5).