Appendix C:

Bessel Functions

In the solution of problems involving cylindrical geometries, we will

repeatedly encounter Bessel functions. In this appendix, we highlight
their most useful properties.

C.1 BESSEL FUNCTIONS AND THEIR PROPERTIES

Bessel’s equation1 is

x2 y

+ xy
+ (µ2 x2 − n2 )y = 0, (C.1.1)

or    
d dy n2
x + µ x−
2
y = 0. (C.1.2)
dx dx x
Because we cannot write down the solution to (C.1.1) or (C.1.2) in a
simple closed form and since Bessel’s equation is singular at x = 0, the
solution is sought as a power series given by the method of Frobenius.

1 The classic reference on Bessel functions is Watson, G. N., 1966: A Treatise on

the Theory of Bessel Functions. Cambridge University Press, 804 pp.

© 2001 by Chapman & Hall/CRC

 1.0

J0 (x)
J1(x)
J2(x)
0.5 J3(x)

0.0

-0.5
0.0 2.0 4.0 6.0 8.0
x
Figure C.1.1: The first four Bessel functions of the first kind over 0 ≤ x ≤ 8.

Since n is usually an integer in this book, we present the results for this
case.
Frobenius’ method yields two linearly independent series that are
written down elsewhere but can be treated as tabulated functions.2 The
first solution to (C.1.1)–(C.1.2) is called a Bessel function of the first
kind of order n, is denoted by Jn (x), and behaves as xn , as x → 0. The
second general solution is Neumann’s Bessel function of the second kind
of order n, Yn (x). The Neumann functions Y0 (x) and Yn (x) behave as
ln(x) and −x−n , respectively, as x → 0 and n > 0. Consequently, the
general solution to (C.1.1) is

y(x) = AJn (µx) + BYn (µx). (C.1.3)

Figure C.1.1 illustrates the functions J0 (x), J1 (x), J2 (x), and J3 (x)
while Figure C.1.2 gives Y0 (x), Y1 (x), Y2 (x), and Y3 (x).
The Bessel functions Jn and Yn enjoy certain recurrence relations
where functions of order n are related to functions of order n + 1 and
n − 1. If we denote Jn or Yn by Cn , the most useful relationships are

Cn−1 (z) + Cn+1 (z) = 2nCn (z)/z, n = 1, 2, 3, . . . , (C.1.4)

Cn−1 (z) − Cn+1 (z) = 2Cn
(z), n = 1, 2, 3, . . . , (C.1.5)

2 See Press, W. H., B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, 1986:

Numerical Recipes: The Art of Scientific Computing. Cambridge University Press,
§6.4.

© 2001 by Chapman & Hall/CRC

 1.0

Y0 (x) Y1 (x)
0.5 Y2 (x)
Y3(x)

0.0

-0.5

-1.0
0.0 2.0 4.0 6.0 8.0
x
Figure C.1.2: The first four Bessel functions of the second kind over 0 ≤ x ≤ 8.

d n
[z Cn (z)] = z n Cn−1 (z), n = 1, 2, 3, . . . , (C.1.6)
dz
and
d  −n 
z Cn (z) = −z −n Cn+1 (z), n = 1, 2, 3, . . . , (C.1.7)
dz
with the special case of J0
(z) = −J1 (z), and Y0
(z) = −Y1 (z).
An equation that is very similar to (C.1.1) is
d2 y dy
x22
+x − (n2 + x2 )y = 0. (C.1.8)
dx dx
√
If we substitute ix = t (where i = −1 ) into (C.1.8), it becomes Bessel’s
equation:

d2 y dy
2
t2
+t + (t2 − n2 )y = 0. (C.1.9)
dt dt
Consequently, we can immediately write the solution to (C.1.8) as

y(x) = c1 Jn (ix) + c2 Yn (ix), (C.1.10)

if n is an integer. Traditionally the solution to (C.1.10) has been written

y(x) = c1 In (x) + c2 Kn (x) (C.1.11)

© 2001 by Chapman & Hall/CRC

 5.0

4.0

3.0

2.0 I0
I1
I2
1.0
I3

0.0
0.0 1.0 2.0 3.0
x
Figure C.1.3: The first four modified Bessel functions of the first kind over 0 ≤
x ≤ 3.

rather than in terms of Jn (ix) and Yn (ix) The function In (x) is the
modified Bessel function of the first kind, of order n, while Kn (x) is the
modified Bessel function of the second kind, of order n. Figure C.1.3
illustrates I0 (x), I1 (x), I2 (x), and I3 (x), while in Figure C.1.4 K0 (x),
K1 (x), K2 (x), and K3 (x) have been graphed. Again, both In and Kn
can be treated as tabulated functions.3 Note that Kn (x) has no real
zeros while In (x) equals zero only at x = 0 for n > 0.
As our derivation suggests, modified Bessel functions are related to
ordinary Bessel functions via complex variables. In particular, Jn (iz) =
in In (z) and In (iz) = in Jn (z) for z complex.
As in the case of Jn and Yn , there are recurrence relations where
functions of order n are related to functions of order n + 1 and n − 1. If
we denote by Cn either In or enπi Kn , we have

Cn−1 (z) − Cn+1 (z) = 2nCn (z)/z, n = 1, 2, 3, . . . , (C.1.12)

Cn−1 (z) + Cn+1 (z) = 2Cn
(z), n = 1, 2, 3, . . . , (C.1.13)

d n
[z Cn (z)] = z n Cn−1 (z), n = 1, 2, 3, . . . , (C.1.14)
dz
and
d  −n 
z Cn (z) = z −n Cn+1 (z), n = 1, 2, 3, . . . , (C.1.15)
dz

3 Ibid., §6.5.

© 2001 by Chapman & Hall/CRC

 4.0

3.0

2.0

K1 K2 K3

1.0
K0

0.0
0.0 1.0 2.0 3.0
x
Figure C.1.4: The first four modified Bessel functions of the second kind over
0 ≤ x ≤ 3.

with the special case of I0
(z) = I1 (z), and K0 (z) = −K1 (z).

C.2 HANKEL FUNCTIONS

Hankel functions are an alternative to Jn and Yn . Here we define

them and show why they are useful.
One of the fundamental equations studied in any differential equa-
tions course is
y

+ k 2 y = 0. (C.2.1)
Its solution is generally written

y(x) = A cos(kx) + B sin(kx), (C.2.2)

where A and B are arbitrary constants. Another, less often quoted,

solution is
y(x) = Ceikx + De−ikx , (C.2.3)
where C and D are generally complex constants. The advantage of
(C.2.3) over (C.2.2) is the convenience of dealing with exponentials
rather than sines and cosines.
Let us return to (C.1.1):
 
x2 y

+ xy
+ x2 − n2 y = 0, (C.2.4)

© 2001 by Chapman & Hall/CRC

Table C.2.1: Some Useful Recurrence Relations for Hankel Functions.

d  n (p) (p) d  (p) (p)

x Hn (x) = xn Hn−1 (x), n = 1, 2, . . . ; H0 (x) = −H1 (x)
dx dx

d  −n (p)
x Hn (x) = −x−n Hn+1 (x),
(p)
n = 0, 1, 2, 3, . . .
dx

(p) (p) 2n (p)

Hn−1 (x) + Hn+1 (x) = H (x), n = 1, 2, 3, . . .
x n

(p)
(p) (p) dHn (x)
Hn−1 (x) − Hn+1 (x) = 2 , n = 1, 2, 3, . . .
dx

where n is a positive (including zero) integer. Its general solution is

most often written
y(x) = A Jn (x) + B Yn (x). (C.2.5)
Again A and B are arbitrary constants. However, another, equivalent
solution is
y(x) = CHn(1) (x) + DHn(2) (x), (C.2.6)
where
Hn(1) (x) = Jn (x) + iYn (x), (C.2.7)
and
Hn(2) (x) = Jn (x) − iYn (x). (C.2.8)
(1) (2)
These functions Hn (x), Hn (x) are referred to as Bessel functions of
the third kind or Hankel functions, after the German mathematician
Hermann Hankel (1839–1873). Some of their properties are listed in
Table C.2.1.
The similarity between the solution (C.2.3) and (C.2.6) is clarified
by the asymptotic expansions for the Hankel functions:

2 i(z−nπ/2−π/4)
Hn (z) ∼
(1)
e , (C.2.9)
πz
and

2 −i(z−nπ/2−π/4)
Hn(2) (z) ∼ e . (C.2.10)
πz

© 2001 by Chapman & Hall/CRC

C.3 FOURIER-BESSEL SERIES

A Fourier-Bessel series is a Fourier series where Jn ( ) is used in place

of the sine and cosine. The exact form of the expansion depends upon
the boundary condition at x = L. There are three possible cases. One
of them is that y(L) = 0 and results in the condition that Jn (µk L) = 0.
Another condition is y
(L) = 0 and gives Jn
(µk L) = 0. Finally, if
hy(L) + y
(L) = 0, then hJn (µk L) + µk Jn
(µk L) = 0. In all of these
cases, the eigenfunction expansion is the same, namely

∞

f (x) = Ak Jn (µk x), (C.3.1)
k=1

where µk is the kth positive solution of either Jn (µk L) = 0, Jn
(µk L) =

0, or hJn (µk L) + µk Jn
(µk L) = 0.
We now need a mechanism for computing Ak . We begin by multi-
plying (C.3.1) by xJn (µm x) dx and integrate from 0 to L. This yields
∞
  L  L
Ak xJn (µk x)J(µm x) dx = xf (x)Jn (µm x) dx. (C.3.2)
k=1 0 0

From the general orthogonality condition

 L
xJn (µk x)Jn (µm x) dx = 0 (C.3.3)
0

if k = m, (C.3.2) simplifies to
 L  L
Am xJn2 (µm x) dx = xf (x)Jn (µm x) dx, (C.3.4)
0 0

or

 L
1
Ak = xf (x)Jn (µk x) dx, (C.3.5)
Ck 0

where  L
Ck = xJn2 (µk x) dx (C.3.6)
0

© 2001 by Chapman & Hall/CRC

and k has replaced m in (C.3.4).
The factor Ck depends upon the boundary conditions at x = L. It
can be shown that
Ck = 12 L2 Jn+1
2
(µk L), (C.3.7)
if Jn (µk L) = 0. If Jn
(µk L) = 0, then

µ2k L2 − n2 2
Ck = Jn (µk L). (C.3.8)
2µ2k

Finally,
µ2k L2 − n2 + h2 L2 2
Ck = Jn (µk L), (C.3.9)
2µ2k
if µk Jn
(µk L) = −hJn (µk L).
All of the preceding results must be slightly modified when n = 0
and the boundary condition is J0
(µk L) = 0 or µk J1 (µk L) = 0. This
modification arises from the additional eigenvalue µ0 = 0 being present
and we must add the extra term A0 to the expansion. For this case,
(C.3.1) becomes
∞

f (x) = A0 + Ak J0 (µk x), (C.3.10)
k=1

where the equation for finding A0 is

 L
2
A0 = f (x) x dx. (C.3.11)
L2 0

Equations (C.3.5) and (C.3.8) with n = 0 give the remaining coefficients.

• Example C.3.1

Let us expand f (x) = x, 0 < x < 1, in the series

∞

f (x) = Ak J1 (µk x), (C.3.12)
k=1

where µk denotes the kth zero of J1 (µ). From (C.3.5) and (C.3.7),
 1
2
Ak = 2 x2 J1 (µk x) dx. (C.3.13)
J2 (µk ) 0

However, from (C.1.6),

d  2 
x J2 (x) = x2 J1 (x), (C.3.14)
dx

© 2001 by Chapman & Hall/CRC

 1.25

1.00

0.75

0.50

0.25

0.00 one term two terms

-0.25
1.25

1.00

0.75

0.50

0.25

0.00 three terms four terms

-0.25
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0
x x

Figure C.3.1: The Fourier-Bessel series representation (C.3.16) for f (x) = x, 0 ≤

x < 1, when we truncate the series so that it includes only the first, first two, first
three, and first four terms.

if n = 2. Therefore, (C.3.13) becomes

µ
2x2 J2 (x) k 2
Ak = 3 2 = , (C.3.15)
µk J2 (µk ) 0 µk J2 (µk )
and the resulting expansion is
∞
 J1 (µk x)
x=2 , 0 ≤ x < 1. (C.3.16)
µk J2 (µk )
k=1

Figure C.3.1 shows the Fourier-Bessel expansion of f (x) = x in trun-

cated form when we only include one, two, three, and four terms.

• Example C.3.2

Let us expand the function f (x) = x2 , 0 < x < 1, in the series

∞

f (x) = Ak J0 (µk x), (C.3.17)
k=1

© 2001 by Chapman & Hall/CRC

 1.00

0.75

0.50

0.25

0.00

-0.25 one term two terms

-0.50
1.00

0.75

0.50

0.25

0.00

-0.25 three terms four terms

-0.50
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0
x x
Figure C.3.2: The Fourier-Bessel series representation (C.3.27) for f (x) = x2 ,
0 < x < 1, when we truncate the series so that it includes only the first, first two,
first three, and first four terms.

where µk denotes the kth positive zero of J0 (µ). From (C.3.5) and
(C.3.7),
 1
2
Ak = 2 x3 J0 (µk x) dx. (C.3.18)
J1 (µk ) 0
If we let t = µk x, the integration (C.3.18) becomes
 µk
2
Ak = 4 2 t3 J0 (t) dt. (C.3.19)
µk J1 (µk ) 0
We now let u = t2 and dv = tJ0 (t) dt so that integration by parts gives
  µk 
2 µk
Ak = t 3
J1 (t) −2 t 2
J1 (t) dt (C.3.20)
µ4k J12 (µk ) 0
0
  µk 
2
= 4 2 µ3k J1 (µk ) − 2 t2 J1 (t) dt , (C.3.21)
µk J1 (µk ) 0
because v = tJ1 (t) from (C.1.6). If we integrate by parts once more, we
find that
 
2
Ak = 4 2 µ J1 (µk ) − 2µk J2 (µk )
3 2
(C.3.22)
µk J1 (µk ) k

© 2001 by Chapman & Hall/CRC

  
2 J1 (µk ) 2J2 (µk )
= 2 − . (C.3.23)
J1 (µk ) µk µ2k

However, from (C.1.4) with n = 1,

J1 (µk ) = 12 µk [J2 (µk ) + J0 (µk )] , (C.3.24)

or
2J1 (µk )
J2 (µk ) = , (C.3.25)
µk
because J0 (µk ) = 0. Therefore,

2(µ2k − 4)J1 (µk )

Ak = , (C.3.26)
µ3k J12 (µk )

and
∞
 (µ2k − 4)J0 (µk x)
x2 = 2 , 0 < x < 1. (C.3.27)
µ3k J1 (µk )
k=1

Figure C.3.2 shows the representation of x2 by the Fourier-Bessel series

(C.3.27) when we truncate it so that it includes only one, two, three, or
four terms. As we add each additional term in the orthogonal expansion,
the expansion fits f (x) better in the “least squares” sense.

© 2001 by Chapman & Hall/CRC