On the computation and applications of Bessel functions with pure imaginary indices (orders) in physics and specifically in corpuscular optics Текст научной статьи по специальности «Математика»

ISSN 0868-5886

НАУЧНОЕ ПРИБОРОСТРОЕНИЕ, 2014, том 24, № 1, c. 144-151

РАБОТЫ ШКОЛЫ ПРОФ. Ю.К. ГОЛИКОВА: -

РАБОТЫ, ПОСВЯЩЕННЫЕ ПАМЯТИ Ю.К. ГОЛИКОВА

УДК 517.6/8(083.3) © A. A. Matyshev, E. B. Fohtung

ON THE COMPUTATION AND APPLICATIONS OF BESSEL FUNCTIONS WITH PURE IMAGINARY INDICES (ORDERS) IN PHYSICS AND SPECIFICALLY IN CORPUSCULAR OPTICS

Bessel functions with pure imaginary index (order) play an important role in corpuscular optics where they govern the dynamics of charged particles in isotrajectory quadrupoles. Much more important role Bessel functions with pure imaginary index play in mathematical physics as solutions to Laplace's equation by method of variable separation in cylindrical coordinates. Also it was shown recently that such functions appear as solutions to Lame's equations which are used to characterize the displacement fields in semiconductor nanostructures. But there is no concrete algorithm that computes Bessel functions with pure imaginary order despite the fact that they were discovered long ago. In this paper an algorithm which can be used for the computation of the normal and modified Bessel functions with pure imaginary index is proposed. The developed algorithm is very fast to compute and for small arguments converges after a few iterations.

Keywords: Bessel functions with pure imaginary indices (orders), corpuscular optics

INTRODUCTION

Bessel functions occur in many branches of mathematical physics [1, 2, 12, 13] as solutions to differential equations when the Dirichlet and Neumann boundary conditions are applied on various space domains such as cylindrical or spherical. In a great variety of these applications (diffraction, electrical induction, etc.) only Bessel functions of the zeroth and first orders namely J0 (x) and J1 (x) occur. In other physical applications (such as solutions to Kepler's equation), the entire index n of Jn (x) are used. The Bessel functions of complex-valued indexes have had limited number of applications up till now. Presently (to the best of our knowledge) there is no existing mathematical software package or algorithm that computes these functions in their entire domain of definition as they have not been found applicable in areas of physics. The normal and modified Bessel function of pure Imaginary order was shown to be the solution governing the motion of corpuscular particles in isotrajectory focusing and deflecting systems [3]. Recently it was shown that they also appear as solutions to the Lame's equation which are used to characterize the displacement fields in semiconductor nanostructures. Usually the Bessel function of pure imaginary order with real or complex valued argument is defined with the aid of the complex Gamma function [4]. Obtaining the real values of these functions, however, is very vital in mathematical physics and applied sciences. The Bessel functions of order v can be represented as the solutions to the following differential equation [5]:

x2y" + xy' + (x2-v2)y = 0, (1)

where y' represent the first derivative with respect to x.

The resulted function can be expressed as an absolute converging series

Jv (x)=( 21 § (-^n^T) (2 J. (2)

Even though the above function has whole number values of v in the entire complex plane, it is practically non-trivial to compute this function for non real values of the index. Most computations and computational packages nowadays on this basis require that each element of the series is obtained via a complex Gamma function [4]. For instance for the case of a purely imaginary index iv and for a natural number n , we have

n

r(iv + n + 1) = r(iv)n(iv + m) . (3)

m=0

Hence to compute real solutions for equation (1) for purely imaginary order, it is required to sum a series with each coefficient being a complex number and not to forget the case of obtaining the complex value function r(iv), which on its own is very cumbersome to compute. Also to the best of our knowledge, there is no concrete algorithm that computes Bessel functions with pure imaginary order despite the fact that they were discovered decade ago.

The author of one of the vastest treatment of Bes-sel functions [5] thought that these functions had no great interest. However these functions were noticed by Lommel [6] and he defined the function (z):

Jv+lß I-

( - ) =

( X 2)V

r(V2)r(v + +1/2)

x[^ (x) + i^ (x)],

(4)

with K and S being real valued functions. This representation by Lommel, however, does not make the computation of these functions any easier. Lommel attempted to obtain a computable integral representation of the Bessel functions to no avail. He however derived the following differential equation as a motivation for the solutions:

x2/ + (2a- 2ßv +1) xy' + +[a (a - 2ßv) + ß2y2 x2ß ] • y = 0,

(5)

this equation has solution (for real nonzero p ) which can be written in the form

y = xpv-a [j (yxp ) + BJ_v (yxp )] . (6)

Comparing equation (5) above with the equation

x2 y" + axy' + ( b + cx2p ) y = 0, (7)

it can be shown [6] that the coefficients in equation (5) and equation (6) have relations shown below

pv - a = -(a -1)/2, (8)

Pv = ^/[( a -1)/2]2 - b, pv =4~c . (9)

It showed that the real valued coefficients ( a, b, c )

gave in general case solution as equation (6) with complex index v . This was also identical to that suggested [3] as possible solution with complex index. At the end, Lommel obtained a very cumbersome solution to equation (7), with real valued coefficients as multiples of a complex valued function. However, no mechanism for the computation of its real valued solution was proposed or mentioned.

Another mathematician M. Bocher approached the modified Bessel function with pure imaginary index while he proposed solutions to Laplace's equation by method of variable separation in cylindrical coordinates [7]. It is very important role of Bessel function with pure imaginary index in mathematical physics.

As opposed to [6], Bocher proposed real valued solutions with the aid of real valued series (however, the convergence of these series was never studied). Each of the coefficients in the series solution suggested in

[7] was partially a polynomial function in v that grew as n!). This approach was found to be intuitive [7]; however, it needed large numbers of iterations as the terms in the series grew. This gave a major problem for its convergence.

Finally the McDonald's function Kiz (x) with pure

imaginary index was used in the integral representation to solve Laplace and wave equations with different boundary values [8]. The author [8] used the representation

' Kiz (x) = jexp (-x cosh t) cos (zt)dt, x > 0. (10)

0

The representation above was found to be valid only for sufficiently large values of the argument x. Unfortunately for small values of x this integral representation is not effective, where pure imaginary indices Bessel functions become applicable and useful in electron optics dynamics [3] and strain field investigations in nanostructures.

The above summary shows that there is not yet an effective algorithm for the computation of Bessel functions with pure imaginary indices.

METHODS AND CALCULATIONS

George Boole [9] developed an operation method to solve differential equation over a century ago. This approach considered calculus as is the treatment of symbols that stands for the operations (e. g., taking the derivative of a function) as if they stood for numerical quantities. Based on the intuition from [9] we can consider the solutions to the equation (1) to be expressed in form

y = A (x) cos (v ln x) + B (x) sin (v ln x), (11)

where A (x) and B (x) are series with the coefficients as shown:

A ( x ) = a0 + a1x + a2 x2 +—anxn + •• B ( x ) = b0 + b1x + b2x2 + •■■bnxn + •••

(12)

Substituting equation (11) into equation (1) and equating the zeros of the coefficients before the linearly independent functions cos (v ln x) and sin (v ln x), we obtain firstly the relation,

a1 = b = 0, (13)

and secondly Vn > 2 the recurrent relationship

nan-2 - 2vbn-2

n(n2 +4v2) '

b = _ lyan-2 + nbn-2

n(n2+4v2)

(15)

From equation (14) and (15), it is clearly seen that odd term coefficients are identical to zeros while those with even values can be obtained via a recurrent relation using the values a0, b0. The series solution to the Bessel equation (1) was first introduced by Boole [9], however, he failed to provide analysis regarding the convergence of the series.

The Boole's recurrent relation thus had to be modified to serve our needs. They can also be readily computed since they seem to exist only for even indices. Let us re-define the function A (x), B (x) such that

,2n m / ^n x ) -A , I x

A (x) = ! a2n 22n I 2 1 = £ A2n I 2

n=0 V 2 I n=0 V 2

Я ( x ) = ! b2n 22n [X] =Z ^2n IX

(16)

(17)

Now for n > 1 the coefficients A2n, B2n can be re-

duced to

A =_ nA2n_2 _vB2n_2

2n n ( n2 +v2 )

B2n ="

vA2n_2 + nB2

n ( n2 +v2 )

(18) (19)

PROOF OF CONVERGENCE AND CONSTRUCTION OF COMPUTABLE SOLUTIONS

In order to compute the strain fields in nanostructures or to implement the Bessel function in other areas of physics, we need to be able to show that the solutions to this differential equation exist and most of all converges in our domain of interest. We now show that the series (14, 15) converges and does so absolutely for any complex value argument and with any value of the orderv . According to [10, 11], a series will converge absolutely if and only if the absolute value of the n -th term (which shall be referred to as the majorant) converges. Let the majorant M2n be defined such that

M 2n = A2n

+ B„

(20)

from (18), (19) and (20) it can follow that

M2n <

+

f

<

v

n2 + v2 ' n (n2 + v2 )

1 H1

— + 4 IM2n_2.

n n I

M2n_2 <

The inequality (20) allows us to be able to test and study the rate and order of which the majorant M2n decays and it is trivial to show that it obeys:

M2n < C

( n!)

(22)

where C is a positive constant defining the zero term of the sequence. Using the d'Alembert's test for convergence [11] the proof of absolute convergence for all values of x and v is completed. We can now without lost of generality form two linearly independent solutions for (1).

Let us choose two pairs of the values for A0 , B0 that are generated from the sequences A2n , B2n with corresponding functions defined by (16) and (17). Under these conditions there exist two solutions for (1). Based on the rate of convergence of the majorant (22), the linearly independent solutions to (1) are shown to be easily computable with algorithm that contains a very limited number of iterations. The solution spawned from the pair

( A,, Bo ) = ( 0,1),

(23)

was notated in [3] as Sfv (x) and Cfv (x) as that spawned from the pair

(A,,Bo) = (1,0) , (24)

then we can easily obtain the analytical expression of these new functions

f ( x ) =

1_

1

1 + v2 V 2

v 12 +.

2

x

-1 +

+

Cfv ( x ) =

1 + v2 V 2

1_

sin ( v ln x ) + cos ( v ln x ),

(25)

1 + v2 V 2

x

- I + ...

cos ( v ln x ) +

+

1 + v2 V 2

x

-1 +...

sin ( v ln x ).

(26)

For small values of the argument x ^ 0, we have two equivalent relationships

1

v

У

n

2

1

2

v

Cfv ( x) ~ cos (v ln x), Sf ( x) ~ sin (v ln x).

(27)

(28)

Hence we see that the symbol used to define this new function are also linked by the equivalent relation and clearly defines physical models earlier proposed [3] and strain field characterization for some nano-structures. Now the general solution of the differential equation (1) can be explicitly written out in the form

y = C1Sf ( x) + C2Cfv ( x).

(30)

x2 y" + xy' + (-x2 + v2 ) y = 0.

Doing this systematically,

cos ( v ln ix) = cosh J ^ J cos (v ln x)-

- sinh |y|sin (v ln x ), (32)

(31)

for (31). Now we can construct two linearly independent solutions for this case. We observe that by choosing two pairs of values C0,D0 obtained from their respective series expressions, we can easily construct two real valued functions being solutions

y1 = C ( x ) cos ( vlnx ) + D ( x ) sin (vlnx ),

(38)

This solution also allows us to easily construct real valued solutions and functions for the modified Bessel function with pure imaginary index, i.e. to solve the equation

y2 = D ( x ) cos (vlnx)- C ( x ) sin (vlnx ). (39)

It is clear that to form two linearly independent solutions to (31) it suffices to take one pair of the spawned sequence

(C0, D0 ) = ( 0,1)

(40)

substituting (34) and (35) into (38) and (39), we obtain two functions denoted [3] as Sdv (x) and Cdv (x):

Sdv ( x ) =

1+

1

1 + v V 2

x

— | +•••

+

1 + vz I 2

2

x

- I +•

sin ( v ln x ) + cos ( v ln x ),

(41)

sin (v ln ix) = cosh ^ J sin (v ln x) +

+ sinh [ ^ ] cos ( v ln x), (33)

where sinh and cosh are the hyperbolic sinus and cosine functions respectively. We can also define functions C ( x ), D (x )

C ( x ) = A ( ix ) = X C2n

V 2 J

n=0 V z J

D ( x ) = B ( ix ) = X D2n (x

(34)

(35)

C = nC2n-2 -vD2n-2 2n (2 2 \ '

n(n +v )

vC2 n-2 + nD2n-2

D2n =

n ( n2 +v2 )

(36)

(37)

Cdv ( x ) =

1+

1 + v2 V 2

2

x

- I +•

cos ( v ln x )-

1 + v2 V 2

2

x

-1 +

sin ( v ln x ), (42)

for small values of the argument x ^ 0, we also have two equivalent relationships

Cdv ( x) ~ cos (v ln x), Sdv ( x) ~ sin ( v ln x),

(43)

(44)

And as in the earlier section, for n > 1 the recurrent equation for C2n, D2n takes the expected form

And it is very trivial to show that (36) and (37) converges absolutely for any x and v. One real valued solution for (1) spawns two real valued solutions

(43) and (44) shows again that these functions are linearly independent and also solutions to (31).

ALGORITHM FOR THE COMPUTATION OF FUNCTION

Cdv ( x ), Sdv (x ), Cfv ( x ), Sfv ( x )

An important property that has to be analyzed for easy computation of these functions is the Wronskian. In mathematics, the Wronskian is a function named after the Polish mathematician Jozef Hoene-Wronski., where it can be used to determine whether a set of solutions are linearly independent [12, 13]. It can be easily shown that the Wronskian W of two linearly independent solutions of the equation

2

v

1

v

/ + Xy' + f ( x ) y = 0,

(45)

CL (x)[Sf (x)]' - Sfv (x)[Cf (x)]' = x, (46) Cdv (x)[Sdv (x)]' - Sdv (x)[Cdv (x)]' = V. (47)

Choosing more terms of the Taylor's series, we can further give the possibility of increasing the accuracy

has the form W = —, with the constant a depending

x

on the concrete solution. The equivalent relationship (27), (28) and (43), (44) allows us to be able to obtain a

m2„ , v2 F

—— < 1 + — + — „2 3 '

m

2n-2

n n

(50)

Taking the natural logarithm of both sides of (50), we obtain:

ln ( m2n )-ln ( m2n-2 )< ln(1 + V2 + F) < V2 + Hf.

n n n n

Using the well known sums below,

It clearly shows that the Wronskian (46) and (47) vanishes for v = 0 . Leading to the following relationship

v ^ 0 ^ Cfv (x J0 ( x),

Sfv (x0.

w 1

I Jr = - -1,

n2 6

1

I— = C(3) -1 = 0.2021,

n"2 n

and also the value for

ACCURACY OF COMPUTATIONAL ALGORITHM

To obtain a measure of the computational error, we required a detailed study of the majorant (remainder) M2n. To understand this, let us consider the following series

M 2n =

( n!)

2 m2n.

(48)

m

<

m

( IvP 1+ U

n

V

/

1 —I .

n

Using Lagrange's formula for remainders, the Taylor's series [10, 11] gives for n > 2

< 1 - +lvl 1 lvl

m2n-2 n2 6n3

(fcj) f1

01v\ 2 --I X

n )

d

3|v| + (|v-2)1 1 + |v|- II 1 --

d

(49)

where 0 < 6 < 1.

We can thus reduce the inequality (49) by removing from the right hand side 6

< 1 - v!+F

~ 2 3 '

m2n-2 n n

mo = I

1 + |v|

finally, we obtain

injmn ^^.-2

m„„

< 0.6449v + 0.2021F

(51)

with

Now we obtain an estimate for remainder m2n for

n>2

1 + v

-exp(0.6449v2 + 0.2021F) = m .

2

1 + v2

Now, the accuracy for M2n takes the form

M2n <

(n!)2

(53)

We can conclude by saying that the sum of the remaining terms for the series functions A (x), B (x), C (x), D (x) cannot exceed the value

= m I

lvl / \2n

nn ( x

n=N(n!)212,

Since the remainder will depend on the value of v, we thus study the case |v| < 2 , then

x

< mT" I

1

[(n -1)!]2

= mT" I

1

4 tN(n!)212)

x 1

(54)

m2

v

n

v

v

v -3

X

2n-2

x

2

n=N+1

l

0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1

A Cfv (x) - Vsfv(x) 1 1 V\ /Xx

f . Vy^ v> 1 1

10

15

20

25

30

-1

1 1 Cd„ [xy

-—

\ ^^^^^ ~~ 1 1 1

0.5

1.5

2.5

Fig. 1. Functions Sfy2(x) and Cf1/2 (x) for 0 < x < 30 Fig. 2. Functions Sd1/2(x) and Cdl/2(x) for 0 < x < 3

The sum at the RHS of this series is simply the remainder of the series for the modified Bessel function of zero 10 (x). This remainder can also be evaluated further as shows below

Z

1

N ( n!) 1

'x^2n 2

/ \ 2 N

x

(N !)2

(N!)212 J Z [(„ + n)!]2 V 2J

(N !)2

'x^2 N

V^y

1+-

x2

(1 + N )2

+ ••• +

v^y

'x^2n

v 2 y

[(1 + N)(2 + N )•••( n + N)]

> <

1

(N!) 1 +

1

( x Y N

v 2 y 1

1 4 • 2

^2 2

n (n +1)!

^2 N

( N !)2

vy

2 I (x ) = —^

x (N!)2

'x^2n

/ N 2 N-1

x

I1 (x).

vy

Now for any x and |v| we have

1 + |v| / 2 \

sN <-4-exp(0.6449v + 0.202LF):

N 1+v2 V }

( N !)2

/ \2 N+1

x

v2y

I1 (x).

(55)

Eliminating v we obtain estimates for x and |v| 1

^n < 15

2 N +1

x

(N!)2

2

I1 (x).

Partial for x < 2 , the equation above reduces to

24

^ N <

(N !)2

(56)

This shows that for the functions A (x), B (x), C (x), D (x) and for small values of the arguments and the order x < 2 and |v| < 2 we have a computational error s8 < 1.5 x 10-16 after the computing only eight members (elements) of the series.

Graphs of Cdv (x), Sdv (x), Cfv (x), Sfv (x) for two different v are shown in Fig. 1-4.

CONCLUSION

We have systematically introduced and define a new algorithm to compute the real valued solutions of some special functions namely the Gamma function with complex order, the Bessel (and modified Bessel) function of pure imaginary orders. This can provide an alternate approach to solving the ground state of the screened coulombs potential problem. It is also promising in solving the inverse scattering problem which is mostly approached via the inverse green

1

x

2 n

1

1

+

+

<

X

1

+

+

2

Fig. 3. Functions Sf3/2 (x) and Cf3/2 (x) for 0 < x < 30

function. In the above applied cases, the value of the arguments of the Bessel functions lies in the vicinity of unity. The algorithm for their computations is shown to be rapidly converging and gives a very favorable accuracy for as few as eight iterations. Further studies of the asymptotic behavior of these functions are still warranted as they play a great role in boundary valued problems, and integral and differential problems of wave scattering.

REFERENCES

1. Chandrasekhar S., Elbert O. // Proc. Cambr. Phil. SOC. 1954. N 5_0. P. 266.

2. Cochranr J.P. // Radio Sciences. 1966. N 1. P. 679.

3. Matyshev A.A. Isotrajectory corpuscular optics. SPb: Science, 2000. 376 p. (in Russian).

4. Abramowitz M., Segun I.A. Handbook of mathematical functions. N.Y.: Dover, 1968.

5. Whittaker E.T., Watson G.N. A course of modern analysis. Cambridge: Cambridge University Press, 1927. 608 p.

6. Lommel E. Zur theorie der bessel'schen funktionen // Math. Ann. 1871. Bd 3. S. 475-487.

7. Bocher M. On some applications of Bessel's function with pure imaginary index // Annals of Math. 1892. V. 6, N 6. P.137-160.

8. Greenberg G.A. Selected problems of mathematical theory of electromagnetic phenomena. M.-L.: Acade-

1.5 -

Cd„ (x)

0.5 -r\

Fig. 4. Functions Sd3,2 (x) and Cd3/2 (x) for 0 < x<3

my of Scienses Publishing, 1948. 728 p. (in Russian).

9. Boole G. On a general method in analysis // Phil. Trans. of the Roy. Soc. of London. 1844. V. 134, Part II. P. 225-282.

10. Dvoretzky A., Rogers A., Ambrose C. Absolute and unconditional convergence in normed linear spaces // Proc. Nat. Acad. Sci. USA. 1950. V. 36. P. 192-197.

11. Rudin W. Principles of mathematical analysis. N.Y.: McGraw-Hill, 1964. 270 p.

12. Watson G.N. A treatise on the theory of Bessel functions. Cambridge University Press, 1945.

13. Gray A., Mathews G.B. A treatise on Bessel functions and their applications to physics. N.Y.: Macmillan and Company, 1895.

St. Petersburg State Polytechnic University, RF (Matyshev A.A.)

New Mexico State University, Department of Physics, USA (Fohtung E.B.)

Контакты: Матышев Александр Александрович., [email protected]

Материал поступил в редакцию 21.01.2014

О ВЫЧИСЛЕНИИ И ПРИМЕНЕНИЯХ БЕССЕЛЕВЫХ ФУНКЦИЙ ЧИСТО МНИМОГО ИНДЕКСА (ПОРЯДКА) В ФИЗИКЕ И ОСОБЕННО В КОРПУСКУЛЯРНОЙ ОПТИКЕ

А. А. Матышев1, E. B. Fohtung2

1 Санкт-Петербургский государственный политехнический университет 2New Mexico State University, Department of Physics, USA

Предложен простой и эффективный алгоритм расчета бесселевых функций чисто мнимого индекса (порядка). Указаны физические проблемы, особенно в области корпускулярной оптики, которые аналитически разрешаются с помощью подобных функций.

Кл. сл.: бесселевы функции чисто мнимого индекса (порядка), корпускулярная оптика