The quantum-mechanical Coulomb propagator in an L2 function representation

The quantum-mechanical Coulomb propagator is represented in a square-integrable basis of Sturmian functions. Herein, the Stieltjes integral containing the Coulomb spectral function as a weight is evaluated. The Coulomb propagator generally consists of two parts. The sum of the discrete part of the spectrum is extrapolated numerically, while three integration procedures are applied to the continuum part of the oscillating integral: the Gauss–Pollaczek quadrature, the Gauss–Legendre quadrature along the real axis, and a transformation into a contour integral in the complex plane with the subsequent Gauss–Legendre quadrature. Using the contour integral, the Coulomb propagator can be calculated very accurately from an L\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^2$$\end{document}2 basis. Using the three-term recursion relation of the Pollaczek polynomials, an effective algorithm is herein presented to reduce the number of integrations. Numerical results are presented and discussed for all procedures.

The long-range Coulomb interaction plays a central role in quantum-mechanical scattering processes with charged particles. Ultra-short laser pulses in the range of femto-and attoseconds allow processes to be monitored in atomic-scales time domains; in particular the final state can be monitored after turning off the laser field, such that the movement of the photoelectron in the field of the charged ion is considered. Phenomena such as above treshold ionization, effects due to carrier envelope phase, and multi-photon ionization, in which an atom or molecule is exposed to an intensive laser pulse for a very short time, require theoretical provision of respective methods to describe the time development of a quantum-mechanical system incorporating Coulomb interactions.
A theoretical treatment of atomic and molecular Coulomb scattering processes, based exclusively on squareintegrable functions (L 2 ), is presented in J-Matrix theory. This can be interpreted as a square-integrable analog to the time-independent R-Matrix theory, in which the scattering wave function is expanded in an L 2 function space. In recent years, the J-Matrix theory has found widespread significance and applications. Two examples are Konovalov and Bray 1,2 , who calculated electron-impact ionization, as well as the extension to relativistic scattering processes by Alhaidari et al. 3 . For an overview of current applications, we direct the reader to Alhaidari et al. 4 .
The basic functions of the J-Matrix scattering method are the Coulomb Sturmian functions, which have been applied to solve ab initio the time-dependent Schrödinger equation in the atto-and femtosecond range via close-coupling equation systems. Madronero and Piraux 5 and Hamido et al. 6 used this approach to calculate the ionization of atomic hydrogen and helium in a femtosecond laser field.
The aforementioned works show that L 2 functions can be successfully used to determine the time development of a system under the influence of a Coulomb field.
For this purpose, we consider the L 2 representation of the Feynman propagator 7,8 for a two-particle Coulomb system: the propagator contains all necessary information to determine the time evolution of a particle in a Coulomb field, i.e., its wave function at time t ′ can be calculated using knowledge of the wave function at time t.
The time propagator can be written formally as an integral and a sum over the eigenstates of the stationary Schrödinger equation and the Coulomb spectral function, bearing in mind that the integrand contains the oscillating function e −iE(t ′ −t) .
(1) K + (r ′ , t ′ ; r, t) = θ(t ′ − t) www.nature.com/scientificreports/ E characterizes an eigenfunction of positive energy, E n characterizes a bound state, and θ(t ′ − t) is the Theta function taking a value of 1 for t ′ > t ; the superior + symbol shows the retardation realized by the Theta function and the superior star symbolizes a complex conjugation.
For the Coulomb propagator, a range of approximations valid in the limit case of ℏ → 0 ( ℏ represents the Planck constant) were developed in the semiclassical limit; representative works on this subject include Gutzwiller 9 , Tomsovic and Heller 10 and Manning and Ezra 11 . Unlike these semiclassical approaches, the complete quantum-mechanical Coulomb propagator without approximations is considered here.
In a number of papers [12][13][14][15] , it has been shown that not only the bound states (with negative energy), but also the continuum states, can be developed using a single function system of the L 2 functions. These solutions to the stationary Schrödinger equation in a L 2 basis can be used for the L 2 representation of the Coulomb propagator.
Although the Coulomb-Greens operator-which is linked to the Coulomb propagator via a Fourier transformation-can be applied by default in scattering computations, the square-integrable analog of the Coulomb propagator has not been studied. To the best of our knowledge, there is no literature on dealing with the Coulomb propagator on an L 2 basis.
The goal of this paper is to evaluate numerical methods for the accurate calculation of the Coulomb propagator in an L 2 basis and to understand its convergence behavior. J-Matrix-based applications in atomic physics, such as electron-atom-scattering and photoionization or nucleon-nucleon scattering and nuclear decay, can benefit from the results presented herein 4 .
For the expansion of the Coulomb propagator, hydrogen-like Sturmian functions 16 are used; these contain an additional, often freely available, parameter that increases the flexibility of the system, in particular for applications where various angular momentum states are involved. Closer consideration of (1) reveals that integration over the continuous spectrum poses a great challenge, since a strongly oscillating function occurs in the integrand. The situation is more complicated due to semi-infinite integration intervals and the infinite sum over bound states.
In previous literature, various methods and libraries have been proposed to numerically solve such integrals with oscillating functions. However, if the integrand falls sufficiently "quickly", standard packages from Matlab, IMSL, Maple, and NAG can be used. Many of these packages trace back to Quadpack 17 , applying a Clenshaw-Curtis algorithm with adaptive partitioning, such that the line of partial sums can be extrapolated with the ǫ algorithm 18 . Evans 19 , Haider and Liu 20 and Sauter 21 also conducted partitioning of the integration interval, as well as Sidi 22 who partitioned according to the zeros of the integrand before calculating the line of partial sums using a convergence accelerator.
Oscillatory integrals often appear in quantum scattering and in wave propagation in classical electrodynamics. There is a great interest and need for techniques that can cope with such integrals. Besides the aforementioned partitioning approach, methods for evaluating such integrals can be divided into three categories: asymptotic expansions of the integrand and Filon-and Levin-type approaches.
In the asymptotic expansion method the integral is approximated by a series either through repeated integration by parts or by employing a formal expansion of the integrand and then evaluating the integrals term by term. In the Filon method the integral dxf(x) e ig(x) is partitioned into N parts and in each part f(x) is approximated by a polynomial. The resulting parts can be integrated either analytically or by standard procedures. The Levin method searches for a function F such that d dx [F(x) e ig(x) ] = f(x) e ig(x) . This leads to dxf(x) e ig(x) = dx d dx [F(x) e ig(x) ] , which can be evaluated trivially. The task to work out the integral is shifted to solve a differential equation for the unknown function F.
A survey of these methods can be found in Deaño, Huybrech and Iserles 23 , where further work on the topic is also referenced, up to and including 2017. Oscillatory integrals are the subject of ongoing research: Yang and Ma 24 , Zaman et al. 25 used Levin-based approaches to calculate highly oscillatory Fourier integrals in one-and two-dimensional domains, whereas Wang und Xiang 26 applied a Levin method to singular integrands. Recent articles from Kayijuka et al. 27 dealt with oscillatory Fourier integrals with singular integrands, while Zaman et al. 28 used a Levin based approach to evaluate the integrals over Bessel functions. For integrals containing Bessel functions 29 , alternative methods are sometimes expedient: for example the transformation to a contour integral in the complex plane, which is then exponentially damped in an asymptotic manner. The previously mentioned path is pursued by this work in order to numerically calculate the oscillating integral in the Coulomb propagator.
In the following method sections the Coulomb propagator is calculated according to three procedures: 1. Gauss-Pollaczek quadrature along the real axis 2. Gauss-Legendre quadrature along the real axis 3. Transformation into a complex contour integral with exponential damping and numerical calculation using a Gauss-Legendre quadrature.
The paper is structured as follows: first the Coulomb propagator in the L 2 representation of Sturmian functions is reviewed. Then the the Gauss-Pollaczek quadrature is explained and numerically stable methods for the calculation of the Pollaczek polynomials are shown. Thereafter the transformation of the continuum part into a contour integral in the complex plane is explained and the calculation of the discrete part of the Coulomb propagator is presented. The result section then demonstrates a comparative analysis of the three integration methods. In the last section potential applications are discussed.

Coulomb propagator in the L 2 representation
The Coulomb propagator K(r ′ , t ′ ; r, t) solves the time-dependent Schrödinger equation using a δ-function as the inhomogeneity, while r symbolizes the space coordinate and t time. where z symbolizes the effective nuclear charge. Since the Hamilton operator H 0 is time-independent, the Coulomb propagator is only a function of t ′ − t in the time part. Thus, K(r ′ , t ′ ; r, t) is also a solution for the following: In the time domain K(r ′ , t ′ ; r, t) satisfies the initial condition: The probability of the system to be in state �(r ′ , t ′ ) -under the condition described by �(r, t i ) at time t i -can therefore be depicted as an integral through all possible paths throughout the configuration space.
For the Coulomb propagator, a partial wave expansion can be conducted in accordance with angular momentum eigenstates; the physical solution of the stationary radial Schrödinger equation to the angular momentum l and energy E = k 2 /2 is as follows: Here, M(a, b, x) is a confluent hypergeometric function and J + l (k) is the Jost function, which is given by: A factor ( 1 2kπ ) 1 2 occurs due to normalization to an energy delta function as follows: The root of the Jost function is given by the Gamma function if the argument l + 1 − i z/k results in a negative integer.
The Coulomb spectral function is defined as the inverse of the square of the absolute value of the Jost function: Consequently, the Coulomb propagator can be depicted as a Stieltjes integral via the Coulomb spectral function.
The Y lm within are spherical harmonics that encapsulate the angular part. Here, the bound states occur at poles in the spectral density at E n b = −z 2 /2n 2 b , where n b = n r + l + 1 is the full quantum number and n r ∈ [0, 1, . . .] is the radial quantum number. Then the residue of the integrand at each bound-state energy is given by with L 2l+1 n r a Laguerre polynomial. As shown in 15 , the origin regular solution � + l (r, E) to positive energy can be developed on an L 2 basis as follows: with square-integrable functions φ l n given by: The variable is a scaling parameter obeying the restriction > 2z/(l + 1) . Larger values of result in maxima of φ l n that are further from the origin. The great advantage of the expansion (16) lies in the factorization of + l in an energy independent spatial part φ l n and a kinetic part, exclusively contained in the expansion coefficients l+ 0 l n . These are representative in an analytical way and consist of the Pollaczek polynomials p l n .
The function l+ 0 is independent of n and given by: and γ , ξ and x are defined by An expansion of the irregular Coulomb function scaled at the origin similar to r −l−1 cannot be realized using the L 2 functions φ l n . However, an expansion according to φ l n and showing the same asymptotic behavior as the irregular Coulomb function in the limit r → ∞ , can be realized. This finite solution at the origin-behaving like the irregular Coulomb function for large r and being normalized to an energy delta function-is given by: with the expansion coefficients Q l+ n , which themselves consist of the Pollaczek functions q l+ n 15 : ] and the insertion of (16), (18), (19), and (10) into (13), results in the L 2 expansion of the Coulomb propagator according to the basis functions φ l n : Upon closer inspection, the endpoint singularity occurring in the integrand at x = 1 is of shape (1 − x) l−1/2 and can be eliminated through transformation to variable k, which is linked to x = (k 2 − 2 4 )/(k 2 + 2 4 ) . (In the following the integration bounds for the x integration will be omitted, in all formulas it is implicitely assumed For a wave function �(r, t) , an expansion into an orthonormalized function system can be pursued (for simplicity, index m is omitted in the following equations): Thus, �(r, t) can be expanded in F l n with time-dependent expansion coefficients α l n (t): The insertion of (25) into the conditional equation (6) and a projection with F l n (r ) then leads to In (26), the integration over x also includes the sum over the bound-state spectrum. Weight ρ l (x) is given by with poles when the argument of the Ŵ function is either zero or a negative integer, and N l n by The integration in the configuration space of matrix elements �φ l n ′′ |f l n � and �φ l n ′ |f l k � can be achieved analytically and is only non zero for n ′′ = n, n + 1 and n ′ = k, k + 1 . Consequently, the sum over n ′′ , n ′ reduces to only a few terms when using the expansion (25). A remaining challenge is to solve the oscillating integrals: and to determine the infinite sum over the bound states. In the second line of (29) the first term describes the continuum part. The discrete part over the bound states can be written in accordance with (14) as the negative sum over all residues of the weight function ρ l evaluated at the poles of the Ŵ function.
For n ′′ , n ′ ∈ [0, N − 1] the number of integrations scales with N 2 for each angular momentum l, and can be drastically reduced by exploiting the three-term recursion relation (see 32) for the Pollaczek polynomials: with b l n = −2(n + l + 1 − 2z/ ) This recursion directly applies to (29). For calculation of the square matrix I l n ′ n for n ′ , n ∈ [0, N − 1] , it is necessary to compute the 2 N-1 start values I l 0n in n, which recur with (30) upwards in n ′ . First test cases give evidence that the recursion is stable. This is further supported by the fact that the diagonal elements I l nn , which are greater than the off-diagonal elements, are constructed from much smaller values I l 0n . In the following sections, we investigate integrals I l n ′ n without restriction for integer values n, n ′ ≥ 0 . Integration via the continuum is examined through three different methods: the Gauss-Pollaczek and Gauss-Legendre quadratures along the real axis and calculations with a rotation into the complex plane, resulting in exponential damping. For this purpose, the Pollaczek polynomials in the Coulomb propagator are discussed and their numerical calculation is introduced in greater detail in the following section.

Gauss-Pollaczek quadrature and calculation of the Pollaczek polynomials
The Pollaczek polynomials are real for real arguments and form an orthogonal function system with regard to the weight function ρ l (x) specified in (27):  (29) differs from (31) only in terms of the factor e −iE(t ′ −t) /(1 − x) ; therefore, it can be assumed that a Gauss-Pollaczek quadrature is well-suited for calculating (29); this hypothesis is examined in more detail in the following sections.
The Pollaczek polynomials obey, similarly to orthogonal polynomials, a three-term recursion relationship as follows: The orthogonality relationship of the Pollaczek polynomials enables the development of an N-Point Gauss quadrature that exactly determines integrals whose integrands contains the Coulomb spectral function as follows: x j , j = 1, . . . , N are the N roots of the Pollaczek polynomial p l N and w j , j = 1, . . . , N , with corresponding weights. The great advantage of the Gauss Pollaczek quadrature is its inherent feature to include the infinite bound state sum directly in the quadrature method, thus avoiding the separate evaluation. This is also reflected by the fact, that the zeros of the Pollaczek polynomials occur at positive as well as at negative energies. Therefore the sum on the right side of (33) approximates both parts-discrete and continuum.
Both the weights and the roots can be calculated very easily via diagonalization of a tridiagonal matrix whose coefficients are given by the monic recursion relation of the (orthonormal) Pollaczek polynomials: The weights in (33) result from the first component of the ith eigenvector v 1i of a tridiagonal matrix together with the zeroth moment, which are thus given by: At this point, the following question arises: how well does this Gauss-Pollaczek quadrature work for integrals of the form (29) containing terms e −iE(t ′ −t) /(1 − x) , which oscillate strongly? This question is examined and answered in more detail in the following section.
For the analysis of (29), the Pollaczek polynomials must be calculated for both negative discrete and positive values of E in the interval [0, ∞].
The sum over the bound-state spectrum-here, the sum over all residues of the integrand-occurs only when parameter z-the atomic number-is positive. Residues result from the poles of the Gamma function when the argument l + 1 − iz/k is a negative integer and characterize the bound states at E n b = −z 2 /2n b 2 . If parameter is chosen in such a way that = 2z/(n r + l + 1) is valid, then φ l n r (r, ) is directly proportional to a bound eigenstate. This value of can thus be used to optimize the wave function of a particular eigenstate.
As stated in 15 , the three-term recursion relationship (32) can be used as a starting point for calculation. The division of all terms in (32) by p l n setting p l n+1 = r l n+1 p l n with 1/r l 0 := 0 leads to the following recursion: As shown by calculations, this algorithm is numerically stable for positive energy values x ∈ [−1, +1] and any values of .
For the discrete part corresponding to negative energies E n b = −z 2 /2n 2 b and → 2z/(n r + l + 1) , however, the recursion (36) becomes numerically unstable for r l n+1 ; these values, however, are needed for (29). At this point, it is best to express the Pollaczek polynomials via a hypergeometric function and implement their evaluation via the Horner scheme as follows: The 2 F 1 reduces to a polynomial in 1 − 1/ξ 2 due to a negative integer in the first argument. The Horner scheme is a fundamental algorithm to evaluate polynomials with a minimum of arithmetic operations and is numerically stable. This makes the algorithm best suited for computer implementation.
For the bound states E n b = −z 2 /2n b 2 , k causes purely imaginary and p l n results for the following finite sum:  Table 1, the values of p 12 14 (x) with E n b = − 1 2n b 2 , n b = 13 were calculated in an exemplary way for different scaling parameter and with varying accuracy. The values are given in columns 3 and 4: this reflects the calculation with the three-term recursion formula (36) and by calculating the sum (38) according to the Horner scheme.
Implementation is achieved using a numerical package (Java Apfloat) that allows for calculating with userdefined accuracy.
In particular, for small values of (for → 2z/(l + 1) ), it can be seen that the three-term recursion (36) becomes numerically unstable; the number of digits required for calculation is significantly higher if the same result is sought as with (38). In contrast, calculations using the Horner scheme are stable for all tested values of the scaling parameter.
Consequently, the following assertions can be made: 1. For positive values of E or equivalent x ∈ [−1, +1] , the three-term recursion formula (36) is numerically stable. 2. The same applies to bound energies E n b = − z 2 2n b 2 under the condition that differs from 2z/(l + 1) in the first two leading digits. 3. For → 2z/(l + 1) and E n b = − z 2 2n b 2 , the calculation with the Horner scheme is numerically stable.

Calculation of the continuum part-contour representation
In this section of the paper, the Stieltjes integral part is considered across the positive energy spectrum. To make a clear dsitinction between the whole Stieltjes integral (29)-which includes the bound state sum-the symbol A l nn ′ is used to indicate the continuum part: www.nature.com/scientificreports/ For this purpose, a transformation is conducted in the following subsection, resulting in exponential damping. Furthermore a contour representation based on a conformal mapping to the unit circle is explained. Exponential damping. Using the connection between the kinetic energy and the absolut value of the impuls, E = k 2 /2 , and transformation to the variable k given by x = (E− 2 /8)/(E+ 2 /8) = (k 2 − 2 4 )/(k 2 + 2 4 ) then leads to the following integral: Here, too, the integrand strongly oscillates. An asymptotic expansion of the integrand identifies the following behavior: This demonstrates that the integrand for large angular momentum numbers l is damped to a much greater extent than for small values of l. The integrand oscillates widely, especially for large values of k. If traditional integration routines are applied, a more accurate numerical evaluation can only be expected with the use of many abscissae, especially for smaller l values.
Consequently, it is necessary to reexpress the integral (39) using a contour integral in the complex plane such that, for the resulting path, the integral is much more strongly damped for large values of k than for integration along the real axis. This is achieved by transformation to the variable y: Thus, the integrand is "practically" free of oscillation with exponential damping. The integration path in the complex plane is shown in Fig. 1 and is characterized by a closed path in the lower complex half plane.
In the limiting case of k → ∞ , the contribution to the integral vanishes along the path C 2 , resulting in the integral I l nn ′ being determined by the path along C 3 . The integral can now be evaluated, for example, with a Gauss-Legendre quadrature, although due to exponential damping, fewer abscissae are necessary than in the integration along the real axis. It is expected that the number of abscissae decreases at higher values of (t ′ − t).
Vice versa, it can be assumed that, for small values of (t ′ − t) ≈ 0 , the number of abscissae will increase in order to maintain a certain accuracy. However, the numerical input for the calculation with the contour representation in the complex plane is significantly higher since, in such cases, the integrand is evaluated using complex arguments.
Transformation to unit circle. By transforming to the variable ξ = ( + 2ik)/( − 2ik) , a conformal mapping is defined with ξξ * = 1 . Consequently, previous work has demonstrated the orthogonality of the Pollaczek polynomials. When using their analytical properties 15 and carrying out the transformation onto the variable ξ , the integral A l nn ′ results in the following: In the above, q l+ n is the coefficient of the irregular Coulomb function-as shown in (21)-and the integration path is a closed curve over the unit circle in the complex ξ plane, as shown in Fig. 2. Within this closed curve lie the bound states with energies E = −z 2 /(2n b 2 ) on the real ξ axis; here, the function q l+ n has poles, while p l n is free of poles. The integrand has another, higher-order pole at digit ξ = −1 due to the argument of the exponential function appearing in the integrand, lying exactly on the integration path, while the poles of q l+ n are located within the unit circle. This contour representation can therefore only be consulted for analysis in the case of (t ′ − t) = 0 . According to this behavior, use of the Cauchy principle value and application of the residual theorem enables the integral over the continuum part to be calculated analytically: where n < specifies the min ( n, n ′ ). The first term arises from the pole at ξ = −1 and the infinite sum over the residues result from the simple bound state poles lying within the closed contour of the unit circle. This sum cancels exactly with the bound states sum in the Coulomb Propagator (29), such that the result for (t ′ − t) = 0 is as follows:  (44) only appears in the case of z > 0 and extends across all discrete eigenvalues with a negative value. The sum is convergent, whereas the individual terms are scaled with 1/n b 3 , which is investigated in more detail in the following section.

Computation of the sum via the discrete spectrum
With analogy to (14), the result for the sum occurring via the discrete spectrum in (29) with (27) is as follows: The function ρ l (x) has poles that are given by the Ŵ function Ŵ(l+1−iz/k) ; these occur if the argument results in 0 or a negative number: Using the relationships Ŵ(x)Ŵ(1 − x) = π/ sin(πx) 31 and sin γ = (ξ − 1/ξ )/2i , the residuum occurring in (46) can be calculated. If transformed to the variable β = iz/k , this results in: and evaluation at β = n b in: where the Pollaczek polynomials are evaluated at argument The sum then scales following 1/n b 3 and exhibits convergence. Since the limit of the function for n b → ∞ cannot be crossed, computers for numerical evaluation aborted the sum at a certain value. The best choice to yield higher accuracy is to calculate a sequence of partial sums (S 1 , S 2 , . . . S N ) and extrapolate these to a limit value S ∞ . This is the basic idea of the Aitken-Neville method: It is an algorithm which uses the partial sum values as input and thus realizes a convergence acceleration.
The algorithm itself is defined by the Aitken-Neville tableau: When defining partial sums and the variable h k as the Aitken-Neville tableau is given as follows: Here, represents a user-defined integer value; Table 2 provides an example of the approximate value of the sum (48), with a specified p l n ≡ p 4 5 and p l n ′ ≡ p 4 2 and t ′ − t = 1 -once calculated directly, whereas the sum was aborted at a upper limit ( n max ) (column 1), the second time was extrapolated via an Aitken-Neville algorithm (column 5). For the extrapolation in column 5, the same upper summation limit was chosen such that columns 2 and 5 are directly comparable ( N = n max ). Our calculations were carried out using three procedures discussed in previous sections: 1. Gauss-Pollaczek quadrature along the real axis. The respective abscissae and weights result from diagonalization of the coefficient matrix with a tridiagonal structure as shown by (34). It is not required for the weight function ρ l (x) to be calculated explicitly, since it is implicitly part of the weight w i ; likewise, this incorporates the sum over the discrete spectrum into the Gauss quadrature, making its separate calculation obsolete. This results in 2. Gauss-Legendre quadrature along the real axis. For the calculation of the abscissae and weights, standard routines are used to calculate the roots of the Legendre polynomial with the Newton procedure. 3. Complex contour integration. In order to enable direct comparison with real axis integration, we chose the same Gauss-Legendre quadrature method.
Specifically, only 2 and 3 can be directly compared to each other, because the Gauss-Pollaczek quadrature approximates the complete term (29), including the sum over the discrete negative bound energies. Thus, for the calculation along the real axis and complex contour, the sum in (52) was calculated following the procedure described in the previous section and added to the integral over the continuum part in all Tables 3, 4 and 5. For this calculation, a Java package (Apfloat) that allows calculations with any user-defined number of digits was used. Calculations were carried out with 40 digits to prevent rounding and truncation errors influencing the result. A particularly interesting aspect of these three integration procedures is the case l = 0 . In this case, the convergence behavior should be least pronounced.
As a result, the following considerations can be given.
1. Convergence behavior improves in all three procedures with increasing l values (see Tables 3,4,5). This is to be expected due to the factor 1/k 2l+2 , and substantiates the qualitative statement from the previous section. 2. The convergence behavior of calculations with complex contour representations improves at higher values of (t ′ − t) ; in the other two procedures, the opposite is valid. 3. If considering the value l = 0 , the Gauss-Pollaczek procedure even converges poorly at t ′ − t ∼ 0 , although no oscillating behavior occurs and the integrand differs from the orthogonality relationship of the Pollaczek polynomials only by the factor 1/ (1 − x) . Even with 1200 abscissae, the result is accurate for only around (53)   www.nature.com/scientificreports/ three digits, i.e., an unexpected finding. At higher values of t ′ − t , the convergence behavior worsens-at this point, the oscillating character of the exponential function comes into play. 4. Gauss-Legendre quadrature along the real axis: Here, an improved convergence behavior occurs, with fewer abscissae than in the Gauss-Pollaczek quadrature. The reason for this probably lies in the distribution of the abscissae. In the Gauss-Pollaczek procedure, they have an accumulation point at E = 0 , while they are more evenly distributed in the Gauss-Legendre procedure. 5. When considering the Gauss-Legendre quadrature of the complex contour integral, the number of abscissae is highest for t ′ − t ∼ 0 to achieve a satisfactory result, whereas that for the value t ′ − t ∼ 0 observed is only slightly worse than in the integration along the real axis. 6. In general, complex contour representation, combined with a Gauss-Legendre quadrature, yields the best result; however, the numerical input is high and the integrand must be evaluated using complex arguments.

Conclusion and prospects
The evaluation of the Coulomb propagator on a L 2 basis can be accurately computed using extrapolation techniques for the bound state sum and mapping of the continuous part to a contour integral in the complex plane with subsequent Gauss quadrature. Direct comparison of integrations along the real axis and contour, requires that the same quadrature method is applied. Test cases give evidence that, for contour integration, the numerical method can be further improved by using a Gauss Hermite quadrature, thereby reducing the number of abscissae. Numerical efforts to evaluate contour integrals decrease with increasing time t − t ′ and angular momentum l. Increasing t − t ′ results in increasing damping, thus lowering the number of abscissae and the number of digits to achieve numerical accuracy. The same is true for increasing angular values l, there the integrand behaves like 1/k 2l+2 . For practical applications, N 2 /2 integrals were evaluated, where N is the maximum number of basis functions included in the expansion. Using the three-term recursion relation for the Pollaczek polynomials, the number of evaluations reduces to 2N.
The methods shown in this article can be used, for example, to propagate an electronic wavepacket released after a few cycles light pulse or in time-dependent perturbation theory. The computational effort is simplified when the initial state is a free state, i.e., it can be decomposed in a superposition of positive energy components, where the sum over bound states vanishes.
The L 2 -Coulomb propagator can also be applied to draft a time-dependent approach of the J-Matrix scattering method. This leads to a L 2 analog of time-dependent R-Matrix theory, wherein the function space is split into an inner and an outer space. As such, time-dependent calculations on the femto-and attosecond photoionization range could be realized for one-electron systems.  Table 3. Results for angular momentum l = 0. Numerical values of the Integral (52) dx ρ l (x) 1 − x p l 5 p l 2 e −iE(x)(t ′ −t) with l = 0 for different values of t ′ − t , (z=1, = 2.01), evaluated using theGauss-Pollaczek quadrature (column 3), Gauss-Legendre quadrature (column 4) along the real axis and as a contour integral with the Gauss-Legendre quadrature (column 5) and, in column 4 and 5, the value of the sum over the bound states is added to enable comparison with column 3.  Table 4. Results for angular momentum l = 1. Same as Table 3 but with l = 1. Table 5. Results for angular momentum l = 4. Same as Table 3 but with l = 4 and ( z = 1, = 0.41).