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

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 summaryProgram title: SCID-TDSE: Time-dependent solution of 1-electron atomic Schrödinger equation in strong laser fields.Catalogue identifier: AEYM_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AEYM_v1_0.htmlProgram obtainable from: CPC Program Library, Queen’s University, Belfast, N. IrelandLicensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.htmlNo. of lines in distributed program, including test data, etc.: 334254No. of bytes in distributed program, including test data, etc.: 20005596Distribution format: tar.gzProgramming 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−MMIN+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∗NR2∗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.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 O(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 O(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 O(NR3) and O(NR3∗(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).