spheroidal wave functions 2

random noise, orthogonal frequency di- vision multiplexing, and anisotropy ofthe cosmic microwave background radi-ation. Therefore, visualizing these func-tions and computing them reliably canbe useful and interesting. What are spheroidal wavefunctions? Spheroidal wave functions are general-izations of Legendre functions andspherical Bessel functions for spheroidalcoordinates rather than for the sphericalpolar coordinates in which the latterfunctions usually occur. The literatureon spheroidal wave functions is often inthe context of specialized applications,but two applied-mathematics mono-graphs are especially useful: those byJulius Stratton and colleagues 1and by Carson Flammer.2Josef Meixner and colleagues have derived analytical re-sults, 3,4and summaries of many results appear in Higher Transcendental Func- tions, edited by Arthur Erdélyi,5and in the Handbook of Mathematical Functions , edited by Milton Abramowitz and IreneStegun. 6For extensive discussions of computational methods—with programsin Mathematica, C, and Fortran—seemy atlas of mathematical functions. 7 Spheroidal coordinates I relate spheroidal coordinates—d, h,x, and f—to Cartesian coordinates—x, y, and z—as in Flammer’s monograph:2For prolate coordinates, , z= dhx/2, -1 £h£1, x‡1, and 0 £f£ 2p. For oblate coordinates, , z= dhx/2, -1 £h£1, x‡0, and 0 £f£2p. The limits xfi¥, dfi0, dx/2= r, and h= cosqproduce spherical polar coordi- nates. Many of the spheroidal coordinatesystems used by other authors, includingAbramowitz and Stegun, do not have thislimit property. Figure 1 illustrates sur-faces corresponding to constant parame-ter values— h, x, or f. The parameter d provides an overall scale factor, as does r in spherical polar coordinates. For detailsof the geometry of spheroidal coordi-nates, see the references by Flammer, 2 Abramowitz and Stegun,6and Parry Moon and Domina Spencer.8Beware, the choice of coordinate systems and no-tations is quite variable!The scalar wave equation inspheroidal coordinates As a context for spheroidal wave func-tions, consider solving the scalar waveequation in spheroidal coordinates. Therelated, but more complicated, vectorwave equation for Maxwell’s equations iscovered by Flammer, by Moon andSpencer, and in technical papers on an-tenna theory and wave scattering byspheroids. The scalar wave equation for wave number k—namely, Ñ2y+ k2y= 0—is separable in spheroidal coordinates bywriting, for prolate coordinates, , in which mis an integer if the fdepen- dence has a period of 2 pand if nis an in- teger. This separability is analogous to thatfor solving the Laplace equation ( k= 0) in spherical polar coordinates. Function S mn is the prolate spheroidal angular function when kis real, because in the limit of small nonsphericity, hbecomes the polar angle q. Parameter cis given by c”kd/2 = pd/l, where lis the wavelength corresponding to wave number k. Thus, cscales as the ra- tio of distance to wavelength. FunctionR mnis the prolate spheroidal radial function , which becomes a spherical Bessel functionin the limit of zero c. Prolate spheroidal functions satisfy Equation 1 for the angular function andEquation 2 for the radial function. (Seethe sidebar for all numbered equations.)In these equations, lmn(c) is the prolate spheroidal eigenvalue , with lmn(0) = n(n+ 1). Oblate coordinates have a similar sep-aration of the wave equation.ψη ξφ φmn mn mnScRcm m=()(),,cos sin xd yd=− + =− +() () () ()211 21122 22ηξ φ ηξ φcos sin xd yd=− − =− −() () () ()211 21122 22ηξ φ ηξ φcos sin 84 COMPUTING IN SCIENCE & ENGINEERINGEditors: Francis Sullivan, [email protected] William J. Thompson, [email protected] PRESCRIPTIONS SPHEROIDAL WAVEFUNCTIONS William J. Thompson, University of North Carolina at Chapel Hill SPHEROIDAL WAVE FUNCTIONS OCCUR IN MANY SCIENTIFIC AND ENGINEERING CONTEXTS, FROM ATOMIC NUCLEI TO THE COSMOS—SCATTERING BY NONSPHERICAL NUCLEI, WAVEFUNCTIONS OF DIATOMIC MOLECULES, ANALYSIS OF BAND-LIMITED 84 COMPUTING IN SCIENCE & ENGINEERING 84 COMPUTING IN SCIENCE & ENGINEERING 84 COMPUTING IN SCIENCE & ENGINEERING 84 COMPUTING IN SCIENCE & ENGINEERING MAY–JUNE1999 85Analytical r esults for spher oidal wave functions ar e commonly pr esented in terms of those for pr olate functions. The transition to the oblate functions (angularor radial functions, or eigenvalues) followsthis r ule: pr olate «oblate by c«–ic, c 2« -c2. However , all spher oidal functions ar e real-valued, in spite of this r ule’s appear - ance. In the following, I give r esults for prolate functions, with the understanding that the r ule is used in analysis (but not in numerical computations!) to obtain r esults for oblate functions. In a dispersive (lossy)medium, the wave number kis a complex number , so cis complex. Le-W ei Li and his colleagues have consider ed this case. 9 Spher oidal wave functions ar e usually expanded in a basis of cor responding spherical functions, with the magnitudeof ccontr olling the range of basis func - tions needed for accurate r esults. Eigenvalues for spheroidal equations The eigenvalues, lmn, are tricky , tedious, and time-consuming to compute accu - rately . I pr esent the necessar y for mulas; for their derivation, see Flammer’ smonograph.2Here I describe a br ute- force method. First, we defi ne two functions that de - pend on m, n, and c(but not on the eigen - value), shown in Equations 3 and 4.Then, we combine these functions to de - fine two functions of lmn, shown in Equa - tions 5 and 6. The fi rst continued fraction terminates with either the ter m contain - ing gm 0or the ter m with gm 1, depending on whether n- mis even or odd, while the second fraction is nonter mi- nating (in principle). As thesecond fraction’ s upper limit increases, the accuracy with which lmncan be deter - mined incr eases. The eigenvalue lmnis the root of the transcendental equation U(lmn) =– U1(lmn) + U2(lmn) = 0. This equation has no closed-for m solutions, except if c= 0 when lmn= n(n + 1) or if the defor mation is large. Accurate eigenvalues are the essential fi rst step for deter mining spher oidal wavefunctions. Appr oximate eigenvalues can be estimated by expanding them as power se - ries in cor as asymptotic series in cand its in - verse powers. The r esulting cumbersome formulas ar e accurate to better than par ts per million only for ver y small or ver y large c. If the inter focal distance of the spher oidal coordinates d»l, a condition that is often interesting, then c»3. For such cvalues, the series for mulas give an accuracy of only a few per cent for most values of mand n. For mand nfrom 0 to 6, and for cwith a magnitude less than 3, the eigenvalueshave a unifor m, slow dependence on de - formation parameter c. For oblate and prolate cases, lmn(c) deviates fr om the spherical coor dinates value, n(n+ 1), in opposite dir ections. This deviation is con - sistent with the leading ter m in a power series expansion in cbeing quadratic. When the magnitude of cis small, the ef - fects of nonsphericity generally becomesmaller as nincreases. To compute the eigenvalues numeri - cally, we can star t with appr oximate solu - tions derived fr om power series or as - ymptotic expansions, then r efine these solutions by a r obust r oot-finding algo - rithm. Using either star ting method givesEquations (1) (2)(3)(4) (5) (6) (7) Rc acjcmnm rmn mr r1 22 0111() + =∞ ()=−() ()() ∑ , ,ξ ξ ξ Umnnmm nmm mnnmm nmm mn22 24 4λβ γ λβ γ λ()≡− −− −−−+ −+−+ −+ K Umn nmm mnnmm nmm mnnmm nmm mn1 22 4λγλβ γ λβ γ λ()≡−− −− −−−− −−−− −− K βrmrr mrmrc m r m r m r≡−()+()+−() +−( )+−( )++( )12 2 1 2 212 232 214 2 γrmmrmrc m m r m r≡+()++()+−− +−( )++( )    1214 1 2 212 232 2 d ddR c dcmRcmn mn mnξξξ ξλξ ξξ2 222 21 10 −()()     −−+ −  ()=,, d ddS c dcmScmn mn mnηηη ηλη ηη 1 102 222 2−()()     +−− −  ()=,, Figure 1. Spheroidal coordinates: (a) prolate; (b) ob - late. h, x, and fare constant parameter values. Coor - dinate surfaces are hyperboloids of revolution for h, and half planes for f.z z j h y x xy x (a) (b). 86 COMPUTING INSCIENCE & ENGINEERINGeigenvalue estimates that ar e usually within 0.2 of the fi nal eigenvalue. W e can therefore use a simple and r obust r oot finder, such as the bisectional method, to locate the r oots. T ypically , a dozen bisec - tions pr oduce par t-per -thousand accu - racy in lmn(c), and appr oximately 30 bi - sections r esults in 10-digit accuracy , the goal for functions in my atlas.7 When cis small in magnitude, the eigenvalue depar ts steadily fr om the spherical-coor dinates value. Y ou might therefore expect that the eigenvalue equation’ s roots ar e unique, as some pr e- vious investigations have assumed. How - ever, as mand nincrease, this is not nec - essarily so.3,4,7Spheroidalangularfunctions Spher oidal angu - lar functions ar e usually expandedinto spherical Le - gendr e functions of the fi rst kind, P m m+r(h), or the sec - ond kind, Qm m+r(h). For functions of the first kind, which are regular at h= ±1, we write , with summation star ting at r= 0 if n-m is even but at r= 1 if n-mis odd. In ei - ther case, rgoes by steps of two. As cfi0, the spher oidal angular function collapses to Pm m+r(h) with the same mand nvalues. The only nonzer o angular coef ficient is then dmn n–m(c), cor responding to r= n-m. For functions of the second kind, whichare irregular at h= ±1, we have . As cfi0, this collapses to Qm n(h), so that the only nonzer o angular coef ficient is dmn n–m(c)—that is, r= n-m. We can compute spher oidal angular coefficients dmn r(c) from the r ecurrence relationwith arand grgiven by and . Recur rence can pr ocede in the dir ection of incr easing or decr easing r. For modest values of n-m, the latter gives mor e ac- curate r esults. (Computing has cer tainly pr ogressed over the last 40 years. The 75,000 nu - merical values that Stratton and col - leagues used r equir ed “about six months of fairly intensive ef fort” from two pr o- grammers and 10 hours of pr oduction time on MIT’ s Whirlwind I computer .1 The output was mor e than fi ve kilometers of paper tape, fr om which they pr epared tables on an electric typewriter . A moder n desktop computer r educes both the exe - cution time and the computer’ s volume by factors of appr oximately 1,000.) Expansion coef ficients for the angular part of spher oidal wave functions depend on the or der n, the degr ee m, and the pa - rameter c, and on whether you ar e using prolate or oblate coor dinates. Figur e 2 displays the drmn(c) values as sur faces made from plaquettes whose ver tices ar e the co - efficient values. The coef ficients for oblate coor dinates behave similarly to those for pr olate coor dinates. If c2is much larger than shown her e, however , the be - havior of the drmn(c) becomes complicated.γrrrc mrmr=−() +−( )+−( )1 2232212αrmrmrc mrmr=++( )++( ) ++( )++( )2221 2232252 α βλ γrrmn r mn rmn rrmnd d d+ −+−() +=2 20,Sc dcmn rmn mrm r2 01() + =∞ ()=()() ∑, ,η ηQSc dcmn rmn mrm r1 01() + =∞ ()=()() ∑, ,η ηPCOMPUTING PRESCRIPTIONS 1 0 02 424 rc 6 8–1dmn r 1 0 02424 rc 6 108–1dmn r1 0 0 2 424 rc 6–1dmn r 1 0 02 424 rc 6 8(a) (c)(b) (d)–1dmn r Figure 2. Coeffi cients for expanding spheroidal wave functions in a basis ofLegendre functions, for mand nvalues of (a) 0, 2; (b) 2, 2; (c) 0, 4; and (d) 2, 4.The coeffi cients peak at r= n-m, which is the unique value when the prolate - ness parameter c= 0. 2 1 0 0.5 1.0η S(1)(c,00)η S(1)(c,00)η S(1)(c,01)η S(1)(c,22 ) c = –i2 c = 2 c = 0(a)3 12 0 0.5 1.0 ηη(b)1 0 0.5 1.0 (c)ηFigure 3. Spheroidal angu - lar functions of the fi rst kind, for prolate ( c= 2), spherical ( c= 0), and oblate (c= -i2) coordinates, for m and nvalues of (a) 0, 0; (b) 2, 2; and (c) 0, 1. MAY–JUNE1999 87Ther e are four ar guments for each spher oidal angular function of the first and second kind: m, n, c, and h. To visual - ize S(1) mn(c,h), we choose çc ç= 2 and super - impose thr ee cur ves for each mand n: c= 2 (prolate case), c= 0 (spherical case), and çc ç= 2 (oblate case). For c= 0, we have the spherical Legendr e functions, P nm(h). For the r egular functions, Figur e 3 shows S(1) mn(c,h)for hover [0, 1]. The spherical angular function is nearly the average ofthe values for pr olate and oblate coor di- nates, indicating that the d rmn(c) are ap- proximately even functions of c. Spheroidal radial functions The spher oidal radial functions, Rmn(c, x), are usually expanded in a basis of spherical Bessel functions.7The expansions ar e quite simple, because radial expansion coef fi- cients ar e proportional to angular expan - sion coef ficients. A given set of spherical- basis radial functions (Bessel, Neumann, orHankel) has cor responding spher oidal ra - dial functions. I discuss the spher oidal ra - dial function called R (1) mn(c, x) by Flammer ,2 by Abramowitz and Stegun,6and in my at - las,7but called jem`(b, x) by Stratton and colleagues,1wher e `= nand h= cin our notation. When c= 0, R(1) mn(c, x) collapses to jn(cx), the spherical Bessel function. The pr ototype spher oidal radial function is regular at x= ±1 and expands in ter ms of regular spherical Bessel functions as shown in Equation 7. In this equation, summationstarts at r= 0 if n -mis even but at r= 1 if n-mis odd. In either case, rgoes by steps of two. The coef ficients armn(c) are radial ex -pansion coef ficients. As cfi0, the spher oidal radial function collapsesto j n(cx), and the only nonzer o ra- dial coef ficient is then amn n–m(0), cor - responding to r= n-m. We can r eadily compute the spher oidal radial coef fi- cients in ter ms of the angular coef ficients drmn. Although the nor malization of the angular coef ficients is dif ferent between Stratton and colleagues1and Flammer ,2 the radial coef ficients ar e the same. Fig - ure 4 displays the armn(c) values as sur faces made fr om plaquettes whose ver tices ar e the radial coef ficient values. T o visualize R(1) mn(c, x) , I choose çc ç= 2 and show two curves for each mand n: çc ç= 2 (pr olate coordinates) and c = -i2 (oblate coor di- nates), as in Figur e 5. After we have computed the angular and radial spher oidal wave functions, we can compute the complete spher oidal wave function, which can help solvemany pr oblems of inter est to scientists and engineers. References 1.J.A. Stratton et al., Spheroidal Wave Functions , Technology Press of M.I.T. and John Wiley &Sons, New York, 1956. 2.C. Flammer, S pheroidal Wave Functions , Stan - ford Univ. Press, Stanford, Calif., 1957. 3.J. Meixner and F.W. Schäfke, Mathieusche Funktionen und Sphäroidfunktionen (MathieuFunctions and Spheroidal Functions), Springer- Verlag, Berlin, 1954. 4.J. Meixner, F.W. Schäfke, and G. Wolf, Math - ieu Functions and Spheroidal Functions andTheir Mathematical Foundations , Springer- Verlag, 1980. 5.A. Erdélyi et al., Higher Transcendental Func - tions, Vol. 3, McGraw-Hill, New York, 1953; reprint edition, Krieger Publishing Co., Mal - abar, Fla., 1981. 6.M. Abramowitz and I.A. Stegun, eds., Hand - book of Mathematical Functions , Dover, New York, 1964. 7.W.J. Thompson, Atlas for Computing Mathe - matical Functions , John Wiley & Sons, 1997, Ch. 13. 8.P. Moon and D.E. Spencer, Field Theory Hand - book, 2nd ed., Springer-Verlag, 1971. 9.L.W. Li et al., “Computations of Spheroidal Harmonics with Complex Arguments: A Re - view with an Algorithm,” Physical Rev. E: Sta - tistical Physics, Plasmas, Fluids, and Related In - terdisciplinary Topics , Vol. 58, No. 5, Nov. 1998, pp. 6792–6806.1 0 02 424 rc 6 8–1amn r 1 0 24 rc –1amn r1 0 24 rc–1amn r 1 0 02 424 rc 6 8(a) (c)(b) (d)–1amn r 0246 1080 2 4 6 Figure 4. Radical coeffi cients for expan - sion into spherical Bessel functions, for m and nvalues of (a) 0, 2; (b) 2, 2; (c) 0, 4; and (d) 2, 4. The coeffi cients peak at r= n-m, the unique value when c= 0. 0.6 2.5 3.0 R(1)(c, )00 c = 2(a)0.2 2.5 3.0R(1)(c, )220.4 0.2 –0.2 (b)0.4 2.5 3.0R(1)(c, )01ξ ξ ξ 0.2 (c) c = –i2ξ ξ ξFigure 5. Spheroidal radial functions of the fi rst kind for prolate ( c= 2) and oblate ( c= -i2) coordinates for mand nvalues of (a) 0, 0; (b) 2, 2; and (c) 0, 1.