conical functions app paper

Methods for calculating electrostatic quantities due to a free charge in a nanoscale three-dimensional tip/base junction Vallorie J. Peridie?) Department of Mechanical Engineering, Temple Universiry Philadelphia, Pennsylvania 19122 Li-Hong Pan and Thomas E. Sullivan Depament of Electrical Engineering, Temple University Philadelphia, Pennsylvania 19122 (Received 27 October 1994; accepted for publication 5 July 1995) The exact solutions for fully three-dimensional electrostatic quantities due to a free charge in a nanoscale tip/base junction, such as the scanning-tunneling microscope, are given for the problem as modeled in the prolate-spheroidal coordinate system. These exact solutions consist of summation series of integrals of conical functions, and methods for calculating these exact solutions with high accuracy are described in detail in this paper. Calculated results for the surface-charge distribution on the tip and the base due to a free electron in a nanoscale tip/base junction are also presented. This analysis has important implications for the ultimate objective of quantifying electron tunneling, from first principles, in nonplanar geometries such as sharp vacuum field emitters. 0 1995 American Institute of Physics. I. INTRODUCTION Fully three-dimensional nanoscale microelectronic struc- tures have significant and practical applicability in the areas of very high-speed switching devices, intense electron/ion sources, and advanced image-display devices (e.g., HDTV), among other applications.‘T2 Consequently, it is desirable to develop the capability to calculate the tunneling current of three-dimensional nanoscale electronic structures from first principles. However, devices such as these, utilizing atomi- cally sharp field-emission tips, are inherently nonplanar ge- ometries and they exhibit tunnel currents which are substan- tially larger than predicted by planar or one-dimensional models.3 An important first step in extending the analysis of elec- tron tunneling beyond the often-used one-dimensional analy- sis of the Fowler-Nordheim field-emission theory4 is the re- cent result reported by Pan et al. (1994)5 for an exact expression, in the prolate-spheroidal coordinate system, of the fully three-dimensional potential-energy barrier associ- ated with emitter tips. The consequent detailed understanding of the three-dimensional electrostatic potential and related quantities is an essential theoretical element for the ultimate objective of quantifying axial (symmetric) and off-axis tun- nel currents. The exact solution for several quantities related to the electrostatic interaction of a free electron in a biased tip/base junction, modeled in the prolate-spheroidal coordinate sys- tem, has been derived by two independent methods.5’6 The three-dimensional electrostatic potential V due to a free charge in a tip/base junction was derived using a classical separation-of-variables approach in Pan et aL5 where com- puted results for this quantity are given. The Green’s func- tion for the region between the tip and the base, modeled in the prolate-spheroidal coordinate system, is derived in Peri- dier et al.;6 this result might enable one, for example, to “Electronic mail: [email protected] model electron tunneling in a region with a charge distribu- tion such as that which would occur in the situation of tun- neling through physical medium as opposed to vacuum. This second paper6 also includes calculations for the potential- energy barrier @ at on-axis and several off-axis sites in the tip/base geometry. Both analyses account properly for the effects of the image and multiple-image interactions that would occur in the limiting cases when the tip/anode spacing is O(fO0 A), such as in the situation of scanning tunneling microscopy (STM). The derived expressions for every electrostatic quantity in both papers5Y6 include a summation series of integrals of conical functions K’J, and while computed results were given in each paper the numerical methods were not described. The intention of this paper is to describe the specific numerical methods required to accurately calculate the electrostatic quantities derived in Refs. 5 and 6. In addition, computed results for the surface-charge distributions, induced on both the tip and base surfaces by a free electron in the tip/base junction, are presented here. II. PRIOR STUDIES He et al7 recognized that the potential-energy barrier for electron tunneling in a tip/base junction might vary consid- erably from a potential-energy barrier determined from pla- nar assumptions. He et al7 obtained analytic expressions for the induced image potential due to a free charge on the axis of symmetry of an STM tip in free space, modeled variously as a cone, hyperboloid, paraboloid, and as a sphere on the tip of a cone. However, it is worthwhile to note that this afore- mentioned analysis7 is fundamentally a one-dimensional study because only free-charge locations on the symmetry axis of the tip/base junction were considered. Furthermore, the result of He et aL7 is somewhat difficult to interpret given that the induced surface charge was considered on the tip surface only of the tip/base junction, but not on the base. (Physically, of course, the surface-charge distributions on the 4888 J. Appl. Phys. 78 (8), 15 October 1995 0021-8979/95/78(8)/4888/7/$6.00 Q 1995 American Institute of Physics hpws(0)-1 I- lP=elwJ~ ‘\ ‘\ ,/ . . . .. .. . . . . . . . . . . . , . . . . . . . ‘\ / + : ,v ); ‘.. -,. -.....,,,.,,,,,, /’ ..” \\ . . . . . . . . . . . . . f . . . . . . . . . . I . . . . J” . . I . . . . . ..’ . ..*..q..* . . . . -.. ,’ ‘... 5. *..*i ..’ i \, ,,,.. * .,I. I I.., . . . . I ‘.., 5 . . . ., /’ ‘*.. %.. c :“. . . .J i‘; f’ * . “,.. ,q ..,, $.g ~~$+~<~..~ J./y **+- */* &z,,-j +* : \ b, i’ ; : -f”-+~‘t-~~ \ ; . i ii i : :::::: i i i .l--?C ; i i/Tt : f : -:;--i i : :. *- : /+ ‘., L&+y~. : -- i i- : *- : cc-c ‘:: ..,. . . . . . i.$.( ,..y ‘x. .*I ‘L., I. . ..I... “, ,. -.;--..! : . . . . . . . . . . . *;” \ ‘... A . . \ / i .:’ --.* . . . N . . . \ “SC’ ‘.... ,’ ‘;” ‘“...,..,.,... . . . . +.’ I /’ \\ ..: 8 . . . . x., 8’ / $ ,/ I \ \ ‘\ a’ x. , b..... I \../” ‘\ t’ . . . . . . . . . r . . . . . . . . . . . ...‘\ \ ‘. ,’ ,’ ;~@J\\ FIG. 1. The prolate-spheroidal coijrdinate system. Hyperboloid surfaces ap- propriate for modeling a tip/base junction are shown in dark print. tip and the base are influenced both by: (i) the free electron in the junction, and (ii) by the surface-charge distribution on the opposing surface.) Ill. THE THREE-DIMENSIONAL PROLATE-SPHEROIDAL MODEL OF A TIP/BASE JUNCTION In our model the junction of the tip and the base is mod- eled as a region bounded by two hyperboloids in the prolate- spheroidal coordinate system, which is given by (Ref. 8, p. 103-104) x=a sinh(u)sin(u)cos($), y=a sinh(lc)sin(u)sin(4), (1) z=a cosh(lc)cos(u). In this coordinate system a is one-half the distance between the foci for each hyperboloid surface; curves corresponding to constant u are hyperboloid surfaces, curves corresponding to constant u are prolate-spheroidal surfaces (see Fig. 1). The variable transformation a$= cash(u). 7,7=cos(u), generate a coordinate system given by x=a@TJm cos(qb), y=a@TJw sin(+), (2) z=atq, with the corresponding coordinate metrics: (3) FIG. 2. The geometry for the model problem h,=&2- v2M1 - v2>, h,=aJ(t2- l)( 1- v2). (4) To model a tip/base junction in this prolate-spheroidal coordinate system two level hyperboloid surfaces of 7 are taken, one to represent the tip (say, at v1=0.9) and the sec- ond to represent the base. (A flat anode would be represented by v2=0.) In the new coordinate system, the range of the variables which define the junction or region between these two hyperboloid surfaces is 8E[7h77,721. (5) This tip/base model in the prolate-spheroidal geometry is shown in Fig. 2, with lines of constant t given for ~={1.0(0.1)1.5}. It is worthwhile to note that the radius of curvature rl of the tip modeled in this fashion is rl =a sin( 8t)tan( et), (6) while the cone angle of the tip 8, (defined as one-half of the angle subtended by the asymptotes of the hyperboloid sur- face) may be determined from 8, =cos-‘( 7,). (7) Furthermore, in this model the separation distance dsep of the tip from the flat surface is dsep=a 771. (8) Consequently, it is worthwhile to note that in this prolate- spheroidal model only two among the three geometric pa- rameters radius of curvature (ri), tip cone angle (e,), and separation distance (dsep), are independent, because fixing two of these three parameters automatically determines the third. J. Appl. Phys., Vol. 78, No. 8, 15 October 1995 Peridier, Pan, and Sullivan 4889 If a charge of strength 4 is placed at an arbitrary location x, in a voltage-biased tip/base junction, the solution for the electrostatic potential V may be written as the sum of two subsidiary potentials v=v,+v,, (9) where V, is the potential due to the bias voltage V, applied across the junction, and Vc incorporates the effects of the free charge’s interaction with the surface charges on the tip and base surfaces. The solutions for V, , Vc are derived in Pan et al5 and are ln{[t 1 + v2M 1 - 77dlH I- d/t I+ d]} “= ‘O in-X( 1+ aaM 1 - v2)][( 1 - vi)/( 1 + 77i)]} ’ vc= 4 47reolx-x,I +v,, with V, being the interaction potential given by x I om[B:m 7) + c:a - 77)lcY5)d~. In Eq. (12) the coefficients BT ,CF are given by By= x K3%)Kr=(- 1?1)-K3- %)K’=(QQ) i i cxmm- 772)-a- 771K3772) ’ CY= & q~~(Ecw'=(771)1 x i cc - %KY 772) -q ?IJK!J( - 772) 1 K~(77*)K’=(-772)-K~(-711)Kr=(~);/2) . (10) (11) (12) (13) (14) For expressions associated with Vc, the charge is located (without loss of generality) at x,=( 5, ,v, ,0), and the coef- ficient term fl in expressions (13) and (14) is: (1% Here, F is the gamma function, and the method for calculat- ing this quantity is given in Appendix A. The expression KF in expressions (12), (13), and (14) is the “conical function” defined by KT(P)sP:++iT(P); (16) conical functions are a special case of generalized Legendre functions which have an integer order m and a continuum of complex degrees of the form (--ifiT), for r~[O,a). The second solution for KY(p) is Ky(-p). Conical functions arise in the solution of the Laplace equation in the prolate- spheroidal coordinate system, and the numerical evaluation of these functions is described in Appendix B. The potential-energy barrier @ for a charge tunneling from the tip to the base, located at x, is derived in Pan et aL5 and is given by @=q(Vg+V,). (17) Finally, the induced surface-charge distributions ‘T, on the tip, and a2 on the base, may be determined from (18) where the metric h ‘I is given in Eq. (4). Taking the derivative of Eq. (9) with respect to 7, evaluating the result at the tip p and using the Wronskian identity” WK%LK:(-~11 (- 1)” 2 cosh(rrr) r($+ir+m) =(l-$) 7r lY($+i7-m) gives (T, , the tip surface-charge distribution: (19) cc -2EoVo u’=~~~~~~ln[(l+~2)~(l-~2~(l-~,)~~l+~~)~+~ ~~~~~‘~~o’-~)~‘~-‘~.“) m X cos( m 4) 1 i d7 r tanh( ~7) r(++iT-m) [K:t&)K:t5)1 0 r($+iT+m) ( K:( - rlc)K:( ~2) - KI=( v,)K)=t - 72) KI=(n)K:(- 772)-K3 - v,)K:(vd )1 . cm Similarly, the surface-charge distribution a2 on the anode surface at 772 is m -2q)vo (+2=- a ~~~~ln~(l+‘li),(l-~~~~l-~~),(~+~~),+~ ~~J~~‘~~o(-~)~‘~-~~~~) 1 I omdr r tanh( nr) r(i+iT-m) X cos( m q5) r(i+iT+m) [K:(tXT(S)I K:t- ‘idK%+K%dK’=(- VI) Kt’(m)K:(- 712)-K3- 77JK%72) 11 ’ (21) 4890 J. Appl. Phys., Vol. 78, No. a, 15 October 1995 Peridier, Pan, and Sullivan In expressions (20) and (21) it is evident that the charge distributions on each surface depend on the applied bias- voltage V,, and on the free charge in the junction at &=( 5, , vc ,O), and on the opposing surface. The above expressions for the electrostatic field V [Eq. (9)], the electron-tunneling potential barrier Q, [Eq. (17)], and the surface-charge distributions on the tip and the base uI ,a2 [Eqs. (20) and (21)] are exact. Computed results for V are given in Pan et al.,’ and computed results for @ are given in Peridier et aL6 where the more generalized Green’s- function result is described. Methods for calculating these quantities, as well as computed results for the surface-charge distributions o, ,a2 induced by a free electron in the tip/base junction, are described below. IV. NUMERICAL METHOD FOR CALCULATING ELECTROSTATIC ,QUANTlTlES IN THE PROLATE- SPHEROIDAL COORDINATE SYSTEM The expressions above for the electrostatic potential V [Eq. (9)3, for th e t unneling potential-energy barrier @ [Eq. (17)], and for the induced surface-charge distributions o1 ,a2 [Eqs. (20) and (21) are exact, and can in principle be calcu- lated to arbitrary accuracy. However, each expression in- cludes a term, similar in form for each expression, that is a summation series (over m =O,l,Z,...,) of integrals (taken over r~[O,w)) of gamma functions r of complex arguments mul- tiplied by conical functions KT. As a practical matter, the numerical evaluation of this term is not trivial and thus will be described in this section, using the expression V, [Eq. (12)] as an example. First, experience has demonstrated that, for evaluating an expression such as V, [Eq. (12)], the result will converge faster if the summation operator is interchanged with the integral, i.e., V 4 m -- r.cdc-4 Eoa s d+( 4, 0 (22) with the argument of Eq. (22) given by the summation series N (23) and each term s,(r) in the summation series F(r) given by [see Eqs. (12), (13), (14), and (15)] The above expression for V,,calc-+VI in the limit N+w, where N is the number of terms taken in Eq. (23). Note that the field point is x=(&v,@) and the location of the free charge is x,=( .$, vc ,O) and that both are fixed constants; the dummy variable of integration is 7: lated for m=0,1,2,...,30 and then, in the numerical proce- dure, each term s,(r) was compared to the s,-,(i) and s,+ i (7) terms. The numerical behavior of these terms, illus- trated schematically in Fig. 3, was as follows. For small or- To calculate V, from Eq. (22), Vl,calc, F( 7) was inte- grated using Simpson’s Rule in steps of Ar=0.05; (25) the integration was carried out over successive intervals in rs[O,l], 7~ [ 1,2], rE [2,3] ,..., until the magnitude of the last interval calculated was less than lop9 of the magnitude of the total integrated result for Vl,calc. In practice the integra- tion was typically terminated for r--0(100). To calculate F(T) at a specific value of r the summation series given by expression (23) was summed for Nc30 terms; the mecha- nism for selecting N is described below. Each term ~~(7) of expression (23) is the product of gamma functions F and conical function KY; the method used to calculate the gamma function F is described in Appendix A, and proce- dures for calculating the conical function KY is discussed in Appendix B. S&l i . 1.5 . . 1.0 -1.0 t . . The number of terms N taken in the series to compute the integral-argument F(r), defined by Eq. (23), was deter- mined as follows. The summation terms sm(r) were calcu- FIG. 3. Schematic showing the numerical character of subsequent terms S,,,(T) for increasing order m. In the schematic example illustrated here, the series solution for F(T), ECq. (22). would be taken to N=25 terms. J. Appl. Phys., Vol. 78, No. a, 15 October 1995 Peridier, Pan, and Sullivan I 4891 der m the value of each term S,(T) was finite, with the mag- nitude of this term ultimately falling off for larger orders of m. However, for some values of x and x, , (particularly for x very close to x,, or x, very close to a surface) after some term N subsequent terms of s,(r) alternated in sign and grew exponentially in magnitude. Consequently, each time F(T) was calculated n was chosen to be the largest value of m prior to the onset of the exponential growth (or to 30 if no such growth occurred). Note that taking m>30 results in a floating-point-overflow error when calculating the gamma function so it was convenient to terminate the series for a choice of NG30. This numerical approach guarantees three to five significant digits of accuracy which is satisfactory for our purposes. The method for calculating the electrostatic quantities that is described above and in Appendices A and B is esti- mated to yield three to five significant digits, with accuracy increasing as the distance of the field-point x from either the charge-location x, or the boundary surfaces is increased. This estimate was obtained by numerical study of an exact expres- sion given in van Nostrand1o for l/h-l as integrals over associated conical functions: &=- z mlo (2- ~O,mb(m+) x I md~4ww35c) 0 an expansion of CY 77)fcY - %I 77,< 77 ~~(-~7)~IS(%) VP?1 . (26) The form of this expression is quite similar to that of the electrostatic quantities described above and, given the fact that l/lrl can be evaluated exactly, it was possible to directly evaluate the accuracy of the numerical-integration scheme by calculating expression (26) using the method described above and comparing the result to the exact solution. V. CALCULATED RESULTS FOR THE SURFACE-CHARGE DISTRIBUTIONS aI ,a2 Equation (20) was used to calculate the surface-charge distribution on the surface of the tip ‘+, , and Eq. (21) for the charge distribution a, on the surface of the anode. Figure 4 illustrates both distributions for the case of zero applied volt- age V. and thus the two surface-charge distributions are in- duced solely by the free electron. The tip in Fig. 4 is modeled as a level hyperboloid surface 0, =0.9 (which corresponds to a tip half angle of about 26”) and the charge is located at a point x,=(5,, vc ,+c) =(13, 5, 0) in the tip/base junction. In Fig. 4 the charge is located in the plane of the paper at the dot, the topological surfaces of the tip and base are repre- sented by the two inner grids, and the corresponding surface- charge distributions are represented by the topmost (c,) and lowest (a& grid surfaces, respectively. (All surfaces are symmetric in the plane of the paper.) Figure 5 provides more detail of the surface-charge distribution on the tip. In both figures it is worthwhile to note that the max- imum (minimum) values of oi (~~2) occur at the point in 4892 J. Appl. Phys., Vol. 78, No. 8, 15 October 1995 FIG. 4. The surface-charge distributions ((T,,u~) on the tip, base surfaces (v,,n). The charge is located at ~~=(~~,~,,~~)=(1.3,0.5,0.0), and the cone angle (for 7, =0.9) is about 25.8”. the surface closest to the charge. The maximum value on the tip is o, =3.523(4/a’), and the minimum value on the base is a,=-2.215(q/a2), where the electron charge q=-1.6x10-‘9 C, and the length-scale a is determined from Eq. (6). For example, if the tip radius of curvature is 20 8, and the tip cone angle is g,=O.9, then a-94.7 A, and the separation distance of the tip and anode is d,,,-85.3 A. FIG. 5. The surface-charge distribution (u,) on the tip (detail of Fig. 4). Peridier, Pan, and Sullivan VI. SUMMARY Exact fully three-dimensional solutions for the electro- static quantities V, a’, and g, ,~z for an electron in a tip/base junction are given in this paper in the prolate-spheroidal co- ordinate system. Very accurate numerical methods for calcu- lating these quantities are described in detail. This result is being used to investigate quantum-size effects in nanoscale electronic devices. APPENDIX A: CALCULATING THE GAMMA FUNCTION l-‘(p) A strategy described in Atfken (Ref. 8, p. 556) was adapted to create a very accurate scheme for calculating the gamma function T(p) for complex argument p. Calculating the gamma function from its series definition can produce large round-off errors; however, if the argument ,U is very large a very accurate result can be obtained by using Stirling’s asymptotic approximation. In this study a variation (Ref. 11, p. 257, Eq. 6.1.40) of Stirling’s asymptotic series, ln(27r) In[T(z)]= 2 + i 1 z -k In(z)-z-t 2 B 211 +**- (2n)(2n- l)(zz”-1) +.** ’ (Al) was used to calculate the gamma function for large z. In Eq. (A I) the coefficients B2,, are Bernoulli numbers, tabulated in Ref. 11 (p. 810, Table 23.2), and the sum was carried out to the z-Z3 term, inclusive. In this study, if lp[>20, then the gamma function was computed directly from Eq. (Al) by Upu)=ew(W(~u)l). 642) On the other hand, if /,~[~20, an intermediate variable z” was created by i=fl-into(L)+20, (A3) where the function int(p) returns the largest integer that is less than or equal to ,u, generating a variable i>20. The gamma function of ,$ was obtained from I72)=exp(lnCUz^)l), (A4) and then I’(p) was obtained from T(t) by recursive division, viz.: r(p)= r(f) (p)(p+ l)(p+2)***(2- 1)’ (.45) This calculation was carried out in FORTRAN using the COMPLEX* 16 variable data types. APPENDIX B: CALCULATING THE ASSOCIATED CONICAL FUNCTIONS K:(p) This appendix describes how to calculate associated conical function KY(p), which are defined by @I) thus the associated conical function Ky(,u) is a special case of the generalized Legendre function P,“(p) with a complex degree n=( - i+ ir). In generalized Legendre functions the order m must be real, and the degree n and the argument p may be complex; but for the problem considered in this study the order m is always integer and the argument ,X always real, with ---cO<~<~. A comprehensive text on Legendre functions is Hobson,’ which gives numerous series representations for P,“(p), and the appropriate choice depends principally on the value of the argument ,u. For example, in the range - I<@ 1 the associated conical function has an exponential character with ‘mf+iT (CL)* 1 as ~31, m=O; Pm++iT(p)+O as P-1, mf0; ‘rIli+iT (Lc)+a as p-+--I, all m. 032) On the other hand, for values of I,u]> 1 the associated conical functions have an oscillatory behavior. Because the funda- mental character of the solution changes dramatically for I~131 + this region of the function argument is the most difficult to evaluate with a high degree of numerical preci- sion. After a fair degree of investigation and experimentation with various expressions, the following approach seemed to provide the greatest accuracy. First, the associated conical functions for order m> 1 were determined recursively using the following two expres- sions from Hobson (Ref. 9, p. 290, Eq. 166): Pz’2(p)= +(n-m)(n+m+ l)Py-2(m+ 1) Prf2(p)= -(n-m)(n+m+ l)Pr-%(m+ 1) x~&~P:“(,u) for ll.~lsl. 033) Because the conical functions are specifically of interest here, n = ( - i+ i r) and the recursion relations are K~‘2(p)=-[~+(m+~)2]K~(p)-2(m+ 1) Xd~K:“(p) for I,ul~l. (B4) Thus, it is very straightforward to obtain KY(p) for m> 1 once K:(p) and K:(g) are evaluated. To calculate K:(p) and K&d the following strategy was utilized. If l&1.1 the following equation (Ref. 9, p. 234, Eq. 68): J. Appl. Phys., Vol. 78, No. 8, 15 October 1995 Peridier, Pan, and Sullivan 4893 P,m(p)=2rn sin((n+m)n) I?(n+m+ 1) 2-c n+m+l) cos(nfl) r(n++)r(;) X(p2- l)m’2 .F(m+$n+n+1,n+&1/22) +2m r(n+$ ryn+ i -m)r(f) ,hq/& l)rn’2 +F(m+ i,m-n,)-n,llz2) 035) was used. In Eq. (B5) z = ,U + dm, F( ) is the hypergeo- metric function of the first kind, and since the conical func- tions are specifically of interest, n=( - 3+ ir). If the argument satisfied the condition that l.O<]p]< 1.1 the associated conical function was evaluated as follows. Start with the relation (Ref. 9, p. 235, Eq. 73) and note that, for integer m and \p]> 1, Pr , Pi” are related by (Ref. 9, p. 205) T(n+m+ 1) --m piY(P)= rtnmrn+ 1) pn (PI* (B7) If in Eq. (B6) the order j+ - m and the resulting expression is incorporated into Eq. (B7), we get p;(P)= r(ft+m+ 1) 2-m r(n-m+ 1) rtrn+ 1) (~2-1)jm’2Zn-m .F f+m,-n+m,l+2m,l-f . i 038) To calculate the conical function from Eq. (B8), again the degree n=( - $+ir). Finally, if the argument ,X satisfies the condition ],~]<l the following expression (Ref. 9, p. 227, Eq. 55) 1 I + p jl2 P!(p)=- - i 1 T(l-j) 1-p *F[-n,nSl,l-j,(l-~)/2] (B9) was utilized. On the crosscut of the Legendre function (p real, and I,uI< 1) PF , Pi” are related by (Ref. 9, p. 230, Eq. 61) I++m+ 1) f?‘tP)=t- 1Y’ r(n-m+ 1) P,“(p). @lo) So, letting j-1 -m in Eq. (40), and substitution of the result into Eq. (BlO) gives p:(P)= (-1)” T(n+m+l) l+p -14~ i 1 r(l+m) r(n-m+l) 1-p .F[--n,n+l,l+m,(l-,u)/2]. 0311) Once again, setting n= ( - $+ ir) results in the associated conical function K:(p). In this study these conical functions KY, the gamma functions I?, and the hypergeometric functions F, were all calculated using the COMPLEXITY data type in FORTRAN, with the gamma function calculated recursively as described in Appendix A, and the hypergeometric function evaluated using a nested, synthetic-division algorithm. ’ T. E. Sullivan, Y. Kuk, and P. H. Cutler, IEEE Trans. Electron Devices 31, 2659 (1989). ‘T. Utsumi, IEEE Trans. Electron Devices 38, 2276 (1991). 3C. A. Spindt, I. Brodie, L. Humphrey, and E. R. Westerberg, J. Appl. Phys. 47, 5248 (1976). 4R. H. Good and E. W. Miiller, Electron-Emission Gas Discharges I, Vol. XX1 of Encyclopedia of Physics (Springer, Berlin, 1956), pp. 176-231. 5L.-H. Pan, T. E. Sullivan, V. J. Peridier, P. H. Cutler, and N. M. Misk- ovsky, Appl. Phys. Lett. 65, 2151 (1994). 6V. J. Peridier, L.-H. Pan, T. E. Sullivan, and P. H. Cutler (unpublished). ‘3. He, P. H. Cutler, N. M. Miskovsky, T. E. Feuchtwang, T. E. Sullivan, and M. Chung, Surf. Sci. 246, 348 (1991). ‘G. Arfken, Mathematical Methods for Physicists (Academic, New York, 1970) pp. 103-104, 556. 9E. W. Hobson, Theory of Spherical and Ellipsoidal Harmonics (Chelsea, New York, 1965), 2nd reprint. ‘OR G van Nostrand, J. Math. Phys. 33, 276 (1954) (Note: this journal is now entitled: Stud. Appl. Math.). “M Abramovitz and I. Stegun, Handbook of Mathematical Functions (Do- ve, New York, 1964). 4894 J. Appt. Phys., Vol. 78, No. 8, 15 October 1995 Peridier, Pan, and Sullivan Journal ofApplied Physics iscopyrighted bytheAmerican Institute ofPhysics (AIP). Redistribution ofjournal material issubject totheAIP online journal license and/or AIP copyright. Formore information, see http://ojps.aip.orgfapofjapcrjsp Copyright ofJournal ofApplied Physics istheproperty ofAmerican Institute of Physics and itscontent may notbecopied oremailed tomultiple sites orposted toa listserv without thecopyright holder's express written permission. However, users may print, download, oremail articles forindividual use.