spheroidal wave functions 1

Theory and Computation of the Spheroidal Wave FunctionsP. E. Falloon, P. C. Abbott, and J. B. WangSchool of Physics, The University of Western Australia35 Stirling Hwy, Crawley WA 6009 AUSTRALIA.AbstractIn this paper we report on a package, written in the Mathematica computer algebrasystem, which has been developed to compute the spheroidal wave functions ofMeixner [J. Meixner and R.W. Schäfke, Mathieusche Funktionen undSphäroidfunktionen, 1954] and is availlable online(www.physics.uwa.edu.au/~falloon/spheroidal/spheroidal.html). This packagerepresents a substantial contribution to the existing software, since it computes thespheroidal wave functions to arbitrary precision for general complex parameters m, n, g and argument z; existing software can only handle integer m,n and does not givearbitrary precision. The package also incorporates various special cases andcomputes analytic power series and asymptotic expansions in the parameter g. Thespheroidal wave functions of Flammer [C. Flammer, Spheroidal Wave Functions,1957] are included as a special case of Meixner’s more general functions. This paperpresents a concise review of the general theory of spheroidal wave functions and adescription of the formulas and algorithms used in their computation, and giveshigh-precision numerical examples.PACS: 02.30.Gp, 02.70.Wz 1 I. IntroductionSpheroidal wave functions are a class of special functions with many applications in physics andapplied mathematics. They satisfy the differential equation(1)dÅÅÅÅÅÅÅÅdzJH1-z2L dfÅÅÅÅÅÅÅÅÅdzN+ikjjjlnmHgL+g2H1-z2L-m2ÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅ1-z2y{zzzfHzL=0,where m,n,g are arbitrary complex parameters. For many applications m,n take on integer values,in which case they are denoted m,n. Although the general theory of spheroidal wave functions hasbeen known for a long time [1-3], they are still regarded as difficult to compute, and at presentthere are few readily available computer packages available for their computation. Solutions of Eq.(1) were first studied by Niven [4] in connection with a problem involving heatconduction in spheroidal bodies, and were subsequently investigated by a number of authors (see[5] and references therein). Early applications included the quantum mechanical two-centreproblem and various electromagnetic boundary-value problems. The general theory andbackground on spheroidal wave functions is contained in the monograph by Flammer [2]. Othernotable monographs are Stratton et al. [3], which has extensive tables of numerical values, andKomarov et al. [6]. These works focus exclusively on spheroidal wave functions with integerparameters m,n. Meixner developed the theory of spheroidal wave functions with arbitrarycomplex parameters m,n (see [1, 7] and references therein). Several packages have been developed recently to compute spheroidal wave functions: Thompson[8] (which uses an incorrect expansion for the angular functions of the second kind) and Li et al.[9] are two of the most recent. Both of these packages are only useful for small values of g, do notprovide arbitrary precision computation, and are limited to integer parameters m,n. Furthermore,no package currently in existence computes the power series and asymptotic expansions for thespheroidal wave functions. The choice of a suitable notation and normalization for spheroidal wave functions presents asignificant challenge, due to the large number of conventions in existence. The two main ones arethose of Meixner [1] and Flammer [2] (also used in [10]). The latter is more commonly used in theliterature, however its use is usually limited to integer parameters. Indeed, Flammer’snormalization scheme cannot readily be generalized to noninteger parameters, because itnormalizes functions to a different constant depending on the parity of n-m. Meixner’s notation,though less commonly used, has the fundamental advantage that it is suitable for (and was in factdeveloped specifically to handle) the case of general complex parameters.2 The choice of a suitable notation and normalization for spheroidal wave functions presents asignificant challenge, due to the large number of conventions in existence. The two main ones arethose of Meixner [1] and Flammer [2] (also used in [10]). The latter is more commonly used in theliterature, however its use is usually limited to integer parameters. Indeed, Flammer’snormalization scheme cannot readily be generalized to noninteger parameters, because itnormalizes functions to a different constant depending on the parity of n-m. Meixner’s notation,though less commonly used, has the fundamental advantage that it is suitable for (and was in factdeveloped specifically to handle) the case of general complex parameters.The purpose of this paper is to describe a package which has recently been developed [12] tocompute spheroidal wave functions with arbitrary parameter values, which overcomes theshortcomings of existing software mentioned above. We have chosen to use the Mathematicacomputer algebra system [11], which is ideal due to its symbolic and high-precision numericalcapabilities, as well as its large library of built-in special functions. We have taken a uniqueapproach to the notation for the spheroidal functions, in that our package computes the functions ofFlammer and Meixner as two distinct sets of functions. In this way, it is hoped that the packagewill be useful to as wide an audience as possible.The layout of this paper is as follows. In Section II we present a concise review of the theory of thespheroidal wave functions, starting with their definition as series of Legendre and spherical Besselfunctions. We then discuss the important special case for the angular functions of the second kindwhen m+n is an integer, before discussing the spheroidal joining factor which relates the angularand radial functions. In Section III we describe the computation of the spheroidal wave functions,beginning with a discussion of the continued fraction and tridiagonal matrix methods used tocompute the spheroidal eigenvalues. We then discuss the general numerical implementation of thespheroidal functions, and finally describe the approach used to generate the asymptotic and powerseries coefficients. We conclude with a discussion of some of the ways in which the accuracy ofour numerical values can be verified. Four appendices are also included: Appendix A containsdefinitions for the Legendre and spherical Bessel functions which are used in the package;Appendix B contains formulas relating Flammer’s spheroidal functions to those of Meixner;Appendix C contains a summary of some important mathematical identities satisfied by thespheroidal wave functions; finally, in Appendix D we present tables of sample function values tohigh-precision, for the purposes of comparison with other packages.II. Theory2.1 The angular spheroidal wave functionsWhen g=0, Eq.(1) reduces to Legendre’s differential equation(2)dÅÅÅÅÅÅÅÅdzJH1-z2LdfÅÅÅÅÅÅÅÅÅdzN+ikjjjlnmH0L-m2ÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅ1-z2y{zzzfHzL=0.The solutions to this equation are the (associated) Legendre functions of the first and second kind,PnmHzL and QnmHzL, with eigenvalue lnmH0L=nHn+1L. Traditionally, these functions are defineddifferently depending on whether or not z lies on the branch cut zœH-1,1L (e.g. [10], §8.1, 8.3).For many applications—particularly those involving computer algebra—this convention isinconvenient, since it precludes the use of identities which are valid for all z. An alternativeapproach, which avoids this difficulty is to define two distinct types of Legendre function, each ofwhich is valid for all values of z. The functions of Type I, PnmHzL and QnmHzL, are equivalent to thoseusually defined on the H-1,1L cut. The functions of Type II, which we denote /GothicCapPnmHzL and /GothicCapQnmHzL, areequivalent to the functions usually defined for z–H-1,1L. In Appendix A we give the definitionsand some important properties of these functions.3 The solutions to this equation are the (associated) Legendre functions of the first and second kind,PnmHzL and QnmHzL, with eigenvalue lnmH0L=nHn+1L. Traditionally, these functions are defineddifferently depending on whether or not z lies on the branch cut zœH-1,1L (e.g. [10], §8.1, 8.3).For many applications—particularly those involving computer algebra—this convention isinconvenient, since it precludes the use of identities which are valid for all z. An alternativeapproach, which avoids this difficulty is to define two distinct types of Legendre function, each ofwhich is valid for all values of z. The functions of Type I, PnmHzL and QnmHzL, are equivalent to thoseusually defined on the H-1,1L cut. The functions of Type II, which we denote /GothicCapPnmHzL and /GothicCapQnmHzL, areequivalent to the functions usually defined for z–H-1,1L. In Appendix A we give the definitionsand some important properties of these functions.For nonzero g, the angular spheroidal wave functions are defined as infinite series of thecorresponding Legendre functions:(3)FnmHz;gL=‚k=-¶¶ H-1Lkan,kmHgLfn+2 kmHzL.Here we follow Meixner’s notation so that F=ps,qs,Ps,Qs and f=P,Q,/GothicCapP,/GothicCapQ respectively.It can readily be shown [1] that the series coefficients an,kmHgL satisfy the three-term recurrencerelation(4a)An,kmHgL an,k-1mHgL+HBn,kmHgL-lnmHgLLan,kmHgL+Cn,kmHgLan,k+1mHgL=0,where(4b)An,kmHgL=-g2Hn-m+2 k-1L Hn-m+2 kLÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅH2n+4 k-3LH2n+4 k-1L,Bn,kmHgL=Hn+2 kL Hn+2 k+1L-2 g2Hn+2 kL Hn+2 k+1L+m2-1ÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅH2 n+4 k-1L H2 n+4 k+3L,Cn,kmHgL=-g2 Hn+m+2 k+1LHn+m+2 k+2LÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅH2 n+4 k+3LH2 n+4 k+5L.The series expansion in Eq.(3) is convergent only when the coefficients an,kmHgL form a minimalsolution to Eq.(4), (i.e. a solution with the property that an,kmHgLêan,k°1mHgLØ0 as kØ≤¶ [16])—in which case it converges for all z. There is a countably infinite set of values for l that correspondto minimal solutions. The spheroidal eigenvalue lnmHgL is defined as a function of m, n and g bychoosing the l value that reduces to nHn+1L continuously as gØ0. In practice, this assignment isnontrivial and must be made numerically.4 For integers m,n with n¥†m§, the coefficients an,kmHgL are normalized so that‡-11 psnmHt;gL2 dt=2ÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅ2 n+1Hn+mL!ÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅHn-mL!,which can be generalized in a natural way to give the relation for general m,n:(5)„k=-¶¶ an,kmHgL2 2 n+1ÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅ2 n+4 k+1 Hn+m+1L2 kÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅHn-m+1L2 k=1.The sign of an,kmHgL is determined by the condition that psnmHz;gLØPnmHzL continuously as gØ0.2.2 Angular functions of the second kind for integer m+nThe functions QnmHzL and /GothicCapQnmHzL diverge when m+n is a negative integer (Eq.A9), from which itimmediately follows that (6)†qsnmHz;gL§,†QsnmHz;gL§=¶,m+n=-1,-2,…For m+n a non-negative integer, we have Cn,H-m-n-2+dLê2mHgL=0, where d=Hm+nLmod2. FromEq.(4) we then find an,kmHgL=0 for k<-Hm+nLê2, and hence the series (6) becomes indeterminatefor qsnmHz;gL and QsnmHz;gL, since the infinite basis functions are multiplied by zero seriescoefficients. The procedure for recovering a valid series representation in this case is reasonablystraightforward, and is described (for integer parameters m,n) by Flammer [2].The essential step is to use the transformation relations Eqs.(A10-11). Taking the case of /GothicCapQnmHzL fordefiniteness, we substitute nØn+e, where m+n is an integer and eØ0, into Eq.(A11) andimmediately findsinHepL /GothicCapQn+emHzL=H-1Lm+n HpeimpcosHHn+eLpL/GothicCapP-n-e-1mHzL-sinHHm-n-eLpL/GothicCapQ-n-e-1mHzLL.Multiplying by the series coefficient an+e,kmHgL and taking the limit eØ0 we have(7)limeØ0an+e,kmHgL /GothicCapQn+2 k+emHzL=aèn,kmHgL H-1Lm+n JeimpcosHnpL/GothicCapP-n-2 k-1mHzL-1ÅÅÅÅÅp sinHHm-nLpL/GothicCapQ-n-2 k-1mHzLN,where we define(8)aèn,kmHgL=limeØ0an+e,kmHgLÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅe,k§k0-1.For k§k0-2, the coefficients aèn,kmHgL can be computed using5 (9)aèn,kmHgLÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅaèn,k+1mHgL=-CkÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅBk-lnmHgL+Ak aèn,k-1mHgLÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅaèn,kmHgL,while for k=k0-1, we have(10)aèn,k0-1mHgL=-CènmHgL an,k0mHgLÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅBk0-1-lnmHgL+Ak0-1 aèn,k0-2mHgLÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅaèn,k0-1mHgL,where(11)CènmHgL=limrØ0Cn+r,k0-1mHgLÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅr=H-1Ldg2ÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅH2 m-2 d-1LH2 m-2 d+1L.For m+n=0,1,2,… the series for QsnmHz;gL therefore reads(12a)QsnmHz;gL=‚k=k0¶ H-1Lk an,kmHgL/GothicCapQn+2 kmHzL+H-1Lm+n ‚k=-¶k0-1H-1Lk aèn,kmHgLµJ‰ÂmpcosHnpL/GothicCapP-n-2 k-1mHzL-sinHHm-nLpLÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅp /GothicCapQ-n-2 k-1mHzLy{zz.For qsnmHz;gL we follow a directly anagolous argument starting from Eq.(A10) in place of Eq.(A11), and obtain(12b)qsnmHz;gL=‚k=k0¶ H-1Lk an,kmHgLQn+2 kmHzL+H-1Lm+n ‚k=-¶k0-1H-1Lk aèn,kmHgLµJcosHmpLcosHnpLP-n-2 k-1mHzL-sinHHm-nLpLÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅp Q-n-2 k-1mHzLy{zz.For integers m,n these expressions reduce to(13a)QsnmHz;gL=‚k=k0¶ H-1Lk an,kmHgL/GothicCapQn+2 kmHzL+‚k=-¶k0-1 H-1Lk aèn,kmHgL/GothicCapP-n-2 k-1mHzL,(13b)qsnmHz;gL=‚k=k0¶ H-1Lk an,kmHgLQn+2 kmHzL+‚k=-¶k0-1 H-1Lk aèn,kmHgLP-n-2 k-1mHzL.2.3 The radial spheroidal wave functionsIn the limit gØ0 and zض such that gz=constant, the two regular singularities of Eq.(1) (atz=≤1) coalesce. This can be seen by changing variables to z=gz, substitutingfHzL=H1-1êz2Lmê2 gHzL, and letting gØ0. Eq.(1) then becomes6 (14)z2d2gÅÅÅÅÅÅÅÅÅÅÅÅdz2+2 zdgÅÅÅÅÅÅÅÅÅdz+Hz2-lLgHzL=0,which is satisfied by the spherical Bessel functions jnHzL and ynHzL ([10], Ch.10 and Appendix A).The radial spheroidal wave functions are defined in terms of these by(15)SnmHkLHz;gL=H1-1êz2Lmê2ÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅAn-mHgL ‚k=-¶¶an,k-mHgLfn+2 kHgzL,where k=1,2 and f=j,y respectively, and(16)AnmHgL=‚k=-¶¶H-1Lk an,kmHgL.It has been shown [13] that the functions defined in Eq.(15) are indeed solutions to Eq.(1) whichare absolutely convergent for †z§>1. The normalization factor AnmHgL is chosen so that thefollowing limits are satisfied:(17)SnmH1LHz;gLöøøøgzض1ÅÅÅÅÅÅÅÅÅgzsinJgz-npÅÅÅÅÅÅÅÅÅÅ2N,SnmH2LHz;gLöøøøgzض-1ÅÅÅÅÅÅÅÅÅgz cosJgz-npÅÅÅÅÅÅÅÅÅÅ2N.The functions SnmHkLHz;gL have branch cuts in the complex z-plane along the line H-1êg¶,0L, forn–/DoubleCapZ, and on the interval H-1,1L, for mê2–/DoubleCapZ.2.4 Joining relations between angular and radial functionsThe angular and radial spheroidal wave functions can both be considered as functions over theentire complex z-plane. From a computational point of view, however, their series expansions areonly useful over a restricted subset of the complex plane. For the angular functions, the series (3) isconvergent over the entire z-plane, but for †z§>1 it becomes too slowly convergent to be of anypractical use. For the radial functions the situation is even worse—the series (15) is in general notconvergent inside the unit circle †z§<1. To allow computation of the functions over the entirecomplex plane, Meixner and Schäfke [1] introduced a joining factor that relates the angular andradial functions.Using well-known series expansions for /GothicCapQnmHzL and jnHzL, it is possible to find the following seriesexpansions for QsnmHz;gL and SnmH1LHz;gL [12]:7 (18a)Qsnm Hz;gL=2-n-1 è!!!p eimpz-m-n-1Hz-1Lmê2Hz+1Lmê2µ„j-¶¶ 1ÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅH2 zL2 j„k=0¶ H-1Lj-k GH2 j+m+n+1LÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅk! GH2 j-k+n+3ÅÅÅÅ2L an,2 j-2 kmHgL,(18b)SnmH1LHz;gL=è!!!p H1-1êz2Lmê2 HgzLnÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅ2n+1 An-mHgL „j=-¶¶1ÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅH2 zL2 j „k=0¶H-1Lk 24 j an,-2 j-2 k-mHgLÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅGH-2 j-k+n+3ÅÅÅÅ2Lk! g2 j.Comparison of Eqs.(18a) and (18b) reveals that the expansions of Qs-n-1mHz;gL and SnmH1LHz;gLinvolve identical powers of z. Furthermore, it is obvious from Eq.(A13) that the expansion forSnmH2LHz;gL will not involve the same powers of z. Now, Qs-n-1mHz;gL must be expressible as alinear combination of the functions SnmH1,2LHz;gL, since it satisfies the same differential equation, soit follows that Qs-n-1mHz;gL and SnmH1LHz;gL must be equal up to a constant factor. Comparing (18a)and (18b), we have then that the ratioikjjjjjjjj„k=0¶H-1Lk 24 j an,-2 j-2 k-mHgLÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅGH-2 j-k+n+3ÅÅÅÅ2Lk! g2 jy{zzzzzzzzìikjjjjjjjj„k=0¶ H-1Lj-k GH2 j+m-nL an,-2 j+2 kmHgLÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅGH2 j-k-n+1ÅÅÅÅ2Lk!y{zzzzzzzzmust be independent of j.In light of the above result, we define the spheroidal joining factor KnmHgL by the relation(19)SnmH1LHz;gL=KnmHgL sinHHm-nLpLÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅp e-iHm+nLp H1-1êz2Lmê2 HgzLnÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅgnzn-m Hz-1Lmê2 Hz+1Lmê2 Qs-n-1mHz;gL.The trigonometric and exponential factors are included so that the joining factor relations reduce toa simple form for integers m,n. Note also the factor involving various powers of z, which includesthe branch cut information for the two functions. Comparing terms in Eqs.(18a) and (18b) for anyparticular j we can obtain an explicit definition for KnmHgL. For definiteness we choose j=0 andobtain(20)Knm HgL=einp 2-2 n-1 GHn-m+1L gnÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅAn-mHgLµikjjjjjjjj„k=0¶H-1Lk an,-2 k-mHgLÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅGH-k+n+3ÅÅÅÅ2Lk!y{zzzzzzzzìikjjjjjjjj„k=0¶ an,2 kmHgLÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅGH-k-n+1ÅÅÅÅ2Lk!y{zzzzzzzz.Using this equation and the symmetry relations given in Appendix C it is possible to obtain(21)An-mHgL K-n-1mHgL AnmHgL Kn-mHgL=pÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅgsinHHm+nLpL,8 which is useful in constructing joining relations for integers m,n.III. Numerical computation3.1 The spheroidal eigenvalues lnmHgLAs we mentioned in Section 2.1, the spheroidal eigenvalues lnmHgL are minimal solutions of thethree-term recurrence (4). There are two standard procedures for finding such solutions. The firstwas developed independently by Bouwkamp [14] and Blanch [15], and makes use of afundamental equivalence between three-term recurrences and continued fractions [16]. Thisprovides a method for determining the eigenvalues numerically to high precision, although it relieson the availability of a sufficiently accurate starting estimate for the eigenvalue. The secondmethod, due to Hodge [17], involves expressing the three-term recurrence as an infinite tridiagonalmatrix equation. It is complementary to the first in the sense that it provides an excellent methodfor generating accurate starting estimates for the eigenvalues, but is inefficient for obtaining highprecision eigenvalues. We now discuss both methods in turn.Continued fraction methodDefining(22)an,kmHgL=An,kmHgLCn,k-1mHgL,bn,kmHgL=Bn,kmHgL,Nn,kmHgL=Cn,k-1mHgL an,kmHgLÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅan,k-1mHgL,the three-term recurrence (4a) can be rewritten in ascending and descending form asNn,k+1mHgL=an,k+1mHgLÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅbn,k+1mHgL-lnmHgL-Nn,k+2mHgL,Nn,k+1mHgL=bn,kmHgL-lnmHgL-an,kmHgLÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅNn,kmHgL.Setting k=0 and iterating these relations we obtain(23a)/ScriptCapUnmH1LHg,lL+/ScriptCapUnmH2LHg,lL=0,where we have defined9 (23b)/ScriptCapUnmH1LHg,lL=b0-l-a0ÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅb-1-l- a-1ÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅb-2-l- ∫/ScriptCapUnmH2LHg,lL=-a1ÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅb1-l- a2ÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅb2-l- a3ÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅb3-l- ∫Here we are using a standard notational convention for continued fractions [18]:a1ÅÅÅÅÅÅÅÅÅÅÅÅÅÅb1+ a2ÅÅÅÅÅÅÅÅÅÅÅÅÅÅb2+ a3ÅÅÅÅÅÅÅÅÅÅÅÅÅÅb3+ ∫=a1ÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅb1+a2ÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅb2+a3ÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅb3+∫ .Eq. (23a) is a transcendental equation in l, whose roots are the spheroidal eigenvalues ln+2 kmHgL.The method of Bouwkamp and Blanch consists of differentiating the left side of Eq.(23a) andusing Newton’s method. Differentiating Eq. (23b) with respect to l we obtain(24)∑/ScriptCapUnmH1LHg,lLÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅ∑l=-ikjjj1+a0ÅÅÅÅÅÅÅÅÅÅN02+a0ÅÅÅÅÅÅÅÅÅÅN02 a-1ÅÅÅÅÅÅÅÅÅÅÅÅÅN-12+a0ÅÅÅÅÅÅÅÅÅÅN02 a-1ÅÅÅÅÅÅÅÅÅÅÅÅÅN-12 a-2ÅÅÅÅÅÅÅÅÅÅÅÅÅN-22+∫y{zzz,∑/ScriptCapUnmH2LHg,lLÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅ∑l=-ikjjjN12ÅÅÅÅÅÅÅÅÅÅa1+N12ÅÅÅÅÅÅÅÅÅÅa1 N22ÅÅÅÅÅÅÅÅÅÅa2+N12ÅÅÅÅÅÅÅÅÅÅa1 N22ÅÅÅÅÅÅÅÅÅÅa2 N32ÅÅÅÅÅÅÅÅÅÅa3+∫y{zzz.To apply Newton’s method to Eq.(23) we begin with a starting estimate l0 and iterate Newton’sformula,(25)dli=-/ScriptCapUnmH1LHg,li-1L+/ScriptCapUnmH2LHg,li-1LÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅ∑/ScriptCapUnmH1LHg,li-1LÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅ∑l+∑/ScriptCapUnmH2LHg,li-1LÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅ∑l,starting with i=1, until the size of the ith iterate dli is beneath the desired level of precision.Provided sufficient precision is used throughout the application of this algorithm, it can be used toobtain lnmHgL to arbitrarily high accuracy. Usually the intermediate working precision needs to besignificantly higher than the desired final precision, as precision is almost always lost in theprocess of numerical computation.Tridiagonal matrix methodIn order to use the continued fraction method just described, it is necessary to begin with areasonably accurate starting value. Traditionally, most authors have used a power series expansionof lnmHgL in powers of g [2], which has a very small radius of convergence (approximately between4 and 5). Hodge [17] appears to have been the first to solve the recurrence Eq.(4) by recasting it asa tridiagonal matrix equation:10 (26)ikjjjjjjjjjjjjjjjjjjjjjjjjjjj.....A-2B-2C-2A0B0C0A2B2C2.....y{zzzzzzzzzzzzzzzzzzzzzzzzzzz ikjjjjjjjjjjjjjjjjjjjjjjjjjjj..a-2a0a2..y{zzzzzzzzzzzzzzzzzzzzzzzzzzz=lnmHgLikjjjjjjjjjjjjjjjjjjjjjjjjjjj..a-2a0a2..y{zzzzzzzzzzzzzzzzzzzzzzzzzzz.Truncating this equation in both directions yields a finite matrix, the eigenvalues of which will beapproximations to the actual eigenvalues lnmHgL. Since eigenvalues of tridiagonal matrices can beevaluated very efficiently, this method allows approximate eigenvalues to be generated rapidly.These matrix eigenvalues are ideal starting estimates for the continued fraction method describedabove. Surprisingly, to our knowledge no author seems to have previously used this hybridapproach.One issue which presents a practical challenge to the computation of the spheroidal eigenvalues isactually deciding which matrix eigenvalue corresponds to which value of n: all that can be said ingeneral is that each eigenvalue of the matrix corresponds to ln+2 kmHgL for some integer k. When g2is real and m,n are integers, the eigenvalues are strictly ordered and hence the correspondence istrivial. However, in the complex g plane the eigenvalues have a complicated branch cut structure,and the ordering relation does not apply. The same is true for non-integer m,n. One rathercumbersome solution to the problem is to start with an eigenvalue lnmHg0L for which the orderingrelation does hold (i.e. m,n are integers g0 is on the real or imaginary axis), and then follow acurve in the m,n,g parameter space towards the desired eigenvalue lnmHgL. Provided the step sizesare small enough, each eigenvalue can be used to choose the correct starting value from the matrixat the subsequent step. However, this approach leads to an intractable amount of computationunless parameter values are very small. In the package, therefore, we simply include an optionalextra argument to the eigenvalue function, which allows a starting estimate to be given by the user.Using this, it is straightforward to implement the above iterative procedure for particular cases.3.2 The spheroidal wave functionsOnce the eigenvalue lnmHgL has been found, it is relatively straightforward to compute the actualspheroidal wave functions using Eqs.(3) and (15). There are two main parts to this computation:generating the series coefficients an,kmHgL and generating the basis functions (Legendre or sphericalBessel functions). We now discuss each of these in turn.11 Series coefficientsTo generate the series coefficients an,kmHgL we follow the standard approach of rewriting therecurrence relation (4) in terms of ascending and descending ratios:(27)an,kmHgLÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅan,k+1mHgL=-Cn,kmHgLÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅBn,kmHgL-lnmHgL+An,kmHgLan,k-1mHgLÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅan,kmHgL,an,kmHgLÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅan,k-1mHgL=-An,kmHgLÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅBn,kmHgL-lnmHgL+Cn,kmHgLan,k+1mHgLÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅan,kmHgL.These ratios converge like 1êk2 as kØ≤¶, so we can set them to approximately to zero for somelarge value †k§=kmax. We then iterate Eq.(27) from k=≤kmax to k=0 to obtain an,kmHgLêan,0mHgLfor -kmax§k§kmax. Once again, in this recursive process precision is usually lost with each step,and hence it is necessary to work with an intermediate precision substantially greater than the finaldesired precision. In practice, our package begins with a working precision of 100 extra digits, thentests at the end whether the final precision is high enough. If it is not, the working precision isincreased by 100 and the coefficients computed again.The normalization relation (5) can be rewritten in the forman,0mHgL>ikjjjjjjjjj„k=-kmaxkmax ikjjjjan,kmHgLÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅan,0mHgLy{zzzz2 2 n+1ÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅ2 n+2 k+1 Hn+m+1LkÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅHn-m+1Lky{zzzzzzzzz-1ê2,allowing us to obtain an,0mHgL and hence the correctly normalized list of coefficients an,kmHgL for-kmax§k§kmax.The question is now whether or not the chosen value of kmax is large enough that the omitted “tail”of the series is insignificant to the given level of precision. To answer this question rigorously isnot easy, and in the package we settle on simply checking that the magnitude of an,≤kmHgL isnegligible compared to the largest series coefficient, i.e.†an,≤kmaxmHgL§ÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅmaxH8†an,kmHgL§»-kmax§k§kmax<L<10-prec.Although not mathematically rigorous, extensive numerical testing shows that this is an eminentlyreasonable condition.12 Basis functionsIn principle, the basis functions could be computed using Mathematica’s built-in functions.However, it is much more efficient to start with the basis functions with k=0,1 and then userecurrence relations satisfied by the basis functions to generate the basis functions for all othervalues of k. The required relations are(28)Hm-n-2LHm-n-1LÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅH2n+1LH2n+3L fn+2mHzL+ikjjjH2nHn+1L-2m2-1LÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅH2n-1LH2n+3L-z2y{zzzfnmHzL+Hm+n-1LHm+nLÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅH2n-1LH2n+1L fn-2mHzL=0,for f=/GothicCapP,/GothicCapQ,P,Q, and(29)fn-2HzLÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅ2 n-1+H2 n+1LJ2ÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅH2 n-1L H2 n+3L-1ÅÅÅÅÅÅÅÅz2N fnHzL+fn+2HzLÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅ2 n+3=0,for f=j,y. These relations are readily obtained from the standard recurrence formulas §8.5.3 and§10.1.19 in [10]. Once again, these computations must be carried out at a much higher precisionthan is required at the end.3.3 Series expansionsThe eigenvalues lnmHgL and the ratios an,kmHgLêan,0mHgL can be expanded in powers of g2 [2], leadingto approximations which are useful for †g§d4 -5. However, for these values of g the tridiagonalmatrix approach of Section 3.1 is much more useful for obtaining numerical values, so in practicethe power series is mainly of theoretical interest. We now describe how the power seriesexpansions are computed in our package.We begin with the following ansätze:(30)lnmHgL=‚j=0¶/ScriptAltLjmng2 j,an,kmHgLÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅan,0mHgL=‚j=0¶ajkmng2 j,k=0,≤1,≤2,…By inspection we immediately have(31)a0,0mn=1,aj,0mn=a0,kmn=0,forj,k∫0.Substituting the expansions in (30) into the recurrence relation (4) with k=0, and defining13 Amnk=Hm+n+k+1LHm+n+k+2LÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅH2 n+2 k+3LH2 n+2 k+5L,Bmnk=1ÅÅÅÅÅ2ikjjj1-4 m2-1ÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅH2 n+2 k-1LH2 n+2 k+3Ly{zzz,Cmnk=Hn-m+kLHn-m+k-1LÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅH2 n+2 k-3LH2 n+2 k-1L,we find(32)nHn+1L-/ScriptAltL0mnÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅg2+bm,n,0-/ScriptAltL1mn+‚j=1¶HAm,n,0aj,1mn+Cm,n,0aj,-1mn-/ScriptAltLj+1mnL g2 j=0.Since this is true for all values of g it follows that(33)/ScriptAltL0mn=nHn+1L,/ScriptAltL1mn=Bm,n,0=1ÅÅÅÅÅ2ikjjj1-4 m2-1ÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅH2 n-1LH2 n+3Ly{zzz,/ScriptAltLjmn=Am,n,0aj-1,1mn+Cm,n,0aj-1,-1mn,j=2,3,…This provides a recursive method for determining /ScriptAltLjmn when the coefficients aj-1,≤1mn are known. Tofind these coefficients we substitute Eq. (30) into Eq. (4) for general k. After some manipulationwe find(34)ajkmn=1ÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅHn+kLHn+k+1L-/ScriptAltL0mnµikjjjjjj‚i=0j-1/ScriptAltLi+1mnaj-i-1,kmn-HAmnkaj-1,k+2mn+Bmnkaj-1,kmn+Cmnkaj-1,k-2mnLy{zzzzzz.Now let k be a positive integer and suppose that aj-1,k-2mn=0 and arsmn=0 for all r§j-1 ands>k. Then Eq.(34) shows that ajkmn=0. This fact, together with the initial conditions a0,kmn=0 fork=≤2,≤4,…, proves by induction that(35)ajkmn=0,j<†k§ê2,and hence an,≤kmHgLêan,0mHgL=OHg2 kL for k=0,1,2,… This result could also have been deduceddirectly from the recurrence Eq.(4). Eqs.(33-35) constitute the recursive scheme by which wecompute the power series expansions (30a-b). 14 The key to an efficient implementation of recursive algorithms of this kind is to use dynamicprogramming, whereby coefficients are “cached” once they are generated, allowing subsequentcoefficients to be generated more quickly. In Mathematica, this is achieved with a definition of theform fHx_L:=fHxL=…With the expansions for the ratios an,kmHgLêan,0mHgL computed, it is not difficult to obtain thecorresponding expansions for the angular functions themselves. One extra step is required,however, since the coefficient an,0mHgL is itself a function of g and so must be expanded as well.This is readily accomplished using the relationan,0mHgL=ikjjjjjjjj„k=-¶¶ ikjjjjan,kmHgLÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅan,0mHgLy{zzzz2 2 n+1ÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅ2 n+2 k+1 Hn+m+1LkÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅHn-m+1Lky{zzzzzzzz-1ê2.For large g, asymptotic series expansions are also known for the spheroidal wave functions. Forinteger m,n the angular functions reduce to Hermite/Laguerre polynomials as gض/gØi¶, andasymptotic expansions in descending powers of g can be found [2]. In our package, theseasymptotic expansions are computed using the method just described for the power series.IV. Results and DiscussionThe most comprehensive sources of tabulated values of spheroidal functions are Flammer [2],Stratton et al. [3] and Van Buren et al. [19]. Some of the tables from Flammer are reproduced inTables 21.1-4 of Abramowitz and Stegun [10]. Some minor errors in these works have alreadybeen pointed out by Li et al. [9]. We have compared the output of our package with all of thesesources, and found agreement up to the precision to which the tabulated values are given. With the availability of packages such as the one we have developed, there is clearly little need forexhaustive tabulations of numerical values. However, it is useful to give some high-precisionnumerical values to provide a benchmark for comparison with other programs. In Appendix D weprovide a set of such values. They are presented to 25 digits and are intended only to represent anillustrative sample.Because no tables or software packages presently available are capable of generating results toarbitrary precision, the most reliable way to check the validity of our numerical functions is toperform self-consistency tests. Although there are relatively few analytic results available for thespheroidal wave functions, there are several important tests (each of which we have applied to ourpackage with perfect results):15 Because no tables or software packages presently available are capable of generating results toarbitrary precision, the most reliable way to check the validity of our numerical functions is toperform self-consistency tests. Although there are relatively few analytic results available for thespheroidal wave functions, there are several important tests (each of which we have applied to ourpackage with perfect results):ËExact solutions: for n=1,2,… we have ln1Hnpê2L=0 [2]. Also, for m=1ê2, the spheroidalfunctions reduce to the Mathieu functions ([10], Ch.20), and the eigenvalues are related to theMathieu characteristic values arHqL by ln1ê2HgL=an+1ê2Hg2ê4L-g2ê2-1ê4 [12]. We cantherefore compare the numerical eigenvalues generated by our package to the built-in Mathieufunctions in Mathematica.ËWronskian: for the spheroidal functions the Wronskian is proportional to Hz2-1L-1, and forthe radial functions it is equal to HgHz2-1LL-1. This test is the most generally useful, since itcan be used for all parameter values.ËSubstitution into the differential equation: it is straightforward to simply substitute thefunctions back into the spheroidal differential equation (1) and verify that it is satisfied to theprecision of the numerical functions. However, because the second derivative has to becomputed numerically, this is not very convenient for testing to a very high level of precision.ËLastly, a simple check that can always be performed is to generate a certain function value totwo different levels of precision (for example 50 and 100 digits). If the two results do notagree up to the precision of the least precise of the two, this indicates there is a problem in themethod of calculation. If they do agree, however, it does not guarantee anything—since, forexample, the truncated series from which the function is being computed may contain too fewterms—but it is a useful guide.Our numerical package, Spheroidal.m, is available online at the URLwww.physics.uwa.edu.au/~falloon/spheroidal/spheroidal.html, along with further documentationregarding its use.AcknowledgmentsPEF is grateful for the support of University Postgraduate Award from the University of WesternAustralia. Michael Trott and Oleg Marichev of Wolfram Research Inc. provided useful informationabout general issues concerning the numerical implementation of special functions in Mathematica.16 References[1] J. Meixner and R.W. Schäfke, Mathieusche Funktionen und Sphäroidfunktionen (Springer-Verlag, Berlin, 1954) [In German].[2] C. Flammer, Spheroidal Wave Functions (Stanford University Press, Stanford, 1957).[3] J.A. Stratton, P.M. Morse, L.J. Chu and R.A. Hutner, Elliptic Cylinder and Spheroidal WaveFunctions (John Wiley and Sons, New York, 1941).[4] C. Niven, Philos. Trans. Roy. Soc. (London) 171, 117 (1880).[5] M.J.O. Strutt, Lamésche, Mathiéusche und verwandte Funktionen in Physik und Technik(Ergebnisse der Mathematik und ihrer Grenzgebiete Vol. 1 No. 3) (Verlag Julius Springer, Berlin,1932) [In German].[6] I.V. Komarov, L.I. Ponomarev and S.Y. Slavyanov, Spheroidal and Coulomb SpheroidalFunctions (Nauka, Moscow, 1976) [In Russian].[7] J. Meixner, R.W. Schäfke, G. Wolf, Mathieu Functions and Spheroidal Functions and TheirMathematical Foundations (further studies), Lecture Notes in Mathematics 837 (Springer-Verlag,Berlin, 1980).[8] W.J. Thompson, Spheroidal Wave Functions, Computing in Science and Eng. 1(3), 84 (1999).[9] L.W. Li, M.S. Leong, T.S. Yeo, P.S. Kooi, and K.Y. Tan, Computations of spheroidalharmonics with complex arguments: A review with an algorithm, Phys. Rev. E 58, 6792 (1998).[10] M. Abramowitz and I. Stegun, Handbook of Mathematical Functions (Dover, New York,1965).[11] S. Wolfram, The Mathematica Book (Wolfram Media/Cambridge University Press, 1999).[12] P.E. Falloon, Theory and Computation of Spheroidal Harmonics with General ComplexParameters, Masters Thesis, The University of Western Australia, 2001. [13] E.W. Leaver, Solutions to a generalized spheroidal wave equation, J. Math. Phys. 27, 1238(1986).[14] C.J. Bouwkamp, On Spheroidal Wave Functions of Order Zero, J. Math. Phys. 26, 79 (1947).17 [15] G. Blanch, On the computation of Mathieu functions, J. Math. Phys. 25, 1 (1946).[16] W. Gautschi, Computational aspects of three-term recurrence relations, SIAM Review 9, 24(1967).[17] D.B. Hodge, Eigenvalues and Eigenfunctions of the Spheroidal Wave Equation, J. Math. Phys.11, 2308 (1970).[18] H.S. Wall, Analytic Theory of Continued Fractions (Van Nostrand, New York, 1948).[19] A.L. Van Buren, B.J. King, R.V. Baier, and S. Hanish, Tables of angular spheroidal wavefunctions Vol. 1-8, Naval Res. Lab. Reports, June 30, 1975; Washington D.C.[20] Wolfram Research Inc., www.specialfunctions.com (2002).Appendix A—Legendre and spherical Bessel functionsIn this appendix we give the definitions of the Legendre and spherical Bessel functions used to thedefine the spheroidal wave functions. The definitions which we use here differ slightly from thosefound in [10], and are based on the approach taken in [20]. In particular, we define two Types ofLegendre function: the functions of Type I, PnmHzL and QnmHzL, are equal to those conventionallyused on the interval zœH-1,1L; the functions of Type II, /GothicCapPnmHzL and /GothicCapQnmHzL, are equal to thoseconventionally used for z–H-1,1L. The practice of using gothic characters /GothicCapP, /GothicCapQ to denote thefunctions of Type II follows Meixner and Schäfke [1].Legendre functionsThe Legendre functions of the first kind are defined by:(A1)PnmHzL=1ÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅGH1-mL H1+zLmê2ÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅH1-zLmê22F1J-n,n+1;1-m;1-zÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅ2N(Type I),(A2)/GothicCapPnmHzL=1ÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅGH1-mL Hz+1Lmê2ÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅHz-1Lmê22F1J-n,n+1;1-m;1-zÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅ2N(Type II),where GHzL is the gamma function ([10], Ch.6) and 2F1Ha,b;c;zL is the Gaussian hypergeometricfunction ([10], Ch.15). Note that these two definitions differ only in their phase, and are triviallyrelated:(A3)/GothicCapPnmHzL=H1-zLmê2ÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅHz-1Lmê2PnmHzL.18 For noninteger m, the functions of the second kind are defined by(A4)QnmHzL=pcscHmpLÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅ2IcosHmpLPnmHzL-Hn-m+1L2mPn-mHzLM(Type I),(A5)/GothicCapQnmHzL=pcscHmpLÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅ2 eimpI/GothicCapPnmHzL-Hn-m+1L2m/GothicCapPn-mHzLM(Type II).The Type I and II functions of the second kind are related by:(A6)/GothicCapQnmHzL=eimpHz-1Lmê2ÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅH1-zLmê2 JQnmHzL+pcscHmpLÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅ2JH1-zLmÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅHz-1Lm-cosHmpLNPnmHzLN,m–/DoubleCapZ,(A7)/GothicCapQnmHzL=H-1LmHz-1Lmê2ÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅH1-zLmê2ikjjjjQnmHzL+pÅÅÅÅÅ2è!!!!!!!!!!1-zÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅè!!!!!!!!!!z-1PnmHzLy{zzzz,mœ/DoubleCapZ.The Type II functions have the following important hypergeometric representation:(A8)/GothicCapQnmHzL=2-n-1eimpè!!!pGHm+n+1LÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅGHn+3ê2Lµz-m-n-1Hz+1Lmê2 Hz-1Lmê22F1Jm+n+1ÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅ2,m+nÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅ2+1;n+3ÅÅÅÅÅ2;1ÅÅÅÅÅÅÅÅz2N.Because of the factor GHm+n+1L, this expansion diverges when m+n=-1,-2,… and hence(A9)†QnmHzL§,†/GothicCapQnmHzL§Ø¶,form+n=-1,-2,…The following relations for nØ-n-1 are used in Section 2.2:(A10)Q-n-1mHzL=cscHpHm-nLLHpcosHmpLcosHnpLPnmHzL-sinHHm+nLpLQnmHzLL,(A11)/GothicCapQ-n-1mHzL=cscHpHm-nLLHp‰ÂmpcosHnpL/GothicCapPnmHzL-sinHHm+nLpL/GothicCapQnmHzLL.In the Mathematica system, the Legendre functions of Type I and II are implemented as “type 2”and “type 3” functions (“type 1” is a redundant variant of “type 2” which is only defined for†z§§1):PnmHzLõLegendreP@n,m,2,zDQnmHzLõLegendreQ@n,m,2,zD/GothicCapPnm HzLõLegendreP@n,m,3,zD/GothicCapQnm HzLõLegendreQ@n,m,3,zDSpherical Bessel functionsThe spherical Bessel function of the first kind may be defined by19 (A12)jnHzL=è!!!pÅÅÅÅÅÅÅÅÅÅÅÅ2 „k=0¶H-1LkÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅGHk+n+3ê2L k! JzÅÅÅÅÅ2N2 k+n.The function of the second kind is defined by(A13)ynHzL=-secHnpLHsinHnpLjnHzL+j-n-1HzLL.These functions are not implemented directly in the Mathematica system, so we compute themusing their relation to the regular Bessel functions JnHzL and YnHzL ([10], Ch. 9):jnHzLõè!!!!!!!!!!pê2 BesselJ@n+1ê2,zDëè!!!zynHzLõè!!!!!!!!!!pê2 BesselY@n+1ê2,zDëè!!!zAppendix B—Flammer’s spheroidal functionsIn this appendix we give the essential relations between the spheroidal functions of Flammer [2]and those of Meixner [1]. Note that Flammer’s functions are only defined for integer parametersm,n with n¥m¥0.Eigenvalues:(B1)lmnHcL=lnmHcL+c2.Angular functions:(B2)SmnHc,hL=wmnHcLpsnmHh;cL,(B3)SmnH2LHc,hL=wmnHcL qsnmHh;cL,where(B4)wmnHcL=loooooooomnooooooooH-1Ln-mÅÅÅÅÅÅÅÅÅÅÅÅ2 Hm+nL!ÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅ2n Hn-mÅÅÅÅÅÅÅÅÅÅÅ2L! Hm+nÅÅÅÅÅÅÅÅÅÅÅ2L! 1ÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅpsnmH0;cL,n-meven,H-1Ln-m-1ÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅ2 Hm+n+1L!ÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅ2n Hn-m-1ÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅ2L! Hm+n+1ÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅ2L! 1ÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅpsnm£H0;cL,n-modd.Radial functions:(B5)RmnH1,2LHc,xL=SnmH1,2LHx;cL. 20 Appendix C—Symmetry relations for the spheroidal functionsThe spheroidal wave functions satisfy a number of useful mathematical identities, which they“inherit” from properties the Legendre and spherical Bessel functions. Some of these wereoriginally derived in a number of papers by Meixner (almost all of whose work was published inGerman—see [1] and references therein), but do not appear in the most popular references onspheroidal wave functions (e.g. [2, 10]), and are therefore effectively unavailable to the majority ofcontemporary readers. The situation is further hampered by the different notations andnormalizations in existence: deriving identities valid for Flammer’s functions from those inMeixner’s is nontrivial. In this appendix we therefore present a concise summary of these identitiesfor the spheroidal wave functions, including all special cases which are not trivially obtained fromthe general ones. All relations are valid throughout the complex plane, and in particular alongbranch cuts. A more detailed discussion of the derivation of these identities can be found in [12].The transformation nÆ-n-1Eigenvalues:(C1)l-n-1mHgL=lnmHgL.Radial normalization factor:(C2)A-n-1mHgL=AnmHgL.Angular functions, general m:(C3)f-n-1mHz;gL=fnmHz;gL,f=Ps,ps,(C4)qs-n-1mHz;gL=cscHHm-nLpLHpcosHmpL cosHnpLpsnmHz;gL-sinHHm+nLpL qsnmHz;gLL,(C5)Qs-n-1mHz;gL=cscHHm-nLpLHpeimpcosHnpLPsnmHz;gL-sinHHm+nLpLQsnmHz;gLL.Angular functions of the second kind, integer m:(C6)qs-n-1mHz;gL=qsnmHz;gL-pcotHnpLpsnmHz;gL,(C7)Qs-n-1mHz;gL=QsnmHz;gL-pcotHnpLPsnmHz;gL.Radial functions, general n:(C8)S-n-1mH1LHz;gL=-sinHnpLSnmH1LHz;gL-cosHnpLSnmH2LHz;gL,(C9)S-n-1mH2LHz;gL=cosHnpLSnmH1LHz;gL-sinHnpLSnmH2LHz;gL.21 Radial functions, integer n:(C10)S-n-1mH1LHz;gL=H-1Ln+1SnmH2LHz;gL,(C11)S-n-1mH2LHz;gL=H-1LnSnmH1LHz;gL.The transformation mÆ-mEigenvalues:(C12)ln-mHgL=lnmHgL.Angular functions, general m:(C13)psn-mHz;gL=GHn-m+1LÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅGHn+m+1LJcosHmpLpsnmHz;gL-2ÅÅÅÅÅpsinHmpLqsnmHz;gLN,(C14)qsn-mHz;gL=GHn-m+1LÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅGHn+m+1LJcosHmpLqsnmHz;gL+pÅÅÅÅÅ2sinHmpLpsnmHz;gLN,(C15)Psn-mHz;gL=GHn-m+1LÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅGHn+m+1LJPsnmHz;gL-2ÅÅÅÅÅpe-impsinHmpLQsnmHz;gLN,(C16)Qsn-mHz;gL=GHn-m+1LÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅGHn+m+1L e-2 impQsnmHz;gL.Angular functions, integer m:(C17)fn-mHz;gL=H-1Lm GHn-m+1LÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅGHn+m+1LfnmHz;gL,f=ps,qs,(C18)fn-mHz;gL=GHn-m+1LÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅGHn+m+1LfnmHz;gL,f=Ps,Qs.Radial functions:(C19)Sn-mHkLHz;gL=SnmHkLHz;gL,k=1,2,3,4.The transformation zÆ-zAngular functions of Type I, general m,n:(C20)psnmH-z;gL=cosHHm+nLpL psnmHz;gL-2ÅÅÅÅÅpsinHHm+nLpL qsnmHz;gL,(C21)qsnmH-z;gL=-cosHHm+nLpL qsnmHz;gL-pÅÅÅÅÅ2sinHHm+nLpL psnmHz;gL.Angular functions of Type II, general m,n and z–H-1,1L:22 (C22)PsnmH-z;gL=expJpn"########-z2ízN PsnmHz;gL-2ÅÅÅÅÅp e-imp sinHpHm+nLL QsnmHz;gL,(C23)QsnmH-z;gL=-expJ-pn"########-z2ízN QsnmHz;gL.Angular functions, integer m,n:(C24)psnmH-z;gL=H-1Lm+n psnmHz;gL,(C25)qsnmH-z;gL=H-1Lm+n+1 qsnmHz;gL,(C26)PsnmH-z;gL=H-1Ln PsnmHz;gL,z–H-1,1L,(C27)QsnmH-z;gL=H-1Ln+1 QsnmHz;gL,z–H-1,1L.Radial functions, general n:(C28)SnmH1LH-z;gL=H-gzLnHgzL-n SnmH1LHz;gL,(C29)SnmH2LH-z;gL=-H-gzL-nHgzLnHSnmH2LHz;gL+H1+H-gzL2nHgzL-2nL tanHnpL SnmH1LHz;gLL.Radial functions, integer n:(C30)SnmH1LH-z;gL=H-1LnSnmH1LHz;gL,(C31)SnmH2LH-z;gL=H-1Ln+1 SnmH2LHz;gL.Relations between angular functions of Type I and II(C32)PsnmHz;gL=H1-zLmê2ÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅHz-1Lmê2psnmHz;gL,(C33)QsnmHz;gL=eimpHz-1Lmê2ÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅH1-zLmê2 JqsnmHz;gL+pcscHmpLÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅ2JH1-zLmÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅHz-1Lm-cosHmpLNpsnmHz;gLN.The second of these must be treated specially for integer m:(C34)QsnmHz;gL=H-1LmHz-1Lmê2ÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅH1-zLmê2ikjjjjqsnmHz;gL+pÅÅÅÅÅ2è!!!!!!!!!!1-zÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅè!!!!!!!!!!z-1psnmHz;gLy{zzzz.Relations between angular and radial functionsIn this section we present the set of relations between the angular and radial functions which weuse to compute the functions throughout the complex z-plane. They can be derived using Eq.(19)and the identities given in this appendix. The forms presented here, which are completely generaland valid along branch cuts, have not previously appeared in the literature.23 Radial functions, general m,n:(C35)SnmH1LHz;gL=KnmHgL sinHHm-nLpLÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅp e-iHm+nLp H1-1êz2Lmê2 HgzLnÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅgnzn-m Hz-1Lmê2 Hz+1Lmê2 Qs-n-1mHz;gL,(C36)SnmH2LHz;gL=secHnpLIS-n-1mH1LHz;gL-sinHnpLSnmH1LHz;gLM.Radial functions, integer m,n:(C37)SnmH1LHz;gL=KnmHgL H1-1êz2Lmê2 zmÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅHz-1Lmê2 Hz+1Lmê2 PsnmHz;gL,(C38)SnmH2LHz;gL=H-1Lm+1ÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅgKn-mHgLAnmHgLAn-mHgL H1-1êz2Lmê2 zmÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅHz-1Lmê2 Hz+1Lmê2 QsnmHz;gL.Type II angular functions, general m,n:(C39)QsnmHz;gL=pcscHHm+nLpLeiHm+nLpÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅK-n-1mHgL Hz-1Lmê2 Hz+1Lmê2 HgzLn+1ÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅH1-1êz2Lmê2 gn+1 zm+n+1 S-n-1mH1LHz;gL,(C40)PsnmHz;gL=secHnpLÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅp ‰-ÂmpHsinHpHm+nLL QsnmHz;gL-sinHHn-mLpLQs-n-1mHz;gLL.Type II angular functions, integer m,n:(C41)PsnmHz;gL=1ÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅKnmHgL H1-1êz2Lmê2 zmÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅHz-1Lmê2 Hz+1Lmê2 SnmH1LHz;gL,(C42)QsnmHz;gL=H-1Lm+1 gKn-mHgLAnmHgLAn-mHgL Hz-1Lmê2 Hz+1Lmê2ÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅH1-1êz2Lmê2 zm SnmH2LHz;gL.Relations involving the Type I angular functions, psnmHz;gL and qsnmHz;gL, can be obtained from thelast four equations by using Eqs.(C32-34).Appendix D—Tables of numerical valuesIn this appendix we present a set of numerical values, to 25 digits of precision, for all of Meixner’sspheroidal functions. The intention is to provide enough values and to high enough precision tofacilitate comparison with any future software implementation. Tables 1 and 2 contain theeigenvalues lnmHgL for integer and complex parameters respectively; Table 3 contains the joiningand normalization factors KnmHgL and AnmHgL; finally, Tables 4 and 5 contain the angular and radialfunctions and their derivatives. Values for the Type II angular functions, as well as all ofFlammer’s functions, are not given, since they can easily be found from those presented here.24 Table 1 Eigenvalues lnmHgL for integer m,n and real g2.mnglnmHgL lnmHigL01011210100101001010010100-90.7716957027500548489877312-9900.7518988910167474495421523-71.8665362671732721853810250-9701.7595433440823666225640610-89.7122312326085318292420084-9899.7468223865850616234724355-70.6610819583855185299419784-9700.744156595858817379153742618.9720560550422438139109191198.997474634082548135724810318.9720619762544159268471575198.997474634082548135724810337.8806498956194532262871049397.989846793913121459744012537.8808487977730112048164244397.9898467939131214597440125Table 2 Eigenvalues lnmHgL for complex m,n,g. Here we use the abbrevation a=1+i.mncReHlnmHgLL ImHlnmHgLL00a10 a0.05947276973503126247061569.2407662146346033515957443-1.3371748778053999710372379-189.98934859565755367515086960a10 a10.50186776246700453075162679.50003165123420466677881692.9507369925182112070617898209.9992573181593545006418858a10 a01-0.9078192346934944943133571-13.78249204145363996320697930.937476128194795842364958017.0373891416686511344181798a10 aa10 aa10 a1.146173558736254250502993213.77544665374287955398693001.331825843494567670634608314.1334443105191566448899153Table 3 Joining factor, KnmHgL, and radial normalization factor AnmHgLmngKnmHgL AnmHgL0101121010 i1010 i1010 i1010 i428.0069932832745874608231771129.955591534821883488517299989.1933715984633653549677126225.0893949188277187831526217i494.806600890673485535027288427.4621485726353680526646286i43.8069890817468723895524520-20.46826393773882905641784460.00092599590016865734973774.35228568796845942426840860.00444351505859583160084892.5127949340421379580116552-0.00008456655430912744771420.5598962979485962334204584-0.00107628476364161970540260.7511925147800865086125805Table 4 Type I angular functions and their derivatives at z=0.mngpsnm H0;gLpsn+1m£ H0;gL01011010i1010i1.8695013198832203237866070µ1008.1392106153914773135592685µ10-4-1.5290337582543180975733869µ100-4.1071723604572527466632257µ10-34.6221868979445343185957783µ1004.2001780506231961222071385µ10-3-8.8274907181871032109649776µ100-4.3315286911297506025068055µ10-2mngqsn+1mH0;gLqsnm£H0;gL01011010i1010i-4.2717498257693456557494192µ10-6-1.5033025515694459977003079µ103-9.2725118702472982516514180µ10-64.9349301484865713534797135µ1024.5866156819976162315752064µ10-72.3273007180665026590360560µ1044.1959649830139821978226061µ10-7-2.8912768871677188308743378µ10325 Table 5 Radial functions and their derivatives at z=1.005.mngSnmH1L Hz;gL SnmH1L£Hz;gL2222223312346.6119132248515374422725009µ10-42.5659296586989964008140566µ10-32.2065345978824180503885691µ10-34.6827642681955017561952436µ10-31.3247288100076832070527852µ10-15.1297872006118942981483008µ10-14.4231954640285939420530600µ10-19.3475721512114037868171462µ10-1mngSnmH2L Hz;gL SnmH2L£Hz;gL222222331234-3.7497722396542435481278539µ102-4.8522267972282203610936955µ101-3.7428718891971076782275646µ101-1.3339979013106281309007387µ1017.5736490437910731355302702µ1049.7369858589493594357303506µ1037.5660512493589672475730118µ1032.6625329643356096410107459µ103 26