THE 1958 PEKERIS-ACCAD-WEIZAC GROUND-BREAKING COLLABORATION THAT COMPUTED GROUND STATES OF TWO-ELECTRON ATOMS (AND ITS 2010 REDUX)

CHRISTOPH KOUTSCHAN^{∗} AND DORON ZEILBERGER^{∗∗}

We have come a long way

In order to appreciate how good we as mathematicians and scientists have it today, with extremely fast hardware and lots and lots of memory, as well as with readily available high-level software, both for numeric and symbolic computation, it may be a good idea to go back to the early days of electronic computers and carefully examine, as a case study, a problem that was considered a huge challenge back then, and compare notes. We chose C.L. Pekeris’ [9] 1958 seminal work^{1} on the ground state energies of two-electron atoms. In particular, we will do all computations ab initio with today’s software and hardware.

Schrödinger

Let’s recall the (time-independent) Schrödinger equation for the state function (alias wave function) ψ(x, y, z) of a one-electron atom with a stationary nucleus (see, for example, [8] Eq. (30-1) with N = 1), in atomic units:

(\frac{∂^{2}}{∂x^{2}} + \frac{∂^{2}}{∂y^{2}} + \frac{∂^{2}}{∂z^{2}} + 2 (E + \frac{Z}{r} )) ψ(x, y, z) = 0,

where Z denotes the nuclear charge, E the energy of the system, and r = \sqrt{x^{2} + y^{2} + z^{2}} is the distance of the electron to the nucleus. Schrödinger’s solution of this eigenvalue problem is one of the greatest classics of modern physics, familiar to all physics (and chemistry, but unfortunately not math) students, using separation of (dependent) variables, and getting explicit and exact results for the eigenvalues (the possible energy levels E) and even for the corresponding eigenfunctions ψ. Since the latter (or more precisely: their squares) are interpreted as probability distributions, certain restrictions have to be imposed on ψ; in particular, the integral of ψ^{2} over the whole domain must be finite. The eigenvalues then are exactly those values of E for which the Schrödinger equation admits such a solution. It turns out that these eigenfunctions are expressible in terms of the venerable special functions of mathematical physics, namely (associated) Legendre and (associated) Laguerre polynomials. Since exactly the same predictions (about the energy levels) were already made by the “old”, ad hoc, Bohr-Sommerfeld, quantum mechanics, the “new” wave- and matrix- quantum theories needed to predict facts that were beyond the scope of the old theory, thereby offering a crucial confirmation. That’s why Schrödinger himself, Hylleraas, and many other physicists tried to derive the energy levels (alias eigenvalues) for two-electron atoms, whose Schrödinger equation, for the wave function ψ = ψ(x_{1}, y_{1}, z_{1}, x_{2}, y_{2}, z_{2}) is

(\frac{∂^{2}}{∂x_{1}^{2}} + \frac{∂^{2}}{∂y_{1}^{2}} + \frac{∂^{2}}{∂z_{1}^{2}} + \frac{∂^{2}}{∂x_{2}^{2}} + \frac{∂^{2}}{∂y_{2}^{2}} + \frac{∂^{2}}{∂z_{2}^{2}} + 2 (E + \frac{Z}{r_{1}} + \frac{Z}{r_{2}} − \frac{1}{r_{12}} )) ψ = 0,

Date: May 26, 2010.

∗ Supported by NFS-DMS 0070567 as a postdoctoral fellow, and by grant P20162 of the Austrian FWF.

∗∗ Supported in part by the United States of America National Science Foundation.

1available on-line from http://www.math.rutgers.edu/˜ zeilberg/akherim/pekeris1958.pdf (viewed June 12, 2016)

where E and Z are as above, while r_{1}, r_{2} are the distances of the electrons from the nucleus, and r_{12} is their mutual distance.

There were some crude attempts to use perturbation theory, but none of their predictions came close to the experimental spectra already known then. It was a major challenge to vindicate the new quantum mechanics by computation. It was one of the rare times in the history of science where the experimenters were ahead, and the theorists had to catch up.

Pekeris

Chaim Leib Pekeris (1908–1993) had a brilliant idea how to catch up. With a computer, of course! He had a carefully laid out approach, to be described soon, that would indeed give a very accurate prediction of the Helium spectra, with a powerful enough computer, and a clever enough programmer.

Except that when he first had that idea, computers didn’t yet exist, and when finally he had access to the JOHNNIAC, during his frequent long visits to the Institute for Advanced Study until von Neumann’s death (in 1957), it was not quite powerful enough, and at any rate too busy, to pursue Pekeris’ plan.

In addition to being a brilliant scientist, Pekeris was also an ardent Zionist. His good friend, (another Chaim, and another scientist), Chaim Weizmann (1874–1952), invited him, already in 1947, to head the department of applied mathematics at the Ziv Institute (that later was renamed the Weizmann Institute of Science), and Pekeris agreed, in principle, but only on condition that they build a computer similar to the JOHNNIAC. A committee was formed, including no lesser figures than Albert Einstein and John von Neumann, to decide whether such a plan was a good idea. Einstein believed not. In those days computers were very expensive, and he thought that such a poor, developing, country can make better use of such a big chunk of money, but von Neumann managed to win Einstein over and the plan was approved. It took a few years to materialize, and finally they recruited one of the members of von Neumann’s team, a visionary electrical engineer by the name of Gerald Estrin (b. 1921) [5]. Estrin narrates ([5], p. 319) that on one short conversation with von Neumann, shortly before his departure, he asked “What will that tiny country do with an electronic computer?”. John von Neumann responded: “Don’t worry about that problem. If nobody else uses the computer, Pekeris will use it full time!”. Estrin goes on to comment that this turned out to be an important prophecy that he often recalled.

C. L. Pekeris^{2}  WEIZAC^{3}

2 3

2Photo courtesy of Rutz, St. Moritz, Switzerland

3Photo courtesy of Yuval Madar under the Creative Commons Attribution ShareAlike 3.0 License, taken from http://en.wikipedia.org/wiki/WEIZAC

Pekeris’ Crazy Plan

The first step was standard. It turns out that the ground state described by ψ only depends on r_{1}, r_{2}, r_{12}, so one “merely” has to deal with functions of three variables, rather than six. The new partial differential equation, phrased in terms of r_{1}, r_{2}, r_{12} is easily derived (Eq. (5) of [9]). The next step (first suggested by H.M. James and A.S. Coolidge, see ref. 4 of [9]) was to make another change of variables, this time a linear one. After substituting E = −ε^{2}, introduce perimetric coordinates:

u = ε(r_{2} + r_{12} − r_{1}),
v = ε(r_{1} + r_{12} − r_{2}),
w = 2ε(r_{1} + r_{2} − r_{12}).

These new variables have the advantage that they range freely and independently from 0 to ∞. In contrast, r_{1}, r_{2}, and r_{12} are the lengths of the sides of a triangle (whose vertices are the two electrons and the nucleus), that must obey the triangle inequalities. In addition, the expected asymptotic behavior of ψ, deduced from the Hydrogen (one-electron) case, suggested to write ([9], Eq. (13))

ψ = e^{− \frac{1}{2} (u+v+w)}F (u, v, w),

and make F (u, v, w) the looked-for function. Pekeris performed this change of variables purely by hand and derived a fairly hairy linear partial differential equation with polynomial coefficients, satisfied by F, that we do not reproduce here, but that the curious reader can either look up ([9], Eq. (14)), or look at the computer output that is available from our webpages. The next step was to express F (u, v, w) as a series expansion of products of (simple) Laguerre polynomials (Eq. (16) of [9]):

F = ∑_{l,m,n=0}^{∞} A(l, m, n)L_{l}(u)L_{m}(v)L_{n}(w),

where L_{n}(x) denotes the Laguerre polynomial

L_{n}(x) = ∑_{k=0}^{n} \binom{n}{k}\frac{(−x)^{k}}{k!} .

Like all families of classical orthogonal polynomials, the Laguerre polynomials satisfy a pure (linear) differential equation, a pure (linear) recurrence equation, and a mixed differential-recurrence relation:

xL′′_{n}(x) = (x − 1)L′_{n}(x) − nL_{n}(x),
xL_{n}(x) = −(n + 1)L_{n+1}(x) + (2n + 1)L_{n}(x) − nL_{n−1}(x),
xL′_{n}(x) = nL_{n}(x) − nL_{n−1}(x),

the primes denoting differentiation with respect to x. Now came an astounding feat! Pekeris substituted the expansion for F (u, v, w), in terms of the yetto-be-determined A(l, m, n), into the above-mentioned linear differential equation (Eq. (14) of [9], politely not shown here), and successively used the above relations for the Laguerre polynomials to first get rid of all differentiations, and then, by using the pure recurrence, of any monomials in u, v, w. Then he collected terms, and got, purely by hand, a huge monster, a 33-term linear partial recurrence equation with polynomial coefficients satisfied by the A(l, m, n). Each of the coefficients of the 33 shifs A(l +α, m+β, n+γ) that showed up was polynomial in l, m, n of degree 3, and of degree 1 in in the charge Z and the yet-to-be found ε. We will kindly spare the reader this recurrence (and more important, us from typing it!), but the curious reader can glance at Eq. (22) of [9] and get horrified. We shudder to think of the great pain that the poor typist who keyed this from Pekeris’ hand-written manuscript, and the type-setter for Physical Review, not to mention Pekeris himself, had to go through. They all deserve lots of credit! In his wonderful essay [5] (p. 331), Estrin understatedly comments that the “appearance of this ugly 33-term recurrence would be enough to discourage most analysts”. The recurrence yielded a homogeneous linear system of equations with ∞^{3} equations and ∞^{3} unknowns, that usually has no non-trivial solutions, but for those ε for which the “determinant” vanishes (i.e., the “eigenvalues”), there are solutions. The largest eigenvalue is of physical relevance since it corresponds to the ground state energy of the atom. But even the most powerful computers can only handle finite systems! Hence the next step consisted in focusing on a finite, truncated version of the system by considering only those l, m, n ≥ 0 for which l + m + n ≤ ω, for some finite ω and set all the A(l, m, n) with l + m + n > ω equal to 0. In addition, the system could be cut approximately in half by using either the symmetry A(l, m, n) = A(m, l, n) in the so-called para states, or the antisymmetry A(l, m, n) = −A(m, l, n) in the so-called ortho states. But in order to make it into a matrix that could be handled on a (at that time, future) computer, one needed a convenient way to order linearly all the triplets of integers (l, m, n) with l+m+n ≤ ω and l ≤ m in the symmetrical case (resp. l < m in the antisymmetrical case). For this Pekeris devised a fairly complicated bijective map k: {(l, m, n) ∈ ℕ_{0}^{3} | l ≤ m} → ℕ that once again we spare the reader, but could be found in Eqs. (27-29) of [9] (by the way, Eq. (28) contains a very rare misprint, there should be \frac{1}{2}(l + m) added to it). It is not known when Pekeris devised this plan, but it is probably several years before he had access to a computer, so he waited eagerly until Chaim Weizmann’s promised computer would materialize, and the recommendation of the above-mentioned committee of Einstein, von Neumann et. al. would be carried out.

WEIZAC

As we have already mentioned above, the person chosen to head the team that would build, from scratch, the first Israeli electronic computer was Estrin, and his vivid account [5] is hereby highly recommended for anyone interested in the history of computing. The WEIZAC team consisted of a cadre of young and talented electrical engineers (including Aviezri Fraenkel (b. 1929) who later went on to do a Ph.D. in number theory, and became, inter alia, an authority on combinatorial games and pioneered the use of computers in religious studies). Finally the computer was ready, and Pekeris was itching to use it on his many problems, including the spectra of Helium, but he needed a programmer (what today we would call a “software engineer”, except there was no such thing as software in those days): of course, neither a Java programmer, nor a Fortran programmer, and not even an Assembly-language programmer. Back in 1957 these were yet to be invented. The only language that WEIZAC understood then was machine language and the alphabet consisted of two letters only, 0 and 1 (via the 16-letter alphabet of hexadecimals). But how to find such programmers? Definitely not among graduates of computer science departments, that didn’t yet exist. What Pekeris did was to ask his secretary to place classified ads in the daily newspapers, asking for high school graduates, after their military service, who attended the megama re’alit (math/science track).

Accad

Yigal Accad (b. 1936), fresh out of his military service, answered such an ad. In a recent e-mail message, dated May 7, 2010, Accad recalls:

“On a 1957 Friday (or was it a Holiday Eve) that happened to be a non-working day at the Weizmann Institute, Prof. Pekeris unexpectedly drove his 1948 Studebaker to our residence at the southern edge of Rehovot. He invited me to join him in his office. Over there he pulled out a pile of handwritten papers and went with me through many of the equations you can find in the 1958 paper, including Eq. (22). As I remember, this tour took at least 2 hours. At the end Prof. Pekeris asked me if I can handle this problem. There were only 2 possible answers to this question and the rest is history. This may have been the best risk I have taken.”

Estrin goes on to state the following accolades ([5], p. 330):

“There is a clear testimony to the fact that Yigal Accad had unusual ability to use WEIZAC as a tool with very little software between him and the machine semantics. That ability, when combined with his talents as an applied mathematician, was a significant factor in the ensuing problem-solving successes at the Weizmann Institute.”

Accad became Pekeris’ right-hand man for many years, and it is hard to imagine what Pekeris would have done without him. Pekeris appreciated Accad’s invaluable work, and suggested that parallel to working full-time as a software engineer, Yigal would enroll in the graduate school (after completing his undergraduate studies at the Hebrew University) and incorporate some of the research into, first a masters’ thesis, in 1969, and a Ph.D. in 1973 (which was a far reaching extension of the work narrated here). Accad stayed at the Weizmann Institute from 1956 until 1989. Between 1977 and 1989 he also served as a consultant to the pioneering Israeli Hi-Tech company Scitex. In 1989 he moved to California and joined Electronics for Imaging (EFI), working there until 2008, ultimately becoming chief scientist.

The Pekeris-Accad-WEIZAC collaboration

And indeed Accad was the perfect person to tame Pekeris’ monster recurrence and to write (machine-language!) programs to generate the truncated matrices, and to implement the iterative algorithm for estimating the largest eigenvalue. The impressive (for its time) WEIZAC output is displayed in Table III of [9] for values of the charge Z ranging from Z = 1 to Z = 10. We are happy to report that our 2010 computations (on three different platforms) completely agree with that table, all the way to the last decimal digit. In a follow-up paper, published a year later, Pekeris [10] (and of course, Accad and WEIZAC—but it took more than 30 years, before a computer, Shalosh B. Ekhad, became co-author!) treat the important special case of Helium (Z = 2) with a greater accuracy, and also consider the ortho state 2 ^{3}S. Our computations agree with that paper, too.

2010

Of course, thanks to Moore’s Law, all these computations can be made so much faster today, and there is no reason for us to be proud of the fact that we can compute the eigenvalues within seconds with nowaday’s hardware and software, a task that kept WEIZAC busy round-the-clock for months: for example, a fixed-point multiplication took 1 millisecond on this early computer and the capacity of its memory was 4096 words (40 bits per word). But what is still remarkable and probably not so obvious to many people: not only the WEIZAC part, the numeric computation that can now be done on every laptop, and the Accad part, the challenging machine-language programming that today becomes an easy exercise with high-level programming languages, but also, and especially, the Pekeris part can be done much faster and mostly automatically, using computer algebra. Even more: in view of the gigabyte-sized recurrences that we can currently deal with (and in fact we do! see for example [7]) with symbolic software, the “monster recurrence” looks rather dwarfish. We don’t know exactly how long it took Pekeris to derive the differential equation and the recurrence, but let’s say 20 person-hours (including checking and rechecking); our program needs 0.108 seconds. To be honest, it took us a couple of hours to program Maple and Mathematica to follow Pekeris’ plan, but with almost the same effort, one could (and we did) program the general problem, that could be used again and again for many other differential equations in future problems. Our programs PEKERIS (for Maple, by DZ) and Pekeris.nb (for Mathematica, by CK) are indeed very general: they basically can input any linear differential equation, in any number of variables, and any series of substitutions, and output the transformed differential operator. Also the recurrence for a Laguerre polynomial expansion is achieved completely automatically. Using the widely-used concept of Gröbner bases (that was invented by Bruno Buchberger in 1965 and hence not yet available for Pekeris!) it is also possible to perform the series expansion for any set of orthogonal polynomials of hypergeometric type. For this purpose, the defining equations for the family of polynomials are represented as a Gröbner basis, which makes sense when the relations are rewritten, in operator notation, as (noncommutative) polynomials. Having chosen an appropriate monomial order, the elimination of the differentials can be achieved by a simple reduction modulo the Gröbner basis. Similarly, by changing the underlying polynomial ring, the elimination of the continuous variables u, v, w can be done. Let us also remark that you don’t need to be a Laguerre or a Pekeris to generate the relations for the Laguerre (and other orthogonal) polynomials. They are all routinely derivable (and provable) by the so-called Wilf-Zeilberger method [12], as implemented, e.g., in the Mathematica package HolonomicFunctions [6] that we employ in our program.

Modular techniques using Chinese remaindering and polynomial interpolation allow for computing the determinant symbolically up to quite large dimensions: for example, the determinant of the 161×161 matrix (para case with ω = 10) is obtained in less than five minutes, yielding a polynomial in ε of degree 161 having integer coefficients with about 500 digits! It is clear that this strategy produces a lot of overhead, so that an alternative way is desirable. We reformulate the problem of finding the largest ε for which the determinant of M ∈ ℤ[ε]^{n×n} vanishes, as a generalized eigenvalue problem:

Av = εBv, M = A − εB with A, B ∈ ℤ^{n×n}.

Although Maple and Mathematica are symbolic programs in the first place, they offer quite some functionality for numerical evaluation, in particular for the above problem. But since we were not 100% satisfied with both—Maple was rather slow for the desired precision and Mathematica didn’t allow higher precision than machine reals (6 decimal digits)—we tried with MATLAB, a software designated for numeric computations, especially in linear algebra. Notably, the program code for building the (sparse) matrices is itself computer-generated! It contains the 33 terms of the recurrence hard-coded to produce the matrix entries, and therefore certainly comes closer to Accad’s machine-code program. We ourselves were very much impressed by MATLAB’s speed and accuracy. Computing all entries of Table III of [9] takes less than a second, and without much effort ω can be increased to 60, corresponding to a 20336 × 20336 matrix.

Software and Sample Output

This article is accompanied by the Maple package PEKERIS, available from

    http://www.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/pekeris.html ,

where the readers can also find lots of output files (and input files if they want to modify them to get more output) that reproduce and far extend the seminal 1958 computations of Pekeris, Accad, and WEIZAC. Further we provide the Mathematica notebook Pekeris.nb (for which the package HolonomicFunctions is required), and the MATLAB programs PekerisPara.m and PekerisOrtho.m, all available from

               http://www.risc.jku.at/people/ckoutsch/pekeris/ .

Our maplephone readers are welcome to play with the first package while the mathematicaphones would probably prefer the latter one. However, even people (shame on you!) who speak neither Maple nor Mathematica can appreciate the output files, written in plain humaneze. The secondnamed author is particularly proud of the procedure PaperPara that fully automatically and seamlessly generates a whole article, ready to be submitted to Physical Review, without any human touch. Changing the parameters can produce many similar papers, see

           http://www.math.rutgers.edu/~zeilberg/tokhniot/oPEKERIS1 .

Conclusion

This article is first and foremost an ode to the vision and ingenuity of computing pioneers, but it also makes the point that there are lots of hidden treasures in the “old” scientific literature, that can be revisited with today’s powerful symbolic computation software. We are not the first to advocate using symbolic computations in scientific computing, see for example [3] (unfortunately he was unaware of [13]), and the current impressive application to high-energy physics [4], but we believe that there is a huge potential for exploiting symbolic computation on problems that previously seemed intractable. This would complement the extensive use (and according to Nobelist Philip Anderson, excessive, and sometimes abusive [1]) of Monte Carlo methods. In particular, the Wilf-Zeilberger algorithmic proof theory [12] (and more importantly the subsequent generalizations to multi-summation and multi-integration [13, 2]), should be taught to all scientists. We would be more than happy if this article is the seed for future collaborations between symbolic computation and physics, chemistry or other sciences.

Encore

Many people, even today, are not comfortable with computer-generated or even computer-assisted proofs, like the four color theorem or the Kepler conjecture, since they are uncomfortable trusting the computer. While the “monster recurrence” discussed above was still derived purely by hand, Pekeris must have started using his own “symbolic” computation when he tackled seemingly intractable problems. Let us end with his prophetic words ([11], quoted in [5], p. 333)

“Here we are confronted with problems where the computer writes the formulae as well as evaluates them. By the nature of their origin such formulae are very long— in many cases too long to be published. We shall therefore be dealing in the future with equations which only the computer will see. The prospect of operating with invisible equations is a frightening one, but the alternative is to accept the situation of the past, where problems have been staring at the applied mathematician for decades, and even more for centuries, without a practical solution being reached. A problem, like the tides of the oceans, for example, is not necessarily insoluble just because it had remained in the books for 184 years.”

References

[1] Philip W. Anderson. What is wrong with QMC?, talk delivered at the 103rd Statistical Mechanics Conference,

    May 10, 2010, Rutgers University, http://www.math.rutgers.edu/events/smm/smm103-invitedtalks.html.

[2] Moa Apagodu and Doron Zeilberger. Multi-Variable Zeilberger and Almkvist-Zeilberger Algorithms and the Sharpening of Wilf-Zeilberger Theory, Adv. Appl. Math. 37 (Special Regev issue), pp. 139–152, 2006.

[3] Michael P. Barnett. Symbolic Computation of Integrals by Recurrence, ACS SIGSAM Bulletin, 37(2), pp. 49–63, 2003.

[4] I. Bierenbaum, J. Blümlein, S. Klein, C. Schneider. Two-Loop Massive Operator Matrix Elements for Unpolarized Heavy Flavor Production to O(ϵ), Nucl.Phys. B803(1-2), pp. 1–41, 2008.

[5] Gerald Estrin. The WEIZAC years, IEEE Trans. for the Annals of the History of Computing 13(4), pp. 317–339, 1991.

[6] Christoph Koutschan. Advanced Applications of the Holonomic Systems Approach, Ph.D. thesis, RISC, Johannes Kepler University, Linz, Austria, 2009.

[7] Christoph Koutschan, Manuel Kauers, and Doron Zeilberger. A proof of George Andrews’ and David Robbins’ q-TSPP conjecture, Technical Report arXiv:1002.4384, 2010.

[8] Linus Pauling and E. Bright Wilson, Jr. “Introduction to Quantum Mechanics with Applications to Chemistry”, MacGraw-Hill, New York 1935 (reprinted by Dover, 1985).

[9] C. L. Pekeris. Ground State of Two-Electron Atoms, Phys. Rev. 112(5) pp. 1649-1658, 1958.

[10] C. L. Pekeris. 1 ^{1}S and 2 ^{3}S States of Helium, Phys. Rev. 115(5), pp. 1216-1221, 1959.

[11] C. L. Pekeris. Propagation of seismic pulses in layered liquids and solids, In Norman David (ed.) “International symposium on stress wave propagation in materials”, New York, Interscience Publishers, Inc., 1960.

[12] Marko Petkovšek, Herbert S. Wilf, Doron Zeilberger. “A=B”, A.K. Peters, 1996, available on-line.

[13] Herbert S. Wilf, Doron Zeilberger. An algorithmic proof theory for hypergeometric (ordinary and “q”) multisum/integral identities, Invent. Math. 108, pp. 575-633, 1992.

Christoph Koutschan, Research Institute for Symbolic Computation, Johannes Kepler University Linz, Austria and (until Aug. 2010) Tulane University

E-mail address: ckoutsch@risc.uni-linz.ac.at

Doron Zeilberger, Mathematics Department, Rutgers University (New Brunswick), Piscataway, NJ, USA.

E-mail address: zeilberg@math.rutgers.edu