How Novel Algorithms and Access to High Performance Computing Platforms are Enabling Scientific Progress in Atomic and Molecular Physics

Over the past 40 years there has been remarkable progress in the quantitative treatment of complex many-body problems in atomic and molecular physics (AMP). This has happened as a consequence of the development of new and powerful numerical methods, translating these algorithms into practical software and the associated evolution of powerful computing platforms ranging from desktops to high performance computational instruments capable of massively parallel computation. We are taking the opportunity afforded by this CCP2015 to review computational progress in scattering theory and the interaction of strong electromagnetic fields with atomic and molecular systems from the early 1960’s until the present time to show how these advances have revealed a remarkable array of interesting and in many cases unexpected features. The article is by no means complete and certainly reflects the views and experiences of the author.


Introduction
The time independent and time dependent many-body Schrodinger equations show a remarkable array of features which often only emerge subsequent to laborious and time consuming calculations. Progress in the field has been marked by the development of novel theoretical and numerical methods, backed by solid algorithms, software development and the continuing evolution of computing platforms ranging from powerful desktops to supercomputers. In this article, I outline historical progress in the treatment of atomic and molecular continua in scattering processes and how these ideas have been extended to the time-dependent interaction of short, intense laser pulses with such targets.

Early History: The 1960's and 1970's
I begin with a short history of early developments in electron-atom/molecule collisions.

i.
Computers advance sufficiently that quite high quality calculations could be performed on the scattering of electrons and photoionization of simple atoms such as H, He and Ne. ii.
Significant experimental advances; A number of experimental groups in AMO physics develop instruments able to extract high quality total and differential electron scattering cross sections and photoabsorbtion cross sections ( See Figure  1 [1,3]) Densitometer trace of the absorption spectrum of He. The absorption is enhanced on the low-frequency side and reduced on the high-frequency side. See [1] Photoabsorption of H − . Comparison of theory (solid curve) [See [2]] and experiment (with error bars) [See [3]]. Data from experiment normalized to theory at 10.90 eV.
(iii) Computation (See 3.1 for discussion of close-coupling (CC)) and experiment discover interesting and in some cases unexpected features such as Feshbach and shape resonances in an otherwise smooth background continuum [4,5]. Such long-lived, quasi-bound states were seen in nuclear scattering earlier.

(iv)
Extension of CC to a number of first row atoms (v) An intuitive but useful approach called the Stabilization method [6]

iii.
Scattering methods were extended to the time-dependent domain where new experimental facilities such as free electron lasers are becoming available to study the interactions of short, intense laser radiation with atomic and molecular targets. Now possible to "see the electrons dance" and even control their motion using optical techniques. Looking directly at electron correlation effects becomes a reality computationally and experimentally

Personal View
In my opinion, the most effective approaches to treat problems in atomic and molecular systems involving a continuous spectrum, have been based on techniques which try to avoid the details of the asymptotic form of the wavefunction. There are a few options. One may divide either physical or function space into two ( or perhaps more ) regions. In the innermost region, one needs to deal with the full interacting set of nuclei and electrons but does not have to worry about the details of asymptopia. This inner region is then quite similar to a bound state problem and it becomes unsurprising that the most successful computational approaches to scattering borrow liberally from those developed by the atomic structure and quantum chemistry communities for discrete states. In the other regions, the Hamilitonian often "simplifies" in some way that makes it computationally easier to solve.
In real space, one might be able to set the interaction to zero or to some problem where an analytic solution is available 3.3. In function space, one could assume the expansion coefficients of the asymptotic functions are known outside the interaction region 3.4. How the different regions are joined, and how the scattering information is extracted becomes the central difference of all the various approaches. There is also the possibility of avoiding the asymptotic region entirely by analytically continuing from the complex plane, rotating the scattering coordinates into the complex plane, or relating a discrete basis set expansion to a numerical quadrature. Among the first class of methods are the R-and-J matrix method, and the Kohn and Schwinger variational principles. The latter category is exemplified by either full or exterior complex scaling and the equivalent quadrature technique. Lets us look briefly at a few of these.

Close Coupling
The close-coupling method [5] proceeds in the most straightforward manner; expand the scattering wavefunction in a set of target states times unknown one-electron scattering functions, insert into the Schrodinger equation, project onto the target states to obtain a set of coupled integro-differential equations and solve those equations by numerical integration. The close coupling or expansion in target states approach became routinely employed over the 1960's-1970's to both neutral and ionized species. The computed cross sections were used in astrophysical and fusion modeling. Computing photoionization cross sections became a simple matter of calculating the dipole matrix element between the electron ion wavefunction and the initial bound state, recognizing the dipole selections rules, a nice benefit once the scattering calculation has been performed. The method was extended to include pseudostates to account for the polarizability of the true target states and to represent the ionization continuum. However, it also became evident that the direct solution of the integrodifferential equations of scattering theory, was not the most efficient way to attack the problem numerically and other approaches were initiated.

Stabilization Method
The "Stabilization Method", was developed by A. Hazi and H. Taylor [6] ( See Figure 2 ) at USC. The success of the approach had a significant influence on much of the work that followed. Taylor correctly and intuitively argued that a discrete basis set expansion would converge, inside the range of the basis, to the true scattering function within an unknown normalization constant. By diagonalizing the Hamiltonian in a larger and larger discrete basis, certain eigenvalues "stabilized" with basis set size. Those eigenvalues that stabilized, were the resonances, and they stabilized simply due to the fact that resonant wavefunctions tend to be localized with rather small asymptotic tails. It was also noticed that the basis could describe the background continuum, but again only to an unknown normalization factor. There are a number of approaches which utilized the stabilized solution to find estimates of the width. Some are motivated by the connection between the stabilized eigenfunction and the variational method proposed by Harris [7]. Another method used the diagonalization to find an approximation to the background continuum function and then computed the width via the Fermi Golden Rule and an approximation to the density of states suggested by Miller [8].

The R-Matrix Method
The R-matrix method was developed in nuclear physics by Wigner and Eisenbud [9]. Phil Burke [10] and his colleagues in the UK, turned the R-matrix method into one of the most successful computational approaches to treat the scattering of electrons from atoms. A few years later, Schneider [11,12] generalized the method to treat electron scattering from molecules and in a significant development [13,14] showed how the theory could be used effectively to treat nuclear motion beyond the usual adiabatic approximation. The central Rmatrix idea, including the relevant equations, is shown in Figure 3. One encloses the strong interaction region in a sphere, imposes a boundary condition on the surface of the sphere to get a discrete eigenvalue problem, solves that problem for its eigenvalues and eigenfunctions and uses that information to extract the scattering matrix by matching to asymptotically known functions. In contrast to nuclei, where one may assume that the potential is zero outside the interaction region, the atomic scattering problem needs to treat the long-range multipolar interaction of the electron and target more precisely. In the atomic R-matrix method this is done using solutions to coupled long range differential equations using these multipolar potentials. One may also take the R-matrix on the inner sphere and directly propagate it to larger distances using the coupled differential equations. The external region is devoid of exchange and electron correlation and this makes for a much simplified treatment. The functions used to expand Ψ k are completely arbitrary. In practice, they are chosen such that the ϰ q vanishes at r = a. This restriction may be relaxed [13,14] to include nuclear motion. One approach which works away from threshold and away from resonances, is to take the fixed nuclei scattering amplitudes at a series of geometries, multiply by a ro-vibrational wavefunction and integrate over the nuclear coordinates. Near a resonance this does not work since the electron is trapped in a quasi-bound state and the molecular cloud wants to adjust to the change in the field but they do so adiabatically. So, the essence of the Born-Oppenheimer approximation is retained. The compound state vibrational energies determine the positions of the resonances and the overlap of the compound state and neutral vibrational eigenfunctions (Franck-Condon factors ) are the strengths. The lower left panel in Figure 4 compares the vibrational excitation cross sections from [15] and [14]. Note that [15] fails to reproduce the details of the compound levels. The lower right panel compares [14] with experiment [16,17]. The calculations agree with experiment quite well and in addition, allowed the absolute normalization of the cross sections to be determined for the first time. Finally, I must mention a new and extremely robust R-matrix approach, using a B-spline, non-orthogonal basis, that has been developed and applied to electron atom collisions by Zatsarinny [18]

The J-Matrix Method
The J-Matrix method was developed by Reinhardt, Heller, Yamani and Fishman [19,20,21]. Applications to both single channel and multichannel problems may be found in the literature. [22,2]. The basic idea is straightforward. For many simple potentials, for example, the free particle or coulomb potential, there are tri-diagonal, analytic expressions available for matrix elements of the appropriate one-particle Hamiltonian in a known L 2 basis. There are also analytic expressions for the coefficients of the regular and irregular solutions to these one particle operators which are exploited in the final solution of the scattering problem. The fundamental idea is quite similar to the R-matrix theory. The matching in space to the unknown linear combination of regular and irregular solutions to the asymptotic Hamiltonian is replaced by the analytically known expansion coefficients when the potential has assumed its asymptotic form. The coefficient joining the two solutions provides the scattering information. All the unknown coefficients are determined by projection onto the complete basis set using the known tridiagonal form of the asymptotic Hamiltonian in the basis.

Complex Scaling
Simon [23,24] extended and made the work of Balslev and Combes [25] on the spectral theory of Schrödinger operators known to atomic physicists. These authors introduced the original concept of complex scaling (let r → r exp(iθ ), an approach that was quickly exploited by a number of researchers, particularly the quantum chemistry community, as a practical tool to compute the positions and widths of scattering resonances. An unfortunate problem was encountered early on in that the complex scaling made the computations of the integrals involving the now, highly oscillatory inner orbitals, quite difficult. Another issue, raised by McCurdy and Rescigno [26] was that complete complex scaling would not work for molecules in a practical sense ( within the BO approximation ) and suggested the use of a mixed set of real and complex basis functions as a substitute. Simon [27] put the idea of using exterior complex scaling (ECS) on a firm mathematical foundation and this method remains one of most robust approaches even today. It has been advantageously used to produce the first truly accurate calculation of impact ionization of H [28] (Figures 5 and 6 )

Equivalent Quadrature
The idea of relating the diagonalization of a Hamiltonian in an L 2 basis to a numerical quadrature was put forward by Heller, Reinhardt and Yamani [29]. This relationship yields a numerical approximation to the Greens's function which properly treats the iε limit required in the scattering problem. The concept was interesting from a historical perspective but never took off as a consequence of the difficulty of having to relate the L 2 basis to a quadrature. This can be done only for a few basis sets. I mention it here for completeness.

The Kohn Variational Methods
Arguably, the Kohn [30,31,32,33,34] (KVP) and the Schwinger [35,36](SVP) variational methods, have, along with the R-matrix method, been the most productive of the methods developed for the electron scattering problem. The Schwinger variational method has the advantage that it does not require the trial function to have the correct asymptotic behavior. A practical disadvantage is that it is based on an integral equation, and requires matrix elements of the Green's function. Here I will concentrate on the Kohn

Time-Dependent (TD) Processes
In the past ten years there has been a significant effort devoted to computational studies of the interaction of intense, short pulse (attosecond) electromagnetic fields with atoms and molecules. The computational effort that is required is substantial, even for relatively small atomic and molecular targets. As a consequence, it is not surprising that the use of supercomputers and parallel processing techniques have become the norm. In the TD approach, one works with wavepackets [ having an energy spread ], needs to [ often ] propagate for long distances and times and then extracts the scattering information by projection onto asymptotic states. The plus side is that there is no need to impose boundary conditions on the wavefunction, which can be difficult to accomplish in ionizing collisions. Space limitations prevent only a brief look at this fascinating and growing field of physics. Arguably, the most effective approach to the space discretization in one and two-electron atomic and molecular systems is the FEDVR [39]. The advantage of this approach in practical calculations is that it leads to a sparse and structured Hamiltonian matrix which is amenable to efficient parallelization. The FEDVR discretization procedure is briefly shown in Equation (7). When coupled to the Short Iterative Lanczos (SIL) process as outlined in Equations (8), (9) and (10), one obtains a powerful computational approach to the TD Schrödinger equation. As can be seen in Figure 10, the method has been successfully applied to the double ionization of the He atom from the energy region where only direct double ionization is possible up to the sequential threshold. The scaling of the method on a variety of parallel processors is shown in Figure 11. One should note that after an initial erratic beginning, the curve flattens indicating that each core is now performing 1/n of the total work. As one can see, only very sophisticated numerical methods are able to reproduce the experimental results. The large rise of the cross section approaching the sequential threshold has been called virtual sequential double ionization in analogy to what is seen from the pole structure due to a resonance. There are two important points to emphasize. The first is that for many problems only very refined physical models, requiring intensive computation, can be relied upon to be predictive. The second is that it is critical to develop algorithms and software that scale well with processor count. A nice summary of what has been happening in the area of attosecond physics may be found in [41] Summary of FEDVR Discretization

i.
Within each finite element: where ϕ i (x) are members of a classical set of orthogonal functions and x i and w i are the corresponding Gauss quadrature points and weights. ii.

Connection between elements
The elements are connected by ensuring that the functions at the element boundaries are continuous.
The structure of the Hamiltonian in one dimension is,

Summary of TD Propagation
Short Time Lanczos (SIL) propagator The Three Term Lanczos Recursion

Conclusion
In this article I have tried to give the reader a flavor of many of the exciting theoretical and computational developments in atomic and molecular physics since it became possible to perform large scale calculations on computers. The early pioneers has to contend with slow processors and tiny memories by todays standards. Luckily, their processors and memories were quite sophisticated. It is a tribute to them that even with substantial digital obstacles they were able to produce many new and interesting results. Today, we are looking towards computers that will operate at the exascale level. The challenges to use that computer power effectively are formidable. If we can meet these challenges by developing the algorithms and software to use these new devices, I am sure the rewards will be significant and we will the truly be able to control atoms and molecules at their most fundamental level.    [14] vs Vibrational Close-Coupling [15] and (R) Theory [14] vs Experiment. [16,17].    Ramsauer-Townsend minimum. [38]