Contents lists available at ScienceDirect

Computer Physics Communications

journal homepage: www.elsevier.com/locate/cpc

Simple, accurate, and efficient implementation of 1-electron atomic time-dependent Schrödinger equation in spherical coordinates*

Serguei Patchkovskii *, H.G. Muller

Max-Born Institute, Max-Born-Straße 2A, 12489 Berlin, Germany

CrossMark

article info

abstract

Modelling atomic processes in intense laser fields often relies on solving the time-dependent Schrödinger equation (TDSE). For processes involving ionisation, such as above-threshold ionisation (ATI) and high-harmonic generation (HHG), this is a formidable task even if only one electron is active. Several powerful ideas for efficient implementation of atomic TDSE were introduced by H.G. Muller some time ago (Muller, 1999), including: separation of Hamiltonian terms into tri-diagonal parts; implicit representation of the spatial derivatives; and use of a rotating reference frame. Here, we extend these techniques to allow for non-uniform radial grids, arbitrary laser field polarisation, and non-Hermitian terms in the Hamiltonian due to the implicit form of the derivatives (previously neglected). We implement the resulting propagator in a parallel Fortran program, adapted for multi-core execution. Cost of TDSE propagation scales linearly with the problem size, enabling full-dimensional calculations of strong-field ATI and HHG spectra for arbitrary field polarisations on a standard desktop PC.

Program summary

Program title: SCID-TDSE: Time-dependent solution of 1-electron atomic Schrödinger equation in strong laser fields.

Catalogue identifier: AEYM_v1_0

Program summary URL: http://cpc.cs.qub.ac.uk/summaries/AEYM_v1_0.html Program obtainable from: CPC Program Library, Queen's University, Belfast, N. Ireland Licensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.html No. of lines in distributed program, including test data, etc.: 334254 No. of bytes in distributed program, including test data, etc.: 20005596 Distribution format: tar.gz

Programming language: Fortran-2003 with OpenMP extensions. Computer: Portable code. Tested on x86_64 Linux. Operating system: Portable code. Tested on x86_64 Linux.

RAM: Memory requirements depend on the input parameters. Time propagation of a given initial state requires O(32 * NR * (LMAX +1) * (MMAX — MMjN + 1)) bytes of memory (double precision), where NR is the number of the radial grid points; LMAX is the maximum desired angular momentum; MMiN and MMAX are respectively the smallest and largest desired angular momentum projections. Preparation of the initial atomic states and analysis of the final wavefunction in terms of the field-free atomic states may require O(64 * N22 * NCPU) bytes of RAM (double precision), where NCPU is the number of processing threads used.

Classification: 2.5, 4.3.

External routines: BLAS and LAPACK (required); libhugetlbfs (optional), DGEFA and DGEDI (LINPACK); routines included with the code.

Article history: Received 3 June 2015 Received in revised form 12 October 2015 Accepted 17 October 2015 Available online 31 October 2015

Keywords:

Time-dependent Schrödinger equation Implicit derivative operators Muller's split propagator Non-Hermitian representation Strong laser fields Linear scaling

* This paper and its associated computer program are available via the Computer Physics Communication homepage on ScienceDirect (http://www.sciencedirect.com/ science/journal/00104655).

* Corresponding author.

E-mail address: Serguei.Patchkovskii@mbi-berlin.de (S. Patchkovskii).

http://dx.doi.org/10.1016/jxpc.2015.10.014

0010-4655/© 2015 The Authors. Published by Elsevier B.V. This is an open access article underthe CC BY license (http://creativecommons.org/licenses/by/4.0/).

Nature of problem: Time propagation of non-relativistic 1-electron Schrödinger equation for a central potential, under the influence of a long-wavelength laser field treated in the velocity-gauge dipole approximation.

Solution method: The propagator is constructed by separating the Hamiltonian into a large number of non-commuting terms, where each term can be handled simply and computationally efficiently (linear scaling). Time-reversibility of the propagator is ensured by combining the individual terms symmetrically around time midpoint (See ref. [1] and the text). The numerical accuracy is achieved through implicit representation of derivatives (accurate to 0(d4) for a uniform grid), combined with variable grid spacing. Restrictions: Ill-conditioned Hamiltonians can occur for some choices or radial grids. The propagator is only approximately norm-conserving; small time steps may be necessary to achieve stable propagation.

Unusual features: Due to the implicit representation of the spatial derivative operators, the overall Hamiltonian is not Hermitian. As the result, the left wavefunction is no longer given by a complex conjugate of the right wavefunction, and must be propagated explicitly.

The code makes no assumptions on the accuracy of numerical types, and can be built for any real or integer kinds supported by the compiler. Detailed instructions for building the code in single, double, and quadruple precision are included.

Running time: Running time is input dependent. Time propagation of the Schrödinger equation scales as 0(NR * (LMAX + 1) * (MMAx — Mmin + 1)) per time step. Preparation of the initial atomic state and analysis of the results in terms of the field-free atomic states may require 0(NR) and 0(NR * (LMAX + 1)) work, respectively. On a 3.6 GHz core i7 desktop with 4 CPU cores available, run times are from 2 s (probability of a perturbative bound- to-bound transition in hydrogen atom) to 1 h (high-harmonic spectrum of a hydrogen atom driven by intense elliptically-polarised IR field).

References:

[1] H.G. Muller, Laser Physics 9 (1999), 138-148.

© 2015 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.Org/licenses/by/4.0/).

1. Introduction

Numerical solution of the Schrodinger equation (SE) is a very old problem, with the first examples dating more than 80 years ago [1]. Solving 1-electron stationary SE in an arbitrary central potential is by now a nearly trivial task. On the other hand, propagating a time-dependent Schrodinger equation (TDSE) in an arbitrary laser field is still a formidable computational challenge, even for a single electron in three-dimensional space. Results of one-electron TDSE simulations are routinely published in leading physics journals [2-6], contribute to the discovery of previously unsuspected phenomena [7,8], and cause intense controversies [9-12]. Despite the apparent simplicity of the problem, development of new, increasingly efficient numerical approaches for solving 1-electron TDSEs in central potentials remains an active area of research [13-26]. The pursuit of numerical efficiency is not entirely surprising: addressing some open research questions in high-harmonic generation (HHG) in optically-dense media and laser filamentation requires simultaneous propagation of 103-106 coupled atomic TDSEs [27-30]. Any improvement in the efficiency of TDSE propagation is therefore welcome.

One of the most elegant and efficient approaches to propagating an atomic TDSE in spherical coordinates has been proposed by H. G. Muller some time ago [15]. The technique relies on three main ideas:

1. Spatial derivative operators, which appear in the field-free and laser-coupling Hamiltonians, are discretised in an implicit form:

f(o) ^ M-1A0f (1)

where f and f(o) are column vectors of the function and its oth derivative (o = 1, 2) on a grid. Tridiagonal matrices M0 and A0 are chosen to minimise the error in the desired derivative. Remarkably, for a uniform grid with spacing h, this form allows derivatives accurate to 0(h4)—despite never explicitly looking beyond the nearest-neighbour grid point.

2. The total time-dependent Hamiltonian is split into a large number of terms—(5LMAX + 2) (2LMAX + 1) in the most general case, where LMAX is the maximum angular momentum considered in the simulation. Each of the individual terms can be implemented using mostly stride-1 memory accesses, and at a linear cost. Furthermore, splitting the Hamiltonian into many terms improves error propagation properties for large grids, and provides ample opportunities for parallel execution.

3. The wavefunction is propagated in a rotating reference frame, chosen such that the vector-potential of the laser field is always along the local Z axis. This choice allows for a substantial reduction in the number of memory accesses and operations needed to evaluate the laser interaction Hamiltonian.

The approach taken by Ref. [15] has been remarkably useful in addressing a number of key problems in atomic strong-field physics [31-37,12]. At the same time, it is subject to a number of limitations. Nearly all of these constraints arise from the ad hoc requirement of Hermiticity imposed on the Hamiltonian matrix. As will be seen below, the implicit representation of the derivative operators (Eq. (1)) on a finite domain necessarily leads to non-Hermitian (non-anti-Hermitian) matrices M-1A2 (M-1 A1).

The non-Hermitian part is small for uniform grids, and for the Laplacian M-1 A2 can be eliminated altogether at the cost of reducing the accuracy to 0(h3) [15]. However, the implicit gradient operator M-1 A1 cannot be made anti-Hermitian without introducing 0(h0) error at a few grid points closest to the origin, even for a uniform grid. Both matrices become essentially non-Hermitian on non-uniform grids.

The goal of this work is to extend the ideas of Ref. [15] to general radial grids, and to present an efficient numerical implementation of the resulting propagator. The rest of this work is organised as follows: Section 2 gives the necessary definitions and derives the propagator expression. Section 3 illustrates the utility of the technique, using the examples of perturbative photo-ionisation, strong-field ionisation, and high-harmonic generation. Finally, Section 4 speculates on possible applications and extensions of the approach. Unless stated otherwise, all expressions in the manuscript use atomic units (h = m = |e| = 1) throughout. We nonetheless retain h, m, and e in most expressions as a dimensionality-analysis aid.

2. Theory

We anticipate dealing with non-Hermitian Hamiltonian matrices below. The algebra of non-self-adjoint matrices is less commonly encountered in chemical physics than its Hermitian counterpart. In order to establish the notation, we would like to restate some key relations here, restricting ourselves to the relevant case of all eigenvectors of the Hamiltonian being genuine [38].

For a non-Hermitian square matrix H of dimension n, the right and left eigenvectors Ri and Li belonging to an eigenvalue ei are column-vectors, defined by:

HRi = e;R; (2)

L] H = eiL[.

For a non-defective matrix H, it is possible to choose n linearly-independent, bi-orthogonal left-/right-eigenvector pairs L,, R;:

L] Rj = Sij. (3)

The condition (3) replaces the familiar normalisation requirement (RHRj = 5j) on eigenfunctions of a Hermitian Hamiltonian. Evaluation of an expectation value (O) of an operator O now requires knowledge of both the left and the right wavefunctions:

<O>i = L] ORi. (4)

We would like to retain the usual time-dependent Schrodinger equation for propagation of the right wavefunctions:

ih dt R = HR. (5)

In the absence of an absorbing boundary condition, we would also like to maintain the total norm of the wavefunction:

dt (L]R) = 0. (6)

It is easy to verify that Eq. (6) is satisfied if the left wavefunction is propagated as follows:

ih dt L = -H] L. (7)

For a Hermitian matrix H, Eq. (7) is simply the complex conjugate of Eq. (5).

If all eigenvalues of the matrix H are real (as should be expected for a physical Hamiltonian), and matrix H is not defective, the right eigenvectors R of Eq. (2) also satisfy a generalised eigenvalue problem with the inner product metric S:

H R = SRE (8)

RtSR = 1 (9)

where R is a column-matrix of the right eigenvectors and E is a diagonal matrix containing the eigenvalues. Hermitian matrices H and S are given by:

H = (r-1)1 er-1 (10)

S = (R-1)1 R-1. (11)

The Hermitian Hamiltonian H may be expensive to evaluate, and therefore unappealing numerically. Nonetheless, its existence is important: it guarantees that the TDSE can be propagated without suffering from exponentially-growing solutions. Furthermore, Eqs. (8) and (9) allow a transparent physical interpretation of the left eigenfunctions L. Comparing Eqs. (3) and (9), the left eigenvectors of Eq. (2) satisfy:

L] =R-1S. (12)

Therefore, the left eigenfunctions L are not independent. Instead, they should be considered as a compact, convenient representation of the inner-product metric S of the Hermitian physical Hamiltonian H.

2.1. Wavefunction ansatz

We are interested in solving a 1-particle, spin-free TDSE in the presence of a time-independent external potential v (?). The time-dependent laser field, treated in the velocity gauge and dipole approximation is specified by the vector-potential A (t). In the Cartesian

laboratory coordinate system, the corresponding Hamiltonian is given by the standard expression: 1 e e2 ( )

Hlab (t) = — P2 "A (t) ■ p + — A2 (t) + v (?) . (13)

2m m 2m

For the reasons comprehensively discussed in Ref. [15], it is computationally advantageous to represent the wavefunction in a rotating reference frame. The z-axis of this frame is chosen along the instantaneous direction of the laser polarisation, so that:

R (t) ■ A (t) =( 0 \ (14)

where R (t) is a 3 x 3 unitary matrix, describing transformation from the laboratory to the rotating frame. The left and right wavefunctions in the laboratory frame are given respectively by:

¿max 1 l

^max 1 L / 1 \

*R (?, t) = E 1 E (r, t) YLm( -R (t) ■ r (15)

L=0 M=-L \ '

lmax 1 L / 1 \

*L (r, t ) = E - E (r, t ) Ylm(-R (t) ■ r (16)

L=0 M=-L \ '

where r = |r |. In Eq. (16), the angular part of the left-wavefunction expansion is conjugated to simplify the book-keeping in evaluation of angular integrals. The phase conventions of spherical harmonics YLM follow Refs. [39,40]. The quantities VjRM (r, t) and (r, t) are the radial wavefunctions in the rotating frame.

Substitution of the right wavefunction (15) into the time-dependent Schrödinger equation with the Hamiltonian (13) yields the rotating-frame Hamiltonian Hrt:

Hrt (t) = T + Hang (t) + Hmix (t) + HA2 (t) + HnAD (t) + V. (17)

In Eq. (17), the kinetic-energy (T), laser-coupling (Hang, Hmix (t), and HA2), non-adiabatic coupling due to frame rotation (HNAD), and the effective-potential (V) terms are defined by:

1 R h 2 9 2 R - *lmYlm = -—Ylm ^ *LM

1 „ ieh L +1 „ ieh L „

rHang (t)-*—YIM =--Art (t) -Cl+1M *LRMyL+1,M + —Art (t )~ Clm *ImYL-1,M (19)

r m r m r

1 „ ieh d ieh d

rH mix (t ) *LMYLM = — Art (t) Cl+1,mYl+1,M^~ *LM + —Art (t ) ClmYl-1,^ — *lm (20)

r m d r m dr

rHA2 (t)-*rmYlm = — A2t (t) *lmYlm (21)

rHNad (t) 1 *lrmYlm = -i*M E Ylm'EMM' (t) (22)

-1 « i L (L + 1)\ R

rv-*RmYlm = {v (r) + (2+2 ) ) *RmYlm. (23)

Coefficients CLm and E—M' are given by:

L2 - M2

Clm = J (24)

EMM' (t) = E i^M'M" (t))* YtDMM» (t) (25)

where DMM, are the Wigner rotation matrices [41,40], corresponding to coordinate transformation R (t):

YLM ( 1r (t) A = E DMM'(t) YLM' (J?) . (26)

r / M'

Eqs. (17)-(25) in principle fully define the Hamiltonian in the rotating frame. Before we proceed to the discussion of the propagator based on this Hamiltonian, it is necessary to define the radial discretisation and implementation of radial derivative operators in Eqs. (18) and (20).

2.2. Implicit representation of spatial derivatives

Evaluation of the kinetic energy operator and the "mixing" term in the laser-interaction operator (Eqs. (18) and (20)) requires discretisation of the first and second derivatives with respect to the radial coordinate r. We follow the approach taken by Ref. [15], and seek the desired discretisation in the form of Eq. (1). However, we do not require matrices Mo and Ao to remain Hermitian. At the interior grid points, Eq. (1) leads to the condition Mof(o) * Aof:

Mu-1f(o) (r-1) + Mi,f(o) (ri) + Mi,i+1f(o) (r+1) * AUi-f (r-1) + Ai,f (n) + Ai,i+1f (ri+1). (27)

After expanding both sides of Eq. (27) in Taylor series around r = ri, we eliminate as many higher-order derivatives of f as possible. Eq. (27) defines Mo and Ao up to an arbitrary overall scale factor. In order to exclude the trivial solution (M = A = 0), we require that the interior rows of M sum up to unity (Mi,i-1 + Mi,i + Mi i+1 = 1,2 < i < NR — 1). If we additionally define f (rNR+1) = f(o) (rNR+^ = 0, the same expression applies to the large-r edge of the radial grid, where it is equivalent to the presence of a hard boundary.

Coulomb singularity at the r = 0 edge of the grid, requires special treatment. For r ^ 0, the implicit-derivative coefficients at the point closest to the origin then must satisfy the equation:

Muf(o) (rO + M1,f(o) (r2) ^ Auf (rO + (r2) (28)

subject to the conditions:

f (0) = 0

f' (0) ={— 2Zf" (0) ifL = 0 (29)

[0 ifL > 0

where f is the effective nuclear charge.

As before, Taylor expansion around r = r1 is followed by elimination of low-order terms in Eq. (28). Similar to Eq. (27), this procedure yields M11, M12,411, and 412 coefficients, which are defined up to an overall multiplicative factor. We choose to fix this factor such that 41,2 matches the expressions for an interior point.

Following tedious, but straightforward algebra, the general-case non-zero elements of matrices M1 and A1 (first derivative) are given

1 di+1 Mi,i—1 = - i+1

2 d2 + didi+1 + d2 =12 (di+di+1)22 i,i 2 d? + didi+1 + d+ 1 d2 + Mi,i+1 = -

^i,i—1 = —

2 d22 + didi+1 + d+

di+1 (2di + di+1)

di (di + di+1) (d2 + didi+1 + d,^) (di+1 — di) (di + di+1)2

didi+1 (d2 + didi+1 + df+1)

d2 (di + 2di+1)

Ai,i+1 = Z-TZ-Z—-^-72—\ • (30)

di+1 (di + di+1) (d2 + didi+1 + di+1)

For L = 0, elements in the upper left corner of M1 and A1 should be replaced by:

(d1 + d2)2 (d1 + 2d2) ((d1 + d2) d1? — 3d1 — d2)

M1,1 =

1.1 (d2 + d2d1 + d2) (2 (d1 + d2) (d1 + 2d2) d1? — 6d2 — 15d2d1 — 8d2) M =_d1 d + 2d2) ((d1 + d2) d1Z — 3d1 — 2d2)_

1.2 (d2 + d2d1 + d2) (2 (d1 + d2) (d1 + 2d2) d1? — 6d2 — 15d2d1 — 8d2)

A __ (d1 + d2)2 (d1 2d2) (2d3z + d2d1 (3 — 2djZ) — 6df + d2) (31)

1,1 = d1d2 (d2 + d2d1 + d2)(2 (d1 + d2)(d1 + 2d2) d1? — 6d2 — 15d2d1 — 8d2) • ( )

In Eqs. (30)-(31), grid spacing di to the left of a point i is given by:

dt = !ri f r1 (32)

i I ri — ri—1 if 1 < i. y '

Note that Eq. (30) require a coordinate of a fictitious grid point rNR+1, which is beyond the rightmost edge of the grid. This value is arbitrary, and can be chosen for computational convenience. ( )

Gradient operator M—1A1 defined by Eqs. (30)-(31) is accurate to O (d3f(5)), both at the grid interior and at the origin. For uniform grids, interior coefficients (Eq. (30)) coincide with the result of Ref. [15].

For the second derivative, the general-case non-zero elements of M2 and A2 are:

m 1 d,2 + didi+1 — d,+1

Mi,i—1 =--

Mi i =-

Mi,i+1 =-

6 di (di + di+1) 1 d2 + 3didi+1 + df+1

6 didi+1

1 — df + didi+1 + d+1

6 di+1 (di + di+1)

A,i_i =

di (di + di+i) 2

Au =--

didi+1

Ai,i+i = d (d \ d ). (33)

di+1 (di + di+1)

For L = 0, the upper left corner of M2 and A2 should be replaced by:

(di + d2) (di (di + d2) (di + 3d2) z - 3 (d2 + 3d2di + d2))

Mi,i = Mi., =

3did2 (di (3di + 4d2) ? - 6 (di + d2)) d3(-z) + d2 (d2Z + 3) + d2di (2d2Z - 3) - 3d2

3d2 (di (3di + 4d2) ? - 6 (di + d2)) 2 (di + d2) (6di - (3di - d2) (di + d2) Z)

A„ = v 2 " v '—^-' ^. (34)

. d2d2 (di (3di + 4d2)z - 6 (di + d2))

For L > 0, the upper left corner becomes:

(di + d2)2 (di + 3d2)

Mi. i =-

Mi 2 = -

3did2 (3di + 4d2) (di - 2d2) (di + d2)

3d2 (3di + 4d2)

2 (3di - d2) (di + d2)2

Ai i =--^-2)( i + 2) . (35)

. d3d2 (3di + 4d2)

The Laplacian operator defined by Eqs. (33)-(35) is accurate to O (d3f(5)) for a general grid. For uniform grids, where coefficients of Eq. (33) reduce to the result of Ref. [i5], accuracy of the derivative for the interior (but not the edge) points improves to O (dAf(6)).

We emphasise that matrices M-i A2 (M-i Ai) defined by Eqs. (30)-(35) are non-Hermitian (non-anti-Hermitian). This property is the inevitable consequence of simultaneously requiring matrices M, A to be tri-diagonal and the resulting operator to have the highest-possible order. For a uniform grid, Hermiticity of the Laplacian M-i A2 can be restored [i5] by requiring M2 to be a linear function of A2, so that the two matrices have common eigenvectors. This additional requirement comes at the cost of decreasing the accuracy of the Laplacian for the edge point to O (d2f(4)). A similar approach is however no longer possible for a non-uniform grid.

Even for a uniform grid, the gradient operator M-iAi cannot be made anti-Hermitian without sacrificing the sparsity, or severely compromising the accuracy. In Ref. [i5], the anti-Hermiticity was restored approximately, by setting M^ i = A^ i = (instead

of the general-case values Mi, i = Mi, i+i = Mi+i,i = i, Ai, i = 0, Ai, i+i = - Ai+i,i = i). Except for a few grid points in the bottom right corner (^50 points in double precision) the resulting matrix M-i Ai is indeed nearly anti-Hermitian. However, large deviations from anti-Hermiticity are present in the bottom-right corner (corresponding to large values of radial coordinate r). More importantly, the anti-Hermitised operator is no longer a good approximation to the radial derivative for small values of r. For the ri grid point, this operator has the accuracy of O (dif(2)) for L > i, decreasing to O (f(i)) for L = 0, thus spoiling the overall accuracy of the first radial derivative. As the result, the anti-Hermitised gradient operator of Ref. [i5] is reliable only for systems where no appreciable wavefunction amplitude may occur at the first few radial grid points, or close to the outer grid edge.

Having established discretisation of all operators in the Hamiltonian (i7)-(23), we can now turn to the propagator.

2.3. Propagator

The Hamiltonian matrix defined by Eqs. (i7)-(35) has the dimension of NH = NR (LMAX + i) for linearly-polarised light (M quantum number is conserved), increasing to NH = NR (LMAX + i)2 for an arbitrary light polarisation. Formally, this matrix has ^3LMAXN^ (linear polarisation) or ^ |LmaxN^ + |LmaxNr (arbitrary polarisation) non-zero elements. Nonetheless, because of the structure imposed on the derivative operators, multiplication of the Hamiltonian by a vector can be performed in O (NH) time. As the result, it is in principle possible to implement the Crank-Nicolson propagator directly for the total Hamiltonian matrix H:

Ucn (t, t + r) = + 2iHT j ^ 1 - liHr) (36)

in formal O (NH) time. This approach is, however, numerically unappealing: the inverse step in Eq. (36) would require an iterative solution of a non-Hermitian linear problem of dimension NH. Convergence of standard iterative solvers for such problems is by no means guaranteed. Even for diagonally-dominant problems (r ^ max (|H|)), a large number of iterations may be required to solve for the inverse.

Instead, we follow Ref. [i5], and split the propagator into a large number of non-commuting individual terms of low dimensionality. For the right and left wavefunctions, respectively:

Uhrgm (t, t') =UR;reerv (t') URtomURJt (t, t') Uixr (t) (37)

UHgM (t, t') =U,asreerv (t') UaLtomUrLot (t, t') UtomUf (t) (38)

where t' - t = 4t. The individual terms of the right propagator are given by:

mmax lmax

Uatom = 0 II Uatom,LM (39)

M=Mmin L=|M|

mmax /lmax ^ \ / lmax-1 \ / lmax-1 \

Ulaser = 0 ( n UR2,LM O Umix,LM Uang,LM II O Umix,LMUang,LM I (40)

M=Mmin V=|M| ) V=|M| + 1,|M|+3,... ) V=|M|,|M|+2,... )

mmax / lmax-1 \ / lmax-1 \ / lmax ^ \

Ulaser = 0 ( n Uang,LMUmix,LM 11 ]""[ Uang,LM Umix,LM 11 ]""[ Ua2,LM I (41)

M=Mmin \l=|m|,|m|+2,... / V=|M|+1,|M|+3,... / V=|M| /

lmax nR na

U rRot =nn nU Rot,L,k,a. (42)

L=1 k=1 a=1

The left-propagator terms are in the same form, except for UaLtom Lm , UaLng Lm , U mix Lm , and UaL2 Lm appearing in the place of the corresponding right sub-block propagators. In Eq. (39) all individual L, M terms commute, and may be computed in parallel. In Eqs. (40) and (41), all terms within the brackets commute, and may be evaluated in any order. Individual M terms also commute, creating additional opportunities for parallel execution. Finally, in Eq. (42), terms with different L and k indices may be evaluated in any order. The rotation-substep terms (index a) must be applied in the order given.

The sub-block propagators appearing in Eqs. (39)-(42) are each in the Crank-Nicolson form:

Uatom,LM = (1 + iTHat,LM) (1 - iTHat,LM) (43)

Uang,LM = (1 + iTHang,LM) (1 - iTHang,LM) (44)

Umix,LM = (1 + iTHmix,LM) i1 - iTHmix,LM) (45)

Ua2,LM = (1 + iTHa2,LM) (1 - iTHa2,LM) (46)

Urot,L,k,a = (1 + iTHrot,L,k,a) i1 - iTHrot,L,k,a) . (47)

The left sub-block propagators are in the same form, except for each H being replaced by —HT. Matrices Hat,LM and Ha2,LM are of dimension Nr. Matrices Hang LM and Hmjx,LM are of dimension 2NR, and couple (L, M) and (L + 1, M) radial sub-blocks. Finally, matrices Hrot,L,k,a are of dimension 2L + 1, coupling all M values with the same L and radial index k. Evaluation of the ÛR/LLM terms, which involve a multiplicative scalar potential ^A'2t, is trivial. The remaining sub-block propagators (43)-(47) are discussed below.

2.3.1. Field-free terms (Eq. (43))

The field-free atomic Hamiltonian matrix is given by:

Hat,LM = ——M—L A2,l + Vl (48)

, 2m 2,L ,

where tri-diagonal matrices M2,L and A2,L are defined by Eqs. (33)-(35), while the diagonal matrix VL contains values of the multiplicative operator (23) at radial grid points. This matrix of dimension NR appears in four individual terms (Eq. (43) and the corresponding left propagator):

Xl,m = (1 — iTHat,LM) Yl,m (49)

XL,M = (1 + iT Hat,LM ) YL,M (50)

XL,M = i1 + iTHIt,LM) YL,M (51)

XL,M = (1 — iT HIt,LM ) YL,M (52)

where YL M is a vector of dimension NR, containing the L, M radial wavefunction. Eqs. (49) and (51) can be evaluated efficiently as written, keeping in mind that both multiplication by a tri-diagonal matrix and solution of a tri-diagonal linear system are O (NR) operations. Eq. (50) is equivalent to:

M2,L — iT^2,l + iTM2,lVl^ Xl,m = M2,lYl,m. (50a)

This is again a tri-diagonal linear system with respect to XL M. The linear system matrices in Eq. (50a) can be precomputed and stored. Finally, Eq. (52) is equivalent to:

m2,l — y A{l + iTVTM,l) Zl,m = Yl,m (52a)

Xl,m = M2,l Zl,m . (52b)

Again, the tri-diagonal linear system matrices can be precomputed to improve efficiency.

2.3.2. Laser coupling: "Angular" term (Eq. (44))

The explicit matrix form of the right "angular" sub-block propagator (Eq. (44)) in the laser interaction Hamiltonian is:

YL+1,M

XL,M 1 -br-1 " -1 " 1 br-1

XL+1,M br-1 1 -br-1 1

where 1 is a unit matrix, while r is a diagonal matrix containing radial coordinates of the grid points. Both matrices are of dimension NR. Coefficient b is given by:

b = t— Art (t) (L + 1) Cl+i,m m

where CL+1,M is defined by Eq. (24). Eq. (53) is equivalent to NR independent systems of 2 x 2 linear equations, with the solution:

(53a) (53b)

Because the Hang matrix is anti-Hermitian, the left "angular" sub-block propagator coincides with the right propagator. Because coefficient b is time-dependent, factors in Eqs. (53a) and (53b) must be re-evaluated on each time step.

Xl,m = (r2 + b21) 1 ((r2 - b21) Yl,m + 2brYL+u Xl+1,m = (r2 + b21)-1 ((r2 - b21) Yl+1,m - 2brYLM

2.3.3. Laser coupling: "Mixing" term (Eq. (45))

The explicit matrix form of the right (Eq. (45)) and left "mixing" sub-block propagators is:

XL,M 1 -c M-L+1A1,L+1 " -1 1 c M -L+1 A1,L+1

XL+1,M _ -cM-,LAU 1 _ cM-,L A1,l 1

XL,M 1 _cAT M-T " LA1,LM1,L 1 ■ 1 cAT M-T c A1,LM1,L

XL+1,M . -cA1,L+1M-L+1 1 . cAT,L+1M-L+1 1

YL+1,M

YL+1,M

Zl M -D N + E

ZL+1 N -E M + D

where matrices M1 and A1 are given by Eqs. (30)-(31), and time-dependent coefficient c is: eh

c = t— Art (t) Cl+1,m . m

Eq. (55) is equivalent to:

m = ^ (mu + mu+1)

n = 2 (Mu - Mu+1)

D = 2 (AU + AU+1)

WL = YL,M + YL+1,M WL+1 = YL,M - YL+1,M

tL = (M + D) wL + (N - E) wL+1 tL+1 = (M - D) wL+1 + (N + E) wL

. - Jl+1

XL,M = ^ (ZL + ZL+1) 1

XL+1,M = 2 (ZL - ZL+1) .

The key step in evaluating Eqs. (55a) is the solution of the linear system for ZL and ZL+1. Two distinct cases have to be considered. For L > 1, M1,L+1 = M1L and A1,l+1 = A1,l, so that N = E = 0, and the solution reduces to the case considered in Ref. [15]. For L = 0, N and E no longer vanish. However, these matrices contain only two non-zero matrix elements (the first two columns of the first row). The solutions of the linear system are then given by the Sherman-Morrison formula [42], using:

M - D - N + E 0

0 M + D - N - E

as the unperturbed matrix, and:

N -E N + E

N -E N + E

as the correction term.

Evaluation of the Sherman-Morrison expression for L = 0 requires solving two auxiliary tri-diagonal systems of linear equations (one for M — D — N + E; another for M + D — N — E), each for two right-hand sides. For L > 1, where N and E vanish, two tridiagonal linear systems (M — D and M + D) have to be solved, each for a single right-hand side. The remaining operations involved in evaluating Eq. (55)a are the same in both cases. The cost of solving a tri-diagonal linear system is dominated by the factorisation step; solving for an additional right-hand side using the same factorisation adds only a minor cost. The overall cost of the L = 0 and L > 1 terms is therefore essentially identical.

Similar transformation of Eq. (56) yields:

M = 2 (M{,L + M1,

N = 2 (MÎ,L — MÎ,L+1) D = 2 «L + 4l+1)

E = 2 «L — 4

Wl = M—,t+1Yl,m + M—,T Yl+1,m

WL 1 M

TL = Tl+1 =

Y, M - M—! Y,

1, l+1yl , M

1LYL+1, M

(M + D) WL + (N + E) W,+1 (M — D) Wl+1 + (N — E) WL

Zl M — D N — E

ZL+1 N+E M + D

XL, M = ^ M1, L+1 (ZL + ZL+1) 1 T

XL+1, M = ^ M1L (ZL — ZL+1) .

As before, the difference matrices N and E are either zero (L > 1), or contain at most two non-zero elements (L = 0: the first two rows of the first column). Again, these equations can be implemented using only multiplications by tri-diagonal matrices and solutions of tri-diagonal linear systems.

Similar to the case of the "angular" propagator (Eqs. (53a), (53b)), linear system matrices in Eqs. (55)-(56) have to be recomputed on each time step.

2.3.4. Laser coupling: Non-adiabatic "frame-rotation" term (Eq. (47))

The effect of the non-adiabatic term in the propagator (Eq. (47)) is to change the direction of the local z axis from A (t) /Art (t) to A (t') /Art (t').The corresponding propagator is known analytically:

n <L,,,a (t , t') = (D(L) (a' , d', y'))* (D(l) (a , d, y))T (60)

where D(L) (a, d, y) are the Wigner rotation matrices [41] for the angular momentum L and Euler angles a, d, y .The Euler angles at times t and t' are chosen such that Eq. (14) is satisfied, fixing angles a and d:

(a,d,y) = (0,Q, 0),

where Q and 0 are the polar and azimuthal coordinates of the vector-potential in the laboratory frame. The choice of the y angle is, in principle, arbitrary. The explicit form of the corresponding left propagator is clear from comparison of Eqs. (15) and (16):

UkLXa^L = (Urot,L,M (^T)* (61)

so that both propagators can be implemented by essentially the same routine. The explicit, dense propagator expressions (60) and (61) work well for small L (L < 5). It is, however, unsatisfactory for large L, since the general expression for the Wigner matrix [41] rapidly loses significance with increasing L values. Furthermore, the change in Euler angles in a single time step is typically small, so that only M' = M, M ± 1 terms are coupled appreciably, and the propagator becomes sparse. It is therefore advantageous to treat this term as a sequence of small-angle rotations [15].

Transforming the local coordinate system defined by Euler angles (0, Q, 0) into a (nearby) frame determined by Euler angles (0 + S0, Q + SQ, 0) requires a sequence of small rotations around each of the local coordinate axes xt, yi, and z;:

Sax = — sin (0) 8<p (62)

Say = —50

Saz = cos (0) 8<fi.

(For small rotations, the order of operations is immaterial.) From the definition of spherical harmonics it then follows that, for a small rotation, Eq. (22) can be rewritten as1 :

rHnad (t) 1 &lRmYlmSt = (-hMSaz) &Mylm + (ih (Say + iSax) D+m) VlmYl,m+i + (ih (-Soy + iSax) D-) ^lmY,m-1 (63)

D+m = (L - M)(L + M + 1) D-m = 1 V(L + M)(L - M + 1).

Because HNAD forms a tri-diagonal matrix Hrot for each value of L and radial grid point k, evaluation of the small-angle rotational propagator of Eq. (47) is straightforward.

In some situations, such as elliptically-polarised fields near the zero of the major field component, the local frame may undergo rapid rotation within a single time step of 4t. Rather than decreasing the time step or increasing the order of the propagator (as is done in Ref. [15]), it is preferable to sub-divide the time step into NA sub-steps (Eq. (42)). For a given L, the Crank-Nicolson-type infinitesimal propagator (42) then requires ^26 (2L + 1) NANR flops. On the other hand, application of the dense finite-rotation propagator (60) requires ^5 (2L + 1)3 Nr flops, plus ^26 (2L + 1)2 NA flops for computing the finite-rotation matrix from the sequence of small-angle propagators. For sufficiently large NA and NR, it becomes advantageous to use the dense-matrix rotation code for such subdivided rotations.

3. Tests and numerical examples

The implementation, input parameters, and possible issues with numerical accuracy of the present implementation are extensively documented in the accompanying source code and program documentation. It therefore remains only to give a few examples of numerical applications.

3.1. Convergence of state energies in hydrogen atom

Convergence of the bound states of the field-free Hamiltonian is illustrated in Fig. 1. Calculations were performed on a uniform grid, extending to a hard boundary at Rmax = 82 Bohr. Grid spacing varied from 1.28 Bohr to ^3 x 10-4 Bohr. Selected eigenstates of the Hamiltonian were calculated using inverse iterations [42]:

Ri+1 = (Hat - c¡ 1)-1 Ri (66)

Li+i = (Ht - e;1)-1 Li (67)

Ri+1 = R¡+1 / |R¡+1 J (68)

Li+1 = Li+1 /LJ+1 Ri+1 (69)

ei+1 = c¡ + 1/L[ R+1 (70)

where Li, Ri, and c; are successive approximations to the left- and right-eigenvalues and the eigenvectors, respectively. Random starting guess was used for the left and right eigenvectors. In the absence of the absorbing boundary, all eigenvalues of the discretised field-free atomic Hamiltonian (Eq. (48)) are real, even though the Hamiltonian matrix is not Hermitian.

Due to the high order of the kinetic energy operator, convergence with increasing grid density is rapid. In single precision, bound-state energies converge to 10-5 Hartree already for the grid spacing d = 0.2 Bohr. Further increase in grid density in fact decreases the accuracy of the eigenvalues in a single-precision calculation.

In double precision, 1s, 2p, and 3d eigenvalues converge to 10-11-10-12 Hartree for d ^ 10-2 Bohr. Further increase in grid density again decreases the accuracy, while the quad-precision calculations continue to converge to the analytical result. Similar behaviour is seen for the 2s, 3s, 3p, and 3d state energies (not shown). At the same time, the n = 4 states (4s, 4p, 4d) converge to ^ - 32 + cl Hartree (c0 = 5.78 x 10-10; c1 = 3.87 x 10-10; e2 = 1.61 x 10-10), instead of the expected analytical free-atom result (- 32), both in double and in quadruple precision. This is not a numerical artifact: The radial extent of the n = 4 hydrogenic solutions is comparable to the position of the outer grid boundary. Radial confinement of the wavefunction then causes the increase of the total energy.

Linear regression of the linear segment of the 1s (q), 2p (q), and 3d (q) in the range of 0.00125 < d < 0.64 Bohr gives IE — Eexactl ^ d4 00±0 05. This scaling behaviour is consistent with the accuracy of the Laplacian operator in the field-free Hamiltonian (see Section 2.2).

3.2. Perturbative ionisation of hydrogen atom with linearly-polarised light

Convergence of the calculated 1-photon photoionisation cross-sections of 1s, 2s, 2p0, and 2p1 states of the hydrogen atom is illustrated in Fig. 2 and Table 1. Calculations are performed on a uniform grid, extending to Rmax = 180 Bohr. Smaller radial grids (Rmax = 90 Bohr) were found to yield unconverged results for the cross-sections, with at most 4 significant digits regardless of the grid density. A transmission-free complex absorbing potential [43], with parameters kmin = 0.2, S = 0.2, was applied starting at RCAP = 146 Bohr. Grid spacing varied from 1.28 Bohr to 0.01 Bohr.

1 Note that in Eq. (63) we use spherical harmonics phase conventions consistent with the definition used in Mathematica software package [39,40]. Definition used in [41]

differs by an additional phase factor of iL+2M.

1e-02 1e-04 1e-06 1e-08 1e-10 1e-12 1e-14 1e-16 1e-02 1e-04 1e-06 1e-08 1e-10 1e-12 1e-14 1e-16 1e-02 1e-04 1e-06 1e-08 1e-10 1e-12 1e-14 1e-16

-1-1—I—I I I I 11-1-1—I—I I I I I |-1-1—I—I I I I I

+ 1s(s) 4s (d) - -X- 4s (s) -e-1s(q) -a-1s(d) -A- 4s (q)

\—i i i 11111-1—i i i 11 II|-1—i i i 11 II|

2p(s) -■ 4p (d) |- 4p (s) -©- 2p (q) 2p(d) -A- 4p (q)

3d (s) -■ 4d (d) - 4d (s) -e- 3d (q) 3d (d) -A- 4d (q)

0.010 0.100 d, Bohr

Fig. 1. Convergence of the bound states of the field-free Hamiltonian in a hydrogen atom with the radial grid density. Calculations are performed in single (s), double (d), and quadruple (q) precision. Uniform grid with spacing d (Bohr) extends to 82 Bohr; Absolute error in the total energy is in Hartree. Panels: (a) L = 0 states; (b) L = 1 states; (c) L = 2 states. Please note that the 4s (d)/4s (q), 4p (d)/4p (q), and 4d (d)/4d (q) curves are not distinguishable on the logarithmic scale for d > 5 x 10-3 Bohr.

Table 1

Photoionisation cross-sections in Bohr2 for low-lying states of hydrogen. The first digit deviating from the corresponding extrapolated numerical result is marked in bold. Truncated Gaussian pulse with central photon energy 1 Hartree (»27.2 eV) and FWHM 60 au[t] (»1.45 fs). See text for details.

Central3 "Exact"6 Extrap.c

(7 (1s, L = -- 1), x102 3.326 053 3.329 278 3.329281

a (2s, L = = 1), x103 3.437 792 3.440751 3.440754

a(2p0, L = 0), x104 1.041755 1.043371 1.043372

a(2p0, L = 2), x104 4.583 723 4.592 508 4.592 513

a(2pt, L = 2), x 104 3.437 792 3.444381 3.444385

a Analytical 1-photon cross-section, at the central frequency of the pulse. b Convolution of the analytical 1-photon cross-section with the power spectrum of the pulse. The integral is converged to »7 significant digits; see text.

c Extrapolation of the quadruple-precision numerical result for d = 0.04, 0.02, and 0.01 Bohr grid to d = 0. To the number of significant digits given, extrapolated value coincides with the d = 0.04 result; see Fig. 2.

Vector-potential of the Z-polarised laser pulse was chosen in the truncated-Gaussian form:

Az (t) = A0 cos (at - t0) f (t - t0)

f (t) =

exp (-at2)

2 /n |t| - t1

exp ( -a t1 + - (t2 - t1) tan ----

n V 2 t2 - t1

if It| < t1 if t1 < It I < t2 if t2 < It I

d, Bohr

Fig. 2. Convergence of 1-photon photoionisation cross-sections for a hydrogen atom with the radial grid density. Laser field is linearly-polarised along the z direction. The photoelectron angular momentum for each initial state is given in parentheses. Calculations are performed in double (d) and quadruple (q) precision.

with the parameters: a « 3.854 x 10-4 au[t]-2, t0 = 150 au[t], t1 = 90 au[t], t2 = 145 au[t], corresponding to a 60 au[t] («1.45 fs) pulse (FWHM power), centred «3.63 fs after simulation start. The simulation continued for 905 au[t] («21.9 fs) after the pulse end, to allow the continuum part of the wavefunction to reach the absorbing boundary. Photoionisation cross-sections were calculated from the probability P of the atom remaining in the initial state at the end of the simulation:

a = - ln (P)//ph (73)

where iph is the photon fluence of the pulse. Double-precision simulations used A0 = 1.069 x 10-3 au, corresponding to the peak intensity of «4 x 1010 W cm-2. Peak intensity was reduced to «4 x 108 W cm-2 (A0 = 1.069 x 10-4 au) for quadruple-precision runs to reduce multi-photon contributions.The time step was5x 10-3 au[t] (2.5 x 10-3 au[t]) in double (quadruple)-precision simulations, corresponding to «0.121 as («0.060 as).

Analytical expressions for the cross-sections for 1-photon ionisation of hydrogenic states are well known [44]. However, some care must be taken in comparing these cross-sections to the simulation results: The spectral width of the short 1.45 fs Gaussian pulse is «1.26 eV. The effective ionisation cross-sections for this pulse are given by a convolution of the analytical monochromatic cross-section with the power spectrum of the pulse. These effective cross-section are substantially different from the result at the central frequency (see Table 1). Furthermore, pulse envelope truncation adds an slowly-decaying pedestal to the power spectrum of the pulse, at the level of «10-7 of the peak power density. As the result, it has proven difficult to converge the analytical expression for the effective 1-photon cross-section to better than 7 significant digits. Instead, we have chosen to extrapolate the results of the quadruple-precision simulations with the grid spacing d = 0.04, 0.02, and 0.01 Bohr, assuming a power-law dependence of the residual error on the grid density. The extrapolation agrees with the analytical result within the accuracy expected from the latter (see Table 1).

The double- and quadruple-precision results are nearly identical for grid spacings d > 0.1 Bohr. Double-precision calculations do not benefit from higher-density grids, with the residual errors in the cross-sections stabilising at «10-8 Bohr2 for d < 0.1. This apparent convergence plateau is a consequence of the numerical accuracy, photon fluence, and the peak pulse intensity. A 1.45 fs, 4 x 1010 W cm-2 pulse with central photon energy of «27.2 eV has photon fluence of «3.98 x 10-4 Bohr-2. If a surviving population of the initial state in Eq. (73) is represented with a 14-digit double-precision real value, we can expect Eq. (73) to be accurate to at most 10-10 Bohr2. In practice, the accuracy of the calculated 1-photon cross-section is limited by the presence of multi-photon effects. For example, for the 1s initial state, 2-photon ionisation at ipeak = 4 x 1010 W cm-1 accounts for «1.9 x 10-12 ionisation probability in the L = 0, 2 continuum channels (see below). The corresponding contribution in the apparent 1-photon cross-section is «4.9 x 10-9 Bohr2, accounting for nearly the entire discrepancy with the converged result. The peak intensity of «4 x 108 W cm-2 used in the quadruple-precision calculation is expected to reduce the multi-photon error to below 10-10 Bohr2.

Linear regression of the quadruple-precision residual errors in the range of 0.01 < d < 0.64 Bohr (Fig. 2) gives |agrjd - aextrap| « d4 08±0 04. This scaling is somewhat better than could be expected from the properties of the radial-gradient operator alone (see Section 2.2). The improved scaling likely reflects the combined effects of the simultaneous convergence of the laser-interaction Hamiltonian and the field-free eigenspectrum.

It is instructive to calculate the energy-resolved photoelectron spectrum in addition to the total photoionisation cross-sections. The most rigorous approach for extracting the photoelectron spectra from an absorbing-boundary TDSE simulation is provided by the t-SURF method [45]. An implementation of t-SURF is planned, but is not presently available in our code. Instead, for a short pulse used here, it is also possible to construct the photoelectron spectrum as described below.

For a given angular channel L, M, the total electron population within the channel is given by:

PLM =L[M RLM (74)

where Llm and RLM are respectively the left and right discretised wavefunctions within the channel. As long as the field-free Hamiltonian matrix for the L, M channel is not defective, its left and right eigenvectors (ULMi and VLMi) form a complete biorthogonal set, so that:

J^VLMUM =1. (75)

Inserting the unity of Eq. (75) into the right-hand-side of Eq. (74), we obtain:

PLM = aLMi (76)

aLMi = (L[MVLMi) (ULMiRLM) . (77)

It is therefore natural to treat aLMi as the population of the field-free eigenstate L, M, i with the complex energy ELMi - 2rLMi, and to construct the energy- and channel-resolved photoelectron spectrum as:

A / 4 ln 2 / 4ln2 A

Pl,m (E) =2_, aLMi exp I —— (E - ELMi) I . (78)

i V LMi V ' LMi '

It should be emphasised what Eq. (78) is ad hoc, and lacks a rigorous justification. The ''state populations" aLMi are in general complex. As long as the time-dependent wavefunctions have not reached the absorbing boundary, overall channel populations PLM are approximately real, and the sum of the imaginary parts of aLMi nearly vanishes. The non-physical, imaginary part of pL^M (E) is also typically small. We therefore choose to treat the real part of the pL^M (E) as an approximation to the photoelectron spectrum.

Despite its obvious shortcomings, Eq. (78) leads to remarkably good photoelectron spectra, as can be seen from Fig. 3.The spectra were calculated for a truncated Gaussian pulse with the central energy »27.21 eV, peak intensity of »4 x 1010 W cm-2, and FWHM duration of »1.45 fs (this is the same pulse used in double-precision calculations of the 1-photon cross-sections above). The calculation used a variable-density radial grid, with 0.02 Bohr grid spacing for r < 2 Bohr, gradually increasing to 0.1 Bohr for r > 3.74 Bohr. A transmissionfree complex absorbing potential [43] (kmjn = 0.2, S = 0.2) was applied starting at r = 785.5 Bohr. The time step was 0.005 au[t]. Angular momentum channels with L < 6, M = 0 were included in the simulation. The calculation was performed using double-precision arithmetics. The photoelectron spectrum was evaluated at T = 300 au[t] (»7.26 fs), 5 au[t] after the end of the laser pulse. The 1-photon ionisation cross-sections calculated from the resulting photoelectron spectrum agree with the converged result (above) to 10-10 cm2.

The total energy-resolved spectra are shown in the top panel of Fig. 3, for the 1s (red) and 2p0 (green) initial states. In both cases, the clear ATI progression is visible, including 1-, 2- and 3-photon peaks (1s: 13.6,40.8 and 68.0 eV; 2p0: 23.8, 51.0, and 78.2 eV). These peaks are nearly Gaussian in shape. The ATI progression is superimposed on the broad background matching the power spectrum of the pedestal of the short ionising pulse, times the 1-photon ionisation cross-section.

Further details emerge in the channel-resolved spectra for the 1s initial state (middle panel). As expected from the angular momentum conservation rules, the L = 1 continuum channel contains the 1- and 3-photon peaks. The two 3-photon peaks (+2 - 1 at 13.6 eV and +3 at 68.0 eV) are seen in the L = 3 channel. The 4-photon peaks, which were obscured by the 1-photon pedestal in the total ATI spectrum, are clearly visible in the L = 4 channel (+3 - 1 at 40.8 eV and +4 at 95.2 eV). In the L = 0 and L = 2 channels, the most prominent feature is the 2-photon peak at 40.8 eV. The pedestal of the power spectrum of the pulse is also clearly visible in the L = 0, 2 channels; as expected for a 2-photon transition, it appears 10 orders of magnitude below the peak maximum, compared to 6 orders of magnitude for the 1-photon pedestal in the L = 1 channel. The L = 0, 2 channels also exhibit some initially unexpected features, namely an asymmetrical peak at 23.8 eV, and a Gaussian peak at 13.6 eV. The 23.8 eV peak is due to an intermediate resonance with the 2p0 — 1s transition, accessed by the wings of the power spectrum of the pulse. The 13.6 eV peaks in L = 0, 2 channels arise due to the final pulse duration; they would not be present in a CW field with the same spectral content [46]. These peaks are the energy-domain manifestation of the spatial asymmetry of the photoelectron spectrum [46]. Finally, numerical noise is evident in the L = 0 channel at the 10-27 eV-1 level, 22 orders of magnitude below the main 1-photon peak.

The 2p0 channel-resolved spectrum (bottom panel of Fig. 3) contains similar features, adjusted for the difference in the initial-state energy and angular momentum. Additionally, a prominent feature at zero photoelectron energy in the L = 1, 3 channels is due to a Raman transition through the L = 0, 2 1-photon continuum.

3.3. High-harmonic generation from 1s state of hydrogen atom with linearly and elliptically polarised light

Having established the accuracy and convergence properties of the new TDSE propagator, we would like to give a few examples of strong-field calculations, which are the intended application domain for the technique.

The first example is high-harmonic generation by a hydrogen atom, initially in the 1s state, driven by an intense, short infrared pulse (800 nm, 1014 W cm-2,4.84 fs FWHM). This example uses a variable-density grid, with 0.04 Bohr grid spacing for r < 0.4 Bohr, gradually increasing to 0.4 Bohr for r > 1.67 Bohr. The grid extends to »110.0 Bohr, with the absorbing boundary (kmjn = 0.2, S = 0.2) starting at »78 Bohr. Angular-momentum channels with L < 20 were included in the simulation. The simulation continued to »24.2 fs, »12.2 fs after the end of the laser pulse. Time step was 0.0025 au[t] (»60.5 zs). The high-harmonic spectrum was calculated from the Fourier transform d (Q) of the real part of the expectation value of the dipole operator d (t):

„ 12

Ihhg (Q) a Q4 |d (Q)l (79)

<u 1e-14

a. 1e-18

<u 1e-14

<0 1e-14

Q. 1e-18

1s H-2p0 -B-

1s(L=0) H-1s(L=1) -X-1 s(L=2) -B-1s(L=3) *-1 s(L=4) -e-

2p0(L=0) H-2p0(L=1) ^ 2p0(L=2) -B-2p0(L=3) 2p0(L=4) -9-

Fig. 3. Simulated photoelectron spectra for the 1s and 2p0 states of hydrogen atom. The laser pulse is a truncated Gaussian with the peak intensity of 4 x 1010 W cm2, central photon energy of 27.21 eV and full width at half maximum of 1.45 fs, polarised along the Cartesian z direction. The x axis is the photoelectron energy in eV. The y axis is ionisation probability density in eV-1. Panels: (a) Total energy-resolved spectra; (b) Energy- and channel-resolved spectra for the 1s initial state. (c) Energy- and channel-resolved spectra for the 2p0 initial state.

9 13 17 21 25 29 Harmonic order

Fig. 4. High-harmonic spectrum of the 1s state of hydrogen atom, driven by an intense short IR pulse (800 nm, 1014W cm-2, 4.84 fs FWHM). Blue (dark) line: linearly-polarised driving field. Green (light) line: elliptically-polarised driving field; Ay/Az = 0.1. Gray vertical lines are at hQ = ¡p («H9) and hQ = ¡p + 3.17Up («H21). The inset shows the vector-potential of the dominant (Az, black (dark)) and the minor (Ay, read (light)) field component. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

where:

d (t) = j VL (r, t) rVR (?, t) dr.

Simulation of the spectra in Fig. 4 required 5 min and 1 h on desktop PC, respectively for the linearly- and elliptically-polarised driving fields.

^ 0.01 o

ri 0.001 LU

■i1 0.0001 EC

1e-05 1e-06

0 10 20 30 40 50 60 70 Photoelectron energy, eV

Fig. 5. Above-threshold ionisation photoelectron spectrum of the 2p0 state of hydrogen atom, driven by an intense short IR pulse (800 nm, 1014 W cm 2, 4.84 fs FWHM), polarised along z direction. The laser field vector-potential is shown in the inset of Fig. 4. Gray vertical lines are at 2 and 10 Up (»12 and 60 eV). Top horizontal axis gives the number of photons at the nominal carrier frequency (1.55 eV) above Ip + Up (3.40 + 6.00 eV) The inset magnifies the low-energy structure of the spectrum in the 0-5 eV energy range. The dense line structure at the low-energy end of the spectrum (<0.5 eV) is an artifact of the finite energy grid resolution.

Calculated high-harmonic spectra are shown in Fig. 4, with the spectrum for the linearly-polarised driver shown in red. Elliptically-polarised pulse produces harmonic spectrum shown in green. The spectra consist of highly-structured peaks close to the expected odd harmonics of the driving field, extending to the harmonic 25 (»38.7 eV). The complex shape of the harmonics is due to a combination of two factors: short-long interference and inter-cycle interference. The free-oscillation radius of an electron driven by 800 nm, 1014 W cm-2 field is only 17 Bohr. Because the absorbing boundary starts at 78 Bohr, all of short, long, and multiple-return trajectories are supported by the radial grid and contribute to the harmonic emission. Emission due to these trajectories interferes in the harmonic spectra. Furthermore, the 4.84 fs pulse duration supports only 4 significant ionisation events at the maxima of the electric field. Because of the rapid variation of the pulse envelope, each of the ionisation events is associated with a distinct set of trajectories, which create additional opportunities for interference. Due to the high ionisation potential of the 1s hydrogen, 98.9% of the atoms remain in the initial state after the end of the pulse. As the result, the left-/right-symmetry breaking is small, and even harmonics, even though clearly visible (H10, H12, H14, H16), remain suppressed.

It is instructive to compare the HHG spectra calculated for the linear- and elliptical driving fields at the same field intensity. Two regions are clearly evident in the spectrum. For the sub-threshold harmonics below H9 (13.6 eV), the two spectra are virtually indistinguishable. In this energy range, Coulomb focusing overwhelms the steering effect of the minor component of the laser field, and total harmonic emission remains unaffected. Above the Ip threshold, the liberated electron wavepacket can be steered away from the core, leading to the universal suppression of the harmonic emission [47].

3.4. Calculation of the ATI spectrum of hydrogen 2p0 state in strong infrared laser field

Our second example is strong-field ionisation of a hydrogen atom, initially in the 2p0 state (Ip = 3.40 eV) by an intense, linearly polarised laser pulse. We use the same laser pulse as in the HHG example above (800 nm, 1014 W cm-2,4.84 fs FWHM, ellipticity e = 0), as well as the same grid density and time-propagation parameters. In order to accommodate the broader range of angular momenta and maintain the entire wavefunction within the simulation volume, we include angular momenta L < 30, and extend the radial grid to »870 Bohr. The photoelectron spectrum is calculated using Eq. (78). The entire simulation requires »40 min on a desktop PC.

Due to the low ionisation potential, the initial state population is nearly completely depleted by the 1014 W cm-2, 4.84 fs pulse, with only 0.03% of the atoms remaining in the initial 2p0 state. At the same time, approximately 24% of the atoms remain bound [48]. The maximum of the bound-state population is found around n = 5, l = 4 quantum numbers; however many states are populated [48]. The most-populated individual final states are 5f0 (2.8%), 5d0 (1.8%), 6f0 (1.5%), and 5g0,7f0 (each with 1.2%).

Simulated photoelectron spectrum is shown in Fig. 5. As expected [49], three distinct regions are clearly discernible in the spectrum. A regular ATI progression starts at »5 eV, corresponding to absorption of 10 IR photons by the initial 2p0 state. The slope of the progression changes around 12 eV(2Up » 12 eV), indicating transition from direct to rescattered electrons [49,50]. The rescattering plateau continues to »55 eV, close to the expected 10Up cutoff at 60 eV. The low-energy part of the spectrum (Fig. 5 inset) shows a complex line structure, with evidence of a second ATI progression offset by »0.8 eV and asymmetric peaks. These structures arise due to resonantly-assisted ionisation. Finally, the low-energy structure, attributed to soft recollisions [51] is seen below 0.3 eV.

4. Conclusions

We developed a new approach for numerical solution of time-dependent 1-electron atomic Schrodinger equation in intense laser fields. The technique is based on earlier ideas of Ref. [15], including use of implicit representation of derivative operators, use of the rotating reference frame, and the overall structure of the propagator. These ideas are extended to support non-uniform radial grids, arbitrary laser field configurations, and higher-accuracy derivative operators.

Non-uniform radial grids permit an improved description of the wavefunction in parts of space where the potential is varying rapidly (typically close to the nucleus), without incurring the cost of the high-density grid away from the origin where it is not needed. However, this possibility comes at a cost, and must be wielded with care: non-uniform grids reduce the accuracy of the implicit derivatives. The deterioration is especially pronounced if the grid density varies rapidly. Furthermore, stability of the propagator (37)-(38) is determined

Photons

10 15 20 25 30 35 40 45 50

by the densest part of the radial grid. As the result, the numerical advantage of a non-uniform grid in a practical simulation is often not as high as might have been naively hoped for.

The new approach shows excellent convergence properties, with both energy and photoionisation cross-sections improving as O (d4), where d is radial grid resolution. The TDSE solution step scales strictly linearly with the total number of grid points used to represent the wavefunction. Overall memory requirements of the propagator are linear in the size of the wavefunction. The structure of the propagator is amenable to shared-memory parallel execution.

The combined numerical efficiency and accuracy of the new technique allows routine numerical simulation of strong-field processes on modest hardware (such as a desktop PC), without compromising numerical convergence of the results. We illustrate possible applications using examples of a high-harmonic spectrum in an elliptically-polarised field and a strong-field above-threshold ionisation in a linearly-polarised field. In either case, numerically converged results are achieved in less than an hour.

Apart from obvious applications in strong-field atomic dynamics, high numerical efficiency and modest hardware requirements of the present implementation of atomic TDSE opens additional possibilities. For example, several thousands of atomic TDSEs can be propagated in lock-step on a single floating-point accelerator card. This capability allows a straightforward implementation of the Maxwell-Schrödinger equation (MSE) for realistic systems, without resorting to supercomputing resources. In turn, routine availability of an MSE solver would enable accurate microscopic simulations of laser beam propagation, filamentation, and harmonic beam formation in dense atomic gases.

Acknowledgements

S.P. would like to acknowledge many helpful discussions with F. Morales-Moreno, gentle but persistent encouragement by M.Yu. Ivanov and O. Smirnova, and help in assigning photoelectron peaks in Fig. 3 generously offered by W. Becker and M.Yu. Ivanov.

References

[1] G.E. Kimball, G.H. Shortley, The numerical solution of Schrödinger's equation, Phys. Rev. 45 (1934) 815-820.

[2] P. Eckle, A.N. Pfeiffer, C. Cirelli, A. Staudte, R. Dörner, H.G. Muller, M. Büttiker, U. Keller, Attosecond ionization and tunneling delay time measurements in helium, Science 322 (5907) (2008) 1525-1529.

[3] P. Eckle, M. Smolarski, P. Schlup, J. Biegert, A. Staudte, M. Schöffler, H.G. Muller, R. Dörner, U. Keller, Attosecond angular streaking, Nat. Phys. 4 (7) (2008) 565-570.

[4] A.N. Pfeiffer, C. Cirelli, M. Smolarski, D. Dimitrovski, M. Abu-samha, L.B. Madsen, U. Keller, Attoclock reveals natural coordinates of the laser-induced tunnelling current flow in atoms, Nat. Phys. 8 (1) (2012) 76-80.

[5] J. Zhao, M. Lein, Determination of ionization and tunneling times in high-order harmonic generation, Phys. Rev. Lett. 111 (2013) 043901.

[6] I.A. Ivanov, A.S. Kheifets, Strong-field ionization of he by elliptically polarized light in attoclock configuration, Phys. Rev. A 89 (2014) 021402.

[7] M. Schultze, M. Fieß, N. Karpowicz, J. Gagnon, M. Korbman, M. Hofstetter, S. Neppl, A.L. Cavalieri, Y. Komninos, T. Mercouris, C.A. Nicolaides, R. Pazourek, S. Nagele, J. Feist, J. Burgdörfer, A.M. Azzeer, R. Ernstorfer, R. Kienberger, U. Kleineberg, E. Goulielmakis, F. Krausz, V.S. Yakovlev, Delay in photoemission, Science 328 (5986) (2010) 1658-1662.

[8] F. Mauger, A.D. Bandrauk, Electronic dynamics and frequency effects in circularly polarized strong-field physics, J. Phys. B 47 (19) (2014) 191001.

[9] G. Orlando, C.R. McDonald, N.H. Protik, G. Vampa, T. Brabec, Tunnelling time, what does it mean? J. Phys. B 47 (20) (2014) 204002.

[10] A. Maquet, J. Caillat, R. Ta'ieb, Attosecond delays in photoionization: time and quantum mechanics, J. Phys. B 47 (20) (2014) 204004.

[11] A.S. Landsman, U. Keller, Attosecond science and the tunnelling time problem, Phys. Rep. 547 (2015) 1-24.

[12] L. Torlina, F. Morales, J. Kaushal, I. Ivanov, A. Kheifets, A. Zielinski, A. Scrinzi, H.G. Muller, S. Sukiasyan, M. Ivanov, O. Smirnova, Interpreting attoclock measurements of tunnelling times, Nat. Phys. 11 (2015) 503-508.

[13] X.-M. Tong, S.-I. Chu, Theoretical study of multiple high-order harmonic generation by intense ultrashort pulsed laser fields: A new generalized pseudospectral time-dependent method, Chem. Phys. 217 (2-3) (1997) 119-130.

[14] M. Nurhuda, F.H.M. Faisal, Numerical solution of time-dependent Schrödinger equation for multiphoton processes: A matrix iterative method, Phys. Rev. A 60 (1999) 3125-3133.

[15] H.G. Muller, An efficient scheme for the time-dependent Schrödinger equation in the velocity gauge, Laser Phys. 9 (1) (1999) 138-148.

[16] A.G. Borisov, Solution of the radial Schrödinger equation in cylindrical and spherical coordinates by mapped Fourier transform algorithms, J. Chem. Phys. 114(18) (2001) 7770-7777.

[17] D. Bauer, P. Koval, Qprop: A Schrödinger-solver for intense laser-atom interaction, Comp. Phys. Comm. 174 (5) (2006) 396-421.

[18] X. Guan, C.J. Noble, O. Zatsarinny, K. Bartschat, B.I. Schneider, ALTDSE: An Arnoldi-Lanczos program to solve the time-dependent Schrödinger equation, Comp. Phys. Comm. 180 (12) (2009) 2401-2409.

[19] T. Sorevik, T. Birkeland, G. Oksa, Numerical solution of the 3D time dependent Schrödinger equation in spherical coordinates: Spectral basis and effects of split-operator technique, J. Comput. Appl. Math. 225 (1) (2009) 56-67.

[20] P. Botheron, B. Pons, Self-consistent bohmian description of strong field-driven electron dynamics, Phys. Rev. A 82 (2010) 021404R.

[21] A.N. Grum-Grzhimailo, B. Abeln, K. Bartschat, D. Weflen, T. Urness, Ionization of atomic hydrogen in strong infrared laser fields, Phys. Rev. A 81 (2010) 043408.

[22] X.-M. Tong, N. Toshima, Time-dependent method in the laser-atom interactions, Comp. Phys. Comm. 182 (1) (2011) 21-23.

[23] T. Dziubak, J. Matulewski, An object-oriented implementation of a solver of the time-dependent Schrödinger equation using the CUDA technology, Comp. Phys. Comm. 183 (3) (2012) 800-812.

[24] J. Shen, W.E. Sha, Z. Huang, M. Chen, X. Wu, High-order symplectic FDTD scheme for solving a time-dependent Schrödinger equation, Comp. Phys. Comm. 184 (3) (2013) 480-492.

[25] C. Ö Broin, L.A.A. Nikolopoulos, A GPGPU based program to solve the TDSE in intense laser fields through the finite difference approach, Comp. Phys. Comm. 185 (6) (2014)1791-1807.

[26] D.F. Gordon, B. Hafizi, A.S. Landsman, Amplitude flux, probability flux, and gauge invariance in the finite volume scheme for the Schrödinger equation, J. Comput. Phys. 280(2015)457-464.

[27] E. Lorin, S. Chelkowski, A. Bandrauk, A numerical Maxwell-Schrödinger model for intense laser-matter interaction and propagation, Comp. Phys. Comm. 177 (2007) 908-932.

[28] E. Lorin, S. Chelkowski, A.D. Bandrauk, Attosecond pulse generation from aligned molecules—dynamics and propagation in H+, New J. Phys. 10 (2) (2008) 025033.

[29] M.B. Gaarde, J.L. Tate, K.J. Schafer, Macroscopic aspects of attosecond pulse generation, J. Phys. B 41 (13) (2008) 132001.

[30] E. Lorin, S. Chelkowski, E. Zaoui, A. Bandrauk, Maxwell-Schrödinger-Plasma (MASP) model for laser-molecule interactions: Towards an understanding of filamentation with intense ultrashort pulses, Physica D 241 (12) (2012) 1059-1071.

[31] H.G. Muller, Numerical simulation of high-order above-threshold-ionization enhancement in argon, Phys. Rev. A 60 (1999) 1341-1350.

[32] H.G. Muller, Tunneling excitation to resonant states in helium as main source of superponderomotive photoelectrons in the tunneling regime, Phys. Rev. Lett. 83 (1999) 3158-3161.

[33] M. Boca, H.G. Muller, M. Gavrila, Dynamic stabilization of ground-state hydrogen in superintense circularly polarized laser pulses, J. Phys. B 37 (1) (2004) 147.

[34] M. Uiberacker, T. Uphues, M. Schultze, A.J. Verhoef, V. Yakovlev, M.F. Kling, J. Rauschenberger, N.M. Kabachnik, H. Schröder, M. Lezius, K.L. Kompa, H.G. Muller, M.J.J. Vrakking, S. Hendel, U. Kleineberg, U. Heinzmann, M. Drescher, F. Krausz, Attosecond real-time observation of electron tunnelling in atoms, Nature 446 (2007) 627-632.

[35] Y. Huismans, A. Rouzee, A. Gijsbertsen, J.H. Jungmann, A.S. Smolkowska, P.S.W.M. Logman, F. Lepine, C. Cauchy, S. Zamith, T. Marchenko, J.M. Bakker, G. Berden, B. Redlich, A.F.G. van der Meer, H.G. Muller, W. Vermin, K.J. Schafer, M. Spanner, M.Y. Ivanov, O. Smirnova, D. Bauer, S.V. Popruzhenko, M.J.J. Vrakking, Time-resolved holography with photoelectrons, Science 331 (6013) (20l1) 61-64.

[36] T. Marchenko, H.G. Muller, K.J. Schafer, M.J.J. Vrakking, Electron angular distributions in near-threshold atomic ionization, J. Phys. B 43 (9) (2010) 095601.

[37] L Torlina, F. Morales, H.G. Muller, O. Smirnova, Ab initio verification of the analytical R-matrix theory for strong field ionization, J. Phys. B 47 (20) (2014) 204021.

[38] P.D. Lax, Linear Algebra, Wiley-Interscience, New York, 1997.

[39] Wolfram Research, Inc., Mathematica, version 10.0 Edition, Wolfram Research, Inc., Champaign, Illinois, 2014.

[40] G.B. Arfken, H.J. Weber, F.E. Harris, Mathematical Methods for Physicists, seventh ed., Academic Press, 2013.

[41] L Landau, E.M. Lifshitz, Quantum Mechanics: Non-relativistic Theory, fifth ed., Fizmatlit, Moscow, 2002.

[42] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes in Fortran 77: The Art of Scientific Computing, second ed., Cambridge University Press, Cambridge, 1992.

[43] D.E. Manolopoulos, Derivation and reflection properties of a transmission-free absorbing potential, J. Chem. Phys. 117 (2002) 9552-9559.

[44] A. Burgess, Tables of hydrogenic photoionization cross-sections and recombination coefficients, Memoirs Roy. Astr. Soc. 69 (1965) 1-17.

[45] L Tao, A. Scrinzi, Photo-electron momentum spectra from minimal volumes: the time-dependent surface flux method, New J. Phys. 14 (2012) 013021.

[46] E. Cormier, P. Lambropoulos, Effect of the initial phase of the field in ionization by ultrashort laser pulses, Eur. Phys. J. D 2 (1998) 15-20.

[47] M.Y. Ivanov, T. Brabec, N. Burnett, Coulomb corrections and polarization effects in high-intensity high-harmonic, Phys. Rev. A 54 (1996) 742-745.

[48] T. Nubbemeyer, K. Gorling, A. Saenz, U. Eichmann, W. Sandner, Strong-field tunneling without ionization, Phys. Rev. Lett. 101 (2008) 233001.

[49] F. Krausz, M. Ivanov, Attosecond physics, Rev. Modern Phys. 81 (2009) 164-234.

[50] S.V. Popruzhenko, Keldysh theory of strong field ionization: history, applications, difficulties and perspectives, J. Phys. B 47 (2014) 204001.

[51] A. Kästner, U. Saalmann, J.M. Rost, Electron-energy bunching in laser-driven soft recollisions, Phys. Rev. Lett. 108 (2012) 033201.