Fourth order real space solver for the time-dependent Schrödinger equation with singular Coulomb potential

Let us consider the axially symmetric three-dimensional time-dependent Schrödinger equation in cylindrical coordinates $\rho=\sqrt{x^{2}+y^{2}}$ and $z$:

$$i\hbar\frac{\partial}{\partial t}\Psi\left(z,\rho,t\right)=-\frac{\hbar^{2}}{2\mu}\left[\frac{\partial^{2}}{\partial z^{2}}+\frac{\partial^{2}}{\partial\rho^{2}}+\frac{1}{\rho}\frac{\partial}{\partial\rho}\right]\Psi\left(z,\rho,t\right)+V\left(z,\rho,t\right)\Psi\left(z,\rho,t\right)$$

(1)

where $\mu$ is the (reduced) electron mass. For the sake of simplicity, we use atomic units throughout this article and we use the notation $\beta=-\hbar^{2}/2\mu$.

As discussed in the Introduction, such a Hamiltonian plays a fundamental role e.g. in the interaction of atomic systems with strong laser pulses. Then the potential is usually the sum of an atomic binding potential (centered in the origin) and an interaction energy term. For a single active electron and a linearly polarized laser pulse, this interaction term in dipole approximation and in length gauge [47] reads as:

$$V\left(z,\rho,t\right)=V_{\text{atom}}\left(z,\rho\right)-q\cdot F(t)\cdot z,$$

(2)

where $q=-1$ denotes the electron’s electric charge and $F(t)$ is the laser’s electric field, pointing in the $z$ direction.

The wave function is defined within the interval $z\in\left[z_{\min},z_{\max}\right]$, $\rho\in\left[0,\rho_{\max}\right]$ and the problem is treated as an initial value problem as $\Psi\left(z,\rho,0\right)=\psi\left(z,\rho\right)$ is known. Box boundary conditions are posed at three edges of the interval:

$$\Psi\left(z>z_{\max},\rho,t\right)=\Psi\left(z<z_{\min},\rho,t\right)=\Psi\left(z,\left|\rho\right|>\rho_{\text{max}},t\right)=0$$

(3)

and because of the singularity in the curvilinear coordinate $\rho$, Neumann boundary condition should be posed along the line $\rho=0$. In order to find this boundary condition, we multiply (1) by $\rho$ and take the limit $\rho\rightarrow 0$:

$$\left.\frac{\partial\Psi}{\partial\rho}\right|_{\rho=0}=-\frac{1}{\beta}\lim_{\rho\rightarrow 0}\left(\rho V\left(z,\rho,t\right)\Psi\left(z,\rho,t\right)\right)=0$$

(4)

where we have assumed the potential is smooth at $\rho=0$ line. Then $\Psi(z,\rho,t)$ is continuous and continuously differentiable everywhere in the $z-\rho$ plane, and it has extremal value in $\rho=0$. Symmetry condition $\Psi(z,\rho,t)=\Psi(z,-\rho,t)$ also holds in this case.

However, this Neumann boundary condition has to be altered when we introduce a singular potential in any form, e.g. the 3D Coulomb potential $V_{{\rm atom}}\left(z,\rho\right)=\gamma/r$ with $r=\sqrt{z^{2}+\rho^{2}}$ and $\gamma=-q^{2}Z/4\pi\epsilon_{0}$.

Our aim is to represent the singularity of the 3D Coulomb potential by implementing the correct boundary condition in the origin. In order to find this condition, it seems necessary to use spherical polar coordinates, since the Coulomb problem is not separable in cylindrical coordinates. Let us first discuss the properties of the radial solutions $\phi(r)$ of the well known radial equation [48, 49]:

$$\left[\beta\frac{\partial^{2}}{\partial r^{2}}+\frac{2\beta}{r}\frac{\partial}{\partial r}+\beta\frac{l(l+1)}{r^{2}}+\frac{\gamma}{r}\right]\phi(r)=E\phi(r),$$

(5)

where $l$ is the angular momentum quantum number. If we multiply both sides with $r$ then take the $r\rightarrow 0$ limit, we can acquire Robin boundary condition for all of the $l=0$ states as

$$\left.\frac{\partial\phi}{\partial r}\right|_{r=0}=-\frac{\gamma}{2\beta}\phi\left(0\right),$$

(6)

but the bound states with higher angular momenta will have $\phi\left(0\right)=0$ boundary condition at the origin, in accordance with the fact that for $l>0$ the particle physically cannot penetrate to $r=0$, because of the singularity of $l(l+1)/r^{2}$.

Since we need the boundary conditions for the cylindrical equation independent of the $l$ values, we turn to the nonseparated symmetric spherical Schrödinger equation for answers:

$$\left[\beta\frac{\partial^{2}}{\partial r^{2}}+\frac{2\beta}{r}\frac{\partial}{\partial r}+\beta\frac{1}{r^{2}\sin\theta}\frac{\partial}{\partial\theta}\left(\sin\theta\frac{\partial}{\partial\theta}\right)+\frac{\gamma}{r}\right]\Psi(r,\theta)=E\Psi(r,\theta)$$

(7)

which has also the term $1/r^{2}$ and has singularity at $r=0$ and at $\theta=k\pi$, $k=0,1,2\dots$. To acquire the boundary conditions for $r=0$, we again multiply both sides with $r$ and take the $r\rightarrow 0$ limit, yielding

$$\lim_{r\rightarrow 0}\left[2\beta\frac{\partial}{\partial r}+\beta\frac{\cos\theta}{r\sin\theta}\frac{\partial}{\partial\theta}+\frac{\beta}{r}\frac{\partial^{2}}{\partial^{2}\theta}+\gamma\right]\Psi(r,\theta)=0.$$

(8)

Now we transform equation (8) into cylindrical coordinates $z=r\cos\theta$ and $\rho=r\sin\theta$, using following expressions of the partial derivatives:

$$\frac{\partial}{\partial r}=\frac{\partial z}{\partial r}\frac{\partial}{\partial z}+\frac{\partial\rho}{\partial r}\frac{\partial}{\partial\rho}=\sin\theta\frac{\partial}{\partial\rho}+\cos\theta\frac{\partial}{\partial z},$$

(9)

$$\frac{\partial}{\partial\theta}=\frac{\partial z}{\partial\theta}\frac{\partial}{\partial z}+\frac{\partial\rho}{\partial\theta}\frac{\partial}{\partial\rho}=r\cos\theta\frac{\partial}{\partial\rho}-r\sin\theta\frac{\partial}{\partial z}$$

(10)

with $\cos\theta=z/\sqrt{z^{2}+\rho^{2}}$ and $\sin\theta=\rho/\sqrt{z^{2}+\rho^{2}}$. After writing back and taking the limit, the formula (8) becomes

$$\left.\left[2\beta\frac{\rho}{\sqrt{z^{2}+\rho^{2}}}\frac{\partial}{\partial\rho}+\beta\frac{z}{\sqrt{z^{2}+\rho^{2}}}\frac{\partial}{\partial z}+\beta\frac{z^{2}}{\rho\sqrt{z^{2}+\rho^{2}}}\frac{\partial}{\partial\rho}+\gamma\right]\Psi(r,\theta)\right|_{r=0}=0.$$

(11)

Then, by substituting $z=0$ into (11), we obtain the Robin boundary condition for the 3D Coulomb singularity:

$$\left.\frac{\partial\Psi}{\partial\rho}\right|_{\rho,z=0}=-\frac{\gamma}{2\beta}\Psi\left(0,0,t\right).$$

(12)

This can be generalized to include multiple Coulomb-cores rather easily, we just need to impose (12) at multiple $z,\rho$ points along the $\rho=0$ axis. Additionally, any continuously differentiable potential added to this configuration does not change this boundary condition. It is interesting to note also that the boundary condition imposed by a 1D Dirac-delta potential $\gamma\delta(\rho)$ also has the form of (12) for a symmetric wave function [49].

Although we are mainly interested in the solution of the TDSE with an axially symmetric initial state, i.e. an initial state of zero magnetic quantum number $m$, we briefly show in the following that initial states with nonzero $m$ also lead to a TDSE of the form (13), however with different boundary conditions at $\rho=0$.

Let us write out the full 3D cylindrical time-dependent Schrödinger equation, again with the same axially symmetric potential that we used previously:

$$i\frac{\partial}{\partial t}\Xi\left(z,\rho,\phi,t\right)=\left[\beta\frac{\partial^{2}}{\partial z^{2}}+\beta\frac{\partial^{2}}{\partial\rho^{2}}+\frac{\beta}{\rho}\frac{\partial}{\partial\rho}+\frac{\beta}{\rho^{2}}\frac{\partial^{2}}{\partial\phi^{2}}+V\left(z,\rho,t\right)\right]\Xi\left(z,\rho,\phi,t\right).$$

(13)

We can write the quantum state with a given magnetic quantum number $m$ in the following form for any $t$

$$\Xi\left(z,\rho,\phi,t\right)=\Psi_{m}(z,\rho,t)e^{im\phi},$$

(14)

due to the fact that the Hamiltonian of (13) conserves the $m$ quantum number. (We can also use real angular basis functions like $\sin(m\phi)$ or $\cos(m\phi)$ instead of $e^{im\phi}$, as long as that they are the eigenfunctions of $\partial_{\phi}^{2}$.)

After substituting (14) to (13) we perform projection in the form of $\intop_{0}^{2\pi}e^{-im\phi}(\cdot){\rm d}\phi$ . Then we get the TDSE for our wave function $\Psi_{m}(z,\rho,t)$ with the specified magnetic quantum number $m$:

$$i\frac{\partial}{\partial t}\Psi_{m}(z,\rho,t)=\left[\beta\frac{\partial^{2}}{\partial z^{2}}+\beta\frac{\partial^{2}}{\partial\rho^{2}}+\frac{\beta}{\rho}\frac{\partial}{\partial\rho}-\beta\frac{m^{2}}{\rho^{2}}+V\left(z,\rho,t\right)\right]\Psi_{m}(z,\rho,t),$$

(15)

Eq. (15) has the same form as the axially symmetric equation (1) but with a new $m$ dependent potential:

$$V_{m}(z,\rho,t)=V\left(z,\rho,t\right)-\beta m^{2}/\rho^{2}.$$

(16)

However, if $m\neq 0$ then the $1/\rho^{2}$ singularity of (16) alters the boundary condition (4) at the $\rho=0$ axis to the following:

$$\Psi_{m}(z,0,t)=0,$$

(17)

which can by derived by multiplying (15) with $\rho^{2}$ and taking the limit of $\rho\rightarrow 0$. The boundary condition (17) also works with Coulomb potential, and in every case where $\rho^{2}V\left(z,\rho,t\right)=0$ in the limit of $\rho\rightarrow 0$. It is also consistent with the boundary conditions for Coulomb states with nonzero angular momenta ($l\neq 0$) mentioned in the previous section.

The formal solution of (1) in terms of the time evolution operator $U(t,t^{\prime})=\mathcal{T}\exp\left[-\frac{i}{\hbar}\int H(t^{\prime\prime})dt^{\prime\prime}\right]$ is the following

$$|\Psi(t)\rangle=U(t,t^{\prime})|\Psi(t^{\prime})\rangle.$$

(18)

Here, the exponential operator is to be understood as a time-ordered quantity, which is a difficult procedure if the Hamilton operators at different time instants do not commute. However, this is the case for the potential we are considering.

To acquire suitable discretization in the time domain of problem (1) we approximate the $U(t,0)$ time evolution operator directly. First, let us divide the time domain $[0,t]$, into $N_{t}$ equal subintervals, then by the group property of $U(t,0)$ we get

$$U\left(t,0\right)=\prod_{k=0}^{N_{t}-1}U\left(t_{k+1},t_{k}\right)$$

(19)

where $t_{k}=k\Delta t$ and $\Delta t=t/N_{t}$. In one interval, we can write the evolution operator with its short time form of

$$U\left(t_{k+1},t_{k}\right)=e^{-i\Delta t\,H_{k}}$$

(20)

where $H_{k}$ is the effective time-independent Hamiltonian related to the original one $H(t)$ by the Magnus expansion [50, 51]:

$$H_{k}=\frac{1}{\Delta t}\intop_{t_{k}}^{t_{k+1}}H(t^{\prime})\text{d}t^{\prime}+\frac{i}{2\Delta t}\intop_{t_{k}}^{t_{k+1}}\intop_{t_{k}}^{t^{{}^{\prime}}}\left[H(t^{\prime\prime}),H(t^{\prime})\right]\text{d}t^{\prime\prime}\text{d}t^{\prime}+\dots.$$

(21)

This includes infinitely many commutators of the Hamiltonians evaluated at different time points, to be integrated with respect to more and more variables.

Using only the first term of this series, one can directly acquire the well known second order approximation of the Hamiltonian as

$$H_{k}^{(2)}=H(t_{k+\frac{1}{2}}).$$

(22)

In order to have information about the next error term of the time evolution, we need to evaluate the Magnus commutator series to fourth order. According to Puzynin et. al. [50], this improved approximation for a TDSE of the form (1) is

$$H_{k}^{(4)}=\frac{1}{2\mu}\left(-i\nabla+\frac{\Delta t^{2}}{12}\nabla\dot{V}(t_{k+1/2})\right)^{2}+V(t_{k+1/2})+\frac{\Delta t^{2}}{24}\ddot{V}(t_{k+1/2}).$$

(23)

where the top dots are the short hand notation for $\partial_{t}$ time derivatives. This formula is important even if we use only the second order approximation, because the leading order of error depends on characteristics of the first and second time derivatives of $V(z,\rho,t)$.

For the exponential operators (20), first we consider the diagonal Padé-approximation of an exponential function

$$e^{\lambda\cdot x}=\frac{{}_{1}F_{1}(-M,-2M,\lambda\cdot x)}{{}_{1}F_{1}(-M,-2M,-\lambda\cdot x)}+O(\lambda^{2M+1})$$

(24)

where ${}_{1}F_{1}$ is the confluent hypergeometric function, which in this case reduces to a polynomial of degree $M$ with real coefficients. This expression can be used for the exponential operators (20) with $\lambda=-i\Delta t$, and it can be shown that for self-adjoint operators the approximation is unitary [23].

From this, a generalized operator approximation scheme that is $\Delta t^{2M+1}$ accurate can be obtained [44, 52] in a general form of

$$e^{-i\Delta t\,H_{k}}=\prod_{s=1}^{M}\left[\left(1+i\frac{\Delta t}{x_{s}^{*}}H_{k}\right)^{-1}\left(1-i\frac{\Delta t}{x_{s}}H_{k}\right)\right]+O(\Delta t^{2M+1})$$

(25)

where $x_{s}$ for $s=1,\dots,M$ are the roots of the polynomial equation

$${}_{1}F_{1}\left(-M,-2M,-x\right)=0.$$

(26)

If we truncate both the Magnus-series at the first term in (21) and take a single coefficient of the Padé-approximation ($M=1$), then we arrive at the second order accurate implicit Crank-Nicolson scheme:

$$\Psi(t_{k+1})=\left(1+iH_{k}\frac{\Delta t}{2}\right)^{-1}\left(1-iH_{k}\frac{\Delta t}{2}\right)\Psi(t_{k}).$$

(27)

One can straightforwardly construct higher order Crank-Nicolson schemes using (25), yielding multiple implicit Crank-Nicolson substeps [44]. For time dependent Hamiltonians though, a high order accurate scheme should use the corresponding high order effective Hamiltonian (21) for consistency.

Our next step is to discretize the effective Hamiltonian $H_{k}$, and to construct the Hamiltonian matrix. We propose to use the method of fourth order finite differences with $z,\rho$ cylindrical coordinates and the following equidistant 2D spatial grid:

$$z_{i}=z_{\min}+i\cdot\Delta z,\,\Delta z=(z_{\max}-z_{\min})/N_{z},\,i\in\left[0,N_{z}\right],$$

(28)

$$\rho_{j}=j\cdot\Delta\rho,\,\Delta\rho=\rho_{\max}/N_{\rho},\,j\in\left[0,N_{\rho}\right].$$

(29)

The first term of the Hamiltonian $H_{k}$ is the Laplacian $\nabla^{2}$. We denote the finite difference approximation of its $z$- and $\rho$-dependent terms by $L_{z}$ and $L_{\rho}$, respectively, and we use the following fourth order accurate finite difference forms [52] for them:

$$L_{z}\Psi_{i,j}(t)=\frac{-\Psi_{i-2,j}+16\Psi_{i-1,j}-30\Psi_{i,j}+16\Psi_{i+1,j}-\Psi_{i+2,j}}{12\Delta z^{2}},$$

(30)

$$L_{\rho}\Psi_{i,j}(t)=\frac{(-1+1/j)\Psi_{i,j-2}+(16-8/j)\Psi_{i,j-1}-30\Psi_{i,j}+(16+8/j)\Psi_{i,j+1}+(-1-1/j)\Psi_{i,j+2}}{12\Delta\rho^{2}}.$$

(31)

These fourth order formulae are optimal in the sense that symmetric finite differences of more than five points are very complicated to be applied for the Coulomb-problem, because the higher order the finite difference formula is, the higher derivatives of $\Psi(z,\rho,t)$ should be continuous.

To discretize the Neumann and Robin boundary conditions, we need a one-sided finite difference formula for the first derivative. Based on the method of Ref. [52], we derived the following fourth order accurate forward difference operator:

$$D_{\rho}\Psi_{i,j}(t)=\frac{-25\Psi_{i,j}+48\Psi_{i,j+1}-36\Psi_{i,j+2}+16\Psi_{i,j+3}-3\Psi_{i,j+4}}{12\Delta\rho}.$$

(32)

Using standard second order accurate form (27) of the exponential operator with the discretized Laplacian in the Hamilton matrix, we get the following implicit scheme for all $i\in\left[0,N_{z}\right]$ and $j\in\left[1,N_{\rho}\right]$:

$$\left(1+\alpha\beta L_{z}+\alpha\beta L_{\rho}+\alpha V_{i,j}\right)\Psi_{i,j}(t_{k+1})=\left(1-\alpha\beta L_{z}-\alpha\beta L_{\rho}-\alpha V_{i,j}\right)\Psi_{i,j}(t_{k}),$$

(33)

where $\alpha=i\Delta t/2$ and the potential is evaluated at the temporal midpoint $V_{i,j}=V(z_{i},\rho_{j},t_{k+1/2})$. The box boundary conditions $\Psi_{i,j}=0$ is applied if $j>N_{\rho}$ or $i>N_{z}$ or $i<0$.

Assuming that the potential is smooth, the Neumann boundary condition (4) prescribes the following implicit equations at $j=0$ for all $i\in[0,N_{z}]$:

$$D_{\rho}\Psi_{i,0}(t_{k+1})=0.$$

(34)

On the other hand, if a Coulomb-core of strength $\gamma$ is present at the grid point $z_{R}=0$ and $m=0$ then, according to (12), equation (34) will be overridden at $i=R$ by

$$\left(D_{\rho}+\frac{\gamma}{2\beta}\right)\Psi_{R,0}(t_{k+1})=0.$$

(35)

Here we have assumed that the origin is included in (28) with $z_{R}=0$, where $R$ can be anywhere in the interval $[0,N_{z}]$ if it is reasonably far from the box boundaries.

For $m\neq 0$ states, the equations of the boundary conditions at $\rho=0$ are simply $\Psi_{i,j=0}(t_{k+1})=0$ with or without Coulomb potential.

For simplicity, let us introduce the short hand notations $X_{i,j}=\Psi_{i,j}(t_{k+1})$, $\Psi_{i,j}=\Psi_{i,j}(t_{k}),$ then by substituting the finite difference Laplacians (30), (31) into (33), we arrive at the following final form of linear equations for all $i\in\left[0,N_{z}\right]$, $j\in\left[1,N_{\rho}\right]$:

$$(-1+1/j)\beta_{\rho}X_{i,j-2}+(16-8/j)\beta_{\rho}X_{i,j-1}+(16+8/j)\beta_{\rho}X_{i,j+1}+(-1-1/j)\beta_{\rho}X_{i,j+2}$$

$$-\beta_{z}X_{i-2,j}+16\beta_{z}X_{i-1,j}+(1-30\beta_{\rho}-30\beta_{z}+\alpha V_{i,j})X_{i,j}+16\beta_{z}X_{i+1,j}-\beta_{z}X_{i+2,j}$$

$$=(1-1/j)\beta_{\rho}\Psi_{i,j-2}+(-16+8/j)\beta_{\rho}\Psi_{i,j-1}+(-16-8/j)\beta_{\rho}\Psi_{i,j+1}+(1+1/j)\beta_{\rho}\Psi_{i,j+2}$$

$$+\beta_{z}\Psi_{i-2,j}-16\beta_{z}\Psi_{i-1,j}+(1+30\beta_{\rho}+30\beta_{z}-\alpha V_{i,j})\Psi_{i,j}-16\beta_{z}\Psi_{i+1,j}+\beta_{z}\Psi_{i+2,j},$$

(36)

where

$$\alpha=i\Delta t/2,\,\,\,\beta_{\rho}=\alpha\beta/(12\Delta\rho^{2}),\,\,\,\beta_{z}=\alpha\beta/(12\Delta z^{2}).$$

(37)

For the $m=0$ configuration, the equations forced by the Neumann and Robin boundary conditions for all $i\in\left[0,N_{z}\right]$ from (32), (34), (35) are

$$(-25+(6(\gamma/\beta)\Delta\rho\cdot\delta_{R,i})X_{i,0}+48X_{i,1}-36X_{i,2}+16X_{i,3}-3X_{i,4}=0.$$

(38)

For the $m\neq 0$ states we have simply

$$X_{i,0}=0.$$

(39)

The spatial discretization in the cylindrical coordinate system and the Neumann or Robin boundary conditions for the $m=0$ states make the unitarity of the algorithm, and in a broader sense, the accuracy of spatial integrations a more subtle issue than usual. One has to find an appropriate discrete inner product formula that is conserved at least with the accuracy of the finite differences, and which can be evaluated with sufficient accuracy using the cylindrical grid. We give a solution to this auxiliary problem in A.

Although it corresponds to 3D propagation, we call this scheme 2D Crank-Nicolson method because it involves only two spatial coordinates. This is already a complete propagation algorithm by itself, however, it suffers from the numerically inefficient solution of the resulting linear equations: if we combine the $i,j$ indices into a single one (by flattening the two dimensional array) as $l=i\cdot(N_{\rho}+1)+j$, we obtain a a block pentadiagonal matrix of size $(N_{z}+1)^{2}(N_{\rho}+1)^{2}$, with block size $(N_{\rho}+1)^{2}$. Inverting this type of matrix is computationally intense [43, 53] because the width of the diagonal is $4N_{\rho}+1$: despite its apparent simplicity, the numerical cost of this task is $\sim N_{z}N_{\rho}^{3}$ which is extremely high compared to the cost $\sim N_{z}N_{\rho}$ of a pentadiagonal scheme of the same size. These facts inspired us to develop an improved algorithm which has almost all advantages of this 2D Crank-Nicolson method but needs much less numerical effort.

The approach of the operator splitting method is to factorize the exponential operator $e^{\lambda H}$ into multiple easy-to-solve parts. From the Taylor-expansion of the exponential operator

$$e^{\lambda(A+B)}=\sum_{n=0}^{\infty}\frac{\left(A+B\right)^{n}}{n!}\lambda^{n},$$

(40)

follow the two main formulae [54] which form the basis of the operator splitting schemes, namely the Baker-Campbell-Hausdorff formula :

$$e^{\lambda A}e^{\lambda B}=e^{\lambda(A+B)+\lambda^{2}\frac{1}{2}[B,A]+\lambda^{3}\frac{1}{12}[A-B,[A,B]]}+O(\lambda^{4}),$$

(41)

and the Zassenhaus formula

$$e^{\lambda(A+B)}=e^{\lambda A}e^{\lambda B}e^{\lambda^{2}\frac{1}{2}[A,B]}e^{\lambda^{3}\frac{1}{6}[A+2B,[A,B]]}+O(\lambda^{4}).$$

(42)

Both of these contain infinitely many commutators of $A$ and $B$. Of course, if $[A,B]=0$ then $e^{\lambda(A+B)}=e^{\lambda A}e^{\lambda B}$ exactly. As the above formulae suggest, the $O(\lambda^{4})$ terms can be further factorized into the exponents. An extended analysis is available in Refs. [55, 56].

If one uses (41) and (42) to acquire a symmetric decomposition, only odd leading order of $\lambda$ will appear in the formula. This is a requirement for quantum propagation though, because the presence of even order terms would destroy the unitary evolution of the wave function. A well-known example of this is the widely used standard symmetric second order accurate formula (or Strang splitting, after [57]):

$$e^{\lambda(A+B)}=e^{\lambda A/2}e^{\lambda B}e^{\lambda A/2}+C_{3}\lambda^{3}+O(\lambda^{4}),$$

(43)

where $C_{3}$ is a combination of commutators of $A$ and $B$. A direct fourth order splitting scheme was derived by Chin and Suzuki [58, 59]:

$$e^{\lambda(A+B)}=e^{\lambda\frac{1}{6}A}e^{\lambda\frac{1}{2}B}e^{\lambda\frac{2}{3}A+\lambda^{3}\frac{1}{72}[A,[B,A]]}e^{\lambda\frac{1}{2}B}e^{\lambda\frac{1}{6}A}+C_{5}\lambda^{5}+O(\lambda^{6}).$$

(44)

This splitting requires also the evaluation of the $[A,[B,A]]$ commutator, which can rise additional difficulties, depending on the particular form of $A$ and $B$.

We proceed by introducing another kind of higher order operator splitting, based on the work of Bandrauk and Shen [60], who developed an iterative method to improve the accuracy of the (43) scheme. Let us denote the second order accurate form with $S_{2}(\lambda)=e^{\lambda A/2}e^{\lambda B}e^{\lambda A/2}$, then their iteration method reads for $n=4,6,...$ as

$$S_{n}(\lambda)=S_{n-2}(s\lambda)S_{n-2}((1-2s)\lambda)S_{n-2}(s\lambda)+C_{n-1}(2s^{n-1}+(1-2s)^{n-1})\lambda^{n-1}+O(\lambda^{n+1})$$

(45)

where the parameter $s$ must be for each iteration step a real root of the corresponding polynomial equation

$\displaystyle 2s^{n-1}+(1-2s)^{n-1}$

$\displaystyle=$

$\displaystyle 0.$

(46)

In (45) only the odd error terms appear, because of the unitarity and the symmetry of the splitting scheme. So the $S_{n}(\lambda)$ requires $3^{n/2}$ evaluations of $S_{2}$, in the worst case. For completeness we note, that using the complex roots of (46) in (45) is also a viable approach in some numerical applications [59, 61].

This scheme was already generalized for time-dependent Hamiltonians of the form $H(t)=A(t)+B(t)$ in [60, 62] as follows. Inserting the second order effective Hamiltonian (22) into the formula (43) with $\lambda=-i\Delta t$ , the second order accurate splitting of the evolution operator becomes

$$U_{2}(t+\Delta t,\,\,t)=e^{-i\frac{\Delta t}{2}A(t+\Delta t/2)}e^{-i\,\Delta t\,B(t+\Delta t/2)}e^{-i\frac{\Delta t}{2}A(t+\Delta t/2)}.$$

(47)

Then (45) will take the form for $n=4,6...$

$$U_{n}(t+\Delta t,\,\,t)=U_{n-2}(t+\Delta t,\,\,t+(1-s)\Delta t)\,\,U_{n-2}(t+(1-s)\Delta t\,\,,t+s\Delta t)\,\,U_{n-2}(t+s\Delta t\,\,,t)+O(\lambda^{n+1})$$

(48)

and $s$ being the same as in the time-independent case (46). Interestingly, the fourth order approximation $U_{4}$ in (48) decreases the error even if the time evolution governed by a nonlinear time-dependent Schrödinger equations, if the nonlinear error term is corrected in the formulation of the $U_{2}$ propagator [62, 61].

Now we write out a new form of iteration (48) for time-dependent Hamiltonians, by using the same principles, for $n=6,10,...$

	$\displaystyle U_{n}(t+\Delta t,\,\,t)=\,\,$	$\displaystyle U_{n-4}(t+\Delta t,\,\,t+(1-s)\Delta t)\,\,U_{n-4}(t+(1-s)\Delta t,\,\,t+(1-s-p)\Delta t)\cdot\,\,$
		$\displaystyle U_{n-4}(t+(1-s-p)\Delta t,\,\,t+(s+p)\Delta t)\cdot\,\,$
		$\displaystyle U_{n-4}(t+(s+p)\Delta t,\,\,t+s\Delta t)\,\,U_{n-4}(t+s\Delta t,\,\,t)+O(\lambda^{n+1}),$		(49)

where $s$, $p$ must be the simultaneous real roots of equations

$$2s^{n-3}+2p^{n-3}+(1-2s-2p)^{n-3}=0,\,\,\,\,\,\text{and}\,\,\,\,\,\,2s^{n-1}+2p^{n-1}+(1-2s-2p)^{n-1}=0.$$

(50)

However, this formula for $S_{6}$ requires five evaluations of $S_{2}$, compared to nine in the case of the (45).

Although the fourth order formula (44) with the fourth order effective Hamiltonian (23) seems to be superior compared to the iterative propagation (48) (because of the extra information given by the temporal and spatial commutators, and less evaluations), the schemes (48) and (49) do not require the calculation of commutators, they decrease all the $\Delta t$ dependent errors simultaneously and they are easy to implement. However, they involve backward time steps, which means they do not work very well with diffusive problems.

Another class of high order split-operator methods for diffusive partial differential equations was developed in [63], where the exponential operator is found by an ansatz of

$$S_{2n}(\lambda)=\sum_{k=1}^{n}c_{k}S_{2}^{k}\left(\frac{\lambda}{k}\right)+O\left(\lambda^{2n+1}\right)\text{ with }c_{k}=\prod_{j=1(\neq k)}^{n}\frac{c_{k}^{2}}{c_{k}^{2}-c_{j}^{2}}.$$

(51)

Here $S_{2}$ is the same second order formula which is used in (45). The fourth order formula is simply given by $S_{4}(\lambda)=-\frac{1}{3}S_{2}(\lambda)+\frac{4}{3}S_{2}^{2}\left(\frac{\lambda}{2}\right)$. Thus, once $S_{2}$ is properly constructed, all high order formulae (48), (49), (51) can be utilized immediately. The (51) method is suitable for the imaginary time propagation [28, 63, 64] with $\lambda\rightarrow-\Delta t$, i.e. for determination of the lowest energy eigenstates of stationary potentials.

A problem also arises when one wishes to use an imaginary potential [65] with (48) and (49), which is a commonly used numerical technique to avoid artificial reflections from the (box) boundaries of the simulation domain. Adding an imaginary potential $iV_{{\rm im}}$ into the Hamiltonian gives the following splitting of the exponential operator from (43):

$$e^{-i\Delta t\,H+\Delta t\,V_{{\rm im}}}=e^{\Delta t\,V_{{\rm im}}/2}e^{-i\Delta t\,H}e^{\Delta t\,V_{{\rm im}}/2}+O(\Delta t^{3}).$$

(52)

It is clear that $V_{{\rm im}}(z,\rho)<0$ performs the required absorption of the wavefunction in the case of forward time propagation $\Delta t>0$, but it would blow up the wave function during the necessary backward steps in (48) and (49). We circumvented this problem by partially disallowing backward steps, i.e. we replaced $\Delta t\,V_{{\rm im}}$ by $|\Delta t|V_{{\rm im}}$.

The most common way of factorizing the exponential operator $e^{\lambda H}$ is that the different spatial coordinate derivatives decouple into different exponential operators, i.e. to use a directional splitting. Then the propagation can be carried out by solving multiple one dimensional TDSE-s in succession. The most frequent realization of this is the so-called potential-kinetic term splitting, which is mainly used in conjunction with the Fourier-transformation methods in Cartesian coordinates [38, 41, 42]. However, if we apply the potential-kinetic term splitting to problem of (1) then we have to write Hamiltonian as $H=T_{\rho}+T_{z}+V$ where $V$ is the potential, $T_{\rho}=\beta\partial_{\rho}^{2}+\beta\rho^{-1}\partial_{\rho}$ and $T_{z}=\beta\partial_{z}^{2}$ are the kinetic energy operators. The $[A,[B,A]]$ commutator will take a rather simple form of $[V,[T_{\rho}+T_{z},V]]=-2\beta|\nabla V|^{2}$. Thus, the direct second order (43) and fourth order (44) symmetric splitting schemes are written as

$$e^{\lambda(T_{\rho}+T_{z}+V)}=e^{\lambda V/2}e^{\lambda T_{\rho}}e^{\lambda T_{z}}e^{\lambda V/2}+C_{3}\lambda^{3}+O(\lambda^{4}),$$

(53)

$$e^{\lambda(T_{\rho}+T_{z}+V)}=e^{\lambda\frac{1}{6}V}e^{\lambda\frac{1}{2}T_{\rho}}e^{\lambda\frac{1}{2}T_{z}}e^{\lambda\frac{2}{3}V+\lambda^{3}\frac{1}{72}\frac{1}{\mu}|\nabla V|^{2}}e^{\lambda\frac{1}{2}T_{\rho}}e^{\lambda\frac{1}{2}T_{z}}e^{\lambda\frac{1}{6}V}+C_{5}\lambda^{5}+O(\lambda^{6}),$$

(54)

where $\lambda=-i\Delta t$. We also note that these split-operator formulae by themself will introduce error, scaling as $\Delta t^{3}$ and $\Delta t^{5}$, compared to the stationary states of the exact time-independent Hamiltonian. In the case of a time-dependent Hamiltonian, the respective second order (22) or fourth order (23) effective Hamiltonians should be used. These truncations of the Magnus-series (21) will introduce additional $\Delta t^{3}$ and $\Delta t^{5}$ dependent errors in time evolution.

The magnitudes of the aforementioned numerical errors will depend on the derivatives of $V(z,\rho,t)$: in the case of an external electric field the smoothness of $V$ in time is a reasonable assumption. On the other hand, the spatial dependence of the atomic potential should also behave the same way if we wish to use operator splitting. As we suspect, it will not always be the case, especially in atomic or molecular physics. We can deduce from (54) that the leading order of error of the second order scheme must have a dependence of $\Delta t^{3}|\nabla V|^{2}$, e.g. the error characteristics strongly depend on the magnitude of the spatial derivatives of $V(z,\rho)$. In the case of an atomic $1/r$ Coulomb potential this error term will be $\Delta t^{3}r^{-4}$ which becomes significant only in the region $r<1$, where it is increasing rapidly with fourth power of $1/r$, and the split operator scheme completely breaks down at $r=0$. This also illustrates the fact that non-differentiability of either the potential or the wavefunction could potentially make operator splitting like (53), (54) less-than-useful.

Up to now, we have seen that directional spitting related exponential operator factorizations result drastic speed improvement for well behaved potentials: we can directly apply the five-point finite difference Crank-Nicolson method to the $e^{\lambda T_{z}}$ and $e^{\lambda T_{\rho}}$exponents for evaluation, thus we have to solve for locally decoupled one dimensional wave functions. Then, the approximate operations count of evaluating the formula (53) is $\sim N_{z}N_{\rho}$, which is smaller by the factor of $N_{\rho}^{2}$ compared to the full Crank-Nicolson problem of the same size. On the other hand, the full Crank-Nicolson scheme is able to incorporate the Neumann and Robin boundary conditions at $\rho=0$ into the implicit linear equations properly, and it does not suffer from the catastrophic error blow up while we approach the point Coulomb singularity.

In this article we propose a merger of a split-operator method with the 2D Crank-Nicolson scheme, in the form of “hybrid splitting”, to get the best of both of these methods.

Let us split the spatial domain as $G=G_{\text{CN}}+G_{\text{Split}}$ where

$$G=\{z,\rho\in\mathbb{R},\,z_{\min}\leq z\leq z_{\min}+N_{z}\Delta z,\,\,0\leq\rho\leq N_{\rho}\Delta\rho\},$$

(55)

$$G_{\text{CN}}=\{z,\rho\in\mathbb{R},\,z_{\min}\leq z\leq z_{\min}+N_{z}\Delta z,\,\,0\leq\rho\leq L\Delta\rho\},$$

(56)

$$G_{\text{\text{Split}}}=\{z,\rho\in\mathbb{R},\,z_{\min}\leq z\leq z_{\min}+N_{z}\Delta z,\,\,L\Delta\rho<\rho\leq N_{\rho}\Delta\rho\},$$

(57)

then we define the pieces of the Hamiltonian $H=H_{z}+H_{\rho}+H_{{\rm CN}}$ as

$$H_{z}=\beta\partial_{z}^{2}\text{ \,\,\ if }(z,\rho)\in G_{\text{Split}},$$

(58)

$$H_{\rho}=\beta\partial_{\rho}^{2}+\beta\rho^{-1}\partial_{\rho}+V(z,\rho,t)\text{ \,\,\ if }(z,\rho)\in G_{\text{Split}},$$

(59)

$$H_{\text{CN}}=\beta\partial_{z}^{2}+\beta\partial_{\rho}^{2}+\beta\rho^{-1}\partial_{\rho}+V(z,\rho,t)\text{ \,\,\ if }(z,\rho)\in G_{\text{CN}}.$$

(60)

Then the original Hamiltonian can be reconstructed as

$$H=\begin{cases}H_{z}+H_{\rho}&\text{if }(z,\rho)\in G_{\text{Split}}\\ H_{\text{CN}}&\text{if }(z,\rho)\in G_{\text{CN}}\end{cases}.$$

(61)

$H$ will never get evaluated outside $G$, in accordance with the boundary conditions for the wave function $\Psi(z,\rho,t)$.