Fourth order real space solver for the time-dependent Schrödinger equation with singular Coulomb potential
FFourth order real space solver for the time-dependent Schrödinger equationwith singular Coulomb potential
Szilárd Majorosi a, ∗ , Attila Czirják a,b a Department of Theoretical Physics, University of Szeged, Tisza L. krt. 84-86, H-6720 Szeged, Hungary b ELI-ALPS, ELI-HU Non-Profit Ltd., Dugonics tér 13, H-6720 Szeged, Hungary
Abstract
We present a novel numerical method and algorithm for the solution of the 3D axially symmetric time-dependent Schrödinger equation in cylindrical coordinates, involving singular Coulomb potential terms be-sides a smooth time-dependent potential. We use fourth order finite difference real space discretization,with special formulae for the arising Neumann and Robin boundary conditions along the symmetry axis.Our propagation algorithm is based on merging the method of the split-operator approximation of the ex-ponential operator with the implicit equations of second order cylindrical 2D Crank-Nicolson scheme. Wecall this method hybrid splitting scheme because it inherits both the speed of the split step finite differenceschemes and the robustness of the full Crank-Nicolson scheme. Based on a thorough error analysis, we ver-ified both the fourth order accuracy of the spatial discretization in the optimal spatial step size range, andthe fourth order scaling with the time step in the case of proper high order expressions of the split-operator.We demonstrate the performance and high accuracy of our hybrid splitting scheme by simulating opticaltunneling from a hydrogen atom due to a few-cycle laser pulse with linear polarization.
Keywords:
Time-dependent Schrödinger equation, Singular Coulomb potential, High order methods,Operator splitting, Crank-Nicolson scheme, Strong field physics, Optical tunneling
1. Introduction
Experiments with attosecond light pulses [1–6], based on high-order harmonic generation (HHG) in noblegases [7–10], have been revolutionizing our view of fundamental atomic, molecular and solid state processesin this time domain [11–17]. A key step in gas-HHG is the tunnel ionization of a single atom and the returnof the just liberated electron to its parent ion [18, 19], due to the strong, linearly polarized femtosecondlaser pulse driving this process. Recent developments in attosecond physics revealed that the accuratedescription of this single-atom emission is more important than ever, e.g. for the correct interpretationof the data measured in attosecond metrology experiments, especially regarding the problem of the zeroof time [15, 20], and the problem of the exit momentum [13]. Although an intuitive and very successfulapproximate analytical solution [21] and many of its refinements exist [22], the most accurate descriptionof the single-atom response is given by the numerical solution of the time-dependent Schrödinger equation(TDSE).The common model of the single-atom response employs the single active electron approximation andthe dipole approximation. These reduce the problem to the motion of an electron in the time-dependentpotential formed by adding the effect of the time-dependent electric field to the atomic binding potential.The peculiarity of this problem is due to the strength of the electric field, which has its maximum typicallyin the range of 0.05-0.1 atomic units. Thus, this electric field enables the tunneling of the electron through ∗ Corresponding author.
Email address: [email protected] (Szilárd Majorosi)
Preprint submitted to Computer Physics Communications July 11, 2017 a r X i v : . [ phy s i c s . a t o m - ph ] J u l he (time-dependent) potential barrier formed by strongly distorting the atomic potential, but this effect isweak during the whole process. On the other hand, this small part of the wave function outside the barrierextends to large distances and in fact this is the main contribution to the time-dependent dipole moment,which is usually considered as the source of the emitted radiation. Thus, a very weak effect needs to becomputed very accurately, and these requirements get even more severe, if the model goes beyond the singleactive electron or the dipole approximation.The fundamental importance of the singularity of the Coulomb potential regarding HHG has alreadybeen analysed and emphasized e.g. in [23]. As other kind of numerical errors decrease, the artifacts causedby the smoothing of the singularity become more disturbing. Therefore, a high precision numerical methodhas to handle the Coulomb singularity as accurately and effectively as possible.The most widely employed approach to the numerical solution of the TDSE for an atomic electron drivenby a strong laser pulse is to use spherical polar coordinates: then the hydrogen eigenfunctions are analyticand one expands the wavefunction in Y l,m ( θ, ϕ ) spherical harmonics with expansion coefficients φ l,m ( r, t ) acting as radial or reduced radial functions, for which the problem is solved directly. For example, theseradial functions can be further expanded using a suitable basis, like Gaussian-Legendre or B-splines [24–26],in which the Hamiltonian matrix does not have any singularity. Also, refined solutions based on highlyaccurate real space discretization of the radial coordinate r exist using this harmonic expansion [27–31], alsorelying on the fact that a boundary condition for φ l,m (0 , t ) can be derived from the analytic properties ofhydrogenic eigenfunctions and can be built into the implicit equations. The usual drawback is that in orderto calculate physical quantities either the real space reconstruction of the wavefunction is necessary or everyobservable must be given in this harmonic expansion. Also, the spherical grid of these methods is not welladapted to the electron’s motion which is mainly along the polarization direction of the laser field.One can easily stumble upon applications other than HHG where neither spherical coordinates, nor spe-cial basis functions are optimal, and one has to use Cartesian or cylindrical coordinate systems in numericalsimulations, and of course, discretization of real space coordinates. The most basic method is just to usea staggered grid which avoids the singularity of the Coulomb potential [32]: then the scheme will be easilysolvable, although with low accuracy. A somewhat refined version of this is using a smoothed version of thepotential, for example the so called soft Coulomb-form of / (cid:112) α + | r | , where the value of α can be fittedby matching the ground state energy with that of the true hydrogen atom [32–35]. These approximationsare rather crude, but can reproduce important features of the quantum system [23]. Recently, Gordon et.al. [36] provided an improved way to include the Coulomb singularity, by using the general asymptoticbehavior of the hydrogen eigenstates around | r | = 0 . They built an improved discrete Hamiltonian matrixthrough a new discrete potential, using a three point finite difference scheme. This way they increased theaccuracy of the numerical solution by an order of magnitude.In the present work, we propose a novel method to accurately incorporate the effect of the Coulombsingularity in the solution of a spinless, axially symmetric TDSE, in cylindrical coordinates. We accountfor the singular Coulomb potential through analytic boundary condition which resolves the problem ofnondifferentiability of the hydrogen eigenfunctions in the cylindrical coordinate system and yields highorder spatial accuracy.Incorporating this kind of analytic boundary condition into the numerical solution of TDSE will constrainthe applicable propagation methods because the singularity of the Coulomb potential at r = 0 and thesingularity of the cylindrical radial coordinate at ρ = 0 will introduce Robin and Neumann boundaryconditions, overriding the Hamilton operator. As a consequence, typical explicit methods, for examplestaggered leapfrog [36, 37], second order symmetric difference method [38, 39], any polynomial expansionmethod [38, 40], and split-step Fourier transformation methods [38, 41, 42] have been ruled out, leaving onlyimplicit methods.In Section 2, we introduce the cylindrical TDSE to be solved, along with its boundary conditions, andwe describe the way how we incorporate the singular Coulomb potential to the problem. In Section 3,we derive the suitable finite-difference scheme with high order spatial accuracy, using second order Crank-Nicolson method [38, 43, 44]. We briefly go through the standard approximations of the evolution operator,then we write down the detailed form of the resulting implicit linear equations. In Section 4, we overviewthe standard split-operator methods for numerically solving the TDSE, and combine them with the finite-2ifference description, which is known as split-step finite difference method [45, 46]. We show that similarlyto the typical explicit methods, real space split-operator methods do not work well for singular Hamiltonians.Based in these conclusions, we introduce our hybrid splitting method, which retains the robustness of thefull Crank-Nicolson method and the speed of the split-operator method. In Section 5, we analyse the newlyarisen coefficient matrix, we describe the way how to take advantage of its structure, then we derive a specialsolver algorithm which is necessary for fast solution of the TDSE. In Section 6, we discuss the results ofnumerical experiments we carried out with our hybrid splitting method: (i) for the accuracy of the spatialdiscretization of the finite difference scheme, we calculate the eigenenergies of selected eigenstates, (ii) forthe temporal accuracy of our split-step finite difference scheme combined with high order evolution operatorfactorizations, we compare the numerical wave function with the analytic solution of the quantum forcedharmonic oscillator, (iii) for final verification, we simulate the hydrogen atom in a time dependent externalelectric field.
2. Schrödinger equation and boundary conditions
Let us consider the axially symmetric three-dimensional time-dependent Schrödinger equation in cylin-drical coordinates ρ = (cid:112) x + y and z : i (cid:126) ∂∂t Ψ ( z, ρ, t ) = − (cid:126) µ (cid:20) ∂ ∂z + ∂ ∂ρ + 1 ρ ∂∂ρ (cid:21) Ψ ( z, ρ, t ) + V ( z, ρ, t ) Ψ ( z, ρ, t ) (1)where µ is the (reduced) electron mass. For the sake of simplicity, we use atomic units throughout thisarticle and we use the notation β = − (cid:126) / µ .As discussed in the Introduction, such a Hamiltonian plays a fundamental role e.g. in the interactionof atomic systems with strong laser pulses. Then the potential is usually the sum of an atomic bindingpotential (centered in the origin) and an interaction energy term. For a single active electron and a linearlypolarized laser pulse, this interaction term in dipole approximation and in length gauge [47] reads as: V ( z, ρ, t ) = V atom ( z, ρ ) − q · F ( t ) · z, (2)where q = − 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 ∈ [ z min , z max ] , ρ ∈ [0 , ρ max ] and the problem is treatedas an initial value problem as Ψ ( z, ρ,
0) = ψ ( z, ρ ) is known. Box boundary conditions are posed at threeedges of the interval: Ψ ( z > z max , ρ, t ) = Ψ ( z < z min , ρ, t ) = Ψ ( z, | ρ | > ρ max , t ) = 0 (3)and because of the singularity in the curvilinear coordinate ρ , Neumann boundary condition should be posedalong the line ρ = 0 . In order to find this boundary condition, we multiply (1) by ρ and take the limit ρ → : ∂ Ψ ∂ρ (cid:12)(cid:12)(cid:12)(cid:12) ρ =0 = − β lim ρ → ( ρV ( z, ρ, t ) Ψ ( z, ρ, t )) = 0 (4)where we have assumed the potential is smooth at ρ = 0 line. Then Ψ( z, ρ, t ) is continuous and continuouslydifferentiable everywhere in the z − ρ plane, and it has extremal value in ρ = 0 . Symmetry condition Ψ( z, ρ, t ) = Ψ( z, − ρ, t ) also holds in this case.However, this Neumann boundary condition has to be altered when we introduce a singular potential inany form, e.g. the 3D Coulomb potential V atom ( z, ρ ) = γ/r with r = (cid:112) z + ρ and γ = − q Z/ π(cid:15) .3 .2. The boundary condition for the 3D Coulomb potential Our aim is to represent the singularity of the 3D Coulomb potential by implementing the correct boundarycondition 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 ofthe radial solutions φ ( r ) of the well known radial equation [48, 49]: (cid:20) β ∂ ∂r + 2 βr ∂∂r + β l ( l + 1) r + γr (cid:21) φ ( r ) = Eφ ( r ) , (5)where l is the angular momentum quantum number. If we multiply both sides with r then take the r → limit, we can acquire Robin boundary condition for all of the l = 0 states as ∂φ∂r (cid:12)(cid:12)(cid:12)(cid:12) r =0 = − γ β φ (0) , (6)but the bound states with higher angular momenta will have φ (0) = 0 boundary condition at the origin,in accordance with the fact that for l > the particle physically cannot penetrate to r = 0 , because of thesingularity of l ( l + 1) /r .Since we need the boundary conditions for the cylindrical equation independent of the l values, we turnto the nonseparated symmetric spherical Schrödinger equation for answers: (cid:20) β ∂ ∂r + 2 βr ∂∂r + β r sin θ ∂∂θ (cid:18) sin θ ∂∂θ (cid:19) + γr (cid:21) Ψ( r, θ ) = E Ψ( r, θ ) (7)which has also the term /r and has singularity at r = 0 and at θ = kπ , k = 0 , , . . . . To acquire theboundary conditions for r = 0 , we again multiply both sides with r and take the r → limit, yielding lim r → (cid:20) β ∂∂r + β cos θr sin θ ∂∂θ + βr ∂ ∂ θ + γ (cid:21) Ψ( r, θ ) = 0 . (8)Now we transform equation (8) into cylindrical coordinates z = r cos θ and ρ = r sin θ , using followingexpressions of the partial derivatives: ∂∂r = ∂z∂r ∂∂z + ∂ρ∂r ∂∂ρ = sin θ ∂∂ρ + cos θ ∂∂z , (9) ∂∂θ = ∂z∂θ ∂∂z + ∂ρ∂θ ∂∂ρ = r cos θ ∂∂ρ − r sin θ ∂∂z (10)with cos θ = z/ (cid:112) z + ρ and sin θ = ρ/ (cid:112) z + ρ . After writing back and taking the limit, the formula (8)becomes (cid:34) β ρ (cid:112) z + ρ ∂∂ρ + β z (cid:112) z + ρ ∂∂z + β z ρ (cid:112) z + ρ ∂∂ρ + γ (cid:35) Ψ( r, θ ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) r =0 = 0 . (11)Then, by substituting z = 0 into (11), we obtain the Robin boundary condition for the 3D Coulombsingularity: ∂ Ψ ∂ρ (cid:12)(cid:12)(cid:12)(cid:12) ρ,z =0 = − γ β Ψ (0 , , t ) . (12)This can be generalized to include multiple Coulomb-cores rather easily, we just need to impose (12) atmultiple z, ρ points along the ρ = 0 axis. Additionally, any continuously differentiable potential added tothis configuration does not change this boundary condition. It is interesting to note also that the boundarycondition imposed by a 1D Dirac-delta potential γδ ( ρ ) also has the form of (12) for a symmetric wavefunction [49]. 4 .3. The states with nonzero magnetic quantum numbers 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 stateswith nonzero m also lead to a TDSE of the form (13), however with different boundary conditions at ρ = 0 .Let us write out the full 3D cylindrical time-dependent Schrödinger equation, again with the same axiallysymmetric potential that we used previously: i ∂∂t Ξ ( z, ρ, φ, t ) = (cid:20) β ∂ ∂z + β ∂ ∂ρ + βρ ∂∂ρ + βρ ∂ ∂φ + V ( z, ρ, t ) (cid:21) Ξ ( z, ρ, φ, t ) . (13)We can write the quantum state with a given magnetic quantum number m in the following form for any t Ξ ( z, ρ, φ, t ) = Ψ m ( z, ρ, t ) e imφ , (14)due to the fact that the Hamiltonian of (13) conserves the m quantum number. (We can also use real angularbasis functions like sin( mφ ) or cos( mφ ) instead of e imφ , as long as that they are the eigenfunctions of ∂ φ .)After substituting (14) to (13) we perform projection in the form of ´ π e − imφ ( · )d φ . Then we get theTDSE for our wave function Ψ m ( z, ρ, t ) with the specified magnetic quantum number m : i ∂∂t Ψ m ( z, ρ, t ) = (cid:20) β ∂ ∂z + β ∂ ∂ρ + βρ ∂∂ρ − β m ρ + V ( z, ρ, t ) (cid:21) Ψ m ( z, ρ, t ) , (15)Eq. (15) has the same form as the axially symmetric equation (1) but with a new m dependent potential: V m ( z, ρ, t ) = V ( z, ρ, t ) − βm /ρ . (16)However, if m (cid:54) = 0 then the /ρ singularity of (16) alters the boundary condition (4) at the ρ = 0 axis tothe following: Ψ m ( z, , t ) = 0 , (17)which can by derived by multiplying (15) with ρ and taking the limit of ρ → . The boundary condition(17) also works with Coulomb potential, and in every case where ρ V ( z, ρ, t ) = 0 in the limit of ρ → . Itis also consistent with the boundary conditions for Coulomb states with nonzero angular momenta ( l (cid:54) = 0 )mentioned in the previous section.
3. The Crank-Nicolson method
The formal solution of (1) in terms of the time evolution operator U ( t, t (cid:48) ) = T exp (cid:2) − i (cid:126) ´ H ( t (cid:48)(cid:48) ) dt (cid:48)(cid:48) (cid:3) isthe following | Ψ( t ) (cid:105) = U ( t, t (cid:48) ) | Ψ( t (cid:48) ) (cid:105) . (18)Here, the exponential operator is to be understood as a time-ordered quantity, which is a difficult procedureif the Hamilton operators at different time instants do not commute. However, this is the case for thepotential we are considering.To acquire suitable discretization in the time domain of problem (1) we approximate the U ( t, timeevolution operator directly. First, let us divide the time domain [0 , t ] , into N t equal subintervals, then bythe group property of U ( t, we get U ( t,
0) = N t − (cid:89) k =0 U ( t k +1 , t k ) (19)5here t k = k ∆ t and ∆ t = t/N t . In one interval, we can write the evolution operator with its short timeform of U ( t k +1 , t k ) = e − i ∆ t H k (20)where H k is the effective time-independent Hamiltonian related to the original one H ( t ) by the Magnusexpansion [50, 51]: H k = 1∆ t t k +1 ˆ t k H ( t (cid:48) ) d t (cid:48) + i t t k +1 ˆ t k t (cid:48) ˆ t k [ H ( t (cid:48)(cid:48) ) , H ( t (cid:48) )] d t (cid:48)(cid:48) d t (cid:48) + . . . . (21)This includes infinitely many commutators of the Hamiltonians evaluated at different time points, to beintegrated with respect to more and more variables.Using only the first term of this series, one can directly acquire the well known second order approximationof the Hamiltonian as H (2) k = H ( t k + ) . (22)In order to have information about the next error term of the time evolution, we need to evaluate the Magnuscommutator series to fourth order. According to Puzynin et. al. [50], this improved approximation for aTDSE of the form (1) is H (4) k = 12 µ (cid:18) − i ∇ + ∆ t ∇ ˙ V ( t k +1 / ) (cid:19) + V ( t k +1 / ) + ∆ t
24 ¨ V ( t k +1 / ) . (23)where the top dots are the short hand notation for ∂ t time derivatives. This formula is important even ifwe use only the second order approximation, because the leading order of error depends on characteristicsof the first and second time derivatives of V ( z, ρ, t ) .For the exponential operators (20), first we consider the diagonal Padé-approximation of an exponentialfunction e λ · x = F ( − M, − M, λ · x ) F ( − M, − M, − λ · x ) + O ( λ M +1 ) (24)where F 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 λ = − i ∆ t , and itcan be shown that for self-adjoint operators the approximation is unitary [23].From this, a generalized operator approximation scheme that is ∆ t M +1 accurate can be obtained [44, 52]in a general form of e − i ∆ t H k = M (cid:89) s =1 (cid:34)(cid:18) i ∆ tx ∗ s H k (cid:19) − (cid:18) − i ∆ tx s H k (cid:19)(cid:35) + O (∆ t M +1 ) (25)where x s for s = 1 , . . . , M are the roots of the polynomial equation F ( − M, − M, − x ) = 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: Ψ( t k +1 ) = (cid:18) iH k ∆ t (cid:19) − (cid:18) − iH k ∆ t (cid:19) Ψ( t k ) . (27)One can straightforwardly construct higher order Crank-Nicolson schemes using (25), yielding multipleimplicit Crank-Nicolson substeps [44]. For time dependent Hamiltonians though, a high order accuratescheme should use the corresponding high order effective Hamiltonian (21) for consistency.6 .2. The finite difference scheme 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, ρ cylindrical coordinates and thefollowing equidistant 2D spatial grid: z i = z min + i · ∆ z, ∆ z = ( z max − z min ) /N z , i ∈ [0 , N z ] , (28) ρ j = j · ∆ ρ, ∆ ρ = ρ max /N ρ , j ∈ [0 , N ρ ] . (29)The first term of the Hamiltonian H k is the Laplacian ∇ . We denote the finite difference approximationof its z - and ρ -dependent terms by L z and L ρ , respectively, and we use the following fourth order accuratefinite difference forms [52] for them: L z Ψ i,j ( t ) = − Ψ i − ,j + 16Ψ i − ,j − i,j + 16Ψ i +1 ,j − Ψ i +2 ,j z , (30) L ρ Ψ i,j ( t ) = ( − /j )Ψ i,j − + (16 − /j )Ψ i,j − − i,j + (16 + 8 /j )Ψ i,j +1 + ( − − /j )Ψ i,j +2 ρ . (31)These fourth order formulae are optimal in the sense that symmetric finite differences of more than fivepoints are very complicated to be applied for the Coulomb-problem, because the higher order the finitedifference formula is, the higher derivatives of Ψ( z, ρ, t ) should be continuous.To discretize the Neumann and Robin boundary conditions, we need a one-sided finite difference formulafor the first derivative. Based on the method of Ref. [52], we derived the following fourth order accurateforward difference operator: D ρ Ψ i,j ( t ) = − i,j + 48Ψ i,j +1 − i,j +2 + 16Ψ i,j +3 − i,j +4 ρ . (32)Using standard second order accurate form (27) of the exponential operator with the discretized Laplacianin the Hamilton matrix, we get the following implicit scheme for all i ∈ [0 , N z ] and j ∈ [1 , N ρ ] : (1 + αβL z + αβL ρ + αV i,j ) Ψ i,j ( t k +1 ) = (1 − αβL z − αβL ρ − αV i,j ) Ψ i,j ( t k ) , (33)where α = i ∆ t/ and the potential is evaluated at the temporal midpoint V i,j = V ( z i , ρ j , t k +1 / ) . The boxboundary conditions Ψ i,j = 0 is applied if j > N ρ or i > N z or i < .Assuming that the potential is smooth, the Neumann boundary condition (4) prescribes the followingimplicit equations at j = 0 for all i ∈ [0 , N z ] : D ρ Ψ i, ( t k +1 ) = 0 . (34)On the other hand, if a Coulomb-core of strength γ 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 (cid:18) D ρ + γ β (cid:19) Ψ R, ( 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 theinterval [0 , N z ] if it is reasonably far from the box boundaries.For m (cid:54) = 0 states, the equations of the boundary conditions at ρ = 0 are simply Ψ i,j =0 ( t k +1 ) = 0 withor without Coulomb potential.For simplicity, let us introduce the short hand notations X i,j = Ψ i,j ( t k +1 ) , Ψ i,j = Ψ i,j ( t k ) , then bysubstituting the finite difference Laplacians (30), (31) into (33), we arrive at the following final form oflinear equations for all i ∈ [0 , N z ] , j ∈ [1 , N ρ ] : 7 − /j ) β ρ X i,j − + (16 − /j ) β ρ X i,j − + (16 + 8 /j ) β ρ X i,j +1 + ( − − /j ) β ρ X i,j +2 − β z X i − ,j + 16 β z X i − ,j + (1 − β ρ − β z + αV i,j ) X i,j + 16 β z X i +1 ,j − β z X i +2 ,j = (1 − /j ) β ρ Ψ i,j − + ( −
16 + 8 /j ) β ρ Ψ i,j − + ( − − /j ) β ρ Ψ i,j +1 + (1 + 1 /j ) β ρ Ψ i,j +2 + β z Ψ i − ,j − β z Ψ i − ,j + (1 + 30 β ρ + 30 β z − αV i,j )Ψ i,j − β z Ψ i +1 ,j + β z Ψ i +2 ,j , (36)where α = i ∆ t/ , β ρ = αβ/ (12∆ ρ ) , β z = αβ/ (12∆ z ) . (37)For the m = 0 configuration, the equations forced by the Neumann and Robin boundary conditions for all i ∈ [0 , N z ] from (32), (34), (35) are ( −
25 + (6( γ/β )∆ ρ · δ R,i ) X i, + 48 X i, − X i, + 16 X i, − X i, = 0 . (38)For the m (cid:54) = 0 states we have simply X i, = 0 . (39)The spatial discretization in the cylindrical coordinate system and the Neumann or Robin boundaryconditions for the m = 0 states make the unitarity of the algorithm, and in a broader sense, the accuracy ofspatial integrations a more subtle issue than usual. One has to find an appropriate discrete inner productformula that is conserved at least with the accuracy of the finite differences, and which can be evaluated withsufficient accuracy using the cylindrical grid. We give a solution to this auxiliary problem in Appendix A.Although it corresponds to 3D propagation, we call this scheme 2D Crank-Nicolson method because itinvolves 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 · ( N ρ + 1) + j , we obtain a ablock pentadiagonal matrix of size ( N z + 1) ( N ρ + 1) , with block size ( N ρ + 1) . Inverting this type ofmatrix is computationally intense [43, 53] because the width of the diagonal is N ρ + 1 : despite its apparentsimplicity, the numerical cost of this task is ∼ N z N ρ which is extremely high compared to the cost ∼ N z N ρ of a pentadiagonal scheme of the same size. These facts inspired us to develop an improved algorithm whichhas almost all advantages of this 2D Crank-Nicolson method but needs much less numerical effort.
4. Operator splitting schemes
As we have seen in the previous section, the 2D Crank-Nicolson scheme with the boundary condition (4)is a possible but ineffective way of solving the TDSE numerically. Now we are going to discuss how to applythe well-known method of operator splitting to the solution of the problem described in the introduction.
The approach of the operator splitting method is to factorize the exponential operator e λH into multipleeasy-to-solve parts. From the Taylor-expansion of the exponential operator e λ ( A + B ) = ∞ (cid:88) n =0 ( A + B ) n n ! λ n , (40)follow the two main formulae [54] which form the basis of the operator splitting schemes, namely the Baker-Campbell-Hausdorff formula : e λA e λB = e λ ( A + B )+ λ [ B,A ]+ λ [ A − B, [ A,B ]] + O ( λ ) , (41)and the Zassenhaus formula 8 λ ( A + B ) = e λA e λB e λ [ A,B ] e λ [ A +2 B, [ A,B ]] + O ( λ ) . (42)Both of these contain infinitely many commutators of A and B . Of course, if [ A, B ] = 0 then e λ ( A + B ) = e λA e λB exactly. As the above formulae suggest, the O ( λ ) 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 λ will appearin the formula. This is a requirement for quantum propagation though, because the presence of even orderterms would destroy the unitary evolution of the wave function. A well-known example of this is the widelyused standard symmetric second order accurate formula (or Strang splitting, after [57]): e λ ( A + B ) = e λA/ e λB e λA/ + C λ + O ( λ ) , (43)where C is a combination of commutators of A and B . A direct fourth order splitting scheme was derivedby Chin and Suzuki [58, 59]: e λ ( A + B ) = e λ A e λ B e λ A + λ [ A, [ B,A ]] e λ B e λ A + C λ + O ( λ ) . (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 Bandraukand Shen [60], who developed an iterative method to improve the accuracy of the (43) scheme. Let us denotethe second order accurate form with S ( λ ) = e λA/ e λB e λA/ , then their iteration method reads for n = 4 , , ... as S n ( λ ) = S n − ( sλ ) S n − ((1 − s ) λ ) S n − ( sλ ) + C n − (2 s n − + (1 − s ) n − ) λ n − + O ( λ n +1 ) (45)where the parameter s must be for each iteration step a real root of the corresponding polynomial equation s n − + (1 − s ) n − = 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 ( λ ) requires n/ evaluations of S , in the worst case. For completeness we note, that using thecomplex 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 λ = − i ∆ t , the second order accurate splitting of the evolution operator becomes U ( t + ∆ t, t ) = e − i ∆ t A ( t +∆ t/ e − i ∆ t B ( t +∆ t/ e − i ∆ t A ( t +∆ t/ . (47)Then (45) will take the form for n = 4 , ...U n ( t + ∆ t, t ) = U n − ( t + ∆ t, t + (1 − s )∆ t ) U n − ( t + (1 − s )∆ t , t + s ∆ t ) U n − ( t + s ∆ t , t ) + O ( λ n +1 ) (48)and s being the same as in the time-independent case (46). Interestingly, the fourth order approximation U in (48) decreases the error even if the time evolution governed by a nonlinear time-dependent Schrödingerequations, if the nonlinear error term is corrected in the formulation of the U propagator [61, 62].Now we write out a new form of iteration (48) for time-dependent Hamiltonians, by using the sameprinciples, for n = 6 , , ...U n ( t + ∆ t, t ) = U n − ( t + ∆ t, t + (1 − s )∆ t ) U n − ( t + (1 − s )∆ t, t + (1 − s − p )∆ t ) · U n − ( t + (1 − s − p )∆ t, t + ( s + p )∆ t ) · U n − ( t + ( s + p )∆ t, t + s ∆ t ) U n − ( t + s ∆ t, t ) + O ( λ n +1 ) , (49)where s , p must be the simultaneous real roots of equations s n − + 2 p n − + (1 − s − p ) n − = 0 , and s n − + 2 p n − + (1 − s − p ) n − = 0 . (50)9owever, this formula for S requires five evaluations of S , compared to nine in the case of the (45).Although the fourth order formula (44) with the fourth order effective Hamiltonian (23) seems to besuperior compared to the iterative propagation (48) (because of the extra information given by the temporaland spatial commutators, and less evaluations), the schemes (48) and (49) do not require the calculationof commutators, they decrease all the ∆ 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 devel-oped in [63], where the exponential operator is found by an ansatz of S n ( λ ) = n (cid:88) k =1 c k S k (cid:18) λk (cid:19) + O (cid:0) λ n +1 (cid:1) with c k = n (cid:89) j =1( (cid:54) = k ) c k c k − c j . (51)Here S is the same second order formula which is used in (45). The fourth order formula is simply given by S ( λ ) = − S ( λ ) + S (cid:0) λ (cid:1) . Thus, once S 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 λ → − ∆ 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), whichis a commonly used numerical technique to avoid artificial reflections from the (box) boundaries of thesimulation domain. Adding an imaginary potential iV im into the Hamiltonian gives the following splittingof the exponential operator from (43): e − i ∆ t H +∆ t V im = e ∆ t V im / e − i ∆ t H e ∆ t V im / + O (∆ t ) . (52)It is clear that V im ( z, ρ ) < performs the required absorption of the wavefunction in the case of forwardtime propagation ∆ t > , but it would blow up the wave function during the necessary backward stepsin (48) and (49). We circumvented this problem by partially disallowing backward steps, i.e. we replaced ∆ t V im by | ∆ t | V im . The most common way of factorizing the exponential operator e λH is that the different spatial coordi-nate derivatives decouple into different exponential operators, i.e. to use a directional splitting. Then thepropagation can be carried out by solving multiple one dimensional TDSE-s in succession. The most fre-quent realization of this is the so-called potential-kinetic term splitting, which is mainly used in conjunctionwith the Fourier-transformation methods in Cartesian coordinates [38, 41, 42]. However, if we apply thepotential-kinetic term splitting to problem of (1) then we have to write Hamiltonian as H = T ρ + T z + V where V is the potential, T ρ = β∂ ρ + βρ − ∂ ρ and T z = β∂ z are the kinetic energy operators. The [ A, [ B, A ]] commutator will take a rather simple form of [ V, [ T ρ + T z , V ]] = − β |∇ V | . Thus, the direct second order(43) and fourth order (44) symmetric splitting schemes are written as e λ ( T ρ + T z + V ) = e λV/ e λT ρ e λT z e λV/ + C λ + O ( λ ) , (53) e λ ( T ρ + T z + V ) = e λ V e λ T ρ e λ T z e λ V + λ µ |∇ V | e λ T ρ e λ T z e λ V + C λ + O ( λ ) , (54)where λ = − i ∆ t . We also note that these split-operator formulae by themself will introduce error, scaling as ∆ t and ∆ t , compared to the stationary states of the exact time-independent Hamiltonian. In the case ofa time-dependent Hamiltonian, the respective second order (22) or fourth order (23) effective Hamiltoniansshould be used. These truncations of the Magnus-series (21) will introduce additional ∆ t and ∆ t dependenterrors in time evolution.The magnitudes of the aforementioned numerical errors will depend on the derivatives of V ( z, ρ, t ) : inthe case of an external electric field the smoothness of V in time is a reasonable assumption. On the otherhand, the spatial dependence of the atomic potential should also behave the same way if we wish to useoperator splitting. As we suspect, it will not always be the case, especially in atomic or molecular physics.10 z N ρ L ρ z N ρ -L N ρ -LN ρ -L ~N z N ρ Splittingerror~ Δ t³ ~N z (N ρ -L) ~N z L Operations : ~N z (N ρ -L) ~N z (N ρ -L)
2D Crank-Nicolson Hybrid Crank-Nicolson Based Operator Splitting e ≈ e e e -i ∆ tH -i ∆ tH z /2 -i ∆ t ( H CN +H ρ ) -i ∆ tH z /2 Figure 1: Sketch of our hybrid splitting scheme. The costly 2D Crank-Nicolson scheme was replaced with a special secondorder symmetric split operator formula except at the L nearest gridpoints in the neighborhood of the ρ = 0 axis, in order toretain accuracy and stability. The solid lines represent coupling between the gridpoints of the exponential operator evaluations.We also indicate the approximate operations count needed to solve the respective systems of linear equations. We can deduce from (54) that the leading order of error of the second order scheme must have a dependenceof ∆ t |∇ V | , e.g. the error characteristics strongly depend on the magnitude of the spatial derivatives of V ( z, ρ ) . In the case of an atomic /r Coulomb potential this error term will be ∆ t r − which becomessignificant only in the region r < , where it is increasing rapidly with fourth power of /r , and the splitoperator scheme completely breaks down at r = 0 . This also illustrates the fact that non-differentiability ofeither 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 drasticspeed improvement for well behaved potentials: we can directly apply the five-point finite difference Crank-Nicolson method to the e λT z and e λT ρ exponents for evaluation, thus we have to solve for locally decoupledone dimensional wave functions. Then, the approximate operations count of evaluating the formula (53) is ∼ N z N ρ , which is smaller by the factor of N ρ compared to the full Crank-Nicolson problem of the samesize. On the other hand, the full Crank-Nicolson scheme is able to incorporate the Neumann and Robinboundary conditions at ρ = 0 into the implicit linear equations properly, and it does not suffer from thecatastrophic 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, inthe form of “hybrid splitting”, to get the best of both of these methods. Let us split the spatial domain as G = G CN + G Split where G = { z, ρ ∈ R , z min ≤ z ≤ z min + N z ∆ z, ≤ ρ ≤ N ρ ∆ ρ } , (55) G CN = { z, ρ ∈ R , z min ≤ z ≤ z min + N z ∆ z, ≤ ρ ≤ L ∆ ρ } , (56) G Split = { z, ρ ∈ R , z min ≤ z ≤ z min + N z ∆ z, L ∆ ρ < ρ ≤ N ρ ∆ ρ } , (57)then we define the pieces of the Hamiltonian H = H z + H ρ + H CN as H z = β∂ z if ( z, ρ ) ∈ G Split , (58) H ρ = β∂ ρ + βρ − ∂ ρ + V ( z, ρ, t ) if ( z, ρ ) ∈ G Split , (59) H CN = β∂ z + β∂ ρ + βρ − ∂ ρ + V ( z, ρ, t ) if ( z, ρ ) ∈ G CN . (60)11hen the original Hamiltonian can be reconstructed as H = (cid:40) H z + H ρ if ( z, ρ ) ∈ G Split H CN if ( z, ρ ) ∈ G CN . (61) H will never get evaluated outside G , in accordance with the boundary conditions for the wave function Ψ( z, ρ, t ) .Thus, we have partitioned the spatial domain into two regions: G CN where (based on the previoussection) we do not use any split operator approximation and propagate with Hamiltonian H CN , and region G Split where we do use operator splitting as H z + H ρ .In order to merge the directional operator splitting approach with the true 2D Crank-Nicolson method,we introduce our second order hybrid splitting scheme: e − i ∆ t ( H z + H ρ + H CN ) = e − i ∆ t H z / e − i ∆ t ( H ρ + H CN ) e − i ∆ t H z / + O (∆ t ) , (62)where we keep H CN and H ρ in the same exponent, in order that the wave function can “freely flow” betweenthe two regions G CN and G Split without introducing further artifacts. (Note that the exponential operator e − i ∆ t ( H ρ + H CN ) cannot be split further in this sense.) We use this scheme as the second order terms in theiterative formulae (45), (49), (51) to gain higher order accuracy in ∆ t . The e − i ∆ t H z / part can be evaluatedwith any method of choice. We have constructed a Numerov-extended Crank-Nicolson line propagationalgorithm for this, which can be found in Appendix B.In order to evaluate the e − i ∆ t ( H ρ + H CN ) operator, we apply the second order Padé-approximation (27)to arrive again at a second order Crank-Nicolson form of Ψ( t k +1 ) = e − i ∆ t ( H ρ + H CN ) Ψ( t k ) = (1 + α ( H ρ + H CN )) − (1 − α ( H ρ + H CN )) Ψ( t k ) + O (∆ t ) (63)where α = i ∆ t/ and V ( z, ρ, t ) is evaluated at the midpoint t + ∆ t/ . We introduce the spatial grids (28),(29) and the discrete operators L z , L ρ and D ρ , then we get the following equations for the two regions for i ∈ [0 , N z ] : (1 + αβL ρ + αβL z + αV i,j ) Ψ i,j ( t k +1 ) = (1 − αβL ρ − αβL z − αV i,j ) Ψ i,j ( t k ) , if j ∈ [1 , L ] , (64) (1 + αβL ρ + V i,j ) Ψ i,j ( t k +1 ) = (1 − αβL ρ − αV i,j ) Ψ i,j ( t k ) , if j ∈ [ L + 1 , N ρ ] . (65)Using five point finite differences, the first set of equations can be expanded resulting in the form of (36),and the expansion of the second set gives the following similar result: ( − /j ) β ρ X i,j − + (16 − /j ) β ρ X i,j − + (1 − β ρ + αV i,j ) X i,j +(16 + 8 /j ) β ρ X i,j +1 + ( − − /j ) β ρ X i,j +2 = (1 − /j ) β ρ Ψ i,j − + ( −
16 + 8 /j ) β ρ Ψ i,j − +(1 + 30 β ρ − αV i,j )Ψ i,j + ( − − /j ) β ρ Ψ i,j +1 + (1 + 1 /j ) β ρ Ψ i,j +2 . (66)With the presence of a Coulomb-core at z R = 0 , the boundary condition at ρ = 0 will be once again for the m = 0 configuration ( D ρ + γ/ (2 β ) · δ R,i ) Ψ i, ( t k ) = 0 (67)with the same expanded form as (38). For the m (cid:54) = 0 states, we just use the Dirichlet-boundary condition(39) instead of (67). The box boundary conditions also apply for every m : Ψ i,j ( t k +1 ) = 0 , if i / ∈ [0 , N z ] , j / ∈ [ − N ρ , N ρ ] . (68)From these, a mixed 2D - 1D Crank-Nicolson scheme can be constructed in the G CN + G Split domaindepending on m , which is fourth order accurate in space and second order accurate in time.What are the advantages of this scheme? First, if L > / ∆ ρ the accuracy of the directional splittingcan be considerably increased in the presence of the Coulomb potential, because we removed directional12plitting near its core (where its gradient is the largest). Second, if L ≥ (the “width” of D ρ ) then the (67)condition no longer affects into the split operator zone: we will maintain stability. Third, the operation-count is approximately ∼ L N z plus ∼ ( N ρ − L ) N z , meaning if L (cid:28) N ρ then we can regain the speed of thedirectional splitting and large part of the algorithm, corresponding to equations (65) can be parallelized fordifferent j indices. So, if we set L to the smallest value sufficient for accuracy then we can acquire a veryefficient scheme.However, it is not straightforward to solve these linear equations effectively, therefore we present asspecial algorithm for this, which we call “hybrid splitting solver algorithm”.
5. Hybrid splitting solver algorithm
In this Section, we write out the matrix form of the linear equations resulting from the approximationof the exponential operator e − i ∆ t ( H ρ + H CN ) to familiarize ourselves with its structure, which is required fordeveloping an efficient propagation algorithm for our hybrid splitting scheme. Then we outline the solutionalgorithm. Recalling the notations X i,j = Ψ i,j ( t + ∆ t ) , Ψ i,j = Ψ i,j ( t ) , let us construct the column vectors corre-sponding to the i th row of the 2D problem as Ψ i = (cid:0) Ψ i, Ψ i, Ψ i, . . . Ψ i,N ρ (cid:1) T and X i = (cid:0) X i, X i, X i, . . . X i,N ρ (cid:1) T . (69)Then, the joint problem (64), (65) will take the block pentadiagonal form of B C F . . . B C F . . . A B C F . . . ... . . . . . . . . . . . . . . . ... . . . E N z − A N z − B N z − C N z − F N z − . . . N z − A N z − B N z − C N z − . . . N z A N z B N z X X X ... X N z − X N z − X N z = y y y ... y N z − y N z − y N z . (70)Here E i , A i , B i , C i , F i are ( N ρ + 1) × ( N ρ + 1) matrices. Particularly, E i = (cid:110) diag (0 , e z,i, , . . . , e z,i,L , , . . . , if i ∈ [2 , N z ] , (71) C i = (cid:110) diag (0 , c z,i, , . . . , c z,i,L , , . . . , if i ∈ [1 , N z ] , (72) A i = (cid:110) diag (0 , a z,i, , . . . , a z,i,L , , . . . , if i ∈ [0 , N z − , (73) F i = (cid:110) diag (0 , f z,i, , . . . , f z,i,L , , . . . , if i ∈ [0 , N z − , (74)are the diagonal matrices responsible for coupling the adjacent ρ -rows (with different values of coordinate z ), and B i is an almost five-diagonal matrix in the form of B i = d ,i d ,i d ,i d ,i d ,i . . . a i, b i, c i, f i, . . . e i, a i, b i, c i, f i, . . . ... . . . . . . . . . . . . . . . ... . . . e i,N ρ − a i,N ρ − b i,N ρ − c i,N ρ − f i,N ρ − . . . e i,N ρ − a i,N ρ − b i,N ρ − c i,N ρ − . . . e i,N ρ a i,N ρ b i,N ρ . (75)13n the above matrices we have already taken into account the box, the Neumann and Robin boundaryconditions for the m = 0 configuration. The d ,i , d ,i , d ,i , d ,i , d ,i coefficients are given by the expandedequation (38) as d ,i = −
25 + 6( γ/β )∆ ρ · δ R,i , d ,i = 48 , d ,i = − , d ,i = 16 , d ,i = − (76)for all i ∈ [0 , N z ] . These coefficients take the diagonal form of d j,i = δ j, for the m (cid:54) = 0 states.In the Crank-Nicolson region ( j ≤ L ), the coefficients are given by equation (36): b i,j = 1 − β ρ − β z + αV i,j , a z,i,j = c z,i,j = 16 β z e z,i,j = f z,i,j = − β z , (77) e i,j = ( − /j ) β ρ , f i,j = ( − − /j ) β ρ , (78) a i,j = (16 − /j ) β ρ , c i,j = (16 + 8 /j ) β ρ . (79)In the split region ( j ≥ L + 1 ), we inspect the equation (66), and note that only b i,j get modified as b i,j = 1 − β ρ + αV i,j if i ∈ [ L + 1 , N ρ ] , j ∈ [0 , N z ] (80)and the a i,j , c i,j , e i,j , f i,j coefficients are the same as in (78), (79).The right hand sides are given by (36) and (66), which can be written in a somewhat simpler matrixform of y i = 2 Ψ i − E i Ψ i − − A i Ψ i − − B i Ψ i − C i Ψ i +1 − F i Ψ i +2 with y i, = 0 , (81)which also takes into account the boundary conditions. One can see that the directional splitting introduced N ρ − L zeros at the end of the diagonal of the matrices E i , A i , C i , F i which means that the corresponding ρ -lines are not directly coupled. Taking advantage of thiswe significantly increase the computational efficiency by eliminating the improper matrix elements to reducethe effective block size to ( L + 1) .To proceed, we take the equations corresponding to the i th block matrix row (cid:2) . . . E i A i B i C i F i . . . (cid:3) · X = y i (82)and write out their coefficient matrix from rows L − to N ρ : . . . e i,L − a i,L − b i,L − c i,L − f i,L − . . . . . .. . . e i,L a i,L b i,L c i,L f i,L . . . . . .. . . e i,L +1 a i,L +1 b i,L +1 c i,L +1 f i,L +1 . . . . . .. . . ... ... ... . . . . . . . . . . . . . . . ... . . .. . . . . . e i,N ρ − a i,N ρ − b i,N ρ − c i,N ρ − f i,N ρ − . . .. . . . . . e i,N ρ − a i,N ρ − b i,N ρ − c i,N ρ − . . .. . . . . . e i,N ρ a i,N ρ b i,N ρ . . . . (83)Here we remind that rows j ≤ L have extra nonzero entries far away from the diagonal, cf. (71)-(74).These lines cannot be used during the row operations, but the rest of them, with j > L can be used.To reduce Eq. (70) to a smaller block five-diagonal problem, f i,L − , c i,L , f i,L must be eliminated for all i ∈ [0 , N z ] . Then, the solution in the 2D Crank-Nicolson region G CN will no longer depend on the solutionin the directionally split region G Split . The structure of the B i matrix makes it possible to use the followingbackward elimination process from the N ρ th equation on, for all i ∈ [0 , N z ] :14 . . . e i,L − a i,L − ˜ b i,L − ˜ c i,L − . . . . . .. . . e i,L ˜ a i,L ˜ b i,L . . . . . .. . . e i,L +1 ˜ a i,L +1 ˜ b i,L +1 . . . . . .. . . ... ... ... . . . . . . . . . . . . . . . ... . . .. . . . . . e i,N ρ − ˜ a i,N ρ − ˜ b i,N ρ − . . .. . . . . . e i,N ρ − ˜ a i,N ρ − ˜ b i,N ρ − . . .. . . . . . e i,N ρ ˜ a i,N ρ ˜ b i,N ρ . . . (84)with right hand side of (cid:101) y i = (cid:0) y i, . . . y i,L − ˜ y i,L − ˜ y i,L ˜ y i,L +1 . . . ˜ y i,N ρ − ˜ y i,N ρ (cid:1) T , (85)where ˜ c i,j = (cid:40) c i,j if j = N ρ − c i,j − ( f i,j / ˜ b i,j +2 )˜ a i,j +2 if j = N ρ − . . . L − , (86) ˜ a i,j = (cid:40) a i,j if j = N ρ a i,j − (˜ c i,j / ˜ b i,j +1 ) e i,j +1 if j = N ρ − . . . L, (87) ˜ b i,j = b i,j if j = N ρ b i,j − (˜ c i,j / ˜ b i,j +1 )˜ a i,j +1 if j = N ρ − b i,j − (˜ c i,j / ˜ b i,j +1 )˜ a i,j +1 − ( f i,j / ˜ b i,j +2 ) e i,j +2 if j = N ρ − . . . Lb i,j − ( f i,j / ˜ b i,j +2 ) e i,j +2 if j = L − , (88) ˜ y i,j = y i,j if j = N ρ y i,j − (˜ c i,j / ˜ b i,j +1 )˜ y i,j +1 if j = N ρ − y i,j − (˜ c i,j / ˜ b i,j +1 )˜ y i,j +1 − ( f i,j / ˜ b i,j +2 )˜ y i,j +2 if j = N ρ − . . . Ly i,j − ( f i,j / ˜ b i,j +2 )˜ y i,j +2 if j = L − , (89)and the rest of the values remain unchanged. During the process, ˜ c i,j needs to be calculated first, then ˜ a i,j , ˜ b i,j , ˜ y i,j can be evaluated. Additionally, the value of ˜ c i,j is needed only in one step, then it can be discarded.When the reduced equations are solved (cf. the next Section), we can solve for all variables with forwardsubstitution: X i,j = (cid:40) solution of the reduced block five diagonal part if j ≤ L (˜ y i,j − ˜ a i,j X i,j − − e i,j X i,j − ) / ˜ b i,j if j = L + 1 . . . N ρ . (90)These formulae above can be obtained by disassembling existing five-diagonal solvers (for example [66]),however, care must be taken handling the boundary values at the j = L edge and in the forward substitutionafterwards.Because of the special structure of the coefficient matrix in (70), this process of backward eliminationand forward substitution can be parallelized to N z + 1 independent threads. They are depending on asynchronization step though, which consists of the solution of the following block five-diagonal part. After we performed the elimination procedure for every block B i , we obtain a block five-diagonal systemjust like (70) with a drastically reduced block size, in the following form:15 (cid:101) B (cid:101) C (cid:101) F . . . (cid:101) A (cid:101) B (cid:101) C (cid:101) F . . . (cid:101) E (cid:101) A (cid:101) B (cid:101) C (cid:101) F . . . ... . . . . . . . . . . . . . . . ... . . . (cid:101) E N z − (cid:101) A N z − (cid:101) B N z − (cid:101) C N z − (cid:101) F N z − . . . (cid:101) E N z − (cid:101) A N z − (cid:101) B N z − (cid:101) C N z − . . . (cid:101) E N z (cid:101) A N z (cid:101) B N z (cid:102) X (cid:102) X (cid:102) X ... (cid:102) X N z − (cid:102) X N z − (cid:102) X N z = (cid:101) y (cid:101) y (cid:101) y ... (cid:101) y N z − (cid:101) y N z − (cid:101) y N z . (91)These block matrices read as: (cid:101) E i = (cid:110) diag (0 , e z,i, , . . . , e z,i,L ) if i ∈ [2 , N z ] , (92) (cid:101) C i = (cid:110) diag (0 , c z,i, , . . . , c z,i,L ) if i ∈ [1 , N z ] , (93) (cid:101) A i = (cid:110) diag (0 , a z,i, , . . . , a z,i,L ) if i ∈ [0 , N z − , (94) (cid:101) F i = (cid:110) diag (0 , f z,i, , . . . , f z,i,L ) if i ∈ [0 , N z − , (95) (cid:101) B i = d ,i d ,i d ,i d ,i d ,i . . . a i, b i, c i, f i, . . . e i, a i, b i, c i, f i, . . . ... . . . . . . . . . . . . . . . ... . . . e i,L − a i,L − b i,L − c i,L − f i,L − . . . e i,L − a i,L − b i,L − ˜ c i,L − . . . e i,L ˜ a i,L ˜ b i,L , (96) (cid:101) y i = (cid:0) y i, y i, . . . y i,L − ˜ y i,L − ˜ y i,L (cid:1) T , (97) (cid:102) X i = (cid:0) X i, X i, . . . X i,L − X i,L − X i,L (cid:1) T . (98)The method of acquiring the solution to this system can be chosen freely. For example, it can be solvedby applying Gaussian elimination (with forward substitution) to the ( L + 1) matrix blocks or directly tothe L + 1 wide diagonal matrix. In either case, the time to obtain the solution is drastically reduced (if L (cid:28) N ρ ). After acquiring the solution of (91) we complete the hybrid splitting solver algorithm by formula(90), as indicated in the previous Section.
6. Numerical results
Now let us turn our attention to the numerical test results of the hybrid splitting method. First, inSection 6.1, we investigate the errors related to the spatial discretization. Section 6.2 is devoted to the testsof the time propagation errors. In Section 6.3 we demonstrate the performance of our algorithm through areal-world example.
To investigate the errors related to the spatial discretization, we numerically compute the stationarystates of selected test potentials. We calculate energy errors caused by the finite differences by comparingthe numerical and the exact energy eigenvalues for certain energy eigenstates. We have also set ∆ ρ = ∆ z .In these tests we have used our implementation of the true singular Coulomb potential (C). Another testcase, called “Soft-Coulomb” potential (SC), differs from C only in that the boundary condition (67) at theorigin is replaced by (34) at the origin. The results of both the C and the SC test cases are compared to the16arget stationary statesID Potential Quantum Numbers Energy State Normalization ParametersC-1s − γ/r n = 1 , l = 0 , m = 0 − γ µ/ e − µγr π − / ( µγ ) / µ = 1 , γ = 1 C-2s − γ/r n = 2 , l = 0 , m = 0 − γ µ/ − µγr ] e − µγr/ − / π − / ( µγ ) / µ = 1 , γ = 1 C- z − γ/r n = 2 , l = 1 , m = 0 − γ µ/ ze − µγr/ − / π − / ( µγ ) / µ = 1 , γ = 1 C- x − γ/r n = 2 , l = 1 , m = 1 − γ µ/ ρe − µγr/ cos φ − / π − / ( µγ ) / µ = 1 , γ = 1 H-000 µω r n x = n y = 0 , n z = 0 ω e − µωr / π − / ( µω ) / µ = 1 , ω = 1 H-001 µω r n x = n y = 0 , n z = 1 ω ze − µωr / / π − / ( µω ) / µ = 1 , ω = 1 Table 1: Properties of the potentials and their bound states that we used for testing the errors regarding the spatial discretiza-tion. Our main test potentials are the Coulomb potential (C) and the 3D harmonic oscillator potential (H). exact eigenenergies and eigenstates of s , s , p z . We also compute the eigenenergy of the state p x , whereC and SC are identical (see Section 2.3). The energy errors of the 3D quantum harmonic oscillator (H) arealso tested for the ground and a first exited state. The most important features of our test potentials aresummarized in Table 2. SC - 1sSC - 2sH - 000H - 001C - 1sC - 2sC - 2p z SC - 2p z C - 2p x - - - D z @ a.u. D A b s o l u t ee n e r g y e rr o r @ a . u . D Figure 2: Absolute energy errors of stationary states listed in Table 1, on a log-log scale. The Coulomb eigenstates were testedwith (C) and without (SC) applying the condition (35) at the origin. (In the case of the p z the error curves are within linethickness.) Curves for s and s show that our hybrid splitting algorithm decreases the errors with high order even with thesingular Coulomb potential, in the range . (cid:46) ∆ z (cid:46) . . We computed the eigenenergy values using imaginary time propagation with the second order hybridsplitting scheme outlined in Section 4.3, where the parameter ∆ t was chosen suitably small to minimizethe spatial errors related to the operator splitting. These energy values were also verified with real timepropagation, by comparing the time dependence of the phase factor of the eigenstates with the exact timedependence of Ψ n ( z, ρ, t ) = ψ n ( z, ρ ) e − iE n t : the first two significant digits of the energy error were equal. Itis worth noting that, although we used our hybrid splitting for calculations, these energy errors are relatedonly to the spatial finite differences (that is, the 2D Crank-Nicolson scheme). For these results, it wassufficient to set the “Crank-Nicolson width” L just above ∼ / ∆ z , further increase of L did not improve theresults. 17bsolute energy errors of selected stationary states ∆ z SC-1s . × − . × − . × − . × − . × − . − . SC-2s . × − . × − . × − . × − . × − . − . SC- z . × − . × − . × − . × − . × − . − . C-1s . × − . × − . × − . × − . × − . − . C-2s . × − . × − . × − . × − . × − . − . C- z . × − . × − . × − . × − . × − . − . C- x . × − . × − . × − . × − . × − . − . H-000 . × − . × − . × − . × − . × − .
31 +1 . H-001 . × − . × − . × − . × − . × − .
20 +2 . Table 2: A detailed table of the energy errors at specific values of ∆ z denoted by vertical dashed lines in Figure 2. TheOrder column gives the exponent of ∆ z , i.e. the order of error decrease, valid from step size ∆ z = 0 . to ∆ z = 0 . , whichcharacterizes the effective accuracy of the method in the given case. The Coulomb p z state also has the same error data setto the first 3 significant digits regardless whether or not the Robin boundary condition was used. From Figure 2 and Table 2 we can draw several conclusions. In the case of the 3D harmonic oscillator (areasonably smooth potential), the method is fourth order accurate in ∆ z , as expected (the effective orderis higher than 4 due to the fifth order accurate second derivatives). For the Coulomb potential, thingsget considerably more elaborate. In the Soft-Coulomb case (SC), the energies of the l = 0 states convergewith second order (more precisely, with 1.86), meaning the method will not be sufficiently accurate or fast.However, for the true singular Coulomb case (C), the accuracy of the l = 0 states drastically improve bytwo orders of magnitude around ∆ z = 0 . . On the other hand, the complicated step size dependence of theenergy error, shown by the corresponding lines in Figure 2, does not allow for a true effective order validin the inspected range. However, our algorithm does decrease the error with high order (close to 4) when . (cid:46) ∆ z (cid:46) . , which is anyhow that range where the computation runtime is reasonable. The unusualstep size dependence of the energy error below ∆ z (cid:46) . for data sets C-1s, C-2s is due to finite differencingwith high spatial accuracy applied to a non-analytic problem ( Ψ is not continuously differentiable in theorigin), along with the artificially high order boundary condition.The test cases with l > work reasonably well for the m = 0 configuration both with C and SC: theSC- z or the C- z case has better (relative) accuracy than H-000 but its order is only . as shown inFigure 2. A slightly higher order of . is achieved for C- x (the only test case with m (cid:54) = 0 ) which hasthe best accuracy of all computed eigenenergies in this Section.In conclusion, the numerical error of the hybrid splitting algorithm displays a high order scaling withthe spatial step size for the stationary states of the 3D Hydrogen problem. In the previous section we focused on the spatial errors in conjunction with the singular Coulomb po-tential, now we turn our attention to the total error of the hybrid splitting algorithm using a smoothtime-dependent potential. This total error is composed typically of several terms related to the finite dif-ferences, to the Padé-approximation, to the factorization of the exponential operator, and to the short timesplitting of the evolution operator.We use the so called forced harmonic oscillator (FHO) problem as the test case, defined by the potential V ( z, ρ, t ) = µω ( z + ρ ) + z F sin ω F t , having a known Ψ A ( z, ρ, t ) analytical solution, which we summarizein Appendix C. Here we only illustrate this time-dependent analytical solution in Figure 3.For error calculations, we compare the analytical and time-dependent wavefunctions at a fixed timeinstant by calculating the so called mean-square error, which provides information about both the totalphase and the amplitude related errors. We define the mean-square error by err L = 2 π ˆ ˆ ρ (cid:12)(cid:12) Ψ A ( z, ρ ) − Ψ( z, ρ ) (cid:12)(cid:12) d ρ d z ≈ π (cid:88) i,j ρ j (cid:12)(cid:12) Ψ A i,j − Ψ i,j (cid:12)(cid:12) . (99)18 = w = (cid:144)
20 20 40 60 80 100 - - @ a.u. D S p a ti a l c oo r d i n a t ez @ a . u D Time dependentmean value x H t L = X z \ w t (cid:217) L H t’ L (cid:226) t’10 (cid:137)x(cid:160) H t L @ a.u. D P h a s e m a g n it ud e Phase phactors of the analitical solution
Figure 3: Time evolution of the analytical solution of the forced harmonic oscillator (Appendix C): the expectation value of z versus time (left), time dependence of the phase contributions (right). Parameters are the same as in Table 4, if not givenexplicitly. We choose the mean-square error as the reference error type because (according to our experience) it is themost reliable evolution error quantifier at a fixed time instant t .ID Propagation Method Evaluations of S S2 Second order operator hybrid splitting scheme (Sec. 4.3). 1S4E3 S2 based 4th order iterative splitting scheme using eq. (48) with n = 4 . 3S6E5 S2 based 6th order iterative splitting scheme using eq. (49) with n = 6 . 5S6E9 S2 based 6th order iterative splitting scheme using eq. (48) with n = 6 . 9 Table 3: Test case definitions for comparison of the different high order split operator schemes.
We tested our algorithm with different split-operators (listed in Table 3), and several ∆ z and ∆ t config-urations to acquire a top-down view of the error properties of the hybrid splitting method.The mean-square error of the different high order split-operators based on our Crank-Nicolson scheme(CN5) can be found in Table 4, and also in Figure 4: the standard convergence of the S2 scheme is veryslow, and there is a drastic improvement using S4E3 on top of S2. Using either S6E9 or S6E5 does not yieldmuch gain in accuracy compared to the their higher numerical costs. Based on these, we propose to use thefourth order split-operator method S4E3 in accordance with (48), along with our hybrid splitting scheme. (cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230) (cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224) (cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236) Test Cases (cid:230) S2 (cid:224) S4E3 (cid:236)
S6E50.01 0.02 0.05 0.10 0.20 0.500.020.050.100.200.501.00 Time step @ a.u. D M e a n s qu a r ee rr o r s D z = Dr = (cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230) (cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224) (cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236) Test Cases (cid:230) S2 (cid:224) S4E3 (cid:236)
S6E50.01 0.02 0.05 0.10 0.20 0.500.0010.0050.0100.0500.1000.5001.000 Time step @ a.u. D M e a n s qu a r ee rr o r s D z = Dr = (cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230)(cid:230) (cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224)(cid:224) (cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236)(cid:236) Test Cases (cid:230) S2 (cid:224) S4E3 (cid:236)
S6E50.01 0.02 0.05 0.10 0.20 0.5010 - @ a.u. D M e a n s qu a r ee rr o r s D z = Dr = Figure 4: Mean-square errors as the function of the time step, on a log-log scale, computed with different split-operator methodsand spatial discretization steps, corresponding to Table 4. These plots clearly show the existence of a threshold value of ∆ t ,below which the total error is not reduced anymore. The lines corresponding to the different split-operators in Figure 4 should exhibit the expected powerscaling with ∆ t , this is only approximately the case and only above a threshold ∆ t . Below this threshold thetotal error is not reduced by decreasing ∆ t , because the evolution error is dominated by the finite differences:the error magnitudes can even be predicted as the product of the stationary energy error (Table 2) and thepropagation time interval. 19alculated L error values at ∆ z = 0 . t .
37 4 . × − . × − . × − . × − . × − . × − S4E3 . × − . × − . × − . × − . × − . × − . × − S6E5 . × − . × − . × − . × − . × − . × − . × − S6E9 . × − . × − . × − . × − . × − . × − . × − Calculated L error values at ∆ z = 0 . t . × − . × − . × − . × − . × − . × − S4E3 . × − . × − . × − . × − . × − . × − . × − S6E5 . × − . × − . × − . × − . × − . × − . × − S6E9 . × − . × − . × − . × − . × − . × − . × − Calculated L error values at ∆ z = 0 . t .
38 7 . × − . × − . × − . × − . × − . × − S4E3 . × − . × − . × − . × − . × − . × − . × − S6E5 . × − . × − . × − . × − . × − . × − . × − S6E9 . × − . × − . × − . × − . × − . × − . × − Table 4: Mean-square errors of the forced harmonic oscillator with different ∆ z (= ∆ ρ ) , ∆ t and split-operator configurations.All calculations were carried out in a box of − ≤ z ≤ and ≤ ρ ≤ with propagation parameters ω F = 2 π/ , F = 1 , µ = 1 , ω = 1 , launched from the corresponding ground state (H-000). All error values are calculated at the time t = 100 . Calculated L error values of CN3 based S4E3 method ∆ t ∆ z = 0 . . × − . × − . × − . × − . × − . × − . × − ∆ z = 0 . . × − . × − . × − . × − . × − . × − . × − ∆ z = 0 .
05 5 . × − . × − . × − . × − . × − . × − . × − Table 5: Mean-square errors of a commonly used propagation scheme, to be compared with the data of Table 4. All theparameters are the same as in Table 4. Comparison of the last columns shows one, two and three orders of magnitude accuracyincrease in favor of the five point discretization, as ∆ z is decreased. One objective of our research was to design a real space finite difference algorithm that achieves fourthorder accuracy in the spatial step size and, most importantly, that is capable to include both the singularcoordinate ρ = 0 and the singular Coulomb potential directly. Therefore, it is interesting to globallybenchmark the results against a second order 3-point Crank-Nicolson finite difference scheme (CN3), whichmeans we degraded the core algorithm to use three-point finite differences in S2.We repeated the main tests with this CN3 method combined with the split-operator configuration S4E3(discussed above). The CN3 results are shown in Table 5 and are to be compared to the results with CN5 inTable 4. Both the CN3 and CN5 schemes display the expected order scaling of the error with ∆ z . Comparingthe L error at ∆ t = 0 . (which is already below the threshold ∆ t ), our fourth order discretization hasone, two and three orders of magnitudes smaller errors at spatial step sizes ∆ z = 0 . , . , . , respectively.Thus, one should always use (at least) the CN5 scheme unless the exotic nature of the problem preventshigh order discretization.In conclusion, our hybrid splitting propagation method is suitable for numerical simulations with time-dependent Hamiltonians, and it is capable of high order scaling with ∆ t , once it is incorporated into theproper high order evolution operator formulae. 20 .3. Hydrogen atom in an external electric field Finally, we conducted a test similar to real world applications by simulating the hydrogen atom in anexternal laser field. Now, we compare our hybrid splitting based CN5 implementation, including Coulombpotential and the Coulomb boundary condition (CN5-C) against the commonplace CN3 discretization withthe best Soft-Coulomb potential approximation allowed by the spatial grid (CN3-SC), which is the sameapproximation that we tested in Section 6.1. We did not use the full 2D CN3 method, but we applied thehybrid splitting scheme using the CN3 equations just like the case of the FHO.The external field was parametrized as f ( t ) = F sin ωt . The simulated time was t ∈ [0 , , and itconsisted of only one field cycle, by setting ω = 2 π/ in atomic units. The field amplitude was F = 0 . ,with simulation box size z ∈ [ − , , ρ ∈ [0 , to contain the escaping electron waves. The initialstate was the s ground state, found by imaginary time propagation method to remove spurious ionizationcomponents. We use the S4E3 high order split operator formula as indicated in the last section. (cid:45)
20 0 20 40 60 (cid:45) (cid:45) (cid:64) a.u (cid:68) Ρ (cid:64) a . u . (cid:68) Snapshot of Log (cid:200) (cid:89) (cid:72) Ρ ,z (cid:76) at time (cid:61)
100 a.u. (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) Figure 5: Density plot of the logarithm of the absolute square of the wave function, in a plane containing the z axis, at the endof the simulation described in Section 6.3. Note that the white waves are 4-5 orders of magnitude smaller than the sphericalpeak. We show the result of this simulation in Figure 5 by a density plot of the logarithm of the absolute squareof the wave function, in a plane containing the z axis. The white waves on the right of the spherical peak,bound by the Coulomb potential of the nucleus, are ca. 4-5 order smaller in density. These waves have tobe computed very accurately, since they contribute the most to the time dependence of the dipole moment,which in turn is of fundamental importance regarding the HHG and the creation of attosecond light pulses.The tests to compare CN5-C with CN3-SC were run with ∆ ρ = ∆ z and with time step ∆ t = 0 . if notindicated explicitly. We quantified the accuracy of the solutions with the error of the expectation value (cid:104) z (cid:105) ,i.e. the magnitude of the time-dependent dipole moment. Unlike the FHO example, analytic solution is notavailable this time, so we use a converged solution to determine numerical errors, that is we compare theresults with a much more accurate numerical solution obtained using smaller ∆ t , ∆ z discretization steps.The results, shown in Figure 6 are as expected: the error of the CN5-C scheme with ∆ z = 0 . is smallerby two orders of magnitude than that of the CN3-SC scheme with ∆ z = 0 . . Based on this, we estimate theperformance difference as follows. Due to the second order convergence of the CN3-SC scheme, ∆ z ≈ . is needed to achieve the accuracy of the CN5-C with ∆ z = 0 . , which means a factor of in the numberof spatial gridpoints, implying a factor of - in runtime. That is, the CN5-C scheme is ca. 1000times more effective than the CN3-SC. Regarding absolute accuracy, the magnitude of the error of (cid:104) z (cid:105) usingCN5-C with ∆ z = 0 . is smaller than the L error of the FHO at the time instant t = 100 , as can be seenin Table 4.Regarding the pointwise convergence of the solution with the CN5-C method, we compare the densityand the phase along a z -line section at ρ = 1 , computed with two different values of ∆ z . These plots, shown21 N3 - SC;with D z = - C;with D z = @ a.u. D E rr o r s o f t h e d i po l e m o m e n t @ a . u . D Figure 6: Comparison of the errors of (cid:104) z (cid:105) , during the time propagation described in Sec. 6.3, using the method CN5-C with ∆ z = 0 . a.u. (blue line) and using the method CN3-SC with ∆ z = 0 . a.u. (purple line). Despite the larger spatial stepused with CN5-C, its error is two orders of magnitude smaller, therefore it is shown here with a 100 times magnification. in Figure 7, convincingly show that a converged numerical solution is obtained already at ∆ z = 0 . . Notethat the accuracy of the phase is crucial in strong field physics, regarding e.g. exit momentum calculations[13, 67] in optical tunneling or regarding phase space methods [68–73] (where usually the Wigner functionis computed from the wave function and this is then further analysed for the physical interpretation of theprocess). D z = D z = z - - -
50 0 50 100 15010 - - - - @ a.u. D P r ob a b ilit y d e n s it y @ a . u . D Snapshot of (cid:160) Y H L⁄ at time =
100 a.u. - - - @ a.u. D A r gu m e n t Arg Y H L< at time = @ a.u. D Figure 7: Probability density and phase (in the inset) along a z-line section at ρ = 1 , at the end of the time propagationdescribed in Sec. 6.3, computed using the CN5-C method with two different values of ∆ z , as indicated in the figure. Theexcellent fit of the two curves shows that the solution can be considered converged already at ∆ z = 0 . .
7. Summary
In this article we presented an algorithm capable of the direct numerical solution of the three dimensionaltime dependent Schrödinger equation, assuming axial symmetry in cylindrical coordinates. The main featureof the algorithm is that it is capable to accurately handle singular Coulomb potentials besides any smoothpotential of the form V ( z, ρ ) . The axial symmetry enables a two dimensional grid and the availability ofthe Cartesian z -axis makes it easy to investigate reduced electron dynamics.22e choose a high order finite difference representation in the spatial coordinate domain. We implementedall singularities that can be reduced to the form /ρ or /r via Neumann and Robin boundary conditionsat the ρ = 0 axis. The accuracy of the algorithm is fourth order in ∆ z , ∆ ρ for smooth potentials and it isclose to fourth order for Coulomb potentials in a restricted discretization parameter range. We completedthis algorithm by a high order scalar product formula designed for this finite difference representation.We based our algorithm on the split-operator approximation of the evolution operator. Due to thenon-separability of the Coulomb-problem in cylindrical coordinates, we constructed a special second orderoperator splitting scheme called hybrid splitting method. This splits the Hamiltonian matrix direction-wiselike the traditional methods, but the innermost region near ρ = 0 is excluded: here the full 2D Crank-Nicolson equations are used. This means that there are many decoupled 1D Schrödinger equations in the z direction, and a 2D Schrödinger equation with a special block pentadiagonal pattern of the coefficient matrix,which has to be taken advantage of for maximum efficiency, like in our hybrid splitting solver algorithm.Thread based parallelization is supported throughout, and we also gave a way to evaluate the decoupled 1DSchrödinger equations with high spatial accuracy and efficiency.We thoroughly investigated the spatial discretization related errors in an optimal discretization parameterrange, determining detailed accuracy characteristics with or without Coulomb potential. We also verifiedthe considerably increased accuracy for numerical simulations forced oscillator. Testing the performance ofthe 4th and the 6th order (iterative) split-operator factorizations built from second order hybrid splittingmethod, we observed a threshold ∆ t value below which there is no accuracy gain. We concluded that highorder (meaning at least 4th order) split-operator formula should be used in practice, accompanying thehybrid splitting algorithm.In order to demonstrate the accuracy and performance of our hybrid splitting algorithm also in a sim-ulation close to the planned applications, we computed the solution for a hydrogen atom in a strong time-dependent electric field of one sine period. The important waves in the probability density, having anamplitude of − - − relative to the peak value, were obtained accurately and efficiently.The hybrid splitting scheme, with some minor modifications of the algorithm, is also capable to handlesingle or multiple coupled nonlinear time dependent Schrödinger equations, such as those arising e.g. intime-dependent density functional theory [74], time-dependent Hartree-Fock methods [75, 76] or severalother areas of physics. Acknowledgements
The authors thank W. Becker, M. G. Benedict, P. Földi, K. Varjú and S. Varró for stimulating discussions.This research has been granted by the Hungarian Scientific Research Fund OTKA under contract No.T81364. Partial support by the ELI-ALPS project is also acknowledged. The ELI-ALPS project (GOP-1.1.1-12/B-2012-000, GINOP-2.3.6-15-2015-00001) is supported by the European Union and co-financed bythe European Regional Development Fund.
References [1] M. Hentschel, R. Kienberger, C. Spielmann, G. A. Reider, N. Milosevic, T. Brabec, P. Corkum, U. Heinzmann, M. Drescher,F. Krausz, Attosecond metrology, Nature 414 (6863) (2001) 509–513.[2] P. . M. Paul, E. Toma, P. Breger, G. Mullot, F. Augé, P. Balcou, H. Muller, P. Agostini, Observation of a train ofattosecond pulses from high harmonic generation, Science 292 (5522) (2001) 1689–1692.[3] R. Kienberger, M. Hentschel, M. Uiberacker, C. Spielmann, M. Kitzler, A. Scrinzi, M. Wieland, T. Westerwalbesloh,U. Kleineberg, U. Heinzmann, et al., Steering attosecond electron wave packets with light, Science 297 (5584) (2002)1144–1148.[4] M. Drescher, M. Hentschel, R. Kienberger, M. Uiberacker, V. Yakovlev, A. Scrinzi, T. Westerwalbesloh, U. Kleineberg,U. Heinzmann, F. Krausz, Time-resolved atomic inner-shell spectroscopy, Nature 419 (6909) (2002) 803–807.[5] A. Baltuška, T. Udem, M. Uiberacker, M. Hentschel, E. Goulielmakis, C. Gohle, R. Holzwarth, V. Yakovlev, A. Scrinzi,T. Hänsch, et al., Attosecond control of electronic processes by intense light fields, Nature 421 (6923) (2003) 611–615.[6] F. Krausz, M. Ivanov, Attosecond physics, Reviews of Modern Physics 81 (1) (2009) 163–234.
7] A. McPherson, G. Gibson, H. Jara, U. Johann, T. S. Luk, I. McIntyre, K. Boyer, C. K. Rhodes, Studies of multiphotonproduction of vacuum-ultraviolet radiation in the rare gases, Journal of the Optical Society of America B 4 (4) (1987)595–601.[8] M. Ferray, A. L’Huillier, X. Li, L. Lompre, G. Mainfray, C. Manus, Multiple-harmonic conversion of 1064 nm radiation inrare gases, Journal of Physics B: Atomic, Molecular and Optical Physics 21 (3) (1988) L31.[9] G. Farkas, C. Tóth, Proposal for attosecond light pulse generation using laser induced multiple-harmonic conversionprocesses in rare gases, Physics Letters A 168 (5) (1992) 447–450.[10] S. Harris, J. Macklin, T. Hänsch, Atomic scale temporal structure inherent to high-order harmonic generation, Opticscommunications 100 (5-6) (1993) 487–490.[11] M. Uiberacker, T. Uphues, M. Schultze, A. J. Verhoef, V. Yakovlev, M. F. Kling, J. Rauschenberger, N. M. Kabachnik,H. Schröder, M. Lezius, et al., Attosecond real-time observation of electron tunnelling in atoms, Nature 446 (7136) (2007)627–632.[12] S. Haessler, J. Caillat, W. Boutu, C. Giovanetti-Teixeira, T. Ruchon, T. Auguste, Z. Diveki, P. Breger, A. Maquet,B. Carré, et al., Attosecond imaging of molecular electronic wavepackets, Nature Physics 6 (3) (2010) 200–206.[13] A. N. Pfeiffer, C. Cirelli, A. S. Landsman, M. Smolarski, D. Dimitrovski, L. B. Madsen, U. Keller, Probing the longitudinalmomentum spread of the electron wave packet at the tunnel exit, Physical Review Letters 109 (8) (2012) 083002.[14] D. Shafir, H. Soifer, B. D. Bruner, M. Dagan, Y. Mairesse, S. Patchkovskii, M. Y. Ivanov, O. Smirnova, N. Dudovich,Resolving the time when an electron exits a tunnelling barrier, Nature 485 (7398) (2012) 343–346.[15] M. Schultze, M. Fieß, N. Karpowicz, J. Gagnon, M. Korbman, M. Hofstetter, S. Neppl, A. L. Cavalieri, Y. Komninos,T. Mercouris, et al., Delay in photoemission, Science 328 (5986) (2010) 1658–1662.[16] A. L. Cavalieri, N. Müller, T. Uphues, V. S. Yakovlev, A. Baltuška, B. Horvath, B. Schmidt, L. Blümel, R. Holzwarth,S. Hendel, et al., Attosecond spectroscopy in condensed matter, Nature 449 (7165) (2007) 1029–1032.[17] A. D. Bandrauk, M. Ivanov, Quantum Dynamic Imaging: Theoretical and Numerical Methods, Springer Science & BusinessMedia, 2011.[18] L. Keldysh, Ionization in the field of a strong electromagnetic wave, Soviet Physics JETP 20 (5) (1965) 1307–1314.[19] P. B. Corkum, Plasma perspective on strong field multiphoton ionization, Physical Review Letters 71 (13) (1993) 1994.[20] P. Eckle, A. Pfeiffer, C. Cirelli, A. Staudte, R. Dörner, H. Muller, M. Büttiker, U. Keller, Attosecond ionization andtunneling delay time measurements in helium, Science 322 (5907) (2008) 1525–1529.[21] M. Lewenstein, P. Balcou, M. Y. Ivanov, A. L’Huillier, Auillier, P. B. Corkum, Theory of high-harmonic generation bylow-frequency laser fields, Physical Review A 49 (3) (1994) 2117.[22] M. Y. Ivanov, M. Spanner, O. Smirnova, Anatomy of strong field ionization, Journal of Modern Optics 52 (2-3) (2005)165–184.[23] A. Gordon, R. Santra, F. X. Kärtner, Role of the coulomb singularity in high-order harmonic generation, Physical ReviewA 72 (6) (2005) 063411.[24] X.-M. Tong, S.-I. Chu, Theoretical study of multiple high-order harmonic generation by intense ultrashort pulsed laserfields: A new generalized pseudospectral time-dependent method, Chemical Physics 217 (2) (1997) 119–130.[25] H. Bachau, E. Cormier, P. Decleva, J. Hansen, F. Martin, Applications of b-splines in atomic and molecular physics,Reports on Progress in Physics 64 (12) (2001) 1815.[26] E. Cormier, P. Lambropoulos, Above-threshold ionization spectrum of hydrogen using b-spline functions, Journal of PhysicsB: Atomic, Molecular and Optical Physics 30 (1) (1997) 77.[27] H. Muller, An efficient propagation scheme for the time-dependent Schrödinger equation in the velocity gauge, LaserPhysics 9 (1999) 138–148.[28] D. Bauer, P. Koval, Qprop: A Schrödinger-solver for intense laser–atom interaction, Computer Physics Communications174 (5) (2006) 396–421.[29] E. S. Smyth, J. S. Parker, K. T. Taylor, Numerical integration of the time-dependent Schrödinger equation for laser-drivenhelium, Computer Physics Communications 114 (1) (1998) 1–14.[30] L. Argenti, R. Pazourek, J. Feist, S. Nagele, M. Liertzer, E. Persson, J. Burgdörfer, E. Lindroth, Photoionization of heliumby attosecond pulses: Extraction of spectra from correlated wave functions, Physical Review A 87 (5) (2013) 053405.[31] B. I. Schneider, L. A. Collins, S. Hu, Parallel solver for the time-dependent linear and nonlinear Schrödinger equation,Physical Review E 73 (3) (2006) 036708.[32] A. Kołakowska, M. Pindzola, F. Robicheaux, D. Schultz, J. Wells, Excitation and charge transfer in proton-hydrogencollisions, Physical Review A 58 (4) (1998) 2872.[33] I. P. Christov, M. M. Murnane, H. C. Kapteyn, High-harmonic generation of attosecond pulses in the single-cycle regime,Physical Review Letters 78 (7) (1997) 1251.[34] S. Hu, L. Collins, Intense laser-induced recombination: The inverse above-threshold ionization process, Physical ReviewA 70 (1) (2004) 013407.[35] J. Wells, D. Schultz, P. Gavras, M. Pindzola, Numerical solution of the time-dependent Schrödinger equation forintermediate-energy collisions of antiprotons with hydrogen, Physical Review A 54 (1) (1996) 593.[36] A. Gordon, C. Jirauschek, F. X. Kärtner, Numerical solver of the time-dependent Schrödinger equation with coulombsingularities, Physical Review A 73 (4) (2006) 042505.[37] I. Gainullin, M. Sonkin, High-performance parallel solver for 3d time-dependent Schrödinger equation for large-scalenanosystems, Computer Physics Communications 188 (2015) 68–75.[38] C. Leforestier, R. Bisseling, C. Cerjan, M. Feit, R. Friesner, A. Guldberg, A. Hammerich, G. Jolicard, W. Karrlein, H.-D.Meyer, et al., A comparison of different propagation schemes for the time dependent Schrödinger equation, Journal ofComputational Physics 94 (1) (1991) 59–80.
39] T. Iitaka, Solving the time-dependent Schrödinger equation numerically, Physical Review E 49 (5) (1994) 4684.[40] A. Castro, M. A. Marques, A. Rubio, Propagators for the time-dependent Kohn–Sham equations, The Journal of chemicalphysics 121 (8) (2004) 3425–3433.[41] S. Blanes, P. Moan, Splitting methods for the time-dependent Schrödinger equation, Physics Letters A 265 (1) (2000)35–42.[42] S. A. Chin, C.-R. Chen, Fourth order gradient symplectic integrator methods for solving the time-dependent Schrödingerequation, The Journal of Chemical Physics 114 (17) (2001) 7338–7341.[43] W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery, Numerical Recipes 3rd Edition: The Art of ScientificComputing, 3rd Edition, Cambridge University Press, 2007.[44] I. Puzynin, A. Selin, S. Vinitsky, A high-order accuracy method for numerical solving of the time-dependent Schrödingerequation, Computer Physics Communications 123 (1) (1999) 1–6.[45] H. Wang, Numerical studies on the split-step finite difference method for nonlinear Schrödinger equations, Applied Math-ematics and Computation 170 (1) (2005) 17–35.[46] J. Shen, E. Wei, Z. Huang, M. Chen, X. Wu, High-order symplectic fdtd scheme for solving a time-dependent Schrödingerequation, Computer Physics Communications 184 (3) (2013) 480–492.[47] B. H. Bransden, C. J. Joachain, Physics of Atoms and Molecules, Pearson Education India, 2003.[48] C. Cohen-Tannoudji, B. Diu, F. Laloë, Quantum Mechanics, Quantum Mechanics, Wiley, 1977.[49] D. J. Griffiths, Introduction to Quantum Mechanics, Pearson Education India, 2005.[50] I. Puzynin, A. Selin, S. Vinitsky, Magnus-factorized method for numerical solving the time-dependent Schrödinger equa-tion, Computer Physics Communications 126 (1) (2000) 158–161.[51] S. Blanes, F. Casas, J. Oteo, J. Ros, The magnus expansion and some of its applications, Physics Reports 470 (5) (2009)151–238.[52] W. van Dijk, F. Toyama, Accurate numerical solutions of the time-dependent Schrödinger equation, Physical Review E75 (3) (2007) 36707.[53] J. M. Varah, On the solution of block-tridiagonal systems arising from certain finite-difference equations, Mathematics ofComputation 26 (120) (1972) 859–868.[54] R. Wilcox, Exponential operators and parameter differentiation in quantum physics, Journal of Mathematical Physics8 (4) (1967) 962–982.[55] R. I. McLachlan, G. R. W. Quispel, Splitting methods, Acta Numerica 11 (2002) 341–434.[56] Y. Yazici, Operator splitting methods for differential equations, Ph.D. thesis, Izmir Institute of Technology (2010).[57] G. Strang, On the construction and comparison of difference schemes, SIAM Journal on Numerical Analysis 5 (3) (1968)506–517.[58] S. A. Chin, Symplectic integrators from composite operator factorizations, Physics Letters A 226 (6) (1997) 344–348.[59] M. Suzuki, Hybrid exponential product formulas for unbounded operators with possible applications to monte carlosimulations, Physics Letters A 201 (5) (1995) 425–428.[60] A. D. Bandrauk, H. Shen, Higher order exponential split operator method for solving time-dependent Schrödinger equa-tions, Canadian Journal of Chemistry 70 (2) (1992) 555–559.[61] A. D. Bandrauk, H. Lu, Exponential propagators (integrators) for the time-dependent Schrödinger equation, Journal ofTheoretical and Computational Chemistry 12 (06) (2013) 1340001.[62] A. D. Bandrauk, H. Shen, Exponential split operator methods for solving coupled time-dependent Schrödinger equations,The Journal of chemical physics 99 (2) (1993) 1185–1193.[63] S. A. Chin, S. Janecek, E. Krotscheck, Any order imaginary time propagation method for solving the Schrödinger equation,Chemical Physics Letters 470 (4) (2009) 342–346.[64] L. Lehtovaara, J. Toivanen, J. Eloranta, Solution of time-independent Schrödinger equation by the imaginary time prop-agation method, Journal of Computational Physics 221 (1) (2007) 148–157.[65] J. G. Muga, J. Palao, B. Navarro, I. Egusquiza, Complex absorbing potentials, Physics Reports 395 (6) (2004) 357–426.[66] A. Karawia, A computational algorithm for solving periodic penta-diagonal linear systems, Applied Mathematics andComputation 174 (1) (2006) 613–618.[67] X. Sun, M. Li, J. Yu, Y. Deng, Q. Gong, Y. Liu, Calibration of the initial longitudinal momentum spread of tunnelingionization, Physical Review A 89 (4) (2014) 045402.[68] A. Czirják, R. Kopold, W. Becker, M. Kleber, W. Schleich, The Wigner function for tunneling in a uniform static electricfield, Optics communications 179 (1) (2000) 29–38.[69] L. Guo, S. Han, J. Chen, Time-energy analysis of above-threshold ionization in few-cycle laser pulses, Physical Review A86 (5) (2012) 053409.[70] S. Gräfe, J. Doose, J. Burgdörfer, Quantum phase-space analysis of electronic rescattering dynamics in intense few-cyclelaser fields, Journal of Physics B: Atomic, Molecular and Optical Physics 45 (5) (2012) 055002.[71] A. Czirják, S. Majorosi, J. Kovács, M. G. Benedict, Emergence of oscillations in quantum entanglement during rescattering,Physica Scripta 2013 (T153) (2013) 014013.[72] C. Zagoya, J. Wu, M. Ronto, D. Shalashilin, C. F. de Morisson Faria, Quantum and semiclassical phase-space dynamicsof a wave packet in strong fields using initial-value representations, New Journal of Physics 16 (10) (2014) 103040.[73] C. Baumann, H.-J. Kull, G. Fraiman, Wigner representation of ionization and scattering in strong laser fields, PhysicalReview A 92 (6) (2015) 063420.[74] C. A. Ullrich, Z.-h. Yang, A brief compendium of time-dependent density functional theory, Brazilian Journal of Physics44 (1) (2014) 154–188.[75] K. C. Kulander, Time-dependent Hartree-Fock theory of multiphoton ionization: Helium, Physical Review A 36 (6) (1987) Appendix A. Approximating the inner product
In this Section we propose a solution to the problem of the discrete inner product formula exposed inSection 3.2.We seek a discretized representation of the inner product formula in the cylindrical coordinate system as (cid:104) Φ | Ψ (cid:105) = 2 π ˆ + ∞−∞ ˆ ∞ ρ Φ ∗ ( z, ρ )Ψ( z, ρ ) d ρ d z = (cid:88) i,j c i,j Φ ∗ i,j Ψ i,j . (A.1)The naive approach with coefficients c i,j = 2 πρ j ∆ ρ ∆ z causes inaccuracy which originates from the ρ = 0 edge and its neighborhood only, because the formula with this particular c i,j has exponential convergenceat the box boundaries [43].Besides the accuracy, the conservation of a discretized scheme of form (A.1) is also an issue: givena discrete Hamiltonian matrix H in the Padé-approximation then the time evolution will be unitary withrespect to the inner product (cid:80) i,j c i,j Φ ∗ i,j Ψ i,j only if H is self-adjoint. Unfortunately, because of the Neumannand Robin boundary conditions (38) the Hamiltonian matrix H does not exist on the ρ = 0 line, thereforethere is no standard norm of form (cid:80) i,j c i,j Ψ ∗ i,j Ψ i,j which is perfectly conserved by the numerical scheme.Therefore, our aim is to approximate (A.1) with an order that is higher than the order of the finitedifference scheme. To achieve this goal, we used Lagrange [43] interpolating polynomial p i,j ( ρ ) defined onpoints ( ρ j , Q i,j ) , ( ρ j +1 , Q i,j +1 ) , ( ρ j +2 , Q i,j +2 ) , ( ρ j +3 , Q i,j +3 ) , ( ρ j +4 , Q i,j +4 ) , ( ρ j +5 , Q i,j +5 ) , ( ρ j +6 , Q i,j +6 ) for a given z i line, where Q i,j is the integrand in (cid:104) Φ | Ψ (cid:105) = ˜ Q ( z, ρ ) d ρ d z . We use an elementary integralformula between ρ j and ρ j +1 : ρ j +1 ˆ ρ j Q ( z, ρ ) d ρ = (cid:20) Q i,j + 27132520 Q i,j +1 − Q i,j +2 + 586945 Q i,j +3 − Q i,j +4 Q i,j +5 − Q i,j +6 (cid:21) · ∆ ρ + O (cid:0) ∆ ρ (cid:1) . (A.2)We sum up this for all i, j points, and utilize the boundary conditions for j ≥ N ρ , then we arrive at thefollowing integral formula, which is our choice to approximate the scalar product (A.1): (cid:104) Φ | Ψ (cid:105) = N z (cid:88) i =0 (cid:20) ρ Φ ∗ i, Ψ i, + 8419960480 ρ Φ ∗ i, Ψ i, + 1886930240 ρ Φ ∗ i, Ψ i, + 3762130240 ρ Φ ∗ i, Ψ i, +5503160480 ρ Φ ∗ i, Ψ i, + 6134360480 ρ Φ ∗ i, Ψ i, + N ρ (cid:88) j =6 ρ j Φ ∗ i,j Ψ i,j · π ∆ ρ ∆ z. (A.3)This achieves the high integration accuracy (it is exact for a polynomial of ρ up to degree 6), which isneeded: the computed norm variations become proportional to ∆ ρ , which is consistent with the accuracyof the spatial finite differences. The constructed Crank-Nicolson scheme stayed stable in our simulations.26 ppendix B. Numerov z-line propagator algorithm We have constructed an efficient way of evaluating e − i ∆ t H z line-by-line, based on the second order Crank-Nicolson algorithm with Numerov-extension [77], which also provides at least fourth order accuracy in ∆ z ,and it reduces the numerical costs because it uses three point finite differences and tridiagonal equationsinstead of five point differences and pentadiagonal equations. We will also outline an optimization techniquefor 1D Crank-Nicolson schemes which is makes the computations even more efficient when they are appliedto 2D problems.Let us fix the value of coordinate ρ (i.e we choose j = const) and we denote Ψ i ( t k ) = Ψ i,j ( t k ) . Here weconsider the Hamiltonian in the form of H z = β∂ z + V ( z, t ) along this line, since the procedure below allowsto include a potential of this form, which is however not present in our hybrid splitting scheme. Again, westart with the second order approximation of the exponential operator e − i ∆ t H z as in (27) (cid:0) αβ∂ z + αV (cid:1) Ψ( t k +1 ) = (cid:0) − αβ∂ z − αV (cid:1) Ψ( t k ) (B.1)with α = i ∆ t/ , V = V ( z, t k +1 / ) . We will use a discretized Laplacian L z based on the standard threepoint finite difference method [43], however, we will also need the leading term of the error: L z Ψ i = Ψ i − − i + Ψ i +1 ∆ z ≈ ∂ Ψ ∂z (cid:12)(cid:12)(cid:12)(cid:12) z i − ∆ z ∂ Ψ ∂z (cid:12)(cid:12)(cid:12)(cid:12) z i + O (∆ z ) . (B.2)To proceed, we introduce the auxiliary variable Y ( t k ) = Ψ( t k +1 ) + Ψ( t k ) and rewrite the equation (B.1) as (cid:0) αβ∂ z (cid:1) Y ( t k ) = 2Ψ( t k ) − (1 + αV ) Y ( t k ) . (B.3)Then, we discretize the equation (B.3) using (B.2): ( αβL z ) Y i ( t k ) = 2Ψ i ( t k ) − (1 + αV i ) Y i ( t k ) + αβ ∆ z ∂ Y∂z (cid:12)(cid:12)(cid:12)(cid:12) i (B.4)Now, we evaluate the error term with ∂ z Y i ( t k ) = L z ( ∂ z Y i ( t k )) + O (∆ z ) also using the right hand side of(B.3) to get to the result of the form (cid:18) αV i + αβL z + ∆ z · L z (1 + αV i ) (cid:19) Y i ( t k ) = 2Ψ i ( t k ) + ∆ z L z Ψ i ( t k ) . (B.5)Restoring the equations for Ψ i ( t k +1 ) we arrive at the Numerov-extended Crank-Nicolson algorithm as (cid:18) αV i + αβL z + ∆ z L z (1 + αV i ) (cid:19) Ψ i ( t k +1 ) = (cid:18) − αV i − αβL z + ∆ z L z (1 − αV i ) (cid:19) Ψ i ( t k ) . (B.6)It is interesting to note that by reverse engineering from the Padé-form, we get a discrete Hamiltonian ofthe form (cid:101) H z,i = V i + βL z − i ∆ z t L z + ∆ z L z V i . (B.7)The extra terms in (B.7) improve the solution to fifth order in ∆ z , however, one should note that (i) L z ( V i Ψ i ( t k )) is part of (cid:101) H z Ψ( t k ) , requiring the potential to be continuously differentiable, (ii) (cid:101) H z is nolonger strictly self-adjoint.In (B.6) we have arrived at a tridiagonal system of linear equations of the form b c · · · a b c · · · a b c · · · ... ... . . . . . . . . . ... · · · a N z − b N z − c N z − · · · a N z b N z Ψ ( t k +1 )Ψ ( t k +1 )Ψ ( t k +1 ) ... Ψ N z − ( t k +1 )Ψ N z ( t k +1 ) = y y y ... y N z − y N z (B.8)27hich after forward Gaussian-elimination reads ˜ b c · · ·
00 ˜ b c · · ·
00 0 ˜ b c · · · ... ... . . . . . . . . . ... · · · b N z − c N z − · · · b N z Ψ ( t k +1 )Ψ ( t k +1 )Ψ ( t k +1 ) ... Ψ N z − ( t k +1 )Ψ N z ( t k +1 ) = ˜ y ˜ y ˜ y ... ˜ y N z − ˜ y N z . (B.9)In the coefficient matrix of (B.8) and (B.9) we denoted a i = (cid:110) αV i − + 12 αβ/ ∆ z if i = 1 , . . . , N z , (B.10) c i = (cid:110) αV i +1 + 12 αβ/ ∆ z if i = 0 , . . . , N z − , (B.11) b i = (cid:110)
10 + 10 αV i − αβ/ ∆ z if i = 0 , . . . , N z , (B.12) ˜ b i = (cid:40) b i if i = 0 ,b i − ( a i / ˜ b i − ) c i − if i = 1 , . . . , N z . (B.13)The right hand side and the solution of (B.9) are given by the following expressions: ˜ y i = (20 − b i )Ψ + (2 − c i )Ψ if i = 0 , (2 − a i )Ψ i − + (20 − b i )Ψ i + (2 − c i )Ψ i +1 − ( a i / ˜ b i − )˜ y i − if i = 1 , . . . , N z − , (2 − a i )Ψ N z − + (20 − b i )Ψ N z − ( a i / ˜ b i − )˜ y N z − if i = N z , (B.14) Ψ i ( t k +1 ) = (cid:40) ˜ y N z / ˜ b N z if i = N z , (˜ y i − c i Ψ i +1 ( t k +1 )) / ˜ b i if i = N z − , . . . , . (B.15)This completes the 1D propagation method.Now let us return to our original 2D problem. Because the discrete Hamiltonian H z is independent fromthe value of ρ , the corresponding coefficient matrix of (B.9) is the same for each j -line. This means that weneed to do the forward elimination (B.13) only once, then it is sufficient to perform only the forward (B.14)and the backward substitution (B.15) steps for each j -line in order to acquire the solution. This yields afactor of 2 speedup in the evaluation of e − i ∆ t H z with three point finite differences, if we have many linesto propagate and cache the appropriate variables. This latter can also be viewed as an LU decompositionbased optimization [43]. Appendix C. Analytical solution of the forced harmonic oscillator
The TDSE for the axially symmetric (i.e. m = 0 ) 3D forced harmonic oscillator (FHO) in cylindricalcoordinates is the following: i ∂∂t Ψ( z, ρ, t ) = (cid:20) β ∂ ∂z + β ∂ ∂ρ + βρ ∂∂ρ + 12 µω ( z + ρ ) + f ( t ) z (cid:21) Ψ( z, ρ, t ) (C.1)where β = − / µ . Using the separability in z and ρ coordinates, we can reduce this problem to a onedimensional time-dependent one, which was solved by K. Husimi et. al. [78]. Then, the analytical time-dependent wave function is of the form Ψ A ( z, ρ, t ) = χ ( z − ξ ( t ) , ρ, t ) exp (cid:18) iµ ( z − ξ ( t )) ˙ ξ ( t ) + i ˆ t L ( t (cid:48) )d t (cid:48) (cid:19) , (C.2)28here χ ( z, ρ, t ) is a solution of the 3D field-free quantum harmonic oscillator problem with axial symmetry.We define it to be the ground state of the field-free problem (H-000): χ ( z, ρ, t ) = (cid:104) µω π (cid:105) / e − µω ( z + ρ ) e − i ω t . (C.3)In formula (C.2) the symbol L ( t ) denotes the Lagrangian of the corresponding classical system: L ( t ) = 12 µ ˙ ξ ( t ) − µω ξ ( t ) − f ( t ) ξ ( t ) , (C.4)and ξ ( t ) is the solution of the initial value problem ¨ ξ ( t ) = f ( t ) /µ − ω ξ ( t ) , with { ξ (0) = 0 , ˙ ξ (0) = 0 } . (C.5)We set f ( t ) = F sin ω F t , then the previous equation is that of the forced harmonic oscillator has the solution: ξ ( t ) = F/µω − ω (cid:20) sin ω F t − ω F ω sin ω t (cid:21) . (C.6)These formulae define the analytical solution that we use in Section 6.2 as one of our test cases. There,in Figure 3, we also plot ξ ( t ) , which is actually the expectation value of the coordinate zz