# Complete Numerical Solution of the Temkin-Poet Three-Body Problem

## Abstract

Although the convergent close-coupling (CCC) method has achieved unprecedented success in obtaining accurate theoretical cross sections for electron-atom scattering, it generally fails to yield converged energy distributions for ionization. Here we report converged energy distributions for ionization of $`\mathrm{H}(1\mathrm{s})`$ by numerically integrating Schrödinger’s equation subject to correct asymptotic boundary conditions for the Temkin-Poet model collision problem, which neglects angular momentum. Moreover, since the present method is complete, we obtained convergence for all transitions in a single calculation (excluding the very highest Rydberg transitions, which require integrating to infinitely-large distances; these cross sections may be accurately obtained from lower-level Rydberg cross sections using the $`1/n^3`$ scaling law). Complete results, accurate to 1%, are presented for impact energies of 54.4 and 40.8 eV, where CCC results are available for comparison.

The Temkin-Poet (TP) model  of electron-hydrogen scattering is now widely regarded as an ideal testing ground for the development of general methods intended for the full three-body Coulomb problem. Although only $`s`$-states are included for both projectile and atomic electrons, this model problem still contains most of the features that make the real scattering problem hard to solve. Indeed, even in this simplified model, converged energy distributions for ionization can not generally be obtained via the close-coupling formalism . Any general method that can not obtain complete, converged results for this model problem will face similar difficulties when applied to the full electron-hydrogen system. Therefore we believe it is essential to develop a numerical method capable of solving the TP model completely before angular momentum is included. Here we report such a method. Complete, precision results for $`e^{}+\mathrm{H}(1\mathrm{s})`$, accurate to 1%, are presented for total energies of 3 and 2 Rydbergs (Ryd). Atomic units (Ryd energy units) are used throughout this work unless stated otherwise.

Our numerical method may be summerized as follows. The model Schrödinger equation is integrated outwards from the atomic center on a grid of fixed spacing $`h`$. The number of difference equations is reduced each step outwards using an algorithm due to Poet , resulting in a propagating solution of the partial-differential equation. By imposing correct asymptotic boundary conditions on this general, propagating solution, the particular solution that physically corresponds to scattering is obtained along with the scattering amplitudes.

The Schrödinger equation in the TP model is given by

$$\left(\frac{^2}{x^2}+\frac{^2}{y^2}+\frac{2}{\mathrm{min}(x,y)}+E\right)\mathrm{\Psi }(x,y)=0,$$
(1)   
with boundary conditions

$$\mathrm{\Psi }(x,0)=\mathrm{\Psi }(0,y)=0$$
(2)   
and symmetry condition

$$\mathrm{\Psi }(y,x)=\pm \mathrm{\Psi }(x,y),$$
(3)   
depending on whether the two electrons form a singlet ($`+`$) or triplet ($``$) spin state. Eq. (1) is separable in the two regions $`xy`$ and $`xy`$. Because of the symmetry condition (3), we can solve Eq. (1) in just one of these regions and this is sufficient to determine all of the scattering information. For brevity, we do not explicitly indicate the total spin since the singlet and triplet cases require completely separate calculations. For $`xy`$, the wave function may be written

$`\mathrm{\Psi }(x,y)=\psi _{ϵ_i}(y)e^{ik_{ϵ_i}x}+{\displaystyle \underset{j=1}{\overset{\mathrm{}}{}}}C_{ϵ_ji}\psi _{ϵ_j}(y)e^{ik_{ϵ_j}x}`$ (4)   
$`+{\displaystyle _0^{\mathrm{}}}𝑑ϵ_bC_{ϵ_bi}\psi _{ϵ_b}(y)e^{ik_{ϵ_b}x}.`$ (5)   
The $`\psi _ϵ`$ are bound and continuum states of the hydrogen atom with zero angular momentum:

$`\psi _ϵ(y)=ye^{qy}{}_{1}{}^{}F_{1}^{}(11/q,2;2qy).`$ (6)   
Here $`q^2=ϵ`$, where $`ϵ`$ is the inner electron energy, and $`{}_{1}{}^{}F_{1}^{}`$ is the confluent hypergeometric function. The momenta in (5) are fixed by energy conservation according to

$$ϵ_i+k_{ϵ_i}^2=ϵ_j+k_{ϵ_j}^2=ϵ_b+k_{ϵ_b}^2=E,$$
(7)   
where $`E>0`$ is the total energy. The $`C_{ϵi}`$ are related to S-matrix elements by normalization factors:

$$\mathrm{S}_{ϵ_ji}=\left(\frac{k_{ϵ_j}}{k_{ϵ_i}}\right)^{1/2}\left(\frac{j}{i}\right)^{3/2}C_{ϵ_ji}$$
(8)   
for discrete transitions and

$$\mathrm{S}_{ϵ_bi}=\left(\frac{k_{ϵ_b}}{k_{ϵ_i}}\right)^{1/2}\left(\frac{1}{i}\right)^{3/2}\left[\frac{1e^{2\pi /k}}{4k}\right]^{1/2}C_{ϵ_bi}$$
(9)   
for ionization, where $`k=\sqrt{ϵ_b}`$. Cross sections are obtained from S-matrix elements in the usual manner.

To convert the partial-differential equation (1) into difference equations we impose a grid of fixed spacing $`h`$ and approximate derivatives by finite differences. After applying the Numerov scheme in both the $`x`$ and $`y`$ directions, our difference equations have the form 

$$𝐀^{(i)}\overline{\mathrm{\Psi }}^{(i1)}+𝐁^{(i)}\overline{\mathrm{\Psi }}^{(i)}+𝐂^{(i)}\overline{\mathrm{\Psi }}^{(i+1)}=\mathrm{𝟎},$$
(10)   
Here we have collected the various $`\mathrm{\Psi }_j^{(i)},j=1,2,\mathrm{},i`$, where $`\mathrm{\Psi }_j^{(i)}\mathrm{\Psi }(x=ih,y=jh)`$, into a vector $`\overline{\mathrm{\Psi }}^{(i)}`$. The matrices $`𝐀^{(i)}`$, $`𝐁^{(i)}`$ and $`𝐂^{(i)}`$ are completely determined by the formulas given by Poet .

At each value of $`i`$ we can solve our equations if we apply symbolic boundary conditions at $`i+1`$ \[solve for $`\mathrm{\Psi }_j^{(i)}`$ in terms of $`\mathrm{\Psi }_j^{(i+1)}`$ ($`j=1,2,\mathrm{},i`$)\]. This procedure yields a propagation matrix $`𝐃^{(i)}`$:

$$\overline{\mathrm{\Psi }}^{(i)}=𝐃^{(i)}\overline{\mathrm{\Psi }}^{(i+1)}.$$
(11)   
We can obtain a recursion relation for $`𝐃^{(i)}`$ by using (11) to eliminate $`\overline{\mathrm{\Psi }}^{(i1)}`$ from equation (10):

$$\left[𝐁^{(i)}+𝐀^{(i)}𝐃^{(i1)}\right]\overline{\mathrm{\Psi }}^{(i)}=𝐂^{(i)}\overline{\mathrm{\Psi }}^{(i+1)}.$$
(12)   
Comparing (12) with (11),

$$𝐃^{(i)}=\left[𝐁^{(i)}+𝐀^{(i)}𝐃^{(i1)}\right]^1𝐂^{(i)}.$$
(13)   
Thus each $`𝐃^{(i)}`$ is determined from the previous one ($`𝐃^{(1)}`$ can be determined by inspection).

In the asymptotic region, the form of the wave function is known and is given in terms of the $`C_{ϵi}`$ by

$`\overline{\mathrm{\Psi }}^{(i)}𝐈^{(i)}+𝐑^{(i)}𝐂.`$ (14)   
Here the matrix $`𝐈^{(i)}`$ contains the incident part of the asymptotic solution while $`𝐑^{(i)}`$ contains the reflected part. The asymptotic solution is identical to the full solution, Eq. (5), except that the quadrature over the continuum extends only up to the total energy $`E`$. The infinite summation over discrete channels is truncated to some finite integer $`N_d`$ and the quadrature over the two-electron continuum is performed prior to matching by first writing the $`C_{ϵ_bi}`$ as a power series in $`ϵ_b`$:

$$C_{ϵ_bi}\underset{n=1}{\overset{N_c}{}}\mathrm{c}_{ni}ϵ_b^n.$$
(15)   
The matching procedure then determines the (in practice, much smaller set of) coefficients $`\mathrm{c}_{ni}`$, rather than the $`C_{ϵ_bi}`$ directly, which eliminates ill conditioning .

To extract an $`N\times N`$ coefficient matrix, where $`N=N_d+N_c`$, we need only $`N`$ of the $`i`$ equations (11). Alternatively, one may use all $`i`$ equations as in Poet . In this case, the system of equations is overdetermined. Nevertheless, a solution can be found by the standard method of minimizing the sum of the squares of the residuals \[the differences between the left- and right-hand sides of equations (11)\]. Previously , we found that the least-squares method is generally stabler than keeping any subset of just $`N`$ equations (11).

Our numerical method is stable and rapidly convergent. For a given grid spacing $`h`$, we established convergence in propagation distance by performing the matching every 40 a.u. until convergence was obtained. At each matching radius, both the number of discrete channels $`N_d`$ and the number of expansion functions for the continuum $`N_c`$ were varied to obtain convergence. Finally, the entire calculation was repeated for a finer grid (using one-half the original grid spacing $`h`$).

The biggest advantage of having a general, propagating solution is that once the grid spacing is chosen, a “single” calculation is all that is needed to establish convergence for the remaining numerical parameters. This is because the D-matrix, the calculation of which consumes nearly all the computational effort, is independent of asymptotic boundary conditions. Thus, in a typical calculation, the same D-matrix is used for, e.g., $`N_c=0,1,\mathrm{},9`$ while $`N_d`$ runs from 1 to 30. This would have required 300 completely separate calculations (each taking about the same time as our “one” calculation) had we solved the original global matrix equations (10).

We have performed complete calculations for electrons colliding with atomic hydrogen at impact energies of 54.4 and 40.8 eV (total energies of 3 and 2 Ryd, respectively). In Table I, we present our calculated cross sections for $`e^{}+\mathrm{H}(1s)e^{}+\mathrm{H}(ns),n8`$. The grid spacing is $`h=1/5`$ a.u. (results using one-half this spacing differed by less than 0.1% for discrete excitations and 0.5% for elastic scattering). One of the advantages of our direct approach is that we are able to obtain the amplitudes for higher-level (Rydberg) transitions as easily as those for low-level excitations, provided the matching radius is large enough to enclose the final Rydberg state. This is in contrast to some other approaches, such as the CCC, which lose accuracy for higher-level transitions.

In Figures 1-4, we present our results (labeled FDM for finite-difference method) for the single-differential cross section (SDCS). For a total energy of 3 Ryd, 240 a.u. proved to be a sufficient matching radius to get convergence of the SDCS and for E = 2 Ryd, a radius of 360 a.u. was required. The SDCS is more sensitive to the number of expansion functions for the continuum than the other observables, particularly about $`ϵ_b=E/2`$. Nevertheless, convergence to better than 1% was readily obtained using 7-8 functions (the largest discrepancy in the SDCS between $`N_c=7`$ and $`N_c=8`$ was smaller than 0.3%; even using just 6 expansion function gave results accurate to 1%).

Also shown in Figs. 1-4 are the results of convergent close-coupling (CCC) calculations . The CCC method of Bray  employs a “distinguishable electron” prescription, which produces energy distributions that are not symmetric about $`ϵ_b=E/2`$. Stelbovics  has shown that a properly symmetrized CCC amplitude yields SDCS that are symmetric about $`E/2`$ as well as being four times larger at $`ϵ_b=E/2`$ than those assuming distinguishable electrons. (Note that our singlet FDM results at $`ϵ_b=E/2`$ are about four times larger than the corresponding CCC results.) Other than making the energy distributions symmetric, it is clear from the figures that symmetrization (coherent summation of the CCC amplitudes $`C_{ϵ_bi}`$ and $`C_{ϵ_ai}`$, where $`ϵ_a=Eϵ_b`$, which correspond to physically indistinguishable processes) will significantly affect only singlet scattering (and then only near $`ϵ_b=E/2`$), since $`C_{ϵ_bi}`$ is practically zero for $`ϵ_b>E/2`$. For singlet scattering, the CCC oscillates about the true value of SDCS, except near (and beyond) $`ϵ_b=E/2`$. CCC results for triplet scattering, on the other hand, are in very good agreement with our results for $`0ϵ_bE/2`$.

Some very recent results from Baertschy et al.  have also been included in the figures. Baertschy et al. rearrange the Schrödinger equation to solve for the outgoing scattered wave. They use a two-dimensional grid like ours, but scale the coordinates by a complex phase factor beyond a certain radius where the tail of the Coulomb potential is ignored. As a result, the scattered wave decays like an ordinary bound state beyond this cut-off radius, which makes the asymptotic boundary conditions very simple. By computing the outgoing flux directly from the scattered wave at several large cut-off radii, and extrapolating to infinity, they obtain the single-differential ionization cross section without having to use Coulomb three-body boundary conditions. This method, called exterior complex scaling (ECS), has just been extended to the full electron-hydrogen ionization problem . It is seen from Figs. 1-4 that the ECS results are in good agreement with our FDM results except when the energy fraction $`ϵ_b/E`$ approaches 0 or 1. Baertschy et al.  note that their method may be unreliable as $`ϵ_b`$ approaches 0 or $`E`$ due to “contamination” of the ionization flux by contributions from discrete excitations.

We note also the recent work of Miyashita et al. , who have presented SDCS for total energies of 4, 2, and 0.1 Ryd using two different methods. One produces an asymmetric energy distribution similar to that of CCC while the other gives a symmetric distribution. Both contain oscillations. The mean of their symmetric curve at $`E=2`$ Ryd (40.8 eV impact energy) is in reasonable agreement with our calculations.

In conclusion, we have presented complete, precision results for the Temkin-Poet electron-hydrogen scattering problem for impact energies of 54.4 and 40.8 eV. It may be possible to improve the speed of the present method by using a variable-spaced grid, like that used by Botero and Shertzer  in their finite-element analysis (this would greatly reduce storage requirements as well). Once we have optimized our code for this simplified model we will proceed to include angular momentum. When angular momentum is included, the ionization boundary condition is no longer separable and this is the major challenge for generalizing the present approach to the full electron-hydrogen scattering problem.

The authors gratefully acknowledge the financial support of the Australian Research Council for this work.