Two-dimensional calculation of sub-barrier heavy-ion fusion cross sections

Two-dimensional calculation of sub-barrier heavy-ion fusion cross sections

Nuclear Physics @ North-Holland A368 (1981) 352-364 Publishing Company TWO-DIMENSIONAL CALCULATION FUSION CROSS OF SUB-BARRIER SECTIONS HEAVY-ION...

795KB Sizes 0 Downloads 16 Views

Nuclear Physics @ North-Holland

A368 (1981) 352-364 Publishing Company

TWO-DIMENSIONAL

CALCULATION FUSION CROSS

OF SUB-BARRIER SECTIONS

HEAVY-ION

S. LANDOWNE Sektion Physik, Universitiil Miinchen. D-8046 Garching. Federal Republic oJ’German)

J. R. NIX + Sektion Physik, Universitci’tMtinchen, D-8046 Garching, Federal Republic of German) and Theoretical Division, Los Alamos National Laboratory, Los Alamos, New Mexico 87545 Received

21 April

1981

Abstract: Using a two-dimensional potential-energy surface and a diagonal kinetic energy we calculate the sub-barrier heavy-ion fusion cross section for the reaction s*Ni + ‘sNi + ’ 16Ba. We represent the effective potential energy by an expression based on separated nuclei that includes monopole and quadrupole Coulomb interaction energies. a Yukawa-plus-exponential nuclear interaction energy, a parabolic deformation energy and a centrifugal energy. To approximate the penetrability through the resulting two-dimensional barrier surface in the separation and deformation coordinates, we take into account the spheroidal zero-point motion and average dynamical deformation of the colliding nuclei by using a gaussian superposition of straight-line paths centered about a nuclear deformation which is determined by solving classical equations of motion up to the turning point. Relative to the result calculated for a one-dimensional barrier corresponding to spherical nuclei, the sub-barrier fusion cross section is-enhanced by the zeropoint motion but is suppressed by the dynamically induced oblate deformations. The combined effect for 58Ni + ‘sNi gives only a slight enhancement which does not explain the experimental measurements.

1. Introduction

Recent experimental measurements of heavy-ion fusion cross sections at subbarrier energies deviate significantly from calculations based on penetrating onedimensional fusion barriers ’ - 4). In particular, it has been found by Beckerman et af. ‘) for reactions involving Ni projectiles and targets and by Evers et al. ‘) for reactions involving S and Ru, that the experimental fusion cross sections are about lo3 times the calculated values at center-of-mass bombarding energies 5 MeV below the Coulomb barrier, although they agree approximately at higher energies near the barrier top. Similar discrepancies are also present in the earlier + Alexander

von Humboldt

Senior US Scientist

Awardee 352

S. Landowne,

J. R. Nix / Two-dimensional

calculation

353

experiments of Stokstad et al. 3, involving Ar and Sm and, to a lesser extent, in those of Sikora et al. “) for Ca and Ni, after the effects of static target deformations have been taken into account. A possible explanation of these deviations is that the nuclei are not rigid during their approach but instead undergo dynamical deformations that, loosely speaking, consist of quanta1 zero-point oscillations about average dynamical trajectories. Kodama et al. ‘) and Esbensen 6, have shown that quanta1 zero-point oscillations enhance the sub-barrier fusion cross section. The effect of the average dynamical motion ’ - ’ “) is more complicated because the oblate deformations induced by the long-range Coulomb force suppress the sub-barrier fusion cross section, whereas the prolate defo~ations induced by the short-range nuclear force enhance it. Do the combined effects of quanta1 zero-point oscillations and average dynamical deformations explain the enhancement of the experimental sub-barrier fusion cross sections? To study this question, we incorporate these two effects into a calculation of the sub-barrier fusion cross section for the reaction 58Ni+ 58Ni + ’ 16Ba measured by Beckerman et al. ‘). The two-dimensional potential-energy surface and diagonal kinetic energy that we use for this purpose are discussed in sect. 2. The approximation that we use for calculating the penetrability through the twodimensional barrier is discussed in sect. 3. Calculations of the sub-barrier fusion cross section are compared with experimental results in sect. 4. Our conclusions are summarized in sect. 5.

2. Two-dimensional potential-energy surface To reduce the number of degrees of freedom that must be considered, we specialize to collisions with symmetric targets and projectiles which have spherical ground-state shapes and we describe the separated nuclei in terms of two collinear spheroids I’). This leads to two collective coordinates, the distance r between the centers of mass of the two nuclei and the semi-symmetry axis c of one of the spheroids. Since all the quantities that enter our calculation of the sub-barrier fusion cross section correspond to separated nuclei, it is not necessary to specify the precise shapes or energies of the nuclei after they come into contact. For a given angular momentum I our expression for the effective potential energy Vf(r, c) is F(r, c) = Zfe2/r+gZie2RI(c-RR,)/r3+

V,(r,c)f+K(c-R,)*+*l(l+

l)V/(~r2),

(1)

where Z,e is the electronic charge of one of the nuclei, R, is the equivalent sharp-surface radius of one of the spherical nuclei, K is the stiffness against spheroidal deformations of the combined two-nucleus system and ~_lis the reduced mass. The five terms in eq. (1) represent, respectively, the monopole Coulomb interaction

354

S. Landowne, .I. R. Nix / Two-dimensional

calculation

energy, the quadrupole Coulomb interaction energy, the nuclear energy, the deformation energy and the centrifugal energy. The radius R, is calculated from

interaction

R, = r,A$,

where A, is the mass number of one of the nuclei and the equivalent sharp-surface nuclear-radius constant r0 has the value 16) r0 = 1.16fm. We calculate the nuclear interaction energy by means of an approximation to the Yukawa-plus-exponential potential r6* l’) that is exact for two spherical nuclei. The deformation dependence is incorporated in terms of the distance between the inner surfaces of the nuclei and the radius of curvature at the point of closest approach by use of the proximity theorem ‘*). The resulting nuclear interaction energy of two separated symmetric spheroidal nuclei is I/,(r,c) =

-2D[F+ r$]$exp[-(r-Zc)/a],

where a is the range of the Yukawa-plus-exponential D is given by “) D = 2dIs(Rl/~)12 exp(-2R&b,{l

potential. The depth constant

-Q(N1 -~,YA,12}/(r~R1),

where a, is the surface-energy constant, K, is the surface-asymmetry is the neutron number of one of the nuclei and

constant, N,

g(x) = x cash (x) - sinh (x). The constant F is given by l’) F

=

4+2R,la-2

f(Rl14 g(Rl,a),

where f(x)

= x2 sinh(x).

For the values of the constants appearing in these expressions we use 16) a = 0.68 fm,

a, = 21.13 MeV,

K,

=

2.3.

As an indication of the accuracy of our approximation for deformed nuclei, we mention that for two touching symmetric 58Ni spheroidal nuclei with semisymmetry axis c = l.2Rl, our approximation yields a nuclear interaction energy of - 16.4 MeV, which is 9 ‘A smaller in absolute value than the exact Yukawaplus-exponential interaction energy calculated by numerical quadrature ’ 9). The stiffness K against spheroidal deformations of the combined two-nucleus system is related to the stiffness K2 against spherical-harmonic e distortions of a

S. Landowe,

J. R. Nix { Two-dimensional calculation

355

single nucleus by K = 2(4n)r$,/R:. For K2 we use the microscopic value extracted by Wong ‘O) from the experimental energy and half-life of the first excited 2+ level in **Ni, which yields K = 25.9 MeV/fm’.

(2)

For separated nuclei and zero angular momentum our expression (1) reproduces the dominant features of two-dimensional potential-energy surfaces calculated within a macroscopic mode1 by numerical quadrature 21). In particular, as the nuclei approach each other from infinity, the minimum potential energy for a given distance between mass centers initially occurs for oblate deformations because of the quadrupole Coulomb interaction energy. At closer distances this minimum shifts to prolate deformations because of the nuclear interaction energy. The kinetic energy of the system is given by T = +pi2 + +Bk=,

(3)

where a dot denotes differentiation with-respect to time and B is the inertia against spheroidal deformations of the combined two-nucleus system. This is related to the inertia B, against spherical-harmonic Yi distortions of a single nucleus by B = 2($n)B,/R; For B, we again use the microscopic value extracted by Wong “) for **Nil which yields B = 4.76 x 10’ MeV/c&

(4)

where co is the speed of light,

3. Two-dimensional

barrier penetration

Our calculation of the sub-barrier fusion cross section involves the penetrability through the two-dimensional potential-energy surface (1) with the kinetic energy (3). A complete calculation would require solving a two-dimensional Schriidinger equation. Because of the dif~culty of this problem, several approximations have been used over the past few years for calculating the penetration through multidimensional barriers of special types ** 6* 2* - =‘). For the problem of spontaneous-fission half-lives Brack et al. ‘=) and Ledergerber and Pauli 23*24) assumed that the nucleus takes a path from its ground state through the multi-dime~siona1 fission barrier that maximizes the one-dimensionalWKB expression for the penetrability. However, this technique would give an

356

S. Landowne, J. R. Nix / Two-dimensional calculation

anomalously large penetrability for two-dimensional fusion barriers which have a low pass situated on the side of the valley through which the incident and reflected waves travel. Also, it neglects the effect of changes in the frequency of the motion in the transverse direction at different distances through the barrier. Hofmann 25) has applied the WKB approximation and Born approximation to a special class of two-dimensional fission barriers. For a potential of the form I+,

.Y) = V,(x) + 4&.(X)f2,

t5a)

and kinetic energy of the form T = fM,(~)i~

+:M,,(x)$~,

the two-dimensional penetrability is approximately calculated for the effective one-dimensional potential v,,,(X) = V,(x) + f@$(4

Ub)

equal to the penetrability

- e$( - m)l

@a)

and kinetic energy

(6b) where the frequency u,(x) is given by

In other words, the penetrability is calculated for an effective one-dimensional barrier along the fission direction that contains at every point the change in the quanta1 zero-point energy in the transverse direction relative to that at - co from where the incident wave is coming. Unfortunately, this approximation cannot be directly applied to sub-barrier heavy-ion fusion because the fusion valley is not straight, with the potential-energy surface changing into a ridge near the point of contact. For the calculation of chemical reaction rates Miller 26) has developed a multi-dimensional barrier-~netration method that invloves solving classical equations of motion for complex time. This method has been used by Ring et al. 27) to calculate the penetrability through a simple two-dimensional fission barrier of the type (5) at energies well below the barrier top. .However, Miller’s method is difficult to compute with and also does not apply for incident energies near the top of the barrier. Kodama et QI. ‘) approximated the two-dimensional fusion barrier in terms of discontinuous portions of parabolic barriers, for which an exact solution of the two-dimensional Schriidinger equation is possible. However, the discontinuities give rise to spurious peaks in the penetrability. Furthermore, despite the unrealistic simplicity of the barrier, the calculation is very difficult to carry out numericaily. Esbensen 6, approximated the effects of zero-point motion on the ~netrabiIity through a multi-dimensional fusion barrier by taking a superposition of one-

S. Lmdowne, J. R. Nix / Two-dimensional calculation

357

dimensional penetrabilities calculated for an effective angular momentum that simulated the zero-point motion. However, Esbensen’s formulation did not account for the shift in the average dynamical path, which should be important for the lowenergy fusion of heavier systems and may also affect the details of the cross sections at higher energies. In our work, we approximate the two-dimensional penetrability in terms of a gaussian superposition of one-dimensional penetrabilities that includes both the zero-point motion in the transverse direction and the shift in the average dynamical path, which we determine by solving classical equations of motion up to the turning point. The penetrability for a given partial wave I is given in our approximation by p, =

--‘__

s SC8

(27ta2)f ()

exp [ -(c - d2/(2a2)] PI,&,

(7)

which is evaluated in practice by extending the lower integration limit to - CCand using a nine-point gaussian-Hermite quadrature formula 28). The variance (T’ of the zero-point motion in the transverse c direction is determined from the square of the ground-state harmonic-oscillator wave function to be o2 = fh/( KB)+. For a system consisting of two 58Ni nuclei, use of the microscopic values (2) and (4) for K and B leads to fs = 0.168 fm. The quantity c1 is the value of the deformation coordinate c at the classical turning point for the Ith partial wave and Pl,, is the penetrability for the fth partial wave calculated for a straight-line path with deformation c. This penetrability is given in the WKB approximation by 29) P I.

c =

l/Cl +ew(Kl, ,)I.

where 2(2/# r* jQc = --[v,(r, ~)-E,.,.-fK(c-R1)~]~dr, h s r1 with ECm being the center-of-mass bombarding energy and r, and r2 are the points of entrance into and emergence from the barrier at the energy &,,.+~K(c-&)~. When applied to the special two-dimensional barrier-penetration problem (5), our approximation (7) gives to lowest order the same result obtained by Hofmann 25) involving the effective one-dimensional system (6). This is seen by noting first of all that for a straight valley the gaussian appearing in eq. (7) is the square of the

358

S. Landowne, J. R. Nix / Two-dimensional

calculation

ground-state wave function in the transverse direction. Then, by first-order perturbation theory, the term Vl(r, c)-+K(cR,)’ appearing in eq. (8) yields to lowest order the effective potential (6a) containing the change in zero-point energy in the transverse direction. We are able to test the accuracy of our approximation for a particular twodimensional fission barrier of the type (5) which was studied by Ring et al. *‘). In this case the inertias A4, and A4,,are constant while V,(x) = V0exp (- x*/a*), K,(x)

= C[ 1 + a exp (- x*/a*)],

where V,, a, C and a are constants. Ring et al. 27) calculated the exact penetrability for this potential by numerically solving a two-dimensional Schrijdinger equation and compared the results to the maximum one-dimensional-WKB penetrability of Brack et al. *‘) and Ledergerber and Pauli 23*24) and to uniform semi-classical calculations based on Miller’s method 26). Our numerical result for the penetrability using the parameters of ref. *‘) is given in table 1 where we also include the result of using the effective one-dimensional potential (6) developed by Hofmann 25). It is seen from table 1 that the present method gives the best correspondence with the exact calculation. For our calculations of the fusion cross section in sect. 4 we make the further approximation that the potential near its maximum is parabolic, in which case eq. (8) becomes Ul,

E

=

‘W?, c-E,.,. -3m-~,)*l/(~q.

,),

where P,, c is the height of the potential for the partial wave 1and deformation

c and

TABLE 1

Comparison of five methods for calculating the penetrability through the two-dimensional fission barrier considered by Ring et al. “) Method exact quanta1 27) maximum WKB penetrability 2’) uniform semi-classical 27) Hofmann’s effective potential present approximation

Penetrability 1.40x 1.60x 1.44x 1.34x 1.37 x

10-5 10-S 1o-5 10-5 10-s

S. Landowne, J. R. Nix / Two-dimensional calculation with Yt, Ebeing

359

the position of the maximum in V,(r, c). This approximation simplifies the calculations, particularly for energies above the barrier when the integral in eq. (8) must be done in the complex r-plane. It should be noted that because of the long range of the Coulomb potential, the parabolic approximation overestimates the penetrability at energies below the barrier. To determine the value cI of the deformation coordinate at the classical turning point for the Eth partial wave, which we require in eq. (7), we integrate Hamilton’s equations of motion using a modilied predictor-corrector method due to Hamming 30). Thus with the kinetic energy (3) and potential energy (1) we form the hamiltonian I?= T+V. Then introducing the conjugate momenta P, = @ P, = Be, Hamilton’s equations give

E=

5 =p,/B, c

2DR:

- 7 . PC=--

i?H

{[F+(r-2c)/a](l

+a/r)-

1) exp[-(r-2c)/a],

= -gZfe2R,/r3+K(c-R,)

4DR; +p arc2

{[F+(r-2c)/a](l-a/c)-l}exp[-(r-2c)/u].

We initialize these equations by the conditions r = lOR,, c = R,, P, = - 12m,

.m. -Z:e’/r-1(Z+

1)fi2/(2p+2)])+,

PC = 0,

and integrate them until the radial momentum p, vanishes.

360

S. Landoune, J. R. Nix / Two-dimensional calculation 4.

Fusion cross section

The fusion cross section err is calculated as Crf= ;

g

(21+ l)P,,

I-O

where k is the wave number at infinity. To illustrate the separate contributions from the two effects that we consider, we calculate the penetrability P, first by including only the zero-point motion and second by including only the average dynamical deformation, before combining both effects together. These results are presented successively in figs. l-3 for the reaction s8Ni+ “Ni + ’ “jBa, in comparison with the experimental data of Beckerman et al. ‘). We also include in each figure the result calculated for a one-dimensional fusion barrier corresponding to spherical nuclei, for which the barrier height for zero angular momentum is 102.1 MeV. Fig. 1 shows the effect of the zero-point motion alone, which is calculated by setting cI = R, in eq. (7) for all partial waves. At high bombarding energies zero-point motion has little effect on the fusion cross section. At sub-barrier

Center-of-Moss

Energy Gm.

(MeV)

Fig. 1. Effect of target and projectile spheroidal zero-point motion on the fusion cross section for 58Ni+58Ni -+ ‘r6Ba. The solid circles give the experimental data of Beckerman et al. I), the dashed curve gives the result for a one-dimensional barrier corresponding to spherical nuclei, and the solid curve gives the result for a gaussian ,superposition of straight-line paths centered about spherical nuclei with deformations determined by spheroidal zero-point motion.

S. Landowe,

J. R. Nix 1 Two-dimensional

58Ni + WJi

calculation

361

* f’ L@ I

.



/1’

I

I i

l

/I- dlmenstonal I.

Center-of-Mass

I

Energy

E,,,fMeVl

Fig. 2. Effect of target and projectile average dynamical deformation on the fusion cross section for 58Ni+58Ni + ‘r6Ba. The solid circles give the experimental data of Beckerman et al. ‘), the dashed curve gives the result for a one-dimensional barrier corresponding to spherical nuclei, and the solid curve gives the result for a single straight-line path with nuclear deformation determined by solving classical equations of motion up to the turning point.

energies zero-point motion enhances the fusion cross section by a factor of about 50, but this still lies below the experimental data by a comparable amount. Furthermore, the slope of the curve including zero-point motion is somewhat smaller at low energy than the experimental slope. Fig. 2 shows the effect of the average dynamical deformation alone, which is calculated by going to the limit B’ = 0 in eq. (7) so that only the penetrability P,, cI contributes for a given partial wave. This result is calculated only for sub-barrier energies because at higher energies the dynamical trajectory leads to overlapping nuclei, for which we have not specified the shapes and energies. Because of the oblate deformation induced dynamically by the Coulomb interaction energy, the barrier height for each partial wave is increased slightly relative to that for spherical nuclei. This suppresses the sub-barrier fusion cross section by a factor of about 15. Fig. 3 shows the combined effect of both the zero-point motion and the average dynamical deformation, which is calculated by use of our full approximation (7) to the penetrability through a two-dimensional potential-energy surface. The net result is very ‘similar to that calculated for a one-dimensional barrier corresponding

362

S. Lmdowne, J. R. Nix 1 Two-dimensional calculation

58Ni+58Ni

a / 2*

M-L

0.01 90

95

Center-of-Mass

100 Energy

/

5,H’

*



105

110

~=.~~MeVl

Fig. 3. Combined effect of both target and projectile zero-point motion and average dynamical deformation on the fusion cross section for 58Ni + 58Ni + t16Ba. The solid circles give the experimental data of Beckerman ef al. ‘), the dashed curve gives the result for a one-dimensional barrier corresponding to spherical nuclei, and the solid curve gives our final result for a two-dimensional potential-energy surface.

to spherical nuclei! However, there is an enhancement by a factor of about 3 at low bombarding energies and the slope of the two-dimensional curve is slightly smaller than that of the one-dimensional curve. At low energies the experimental data are still about lo3 times our final result.

5. Conclusions We have calculated the sub-barrier fusion cross section for the reaction 58Ni+ s8Ni + l16Ba by taking into account the zero-point motion and average, dynamical deformation for a two-dimensional potential-energy surface and diagonal kinetic energy involving the distance between mass centers and spheroidal deformation of the colliding nuclei. No adjustable parameters were introduced, but instead all of the quantities that enter the calculation were taken from previous work. We found that relative to the result calculated for a one-dimensional barrier corres~nding to spherical nuclei, the sub-barrier fusion cross section is enhanced by the zero-point motion but is suppressed by the dynamically induced oblate

S. Landowe,

J. R. Nix / Tuso-dimensional calculation

363

deformations. Although the combined effect is a slight enhancement, the experimental data at low energies are still approximately lo3 times our tinal result. Therefore, zero-point motion and dynamical deformations of the type considered here are not the explanation for the enhanced sub-barrier fusion cross sections that have been observed experimentally. This leads us to conclude that the most likely alternative explanation is that the attempted exchange of neutrons near the classical turning point leads to the formation of a neck and subsequent fusion through a sequence of more complicated shapes, a possibility that has also been suggested in refs. ‘. 2). We are grateful to D. Evers. H. Hofmann, and H. H. Wolter for stimulating discussions Yukawa-plus-exponential interaction energy nuclei.

H. J. Krappe, P. Ring, S. J. Skorka and to A. J. Sierk for calculating the of two touching prolate spheroidal

This work was supported by the Federal Republic of Germany Bundesministerium fur Forschung und Technologie, the Alexander von Humboldt Foundation and the US Department of Energy.

References 1) M. Beckerman, 2) 3) 4) 5) 6) 7) 8) 9) 10) I I) 12) 13) 14) IS) 16) 17) 18) 19) 20) 21) 22)

M. Salomaa, A. Sperduto, H. Enge, J. Ball, A. DiRienzo, S. Gazes, Yan Chen, J. D. Molitoris and Mao Nai-feng, Phys. Rev. Lett. 45 (1980) 1472 D. Evers. K. Rudolph. 1. Weidl, P. Konrad, U. Quade, K. E. G. Lobner and S. J. Skorka, Proc. 19th Int. Winter Meeting on nuclear physics, Bormio, Italy, 1981 R. G. Stokstad, W. Reisdorf, K. D. Hildenbrand, J. V. Kratz, G. Wirth, R. Lucas and J. Poitou, Z. Phys. A295 (1980) 269 B. Sikora. J. Bisplinghoff, W. Scobel, M. Beckerman and M. Blann, Phys. Rev. C20 (1979) 2219 T. Kodama. R. A. M. S. Nazareth, P. Mbller and J. R. Nix, Phys. Rev. Cl7 (1978) 111 H. Esbensen, Nucl. Phys. A352 (1981) 147 J. Maly and J. R. Nix, Contributions, Int. Conf. on nuclear structure, Tokyo, Japan, 1967, University of Tokyo report (1967) p. 224 H. Holm, W. Scheid and W. Greiner, Phys. Lett. 29B (1969) 473 H. Holm, D. Scharnweber, W. Scheid and W. Greiner, Z. Phys. 231 (1970) 450 H. Holm and W. Greiner, Phys. Rev. Lett. 24 (1970) 404 A. S. Jensen and C. Y. Wong, Phys. Rev. Cl (1970) 1321 A. S. Jensen and C. Y. Wong, Nucl. Phys. A171 (1971) 1 P. W. Riesenfeldt and T. D. Thomas, Phys. Rev. C2 (1970) 71 I P. W. Riesenfeldt and T. D. Thomas, Phys. Rev. C2 (1970) 2448 J. R. Nix and W. J. Swiatecki, Nucl. Phys. 71 (1965) 1 P. Moller and J. R. Nix, Nucl. Phys. A361 (1981) I I7 H. J. Krappe, J. R. Nix and A. J. Sierk, Phys. Rev. C20 (1979) 992 J. Blocki, J. Randrup, W. J. Swiatecki and C. F. Tsang, Ann. of Phys. 105 (1977) 427 A. J. Sierk, private communication (1981) C. Y. Wong, Nucl. Data A4 (1968) 271 P. Moller and J. R. Nix, Nucl. Phys. A272 (1976) 502 M. Brack, J. Damgaard, A. S. Jensen, H. C. Pauli, V. M. Strutinsky and C. Y. Wong, Rev. Mod. Phys. 44 (1972) 320

364

S. Lundowne, J. R. Nix / Two-dimensional

23) T. Ledergerber and H. C. Pauli, Nucl. Phys. A207 (1973)

calculation

1

24) H. C. Pauli and T. Ledergerber, Proc. 3rd IAEA Symp. on the physics and chemistry of fission, Rochester, 1973. vol. I (IAEA, Vienna, 1974) p. 463 25) H. Hofmann, Nucl. Phys. A224 (1974) 116 26) W. H. Miller, Adv. Chem. Phys. 25 (1974) 69 27) P. Ring, H. Massmann and J. 0. Rasmussen, Nucl. Phys. A296 (1978) 50 28) V. I. Krylov, Approximate calculation of integrals (Macmillan, New York, 1962) pp. 129-130 and 343-346 29) N. Froman and P. 0. Froman, JWKB approximation (North-Hohand, Amsterdam, 1965) pp. 92-97 30) A. Ralston, in Mathematical methods for digital computers, edited by A. Ralston and H. S. Wilf (Wiley, New York, 1960) p. 95