Efficient Evaluation of a Class of Integrals Commonly Encountered in the Study of Two-Electron Systems with One or Both Electrons Bound to a Nucleus

A system that consists of two electrons and an infinitely heavy nucleus is considered. At least one electron is bound to a stationary nucleus; the other electron can be bound or unbound. As a consequence of a special broken rotational symmetry, the wavefunction separates into a finite linear combination of terms, each of which factorizes into a function of the external (Euler) angles and a function of the interparticle coordinates r1=|r→1|, r2=|r→2|, and r3=|r→1−r→2|, with r3 the separation of electrons located at distances r1 and r2 from the nucleus. A theoretical treatment often proceeds by introducing a finite basis composed of products of powers of r1, r2, and r3 multiplied by exponentials in r1 and r2. After integrating over the external angles, most matrix elements of interest can be expressed in terms of three-dimensional integrals whose integrands are individual functions of r1, r2, and r3 multiplied by the Legendre polynomial Pl(r^1·r^2), with l a nonnegative integer related to the internal and total angular momenta. Integration over r3 is straightforward. However, since the derivative of r3 with respect to r1 or r2 is undefined at r1=r2, the remaining two-dimensional integral is highly prone to numerical roundoff error, particularly if a large basis is employed. It is shown that roundoff error can be managed by a rearrangement of the integrand, and the use of numerically stable recurrence relations.