Computational treatment of electron and photon collisions with atoms, ions, and molecules: the legacy of Philip G Burke

The son of Henry and Frances Mary (Sprague) Burke, Philip George Burke, shown during the 2011 International Conference on Photonic, Electronic, and Atomic Collisions (ICPEAC) in Belfast in figure 1, was born in London on October 18, 1932. He lived in London during the Second World War Blitz, and memories of the bombings were still fresh in his mind when one of us (NSS) last reminisced with him in September 2018. He attended the University College of the South West of England, later to become the University of Exeter, where he read physics, graduating in 1953. Phil undertook a PhD at University College London (UCL), submitting a thesis in 1956 entitled, 'A theoretical investigation into the scattering of nucleons by light nuclei and the nature of nuclear forces'. He worked under the supervision of Sir Harry Massey FRS and Professor R A Buckingham.

Zoom In Zoom Out Reset image size

In 1956, Phil was appointed to a Research Fellowship at UCL, where he was given the opportunity to use the pilot ACE and DEUCE computers at the National Physical Laboratory (NPL) at Teddington. Here he interacted with Dr H H Robertson and Dr J H Wilkinson, jointly publishing a paper with Robertson (Burke and Robertson 1957), on low-energy elastic scattering of neutrons by deuterons.

In 1957, Phil was appointed Assistant Lecturer at the newly-formed Computer Unit of the University of London, which was headed by Professor R A Buckingham. Here he engaged in computational atomic and nuclear collision research using the newly installed Ferranti Mercury computer. This enabled Phil to develop further his interest and expertise in the emerging field of scientific computing.

In 1959, Phil married Valerie Martin, then a research student at UCL, and they went to work at the Lawrence Radiation Laboratory in Berkley, California. Here Phil had access to IBM 704, 709, and 7090 machines, which were far superior to any machines in the UK, and he was able to develop code in FORTRAN rather than machine code. Phil was initially employed in the preparation and analysis of experiments on the 72 inch hydrogen bubble chamber, under the general direction of Professor Louis Alvarez. After that he joined the Theory Division, where he carried out research in computational atomic and particle physics.

Phil's first paper in atomic physics was in collaboration with Ken Smith at Argonne National Laboratory (ANL). In Smith and Burke (1961) they examined the effect of virtual excitations on the elastic scattering of electrons and positrons by atomic hydrogen. The numbers reported were obtained through the use of an IBM 704 at ANL. This paper marked the start of Phil's enormous contribution to the study of resonances in atomic and molecular scattering. The following comment in the paper also testifies to Phil's early recognition of the need to develop general rather than one-off special-purpose codes: 'A generalized version of the program used to calculate the numbers reported here has been written and is being tested. With this program we hope to calculate...the total elastic cross section in the strong-coupling exchange approximation over the entire energy region and so be able to compare the results with the experimental results of Bederson et al and Brackman et al'.

Phil's second atomics physics paper was with Harry Schey (Burke and Schey 1962). It focused on the topic of low-energy elastic scattering of electrons by atomic hydrogen, a problem introduced to Phil by Dr Mike Seaton at UCL, and involved the use of an IBM 709. The comment in the paper that, 'this work was undertaken to bring to bear upon this problem calculational opportunities placed at the disposal of theorists by modern computing facilities, in the hope that certain discrepancies between theory and experiment may be resolved by the nearly exact solutions now possible', established a rationale that pervaded all of Phil's future work.

In the UK, toward the end of 1956, joint Atomic Energy Authority/Civil Service Commission Boards (Hoff Boards) were established to stem the 'brain drain' and recruit emigré scientists back to the UK. Phil was recruited in 1962 and took up an appointment as Research Fellow, later Senior Principal Scientific Officer, in the Theory Division at the Atomic Energy Research Establishment at Harwell. Phil continued his research in particle physics in collaboration with scientists at the Rutherford Appleton Laboratory and in atomic physics in support of the plasma physics programme at the Culham Laboratory. In this work he used the recently acquired IBM STRETCH computer at AWRE, Aldermaston. This productive time included a series of four papers on low-energy electron scattering by atomic hydrogen, published back-to-back in 1967 (Burke et al 1967a, Taylor and Burke 1967, Burke et al 1967b, Macek and Burke 1967).

Following the departure of Professor Alex Dalgarno FRS from Queen's University Belfast to the Harvard College Observatory and Smithsonian Astrophysical Observatory, Phil was appointed to the Chair of Mathematical Physics at Queen's University Belfast on the 1st November 1967. He delivered his inaugural lecture on 11th December 1968 on the topic entitled 'Particles and resonances in modern physics'. At the start of the lecture he commented that this topic for discussion was chosen because it provided '...some sort of unifying link across several branches of physics...it is also one which I believe enables a common language and a common appreciation to be developed for branches of physics as diverse as the properties of liquids and solids at almost zero temperatures, the behavior of atomic nuclei in the interior of stars, the behavior of gases at high temperatures, and the understanding of the interaction of matter under the acceleration of billions of volts' (Burke 1968).

Phil joined the Department of Applied Mathematics and Theoretical Physics, which was headed by Professor (later Sir) David Bates FRS. Interestingly, both had had their early careers shaped by Professor Sir Harry Massey FRS. Alan Hibbert, who joined Queen's a few weeks after Phil, recounts in Scott et al (2020) that he 'quickly established a small research group...with the aim of constructing robust and general programs to enable both atomic structure and atomic collision processes to be calculated'.

From 1977 to 1982 Phil held a joint appointment with Daresbury Laboratory, in Cheshire, where he was Head of the Theory and Computational Science Division. With typical devotion to duty, as recalled by Alan Hibbert in Scott et al (2020), 'he maintained his teaching and research commitments in Queen's by traveling, every week in term time, on the Liverpool–Belfast boat—on Wednesday night from Liverpool and returning on Friday night. He spent the days in Belfast undertaking his undergraduate teaching and meeting graduate students and colleagues'. One of us (NSS) can attest to this, as he was one of Phil's PhD students at the time.

During his 40+ years at Queen's, Phil developed further the resonance ideas highlighted in his inaugural lecture. Encouraged by Ugo Fano and Mike Seaton, he exploited his knowledge of resonances from nuclear physics to develop a generalized R-matrix theory of atomic collisions and its application to the ab initio study of atomic, molecular, and optical collision processes. Phil and his collaborators developed associated computer codes that were acknowledged as the best computational tools to model electron and photon collisions with atoms, ions, molecules and solids. The authors of this paper are privileged to have collaborated with Phil as students, colleagues, and friends over many years. As his early career was shaped by Harry Massey, so ours was shaped by Phil.

Phil continued to research and collaborate after retiring from the University on 30th September 1998. However, release from formal university duties allowed Phil to set about writing his magnum opus, a monograph on the R-matrix theory of atomic collisions and its application to atomic, molecular, and optical physics (Burke 2011). Phil's well-thumbed notebook, containing a handwritten record of every paper he had read, was put to good use. A subset of these topics and Phil's seminal contribution to their understanding are explored in the present paper.

Phil's contribution was not limited to atomic, molecular, and optical collision processes. He also had a significant influence on computational science nationally and internationally. First, Phil was instrumental in the establishment of the UK Collaborative Computational Projects (CCPs). CCPs bring together leading UK expertise in key fields of computational research to tackle large-scale scientific software development, maintenance, and distribution. Second, Phil was the founding editor of the international journal Computer Physics Communications and Director of its associated Program Library. The journal recently celebrated its 50th anniversary and is a major resource for the computational physics community worldwide. Phil's leadership in these two areas is elucidated in this paper as well.

Phil also influenced the development of supercomputing within the UK by serving on numerous government committees. He was a member of: the UK Science and Engineering Research Council (SERC) Swindon, 1989–1994; the Joint Policy Committee on Advanced Research Computing, 1988–1990; the Advisory Board for the Research Councils Supercomputing Sub-Committee, 1991–1994; and he chaired the SERC Supercomputing Management Committee, 1991–1994.

Phil was a physically tall man, standing head and shoulders above many of his colleagues. This was equally true intellectually, and Phil's eminence in the field is evidenced through the following elections, prizes, and awards. He was a Fellow of the Institute of Physics, the American Physical Society, the Royal Astronomical Society, and the European Physical Society. He was elected to membership of the Royal Irish Academy in 1974 and to Fellowship of the Royal Society of London in 1978; he was awarded the Institute of Physics' Guthrie Medal and Prize in 1994 and the David Bates Prize in 2000; he was recipient of the American Physical Society Will Allis Prize in 2012. Phil was appointed as a CBE in the 1993 Queen's Birthday Honors for services to science.

In the academic work of cut-and-thrust, Phil was a true gentleman, always looking for the best in everyone. As Alan Hibbert recounts in Scott et al (2020), 'Phil was a humble man.... He had a fine mind, and a warm personality. He was a man of integrity.' And, whether helping the children with their homework, or dealing with a complex mathematical problem from a student or colleague, he developed a solution from first principles with patience and understanding. Phil's altruism was also noted by Mike Seaton who commented in Seaton (1994), 'Phil has always given much of his time to the wider interests of the scientific community...If the UK has any world standing in computational science, then Phil Burke is the first person we have to thank'.

This manuscript is dedicated to the memory of Phil and highlights his important contributions in key areas of theoretical and computational physics. It is organized as follows. Section 2 focusses on electron collisions with atoms and ions, as well as weak-field photoionisation with continuous radiation. The final state of the latter process is effectively scattering from the residual ion. This is followed by the extension to molecular targets in section 3. We then move to photon impact with stronger (but still many-cycle) radiation. This enables multi-photon processes, which we discuss in section 4, before treating explicitly time-dependent processes such as short-pulse intense laser-atom and laser-molecule interactions in section 5. Finally, section 6 describes Phil's role in founding the Computer Physics Communications journal and its associated Program Library, as well as CCPs in the United Kingdom.

There are many ways to formulate the electron–atom/ion collision problem, as well as the closely related case of (weak-field) photoionisation by essentially continuous radiation. In the latter case, one needs transition elements from the initial bound state to a continuum state that corresponds to half an electron collision with the residual ion or the neutral atom in photodetachment of negative ions. Due to selection rules (we will limit ourselves to the electric dipole approximation), only a few partial waves will contribute to photoionisation, while the partial-wave expansion for pure electron collisions contains many (in principle, an infinite number of) partial waves. In the latter case, the convergence of the expansion needs to be checked carefully.

Here we will concentrate on just one method, namely the time-independent close-coupling expansion. Phil Burke pioneered many numerical methods for the solution of the resulting equations. We will only list very few here, but point to a number of references that interested readers should consult for details. These include the review with his long-term colleague M J Seaton (Burke and Seaton 1971), a description of numerical methods for the Opacity Project (Berrington et al 1987), his extensive monograph about the use of the R-matrix method for a large number of processes in atomic, molecular, and optical (AMO) physics (Burke 2011), and his chapter in G W F Drake's AMO Handbook (Burke 2006), which Phil entrusted two of us (KB and JT) to update.

To describe the general ideas, we use excerpts of Phil's own words from the Handbook chapter, in which he first describes the general close-coupling expansion, formulated in coordinate space, for electron–atom and electron–ion scattering before concentrating on a particular way of solving the resulting equations using the R-matrix method (Burke et al 1971, Burke and Robb 1975). We then select a few examples to illustrate the developments that took place over the past 50 years (about 1965–2015), when Phil was extremely active.

After introducing the R-matrix method for low-energy electron–atom and electron–ion collisions, Phil initiated a number of major extensions. These include photoionisation (Burke and Taylor 1975), relativistic effects within the Breit–Pauli approximation (Scott and Burke 1980), a hybrid distorted-wave + R-matrix treatment of electron-impact ionisation (Bartschat and Burke 1987), the intermediate-energy R-matrix method (IERM) (Burke et al 1987), and the R-matrix with pseudo-states (RMPS) (Bartschat et al 1996) approach. Even more (molecules, multi-photon and strong-field, time-dependent processes) will be described in other sections of this manuscript.

We begin with Phil's (abbreviated) description of the close-coupling expansion. If all relativistic effects are neglected, the time-independent Schrödinger equation (TISE) describing the scattering of an electron by a target atom or ion containing N electrons and nuclear charge Z is

$$\begin{equation}{H}_{N+1}{\Psi}=E\enspace {\Psi},\end{equation} \tag{ 1 }$$

where E is the total energy of the system. The (N + 1)-electron nonrelativistic Hamiltonian H_N+1 is given in atomic units by

$$\begin{equation}{H}_{N+1}=\sum _{i=1}^{N + 1}\left(-\frac{1}{2}{\nabla }_{i}^{2}-\frac{Z}{{r}_{i}}\right)+\sum _{i{ >}j=1}^{N + 1}\frac{1}{{r}_{ij}},\end{equation} \tag{ 2 }$$

where r_ij = |r_i − r_j| with r_i and r_j denoting the vector coordinates of electrons i and j relative to the origin of coordinates taken to be the target nucleus, which is assumed to have infinite mass.

In order to solve the Schrödinger equation to obtain the scattering amplitude and cross section at low energies, we make a partial-wave expansion of the total wave function

$$\begin{align}\hfill {{\Psi}}_{j}^{{\Gamma}}\left({\boldsymbol{X}}_{N + 1}\right)& =\mathcal{A}\sum _{i=1}^{n}{\overline{{\Phi}}}_{i}^{{\Gamma}}\left({\boldsymbol{x}}_{1},\dots ,{\boldsymbol{x}}_{N};{\hat{\boldsymbol{r}}}_{N + 1}{\sigma }_{N + 1}\right)\frac{{F}_{ij}^{{\Gamma}}\left({r}_{N + 1}\right)}{{r}_{N + 1}}\hfill \\ \hfill & \quad +\sum _{i=1}^{m}{\chi }_{i}^{{\Gamma}}\left({\boldsymbol{x}}_{1},\dots ,{\boldsymbol{x}}_{N + 1}\right){b}_{ij}^{{\Gamma}},\hfill \end{align} \tag{ 3 }$$

where X_N+1 ≡ x₁, x₂...x_N+1 represents the space and spin coordinates of all N + 1 electrons, x_i ≡ r_iσ_i represents the space and spin coordinates of the ith electron, and $\mathcal{A}$ is the operator that antisymmetrizes the first summation with respect to exchange of all pairs of electrons in accordance with the Pauli exclusion principle. The channel functions ${\overline{{\Phi}}}_{i}^{{\Gamma}}$, assumed to be n in number, are obtained by coupling the orbital and spin angular momenta of the target states Φ_i with those of the scattered electron to form eigenstates of the total orbital and spin angular momenta L and S, their z components M_L and M_S, and the parity π, where Γ ≡ LM_LSM_Sπ is conserved in the collision. The square-integrable correlation functions ${\chi }_{i}^{{\Gamma}}$ allow for additional correlation effects not included in the first expansion in equation (3), which runs over a limited number of target eigenstates and possibly pseudo-states.

By substituting equation (3) into the Schrödinger equation, projecting onto the channel functions ${\overline{{\Phi}}}_{i}^{{\Gamma}}$ and the square-integrable functions ${\chi }_{i}^{{\Gamma}}$, and eliminating the coefficients ${b}_{ij}^{{\Gamma}}$, we obtain n coupled integro-differential equations satisfied by the reduced radial functions ${F}_{ij}^{{\Gamma}}$ representing the motion of the scattered electron:

$$\begin{align}\hfill & \left(\frac{{\mathrm{d}}^{2}}{\mathrm{d}{r}^{2}}-\frac{{\ell }_{i}\left({\ell }_{i}+1\right)}{{r}^{2}}+\frac{2\left(Z-N\right)}{r}+{k}_{i}^{2}\right){F}_{ij}^{{\Gamma}}\left(r\right)\hfill \\ \hfill & =2\sum _{\ell }\left\{{V}_{i\ell }^{{\Gamma}}\left(r\right){F}_{\ell j}^{{\Gamma}}\left(r\right)\right.\hfill \\ \hfill & \quad +\left.{\int }_{0}^{\infty }\left[{K}_{i\ell }^{{\Gamma}}\left(r,{r}^{\prime }\right)+{X}_{i\ell }^{{\Gamma}}\left(r,{r}^{\prime }\right)\right]{F}_{\ell j}^{{\Gamma}}\left({r}^{\prime }\right)\mathrm{d}{r}^{\prime }\right\},\hfill \end{align} \tag{ 4 }$$

here ℓ_i is the orbital angular momentum of the scattered electron, while ${V}_{i\ell }^{{\Gamma}}$, ${W}_{i\ell }^{{\Gamma}}$, and ${X}_{i\ell }^{{\Gamma}}$ are the local direct, nonlocal exchange, and nonlocal correlation potentials, respectively. If the correlation potential, which arises from the ${\chi }_{i}^{{\Gamma}}$ terms in equation (3), is not included, equation (4) are called the close-coupling equations. These need to be solved subject to the appropriate boundary conditions for collision processes (Burke 2006).

The direct potential can be written as

$$\begin{align}\hfill {V}_{ij}^{{\Gamma}}\left({r}_{N + 1}\right)& =\left\langle {\overline{{\Phi}}}_{i}^{{\Gamma}}\left({\boldsymbol{x}}_{1},\dots ,{\boldsymbol{x}}_{N};{\hat{\boldsymbol{r}}}_{N + 1}{\sigma }_{N + 1}\right)\left\vert \sum _{i=1}^{N}\frac{1}{{r}_{iN + 1}}-\frac{N}{{r}_{N + 1}}\right\vert \right.\hfill \\ \hfill & \quad {\times}\left.{\overline{{\Phi}}}_{j}^{{\Gamma}}\left({\boldsymbol{x}}_{1},\dots ,{\boldsymbol{x}}_{N};{\hat{\boldsymbol{r}}}_{N + 1}{\sigma }_{N + 1}\right), \right\rangle ,\hfill \end{align} \tag{ 5 }$$

where the integral is taken over all electron space and spin coordinates, except for the radial coordinate of the (N+1)-th electron. This potential has the asymptotic form

$$\begin{equation}{V}_{ij}^{{\Gamma}}\left(r\right)=\sum _{\lambda =1}^{\lambda \mathrm{max}}{a}_{ij}^{\lambda }{r}^{-\lambda -1}\quad \mathrm{f}\mathrm{o}\mathrm{r}\;\;r{\geqslant}a,\end{equation} \tag{ 6 }$$

where a is the range beyond which the orbitals in the target states Φ_i included in the first expansion in equation (3) are negligible.

The exchange and correlation potentials, unlike the direct potential, are both nonlocal, and the exchange potential vanishes exponentially for large r. In practice, all these potentials are evaluated by general computer programs.

The scattering amplitude and cross section can be obtained by solving equation (4) for all relevant conserved quantum numbers Γ subject to the K-matrix asymptotic boundary conditions

$$\begin{equation}{F}_{ij}^{{\Gamma}}\underset{r\to \infty }{\approx }{k}_{i}^{-\frac{1}{2}}\left(\mathrm{sin}\enspace {\theta }_{i}{\delta }_{ij}+\mathrm{cos}\enspace {\theta }_{i}{K}_{ij}^{{\Gamma}}\right)\mathrm{f}\mathrm{o}\mathrm{r}\enspace \mathrm{o}\mathrm{p}\mathrm{e}\mathrm{n}\enspace \mathrm{c}\mathrm{h}\mathrm{a}\mathrm{n}\mathrm{n}\mathrm{e}\mathrm{l}\mathrm{s}\enspace \left({k}_{i}^{2}{\geqslant}0\right);\end{equation} \tag{ 7a }$$

$$\begin{equation}\quad \quad \quad \quad {F}_{ij}^{{\Gamma}}\underset{r\to \infty }{\approx }0\enspace \mathrm{f}\mathrm{o}\mathrm{r}\enspace \mathrm{c}\mathrm{l}\mathrm{o}\mathrm{s}\mathrm{e}\mathrm{d}\enspace \mathrm{c}\mathrm{h}\mathrm{a}\mathrm{n}\mathrm{n}\mathrm{e}\mathrm{l}\mathrm{s}\enspace \left({k}_{i}^{2}{< }0\right),\end{equation} \tag{ 7b }$$

here

$$\begin{equation}{\theta }_{i}={k}_{i}r-\frac{1}{2}{\ell }_{i}\pi +\frac{z}{{k}_{i}}\enspace \mathrm{ln}\left(2{k}_{i}r\right)+{\sigma }_{{\ell }_{i}}\end{equation} \tag{ 8 }$$

with z = Z − N and ${\sigma }_{{\ell }_{i}}$ is the Coulomb phase. The S-matrix and T-matrix are related to the K-matrix defined by equation (7a) through the matrix equations

$$\begin{equation}{\boldsymbol{S}}^{{\Gamma}}=\frac{\boldsymbol{I}+\mathrm{i}{\boldsymbol{K}}^{{\Gamma}}}{\boldsymbol{I}-\mathrm{i}{\boldsymbol{K}}^{{\Gamma}}},\quad {\boldsymbol{T}}^{{\Gamma}}={\boldsymbol{S}}^{{\Gamma}}-\boldsymbol{I}=\frac{2\mathrm{i}{\boldsymbol{K}}^{{\Gamma}}}{\boldsymbol{I}-\mathrm{i}{\boldsymbol{K}}^{{\Gamma}}}.\end{equation} \tag{ 9 }$$

The dimensions of the matrices in these equations are n_a × n_a, where n_a is the number of open channels at the energy under consideration for the given Γ. The hermiticity and time-reversal invariance of the Hamiltonian ensures that K^Γ is real and symmetric, while S^Γ is unitary and symmetric.

Early calculations using the close-coupling equation by solving the resulting integro-differential equations were carried out by Phil and a number of co-workers, one of them being his wife Val, for simple targets in few-state approximations. These included the inclusion of the 1s, 2s, and 2p states in atomic hydrogen (Burke et al 1962) or He⁺, including calculations of the photoionisation of helium (Burke and McVicar 1965), which reproduced the well-known series of Fano resonances (Fano 1961) that had also been seen experimentally (Madden and Codling 1965). Already in 1963, however, Phil noted problems with the convergence of the close-coupling expansion if only physical discrete target states were included. Adding the physical n = 3 states to the expansion dropped the theoretical predictions for the excitation of the 2s and 2p states dramatically, leading Phil to the conclusion that '...the experimental curves could well be correct and the outstanding discrepancy due to the slow convergence of the close-coupling expansion...' (Burke 1963).

In 1970, he decided to use so-called 'pseudo-states', which are built from square-integrable functions that are, however, not eigenstates of the Hamiltonian. They have the correct number of nodes to maintain orthogonality, but their range is shorter than that of the physical orbitals. In fact, one can choose the range to optimize their effect on the predictions for the physical states of interest. As a result, the pseudo-orbitals have higher energies than the physical bound orbitals, including positive energies that approximate coupling to the ionisation continuum.

Figure 2 shows results of Burke and Webb (1970) for the excitation cross section of the 1s → 2p transition atomic hydrogen. With just two pseudo-states corresponding to $3\overline{\mathrm{s}}$ and $3\overline{\mathrm{p}}$, the predicted cross section drops dramatically, and much better agreement with experiment is achieved. The authors conclude '...with new techniques of solving coupled integro-differential equations (Burke et al 1971) there is no reason why more pseudo-states should not be included and the convergence of the pseudo-state expansion studied in greater detail' (Burke and Webb 1970). Interestingly, it took another 25 years until this was finally done in a systematic way in Belfast to produce the RMPS (Bartschat et al 1996) approach, a few years after our Australian colleagues developed the 'convergent close-coupling' (CCC) method (Bray and Stelbovics 1992) using the very same basic idea but a completely different numerical implementation.

Zoom In Zoom Out Reset image size

The reference to 'Burke et al (1971)' points to a key paper of Phil's career. It describes the introduction of the R-matrix method for atomic collisions (see also Burke and Robb (1975)). It soon became, and still remains, the method of choice for low- and intermediate-energy electron and photon collisions with complex atoms, ions, and molecules.

Once again, we use an excerpt of Phil's own writing for the AMO Handbook. The R-matrix method was first introduced in nuclear physics (Wigner 1946, Wigner and Eisenbud 1947) in a study of resonance reactions. It has since been applied to a wide range of atomic, molecular, and optical processes. A comprehensive overview of the method and its many applications can be found in Burke (2011).

The method starts by partitioning configuration space into two regions by a sphere of radius a, chosen so that the direct potential has achieved its asymptotic form given by (6) and the exchange and correlation potentials are negligible for r ⩾ a. The objective is then to calculate the R-matrix elements ${R}_{ij}^{{\Gamma}}\left(E\right)$, which are defined by

$$\begin{equation}{F}_{ij}^{{\Gamma}}\left(a\right)=\sum _{\ell =1}^{n}{R}_{i\ell }^{{\Gamma}}\left(E\right){\left(a\frac{\mathrm{d}{F}_{\ell j}^{{\Gamma}}}{\mathrm{d}r}-{b}_{\ell }{F}_{\ell j}^{{\Gamma}}\right)}_{r=a},\end{equation} \tag{ 10 }$$

by solving (4) in the internal region.

The collision problem is solved for r ⩽ a by expanding the wave function, in analogy with (3), in the form

$$\begin{align}\hfill {{\Psi}}_{k}^{{\Gamma}}\left({\boldsymbol{X}}_{N + 1}\right)& =\mathcal{A}\sum _{ij}{\overline{{\Phi}}}_{i}^{{\Gamma}}\left({\boldsymbol{x}}_{1},\dots ,{\boldsymbol{x}}_{N};{\hat{\boldsymbol{r}}}_{N + 1}{\sigma }_{N + 1}\right)\hfill \\ \hfill & \quad {\times}{r}_{N + 1}^{-1}{u}_{j}\left({r}_{N + 1}\right){a}_{ijk}^{{\Gamma}}\hfill \\ \hfill & \quad +\sum _{i}{\chi }_{i}^{{\Gamma}}\left({\boldsymbol{x}}_{1},\dots ,{\boldsymbol{x}}_{N + 1}\right){b}_{ik}^{{\Gamma}},\hfill \end{align} \tag{ 11 }$$

where the u_j are radial basis functions defined over the range 0 ⩽ r ⩽ a. For radial basis functions u_j satisfying arbitrary boundary conditions at r = a, the Hamiltonian H_N+1 defined by equation (2) is not hermitian in the internal region due to the kinetic energy operators. It can, however, be made hermitian by adding the Bloch operator (Bloch 1957)

$$\begin{equation}{L}_{\text{b}}=\sum _{i=1}^{n}\vert {\overline{{\Phi}}}_{i}^{{\Gamma}}\rangle \frac{1}{2}\delta \left(r-a\right)\left(\frac{\mathrm{d}}{\mathrm{d}r}-\frac{b-1}{r}\right)\langle {\overline{{\Phi}}}_{i}^{{\Gamma}}\vert \end{equation} \tag{ 12 }$$

to H_N+1, where b is an arbitrary parameter. The Schrödinger equation (1) then becomes

$$\begin{equation}\left({H}_{N+1}+{L}_{\text{b}}-E\right){{\Psi}}^{{\Gamma}}={L}_{\text{b}}{{\Psi}}^{{\Gamma}},\end{equation} \tag{ 13 }$$

which can be formally solved as

$$\begin{equation}{{\Psi}}^{{\Gamma}}={\left({H}_{N+1}+{L}_{\text{b}}-E\right)}^{-1}{L}_{\text{b}}{{\Psi}}^{{\Gamma}}.\end{equation} \tag{ 14 }$$

We now expand the Green's function ${\left({H}_{N+1}+{L}_{\text{b}}-E\right)}^{-1}$ in terms of the basis ${{\Psi}}_{k}^{{\Gamma}}$, where the coefficients ${a}_{ijk}^{{\Gamma}}$ and ${b}_{ik}^{{\Gamma}}$ in (11) are chosen to diagonalize H_N+1 + L_b according to

$$\begin{equation}\left\langle {{\Psi}}_{k}^{{\Gamma}}\left\vert {H}_{N+1}+{L}_{\text{b}}\right\vert {{\Psi}}_{{k}^{\prime }}^{{\Gamma}}\right\rangle ={E}_{k}^{{\Gamma}}{\delta }_{k{k}^{\prime }}.\end{equation} \tag{ 15 }$$

Equation (14) can then be written as

$$\begin{equation}\left\vert {{\Psi}}^{{\Gamma}}\right\rangle =\sum _{k}\frac{\left\vert {{\Psi}}_{k}^{{\Gamma}}\right\rangle \left\langle {{\Psi}}_{k}^{{\Gamma}}\right\vert }{{E}_{k}^{{\Gamma}}-E}{L}_{\text{b}}\left\vert {{\Psi}}^{{\Gamma}}\right\rangle .\end{equation} \tag{ 16 }$$

Finally, we project this equation onto the channel functions ${\overline{{\Phi}}}_{i}^{{\Gamma}}$ and evaluate it at r = a. Assuming Ψ^Γ is given by (3), we retrieve (10), where the R-matrix elements are calculated from the expansion

$$\begin{equation}{R}_{ij}^{{\Gamma}}\left(E\right)=\frac{1}{2a}\sum _{k}\frac{{w}_{ik}^{{\Gamma}}{w}_{jk}^{{\Gamma}}}{{E}_{k}^{{\Gamma}}-E},\end{equation} \tag{ 17 }$$

here we have introduced the surface amplitudes

$$\begin{equation}{w}_{ik}^{{\Gamma}}=\sum _{j}{u}_{j}\left(a\right){a}_{ijk}^{{\Gamma}}.\end{equation} \tag{ 18 }$$

The main part of the calculation involves setting up and diagonalizing the matrix given by equation (15). This has to be carried out only once to determine the R-matrix for all energies E, thereby making the method extremely efficient if results for many energies are required, e.g., for a detailed tracing of resonance structures.

In the external region r ⩾ a, equation (4) reduce to ordinary differential equations coupled by the potential ${V}_{ij}^{{\Gamma}}\left(r\right)$, which has achieved its asymptotic form (6). There are many numerical implementations of the R-matrix method in atomic, molecular, and optical physics. Some of the packages will be mentioned in section 6.

Another very important step in the development of the R-matrix method was the extension to photoionisation processes. In this case, the essential problem is to calculate the dipole matrix element

$$\begin{equation}\langle {{\Psi}}_{f}^{-}\left(\hat{\boldsymbol{k}}\right)\vert \sum _{j=1}^{N+1}\hat{\boldsymbol{\varepsilon }}\cdot {\boldsymbol{r}}_{j}\enspace \vert {{\Psi}}_{i}\rangle \end{equation} \tag{ 19 }$$

in the length form of the electric dipole operator or the corresponding matrix element if the velocity form is used. In (19), ${{\Psi}}_{f}^{-}\left(\hat{\boldsymbol{k}}\right)$ describes the photoelectron emitted along the direction $\hat{\boldsymbol{k}}$, i.e., effectively a scattering state from the residual N-electron ion with appropriate boundary conditions indicated by the superscript −, $\hat{\boldsymbol{\varepsilon }}$ is the unit vector describing the polarization of the electric field, and Ψ_i is the initial (N + 1)-electron bound state.

The general theory was outlined by Burke and Taylor (1975) and immediately applied to Ne, Ar, and Al. Figure 3 shows the results obtained by LeDourneuf et al (1975) for photoionisation of $\mathrm{A}\mathrm{l}{\left({\mathrm{3}\mathrm{s}}^{2}3\mathrm{p}\right)}^{2}{\mathrm{P}}^{\text{o}}$ leading to any of the final states ${\mathrm{A}\mathrm{l}}^{+}{\left(3{\mathrm{s}}^{2}\right)}^{1}{\mathrm{S}}^{\text{e}}$, Al⁺(3s3p)³P^o, and Al⁺(3s3p)¹P^o, respectively, as a function of the wavelength. Note the series of Rydberg resonances that converge to the Al⁺(3s3p)³P^o and Al⁺(3s3p)¹P^o thresholds when these final ionic states are still energetically closed.

Zoom In Zoom Out Reset image size

Burke and Taylor (1975) end their abstract with '...showing that the method has a wide range of applicability'. Looking of at the developments over the past 45 years, even this seemingly optimistic outlook was an understatement. The R-matrix approach for photoionisation as outlined in the above paper has been, and continues to be, the method of choice for numerous applications, especially those relevant for the interpretation of astrophysical observations.

The ground work for the next major development of the R-matrix method for electron collisions with atoms and ions was laid by Scott and Burke (1980) who outlined the extension to include terms of the Breit–Pauli Hamiltonian. This enabled more reliable calculations for targets with intermediate and even high nuclear charge Z, such as Fe (Z = 26) and its many ions that are also highly relevant for astrophysics. In particular, transitions between fine-structure levels could now be described ab initio by including some dynamical effects, rather than being limited to algebraic recoupling of nonrelativistic results. Burke and Scott end their conclusions with the statement: 'we also intend to apply the package to study electron scattering by heavy neutral atoms where there is beginning to be reliable low-energy experimental data for comparison.'

In retrospect, that statement described perfectly my situation in 1982. A year before I had obtained my Diploma (Master) degree at the University of Münster in Germany, participating in experiments on spin-polarization effects in e-Hg scattering (Bartschat et al 1981, Bartschat et al 1982) in the group of Joachim Kessler and Friedrich Hanne, but simultaneously carrying out essentially algebraic calculations under the guidance of Karl Blum. I learned a lot in both of those settings, but I noticed that all dynamical parameters were just given names, such as 'cross section', 'asymmetry function', 'Sherman function', 'state multipoles', etc. When I asked Karl if we could ever calculate these things explicitly (beyond saying that they might be zero in certain situations (Bartschat and Blum 1982)), he said: 'you need to go to Belfast and work with Phil Burke.' I was fortunate to obtain a stipend from the British Council to do just that and arrived in Belfast in October of 1982. Phil was a very busy man, and so Stan Scott was assigned to take care of me. Looking at the history of Northern Ireland, this was certainly not the easiest time to be in Belfast. But I still remember Stan saying to me at our first coffee break: 'there is nothing to worry about here when it comes to physics. Phil knows everything!' I found out over the next 30 years that Stan was absolutely correct.

We then set out to perform a five-state Breit–Pauli calculation on e-Hg collisions, with special emphasis on excitation of the (6s6p)³P₁ and (6s6p)³P₂ states, which I had worked on the measurements at Münster. The results for the cross sections for elastic scattering and excitation of the (6s6p)³P_0,1,2 states are shown in figure 4 (Scott et al 1983). Note, in particular, the very sharp near-threshold resonance in the excitation of the (6s6p)³P₁ state. We (and some others) had seen this resonance before, and Heddle (1975) actually suspected it to be of predominantly ${\left({\mathrm{6}\mathrm{s}\mathrm{6}\mathrm{p}}^{2}\right)}^{4}{\mathrm{P}}_{5/2}$ character. Looking at the various partial-wave contributions, we could confirm this by our ab initio calculation, which stood as the benchmark for this collision system until the mid 1990s and was used to interpret a number of follow-up experiments. This particular resonance (strength and position) was actually critical for the success of the Franck–Hertz experiment (Franck and Hertz 1914), as described in detail by Hanne (1988).

Zoom In Zoom Out Reset image size

Recall that weak-field photoionisation is effective electron scattering from the residual ion, and hence there will be numerous resonances that appear as Rydberg series converging to the thresholds of closed excited ionic states. A similar feature occurs in another development that carries Phil's signature. When I was in Belfast on a post-doctoral stipend from the German Research Council in 1985, he suggested developing a hybrid method for electron impact ionisation, in which a 'fast' projectile was described by a distorted wave while the initial state and the 'slow' ejected-electron–residual-ion interaction would be treated by an R-matrix expansion. This required replacing the single energy-independent dipole operator in equation (19) by an energy-dependent operator ${\sum }_{j=1}^{N+1}{R}_{{\ell }_{0}{\ell }_{1}}^{\lambda }\left({E}_{0},{E}_{1},{r}_{j}\right)$, which depends on the partial-wave angular momenta (ℓ₀, ℓ₁), the initial and final energies (E₀, E₁) of the projectile, and the multipolarity λ after the expansion of the Coulomb interaction between the projectile and the N + 1 target electrons (Bartschat and Burke 1987).

I still remember sitting in a train from my home town in Germany to Amsterdam (from which I would fly back to Belfast) in January of 1985, working out the number of integrals one would need for this hybrid distorted-wave + R-matrix method for electron-impact ionisation. Back in Belfast, I told Phil that it was pretty hopeless, and we should just drop the idea altogether. His answer was: 'do not worry. By the time you get your program to work, we will have a computer to run it on.' He was (of course) correct, and we published the first paper with actual results in (Bartschat and Burke 1987). It was one of several cases in which Phil seemed to have perfect timing in anticipating how computers would develop and how one would be able to put his ideas into practice by using the latest technology.

Figure 5 shows the single-differential (with respect to the energy loss) cross section for electron impact ionisation of neutral argon atoms for an incident electron energy of 136 eV (Bartschat and Burke 1988). For the 'dipole' component (λ = 1), we see the same type of window resonances that would appear in photoionisation, but the other partial waves of the ejected electron also show resonances, particularly for the monopole (λ = 0) and quadrupole (λ = 2) contributions. As Phil predicted, it took me quite a while to develop the program, but it was ultimately used (and is still being used in appropriate kinematical situations (deHarak et al 2008, Martin et al 2018)) with significant success. Phil encouraged me to write it up for Computer Physics Communications (Bartschat 1993), and we also decided to convert Phil's original many-channel resonance fitting program into a publicly available version (Bartschat and Burke 1986). Once again I learned a lot from him in the process.

Zoom In Zoom Out Reset image size

Clearly, reliable calculations for electron impact ionisation became one of Phil's priorities at the same time. In addition to the hybrid method, he suggested the IERM approach (Burke et al 1987). The main idea was to modify equation (11) to have two continuum functions in the inner region of the R-matrix sphere. While this allows, in principle, to handle ionisation by then allowing two electrons also in the continuum, it leads to difficulties in the matching procedure.

Figure 6 shows results for one partial-wave T-matrix element in elastic scattering from atomic hydrogen (Scott et al 1989). The very first resonance, ${\left(2{\mathrm{p}}^{2}\right)}^{1}{\mathrm{D}}^{\text{e}}$, is indeed a physical resonance, while all the other structures are pseudo-resonances. Phil and his collaborators then developed an averaging procedure over those unphysical resonances and obtained excellent results for the cross sections. An example is shown in figure 7 for excitation of the 1s–2s transition in atomic hydrogen. In this case, coupling to the continuum is very important, as can be seen by the failure of a simple three-state (1s, 2s, 2p) close-coupling model (Kingston et al 1976).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The last item in my list of Phil's contributions is the development of the RMPS approach. As mentioned above, he considered the idea already in 1970, but in those days the computational facilities were simply insufficient to push things to convergence. I was fortunate to spend part of a sabbatical in Belfast during the fall of 1995. At that time, we were impressed by the success of the CCC method developed by Bray and Stelbovics (1992) in Australia a few years earlier. They used a Laguerre basis to obtain a large number of pseudo-states and were able to convincingly demonstrate the convergence of the close-coupling expansion—provided a sufficient number of positive-energy pseudo-states were included in the close-coupling expansion—by solving the Lippman–Schwinger equations for the T-matrix in momentum space.

Due to the structure of the R-matrix codes, it seemed straightforward to do the same using a Sturmian basis, which should yield the same final results as long as one includes a sufficient number of functions. Without going into details, once again some numerical issues had to be overcome. On some week-end in October or November of 1995, I had identified one of those problems, immediately contacted Val Burke to let Phil know (Phil did not do email in those days!), and even came up with what I thought was a complicated but viable solution. When we met on the following Monday morning in Phil's office, he suggested 'you go first'. I outlined my procedure, and then he said something like this: 'very good, but I also have an idea. How about this?' Not surprisingly, we ended up using Phil's suggestion instead of mine. The theory was outlined in Bartschat et al (1996).

Figure 8 shows a benchmark calculation for e-H collisions, in which we compared IERM, CCC, and RMPS results for elastic scattering and excitation of the 2s and 2p states of atomic hydrogen between the n = 2 and n = 3 thresholds. The curves are almost indistinguishable, except when the numbers become very small (e.g., for L = 3) and are more sensitive to small numerical issues. Today the RMPS method is being used extensively in several versions of the R-matrix approach, including the B-spline R-matrix (BSR) method (Zatsarinny 2006, Zatsarinny and Bartschat 2013) developed by Oleg Zatsarinny in my own group. As demonstrated in figures 9 and 10, RMPS can handle not only total ionisation cross sections, but with appropriate further processing of the raw output, one can reproduce fully-differential ionisation cross sections even for complex atoms like Ar.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Let me finish this part describing Phil's contributions to our field with one more personal impression. Developments like CCC and BSR could certainly be seen as stepping directly into Phil's territory. Instead of choosing 'territorial defense', Phil was always eager to collaborate and to lend his full support even to competitive efforts. He wrote numerous letters of support for my grant proposals, which undoubtedly were extremely helpful for my career.

Phil's first foray into the physics of electron–molecule collisions was a study of low-energy collisions with the nitrogen molecule performed in collaboration with SinfaiLam (Burke and Sinfailam 1970) in which they recovered the now well-known low-lying ²Π_g shape resonance. This calculation (which was styled as paper II) was performed using the single-centre expansion method proposed by Faisal (1970) in what was billed as paper I in the series. In paper III, Burke and Chandra (1972) developed a pseudo-potential representation of exchange, which was again applied to the N₂   ²Π_g shape resonance; indeed calculations on electron collisions with the nitrogen molecule represent something of a theme in Phil's molecular work, as he returned to this benchmark system to test each new theoretical development as it occurred (Buckley and Burke 1977, Noble et al 1982, Burke et al 1983, Gillan et al 1987, Gillan et al 1988, Gillan et al 1990, Gillan et al 1996, Burke et al 1999). Paper IV by Burke et al (1972) extended the use of single-centre expansions and (pseudo)-potentials to electron collisions with polyatomic molecules. While this paper laid the foundation for work on the use of the single-centre expansions calculations performed by Gianturco and others (e.g. (Curik et al 2001)) over several decades, Phil's work moved onto more rigorous formulations of the electron–molecule collision problem.

The 1977 paper by Burke, Mackey, and Shimamura (Burke et al 1977) developed a general R-matrix theory of low-energy electron scattering by diatomic molecules. The R-matrix equations for electron collisions by molecules are relatively straightforward extensions of those laid out in the previous section. In particular, the inner-region wavefunction expansion equation (3), asymptotic boundary condition, equation (7), and the fundamental equation of the R-matrix on the boundary, equation (17), all remain unchanged. The major difficulties are numerical and computational, as the introduction of a second atom into the scattering problem leads to both a loss of symmetry and the need to consider much larger partial-wave expansions.

At about this time Phil took a secondment to Daresbury Laboratory to become Head of the Theory Group. At Daresbury he assembled a small team to work on the electron scattering from diatomic molecules. In due course this team evolved to comprise Cliff Noble and Stefano Salvini, a PhD student at Queen's, and eventually me. While the formalism for the electron–molecule collision problem resembles the atomic problem, the numerical issues are much more challenging. The initial developments focused on suitable basis functions with which to represent both the target and the continuum wavefunctions. After a number of trial solutions, which largely did not make it into the scientific literature, Noble, Burke, and Salvini published their first electron–molecule R-matrix calculation (Noble et al 1982). This was built around the use of the Alchemy I diatomic electronic structure code originally developed at IBM San Jose (McLean 1971). Alchemy I used Slater-type orbitals (STOs) to represent the (target) wavefunction. This first work, which featured static exchange (SE) calculations for electron collisions with both H₂ and N₂, also used STOs to represent the continuum but concluded that the use of STO continuum functions inherently limited the range of energies that could be studied. While mentioning the possibility of using Gaussian-type orbitals (GTOs) to represent the continuum, the authors concluded that the use of numerically defined continuum orbitals was likely to be more successful. The subsequent paper by the same authors (Burke et al 1983) again considered electron collisions with N₂, this time at both the SE and the static exchange plus polarisation (SEP) level, and it included explicit consideration of vibrational effects within a Born–Oppenheimer treatment. This work used numerical basis functions, which were obtained as the solutions of appropriate radial equations. For these equations they used the isotropic part of a single-centre expansion of the N₂ target potential, V₀:

$$\begin{equation}\left(\frac{{\mathrm{d}}^{2}}{\mathrm{d}{r}^{2}}-\frac{{\ell }_{i}\left({\ell }_{i}+1\right)}{{r}^{2}}+{V}_{0}+{k}_{j}^{2}\right){u}_{i,j}\left(r\right)=\sum _{n}{\lambda }_{j,n}{P}_{i,n}\left(r\right),\end{equation} \tag{ 20 }$$

where the right-hand side ensured strict (Lagrange) orthogonalisation to the occupied target orbitals, which themselves were represented by single-centre expansion coefficients ${P}_{i,n}\left(r\right)$. Equation (20) was solved subject to the boundary conditions that the radial functions ${u}_{i,j}\left(r\right)$ were zero at the origin (r = 0) and had a zero boundary condition at the R-matrix boundary (r = a); this latter constraint necessitated the use of a Buttle (1967) correction. These continuum functions were found to provide a good basis set but proved to be rather numerous, as the Lagrange orthogonalisation constraint led to a new set being required for each ℓ.

The third paper by the same authors (Salvini et al 1984) reported SE and SEP calculations on the polar molecules LiH and CO. In this work numerical continuum basis functions for all ℓ were generated as simple solutions of equation (20) without an orthogonalisation constraint, i.e., with the right-hand side set to zero. This procedure provided a significantly smaller set of continuum functions, which were found to be sufficient for practical calculations. These functions, with the later addition of an optional, partial Lagrange orthogonalisation in case of severe linear dependence (Tennyson et al 1987), became the backbone of the UK diatomic R-matrix code (Gillan et al 1995, Morgan et al 1998) which is still in use today (Chakrabarti et al 2019).

Electron–molecule collisions lead to more complicated outer-region problems than electron–atom problems because of the many degenerate partial waves associated with each asymptotic state, as well as the presence of long-range potentials associated with the permanent moments of the target such as dipoles and quadrupoles. The original electron–molecule R-matrix calculations used the Burke and Schey (1962) propagator, but this is not particularly efficient for treating coupled many-channel problems. Baluja, Burke, and Morgan (Baluja et al 1982) reformulated the outer-region propagation problem by dividing the problems into sectors, which are then treated using a basis set of Legendre polynomials. This has proved to be a highly efficient method for solving the outer-region problem. The subsequent implementation by Morgan (1984) still forms the backbone of the UK molecular R-matrix (UKRmol) codes (Carr et al 2012, Mašín et al 2020).

I arrived at Daresbury in late 1982, just as Phil returned full-time to Queen's. While consultations with Phil were regular, the development of the electron–molecule collision code took an independent existence under the supervision of Cliff Noble, with help from myself and Lesley Morgan from Royal Holloway College who also joined the project at about the same time. Over the following decade, Cliff and Phil used the diatomic codes for a range of different calculations, including a series of increasingly sophisticated calculations on electron collisions with molecular oxygen (Noble and Burke 1992, Nordbeck et al 1993, Higgins et al 1994, Middleton et al 1994, Higgins et al 1995, Woste et al 1995, Noble et al 1996). Although Cliff and I spent most of our efforts developing and using the electron–molecule collision code, for instance performing the first R-matrix electron–molecule close-coupling calculations with Salvini (Tennyson et al 1984), this was not my job specification.

My position at Daresbury was actually associated with synchrotron science, as Daresbury then ran the SRS (synchrotron radiation source), and hence I was meant to study near-threshold photoionisation. With Phil and Cliff I therefore adapted the R-matrix photoionisation formalism of Burke and Taylor (1975) to the case of diatomic molecules and demonstrated its use for photoionisation of H₂ (Tennyson et al 1986) with a particular focus on the behavior at resonances.

In the 1970s Schneider, LeDourneuf, and Burke (Schneider et al 1979) had developed a sophisticated theory, which allows for the treatment of beyond Born–Oppenheimer effects in electron–diatomic molecule collisions by using an R-matrix representation of the nuclear coordinate as well as the electron one. This method allows for the treatment of a variety of processes, shown in figure 11, as well as vibrational excitation as part of a single calculation that treats the coupled electron and nuclear motion within the R-matrix hypersphere at the same time.

Zoom In Zoom Out Reset image size

In the late 1980s, Erhardt's group in Kaiserslautern reported a series of high-resolution measurements of near-threshold vibrational and rotational excitation cross sections for hydrogen halides such as HCl (Knoth et al 1989, Radle et al 1989), which built on similar earlier experimental studies by Rohr and Linder (1975, 1976). These experiments revealed cross sections that showed a series of complicated resonance-like structures in the region of each threshold to vibrational excitation; see figure 12. Phil, in collaboration with Lesley Morgan and a number of students, used the theory of Schneider, LeDourneuf, and Burke (Schneider et al 1979) to perform calculations on electron collisions with alkali halides, namely HF (Morgan and Burke 1988), HCl (Morgan et al 1990, McCartney et al 1990), and HBr (Fandreyer et al 1993). These studies showed that the unifying characteristic of these systems is the presence of a weakly bound anion state whose vibrationally excited states manifest themselves as resonances when nuclear motion is considered; see figure 13. The effect of these nuclear-excited resonances can only be properly considered by accounting for non-Born–Oppenheimer or non-adiabatic coupling between the electronic and nuclear motions. The double-R-matrix method (Schneider et al 1979) proved to be highly successful at doing this, as can be seen from figure 12.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The success of the diatomic code led Phil to consider the development of a general code capable of treating electron collisions with polyatomic molecules. Such a code would necessarily use GTOs to represent both the target and the continuum in the inner region, as these had proved outstandingly successful for treating the molecular electronic structure problem and thus formed the basis of essentially all available quantum chemistry codes upon which any polyatomic R-matrix code would be built. Phil proposed to the research council (SERC at that point) that his team in Belfast start on this enterprise, characterising it as a ten-year project. To his dismay, and I suspect his disgust (although he was always too gentlemanly to express such an opinion in public), this project was rejected as being too long-term. Eventually Lesley Morgan, in collaboration with Phil and myself, secured funding from the newly formed EPSRC (Engineering and Physical Sciences Research Council) to start building a GTO-based electron-polyatomic scattering code. However, this was not before, in the lacuna, the Bonn group had already launched such a code (Pfingst et al 1994). The new UKRmol code was built about the Sweden-molecule quantum chemistry code of Almlöf and Taylor (1984). Charles Gillan, Phil's student and then post-doc, was dispatched to Minnesota to work with Almlöf on this project. Lesley and Charles became the lead developers of this project (Morgan et al 1997, Morgan et al 1998). Phil, despite remaining interested in the development of a polyatomic code (Gillan et al 1995, Burke and Tennyson 2005), took something of a back seat. This was at least partly because Phil's interests in molecular processes had moved on in two different directions. One was the study of photoionisation by time-dependent fields, which is discussed the following section, and the other was the study of electron collisions with molecules weakly bounded (physisorbed) to surfaces.

Around 1995 Phil started a collaboration with John Inglesfield, his successor as Head of the Theory Group at Daresbury Laboratory but now at the University of Cardiff. John was a theoretical condensed matter physicist, and with Phil's student and then post-doc Katrina Higgins, they set about adapting the polyatomic R-matrix code to treat the problem of molecules on surfaces. This involved considering oriented molecules, something Phil had already done in collaboration with Karl Blum (Nordbeck et al 1993, Middleton et al 1994), and the inclusion of extra potentials, particularly the image potential due to the scattering electron (Burke et al 1999). This method was also used to develop an atomic theory of electron energy loss from transition metal oxides (Michiels et al 1997) and hence simulate spectra of CaO and NiO crystals (Jones et al 2000).

Phil became progressively less involved in the day-to-day developments of the UK molecular R-matrix codes, but his influence remained profound. One reason for this is that there is a transparent pathway in that developments made originally for electron–atom scattering find their way into electron–molecule problems, often after a decade or more. Thus, for example, the bound-state finding methodology Phil developed with Mike Seaton (Burke and Seaton 1984) was duly implemented in the molecular codes (Sarpal et al 1991) and remains in regular use (Little and Tennyson 2013). Similarly, Phil developed the use of polarised pseudostates to represent excited target states, including target continuum states neglected in the close-coupling expansion for electron–atom scattering problems (Burke and Mitchell 1974, Lan et al 1976). More than two decades later, he applied it to electron–molecule scattering by again returning to the electron–N₂ problem (Gillan et al 1988). Subsequently, in 1996, with Klaus Bartschat and others he developed the more powerful RMPS method (Bartschat et al 1996; see the previous section). The RMPS method was implemented in the molecular codes about eight years later by Jimena Gorfinkiel and myself (Gorfinkiel and Tennyson 2004, Gorfinkiel and Tennyson 2005). Finally, as an algorithmic example, Phil, with his wife Val and Kevin Dunseath, developed a very much more efficient algorithm for constructing the inner-region Hamiltonian matrix based on expanding the total wavefunction in target configurations rather than in individual target states and taking advantage of the special status of the scattered electron in the collisional wavefunction. This development became the heart of the RMATRX II atomic code (Burke et al 1994). Inspired by this and after extensive discussions with Phil, I developed a similar algorithm for the molecular code known as SCATCI (Tennyson 1996), recently updated for MPI computers (Al-Refaie and Tennyson 2017). This has become the bedrock of both the UKRMol (Carr et al 2012) and UKRMol+ (Mašín et al 2020) electron–molecule scattering codes and hence also their easy-to-use commercially available versions known as Quantemol-N (Tennyson et al 2007) and Quantemol electron scattering (QEC) (Cooper et al 2019), respectively.

In the late 1980s, Phil realized that R-matrix theory could be combined with the Floquet-Fourier ansatz to treat atomic multiphoton processes. This can be thought to include multi-photon ionisation (probably the main application), laser-assisted electron–atom collisions, and harmonic generation. This development was well-timed to help in the analysis of a range of experiments that took advantage of the new laser technologies that were (and continue to be) driving many advances in atomic and molecular scattering processes. This major advance in application of the R-matrix theory came about through a collaboration that Phil fostered with the Belgian groups of Professor Charles Joachain (Université Catholique de Louvain) and Professor Francken (Université Libre de Bruxelles). It was around this time that Phil also co-authored a book titled 'theory of electron-atomic collisions' (Joachain and Burke 1994) with Joachain.

The details of what became quickly known as R-matrix-Floquet theory were discussed in Burke et al (1990) as well as in other sources (Burke et al 1991, Dorr et al 1992). As in previous sections, here we very briefly discuss only the salient theoretical points, again basing our descriptions on Phil's book (Burke 2011). In R-matrix-Floquet theory, the main assumption made about the interaction of an intense laser field with the atom is that many cycles of the laser field are involved, and hence the method is not appropriate for few-cycle or extremely short pulses. Explorations of the range of validity of this approximation are discussed below. The laser field is treated classically and assumed to be monochromatic, linearly polarized, and spatially homogeneous. The corresponding electric field vector can be written as

$$\begin{equation}\mathcal{E}\left(t\right)=-\frac{1}{c}\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{A}\left(t\right)=\hat{{\epsilon}}{\mathcal{E}}_{0}\enspace \mathrm{cos}\enspace \omega t,\end{equation} \tag{ 21 }$$

where $\hat{{\epsilon}}$ is a unit vector along the laser polarization direction, $\mathcal{E}\left(t\right)$ is the electric field strength, ω is the angular frequency, and A(t) is the vector potential. The atomic system in the presence of the external laser field is described by the time-dependent Schrödinger equation

$$\begin{align}\hfill & \left[{H}_{N+1}+\frac{1}{c}\mathbf{A}\left(t\right)\cdot {\mathbf{P}}_{N+1}+\frac{N+1}{2{c}^{2}}{\mathbf{A}}^{2}\left(t\right)\right]\tilde {{\Psi}}\left({\mathbf{X}}_{N+1},t\right)\hfill \\ \hfill & =\mathrm{i}\frac{\partial }{\partial t}\tilde {{\Psi}}\left({\mathbf{X}}_{N+1},t\right),\hfill \end{align} \tag{ 22 }$$

where H_N+1 is the non-relativistic Hamiltonian (see equation (2)) of the (N + 1)-electron atomic system in the absence of the laser field, and P_N+1 is the total electron momentum operator. The tilde on the time-dependent wavefunction $\tilde {{\Psi}}$ distinguishes this function from the time-independent wave function Ψ, which will be considered later.

The next key step in this theory is the representation of the time-dependent wavefunction $\tilde {{\Psi}}\left({\mathbf{X}}_{N+1},t\right)$ in terms of time-independent wavefunctions Ψ_n(X_N+1) via the use of the Floquet-Fourier expansion given by

$$\begin{equation}\tilde {{\Psi}}\left({\mathbf{X}}_{N+1},t\right)={\mathrm{e}}^{-\mathrm{i}Et}{{\Sigma}}_{n=-\infty }^{\infty }{\mathrm{e}}^{-\mathrm{i}n\omega t}{{\Psi}}_{n}\left({\mathbf{X}}_{N+1}\right).\end{equation} \tag{ 23 }$$

After substituting this expansion into the time-dependent Schrödinger equation, an infinite set of coupled time-independent equations for Ψ_n(X_N+1) results, where the expansion over n is, in practice, truncated at some (relatively) small and finite n value. Physically, the positive and negative terms correspond to absorption or emission of n photons.

In R-matrix-Floquet theory, as in other R-matrix formalisms, space is configured into different regions as shown in figure 14. In the internal region, the dipole length gauge is used in the time-dependent Schrödinger equation via a unitary transformation, since the dipole interaction $\mathcal{E}\left(t\right)\cdot {R}_{N+1}$ tends to zero as R_N+1 → 0. The derivation of the internal-region Hamiltonian is discussed in detail in Burke et al (1991). When written as a matrix equation, the TISE takes the form

$$\begin{equation}\left({\mathbf{H}}_{F}^{L}-E\mathbf{I}\right){\mathbf{\Psi }}^{L}=0,\end{equation} \tag{ 24 }$$

where the Hamiltonian is written as

$$\begin{equation}{\mathbf{H}}_{F}=\left(\begin{matrix}\hfill \ddots \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \\ \hfill {D}_{N+1}\hfill & \hfill {H}_{N+1}-\left(n-1\right)\omega \hfill & \hfill {D}_{N+1}\hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill {D}_{N+1}\hfill & \hfill {H}_{N+1}-n\omega \hfill & \hfill {D}_{N+1}\hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill {D}_{N+1}\hfill & \hfill {H}_{N+1}-\left(n+1\right)\omega \hfill & \hfill {D}_{N+1}\hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \ddots \hfill \end{matrix}\right).\end{equation} \tag{ 25 }$$

Expressed in this form, the Hamiltonian is a block-tridiagonal matrix, where the off-diagonal elements are defined by the dipole length operator ${D}_{N+1}=\frac{1}{2}{\mathcal{E}}_{0}\hat{{\epsilon}}.{\mathbf{R}}_{N+1}$. The solution in the internal region then proceeds in a similar manner to that described in section 2.

Zoom In Zoom Out Reset image size

In the external region, the ejected (or scattered) electron is described using the velocity gauge, while the remaining N electrons are described using the length gauge. This is possible since the outer electron and the N inner electrons occupy different regions of space and are distinguishable. The expansion of the wavefunction and subsequent insertion into the time-dependent Schrödinger equation are again discussed in detail in Burke et al (1991). There, it is shown that one can recast the external-region solution through a set of coupled second-order differential equations that are very similar to equation (4).

In the asymptotic region, the ejected (or scattered) electron may be described in either the velocity gauge or the acceleration frame of reference, using the Kramers–Henneberger transformation. Various numerical methods for the solution in this region have been described and a simplified analysis has been given by Charlo et al (1998) and Terao-Dunseath and Dunseath (2001).

Multiphoton ionisation rates of a number of atomic systems have been systematically explored using R-matrix-Floquet approaches. These include studies of resonance-enhanced multi-photon ionisation, often known as REMPI. As an example, we show in figure 15 the ionisation rate for Ar subjected to 390 nm laser light, as a function of the laser intensity (van der Hart 2006). The various peaks evident in the ionisation rate correspond to resonances such as the 3s²3p⁵5s  ¹P^o state that is indicated. The REMPI rate found by the R-matrix-Floquet calculations is considerably higher than that predicted by the Ammosov–Delone–Krainov (ADK) tunneling model (Ammosov et al 1986).

Zoom In Zoom Out Reset image size

Results from R-matrix-Floquet calculations have also been compared with results from direct propagation of the time-dependent Schrödinger equation for helium by Parker et al (2000). Figure 16 shows the ionisation rate of helium as a function of laser intensity at a laser wavelength of 248.6 nm. The sharp peaks correspond to resonances of various intermediate bound states. At low laser intensities, five photons are sufficient to ionise helium, but at higher intensities, six photons are needed. This is due to the increase in the ionisation threshold by the ponderomotive energy given by ${E}_{\mathrm{p}}={\mathcal{E}}_{0}^{2}/4{\omega }^{2}$. The solid lines in figure 16 correspond to the R-matrix-Floquet calculations while the solid circles are results from the time-dependent Schrodinger equation calculations of Parker et al (2000). Excellent, but not perfect, agreement is found between the predictions from the two methods. Near the resonances, time-dependent methods require propagation in time to hundreds of field periods to obtain highly accurate results. Conversely, the position of some resonances will be sensitive to the close-coupling expansion employed in the R-matrix-Floquet calculations, which can result in slightly different resonance energies and positions. However, the maximum difference between the ionisation rates shown in figure 16 is still less than 20%.

Zoom In Zoom Out Reset image size

R-matrix-Floquet theory has also been used to explore the phenomenon known as laser-induced degenerate states (LIDS), where an autoionising state in the electron–ion collision system and a continuum state corresponding to the atom 'dressed' by the laser field become degenerate at a particular laser energy and intensity. This can be thought of as a 'double pole' in the laser-assisted electron–ion collision S-matrix. Figure 17 shows such a feature in Ar one-photon ionisation that was explored by Latinne et al (1995). The figure shows the trajectories in the complex energy plane as a function of laser intensity. The evolution of each state in this plane is traced as a function of laser intensity. The zero-field position of the Ar ground state is indicated as E_g (at 0.578 16 a.u.) and the circles indicate the zero-field position of the autoionising state (which depends on the laser frequency ω). For each frequency, two curves adiabatically connect the ground state or the autoionising state (the frequencies of each line are indicated on the figure). At some intermediate value of the detuning frequency d = E_a − E_g − ω, structures are visible, about which the curves of the ground state and the autoionising state exchange roles. At this point in the complex plane the complex energies of the two states are exactly degenerate. A much more detailed discussion of the implications of these fascinating phenomena is given in Latinne et al (1995), including the possibility of trapping population in an intermediate state as a function of laser intensity, if the laser parameters can be varied in a precise manner.

Zoom In Zoom Out Reset image size

The R-matrix-Floquet theory has also been applied to other phenomena, such as harmonic generation (Gebarowski et al 1997a, Gebarowski et al 1997b), laser-assisted electron scattering (Dorr et al 1995a, Terao-Dunseath and Dunseath 2001, Terao-Dunseath et al 2001) and above-threshold ionisation (Dorr et al 1993, Vinci et al 2000a), as well as studies of other multiphoton ionisation rates in a number of atomic systems (Dorr et al 1992, Dorr et al 1993, Purvis et al 1993, Cyr et al 1997, Glass et al 1997, Kylstra et al 1998, Glass and Burke 2000b, Glass and Burke 2000a), including negative ions (Vinci et al 2000a, Dorr et al 1995b, Glass et al 1998, Vinci et al 2000b, Vinci et al 2003). One of us (HWvdH) also applied R-matrix-Floquet theory to 'two-colour' multiphoton processes, that is, ionisation by two laser fields of different frequencies (van der Hart 1996). This approach has also found some recent applications (Costa i Bricha et al 2004). Extensive discussions of all these topics are given in Burke (2011).

We close this section by briefly discussing an extension of R-matrix-Floquet theory to molecular multiphoton processes. This was the PhD topic of my thesis (JC). The research was performed in the late 1990s, before and after Phil's official retirement in 1998. Retirement certainly did not slow Phil down, and we suspect that the release from teaching and official university duties was privately welcomed by Phil. The treatment of molecular multiphoton processes was the result of a combination of the R-matrix electron–molecule collision codes and the multiphoton parts of the atomic R-matrix-Floquet codes. Phil had the insight that the internal region electron–molecule codes could be adapted with relatively little change, along with a code that computed the molecular analogue of the dipole matrix elements discussed earlier. Key contributions were also made here by David Glass and Katrina Higgins, at that time both post-doctoral scholars at QUB.

Space is again partitioned in a similar way to that indicated in figure 14, apart from the use of a multi-centre expansion in the internal region. A key component of this development was the need of a transformation from the molecular fixed coordinate system (where the z-axis lies along the internuclear axis) to a laser-field coordinate system where the z'-axis lies along the laser-field polarization direction. This transformation, implemented through the use of Euler angles and standard rotation algebra, is required to enable the computation of various scattering quantities and ionisation rates.

As in our previous sections, we leave the details of the development of molecular R-matrix-Floquet theory to Phil's book (Burke 2011) and the journal papers from that time (Colgan et al 1998, Colgan et al 2000, Burke et al 2000, Colgan et al 2001). This method was used to examine multiphoton ionisation rates in the H₂ molecule. For example, figure 18 shows the ionisation rates of H₂ as a function of the laser frequency for three different laser intensities. One-state (solid line) and two-state (dashed line) approximations are compared. As in atomic systems, we find REMPI peaks, corresponding to three-photon excitation of Rydberg bound states followed by ionisation through absorption of a fourth photon. The position of these peaks changes with laser intensity due to the ponderomotive shift of the ionisation threshold. Further exploration also considered the position and widths of such REMPI peaks as a function of the internuclear separation of H₂.

Zoom In Zoom Out Reset image size

The influence of Phil Burke's research into AMO physics extends beyond those areas, already discussed, where he had a direct involvement. To further demonstrate Phil's lasting legacy we include a brief discussion of important work that has taken place since his retirement, none of which would have been possible without the decades-long development of R-matrix theory and the associated computer codes that he led.

The formalism outlined in section 2 can be extended into the time domain with the explicit inclusion of the time-dependent potential of the laser field. In the semi-relativistic case, the TDSE for an N-electron system may be written as

$$\begin{equation}\mathrm{i}\frac{\partial }{\partial t}{\Psi}\left({\mathbf{X}}_{N},t\right)=\left[{H}_{\mathrm{N}}+{H}_{\mathrm{S}\mathrm{O}}+{H}_{\mathrm{D}}\left(t\right)\right]{\Psi}\left({\mathbf{X}}_{N},t\right),\end{equation} \tag{ 26 }$$

where H_N is the time-independent, non-relativistic, atomic Hamiltonian as defined in equation (2), and H_SO is the Breit–Pauli spin–orbit interaction applied to all electrons. The term H_D(t) is the time-dependent dipole-operator, which describes the interaction of the laser with all electrons,

$$\begin{equation}{H}_{\mathrm{D}}\left(t\right)=\mathbf{E}\left(t\right)\cdot \sum _{i=1}^{N}{\mathbf{r}}_{i},\end{equation} \tag{ 27 }$$

where r_i is the position vector of the ith electron.

Unlike the R-matrix-Floquet approach discussed in section 4, in time-dependent R-matrix theory we take some initial wavefunction (provided by a time-independent R-matrix calculation) and propagate it explicitly in time, applying, at each time step, the atomic and laser potentials, and thereby computing ${\Psi}\left({\mathbf{X}}_{N},t\right)$ for each time, t ⩾ 0.

Phil and Val published a paper detailing time-dependent R-matrix (TDRM) theory in 1997—shortly before Phil's retirement from Queen's—and they demonstrated its application to multiphoton ionisation of a 1D model system (Burke and Burke 1997). It then took a full decade for two independently developed implementations of TDRM to be published (Guan et al 2007, van der Hart et al 2007). There followed a number of publications showcasing the capabilities of TDRM to account for multielectron dynamics in multiphoton processes (Guan et al 2008a, Guan et al 2009, Guan et al 2008b, Lysaght et al 2009a, Hutchinson et al 2010, Brown et al 2012, Lysaght et al 2009b, Lysaght et al 2009c).

The advantages of TDRM theory are offset somewhat, however, by the sheer computational cost of implementation. While applications of TDRM to multiphoton processes are impressive theoretical accomplishments, until quite recently the science they described still lagged some way behind both experiment and the best model approaches. Most recently however, a new implementation, the R-matrix with time-dependence (RMT) approach, has been able to surmount the prohibitive computational barriers to become the only method of its kind capable of describing general, multielectron systems in strong fields (Moore et al 2011, Clarke et al 2018, Wragg et al 2019). RMT has been applied to strong field processes such as high harmonic generation (Hassouneh et al 2014, Brown and van der Hart 2016), attosecond transient absorption spectroscopy (Ding et al 2016) and strong-field rescattering (Hassouneh et al 2015). More recently, the method has been extended to account for processes in arbitrarily polarised laser pulses (Clarke et al 2018)—allowing investigations into electron vortices (Armstrong et al 2019), and attosecond angular streaking (Armstrong et al 2020)—semi-relativistic phenomena (Wragg et al 2019), and molecular targets (Brown et al 2020).

6.1. Computer Physics Communications

6.2. Collaborative Computational Projects

In this manuscript, we have tried to summarize and highlight the extraordinary contributions of Philip George Burke CBE FRS, not only to the computational treatment of electron and photon collisions with atoms, ions, and molecules, but also to the organizational structure of atomic, molecular, and optical physics in the UK and worldwide. Phil's leadership was outstanding, and his legacy will last for many years to come.

The work of KB has been supported for many years by the United States National Science Foundation, most recently under Grant Nos. PHY-1803844 and OAC-1834740, as well as supercomputer support through the eXtreme Science and Engineering Discovery Environment (XSEDE) program under Grant No. PHY-090031. The RMT and UKRMol+ codes are part of the UK-AMOR suite, and can be obtained for free at reference (UK-AMOR (2019) code repository). Parts of this work benefited from computational support by CoSeC, the Computational Science Center for Research Communities, through CCPQ. The authors acknowledge funding from the EPSRC under Grants EP/P022146/1, EP/P013953/1 and EP/R029342/1. Parts of this work relied on the ARCHER UK National Supercomputing Service (www.archer.ac.uk), for which access was obtained via the UK-AMOR consortium funded by EPSRC grant EP/R029342/1. Part of this work was supported by the Physics & Engineering Materials (PEM) program of the US Department of Energy through the Los Alamos National Laboratory. Los Alamos National Laboratory is operated by Triad National Security, LLC, for the National Nuclear Security Administration of U.S. Department of Energy (Contract No. 89233218NCA000001). The authors have confirmed that any identifiable participants in this study have given their consent for publication

All of us, and certainly many of our colleagues, are most grateful for the interactions we have had with Phil Burke over many years. Phil was not only a 'gentleman'—he was also a 'gentle man'. We will miss him, but will do our best to carry on his legacy.