Bessel functions

Bessel functions of all kinds and their properties.

In this chapter 7

Tick a unit to add it to your PDF basket.

Bessel functions

Bessel functions

Definition

Bessel functions of the first kind

Adapted

Add to PDF

Bessel's equation has the Frobenius solution $J_\nu(x) = \sum_m (-1)^m (x/2)^{\nu+2m}/[m!\,\Gamma(\nu+m+1)]$ , convergent for all $x$ and regular at the origin when $\nu \ge 0$ ; for integer order the whole family is packed into the generating function $e^{(x/2)(t-1/t)} = \sum_n J_n(x)t^n$ , and four recurrence relations connect neighbouring orders and derivatives.

Prerequisites: The Frobenius method, The Gamma function, Classifying singular points, including infinity

Definition(Bessel's equation and the function

J_\nu

)AdaptedAdapted from Riley, Hobson & Bence (2006), §18.5-18.5.1, Eqs. (18.73)-(18.79), pp. 602-605; Arfken, Weber & Harris (2013), §14.1, Eqs. (14.1) and (14.6), pp. 643-645; Strauss (2008), §10.5, Eqs. (1) and (4), pp. 282-283

Bessel's equation of order $\nu$ is

x^{2}y'' + x\,y' + (x^{2}-\nu^{2})\,y = 0 ,

with $\nu$ a fixed parameter. Dividing by $x^{2}$ puts it in the normalised form $y'' + x^{-1}y' + (1-\nu^{2}x^{-2})y = 0$ , so $x = 0$ is a regular singular point and $x = \infty$ is an irregular one (Classifying singular points, including infinity). Replacing $\nu$ by $-\nu$ changes nothing, so we may take $\mathrm{Re}\,\nu \ge 0$ without loss.

Substituting a Frobenius series $y = x^{\sigma}\sum_{n\ge0}a_nx^{n}$ (The Frobenius method) gives the indicial equation $\sigma^{2}-\nu^{2} = 0$ , so the exponents are $\sigma = \pm\nu$ ; the coefficient of $x^{\sigma+1}$ gives $\bigl[(\sigma+1)^{2}-\nu^{2}\bigr]a_1 = 0$ , and for $n \ge 2$

\bigl[(\sigma+n)^{2}-\nu^{2}\bigr]a_n + a_{n-2} = 0 .

Take the larger exponent $\sigma = \nu$ . Then the $n = 1$ equation reads $(1+2\nu)a_1 = 0$ and we set $a_1 = 0$ ; the recurrence becomes $n(n+2\nu)\,a_n = -a_{n-2}$ , which kills every odd coefficient and determines the even ones from $a_0$ . The conventional normalisation is

a_0 = \frac{1}{2^{\nu}\,\Gamma(1+\nu)} ,

with $\Gamma$ the Gamma function; iterating the recurrence then gives $a_{2m} = (-1)^{m}\bigl[2^{\nu+2m}\,m!\,\Gamma(\nu+m+1)\bigr]^{-1}$ . The resulting function is the Bessel function of the first kind of order $\nu$ ,

J_\nu(x) \;=\; \sum_{m=0}^{\infty}\frac{(-1)^{m}}{m!\;\Gamma(\nu+m+1)}\left(\frac{x}{2}\right)^{\nu+2m} .

Four facts about this series:

It converges for every $x$ . The ratio of consecutive terms is $-x^{2}/\bigl[4m(\nu+m)\bigr]$ , which tends to $0$ ; the singularity at $x=\infty$ is what the Frobenius series near $x=0$ cannot see.
It is defined for every real order, once one convention is fixed. The convention is that $1/\Gamma$ is read as the entire function it is, so that $1/\Gamma(z) = 0$ at each pole $z = 0,-1,-2,\dots$ of $\Gamma$ (The Gamma function). This handbook adopts that reading throughout, and the identity below and the generating function depend on it. With it, $\nu = -n$ for $n \in \mathbb{N}$ causes no trouble: the first $n$ coefficients are $0$ rather than undefined, so the sum begins at $m = n$ , and shifting $m \to m+n$ gives $J_{-n}(x) = (-1)^{n}J_n(x) \qquad (n \in \mathbb{Z}) ,$ so at integer order the two exponents $\pm n$ produce the same solution up to sign, and a second solution must be built differently. The other reading is also in print: Strauss (Strauss (2008), §10.5, p. 283) declares the series undefined at negative integer order, on the ground that $\Gamma(\nu+m+1)$ is infinite there, and so states his existence claim for every real order except the negative integers. Nothing mathematical separates the two; what follows is written in the $1/\Gamma$ convention, and where a statement is taken from Strauss it is quoted at his orders, not at ours.
Behaviour at the origin. The leading term is $(x/2)^{\nu}/\Gamma(1+\nu)$ , so $J_0(0) = 1$ , $J_\nu(0) = 0$ for $\nu > 0$ , and $J_\nu$ is unbounded at $0$ for $\nu < 0$ and not an integer, the leading power $x^{\nu}$ then being genuinely negative. The exception matters: at $\nu = -n$ the leading term displayed above is $0$ , the bullet above gives $J_{-n} = (-1)^{n}J_n$ , and $J_{-n}$ is as bounded at the origin as $J_n$ is — for instance $J_{-1}(0.001) = -5.0\times10^{-4}$ and $J_{-2}(0.5) = J_2(0.5) = 0.0306040$ . Strauss (Strauss (2008), §10.2, p. 267) records that for integer order $J_n$ is, up to a constant factor, the only solution of Bessel's equation finite at $x = 0$ — which is why it is the one that survives in a problem whose domain contains the axis.
The name. "First kind" is Arfken, Weber and Harris's (Arfken, Weber & Harris (2013), §14.1, p. 643) label for exactly these solutions, the ones regular at the origin for non-negative integer order; every solution independent of $J_\nu$ is irregular there.

For $\nu \notin \mathbb{Z}$ the exponents $\nu$ and $-\nu$ differ by a non-integer, so $J_\nu$ and $J_{-\nu}$ are two independent Frobenius solutions and the general solution of Bessel's equation is $c_1J_\nu + c_2J_{-\nu}$ . The independence is confirmed by a Wronskian computation, carried out in Neumann and Hankel functions.

Definition(Generating function for integer order)AdaptedAdapted from Arfken, Weber & Harris (2013), §14.1, Eqs. (14.2)-(14.5), pp. 644-645; Riley, Hobson & Bence (2006), §18.5.3, Eq. (18.98), pp. 612-613

For every $x \in \mathbb{C}$ and every $t \ne 0$ ,

\exp\!\left[\frac{x}{2}\left(t - \frac{1}{t}\right)\right] \;=\; \sum_{n=-\infty}^{\infty} J_n(x)\,t^{n} ,

the series converging absolutely. That is, the Laurent coefficients in $t$ of the left-hand side are precisely the Bessel functions of integer order, at exactly the scale fixed by the choice of $a_0$ above.

Setting $t = e^{i\theta}$ turns this into the Jacobi-Anger expansion

e^{\,ix\sin\theta} \;=\; \sum_{n=-\infty}^{\infty} J_n(x)\,e^{\,in\theta} ,

whose real and imaginary parts, combined with the orthogonality of $\{e^{in\theta}\}$ , give the integral representation

J_n(x) \;=\; \frac{1}{\pi}\int_0^{\pi}\cos\bigl(x\sin\theta - n\theta\bigr)\,d\theta \qquad (n \in \mathbb{Z}) .

Nothing of the kind is available at non-integer order: there is no Laurent series, and the scale of $J_\nu$ is then fixed by convention alone, chosen to agree with the integer case ((Arfken, Weber & Harris (2013), §14.1, p. 653)). Further representations, valid for general $\nu$ , are the subject of Integral representations of Bessel functions.

Idea of the proofOriginal

The generating function works because $e^{xt/2}$ and $e^{-x/(2t)}$ are each a single geometric-looking series, one in positive powers of $t$ and one in negative powers. Multiplying them and collecting the coefficient of $t^{n}$ pairs an $r$ from the first factor with an $s = r-n$ from the second; the double sum collapses to a single sum over $s$ , and what is left is literally the Frobenius series with $m = s$ . The two recurrence relations then come for free: differentiating the exponential with respect to $t$ and with respect to $x$ reproduces the same exponential times an elementary factor, so equating Laurent coefficients turns calculus into algebra among neighbouring orders.

Proof(The generating function)AdaptedAdapted from Arfken, Weber & Harris (2013), §14.1, Eqs. (14.2)-(14.5), pp. 644-645; Riley, Hobson & Bence (2006), §18.5.3, p. 613

Fix $x$ and $t \ne 0$ and expand each factor of $e^{xt/2}\,e^{-x/(2t)}$ by the exponential series:

e^{xt/2} = \sum_{r=0}^{\infty}\frac{1}{r!}\left(\frac{x}{2}\right)^{r}t^{r} , \qquad e^{-x/(2t)} = \sum_{s=0}^{\infty}\frac{(-1)^{s}}{s!}\left(\frac{x}{2}\right)^{s}t^{-s} .

Both converge absolutely, so their product may be multiplied out and the terms rearranged at will. Doing so,

\exp\!\left[\frac{x}{2}\left(t-\frac1t\right)\right] = \sum_{r=0}^{\infty}\sum_{s=0}^{\infty}\frac{(-1)^{s}}{r!\,s!}\left(\frac{x}{2}\right)^{r+s}t^{\,r-s} .

Now collect the terms with a fixed power $t^{n}$ , that is, those with $r = n+s$ . Since $r \ge 0$ , the admissible $s$ are $s \ge \max(0,-n)$ , and $r+s = n+2s$ . Hence the coefficient of $t^{n}$ is

\sum_{s \ge \max(0,-n)}\frac{(-1)^{s}}{(n+s)!\,s!}\left(\frac{x}{2}\right)^{n+2s} .

Two cases finish the identification.

$n \ge 0$ . The sum runs over all $s \ge 0$ , and since $(n+s)! = \Gamma(n+s+1)$ it is exactly the series defining $J_n(x)$ .
$n < 0$ . Write $n = -k$ with $k > 0$ . The sum runs over $s \ge k$ ; substituting $s = k+m$ turns it into $\sum_{m\ge0}(-1)^{k+m}\bigl[m!\,(k+m)!\bigr]^{-1}(x/2)^{k+2m} = (-1)^{k}J_k(x)$ , which by the identity $J_{-k} = (-1)^{k}J_k$ established in the definition is $J_{-k}(x) = J_n(x)$ .

So in both cases the coefficient of $t^{n}$ is $J_n(x)$ , which is the claim. Absolute convergence of the double series gives absolute convergence of the Laurent series.

For the Jacobi-Anger form put $t = e^{i\theta}$ , so that $t - 1/t = 2i\sin\theta$ and the exponent becomes $ix\sin\theta$ . Multiplying the resulting identity by $e^{-in\theta}$ and integrating over $\theta \in [-\pi,\pi]$ , the orthogonality $\int_{-\pi}^{\pi}e^{i(m-n)\theta}d\theta = 2\pi\delta_{mn}$ leaves

J_n(x) = \frac{1}{2\pi}\int_{-\pi}^{\pi}e^{\,i(x\sin\theta - n\theta)}\,d\theta .

The imaginary part of the integrand, $\sin(x\sin\theta-n\theta)$ , is odd in $\theta$ for integer $n$ and integrates to zero; the real part is even, which halves the range and gives the stated formula. $\square$

Definition(Recurrence relations)AdaptedAdapted from Riley, Hobson & Bence (2006), §18.5.3, Eqs. (18.92)-(18.97), pp. 611-612; Arfken, Weber & Harris (2013), §14.1, Eqs. (14.7)-(14.12), pp. 645-646; Hassani (2013), §15.5, Eqs. (15.47)-(15.48), p. 485

For every $\nu$ for which the functions are defined,

\frac{d}{dx}\bigl[x^{\nu}J_\nu(x)\bigr] = x^{\nu}J_{\nu-1}(x) , \qquad \frac{d}{dx}\bigl[x^{-\nu}J_\nu(x)\bigr] = -x^{-\nu}J_{\nu+1}(x) ,

equivalently

x\,J_\nu'(x) + \nu\,J_\nu(x) = x\,J_{\nu-1}(x) , \qquad x\,J_\nu'(x) - \nu\,J_\nu(x) = -x\,J_{\nu+1}(x) .

Subtracting and adding these two gives the three-term recurrence and the derivative recurrence

J_{\nu-1}(x) + J_{\nu+1}(x) = \frac{2\nu}{x}J_\nu(x) , \qquad J_{\nu-1}(x) - J_{\nu+1}(x) = 2J_\nu'(x) ,

of which the case $\nu = 0$ of the second, with $J_{-1} = -J_1$ , is

J_0'(x) = -J_1(x) .

Two remarks on how far this reaches.

The relations do not single out $J_\nu$ . They hold for $Y_\nu$ and for the Hankel functions as well, and more generally for any family built from solutions of Bessel's equation at the standard scale ((Hassani (2013), §15.5, p. 485)); the pair of first-order relations says that $x^{\pm\nu}\frac{d}{dx}x^{\mp\nu}$ raises and lowers the order.
Conversely they force Bessel's equation. Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §14.1, Eqs. (14.13)-(14.14), pp. 646-647) show that any family $Z_\nu$ satisfying the three-term and derivative recurrences satisfies $x^{2}Z_\nu'' + xZ_\nu' + (x^{2}-\nu^{2})Z_\nu = 0$ . So the recurrences are not a by-product of the equation; they are equivalent to it.

Read backwards, the first pair are antiderivatives:

\int x^{\nu}J_{\nu-1}(x)\,dx = x^{\nu}J_\nu(x) , \qquad \int x^{-\nu}J_{\nu+1}(x)\,dx = -x^{-\nu}J_\nu(x) .

Idea of the proofOriginal

Only two of the four recurrences need proving; the rest is algebra. The two that do are the ones with a prefactor $x^{\pm\nu}$ , and the reason they are the right ones is that $x^{\pm\nu}$ cancels the fractional power in the series and leaves a power series in $x^{2}$ whose exponents are $2\nu+2m$ and $2m$ . Differentiating then multiplies the $m$ -th term by that exponent — by $2(\nu+m)$ in the first case, by $2m$ in the second — and in each case the new factor is exactly what the Gamma function needs to shift its argument by one, because $\Gamma(\nu+m+1) = (\nu+m)\Gamma(\nu+m)$ and $m! = m\,(m-1)!$ . So differentiation shifts the order by one. The two cases differ in one detail with a large consequence: the factor $2m$ annihilates the $m = 0$ term, so the second series has to be relabelled, and the relabelling is where the minus sign comes from.

Proof(The recurrences from the series)AdaptedAdapted from Riley, Hobson & Bence (2006), §18.5.3, Eqs. (18.92)-(18.95), pp. 611-612; Arfken, Weber & Harris (2013), §14.1, pp. 645-646

We prove the two first-order relations directly from the power series; everything else is algebra.

Step 1: raising the order downwards. Multiply the series by $x^{\nu}$ and differentiate term by term, which is legitimate because the series converges locally uniformly:

\frac{d}{dx}\bigl[x^{\nu}J_\nu(x)\bigr] = \frac{d}{dx}\sum_{m\ge0}\frac{(-1)^{m}\,x^{2\nu+2m}}{m!\,\Gamma(\nu+m+1)\,2^{\nu+2m}} = \sum_{m\ge0}\frac{(-1)^{m}\,(2\nu+2m)\,x^{2\nu+2m-1}}{m!\,\Gamma(\nu+m+1)\,2^{\nu+2m}} .

Use $\Gamma(\nu+m+1) = (\nu+m)\Gamma(\nu+m)$ to cancel the factor $2(\nu+m)$ against the denominator:

= x^{\nu}\sum_{m\ge0}\frac{(-1)^{m}}{m!\,\Gamma(\nu+m)}\left(\frac{x}{2}\right)^{\nu-1+2m} = x^{\nu}J_{\nu-1}(x) ,

since $\Gamma(\nu+m) = \Gamma\bigl((\nu-1)+m+1\bigr)$ is exactly the Gamma factor for order $\nu-1$ .

Step 2: lowering the order. The same computation on $x^{-\nu}J_\nu$ gives

\frac{d}{dx}\bigl[x^{-\nu}J_\nu(x)\bigr] = \sum_{m\ge0}\frac{(-1)^{m}\,2m\,x^{2m-1}}{m!\,\Gamma(\nu+m+1)\,2^{\nu+2m}} .

The $m = 0$ term drops out — this is the one asymmetry between the two computations — so the sum starts at $m = 1$ ; cancelling $2m$ against $m! = m\,(m-1)!$ and relabelling $m = k+1$ ,

= -x^{-\nu}\sum_{k\ge0}\frac{(-1)^{k}}{k!\,\Gamma(\nu+k+2)}\left(\frac{x}{2}\right)^{\nu+1+2k} = -x^{-\nu}J_{\nu+1}(x) ,

because $\Gamma(\nu+k+2) = \Gamma\bigl((\nu+1)+k+1\bigr)$ .

Step 3: the other two relations. Expanding the derivative in Step 1 and dividing by $x^{\nu-1}$ gives $xJ_\nu' + \nu J_\nu = xJ_{\nu-1}$ ; expanding Step 2 and multiplying by $x^{\nu+1}$ gives $xJ_\nu' - \nu J_\nu = -xJ_{\nu+1}$ . Subtract the second from the first and divide by $x$ :

\frac{2\nu}{x}J_\nu = J_{\nu-1}+J_{\nu+1} .

Add them and divide by $x$ : $2J_\nu' = J_{\nu-1}-J_{\nu+1}$ . Finally, at $\nu = 0$ the identity $J_{-1} = -J_1$ turns the last relation into $2J_0' = -2J_1$ . $\square$

The same relations from the generating function. For integer order there is a shorter route, and it is worth recording because it generalises to other generating functions. Differentiating $g(x,t) = e^{(x/2)(t-1/t)}$ in $t$ and in $x$ gives

\frac{\partial g}{\partial t} = \frac{x}{2}\left(1+\frac{1}{t^{2}}\right)g , \qquad \frac{\partial g}{\partial x} = \frac{1}{2}\left(t-\frac{1}{t}\right)g .

Replacing $g$ on both sides by $\sum_n J_n(x)t^{n}$ and equating the coefficient of $t^{n-1}$ in the first and of $t^{n}$ in the second yields $nJ_n = \tfrac{x}{2}(J_{n-1}+J_{n+1})$ and $2J_n' = J_{n-1}-J_{n+1}$ at once ((Riley, Hobson & Bence (2006), §18.5.3, p. 613); (Arfken, Weber & Harris (2013), §14.1, p. 646)).

Remark(Where the hypotheses are used)Original

The $1/\Gamma$ convention, in the definition. The series needs $1/\Gamma(\nu+m+1)$ to have a value for every $m \ge 0$ . Off the negative integers this is automatic; at $\nu = -n$ the first $n$ of those reciprocals are $1/\Gamma$ of a non-positive integer, and the convention adopted in the definition reads them as $0$ rather than as undefined, which is legitimate because $1/\Gamma$ is entire (The Gamma function). What comes out is then $\pm J_{|\nu|}$ — a solution, but not a new one, so the restriction that really bites is on finding a second solution, not on defining the first. Strauss (Strauss (2008), §10.5, p. 283) draws the line in the other place, calling the series undefined at negative integer order and stating his existence claim for every real order except those; the difference is one of convention, and the consequence — no second solution at integer order — is the same either way.
$\nu \notin \mathbb{Z}$ , for a basis from $J_{\pm\nu}$ . This is the Frobenius dichotomy (The Frobenius method): exponents differing by a non-integer give two independent series. For $\nu = n$ the exponents differ by $2n \in \mathbb{Z}$ , the smaller one fails, and the second solution acquires a logarithm. Half-integer $\nu$ is the instructive near-miss: the exponents differ by the odd integer $2\nu$ , so the theorem gives no guarantee, yet the second series survives because the recurrence $n(n-2\nu)a_n = -a_{n-2}$ advances $n$ in steps of two and therefore steps over the obstruction at $n = 2\nu$ ((Riley, Hobson & Bence (2006), §18.5.1, pp. 603-604)). That is why $J_{1/2}$ and $J_{-1/2}$ are independent while $J_1$ and $J_{-1}$ are not.
$a_1 = 0$ . For $\nu \ne -\tfrac12$ the $n=1$ equation forces it. At $\nu = -\tfrac12$ it does not, and $a_1$ is free; choosing $a_1 = 0$ is then a choice, and the series obtained is $J_{-1/2}$ , with the admitted odd part being a multiple of $J_{1/2}$ . Nothing is lost, but the step is a convention, not a deduction.
Integer order, for the generating function. A Laurent series in $t$ has integer powers, so no expansion of this shape can produce non-integer orders. The identity $J_{-n} = (-1)^nJ_n$ is what makes the doubly infinite sum consistent; for non-integer $\nu$ there is no such relation and the generating function is silent. Any statement about $J_\nu$ at general $\nu$ proved "by the generating function" is therefore proved only for integer order, and must be extended by the series (or by analytic continuation in $\nu$ ) before it may be used.
Absolute convergence, in the proof. The rearrangement of the double sum into Laurent coefficients is the one analytic step; it is justified by absolute convergence of both exponential series, which holds for all $x$ and all $t \ne 0$ . At $t = 0$ the left-hand side has an essential singularity, so the Laurent series has no constant-radius disc of convergence around $t=0$ — only the punctured plane.
Term-by-term differentiation. Step 1 and Step 2 differentiate a power series inside its disc of convergence, which here is all of $\mathbb{C}$ ; no separate justification is needed. The loss of the $m=0$ term in Step 2 is not a technicality but the source of the minus sign, and hence of the asymmetry between raising and lowering.
A numerical warning. The three-term recurrence is exact but not safe to iterate upwards in $\nu$ : $J_\nu(x)$ decays very fast once $\nu > x$ , so forward recursion amplifies rounding error into a spurious multiple of the second solution $Y_\nu$ , which grows as $\nu$ increases instead of decaying. Recursion downwards in $\nu$ is stable, which is the basis of the standard algorithms.

Method(Working with

J_\nu

)AdaptedAdapted from Riley, Hobson & Bence (2006), §18.5.3, pp. 611-612; Arfken, Weber & Harris (2013), §14.1, Eqs. (14.10)-(14.11), p. 646

To find the order. Put the equation in the form $x^{2}y''+xy'+(x^{2}-\nu^{2})y = 0$ ; the order is $\nu = \sqrt{\text{(the constant subtracted from }x^{2})}$ . If the independent variable appears as $k\rho$ rather than $x$ , the solutions are $J_\nu(k\rho)$ and the equation reads $\rho^{2}\frac{d^2}{d\rho^2} + \rho\frac{d}{d\rho} + (k^{2}\rho^{2}-\nu^{2})$ .
To choose the solution. If the domain includes $\rho = 0$ and the solution must be bounded there, take $J_\nu$ with $\nu \ge 0$ and discard everything else. Otherwise a second solution is needed (Neumann and Hankel functions).
To evaluate for small or moderate $x$ . Sum the defining series; the terms alternate and fall off like $(x/2)^{2m}/[m!\,\Gamma(\nu+m+1)]$ , so a handful suffice for $x \lesssim \nu+5$ . For large $x$ use the asymptotic forms of Asymptotic forms of Bessel functions instead.
To move between orders. Use $J_{\nu+1} = (2\nu/x)J_\nu - J_{\nu-1}$ to go up, or the same relation solved for $J_{\nu-1}$ to go down; in floating point, go down (see the remark above). For half-integer $\nu$ , start from $J_{1/2} = \sqrt{2/\pi x}\,\sin x$ and $J_{-1/2} = \sqrt{2/\pi x}\,\cos x$ , and every $J_{n+1/2}$ comes out in closed form.
To integrate. Look for the pattern $x^{\nu}J_{\nu-1}$ or $x^{-\nu}J_{\nu+1}$ , which integrate exactly to $x^{\nu}J_\nu$ and $-x^{-\nu}J_\nu$ . A power that does not match is reduced by writing $x^{\mu+1} = x^{\mu-\nu}\cdot x^{\nu+1}$ and integrating by parts.
To expand a plane wave or a phase-modulated signal. Apply the Jacobi-Anger expansion: any $e^{\,i\beta\sin\theta}$ is a Fourier series in $\theta$ whose coefficients are $J_n(\beta)$ .

Example(Closed forms at half-integer order)AdaptedAdapted from Riley, Hobson & Bence (2006), §18.5.1, the reduction of $J_{\pm1/2}$ to trigonometric functions, p. 605, and §18.5.3, the reduction of $J_{\pm3/2}$ by the three-term recurrence, p. 612

(a) Closed forms at half-integer order. The route — collapse $J_{\pm1/2}$ to sine and cosine from the series, then climb the three-term recurrence — is Riley, Hobson and Bence's (Riley, Hobson & Bence (2006), §18.5.1, p. 605; §18.5.3, p. 612); the Gamma-factor identity used to collapse the series, the order $\tfrac52$ , and every number below are ours.

With $\nu = \tfrac12$ the Gamma factors telescope. Using $\Gamma(\tfrac32+m) = \frac{(2m+1)!}{m!\,4^{m}}\cdot\frac{\sqrt{\pi}}{2}$ ,

J_{1/2}(x) = \sum_{m\ge0}\frac{(-1)^{m}}{m!\,\Gamma(m+\tfrac32)}\left(\frac{x}{2}\right)^{1/2+2m} = \sqrt{\frac{2}{\pi x}}\sum_{m\ge0}\frac{(-1)^{m}x^{2m+1}}{(2m+1)!} = \sqrt{\frac{2}{\pi x}}\,\sin x ,

and the same computation at $\nu = -\tfrac12$ gives $J_{-1/2}(x) = \sqrt{2/\pi x}\,\cos x$ . Now climb with $J_{\nu+1} = (2\nu/x)J_\nu - J_{\nu-1}$ :

J_{3/2}(x) = \frac{1}{x}J_{1/2}(x) - J_{-1/2}(x) = \sqrt{\frac{2}{\pi x}}\left(\frac{\sin x}{x} - \cos x\right) ,

J_{5/2}(x) = \frac{3}{x}J_{3/2}(x) - J_{1/2}(x) = \sqrt{\frac{2}{\pi x}}\left[\left(\frac{3}{x^{2}}-1\right)\sin x - \frac{3\cos x}{x}\right] .

Check at $x = 3$ , where $\sin 3 = 0.1411200$ and $\cos 3 = -0.9899925$ . Then $\sqrt{2/3\pi} = 0.4606589$ and

J_{5/2}(3) = 0.4606589\left[\left(\tfrac13 - 1\right)(0.1411200) - (-0.9899925)\right] = 0.4606589 \times 0.8959125 = 0.4127100 .

Summing the defining series at $\nu = \tfrac52$ , $x = 3$ term by term gives $0.4127100$ to seven places, as it must.

Example(An integral, and a modulation index)Original

(b) A definite integral. Evaluate $\int_0^1 x^{2}J_1(x)\,dx$ . The integrand matches $x^{\nu}J_{\nu-1}$ with $\nu = 2$ , so by the first recurrence $x^{2}J_1 = \frac{d}{dx}\bigl[x^{2}J_2\bigr]$ and

\int_0^1 x^{2}J_1(x)\,dx = \bigl[x^{2}J_2(x)\bigr]_0^1 = J_2(1) .

To get a number without a table, use the three-term recurrence at $\nu = 1$ : $J_2(1) = 2J_1(1) - J_0(1)$ . Summing the series at $x = 1$ gives $J_0(1) = 0.7651977$ and $J_1(1) = 0.4400506$ , so

\int_0^1 x^{2}J_1(x)\,dx = 2(0.4400506) - 0.7651977 = \boxed{0.1149035} .

Numerical quadrature of the left-hand side returns $0.11490348$ , agreeing to eight figures.

(c) The modulation index that kills a carrier. A signal of carrier frequency $\omega_c$ whose phase is modulated at frequency $\omega_m$ with peak deviation $\beta$ is $s(t) = \mathrm{Re}\,e^{\,i(\omega_ct + \beta\sin\omega_mt)}$ . By Jacobi-Anger with $\theta = \omega_mt$ ,

e^{\,i\beta\sin\omega_mt} = \sum_{n=-\infty}^{\infty}J_n(\beta)\,e^{\,in\omega_mt} , \qquad\text{so}\qquad s(t) = \sum_{n=-\infty}^{\infty}J_n(\beta)\cos\bigl[(\omega_c+n\omega_m)t\bigr] .

The spectrum is a comb of sidebands spaced by $\omega_m$ , the $n$ -th with amplitude $J_n(\beta)$ . The carrier itself ( $n=0$ ) has amplitude $J_0(\beta)$ , so it disappears exactly when $\beta$ is a zero of $J_0$ : the smallest such $\beta$ is $2.404826$ . Setting $\beta = 2.404826$ and summing $\sum_n J_n(\beta)^{2}$ over $|n| \le 20$ returns $1.0000000$ , as the total power must, with $J_0 = 0$ and the first sidebands carrying $J_{\pm1}(\beta)^{2} = (0.5191475)^{2} = 0.2695141$ each — that is, $53.9\%$ of the power in the first pair of sidebands and none at all in the carrier.

ApplicationOriginal

Anything cylindrical. Separating the Helmholtz equation $\nabla^{2}u + k^{2}u = 0$ in cylindrical polars $(\rho,\varphi,z)$ with $u = P(\rho)e^{im\varphi}e^{i\ell z}$ leaves $\rho^{2}P'' + \rho P' + \bigl[(k^{2}-\ell^{2})\rho^{2}-m^{2}\bigr]P = 0$ , which is Bessel's equation of order $m$ in the variable $\sqrt{k^{2}-\ell^{2}}\,\rho$ . Single-valuedness in $\varphi$ forces $m \in \mathbb{Z}$ , and regularity on the axis forces $P \propto J_m$ . This one calculation underlies the vibrating circular membrane, the cylindrical resonant cavity of a particle accelerator, the acoustic modes of an organ pipe of circular cross-section, the temperature in a cooling cylindrical billet, and the TE and TM modes of a circular waveguide; the radial boundary condition at $\rho = a$ then quantises $k$ through the zeros of $J_m$ or $J_m'$ (Orthogonality and Fourier-Bessel series, Separation of variables in standard coordinates).

Diffraction by a circular aperture. The Fraunhofer amplitude from a uniformly illuminated circular hole of radius $a$ is the two-dimensional Fourier transform of a disc. The angular integral turns into $J_0$ by the integral representation above, and the radial integral is then done exactly by $\int_0^a J_0(br)r\,dr = aJ_1(ab)/b$ — the antiderivative $\frac{d}{dr}[rJ_1(br)] = brJ_0(br)$ . The result is the Airy pattern, intensity $\propto [J_1(ka\sin\alpha)/\sin\alpha]^{2}$ , whose first dark ring sits at $ka\sin\alpha = 3.831706$ , the first zero of $J_1$ . That number is the diffraction limit of every circular lens, telescope and microscope objective, and the same formula with $k$ a particle wavenumber describes the diffraction minima in low-energy nuclear scattering.

Frequency modulation and sidebands. Example (c) is how FM radio, laser phase modulators and optical frequency combs are analysed, and the Bessel-zero trick is a standard laboratory method for calibrating a modulator: increase the drive until the carrier vanishes on a spectrum analyser, and the modulation index is known to be $2.404826$ without measuring the drive at all. The same expansion gives the Debye-Waller-like sidebands in acousto-optic deflection and the Rabi sidebands of a two-level atom driven by a frequency-modulated field.

Hankel transforms and axisymmetric problems. The integral $\int_0^\infty f(\rho)J_\nu(\gamma\rho)\rho\,d\rho$ is to an axisymmetric problem what the Fourier transform is to a translation-invariant one: it diagonalises the radial part of the Laplacian. Potential problems for a charged disc, the stress field under a circular punch, and the diffusion of heat from a circular source are all solved this way (Hankel and Mellin transforms).

Where the order is not an integer. Half-integer orders arise whenever a three-dimensional problem is separated in spherical polars, since the radial factor $r^{-1/2}$ converts the spherical Bessel equation into Bessel's equation of order $\ell+\tfrac12$ (Spherical Bessel functions); order $\tfrac13$ appears in the solution of Airy's equation and so in every turning-point problem of wave optics and quantum mechanics (Airy functions); and continuous orders appear in the normal modes of a hanging chain and in the Kelvin functions of eddy-current theory.

Definition

Neumann and Hankel functions

Adapted

Add to PDF

The Neumann function $Y_\nu = (\cos\nu\pi\,J_\nu - J_{-\nu})/\sin\nu\pi$ , defined at integer order by its limit, is a second solution of Bessel's equation, unbounded at the origin — like $x^{-\nu}$ for $\nu > 0$ , like $\ln x$ at $\nu = 0$ ; the Wronskian $J_\nu Y_\nu' - J_\nu' Y_\nu = 2/\pi x$ is independent of $\nu$ and never vanishes, and the combinations $H^{(1,2)}_\nu = J_\nu \pm iY_\nu$ are the cylindrical analogues of $e^{\pm ix}$ , carrying energy outward and inward.

Prerequisites: Bessel functions of the first kind, Reduction of order, Adjoint and self-adjoint differential operators, Digamma and polygamma functions, The Gamma function

Definition(The Neumann function

Y_\nu

)AdaptedAdapted from Arfken, Weber & Harris (2013), §14.3, Eqs. (14.57) and (14.59)-(14.60), pp. 667-668; Riley, Hobson & Bence (2006), §18.5.2, Eqs. (18.81)-(18.83), pp. 606-607; Hassani (2013), §15.5, Eq. (15.44), p. 483

At non-integer order the two Frobenius solutions $J_\nu$ and $J_{-\nu}$ of Bessel's equation are independent (Bessel functions of the first kind), but at integer order $J_{-n} = (-1)^{n}J_n$ and one solution is missing. The standard repair is not to construct the second solution by reduction of order but to choose a particular combination of $J_{\pm\nu}$ which has a finite limit as $\nu$ approaches an integer.

The Neumann function, or Bessel function of the second kind of order $\nu$ , is

Y_\nu(x) \;=\; \frac{\cos(\nu\pi)\,J_\nu(x) - J_{-\nu}(x)}{\sin(\nu\pi)} \qquad (\nu \notin \mathbb{Z}) ,

and at integer order

Y_n(x) \;=\; \lim_{\nu \to n} Y_\nu(x) .

For $\nu \notin \mathbb{Z}$ this is a linear combination of two solutions, hence a solution, and it is independent of $J_\nu$ because $J_{-\nu}$ enters with the non-zero coefficient $-1/\sin\nu\pi$ . The content of the definition is therefore entirely in the limit: that it exists, and that what it produces is independent of $J_n$ .

At $\nu = n$ both numerator and denominator vanish — $\cos n\pi = (-1)^{n}$ and $J_{-n} = (-1)^{n}J_n$ — so l'Hopital's rule in the variable $\nu$ applies and gives

Y_n(x) = \frac{1}{\pi}\left[\frac{\partial J_\nu}{\partial\nu} - (-1)^{n}\frac{\partial J_{-\nu}}{\partial\nu}\right]_{\nu = n} ,

which, reading $\partial J_{-\nu}/\partial\nu$ as the derivative of $\nu \mapsto J_{-\nu}$ , is the same as $\pi^{-1}\bigl[\partial_\nu J_\nu|_{\nu=n} + (-1)^{n}\partial_\nu J_\nu|_{\nu=-n}\bigr]$ . Differentiating the defining series with respect to $\nu$ brings down a factor $\ln(x/2)$ from $(x/2)^{\nu+2m}$ and a digamma function from $1/\Gamma(\nu+m+1)$ (Digamma and polygamma functions), and the result is a convergent series:

Y_n(x) = \frac{2}{\pi}J_n(x)\ln\frac{x}{2} - \frac{1}{\pi}\sum_{k=0}^{n-1}\frac{(n-k-1)!}{k!}\left(\frac{x}{2}\right)^{2k-n} - \frac{1}{\pi}\sum_{k=0}^{\infty}\frac{(-1)^{k}\bigl[\psi(k+1)+\psi(n+k+1)\bigr]}{k!\,(n+k)!}\left(\frac{x}{2}\right)^{2k+n} .

For $n = 0$ the finite sum is empty and, using $\psi(k+1) = -\gamma + H_k$ with $H_k = \sum_{j\le k}j^{-1}$ the harmonic number and $\gamma$ the Euler-Mascheroni constant, this collapses to

Y_0(x) = \frac{2}{\pi}J_0(x)\left(\gamma + \ln\frac{x}{2}\right) - \frac{2}{\pi}\sum_{k=1}^{\infty}\frac{(-1)^{k}H_k}{(k!)^{2}}\left(\frac{x}{2}\right)^{2k} .

Notation. $Y_\nu$ is the notation of the standard tables and of (Arfken, Weber & Harris (2013), §14.1, p. 644, footnote 1), which also records $N_\nu$ as a common alternative; Strauss (Strauss (2008), §10.5, Eq. (13), p. 285) writes $N_s$ . The names Neumann function and Weber function are also in use for the same object ((Riley, Hobson & Bence (2006), §18.5.2, p. 607)). The handbook uses $Y_\nu$ throughout.

Definition(Behaviour at the origin, and independence from

J_\nu

)AdaptedAdapted from Arfken, Weber & Harris (2013), §14.3, Eqs. (14.58) and (14.61)-(14.62), pp. 667-668; Riley, Hobson & Bence (2006), §18.5.3, p. 608; Strauss (2008), §10.5, Eq. (16), p. 286

As $x \downarrow 0$ ,

Y_\nu(x) = -\frac{\Gamma(\nu)}{\pi}\left(\frac{2}{x}\right)^{\nu} + \cdots \quad (\nu > 0) , \qquad Y_0(x) = \frac{2}{\pi}\left(\ln\frac{x}{2}+\gamma\right) + \cdots .

So $Y_\nu$ is unbounded at the origin for every $\nu \ge 0$ : algebraically for $\nu > 0$ , logarithmically for $\nu = 0$ . Consequently

for $\nu \ge 0$ , $Y_\nu$ cannot be a multiple of $J_\nu$ , which is bounded at the origin, so $\{J_\nu, Y_\nu\}$ is a basis of solutions of Bessel's equation and the general solution is $c_1J_\nu + c_2Y_\nu$ . For non-integer $\nu < 0$ this particular argument separates nothing, since $J_\nu$ is then unbounded at the origin as well; at negative integer order it still works, because $J_{-n} = (-1)^{n}J_n$ and $Y_{-n} = (-1)^{n}Y_n$ reduce such an order to a positive one, where $J_n$ is bounded and $Y_n$ is not; the conclusion nevertheless holds for every $\nu$ , integer or not, because the Wronskian $W[J_\nu, Y_\nu] = 2/\pi x$ of the next block never vanishes ((Riley, Hobson & Bence (2006), §18.5.2, p. 607));
in a problem whose domain contains the axis $x = 0$ and whose solution must be bounded there, $c_2 = 0$ and only $J_\nu$ survives ((Riley, Hobson & Bence (2006), §18.5.3, p. 608)). $Y_\nu$ earns its place when the axis is excluded — an annulus, a coaxial line, the exterior of a cylinder — or when the singular behaviour is physically wanted, as at a line source.

Since $Y_\nu$ is a fixed linear combination of $J_{\pm\nu}$ with coefficients independent of $x$ , it satisfies the same recurrence relations as $J_\nu$ ((Riley, Hobson & Bence (2006), §18.5.3, p. 612); (Arfken, Weber & Harris (2013), §14.3, p. 669)):

Y_{\nu-1}+Y_{\nu+1} = \frac{2\nu}{x}Y_\nu , \qquad Y_{\nu-1}-Y_{\nu+1} = 2Y_\nu' , \qquad Y_0' = -Y_1 .

Two identities at special orders, from (Hassani (2013), §15.5, Problem 15.37, p. 489) and (Riley, Hobson & Bence (2006), §18.5.2, p. 607): $Y_{-n} = (-1)^{n}Y_n$ for $n \in \mathbb{Z}$ , and

Y_{n+1/2}(x) = (-1)^{n+1}J_{-n-1/2}(x) ,

so at half-integer order the Neumann functions are elementary as well. Their singularity at the origin is then a pure power and no logarithm: $J_{-n-1/2}$ has leading term $(x/2)^{-n-1/2}/\Gamma(\tfrac12-n)$ , so $Y_{n+1/2}(x) = O(x^{-n-1/2})$ , exactly the power $x^{-\nu}$ of the display above with $\nu = n+\tfrac12$ . The power does grow with the order — $Y_{1/2} = -\sqrt{2/\pi x}\,\cos x$ is only $O(x^{-1/2})$ , but $Y_{3/2}(x) = -\sqrt{2/\pi x}\,\bigl(x^{-1}\cos x + \sin x\bigr)$ is $O(x^{-3/2})$ , with $Y_{3/2}(0.001) = -2.5231\times10^{4}$ . The special feature of half-integer order is the absence of the $\ln x$ that the integer orders carry, not a bound on the power.

Definition(Wronskians)AdaptedAdapted from Arfken, Weber & Harris (2013), §14.3, Eqs. (14.64)-(14.71), pp. 670-671, and §14.4, Eqs. (14.85)-(14.87), p. 675; Hassani (2013), §15.5, Problem 15.36, p. 489

Write $W[u,v] = uv' - u'v$ . For all $x > 0$ ,

W\bigl[J_\nu, J_{-\nu}\bigr](x) = -\frac{2\sin(\nu\pi)}{\pi x} , \qquad W\bigl[J_\nu, Y_\nu\bigr](x) = \frac{2}{\pi x} , \qquad W\bigl[H^{(1)}_\nu, H^{(2)}_\nu\bigr](x) = -\frac{4i}{\pi x} ,

the second and third holding for every $\nu$ including integers. In particular:

$W[J_\nu, J_{-\nu}]$ vanishes identically exactly when $\nu \in \mathbb{Z}$ — a quantitative statement of the collapse $J_{-n} = (-1)^{n}J_n$ , and the reason a second solution must be manufactured.
$W[J_\nu, Y_\nu] = 2/\pi x$ is never zero and does not depend on $\nu$ . This is the cleanest proof that $\{J_\nu, Y_\nu\}$ is a basis at every order.
Equivalent forms follow from the recurrences, for instance
$J_\nu(x)Y_{\nu+1}(x) - J_{\nu+1}(x)Y_\nu(x) = -\frac{2}{\pi x} , \qquad J_{\nu-1}(x)H^{(1)}_\nu(x) - J_\nu(x)H^{(1)}_{\nu-1}(x) = \frac{2}{i\pi x} .$

Idea of the proofOriginal

Both Wronskian formulas come from one observation and one evaluation. The observation is that Bessel's equation in self-adjoint form is $(xy')' + (x - \nu^{2}/x)y = 0$ , so the coefficient of $y''$ is $p(x) = x$ and Abel's identity gives $W = A/x$ with $A$ a constant: the whole $x$ -dependence is settled before any series is written. The evaluation then needs only the leading term of each series, because multiplying $W$ by $x$ makes the limit at the origin finite, and the Gamma-function reflection formula $\Gamma(1+\nu)\Gamma(1-\nu) = \pi\nu/\sin\nu\pi$ converts what comes out into $\sin\nu\pi$ . Everything else — the value for $(J_\nu,Y_\nu)$ , for the Hankel pair, and the independence statements — is then linear algebra on the definitions.

Proof(The Wronskians)AdaptedAdapted from Arfken, Weber & Harris (2013), §14.3, Eqs. (14.64)-(14.67), pp. 670-671; Hassani (2013), §15.5, Problem 15.34, p. 489

Step 1: the shape of the Wronskian. Divide Bessel's equation by $x$ and rewrite it as

(x\,y')' + \left(x - \frac{\nu^{2}}{x}\right)y = 0 ,

which is self-adjoint with $p(x) = x$ (Adjoint and self-adjoint differential operators). If $u$ and $v$ both solve it, then multiplying $u$ 's equation by $v$ , $v$ 's by $u$ and subtracting gives

v\,(xu')' - u\,(xv')' = 0 , \qquad\text{that is}\qquad \frac{d}{dx}\Bigl[x\bigl(u'v - uv'\bigr)\Bigr] = 0 ,

since the two first-order terms cancel in the derivative. Hence $x\,W[u,v](x)$ is constant: $W[u,v] = A/x$ with $A$ depending on $\nu$ and on the pair, not on $x$ .

Step 2: the constant for $(J_\nu, J_{-\nu})$ , at non-integer $\nu$ . Because $xW$ is constant, we may compute $A = \lim_{x \downarrow 0}xW$ . From the defining series,

J_{\pm\nu}(x) = \frac{1}{\Gamma(1\pm\nu)}\left(\frac{x}{2}\right)^{\pm\nu}\bigl[1 + O(x^{2})\bigr] , \qquad J_{\pm\nu}'(x) = \frac{\pm\nu}{2\,\Gamma(1\pm\nu)}\left(\frac{x}{2}\right)^{\pm\nu-1}\bigl[1 + O(x^{2})\bigr] .

The two products in $W = J_\nu J_{-\nu}' - J_\nu' J_{-\nu}$ have leading terms of the same order $x^{-1}$ , the powers $x^{\nu}$ and $x^{-\nu}$ cancelling:

J_\nu J_{-\nu}' = \frac{-\nu}{x\,\Gamma(1+\nu)\Gamma(1-\nu)}\bigl[1+O(x^{2})\bigr], \qquad J_\nu' J_{-\nu} = \frac{+\nu}{x\,\Gamma(1+\nu)\Gamma(1-\nu)}\bigl[1+O(x^{2})\bigr] .

Subtracting and multiplying by $x$ , the corrections are $O(x^{2})$ and drop out in the limit:

A = \frac{-2\nu}{\Gamma(1+\nu)\Gamma(1-\nu)} .

Now $\Gamma(1+\nu) = \nu\Gamma(\nu)$ and the reflection formula $\Gamma(\nu)\Gamma(1-\nu) = \pi/\sin\nu\pi$ (The Gamma function) give $\Gamma(1+\nu)\Gamma(1-\nu) = \pi\nu/\sin\nu\pi$ , so $A = -2\sin(\nu\pi)/\pi$ and

W[J_\nu,J_{-\nu}](x) = -\frac{2\sin(\nu\pi)}{\pi x} .

Step 3: the constant for $(J_\nu, Y_\nu)$ . The Wronskian is bilinear and antisymmetric, and $W[u,u] = 0$ . So for $\nu \notin \mathbb{Z}$ , using the definition of $Y_\nu$ ,

W[J_\nu, Y_\nu] = \frac{1}{\sin\nu\pi}\Bigl(\cos\nu\pi\,W[J_\nu,J_\nu] - W[J_\nu,J_{-\nu}]\Bigr) = \frac{1}{\sin\nu\pi}\cdot\frac{2\sin\nu\pi}{\pi x} = \frac{2}{\pi x} ,

independent of $\nu$ — the factors $\sin\nu\pi$ cancel exactly, which is the point of the normalisation in the definition.

Step 4: integer order. The right-hand side $2/\pi x$ has no $\nu$ in it, and both $Y_\nu(x)$ and $Y_\nu'(x)$ converge as $\nu \to n$ (the series above converges locally uniformly in $x$ and is continuous in $\nu$ , and $\partial_\nu$ of the defining series of $J_\nu$ is again a locally uniformly convergent series). Letting $\nu \to n$ in Step 3 therefore gives

W[J_n,Y_n](x) = \frac{2}{\pi x} \ne 0 ,

so $Y_n$ exists as a genuine second solution and is independent of $J_n$ . (This is an alternative to reading the independence off the logarithmic term in the series: the series is used only for its convergence, never for the shape of the singularity it exhibits.)

Step 5: the Hankel pair. For constants $a,b,c,d$ the bilinearity and antisymmetry of $W$ give $W[aU+bV,\;cU+dV] = (ad-bc)\,W[U,V]$ , because the $W[U,U]$ and $W[V,V]$ terms vanish. With $U = J_\nu$ , $V = Y_\nu$ and $(a,b,c,d) = (1,i,1,-i)$ we have $ad-bc = -i-i = -2i$ , so

W\bigl[H^{(1)}_\nu,H^{(2)}_\nu\bigr] = -2i\,W[J_\nu,Y_\nu] = -\frac{4i}{\pi x} .

Being non-zero, this confirms that $H^{(1)}_\nu$ and $H^{(2)}_\nu$ are also a basis. $\square$

Definition(Hankel functions)AdaptedAdapted from Arfken, Weber & Harris (2013), §14.4, Eqs. (14.76)-(14.78), pp. 674-675; Riley, Hobson & Bence (2006), §18.5.2, p. 607; Hassani (2013), §15.5, Eq. (15.46), p. 484

The Hankel functions of the first and second kind (also, in older usage, Bessel functions of the third kind) are

H^{(1)}_\nu(x) = J_\nu(x) + i\,Y_\nu(x) , \qquad H^{(2)}_\nu(x) = J_\nu(x) - i\,Y_\nu(x) .

They stand to $J_\nu, Y_\nu$ exactly as $e^{\pm i\theta}$ stands to $\cos\theta, \sin\theta$ ((Arfken, Weber & Harris (2013), §14.4, Eq. (14.78), p. 675)), and the analogy is not merely formal: it is reproduced in their large-argument behaviour. For real $x$ and real $\nu$ they are complex conjugates of one another. The recurrences of the statement above transfer to them unchanged, because each is linear in the function with coefficients that do not depend on which solution is substituted, and $\pm i$ does not depend on $x$ ; so

H_{\nu-1}+H_{\nu+1} = \frac{2\nu}{x}H_\nu , \qquad H_{\nu-1}-H_{\nu+1} = 2H_\nu' ,

for either superscript, and both are singular at the origin, with $H^{(1)}_0(x) \approx 1 + (2i/\pi)(\ln(x/2)+\gamma)$ and $H^{(1)}_\nu(x) \approx -i\Gamma(\nu)(2/x)^{\nu}/\pi$ for $\nu > 0$ ((Arfken, Weber & Harris (2013), §14.4, Eqs. (14.79)-(14.82), p. 675)).

Definition(Outgoing and incoming cylindrical waves)AdaptedAdapted from Arfken, Weber & Harris (2013), §14.6, Eqs. (14.138)-(14.139), p. 692, and Example 14.6.1, pp. 694-695; Strauss (2008), §10.5, Eq. (15), p. 285

As $x \to +\infty$ , with $\nu$ fixed and real,

H^{(1)}_\nu(x) = \sqrt{\frac{2}{\pi x}}\;\exp\!\left[\,i\Bigl(x - \frac{\nu\pi}{2} - \frac{\pi}{4}\Bigr)\right]\bigl[1 + O(x^{-1})\bigr] ,

H^{(2)}_\nu(x) = \sqrt{\frac{2}{\pi x}}\;\exp\!\left[-i\Bigl(x - \frac{\nu\pi}{2} - \frac{\pi}{4}\Bigr)\right]\bigl[1 + O(x^{-1})\bigr] ,

and taking real and imaginary parts recovers $J_\nu \sim \sqrt{2/\pi x}\cos(x - \nu\pi/2 - \pi/4)$ and $Y_\nu \sim \sqrt{2/\pi x}\sin(x - \nu\pi/2 - \pi/4)$ (Asymptotic forms of Bessel functions develops the full expansions). Strauss (Strauss (2008), §10.5, p. 285) identifies this behaviour at large argument as the property that singles the Hankel and Neumann functions out among all the solutions of Bessel's equation.

The wave reading. In a time-harmonic problem with the convention $u(\rho,t) = R(\rho)e^{-i\omega t}$ and $k = \omega/c$ , the choice $R = H^{(1)}_\nu(k\rho)$ gives, at large $k\rho$ ,

u \sim \sqrt{\frac{2}{\pi k\rho}}\;e^{\,i(k\rho - \omega t)}\,e^{-i(\nu\pi/2+\pi/4)} ,

a wave whose surfaces of constant phase $k\rho - \omega t = \text{const}$ move outwards at speed $c$ , with the amplitude falling like $\rho^{-1/2}$ as the energy spreads over a circle of circumference $2\pi\rho$ . $H^{(2)}_\nu(k\rho)e^{-i\omega t}$ is the incoming wave. So:

radiation from a source near the axis is $H^{(1)}$ , with the $e^{-i\omega t}$ convention; this is the boundary condition at infinity, and it replaces the regularity condition at the origin that selects $J_\nu$ in an interior problem ((Arfken, Weber & Harris (2013), §14.6, Example 14.6.1, pp. 694-695));
standing waves are real combinations, $J_\nu$ and $Y_\nu$ ;
with the opposite time convention $e^{+i\omega t}$ the roles of $H^{(1)}$ and $H^{(2)}$ are exchanged. The convention must be stated.

At half-integer order the asymptotic expansions terminate: the series $P_\nu$ and $Q_\nu$ that multiply the trigonometric factors become polynomials in $1/x$ , so the asymptotic solution is an exact one ((Arfken, Weber & Harris (2013), §14.6, p. 694)). Be careful what this licenses. The displays above keep only the leading term, and that term alone is exact only at $\nu = \pm\tfrac12$ , where $P_\nu$ reduces to the constant $1$ and $Q_\nu$ to $0$ , so that the bracket $P_\nu + iQ_\nu$ multiplying the exponential in $H^{(1)}_\nu$ is exactly $1$ : there $H^{(1)}_{1/2}(x) = -i\sqrt{2/\pi x}\,e^{ix}$ with no error term at all. At $\nu = \tfrac32$ the terminating expansion has two terms, and dropping the second leaves a genuine $O(x^{-1})$ relative error — as the closed form $Y_{3/2} = -\sqrt{2/\pi x}(x^{-1}\cos x + \sin x)$ shows directly. Example (a) below checks the exact case at $\nu = \tfrac12$ .

Remark(Where the hypotheses are used)Original

$\nu \notin \mathbb{Z}$ , in the formula for $Y_\nu$ . The denominator $\sin\nu\pi$ is the whole difficulty. For $\nu \notin \mathbb{Z}$ it is harmless and $Y_\nu$ is an ordinary linear combination; at $\nu = n$ the numerator vanishes too, and the definition is the limit, not the formula. Anything proved from the formula must therefore be extended to integer order by a continuity argument — which is exactly what Step 4 of the proof does for the Wronskian, and what fails for naive manipulations like "divide by $\sin\nu\pi$ ".
Self-adjoint form, in Step 1. Abel's identity needs the equation written with $p(x) = x$ ; if one works with the normalised form $y'' + x^{-1}y' + \cdots$ instead, the conclusion is $W = A\exp(-\int x^{-1}) = A/x$ , the same thing. The reason the answer is $1/x$ and not, say, a constant is purely that the first-derivative coefficient is $1/x$ .
$x > 0$ and the limit at the origin, in Step 2. The constant $A$ is extracted at $x \downarrow 0$ , where the series are at their most usable. The step needs the leading behaviour of all four of $J_{\pm\nu}, J_{\pm\nu}'$ , and it needs the two leading terms to cancel in the powers of $x$ but not in the coefficients. The corrections are relatively $O(x^{2})$ , not $O(x)$ , because the series contain only even powers beyond the leading one — so there is room to spare. Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §14.3, Exercise 14.3.6, p. 673) suggest reaching the same conclusion without a limit, by uniqueness of power series: the coefficient of $x^{-1}$ in $W$ is already $A$ , and it remains only to check that the next two coefficients — those of $x^{0}$ and of $x^{1}$ in $W$ itself, equivalently of $x^{1}$ and of $x^{2}$ in the constant $xW$ — vanish.
Why not reduction of order? Reduction of order gives a second solution $J_n(x)\int dx/[xJ_n(x)^{2}]$ ((Hassani (2013), §15.5, Problem 15.34, p. 489)). It is a perfectly good solution, but its scale is arbitrary and it is useless at the zeros of $J_n$ ; the point of the $Y_\nu$ normalisation is that it fixes the scale so that the recurrences, the Wronskian and the asymptotic phase all come out in the form stated, uniformly in $\nu$ .
Fixed real $\nu$ , in the asymptotic statement. The forms above are asymptotic in $x$ at fixed order. They are useless in the regime $\nu \gtrsim x$ . Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §14.6, Eq. (14.145), p. 694) give the honest condition for the leading term alone to be accurate, and they print it at integer order and without a modulus: $8x \gg 4n^{2}-1$ . Read as it stands it is the condition that the first correction be small compared with the leading term, so for the fractional orders $|\nu| < \tfrac12$ , where $4\nu^{2}-1$ is negative, what is wanted is $8x \gg |4\nu^{2}-1|$ ; that modulus is ours, not theirs, and is recorded in the project's discrepancy ledger rather than attributed to them. Uniform-in-order asymptotics is a different and harder problem.
The time convention, in the wave reading. "Outgoing" is a statement about $H^{(1)}_\nu(k\rho)e^{-i\omega t}$ , not about $H^{(1)}$ alone. The function $H^{(1)}_\nu$ knows nothing about time; swapping the convention to $e^{+i\omega t}$ makes $H^{(2)}$ the outgoing solution. A paper that fixes no convention cannot be read.
Negative order needs nothing new. Nothing goes wrong below $\nu = 0$ , provided one does not expect $Y_{-\nu}$ to be $\pm Y_\nu$ . For non-integer $\nu < 0$ the defining formula applies verbatim: it is written for any order whose sine is non-zero, and $-\nu$ is such an order whenever $\nu$ is. What it does not do is reproduce $Y_\nu$ up to a sign. Substituting $-\nu$ for $\nu$ in the definition and using $J_{-\nu}$ , $\cos$ and $\sin$ parities gives $Y_{-\nu} = \sin(\nu\pi)\,J_\nu + \cos(\nu\pi)\,Y_\nu$ ((Hassani (2013), §15.5, Problem 15.37, p. 489)), so $Y_{-\nu}$ is a genuine mixture of the two solutions, degenerating to $\pm Y_\nu$ only when $\sin\nu\pi = 0$ . At integer order that limit is the identity $Y_{-n} = (-1)^{n}Y_n$ , which reduces negative integer order to positive; and the finite sum $\sum_{k<n}$ in the integer-order series is simply empty for $n = 0$ , which is how $Y_0$ comes out purely logarithmic.

Method(Choosing the right pair of solutions)AdaptedAdapted from Arfken, Weber & Harris (2013), §14.3, the list of uses of $Y_\nu$, p. 671, and §14.6, Example 14.6.1, pp. 694-695; Riley, Hobson & Bence (2006), §18.5.3, p. 608

Given a radial equation that reduces to Bessel's equation of order $\nu$ in the variable $k\rho$ :

Is $\rho = 0$ in the domain, and must the solution be bounded there? If yes, write the solution as $c\,J_\nu(k\rho)$ and stop. $Y_\nu$ , $H^{(1)}_\nu$ and $H^{(2)}_\nu$ are all excluded by the singularity.
Is the domain an annulus $a \le \rho \le b$ , or the exterior of a cylinder, with a condition at each end? Use $c_1J_\nu(k\rho) + c_2Y_\nu(k\rho)$ . Imposing the two homogeneous boundary conditions gives a $2\times2$ homogeneous system; a non-trivial solution requires its determinant to vanish, which is the transcendental equation that quantises $k$ . The ratio $c_2/c_1$ then comes from either row.
Does the domain extend to infinity with radiation outwards and no incoming wave? Fix the time convention; with $e^{-i\omega t}$ , take $c\,H^{(1)}_\nu(k\rho)$ and discard $H^{(2)}_\nu$ . If the source is on the axis, the resulting singularity at $\rho = 0$ is the source, not an error.
Is the problem a decaying one rather than an oscillatory one (the equation has $-k^{2}\rho^{2}$ where Bessel's has $+x^{2}$ )? Then the relevant pair is $I_\nu, K_\nu$ : see Modified Bessel functions.
To evaluate $Y_\nu$ , use the series above for small $x$ , the recurrences to move in $\nu$ — upwards, which is the opposite direction to $J_\nu$ and for the same reason — and Asymptotic forms of Bessel functions for large $x$ . Of the two solutions of $f_{\nu+1} = (2\nu/x)f_\nu - f_{\nu-1}$ , $J_\nu$ is the one that decays as $\nu$ grows past $x$ and $Y_\nu$ the one that grows; a three-term recurrence is stable only in the direction in which the wanted solution is the growing one. So $J_\nu$ must be recurred downwards (Bessel functions of the first kind) and $Y_\nu$ upwards: climbing from $Y_0,Y_1$ , the rounding error picked up at each step is a small multiple of $J_\nu$ and shrinks against $Y_\nu$ , whereas descending from large $\nu$ that same error component grows like $J_\nu$ and swamps the answer. The failure is not marginal — in double precision at $x=1$ , descending from exact $Y_{19},Y_{20}$ returns $-1.8\times10^{18}$ for $Y_{10}$ (true value $-1.216\times10^{8}$ ) and $-5.3\times10^{27}$ for $Y_0$ (true value $0.08826$ ), while ascending from exact $Y_0,Y_1$ reproduces $Y_{20}$ to a relative error of $2\times10^{-16}$ . (These figures are ours; none of the sources cited here discusses the stability of the recurrence.) The Wronskian $J_\nu Y_\nu' - J_\nu'Y_\nu = 2/\pi x$ is a cheap and very sharp check on any table or routine.

Example(Two numerical checks)Original

(a) Half-integer order: the asymptotic form made exact. Take $\nu = \tfrac12$ . From the definition, $\cos(\pi/2) = 0$ and $\sin(\pi/2) = 1$ , so

Y_{1/2}(x) = -J_{-1/2}(x) = -\sqrt{\frac{2}{\pi x}}\,\cos x ,

consistent with the identity $Y_{n+1/2} = (-1)^{n+1}J_{-n-1/2}$ at $n = 0$ . Hence

H^{(1)}_{1/2}(x) = J_{1/2}+iY_{1/2} = \sqrt{\frac{2}{\pi x}}\bigl(\sin x - i\cos x\bigr) = -i\sqrt{\frac{2}{\pi x}}\;e^{\,ix} ,

using $\sin x - i\cos x = -i(\cos x + i\sin x)$ . Compare the asymptotic formula at $\nu = \tfrac12$ : the phase is $x - \pi/4 - \pi/4 = x - \pi/2$ , and $e^{\,i(x-\pi/2)} = -ie^{\,ix}$ . The two agree exactly, with no error term — the promised termination of the asymptotic series at half-integer order. Numerically at $x = 3$ , with $\sqrt{2/3\pi} = 0.4606589$ and $\cos 3 = -0.9899925$ : $J_{1/2}(3) = 0.0650082$ and $Y_{1/2}(3) = 0.4560488$ , while $-i\sqrt{2/3\pi}\,e^{3i} = 0.0650082 + 0.4560488\,i$ . The real and imaginary parts match to seven figures.

(b) The Wronskian as a numerical check. At $x = 2$ , tables give $J_0(2) = 0.2238908$ , $J_1(2) = 0.5767248$ , $Y_0(2) = 0.5103757$ , $Y_1(2) = -0.1070324$ . Using $J_0' = -J_1$ and $Y_0' = -Y_1$ ,

J_0(2)Y_0'(2) - J_0'(2)Y_0(2) = -J_0(2)Y_1(2) + J_1(2)Y_0(2) = 0.0239636 + 0.2943463 = 0.3183099 ,

while $2/(2\pi) = 1/\pi = 0.3183099$ . Four tabulated numbers, one exact identity, agreement to seven figures: this is how a Bessel routine is validated.

Example(The cutoff of a coaxial line)AdaptedAdapted from Arfken, Weber & Harris (2013), §14.3, Example 14.3.1 and Eqs. (14.74)-(14.75), pp. 671-672, and Exercise 14.3.10, p. 673

The set-up is Arfken, Weber and Harris's Example 14.3.1 — a coaxial guide, $E_z = c_1J_m + c_2Y_m$ , the two perfect-conductor conditions and the vanishing determinant that quantises the cutoff — and their Exercise 14.3.10 sets the same condition for $b/a = 2$ , written as the vanishing of $J_0(2x)/Y_0(2x) - J_0(x)/Y_0(x)$ , which is our determinant at $b/a = 2$ divided by the factor $-Y_0(2x)Y_0(x)$ . That factor is not everywhere nonzero: it vanishes at each zero of $Y_0(x)$ and of $Y_0(2x)$ , which are exactly the poles of his quotients in his Fig. 14.9. The determinant form is the better behaved of the two there, since it stays finite, and it also stays nonzero unless the two $Y_0$ factors vanish together — at a zero of $Y_0(x)$ it equals $J_0(x)Y_0(2x)$ , and $J_0$ and $Y_0$ have no common zero because their Wronskian is $2/\pi x$ . No zero of $Y_0$ among the first twenty is twice another, so over the range that matters here the two conditions have the same roots. The radii, the numerical work and the mode shape are ours: we take $b/a = 2.5$ and carry the calculation through to a frequency.

A coaxial line has inner conductor of radius $a = 1.0$ cm and outer conductor of radius $b = 2.5$ cm. For the lowest transverse-magnetic mode with no azimuthal variation ( $\nu = 0$ ), the axial electric field is $E_z(\rho) = c_1J_0(\gamma\rho) + c_2Y_0(\gamma\rho)$ , and the perfect-conductor condition $E_z = 0$ at $\rho = a$ and $\rho = b$ gives

\begin{pmatrix} J_0(\gamma a) & Y_0(\gamma a) \\ J_0(\gamma b) & Y_0(\gamma b)\end{pmatrix} \begin{pmatrix}c_1\\c_2\end{pmatrix} = 0 , \qquad\text{so}\qquad J_0(\gamma a)\,Y_0(\gamma b) - J_0(\gamma b)\,Y_0(\gamma a) = 0 .

Writing $z = \gamma a$ and $b/a = 2.5$ , the eigenvalue condition is $F(z) = J_0(z)Y_0(2.5z) - J_0(2.5z)Y_0(z) = 0$ . The radial problem is essentially a string of length $b - a$ , so the first root should be near $\pi a/(b-a) = \pi/1.5 = 2.0944$ ; solving numerically from that start gives

z_1 = 2.073228849 , \qquad z_2 = 4.177298828 .

Check $z_1$ : $J_0(z_1) = 0.181856837$ , $Y_0(z_1) = 0.516714999$ , $J_0(2.5z_1) = -0.116093857$ , $Y_0(2.5z_1) = -0.329860777$ , and

(0.181856837)(-0.329860777) - (-0.116093857)(0.516714999) = -0.059987437 + 0.059987437 = 0 ,

both products agreeing to nine figures. The mode shape is $E_z \propto J_0(\gamma\rho) - 0.351948051\,Y_0(\gamma\rho)$ , the ratio being $-J_0(z_1)/Y_0(z_1)$ ; the same ratio computed from the outer row, $-J_0(2.5z_1)/Y_0(2.5z_1) = -0.351948051$ , agrees, which is the consistency check that the determinant vanished.

Finally the cutoff frequency. With $\gamma_1 = z_1/a = 2.0732288\ \mathrm{cm^{-1}}$ and $c = 2.99792458\times10^{10}\ \mathrm{cm\,s^{-1}}$ ,

f_c = \frac{c\,\gamma_1}{2\pi} = \frac{(2.99792458\times10^{10})(2.0732288)}{2\pi}\ \mathrm{s^{-1}} = 9.892\ \mathrm{GHz} .

Below $9.892$ GHz this TM mode cannot propagate. Note that the cutoff is a statement about $Y_0$ as much as about $J_0$ : had we kept only $J_0$ , as regularity on the axis would demand in a hollow pipe, there would be no second constant with which to satisfy the second boundary condition and no mode at all.

ApplicationOriginal

Coaxial lines and annular cavities. The coaxial-line example above is the design calculation for a coaxial cable or a coaxial resonator: the TEM mode of a coaxial line has no cutoff, which is why such cables carry signals from DC upwards, but the first higher mode does, and the usable bandwidth of the cable ends there. The same determinant condition, with $J_\nu'$ and $Y_\nu'$ in place of $J_\nu$ and $Y_\nu$ , gives the TE modes and the acoustic modes of an annular duct; with $\nu = m$ it gives the modes with $m$ azimuthal nodal lines.

Line sources and the two-dimensional Helmholtz Green's function. The outgoing solution of $(\nabla^{2}+k^{2})G = -\delta^{(2)}(\boldsymbol{\rho})$ in the plane is $G = \tfrac{i}{4}H^{(1)}_0(k\rho)$ . Both features of $H^{(1)}_0$ are needed: the logarithmic singularity at $\rho = 0$ reproduces the delta function (it matches the two-dimensional Laplace Green's function $-\ln\rho/2\pi$ near the source), and the outgoing phase at large $\rho$ enforces the Sommerfeld radiation condition. This is the kernel of every two-dimensional scattering calculation and of the boundary-element method in acoustics and electromagnetics, and $H^{(1)}_0$ likewise appears in the Hankel-transform solution of axisymmetric problems on $0 < \rho < \infty$ , where it supplies the Green's function of the singular Bessel operator ((Stakgold & Holst (2011), §7.3, Eq. (7.3.33), p. 450)).

Scattering off a cylinder. Expand the incident plane wave in $J_m(k\rho)e^{im\varphi}$ (Bessel functions of the first kind) and the scattered field in $H^{(1)}_m(k\rho)e^{im\varphi}$ : one is regular everywhere, the other purely outgoing, and the boundary condition on the cylinder surface fixes one coefficient per $m$ . The resulting partial-wave series is the textbook solution for radar cross-sections of wires, for the scattering of sound by a pipe, and — with $\nu$ half-integral and spherical Bessel functions — for Mie scattering off a sphere (Spherical Bessel functions).

Where the singularity is the physics. In a dislocation core, at a crack tip, along a vortex line in a superfluid or at a current-carrying filament, the field genuinely diverges on an axis, and the solution one wants is the singular one. $Y_\nu$ and the Hankel functions are the admissible solutions there, and the coefficient of the singular part is the strength of the source — the Burgers vector, the stress-intensity factor, the circulation, the current.

Quantum scattering and the phase shift. For a radial Schrodinger problem the pair $(J,Y)$ or $(H^{(1)},H^{(2)})$ at large $r$ defines the scattering phase shift: the exterior solution is a combination $\cos\delta_\ell\,j_\ell - \sin\delta_\ell\,y_\ell$ , and $\delta_\ell$ is read off by matching to the interior solution. The constancy of the Wronskian is what makes this matching well posed, and the same quantity appears as the flux normalisation of the continuum states.

Theorem

Orthogonality and Fourier-Bessel series

Adapted

Add to PDF

For fixed order $\nu > -1$ the functions $J_\nu(\alpha_{\nu j}r/a)$ , with $\alpha_{\nu j}$ the positive zeros of $J_\nu$ , are orthogonal on $[0,a]$ with weight $r$ and have squared norm $\tfrac12a^{2}J_{\nu+1}(\alpha_{\nu j})^{2}$ , so any reasonable function on the disc expands in a Fourier-Bessel series with explicitly computable coefficients.

Prerequisites: Bessel functions of the first kind, Sturm-Liouville eigenvalue problems, Adjoint and self-adjoint differential operators, Eigenfunction expansions, Sturm comparison and oscillation theorems

Definition(The zeros of

J_\nu

)AdaptedAdapted from Strauss (2008), §10.5, pp. 283-284; Arfken, Weber & Harris (2013), §14.1, Table 14.1, p. 648, Exercises 14.1.11 and 14.2.4, pp. 656 and 665

Fix $\nu > -1$ and real. The function $J_\nu$ has infinitely many positive zeros

0 < \alpha_{\nu1} < \alpha_{\nu2} < \alpha_{\nu3} < \cdots \longrightarrow \infty ,

each of them simple ( $J_\nu'(\alpha_{\nu j}) \ne 0$ ), and no complex zeros at all ((Arfken, Weber & Harris (2013), §14.2, Exercise 14.2.4, p. 665)). We write $\beta_{\nu j}$ for the $j$ -th positive zero of $J_\nu'$ . Three further facts, all from Strauss (Strauss (2008), §10.5, p. 283) except where noted:

Interlacing. Between two consecutive zeros of $J_\nu$ there lies a zero of $J_{\nu+1}$ , and conversely, so the two sets of zeros separate each other. This is Strauss's (Strauss (2008), §10.5, p. 283) statement, made for every real order; Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §14.1, Exercise 14.1.11, p. 656) set the sharper form — one and only one zero in between — as an exercise, but only at integer order $n$ . More generally the zeros of any two independent solutions of Bessel's equation separate each other, which is Sturm separation applied to this equation.
Spacing. Consecutive zeros approach a separation of $\pi$ . This follows from the leading asymptotic form $J_\nu(x) \sim \sqrt{2/\pi x}\cos\bigl(x - \nu\pi/2 - \pi/4\bigr)$ ((Arfken, Weber & Harris (2013), §14.6, Eq. (14.144), p. 694); Asymptotic forms of Bessel functions), which also gives the useful estimate $\alpha_{\nu j} \;\approx\; \left(j + \frac{\nu}{2} - \frac14\right)\pi .$
No closed formula. The zeros are not expressible in elementary terms and must be found numerically; they are tabulated in (Arfken, Weber & Harris (2013), §14.1, Table 14.1, p. 648).

The first few, to seven decimals:

$j$	$\alpha_{0j}$	$\alpha_{1j}$	$\alpha_{2j}$	$\beta_{1j}$
1	2.4048256	3.8317060	5.1356223	1.8411838
2	5.5200781	7.0155867	8.4172441	5.3314428
3	8.6537279	10.1734681	11.6198412	8.5363164
4	11.7915344	13.3236919	14.7959518	11.7060049

The estimate above gives $\alpha_{03} \approx 2.75\pi = 8.6393798$ against the true $8.6537279$ , an error of $0.17\%$ ; by $j=4$ the error is $0.09\%$ . Two consistency checks on the table: since $J_0' = -J_1$ , the zeros of $J_0'$ are $0$ together with the $\alpha_{1j}$ , so the $\beta_{0j}$ need no separate column; and interlacing requires $\alpha_{0j} < \alpha_{1j} < \alpha_{0,j+1}$ , which the first two columns satisfy at every row.

Theorem(Orthogonality and normalisation)AdaptedAdapted from Arfken, Weber & Harris (2013), §14.2, Eqs. (14.41)-(14.46), pp. 661-663, with Exercises 14.2.1-14.2.3, p. 665; Riley, Hobson & Bence (2006), §18.5.3, Eqs. (18.84)-(18.89), pp. 608-610; Strauss (2008), §10.5, Eqs. (8)-(9), p. 284

Fix $\nu > -1$ and $a > 0$ .

(i) The Lommel identity. For any $\alpha, \beta > 0$ with $\alpha \ne \beta$ ,

\bigl(\alpha^{2}-\beta^{2}\bigr)\int_0^{a} J_\nu(\alpha r)\,J_\nu(\beta r)\,r\,dr = a\Bigl[\beta\,J_\nu(\alpha a)J_\nu'(\beta a) - \alpha\,J_\nu(\beta a)J_\nu'(\alpha a)\Bigr] ,

where $J_\nu'$ is the derivative with respect to the whole argument. The lower endpoint contributes nothing.

(ii) Orthogonality. Consequently the right-hand side, and hence the integral, vanishes whenever the boundary data at $r = a$ are such that the bracket is zero. The two standard cases are

Dirichlet: $\alpha = \alpha_{\nu i}/a$ and $\beta = \alpha_{\nu j}/a$ , so that $J_\nu(\alpha a) = J_\nu(\beta a) = 0$ ;
Neumann: $\alpha = \beta_{\nu i}/a$ and $\beta = \beta_{\nu j}/a$ , so that $J_\nu'(\alpha a) = J_\nu'(\beta a) = 0$ .

In either case, for $i \ne j$ ,

\int_0^{a} J_\nu\!\left(\lambda_{\nu i}\frac{r}{a}\right)J_\nu\!\left(\lambda_{\nu j}\frac{r}{a}\right)r\,dr = 0 , \qquad \lambda = \alpha \text{ or } \beta .

(iii) Normalisation. For every $A > 0$ ,

2\int_0^{A} J_\nu(x)^{2}\,x\,dx = \bigl[A\,J_\nu'(A)\bigr]^{2} + \bigl(A^{2}-\nu^{2}\bigr)J_\nu(A)^{2} ,

and hence, after rescaling,

\int_0^{a} J_\nu\!\left(\alpha_{\nu j}\frac{r}{a}\right)^{2}r\,dr = \frac{a^{2}}{2}J_\nu'(\alpha_{\nu j})^{2} = \frac{a^{2}}{2}J_{\nu+1}(\alpha_{\nu j})^{2}

in the Dirichlet case, and

\int_0^{a} J_\nu\!\left(\beta_{\nu j}\frac{r}{a}\right)^{2}r\,dr = \frac{a^{2}}{2}\left(1 - \frac{\nu^{2}}{\beta_{\nu j}^{2}}\right)J_\nu(\beta_{\nu j})^{2}

in the Neumann case.

Where this sits. Dividing Bessel's equation in the variable $kr$ by $r$ exhibits it as the Sturm-Liouville problem

-\frac{d}{dr}\left(r\frac{du}{dr}\right) + \frac{\nu^{2}}{r}u = k^{2}\,r\,u ,

with $p(r) = r$ , $q(r) = \nu^{2}/r$ , weight $w(r) = r$ and eigenvalue $k^{2}$ . The weight $r$ is not an artifice: on a disc of radius $a$ it is exactly the Jacobian of plane polar coordinates, so this orthogonality is orthogonality in $L^{2}$ of the disc ((Strauss (2008), §10.2, Eq. (23), p. 269)). But the problem is singular: $p(0) = 0$ and $q$ blows up at the origin, so the regular Sturm-Liouville theorem does not apply and conclusions have to be re-earned. Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §18.5.3, p. 608) make the point sharply: because $p(x) = x$ , nothing in the equation itself picks out an interval on which solutions for different eigenvalues are automatically orthogonal; the interval and the boundary conditions must be supplied, and $[0,a]$ works only because the $r = 0$ endpoint takes care of itself.

Idea of the proofOriginal

Both halves are the same trick used twice. For orthogonality, write down the equation for $J_\nu(\alpha r)$ and the equation for $J_\nu(\beta r)$ , multiply each by the other function and subtract: the second-derivative and first-derivative terms assemble into the exact derivative $\frac{d}{dr}\bigl[r(fg'-gf')\bigr]$ , and what is left on the other side is $(\alpha^{2}-\beta^{2})rfg$ . Integrating turns a statement about an integral into a statement about one boundary value. For the normalisation the two functions coincide, so this identity degenerates to $0 = 0$ ; instead multiply the self-adjoint equation by $2ru'$ , which again produces an exact derivative, this time of $(ru')^{2} + (r^{2}-\nu^{2})u^{2}$ . In both cases the endpoint $r = 0$ contributes nothing, because the factor $r$ in $p$ suppresses it.

Proof(Orthogonality and the Lommel identity)AdaptedAdapted from Riley, Hobson & Bence (2006), §18.5.3, Eqs. (18.85)-(18.88), pp. 609-610; Arfken, Weber & Harris (2013), §14.2, Eq. (14.43), p. 661

Step 1: the two equations. Put $f(r) = J_\nu(\alpha r)$ and $g(r) = J_\nu(\beta r)$ . If $w = \alpha r$ then $f' = \alpha J_\nu'(w)$ and $f'' = \alpha^{2}J_\nu''(w)$ , so

r^{2}f'' + rf' + (\alpha^{2}r^{2}-\nu^{2})f = w^{2}J_\nu''(w) + wJ_\nu'(w) + (w^{2}-\nu^{2})J_\nu(w) = 0

by Bessel's equation itself; likewise $r^{2}g'' + rg' + (\beta^{2}r^{2}-\nu^{2})g = 0$ . The rescaling is the only place the parameters $\alpha,\beta$ enter.

Step 2: subtract, and recognise a derivative. Multiply the $f$ -equation by $g$ , the $g$ -equation by $f$ , subtract, and divide by $r$ :

r\bigl(f''g - g''f\bigr) + \bigl(f'g - g'f\bigr) + \bigl(\alpha^{2}-\beta^{2}\bigr)r\,fg = 0 .

The first two groups are exactly $-\frac{d}{dr}\bigl[r(fg'-gf')\bigr]$ , since

\frac{d}{dr}\Bigl[r\bigl(fg'-gf'\bigr)\Bigr] = r\bigl(fg''-gf''\bigr) + \bigl(fg'-gf'\bigr) ,

the terms $f'g'$ cancelling. Hence

\frac{d}{dr}\Bigl[r\bigl(fg'-gf'\bigr)\Bigr] = \bigl(\alpha^{2}-\beta^{2}\bigr)\,r\,fg .

Step 3: integrate from $0$ to $a$ . Integrating,

\bigl(\alpha^{2}-\beta^{2}\bigr)\int_0^{a}r\,f g\,dr = \Bigl[r\bigl(fg'-gf'\bigr)\Bigr]_0^{a} .

At $r = a$ the bracket is $a\bigl[\beta J_\nu(\alpha a)J_\nu'(\beta a) - \alpha J_\nu(\beta a)J_\nu'(\alpha a)\bigr]$ , which is the stated right-hand side.

Step 4: the endpoint $r = 0$ gives nothing, for $\nu > -1$ . Near $0$ the series give $f(r) = c_\alpha r^{\nu}\bigl[1+O(r^{2})\bigr]$ with $c_\alpha = (\alpha/2)^{\nu}/\Gamma(1+\nu)$ , and therefore $f'(r) = \nu f(r)/r\bigl[1+O(r^{2})\bigr]$ , and the same for $g$ . Hence

fg' - gf' = \frac{\nu}{r}fg - \frac{\nu}{r}gf + O\bigl(r^{2\nu+1}\bigr) = O\bigl(r^{2\nu+1}\bigr) ,

the leading terms cancelling identically. Multiplying by $r$ gives $O(r^{2\nu+2})$ , which tends to $0$ precisely when $\nu > -1$ . (For $\nu \ge 0$ no cancellation is even needed: each term separately vanishes. The cancellation is what rescues the range $-1 < \nu < 0$ , where $f$ and $g$ are individually unbounded; see (Arfken, Weber & Harris (2013), §14.2, Exercise 14.2.3, p. 665).) The integral on the left converges for the same reason: $r\,fg \sim r^{2\nu+1}$ is integrable at $0$ exactly when $\nu > -1$ .

Step 5: the two boundary conditions. If $\alpha a$ and $\beta a$ are both zeros of $J_\nu$ , both terms in the bracket at $r=a$ contain a vanishing factor $J_\nu(\alpha a)$ or $J_\nu(\beta a)$ . If instead both are zeros of $J_\nu'$ , both terms contain a vanishing $J_\nu'$ . Either way the right-hand side is zero, and since $\alpha^{2} \ne \beta^{2}$ we may divide:

\int_0^{a}J_\nu(\alpha r)J_\nu(\beta r)\,r\,dr = 0 . \qquad\square

(More general homogeneous conditions $\kappa J_\nu(\lambda a) + \lambda J_\nu'(\lambda a) = 0$ with a fixed $\kappa$ also annihilate the bracket, giving the so-called Dini series; the computation is the same.)

Proof(The normalisation)AdaptedAdapted from Strauss (2008), §10.5, Eqs. (8)-(9), p. 284; Arfken, Weber & Harris (2013), §14.2, Eqs. (14.45)-(14.46), pp. 662-663; Hassani (2013), §15.6, Problem 15.39, p. 490

Step 1: an exact derivative. Let $u(x) = J_\nu(x)$ , which satisfies the self-adjoint form $(xu')' + \bigl(x - \nu^{2}/x\bigr)u = 0$ . Multiply by $2xu'$ and read the two pieces as derivatives:

2(xu')(xu')' = \frac{d}{dx}\bigl[(xu')^{2}\bigr] , \qquad 2xu'\left(x-\frac{\nu^{2}}{x}\right)u = 2\bigl(x^{2}-\nu^{2}\bigr)uu' = \frac{d}{dx}\Bigl[\bigl(x^{2}-\nu^{2}\bigr)u^{2}\Bigr] - 2xu^{2} .

Adding, the equation becomes

\frac{d}{dx}\Bigl[(xu')^{2} + \bigl(x^{2}-\nu^{2}\bigr)u^{2}\Bigr] = 2x\,u^{2} .

Step 2: integrate. Integrating from $0$ to $A$ ,

2\int_0^{A}x\,J_\nu(x)^{2}\,dx = \Bigl[(xJ_\nu'(x))^{2} + (x^{2}-\nu^{2})J_\nu(x)^{2}\Bigr]_0^{A} .

At the lower limit, $J_\nu(x) = cx^{\nu}[1+O(x^{2})]$ and $xJ_\nu'(x) = \nu cx^{\nu}[1+O(x^{2})]$ , so the bracket is $\bigl(\nu^{2} - \nu^{2}\bigr)c^{2}x^{2\nu} + O(x^{2\nu+2}) = O(x^{2\nu+2})$ , which vanishes for $\nu > -1$ — again by cancellation of the two leading terms, not by each vanishing separately. This gives part (iii) of the statement.

Step 3: specialise to a zero. Take $A = \alpha_{\nu j}$ , so $J_\nu(A) = 0$ and only the first term survives:

\int_0^{\alpha_{\nu j}}x\,J_\nu(x)^{2}\,dx = \tfrac12\,\alpha_{\nu j}^{2}\,J_\nu'(\alpha_{\nu j})^{2} .

Step 4: rescale to $[0,a]$ . Substituting $x = \alpha_{\nu j}r/a$ , so $x\,dx = (\alpha_{\nu j}/a)^{2}r\,dr$ ,

\int_0^{a}J_\nu\!\left(\alpha_{\nu j}\frac{r}{a}\right)^{2}r\,dr = \left(\frac{a}{\alpha_{\nu j}}\right)^{2}\int_0^{\alpha_{\nu j}}x\,J_\nu(x)^{2}\,dx = \frac{a^{2}}{2}J_\nu'(\alpha_{\nu j})^{2} .

Step 5: replace $J_\nu'$ by $J_{\nu+1}$ . The recurrence $xJ_\nu'(x) - \nu J_\nu(x) = -xJ_{\nu+1}(x)$ of Bessel functions of the first kind, evaluated at $x = \alpha_{\nu j}$ where $J_\nu$ vanishes, gives $J_\nu'(\alpha_{\nu j}) = -J_{\nu+1}(\alpha_{\nu j})$ . Squaring removes the sign, and the Dirichlet normalisation follows.

Step 6: the Neumann case. Take instead $A = \beta_{\nu j}$ , a zero of $J_\nu'$ . Then the first term of the bracket in Step 2 vanishes and the second survives:

\int_0^{\beta_{\nu j}}x\,J_\nu(x)^{2}\,dx = \tfrac12\bigl(\beta_{\nu j}^{2}-\nu^{2}\bigr)J_\nu(\beta_{\nu j})^{2} ,

and the same rescaling gives $\tfrac12 a^{2}\bigl(1-\nu^{2}/\beta_{\nu j}^{2}\bigr)J_\nu(\beta_{\nu j})^{2}$ , which is Arfken, Weber and Harris's (Arfken, Weber & Harris (2013), §14.2, Exercise 14.2.2(b), p. 665) answer. $\square$

Remark(Where the hypotheses are used)Original

$\nu > -1$ . This single condition does three jobs: it makes $r J_\nu(\alpha r)^{2} \sim r^{2\nu+1}$ integrable at the origin, it makes the boundary term at $r = 0$ vanish, and it is the range in which $J_\nu$ has only real zeros. What fails below it is the first job: for non-integer $\nu$ the integrand near the origin is $\Theta(r^{2\nu+1})$ , and $\int_0 r^{2\nu+1}dr$ diverges once $2\nu+1 \le -1$ , logarithmically in the limiting case. So the condition is sharp approached through non-integer orders — it is $\nu \downarrow -1$ that breaks the integral, not $\nu = -1$ itself. The endpoint is in fact fine, which is why Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §14.2, p. 661, footnote 3) state the vanishing of the lower boundary term for all $\nu \ge -1$ : at $\nu = -1$ the $1/\Gamma$ convention gives $J_{-1} = -J_1$ (Bessel functions of the first kind), the leading power is $r^{3}$ and not $r^{-1}$ , and the formulas of this unit hold verbatim — for instance $\int_0^{1}rJ_{-1}(\alpha_{11}r)^{2}dr = 0.0811075654 = \tfrac12J_0(\alpha_{11})^{2}$ , which is exactly $\tfrac12J_{\nu+1}(\alpha_{\nu 1})^{2}$ at $\nu = -1$ . The same holds at every negative integer order, since $J_{-n} = (-1)^{n}J_n$ there; the unit states the hypothesis as $\nu > -1$ because that is the range in which the proof below runs without case distinctions. For $-1 < \nu < 0$ the result survives only because of a cancellation between two individually divergent boundary terms — a case worth checking explicitly rather than waving at, which is why (Arfken, Weber & Harris (2013), §14.2, Exercise 14.2.3, p. 665) sets it as an exercise.
Same order $\nu$ , different scale. The orthogonal family is $\{J_\nu(\alpha_{\nu j}r/a)\}_j$ at fixed $\nu$ . There is no orthogonality between $J_\nu$ and $J_\mu$ for $\nu \ne \mu$ ; in the two-dimensional problem the different orders are separated instead by the orthogonality of $e^{im\varphi}$ in the angular variable. Mixing the two kinds of index is the commonest error here.
The lower limit must be $0$ . On an annulus $[b,a]$ with $b > 0$ , Step 4 of the first proof is unavailable and the boundary term at $r = b$ must be killed by a condition there. The eigenfunctions are then combinations of $J_\nu$ and $Y_\nu$ (Neumann and Hankel functions) and the eigenvalues solve a transcendental determinant condition, not $J_\nu(\lambda) = 0$ .
Dirichlet or Neumann, consistently. The bracket in the Lommel identity vanishes when both $\alpha a$ and $\beta a$ satisfy the same condition. A set mixing zeros of $J_\nu$ with zeros of $J_\nu'$ is not orthogonal.
Simplicity of the zeros, in the normalisation. Step 3 divides by nothing, but Step 5 would be vacuous if $J_\nu'(\alpha_{\nu j})$ were zero: the normalisation constant would vanish and the expansion would be meaningless. Simplicity of the zeros is therefore not decoration. It follows from uniqueness for the initial value problem: $J_\nu(\alpha) = J_\nu'(\alpha) = 0$ at an interior point $\alpha > 0$ would force $J_\nu \equiv 0$ .
Completeness is not proved here. Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §14.2, p. 663) preface the Bessel series by assuming, in so many words and without argument, that the family is complete; we follow them, and label the hypothesis as a hypothesis in the statement above. The justification belongs to the theory of singular Sturm-Liouville problems (Eigenfunction expansions); orthogonality and the norm, which are what this unit proves, are completely elementary by comparison.

Theorem(Fourier-Bessel series)AdaptedAdapted from Arfken, Weber & Harris (2013), §14.2, Eqs. (14.47)-(14.48), p. 663, and Exercise 14.2.7, p. 666; Riley, Hobson & Bence (2006), §18.5.3, Eqs. (18.90)-(18.91), pp. 610-611

Fix $\nu > -1$ and $a > 0$ , and assume the set $\bigl\{J_\nu(\alpha_{\nu j}r/a)\bigr\}_{j\ge1}$ is complete in the weighted space $L^{2}\bigl([0,a],\,r\,dr\bigr)$ of complex-valued functions, with the inner product $(f,g) = \int_0^{a}f(r)\overline{g(r)}\,r\,dr$ . The basis functions are real for real $\nu$ , so the conjugate plays no part in the coefficient formula; it does appear in Parseval's relation, and both are written below with moduli so that complex $f$ is covered. Then every $f$ in that space has the Fourier-Bessel expansion

f(r) = \sum_{j=1}^{\infty}c_{j}\,J_\nu\!\left(\alpha_{\nu j}\frac{r}{a}\right) , \qquad c_{j} = \frac{2}{a^{2}\,J_{\nu+1}(\alpha_{\nu j})^{2}}\int_0^{a}f(r)\,J_\nu\!\left(\alpha_{\nu j}\frac{r}{a}\right)r\,dr ,

the series converging in mean square with weight $r$ . Parseval's relation reads

\int_0^{a}\bigl|f(r)\bigr|^{2}\,r\,dr = \frac{a^{2}}{2}\sum_{j=1}^{\infty}\bigl|c_{j}\bigr|^{2}\,J_{\nu+1}(\alpha_{\nu j})^{2}

((Arfken, Weber & Harris (2013), §14.2, Exercise 14.2.7, p. 666)). The coefficient formula is obtained exactly as for any orthogonal expansion: multiply by $J_\nu(\alpha_{\nu i}r/a)\,r$ , integrate, and use orthogonality to kill every term but $i = j$ and the normalisation to evaluate that one ((Riley, Hobson & Bence (2006), §18.5.3, p. 611)).

Replacing $\alpha_{\nu j}$ by $\beta_{\nu j}$ throughout, and the norm by its Neumann value, gives the companion expansion appropriate to a flux-free boundary.

Method(Expanding a function in a Fourier-Bessel series)AdaptedAdapted from Arfken, Weber & Harris (2013), §14.2, Eqs. (14.47)-(14.48), p. 663; Riley, Hobson & Bence (2006), §18.5.3, pp. 610-611

Fix the order. $\nu$ is set by the problem — by the azimuthal index $m$ in a disc problem, by $\ell+\tfrac12$ in a spherical one — not by the function being expanded.
Fix the boundary condition at $r = a$ , and with it the eigenvalues: $\lambda_{j} = \alpha_{\nu j}$ for $f(a) = 0$ , $\lambda_{j} = \beta_{\nu j}$ for $f'(a) = 0$ . Look the zeros up or compute them from the estimate $(j+\nu/2-\tfrac14)\pi$ refined by Newton's method, using $J_\nu' = J_{\nu-1} - \nu J_\nu/x$ .
Compute the overlap integrals $\int_0^{a}f(r)J_\nu(\lambda_{j}r/a)\,r\,dr$ . Where $f$ is a power of $r$ , use the antiderivatives $\int x^{\nu}J_{\nu-1}dx = x^{\nu}J_\nu$ and $\int x^{-\nu}J_{\nu+1}dx = -x^{-\nu}J_\nu$ of Bessel functions of the first kind, reducing higher powers by parts.
Divide by the norm, $\tfrac12a^{2}J_{\nu+1}(\alpha_{\nu j})^{2}$ or $\tfrac12a^{2}(1-\nu^{2}/\beta_{\nu j}^{2})J_\nu(\beta_{\nu j})^{2}$ .
Check with Parseval. Compute $\int_0^{a}f^{2}r\,dr$ independently and compare with the sum of squared coefficients times their norms. This catches a wrong normalisation immediately, which a few numerical partial sums may not.
Expect slow convergence if $f$ does not satisfy the boundary condition. The partial sums converge to the boundary value imposed by the eigenfunctions, so forcing $f(a) \ne 0$ into a Dirichlet Bessel series produces Gibbs-like behaviour near $r = a$ .

Example(A Fourier-Bessel series, checked three ways)Original

Expand $f(r) = 1-r^{2}$ on $0 \le r \le 1$ in Bessel functions of order zero with the Dirichlet condition, that is, in $\bigl\{J_0(\alpha_{j}r)\bigr\}$ with $\alpha_{j} = \alpha_{0j}$ . Note first that $f(1) = 0$ , so $f$ obeys the boundary condition and good convergence is expected.

Step 1: the overlap integral. Put $u = \alpha r$ , $\alpha = \alpha_{j}$ . Then

\int_0^{1}(1-r^{2})J_0(\alpha r)\,r\,dr = \frac{1}{\alpha^{2}}\int_0^{\alpha}\left(1-\frac{u^{2}}{\alpha^{2}}\right)J_0(u)\,u\,du .

For the first piece, $\frac{d}{du}[uJ_1(u)] = uJ_0(u)$ , so $\int_0^{\alpha}uJ_0\,du = \alpha J_1(\alpha)$ . For the second, integrate by parts against the same antiderivative and then use $\frac{d}{du}[u^{2}J_2] = u^{2}J_1$ :

\int_0^{\alpha}u^{3}J_0\,du = \Bigl[u^{2}\cdot uJ_1\Bigr]_0^{\alpha} - 2\int_0^{\alpha}u^{2}J_1\,du = \alpha^{3}J_1(\alpha) - 2\alpha^{2}J_2(\alpha) .

Hence

\int_0^{1}(1-r^{2})J_0(\alpha r)\,r\,dr = \frac{1}{\alpha^{2}}\left[\alpha J_1(\alpha) - \frac{1}{\alpha^{2}}\Bigl(\alpha^{3}J_1(\alpha) - 2\alpha^{2}J_2(\alpha)\Bigr)\right] = \frac{1}{\alpha^{2}}\Bigl[\alpha J_1(\alpha) - \alpha J_1(\alpha) + 2J_2(\alpha)\Bigr] = \frac{2J_2(\alpha)}{\alpha^{2}} .

Since $\alpha$ is a zero of $J_0$ , the three-term recurrence gives $J_2(\alpha) = (2/\alpha)J_1(\alpha) - J_0(\alpha) = (2/\alpha)J_1(\alpha)$ , so the integral is $4J_1(\alpha)/\alpha^{3}$ .

Step 2: the coefficients. With $a = 1$ , $\nu = 0$ , the norm is $\tfrac12 J_1(\alpha_{j})^{2}$ , so

c_{j} = \frac{2}{J_1(\alpha_{j})^{2}}\cdot\frac{4J_1(\alpha_{j})}{\alpha_{j}^{3}} = \frac{8}{\alpha_{j}^{3}\,J_1(\alpha_{j})} .

Numerically:

$j$	$\alpha_{j}$	$J_1(\alpha_{j})$	$c_{j}$
1	2.4048256	0.5191475	1.1080223
2	5.5200781	$-0.3402648$	$-0.1397775$
3	8.6537279	0.2714523	0.0454765
4	11.7915344	$-0.2324598$	$-0.0209909$
5	14.9309177	0.2065464	0.0116362

Check 1: the value at the centre. $J_0(0) = 1$ , so $\sum_j c_j$ must equal $f(0) = 1$ . The partial sums are $1.0043666$ ( $5$ terms), $0.9998451$ ( $20$ terms), $0.9999898$ ( $60$ terms) — converging to $1$ , slowly, at the rate the coefficients dictate. That rate is $|c_j| = \Theta(j^{-5/2})$ , not $\Theta(j^{-3})$ : the explicit $\alpha_j^{-3}$ is partly undone by the $J_1(\alpha_j)$ in the denominator, which by the leading asymptotic form satisfies $|J_1(\alpha_j)| = \Theta(\alpha_j^{-1/2})$ , so $|c_j| = \Theta(\alpha_j^{-5/2}) = \Theta(j^{-5/2})$ since $\alpha_j \approx (j-\tfrac14)\pi$ . The products $|c_j|\alpha_j^{5/2}$ are $10.0221$ , $10.0260$ , $10.0265$ , $10.0265$ at $j = 4, 11, 51, 101$ , flat to five figures, whereas $|c_j|j^{3}$ is $1.343$ , $2.013$ , $4.144$ , $5.796$ at the same $j$ and still climbing.

Check 2: a value off centre. At $r = 0.5$ the exact value is $0.75$ . The partial sums are $0.7502500$ ( $5$ terms), $0.7499708$ ( $20$ terms), $0.7499989$ ( $60$ terms).

Check 3: Parseval, which gives a closed-form sum. Directly,

\int_0^{1}(1-r^{2})^{2}r\,dr = \int_0^{1}\bigl(r - 2r^{3} + r^{5}\bigr)dr = \frac12 - \frac12 + \frac16 = \frac16 .

Parseval's relation with $a = 1$ , $\nu = 0$ says this equals $\tfrac12\sum_j c_j^{2}J_1(\alpha_j)^{2}$ — here $f$ and the $c_j$ are real, so the moduli in the general statement are immaterial. Substituting $c_j = 8/(\alpha_j^{3}J_1(\alpha_j))$ , the factors $J_1(\alpha_j)^{2}$ cancel exactly and

\frac16 = \frac12\sum_{j=1}^{\infty}\frac{64}{\alpha_{j}^{6}} = 32\sum_{j=1}^{\infty}\alpha_{0j}^{-6} , \qquad\text{so}\qquad \boxed{\;\sum_{j=1}^{\infty}\frac{1}{\alpha_{0j}^{6}} = \frac{1}{192}\;}

a closed-form value for a sum over the zeros of $J_0$ , obtained purely from the normalisation constant. Summing the first $199$ terms numerically gives $0.0052083333$ against $1/192 = 0.0052083333$ — agreement to ten decimal places, and a decisive check that the normalisation $\tfrac12a^{2}J_{\nu+1}(\alpha_{\nu j})^{2}$ is right. (The same argument applied to $f(r) = 1$ , whose coefficients are $c_j = 2/[\alpha_j J_1(\alpha_j)]$ , gives $\sum_j \alpha_{0j}^{-2} = \tfrac14$ .)

ApplicationOriginal

The drumhead, and what a drum sounds like. A circular membrane of radius $a$ fixed at the rim has normal modes $J_m(\alpha_{mj}r/a)\,e^{\pm im\varphi}$ with angular frequencies $\omega_{mj} = c\,\alpha_{mj}/a$ , and the initial shape and velocity are matched by a Fourier-Bessel series in $j$ inside a Fourier series in $\varphi$ ((Strauss (2008), §10.2, pp. 267-269)). The orthogonality proved here is what makes the matching possible, and the normalisation is what makes the coefficients explicit. The frequencies are in the ratios $\alpha_{01}:\alpha_{11}:\alpha_{21}:\alpha_{02} = 2.4048:3.8317:5.1356:5.5201$ , that is $1:1.593:2.136:2.295$ — not integer multiples. That is exactly why a drum has no definite pitch while a string, whose eigenvalues are $n\pi$ , does; and it is why a timpanist loads the head to push the first few overtones towards harmonic ratios.

Heat and diffusion in a cylinder. For a cylinder of radius $a$ held at zero on its curved surface, the temperature relaxes as $T(r,t) = \sum_j c_j J_0(\alpha_{0j}r/a)\exp\bigl(-\kappa\alpha_{0j}^{2}t/a^{2}\bigr)$ , with $c_j$ the Fourier-Bessel coefficients of the initial profile. After a short time only the first term survives, so the cooling of any axially symmetric billet, wire or reactor fuel rod settles into a universal $J_0$ profile decaying with time constant $a^{2}/(\kappa\alpha_{01}^{2}) = a^{2}/(5.783\,\kappa)$ , independently of how it started. The identical series solves the corresponding Schrodinger problem for a particle in a cylindrical box, under the substitution $\kappa \to \mathrm{i}\hbar/2m$ that turns the diffusion equation $\partial_tT = \kappa\nabla^{2}T$ into $\mathrm{i}\hbar\,\partial_t\psi = -(\hbar^{2}/2m)\nabla^{2}\psi$ : the decaying factors become the phases $\exp(-\mathrm{i}E_jt/\hbar)$ with $E_j = \hbar^{2}\alpha_{0j}^{2}/(2ma^{2})$ , so the energies of the axially symmetric states are fixed by the same zeros. The sign matters, and so does the order of the substitution: writing $-\mathrm{i}\kappa$ in place of $\mathrm{i}\kappa$ runs the evolution backwards in time.

Electrostatics and the cylindrical Dirichlet problem. For a grounded cylindrical can with a prescribed potential on one end cap, separation gives $\psi = \sum_{m,j}c_{mj}J_m(\alpha_{mj}\rho/a)e^{im\varphi}\sinh(\alpha_{mj}z/a)$ , and the $c_{mj}$ come from a Fourier-Bessel integral in $\rho$ inside a Fourier integral in $\varphi$ ((Arfken, Weber & Harris (2013), §14.2, Example 14.2.1, pp. 663-665)). The same expansion with $J_m'$ zeros handles an insulated side wall.

Waveguides and optical fibres. The cutoff wavenumbers of the TM modes of a circular pipe are $\alpha_{mj}/a$ and of the TE modes $\beta_{mj}/a$ ; since $\beta_{11} = 1.8412 < \alpha_{01} = 2.4048$ , the lowest mode of a circular waveguide is TE $_{11}$ , and the single-mode band of the guide is the interval between those two numbers. The same zeros fix the single-mode condition of a step-index optical fibre.

From series to transform. Letting $a \to \infty$ at fixed $r$ turns the discrete spectrum $\alpha_{\nu j}/a$ into a continuum and the Fourier-Bessel series into the Hankel transform $F(\gamma) = \int_0^{\infty}f(r)J_\nu(\gamma r)r\,dr$ with the inversion $f(r) = \int_0^{\infty}F(\gamma)J_\nu(\gamma r)\gamma\,d\gamma$ , the orthogonality relation becoming $\int_0^{\infty}J_\nu(\gamma r)J_\nu(\gamma' r)r\,dr = \delta(\gamma-\gamma')/\gamma$ (Hankel and Mellin transforms). This is the standard tool for axisymmetric problems on an unbounded domain, and in two-dimensional imaging it is the radial part of the Fourier transform — which is why Fourier-Bessel expansions appear in computed tomography, in the point-spread functions of optical systems, and in the analysis of the cosmic microwave background on a disc of sky.

Definition

Modified Bessel functions

Adapted

Add to PDF

Rotating the argument of Bessel's equation through a right angle turns the oscillating solutions into a growing one and a decaying one: $I_\nu(x) = i^{-\nu}J_\nu(ix)$ is finite at the origin and grows like $e^{x}/\sqrt{2\pi x}$ , $K_\nu(x) = \tfrac{\pi}{2}i^{\nu+1}H^{(1)}_\nu(ix)$ is singular at the origin and decays like $\sqrt{\pi/2x}\,e^{-x}$ , and their Wronskian is exactly $-1/x$ .

Prerequisites: Neumann and Hankel functions, The Gamma function, Sturm comparison and oscillation theorems

Definition(The modified Bessel equation and the functions

I_\nu

K_\nu

)AdaptedAdapted from Arfken, Weber & Harris (2013), §14.5, Eq. (14.97), p. 680, Eqs. (14.98)-(14.99), p. 681, and Eqs. (14.106)-(14.107), pp. 682-683; Hassani (2013), §15.5, the equation after Eq. (15.46) and the two definitions that follow it, p. 484

The modified Bessel equation of order $\nu$ is

x^{2}y'' + x\,y' - (x^{2} + \nu^{2})\,y \;=\; 0 , \qquad x > 0 ,

which differs from Bessel's equation only in the sign of the $x^{2}$ term. Equivalently, in a problem separated in cylindrical coordinates the radial equation

\rho^{2}\frac{d^{2}P}{d\rho^{2}} + \rho\frac{dP}{d\rho} - (k^{2}\rho^{2} + \nu^{2})P \;=\; 0

becomes the above with $x = k\rho$ . Because the sign change is produced by replacing $k$ by $ik$ , every solution of the modified equation is a solution of Bessel's equation evaluated at a rotated argument: if $Z_\nu$ solves Bessel's equation then $Z_\nu(ix)$ solves the modified one.

The modified Bessel function of the first kind is the real-valued normalisation of $J_\nu(ix)$ ,

I_\nu(x) \;=\; i^{-\nu}J_\nu(ix) \;=\; e^{-i\nu\pi/2}J_\nu\bigl(xe^{i\pi/2}\bigr) \;=\; \sum_{s=0}^{\infty}\frac{1}{s!\,\Gamma(s+\nu+1)}\Bigl(\frac{x}{2}\Bigr)^{\nu+2s} .

The second solution is written $K_\nu$ and called the modified Bessel function of the second kind. Rather than being pinned to a series, it is pinned to the Hankel function $H^{(1)}_\nu$ :

K_\nu(x) \;=\; \frac{\pi}{2}\,i^{\nu+1}H^{(1)}_\nu(ix) \;=\; \frac{\pi}{2}\,i^{\nu+1}\bigl[J_\nu(ix) + iY_\nu(ix)\bigr] ,

the factor $i^{\nu+1}$ being chosen to make $K_\nu(x)$ real for real $x > 0$ . An equivalent form, derived below, is

K_\nu(x) \;=\; \frac{\pi}{2}\,\frac{I_{-\nu}(x) - I_\nu(x)}{\sin\nu\pi} , \qquad \nu \notin \mathbb{Z} , \qquad\qquad K_n(x) \;=\; \lim_{\nu\to n}K_\nu(x) ,

the limit being needed because $I_n = I_{-n}$ makes the quotient indeterminate at integer order. The pair $\{I_\nu, K_\nu\}$ is the standard basis.

Two names are in circulation. Because $I_\nu$ stands to $J_\nu$ as $\sinh$ stands to $\sin$ , the two functions are sometimes called hyperbolic Bessel functions; $K_\nu$ is sometimes called a Whittaker function (Arfken, Weber & Harris (2013), §14.5, pp. 682-683). We use $I_\nu$ and $K_\nu$ throughout. For the convention followed here, the branch of $K_n$ taken for real $x>0$ is the one that is real there.

Definition(Basic properties)AdaptedAdapted from Arfken, Weber & Harris (2013), §14.5, Eqs. (14.100)-(14.101), p. 681, Eqs. (14.104)-(14.105), pp. 681-682, Eqs. (14.108)-(14.111), p. 683, and the Wronskian value on p. 685 (Exercise 14.5.11, p. 687); Arfken, Weber & Harris (2013), the large-$x$ forms (14.127), p. 690, and (14.143), p. 693; Hassani (2013), §14.4, Corollary 14.4.12 and Example 14.4.13, p. 432 (item 6)

Let $x > 0$ .

Symmetry in the order. For every integer $n$ , $I_n(x) = I_{-n}(x)$ ; for every $\nu$ , $K_\nu(x) = K_{-\nu}(x)$ . For non-integer $\nu$ the functions $I_\nu$ and $I_{-\nu}$ are linearly independent, and $\{I_\nu, I_{-\nu}\}$ is then also a basis; at integer order it collapses, which is the reason for introducing $K_\nu$ .
Recurrences. For all $\nu$ ,
$I_{\nu-1}(x) - I_{\nu+1}(x) = \frac{2\nu}{x}\,I_\nu(x) , \qquad I_{\nu-1}(x) + I_{\nu+1}(x) = 2\,I_\nu'(x) ,$
while the $K$ family satisfies the same two relations with both right-hand sides negated:
$K_{\nu-1}(x) - K_{\nu+1}(x) = -\frac{2\nu}{x}\,K_\nu(x) , \qquad K_{\nu-1}(x) + K_{\nu+1}(x) = -2\,K_\nu'(x) .$
In particular $I_0' = I_1$ and $K_0' = -K_1$ .
Wronskian. For every $\nu$ ,
$I_\nu(x)K_\nu'(x) - I_\nu'(x)K_\nu(x) \;=\; -\frac{1}{x} ,$
equivalently $I_\nu(x)K_{\nu+1}(x) + I_{\nu+1}(x)K_\nu(x) = 1/x$ . The constant in Abel's identity — that is, $x$ times the Wronskian — is exactly $-1$ , with no $\nu$ dependence; so $I_\nu$ and $K_\nu$ are independent for every order, integer or not.
Behaviour at the origin. As $x \to 0^{+}$ ,
$I_\nu(x) = \frac{x^{\nu}}{2^{\nu}\Gamma(\nu+1)} + \cdots , \qquad K_0(x) = -\ln x - \gamma + \ln 2 + \cdots , \qquad K_\nu(x) = 2^{\nu-1}\Gamma(\nu)\,x^{-\nu} + \cdots \ \ (\nu > 0) ,$
with $\gamma$ the Euler-Mascheroni constant. So for $\nu \ge 0$ the function $I_\nu$ is finite at $x = 0$ (and vanishes there unless $\nu = 0$ , where $I_0(0) = 1$ ), while $K_\nu$ is singular: logarithmically at $\nu = 0$ and like $x^{-\nu}$ otherwise.
Behaviour at infinity. As $x \to +\infty$ ,
$I_\nu(x) \sim \frac{e^{x}}{\sqrt{2\pi x}} , \qquad K_\nu(x) \sim \sqrt{\frac{\pi}{2x}}\;e^{-x} ,$
for every fixed $\nu$ : one solution grows exponentially, the other decays exponentially, and the leading factor does not depend on the order. The full asymptotic expansions, and the derivations of these leading terms, are in Asymptotic forms of Bessel functions.
Non-oscillation. For $\nu \ge \tfrac12$ no non-trivial solution of the modified Bessel equation has more than one zero in $(0,\infty)$ . This is the precise content of the informal statement that the modified Bessel functions are not oscillatory, in contrast to $J_\nu$ and $Y_\nu$ , which have infinitely many zeros.

Idea of the proofOriginal

Everything here follows from one move: turn the argument of Bessel's equation through a right angle. Bessel's equation says, in its Liouville normal form, that a solution multiplied by $\sqrt{x}$ obeys $v'' + (1 - \tfrac{4\nu^{2}-1}{4x^{2}})v = 0$ — that is, $v = \sqrt{x}\,y$ , the substitution written as $y = v/\sqrt{x}$ in the remark below, which at large $x$ is the harmonic oscillator; replacing $x$ by $ix$ flips the sign of that $1$ and gives $v'' - (1 + \tfrac{4\nu^{2}-1}{4x^{2}})v = 0$ , which is $v'' = v$ at large $x$ . An oscillator becomes a pair of real exponentials. All the oscillatory structure of $J_\nu$ and $Y_\nu$ — infinitely many zeros, orthogonality on an interval, a $1/\sqrt{x}$ -damped sinusoid at infinity — disappears, and is replaced by a two-dimensional space of solutions in which one direction grows like $e^{x}$ and the other decays like $e^{-x}$ .

What is left to do is bookkeeping, and it is where all the conventions live. The rotated function $J_\nu(ix)$ has a series whose terms all carry the same phase $i^{\nu}$ , so dividing that phase out gives a real function with positive coefficients: that is $I_\nu$ . For the second solution, taking $Y_\nu(ix)$ would be the obvious choice, but it is a bad one, because $Y_\nu(ix)$ contains both exponentials. The decaying solution is the one a physicist wants, and the function that decays at infinity is the rotated Hankel function $H^{(1)}_\nu$ , whose asymptotic form is a single outgoing exponential $e^{ix}$ : rotating that argument makes $e^{-x}$ . So $K_\nu$ is defined from $H^{(1)}_\nu(ix)$ , with a phase factor put in front to make it real. The price of defining $K_\nu$ by its behaviour at infinity rather than by a series is that it does not inherit the recurrences of $I_\nu$ unchanged; the two signs flip.

The Wronskian is then the only quantitative input needed to use the pair, and Abel's identity gives it for free up to a constant: the equation in self-adjoint form is $(xy')' - (x + \nu^{2}/x)y = 0$ , so the Wronskian times $x$ is constant. The constant is fixed by multiplying the leading behaviours at the origin, where the Gamma-function identity $\nu\Gamma(\nu) = \Gamma(\nu+1)$ makes everything cancel and leaves exactly $-1$ .

Proof(The series for

I_\nu

, and the quotient form of

K_\nu

)AdaptedAdapted from Arfken, Weber & Harris (2013), §14.5, Eqs. (14.98)-(14.99), p. 681, and Eqs. (14.106)-(14.107), pp. 682-683

Step 1 (the rotated series). Putting $ix$ for $x$ in the series for $J_\nu$ gives

J_\nu(ix) \;=\; \sum_{s=0}^{\infty}\frac{(-1)^{s}}{s!\,\Gamma(s+\nu+1)}\Bigl(\frac{ix}{2}\Bigr)^{\nu+2s} \;=\; i^{\nu}\sum_{s=0}^{\infty}\frac{(-1)^{s}i^{2s}}{s!\,\Gamma(s+\nu+1)}\Bigl(\frac{x}{2}\Bigr)^{\nu+2s} ,

and $i^{2s} = (-1)^{s}$ cancels the alternating sign term by term. Hence $i^{-\nu}J_\nu(ix)$ has the stated series, all of whose terms are positive for $x>0$ and $\nu > -1$ ; in particular it is real. This is the definition of $I_\nu$ , and it converges for every $x$ because the $J_\nu$ series does.

Step 2 (symmetry at integer order). From $J_\nu(ix) = i^{\nu}I_\nu(x)$ and the integer-order relation $J_{-n} = (-1)^{n}J_n$ ,

I_{-n}(x) \;=\; i^{\,n}J_{-n}(ix) \;=\; i^{\,n}(-1)^{n}J_n(ix) \;=\; i^{\,2n}(-1)^{n}I_n(x) \;=\; I_n(x) ,

since $i^{2n}(-1)^{n} = (-1)^{n}(-1)^{n} = 1$ . So the two solutions $I_\nu, I_{-\nu}$ coincide at integer order, exactly as $J_n$ and $J_{-n}$ do.

Step 3 (from the Hankel definition to the quotient). Write $H^{(1)}_\nu = J_\nu + iY_\nu$ with $Y_\nu = (\cos\nu\pi\,J_\nu - J_{-\nu})/\sin\nu\pi$ , and evaluate at $ix$ with $\nu \notin\mathbb{Z}$ :

J_\nu(ix) + iY_\nu(ix) \;=\; J_\nu(ix)\,\frac{\sin\nu\pi + i\cos\nu\pi}{\sin\nu\pi} \;-\; \frac{i}{\sin\nu\pi}\,J_{-\nu}(ix) .

Now use three identities: $\sin\nu\pi + i\cos\nu\pi = i(\cos\nu\pi - i\sin\nu\pi) = ie^{-i\nu\pi}$ ; $J_\nu(ix) = i^{\nu}I_\nu(x)$ ; and $J_{-\nu}(ix) = i^{-\nu}I_{-\nu}(x)$ . The first factor becomes $i^{\nu}e^{-i\nu\pi} = e^{i\nu\pi/2}e^{-i\nu\pi} = e^{-i\nu\pi/2} = i^{-\nu}$ , so

H^{(1)}_\nu(ix) \;=\; \frac{i\,i^{-\nu}}{\sin\nu\pi}\bigl[I_\nu(x) - I_{-\nu}(x)\bigr] .

Multiplying by $\tfrac{\pi}{2}i^{\nu+1}$ and using $i\cdot i = -1$ gives

K_\nu(x) \;=\; \frac{\pi}{2}\,\frac{I_{-\nu}(x) - I_\nu(x)}{\sin\nu\pi} ,

which is the quotient form. Replacing $\nu$ by $-\nu$ changes the sign of both the numerator and $\sin\nu\pi$ , so $K_{-\nu} = K_\nu$ for non-integer $\nu$ , and by continuity of the limit also at integer order. At integer order the quotient is $0/0$ , and the limit $\nu\to n$ exists by l'Hôpital's rule applied in $\nu$ ; carrying the differentiation out produces a term $(-1)^{n+1}I_n(x)\ln(x/2)$ alongside power series in $x$ , which is the source of the logarithm in $K_0$ (Hassani (2013), §15.5, the series for $K_n$, p. 484). $\square$

Proof(The recurrences)AdaptedAdapted from Arfken, Weber & Harris (2013), §14.5, Eqs. (14.102)-(14.105), pp. 681-682, and Eqs. (14.108)-(14.109), p. 683

Step 1 (the $I$ family). Start from the two Bessel recurrences $J_{\nu-1}(z) + J_{\nu+1}(z) = (2\nu/z)J_\nu(z)$ and $J_{\nu-1}(z) - J_{\nu+1}(z) = 2J_\nu'(z)$ , and set $z = ix$ . Differentiating $J_\nu(ix) = i^{\nu}I_\nu(x)$ in $x$ gives $i\,J_\nu'(ix) = i^{\nu}I_\nu'(x)$ , that is $J_\nu'(ix) = i^{\nu-1}I_\nu'(x)$ .

In the first recurrence, $i^{\nu-1}I_{\nu-1} + i^{\nu+1}I_{\nu+1} = \frac{2\nu}{ix}i^{\nu}I_\nu$ . Dividing by $i^{\nu-1}$ and using $i^{2} = -1$ ,

I_{\nu-1}(x) - I_{\nu+1}(x) \;=\; \frac{2\nu}{x}\,I_\nu(x) .

In the second, $i^{\nu-1}I_{\nu-1} - i^{\nu+1}I_{\nu+1} = 2i^{\nu-1}I_\nu'$ ; dividing by $i^{\nu-1}$ gives $I_{\nu-1} + I_{\nu+1} = 2I_\nu'$ .

So the subtraction and addition have changed places relative to the $J$ family. That swap is bookkeeping, not an argument about zeros: neither relation on its own forbids $I_\nu$ from vanishing. What rules out zeros of $I_\nu$ on $(0,\infty)$ for $\nu\ge0$ is the series of Step 1, every term of which is positive.

Step 2 (the $K$ family). Work from the quotient form, for non-integer $\nu$ ; the integer case then follows by taking $\nu\to n$ in the finished identities. Since $\sin(\nu\pm1)\pi = -\sin\nu\pi$ ,

K_{\nu-1} = \frac{\pi}{2\sin\nu\pi}\bigl(I_{\nu-1} - I_{1-\nu}\bigr) , \qquad K_{\nu+1} = \frac{\pi}{2\sin\nu\pi}\bigl(I_{\nu+1} - I_{-\nu-1}\bigr) .

Subtracting,

K_{\nu-1} - K_{\nu+1} = \frac{\pi}{2\sin\nu\pi}\Bigl[\bigl(I_{\nu-1} - I_{\nu+1}\bigr) - \bigl(I_{1-\nu} - I_{-\nu-1}\bigr)\Bigr] .

The first bracket is $\frac{2\nu}{x}I_\nu$ by Step 1. For the second, apply Step 1 at order $-\nu$ : $I_{-\nu-1} - I_{-\nu+1} = \frac{-2\nu}{x}I_{-\nu}$ , hence $I_{1-\nu} - I_{-\nu-1} = \frac{2\nu}{x}I_{-\nu}$ . Therefore

K_{\nu-1} - K_{\nu+1} = \frac{2\nu}{x}\cdot\frac{\pi}{2\sin\nu\pi}\bigl(I_\nu - I_{-\nu}\bigr) = -\frac{2\nu}{x}\,K_\nu .

Adding instead of subtracting, $K_{\nu-1} + K_{\nu+1} = \frac{\pi}{2\sin\nu\pi}\bigl[(I_{\nu-1}+I_{\nu+1}) - (I_{1-\nu}+I_{-\nu-1})\bigr] = \frac{\pi}{2\sin\nu\pi}\bigl[2I_\nu' - 2I_{-\nu}'\bigr] = -2K_\nu'$ , using $I_{\nu-1}+I_{\nu+1} = 2I_\nu'$ at orders $\nu$ and $-\nu$ . $\square$

Proof(The Wronskian is exactly

-1/x

)Original

Step 1 (the $x$ -dependence, from Abel's identity). Multiplying the modified Bessel equation by $1/x$ puts it in self-adjoint form,

\bigl(x\,y'\bigr)' - \Bigl(x + \frac{\nu^{2}}{x}\Bigr)y \;=\; 0 ,

so the coefficient of $y''$ in self-adjoint form is $p(x) = x$ , with $q(x) = -(x + \nu^{2}/x)$ . For any two solutions $u, v$ of $(py')' + qy = 0$ ,

\frac{d}{dx}\Bigl[p\bigl(uv' - u'v\bigr)\Bigr] = u\,(pv')' - v\,(pu')' = u\,(-qv) - v\,(-qu) = 0 ,

the first equality because the terms $pu'v'$ cancel. Hence

x\bigl[I_\nu(x)K_\nu'(x) - I_\nu'(x)K_\nu(x)\bigr] \;=\; C_\nu

is constant on $(0,\infty)$ .

Step 2 (the constant, for $\nu > 0$ ). Evaluate the constant from the behaviour at the origin, item 4 of the statement: $I_\nu(x) = a\,x^{\nu}(1+o(1))$ with $a = 1/\bigl(2^{\nu}\Gamma(\nu+1)\bigr)$ , and $K_\nu(x) = b\,x^{-\nu}(1+o(1))$ with $b = 2^{\nu-1}\Gamma(\nu)$ . Differentiating the leading terms,

I_\nu K_\nu' - I_\nu' K_\nu \;=\; a x^{\nu}\bigl(-\nu b x^{-\nu-1}\bigr) - \nu a x^{\nu-1}\bigl(b x^{-\nu}\bigr) + o(x^{-1}) \;=\; -\frac{2\nu ab}{x} + o(x^{-1}) .

The product is

2\nu ab \;=\; 2\nu\cdot\frac{2^{\nu-1}\Gamma(\nu)}{2^{\nu}\Gamma(\nu+1)} \;=\; \frac{\nu\,\Gamma(\nu)}{\Gamma(\nu+1)} \;=\; 1 ,

by the functional equation $\Gamma(\nu+1) = \nu\Gamma(\nu)$ of the Gamma function. So $C_\nu = -1$ for every $\nu > 0$ .

Step 3 (the remaining orders). For $\nu = 0$ the same computation runs with $I_0 \to 1$ and $K_0 = -\ln x + O(1)$ , $K_0' = -1/x + o(1/x)$ : then $I_0K_0' - I_0'K_0 = -1/x + o(1/x)$ , so $C_0 = -1$ as well, the $I_0'K_0$ term being $O(x\ln x)$ and negligible. For $\nu < 0$ use $K_{-\nu} = K_\nu$ and $I_{-\nu} = I_\nu + \tfrac{2}{\pi}\sin(\nu\pi)K_\nu$ — which is the quotient form of $K_\nu$ rearranged — so that $I_{-\nu}K_{-\nu}' - I_{-\nu}'K_{-\nu} = I_\nu K_\nu' - I_\nu'K_\nu$ , the added multiple of $K_\nu$ contributing $\tfrac{2}{\pi}\sin(\nu\pi)(K_\nu K_\nu' - K_\nu'K_\nu) = 0$ . Hence $C_\nu = -1$ for all $\nu$ .

The second form in item 3 follows from the first with the recurrences: subtracting the two relations of each family in item 2 gives $K_\nu' = -K_{\nu+1} + \tfrac{\nu}{x}K_\nu$ and $I_\nu' = I_{\nu+1} + \tfrac{\nu}{x}I_\nu$ . Substituting both,

I_\nu K_\nu' - I_\nu' K_\nu = -I_\nu K_{\nu+1} + \frac{\nu}{x}I_\nu K_\nu - I_{\nu+1}K_\nu - \frac{\nu}{x}I_\nu K_\nu = -\bigl(I_\nu K_{\nu+1} + I_{\nu+1}K_\nu\bigr) ,

the two $\tfrac{\nu}{x}I_\nu K_\nu$ terms cancelling against each other with nothing left over; and the left-hand side is $-1/x$ by Step 2. As a numerical check at $x = 1$ : $I_0(1)K_1(1) + I_1(1)K_0(1) = (1.2660659)(0.6019072) + (0.5651591)(0.4210244) = 0.7620542 + 0.2379458 = 1.0000000$ . $\square$

Remark(Where the hypotheses are used, and what the conventions cost)AdaptedAdapted from Arfken, Weber & Harris (2013), §14.5, the footnotes on p. 683 and the alternative $\\cos\\nu\\pi$ convention, p. 683; Hassani (2013), §14.4, Corollary 14.4.12 and Example 14.4.13, p. 432

$x > 0$ throughout. Every statement above is on the open half-line. At $x = 0$ the equation is singular (the coefficient of $y''$ vanishes), $K_\nu$ is not defined there, and the Wronskian identity $-1/x$ blows up — correctly, since two solutions cannot both stay bounded with bounded derivative at a point where the equation degenerates.

Where $\nu\notin\mathbb{Z}$ is needed. Only in the quotient form of $K_\nu$ , whose denominator $\sin\nu\pi$ vanishes at integer order. Every conclusion about $K_\nu$ — the recurrences, the Wronskian, the asymptotics — holds at integer order too, but each is obtained as a limit $\nu\to n$ of an identity proved off the integers, or else (as in the Wronskian proof) directly. What does not survive the limit is the pair $\{I_\nu, I_{-\nu}\}$ as a basis: by item 1 it degenerates, and this is the whole reason $K_\nu$ is introduced.

Where $\nu \ge \tfrac12$ is needed in item 6. Setting $y = v/\sqrt{x}$ removes the first-derivative term and turns the modified equation into its Liouville normal form

v'' - \Bigl(1 + \frac{4\nu^{2}-1}{4x^{2}}\Bigr)v \;=\; 0 ,

the same calculation that turns Bessel's equation into $v'' + (1 - \tfrac{4\nu^{2}-1}{4x^{2}})v = 0$ Hassani (Hassani (2013), Example 14.4.13, p. 432). The coefficient of $v$ is $\le 0$ precisely when $4x^{2} + 4\nu^{2} - 1 \ge 0$ , which for $\nu \ge \tfrac12$ holds at every $x>0$ , and Sturm's comparison theorem then forbids more than one zero. For $\nu < \tfrac12$ the hypothesis fails only on the initial segment $0 < x < \tfrac12\sqrt{1-4\nu^{2}}$ , so the conclusion still holds unchanged on $x \ge \tfrac12\sqrt{1-4\nu^{2}}$ . On the initial segment it survives as well, but a bounded interval is not on its own a reason: a solution can have infinitely many zeros in one. What rules them out is a second comparison. Writing the normal form as $v'' + qv = 0$ with $q(x) = -1 + (1-4\nu^{2})/4x^{2} \le 1/4x^{2}$ , compare with the Euler equation $w'' + w/4x^{2} = 0$ , whose solutions are $\sqrt{x}$ and $\sqrt{x}\ln x$ . If $v$ had two zeros in $(0,\tfrac12\sqrt{1-4\nu^{2}})$ then, since $q \le 1/4x^{2}$ , Sturm's theorem would force every solution of the Euler equation to vanish strictly between them — and $w = \sqrt{x}$ never vanishes on $(0,\infty)$ . So there is at most one zero there too.

The price of defining $K_\nu$ by its behaviour at infinity. Because $K_\nu$ is pinned to $H^{(1)}_\nu$ rather than to a power series, it fails to satisfy the recurrences of $I_\nu$ : both right-hand signs flip. Some authors therefore insert an extra factor $\cos\nu\pi$ into the definition, which restores the $I$ recurrences but makes the function vanish identically at $\nu = \tfrac12, \tfrac32, \tfrac52, \dots$ (Arfken, Weber & Harris (2013), §14.5, p. 683). We follow Arfken's convention, Eq. (14.106), with no such factor; with it, the identities of item 2 are the ones to use. A second consequence of the same choice: for non-integer $\nu$ the function $K_\nu$ has a branch point at the origin from the fractional power, and for integer $n$ from the logarithm, so " $K_n$ " always means the branch that is real on the positive axis.

What the sign change in the equation destroys. Bessel's equation on a finite interval is a regular Sturm-Liouville problem with an infinite sequence of eigenvalues, which is what produces the zeros of $J_\nu$ and the Fourier-Bessel expansions of Orthogonality and Fourier-Bessel series. None of that transfers. The modified equation has no such spectrum on $(0,\infty)$ : by item 6 its solutions do not oscillate, there are no interior zeros to impose a boundary condition at, and there is no analogue of a Fourier-Bessel series in $I_\nu$ .

Method(Choosing the right modified Bessel function in a boundary-value problem)Original

A separated radial problem produces the modified Bessel equation whenever the separation constant enters with the sign opposite to the oscillatory case — for instance a Laplace or Helmholtz problem in cylindrical coordinates whose $z$ -dependence is $\cos kz$ rather than $e^{-kz}$ , or a screened (Yukawa) operator $-\nabla^{2} + m^{2}$ . Given

\frac{1}{\rho}\frac{d}{d\rho}\Bigl(\rho\frac{du}{d\rho}\Bigr) - \Bigl(m^{2} + \frac{\nu^{2}}{\rho^{2}}\Bigr)u \;=\; 0 ,

proceed as follows.

Put it in standard form. Substitute $x = m\rho$ . The general solution is $u = A\,I_\nu(x) + B\,K_\nu(x)$ .
Impose the condition at the inner end. If the domain includes the axis $\rho = 0$ and $u$ must be bounded there, set $B = 0$ by item 4: $K_\nu$ is the singular solution. If instead the axis is excluded (an annulus, a coaxial geometry), keep both constants.
Impose the condition at the outer end. If the domain is unbounded and $u$ must decay, set $A = 0$ by item 5: $I_\nu$ is the growing solution. If both conditions apply — bounded at the axis and decaying at infinity — there is no non-trivial solution, which is the statement that the homogeneous problem has no eigenvalues (see the last paragraph of the remark).
Fix the surviving constants from the data. With two boundary values $u(\rho_1) = u_1$ , $u(\rho_2) = u_2$ on an annulus, solve the $2\times2$ system; the determinant is $I_\nu(x_1)K_\nu(x_2) - I_\nu(x_2)K_\nu(x_1)$ , which is non-zero for $x_1\neq x_2$ because $I_\nu$ is increasing and $K_\nu$ decreasing.
Differentiate with the recurrences, not by hand. Fluxes and normal derivatives are wanted as often as values: use $I_0' = I_1$ , $K_0' = -K_1$ , and in general $I_\nu' = I_{\nu+1} + \tfrac{\nu}{x}I_\nu$ , $K_\nu' = -K_{\nu+1} + \tfrac{\nu}{x}K_\nu$ — note the $+$ in the second, which is what item 2 gives and is easy to get wrong.
Check against the Wronskian. $I_\nu(x)K_{\nu+1}(x) + I_{\nu+1}(x)K_\nu(x) = 1/x$ at any $x$ catches a wrong function value or a mixed-up order immediately.

Example(Steady screened potential in an annulus)Original

Solve

\frac{1}{r}\frac{d}{dr}\Bigl(r\frac{du}{dr}\Bigr) - u \;=\; 0 \quad\text{on } 1 \le r \le 2 , \qquad u(1) = 1 , \qquad u(2) = 0 ,

and report $u(1.5)$ and the flux $-u'(2)$ at the outer wall.

Step 1 (identify the equation). Expanding, $r^{2}u'' + ru' - r^{2}u = 0$ : the modified Bessel equation with $\nu = 0$ and $x = r$ (so $m = 1$ ). Hence

u(r) \;=\; A\,I_0(r) + B\,K_0(r) .

The axis is not in the domain, so both constants survive — step 2 of the method.

Step 2 (apply the boundary conditions). The two conditions read

A\,I_0(1) + B\,K_0(1) = 1 , \qquad A\,I_0(2) + B\,K_0(2) = 0 .

With $I_0(1) = 1.2660659$ , $K_0(1) = 0.4210244$ , $I_0(2) = 2.2795853$ , $K_0(2) = 0.1138939$ , the determinant is

D = I_0(1)K_0(2) - I_0(2)K_0(1) = 0.1441971 - 0.9597611 = -0.8155640 ,

non-zero as step 4 of the method promised. Cramer's rule gives

A = \frac{K_0(2)}{D} = \frac{0.1138939}{-0.8155640} = -0.1396504 , \qquad B = \frac{-I_0(2)}{D} = \frac{-2.2795853}{-0.8155640} = 2.7951030 .

Step 3 (evaluate). With $I_0(1.5) = 1.6467232$ and $K_0(1.5) = 0.2138056$ ,

u(1.5) = (-0.1396504)(1.6467232) + (2.7951030)(0.2138056) = -0.2299656 + 0.5976086 = 0.3676429 .

Step 4 (the flux). By step 5 of the method, $u' = A I_0' + B K_0' = A\,I_1(r) - B\,K_1(r)$ . At $r = 2$ , with $I_1(2) = 1.5906369$ and $K_1(2) = 0.1398659$ ,

u'(2) = (-0.1396504)(1.5906369) - (2.7951030)(0.1398659) = -0.2221331 - 0.3909395 = -0.6130727 ,

so the outward flux is $-u'(2) = 0.6130727$ .

Step 5 (check). Two independent checks. First, the solution must be positive and decreasing on the interval, since $u$ satisfies $(ru')' = ru > 0$ wherever $u>0$ : and indeed $u(1) = 1 > u(1.5) = 0.3676 > u(2) = 0$ with $u'(2) < 0$ . Second, the Wronskian check of step 6 at $x = 1.5$ : $I_0(1.5)K_1(1.5) + I_1(1.5)K_0(1.5) = (1.6467232)(0.2773878) + (0.9816664)(0.2138056) = 0.4567809 + 0.2098857 = 0.6666667 = 1/1.5$ , as required. All function values are rounded from a 15-digit computation; the quoted digits of $u(1.5) = 0.367643$ and $-u'(2) = 0.613073$ are correct.

Application(Where

I_\nu

and

K_\nu

turn up)AdaptedAdapted from Arfken, Weber & Harris (2013), §14.5, Example 14.5.1 and the Wronskian factor, pp. 684-685, and the summary, p. 685; Hassani (2013), §19.6.3, the skin-depth equation (19.35) and the Kelvin equation (19.36), p. 589, and the Kelvin functions ber, bei, p. 590

The cylindrical Green's function. The fundamental solution of the Laplacian in cylindrical coordinates, expanded in the azimuthal angle and Fourier-transformed in $z$ , has radial factor $I_m(k\rho_<)K_m(k\rho_>)$ , where $\rho_<$ and $\rho_>$ are the smaller and larger of the two radii. The reason is exactly items 4 and 5 above: the factor must be bounded on the axis, so $I_m$ is used for the inner point, and must decay at large radius, so $K_m$ is used for the outer. The normalising constant is set by the jump condition, which is the Wronskian $-1/(k\rho)$ of item 3 (Arfken, Weber & Harris (2013), Example 14.5.1, pp. 684-685). This is a worked case of the construction in Green's functions for Laplace and Poisson problems.

The two-dimensional screened Coulomb (Yukawa) potential. For $(-\nabla^{2} + m^{2})G = \delta^{(2)}$ in the plane, $G$ must solve the modified Bessel equation of order zero away from the origin, must decay at infinity, and must have the planar logarithmic singularity $-\tfrac{1}{2\pi}\ln r$ at the origin. By items 4 and 5 the only candidate is a multiple of $K_0(mr)$ , and matching $K_0(x) = -\ln x - \gamma + \ln 2 + \cdots$ against $-\tfrac{1}{2\pi}\ln r$ fixes it: $G = \tfrac{1}{2\pi}K_0(mr)$ . The exponential decay of $K_0$ is the screening length $1/m$ . This is the planar analogue of the three-dimensional $e^{-mr}/4\pi r$ ; the general construction is in Green's functions for the Helmholtz equation.

Heat loss from a fin or a buried pipe. A thin radial fin losing heat to its surroundings at a rate proportional to the local excess temperature satisfies, in the steady state, the modified Bessel equation of order zero with $m^{2}$ the loss coefficient over the conductivity. Bounded temperature on the axis selects $I_0$ ; a prescribed temperature at the rim normalises it, as in the worked example above. The engineering quantity of interest, the total heat flow, is the flux computed with $I_0' = I_1$ .

The skin effect. The current density in a long cylindrical wire carrying alternating current satisfies a Bessel-type equation with purely imaginary parameter, $w'' + w'/x - ik^{2}w = 0$ , whose solutions are Bessel functions of complex argument — the Kelvin functions $\operatorname{ber}$ and $\operatorname{bei}$ , obtained from $J_0$ by rotating the argument through $45^{\circ}$ rather than $90^{\circ}$ Hassani (Hassani (2013), §19.6.3, Eqs. (19.35)-(19.36) and the definition of ber and bei, pp. 589-590). The skin depth is the decay length of the resulting exponential, and the mechanism is the one in the idea above: a sign (here a phase) in front of the $x^{2}$ term decides whether the radial solutions oscillate or decay.

Definition

Spherical Bessel functions

Adapted

Add to PDF

The radial part of the Helmholtz equation in spherical coordinates is Bessel's equation of half-integer order; normalised as $j_n(x) = \sqrt{\pi/2x}\,J_{n+1/2}(x)$ the solutions are elementary — $j_0 = \sin x/x$ and all the rest obtained from it by the Rayleigh formula $f_n = (-1)^nx^n(x^{-1}d/dx)^nf_0$ — and they expand a plane wave as $e^{ikr\cos\theta} = \sum_l(2l+1)i^lj_l(kr)P_l(\cos\theta)$ .

Prerequisites: Neumann and Hankel functions, Legendre polynomials

Definition(The spherical Bessel equation and its four solutions)AdaptedAdapted from Arfken, Weber & Harris (2013), §14.7, Eqs. (14.148)-(14.151), pp. 698-699; Hassani (2013), §19.7.2, Eqs. (19.39)-(19.40), p. 593; Riley, Hobson & Bence (2006), §18.6, Eqs. (18.101)-(18.103), pp. 614-615

Separating the Helmholtz equation $\nabla^{2}\psi + k^{2}\psi = 0$ in spherical coordinates leaves, for the radial factor, the spherical Bessel equation

r^{2}\frac{d^{2}R}{dr^{2}} + 2r\frac{dR}{dr} + \bigl[k^{2}r^{2} - n(n+1)\bigr]R \;=\; 0 ,

in which $n(n+1)$ is the separation constant of the angular equation and $n$ is a non-negative integer whenever the angular factor is required to be single-valued on the whole sphere. Substituting

R(kr) \;=\; \frac{Z(kr)}{(kr)^{1/2}}

turns it into Bessel's equation of order $n + \tfrac12$ (the computation is in the proof below), so the solutions are $(kr)^{-1/2}$ times Bessel functions of half-integer order. Because the equation is homogeneous, the overall scale is a matter of convention; the standard choice fixes the factor $\sqrt{\pi/2x}$ . For integer $n$ (not necessarily non-negative) define

j_n(x) = \sqrt{\frac{\pi}{2x}}\,J_{n+1/2}(x) , \qquad y_n(x) = \sqrt{\frac{\pi}{2x}}\,Y_{n+1/2}(x) ,

h^{(1)}_n(x) = \sqrt{\frac{\pi}{2x}}\,H^{(1)}_{n+1/2}(x) = j_n(x) + i\,y_n(x) , \qquad h^{(2)}_n(x) = \sqrt{\frac{\pi}{2x}}\,H^{(2)}_{n+1/2}(x) = j_n(x) - i\,y_n(x) .

These are the spherical Bessel function $j_n$ , the spherical Neumann function $y_n$ , and the spherical Hankel functions $h^{(1)}_n$ , $h^{(2)}_n$ . All four are covered at once by writing $z_n(x) = \sqrt{\pi/2}\;Z_{n+1/2}(x)/\sqrt{x}$ with $Z$ any solution of Bessel's equation of order $n+\tfrac12$ Hassani (Hassani (2013), Eq. (19.40), p. 593). The general solution of the radial equation is $R(r) = A\,j_n(kr) + B\,y_n(kr)$ .

Throughout, $f_n$ stands for any one of $j_n$ , $y_n$ , $h^{(1)}_n$ , $h^{(2)}_n$ , since all four satisfy the same recurrences. Double factorials are used in the form

(2n+1)!! = 1\cdot3\cdot5\cdots(2n+1) , \qquad (2n-1)!! = 1\cdot3\cdot5\cdots(2n-1) , \qquad (-1)!! = 1 .

Notation. We write $y_n$ for the spherical Neumann function, with Arfken and Hassani. Riley, Hobson and Bence write $n_\ell$ for the same function (Riley, Hobson & Bence (2006), §18.6, Eq. (18.103), p. 615), and in quantum-scattering literature $n_\ell$ is at least as common as $y_\ell$ ; the functions are identical, only the letter differs. The index is written $n$ here and $\ell$ or $l$ in texts that reserve $\ell$ for angular momentum.

Definition(Closed forms, recurrences and the two ends)AdaptedAdapted from Arfken, Weber & Harris (2013), §14.7, Eqs. (14.152)-(14.162), pp. 699-702, Eqs. (14.169)-(14.182), pp. 702-703; Riley, Hobson & Bence (2006), §18.6, Eqs. (18.104)-(18.106) and the low-order list, pp. 615-616; Strauss (2008), §10.5, Eq. (12), p. 285

Let $n$ be an integer and $x > 0$ .

Reflection in the index. $y_n(x) = (-1)^{n+1}j_{-n-1}(x)$ . So the second kind is not an independent object: negative indices of $j$ already supply it.
Series. With $(a)_s = a(a+1)\cdots(a+s-1)$ the Pochhammer symbol,
$j_n(x) = \frac{x^{n}}{(2n+1)!!}\sum_{s=0}^{\infty}\frac{(-1)^{s}}{s!\,\bigl(n+\tfrac32\bigr)_s}\Bigl(\frac{x}{2}\Bigr)^{2s} , \qquad y_n(x) = -\frac{(2n-1)!!}{x^{n+1}}\sum_{s=0}^{\infty}\frac{(-1)^{s}}{s!\,\bigl(\tfrac12-n\bigr)_s}\Bigl(\frac{x}{2}\Bigr)^{2s} .$
Elementary closed forms. Every spherical Bessel function is a finite combination of $\sin x$ , $\cos x$ and powers of $1/x$ :
$j_0(x) = \frac{\sin x}{x} , \qquad y_0(x) = -\frac{\cos x}{x} , \qquad h^{(1)}_0(x) = -\frac{i}{x}e^{ix} , \qquad h^{(2)}_0(x) = \frac{i}{x}e^{-ix} ,$ $j_1(x) = \frac{\sin x}{x^{2}} - \frac{\cos x}{x} , \qquad j_2(x) = \Bigl(\frac{3}{x^{3}} - \frac{1}{x}\Bigr)\sin x - \frac{3}{x^{2}}\cos x ,$ $y_1(x) = -\frac{\cos x}{x^{2}} - \frac{\sin x}{x} , \qquad y_2(x) = -\Bigl(\frac{3}{x^{3}} - \frac{1}{x}\Bigr)\cos x - \frac{3}{x^{2}}\sin x .$
In general the spherical Hankel function is given by a terminating sum,
$h^{(1)}_n(x) = (-i)^{n+1}\,\frac{e^{ix}}{x}\sum_{s=0}^{n}\frac{i^{s}(n+s)!}{s!\,(2x)^{s}\,(n-s)!} ,$
with $h^{(2)}_n = \overline{h^{(1)}_n}$ for real $x$ ; $j_n$ and $y_n$ are the real and imaginary parts of $h^{(1)}_n$ .
Recurrences. For $f_n$ any of the four kinds,
$f_{n-1}(x) + f_{n+1}(x) = \frac{2n+1}{x}f_n(x) , \qquad n\,f_{n-1}(x) - (n+1)f_{n+1}(x) = (2n+1)f_n'(x) ,$ $\frac{d}{dx}\bigl[x^{n+1}f_n(x)\bigr] = x^{n+1}f_{n-1}(x) , \qquad \frac{d}{dx}\bigl[x^{-n}f_n(x)\bigr] = -x^{-n}f_{n+1}(x) .$
Rayleigh formulas. For every $n \ge 0$ ,
$j_n(x) = (-1)^{n}x^{n}\Bigl(\frac{1}{x}\frac{d}{dx}\Bigr)^{n}\frac{\sin x}{x} , \qquad y_n(x) = -(-1)^{n}x^{n}\Bigl(\frac{1}{x}\frac{d}{dx}\Bigr)^{n}\frac{\cos x}{x} ,$ $h^{(1)}_n(x) = -i\,(-1)^{n}x^{n}\Bigl(\frac{1}{x}\frac{d}{dx}\Bigr)^{n}\frac{e^{ix}}{x} , \qquad h^{(2)}_n(x) = i\,(-1)^{n}x^{n}\Bigl(\frac{1}{x}\frac{d}{dx}\Bigr)^{n}\frac{e^{-ix}}{x} .$
Small argument. For $x \ll 1$ ,
$j_n(x) \approx \frac{x^{n}}{(2n+1)!!} , \qquad y_n(x) \approx -\frac{(2n-1)!!}{x^{n+1}} ,$
so $j_n$ is regular at the origin and vanishes there like $x^{n}$ (and $j_0(0) = 1$ ), while $y_n$ diverges like $x^{-n-1}$ . The spherical Hankel functions diverge at the same rate, going as $\pm i\,y_n(x)$ : $h^{(1)}_n \approx i\,y_n$ and $h^{(2)}_n \approx -i\,y_n$ , since the regular part $j_n$ is negligible there.
Large argument. For $x \gg n(n+1)/2$ ,
$j_n(x) \sim \frac{1}{x}\sin\Bigl(x - \frac{n\pi}{2}\Bigr) , \qquad y_n(x) \sim -\frac{1}{x}\cos\Bigl(x - \frac{n\pi}{2}\Bigr) ,$ $h^{(1)}_n(x) \sim -i\,\frac{e^{i(x-n\pi/2)}}{x} , \qquad h^{(2)}_n(x) \sim i\,\frac{e^{-i(x-n\pi/2)}}{x} .$
So $j_n$ and $y_n$ describe standing spherical waves and $h^{(1)}_n$ , $h^{(2)}_n$ travelling ones: with a time factor $e^{-i\omega t}$ , $h^{(1)}_n$ is outgoing and $h^{(2)}_n$ incoming.

Definition(The Rayleigh plane-wave expansion)AdaptedAdapted from Hassani (2013), §19.7.4, Eqs. (19.43)-(19.45), pp. 594-595

Let $k > 0$ , let $\mathbf{k}$ be a vector with $|\mathbf{k}| = k$ , and let $\gamma$ be the angle between $\mathbf{k}$ and $\mathbf{r}$ , so that $\mathbf{k}\cdot\mathbf{r} = kr\cos\gamma$ . Then

e^{i\mathbf{k}\cdot\mathbf{r}} \;=\; \sum_{l=0}^{\infty}(2l+1)\,i^{\,l}\,j_l(kr)\,P_l(\cos\gamma) ,

with $P_l$ the Legendre polynomials. In particular, taking $\mathbf{k}$ along the polar axis,

e^{ikr\cos\theta} \;=\; \sum_{l=0}^{\infty}(2l+1)\,i^{\,l}\,j_l(kr)\,P_l(\cos\theta) .

Applying the addition theorem for spherical harmonics to $P_l(\cos\gamma)$ turns this into the expansion of a plane wave in spherical harmonics, which is the form used in scattering theory.

Idea of the proofOriginal

Three facts are doing all the work, and they are worth separating.

Why the functions are elementary. Bessel's equation at order $\nu$ becomes, after the substitution that removes the first-derivative term, $v'' + \bigl(1 - \tfrac{4\nu^{2}-1}{4x^{2}}\bigr)v = 0$ . At $\nu = \tfrac12$ the bracket is exactly $1$ : the equation is $v'' + v = 0$ , whose solutions are $\sin x$ and $\cos x$ . Half-integer order is the one place where Bessel's equation degenerates to the harmonic oscillator, and that is the whole reason $j_0 = \sin x/x$ rather than a transcendental new function. Everything else in the family descends from that one coincidence.

Why the Rayleigh formula is the right way down the family. Among the recurrences, the one that matters is $\frac{d}{dx}[x^{-n}f_n] = -x^{-n}f_{n+1}$ , because it expresses $f_{n+1}$ in terms of a single differentiation of $f_n$ — no second function, no second index. Rearranged it reads $f_{n+1} = -x^{n}\frac{d}{dx}(x^{-n}f_n)$ , and iterating from $n = 0$ collapses the whole chain of operators $-x^{n}\frac{d}{dx}x^{-n}$ into one clean expression $(-1)^{n}x^{n}(x^{-1}d/dx)^{n}$ . This is the same mechanism as a Rodrigues formula for orthogonal polynomials: a first-order ladder operator applied $n$ times to the seed.

Why the plane-wave expansion has the coefficients it has. A plane wave solves the Helmholtz equation, so it must be a superposition of the regular solutions $j_l(kr)$ times angular functions; the exponential is finite everywhere, which excludes $y_l$ , and it is axially symmetric about $\mathbf{k}$ , which leaves only the $m = 0$ harmonics, that is the Legendre polynomials. So the form $\sum_lA_lj_l(kr)P_l(\cos\theta)$ is forced before any computation. To find $A_l$ , project with $P_l$ and then look at small $kr$ : on the left the lowest power of $kr$ that survives the projection is $(kr)^{l}$ , because $\int_{-1}^{1}P_l(u)u^{m}\,du$ vanishes for $m < l$ by orthogonality; on the right $j_l(kr)$ also starts at $(kr)^{l}$ . Matching the two coefficients of $(kr)^{l}$ is one line of factorials, and the double factorial identity $(2l+1)! = (2l+1)!!\,2^{l}l!$ makes everything cancel down to $i^{l}(2l+1)$ .

Proof(Reduction to half-integer order, and

j_0

y_0

)AdaptedAdapted from Riley, Hobson & Bence (2006), §18.6, p. 614; Arfken, Weber & Harris (2013), §14.7, Eqs. (14.149)-(14.151), pp. 698-699; Strauss (2008), §10.5, Eq. (11) and the line after it, p. 285

Step 1 (change of dependent variable). Put $x = kr$ and write $R = x^{-1/2}Z(x)$ . Then

R' = x^{-1/2}Z' - \tfrac12 x^{-3/2}Z , \qquad R'' = x^{-1/2}Z'' - x^{-3/2}Z' + \tfrac34 x^{-5/2}Z ,

where the prime is $d/dx$ . Substituting into $x^{2}R'' + 2xR' + [x^{2}-n(n+1)]R = 0$ and collecting,

x^{3/2}Z'' + \bigl(-x^{1/2} + 2x^{1/2}\bigr)Z' + \Bigl(\tfrac34x^{-1/2} - x^{-1/2} + x^{3/2} - n(n+1)x^{-1/2}\Bigr)Z = 0 .

Multiplying by $x^{1/2}$ ,

x^{2}Z'' + xZ' + \Bigl(x^{2} - n(n+1) - \tfrac14\Bigr)Z = 0 .

Step 2 (identify the order). Since $n(n+1) + \tfrac14 = \bigl(n+\tfrac12\bigr)^{2}$ , this is Bessel's equation of order $\nu = n+\tfrac12$ . Hence $Z$ is a combination of $J_{n+1/2}$ and $Y_{n+1/2}$ , and $R = x^{-1/2}Z$ is the corresponding combination of $j_n$ and $y_n$ up to the constant $\sqrt{\pi/2}$ built into the definitions.

Step 3 (the seed functions). Half-integer order is the one case in which Bessel's equation is elementary. Removing the first-derivative term from $z^{2}u'' + zu' + (z^{2}-s^{2})u = 0$ by $u = z^{-1/2}v$ gives

v'' + \Bigl(1 - \frac{s^{2}-\tfrac14}{z^{2}}\Bigr)v \;=\; 0

(Strauss (2008), §10.5, Eq. (10), p. 284), and at $s = \tfrac12$ the bracket is exactly $1$ : the equation is $v'' + v = 0$ . So $J_{1/2}$ and $J_{-1/2}$ are $z^{-1/2}\sin z$ and $z^{-1/2}\cos z$ up to normalisation, and the standard normalisation makes this

J_{1/2}(x) = \sqrt{\frac{2}{\pi x}}\,\sin x , \qquad J_{-1/2}(x) = \sqrt{\frac{2}{\pi x}}\,\cos x

(Strauss (2008), §10.5, Eq. (11) and the sentence after it, p. 285). Therefore

j_0(x) = \sqrt{\frac{\pi}{2x}}\cdot\sqrt{\frac{2}{\pi x}}\sin x = \frac{\sin x}{x} .

For $y_0$ , note that at half-integer order the definition of $Y_\nu$ simplifies: with $\nu = n+\tfrac12$ one has $\cos\nu\pi = 0$ and $\sin\nu\pi = (-1)^{n}$ , so

Y_{n+1/2}(x) = \frac{\cos\nu\pi\,J_\nu(x) - J_{-\nu}(x)}{\sin\nu\pi}\bigg|_{\nu = n+1/2} = (-1)^{n+1}J_{-n-1/2}(x) ,

which is item 1 of the statement after multiplying by $\sqrt{\pi/2x}$ . At $n = 0$ it gives $Y_{1/2} = -J_{-1/2}$ and hence $y_0(x) = -\cos x/x$ . $\square$

Proof(The recurrences and the Rayleigh formulas)AdaptedAdapted from Arfken, Weber & Harris (2013), §14.7, Eqs. (14.169)-(14.176), pp. 702-703; Riley, Hobson & Bence (2006), §18.6, the derivation of Eq. (18.106), p. 616

Step 1 (transfer the Bessel recurrences). The Bessel functions of the first kind satisfy

J_{\nu-1} + J_{\nu+1} = \frac{2\nu}{x}J_\nu , \qquad \frac{d}{dx}\bigl[x^{\nu}J_\nu\bigr] = x^{\nu}J_{\nu-1} , \qquad \frac{d}{dx}\bigl[x^{-\nu}J_\nu\bigr] = -x^{-\nu}J_{\nu+1} ,

and so does $Y_\nu$ , hence also $H^{(1)}_\nu$ and $H^{(2)}_\nu$ , which are constant-coefficient combinations of the two. Put $\nu = n+\tfrac12$ .

For the first, multiply through by $\sqrt{\pi/2x}$ : the factor is the same at every order, and $2\nu = 2n+1$ , so

f_{n-1}(x) + f_{n+1}(x) = \frac{2n+1}{x}f_n(x) .

For the third, observe that $x^{-n}f_n = \sqrt{\pi/2}\;x^{-n-1/2}J_{n+1/2} = \sqrt{\pi/2}\;x^{-\nu}J_\nu$ — the half-power in the definition of $f_n$ is exactly what is needed to turn $x^{-n}$ into $x^{-\nu}$ . Hence

\frac{d}{dx}\bigl[x^{-n}f_n\bigr] = \sqrt{\pi/2}\;\frac{d}{dx}\bigl[x^{-\nu}J_\nu\bigr] = -\sqrt{\pi/2}\;x^{-\nu}J_{\nu+1} = -x^{-n}\sqrt{\frac{\pi}{2x}}J_{n+3/2} = -x^{-n}f_{n+1} .

The same device on the second gives $\frac{d}{dx}[x^{n+1}f_n] = x^{n+1}f_{n-1}$ , since $x^{n+1}f_n = \sqrt{\pi/2}\,x^{\nu}J_\nu$ . Carrying out the two differentiations by the product rule turns them into

f_{n-1} = f_n' + \frac{n+1}{x}f_n , \qquad f_{n+1} = -f_n' + \frac{n}{x}f_n ,

whose sum reproduces $f_{n-1}+f_{n+1} = \tfrac{2n+1}{x}f_n$ once more (the $f_n'$ cancels), and whose combination $n\times(\text{first}) - (n+1)\times(\text{second})$ cancels the $f_n$ terms instead and leaves the remaining recurrence $nf_{n-1} - (n+1)f_{n+1} = (2n+1)f_n'$ .

Step 2 (one step of the ladder). Rearranging the third identity,

f_{n+1}(x) = -x^{n}\frac{d}{dx}\bigl[x^{-n}f_n(x)\bigr] .

Step 3 (induction). Claim: $f_n(x) = (-1)^{n}x^{n}\bigl(\tfrac1x\tfrac{d}{dx}\bigr)^{n}f_0(x)$ for every $n \ge 0$ . For $n = 0$ both sides are $f_0$ . Assume it for $n$ . Then by Step 2,

f_{n+1} = -x^{n}\frac{d}{dx}\Bigl[x^{-n}\cdot(-1)^{n}x^{n}\bigl(\tfrac1x\tfrac{d}{dx}\bigr)^{n}f_0\Bigr] = (-1)^{n+1}x^{n}\frac{d}{dx}\bigl(\tfrac1x\tfrac{d}{dx}\bigr)^{n}f_0 ,

the $x^{-n}$ and $x^{n}$ cancelling inside the bracket. Writing $x^{n}\frac{d}{dx} = x^{n+1}\cdot\tfrac1x\tfrac{d}{dx}$ gives

f_{n+1} = (-1)^{n+1}x^{n+1}\bigl(\tfrac1x\tfrac{d}{dx}\bigr)^{n+1}f_0 ,

which is the claim at $n+1$ . Inserting the seeds $j_0 = \sin x/x$ , $y_0 = -\cos x/x$ , $h^{(1)}_0 = -ie^{ix}/x$ and $h^{(2)}_0 = ie^{-ix}/x$ from the previous proof gives the four Rayleigh formulas of item 5. Since each application of $\tfrac1x\tfrac{d}{dx}$ to a combination of $\sin x$ , $\cos x$ and inverse powers of $x$ returns another such combination, item 3 follows: every $j_n$ and $y_n$ is elementary. $\square$

Proof(The plane-wave expansion)AdaptedAdapted from Hassani (2013), §19.7.4, Eqs. (19.43)-(19.45) and the two integrals quoted on p. 595, pp. 594-595

Take $\mathbf{k} = k\hat{\mathbf{z}}$ , so that $\mathbf{k}\cdot\mathbf{r} = kr\cos\theta$ ; the general case follows by rotating the coordinate frame, which changes $\theta$ into the angle $\gamma$ between $\mathbf{k}$ and $\mathbf{r}$ and leaves the right-hand side unchanged in form. Write $u = \cos\theta$ .

Step 1 (the form of the expansion). The function $e^{i\mathbf{k}\cdot\mathbf{r}}$ satisfies $\nabla^{2}\psi + k^{2}\psi = 0$ , since $\nabla^{2}e^{i\mathbf{k}\cdot\mathbf{r}} = -|\mathbf{k}|^{2}e^{i\mathbf{k}\cdot\mathbf{r}}$ . Hence it is a superposition of the separated solutions $j_l(kr)Y_{lm}(\theta,\varphi)$ and $y_l(kr)Y_{lm}(\theta,\varphi)$ . Two restrictions apply. The plane wave is finite at $r = 0$ , so by item 6 the $y_l$ terms are absent. And it is independent of $\varphi$ with this choice of $\mathbf{k}$ , so only $m = 0$ survives, and $Y_{l0}$ is a multiple of $P_l(u)$ . Therefore

e^{ikru} \;=\; \sum_{l=0}^{\infty}A_l\,j_l(kr)\,P_l(u)

for some constants $A_l$ , independent of $r$ and $u$ .

Step 2 (project with $P_n$ ). Multiply by $P_n(u)$ and integrate over $u \in [-1,1]$ , using the orthogonality $\int_{-1}^{1}P_nP_l\,du = \frac{2}{2n+1}\delta_{nl}$ of the Legendre polynomials:

\int_{-1}^{1}P_n(u)\,e^{ikru}\,du \;=\; A_n\,j_n(kr)\,\frac{2}{2n+1} , \qquad\text{so}\qquad A_n\,j_n(kr) = \frac{2n+1}{2}\int_{-1}^{1}P_n(u)\,e^{ikru}\,du .

Step 3 (expand the exponential and use orthogonality again). Insert $e^{ikru} = \sum_{m\ge0}(ikr)^{m}u^{m}/m!$ and integrate term by term (legitimate because the series converges uniformly in $u$ on $[-1,1]$ for each fixed $r$ ):

A_n\,j_n(kr) = \frac{2n+1}{2}\sum_{m=0}^{\infty}\frac{(ikr)^{m}}{m!}\int_{-1}^{1}P_n(u)\,u^{m}\,du .

Now $u^{m}$ is a polynomial of degree $m$ , and $P_n$ is orthogonal to every polynomial of degree less than $n$ ; so all the terms with $m < n$ vanish, and the series begins at $m = n$ with

\int_{-1}^{1}P_n(u)\,u^{n}\,du = \frac{2^{n+1}(n!)^{2}}{(2n+1)!}

(Hassani (2013), §19.7.4, the second integral quoted on p. 595).

Step 4 (match the lowest power of $kr$ ). Both sides of Step 3 are power series in $kr$ beginning at order $n$ : on the left, by item 6, $j_n(kr) = (kr)^{n}/(2n+1)!! + O\bigl((kr)^{n+2}\bigr)$ . Equating the coefficients of $(kr)^{n}$ ,

\frac{A_n}{(2n+1)!!} \;=\; \frac{2n+1}{2}\cdot\frac{i^{\,n}}{n!}\cdot\frac{2^{n+1}(n!)^{2}}{(2n+1)!} \;=\; \frac{(2n+1)\,i^{\,n}\,2^{n}\,n!}{(2n+1)!} .

Step 5 (simplify with the double factorial). Splitting $(2n+1)!$ into its odd and even factors, $(2n+1)! = (2n+1)!!\,(2n)!! = (2n+1)!!\,2^{n}n!$ . Hence $2^{n}n!/(2n+1)! = 1/(2n+1)!!$ and

A_n = (2n+1)!!\cdot\frac{(2n+1)\,i^{\,n}}{(2n+1)!!} = (2n+1)\,i^{\,n} .

Substituting back into Step 1 gives the expansion. $\square$

Remark(Where the hypotheses are used, and what the conventions cost)Original

Integer $n$ . The definitions make sense for any index, but two things depend on $n$ being an integer. Item 1, $y_n = (-1)^{n+1}j_{-n-1}$ , uses $\cos(n+\tfrac12)\pi = 0$ and $\sin(n+\tfrac12)\pi = (-1)^{n}$ , which fail off the integers. And the terminating sum of item 3 exists only because the asymptotic series for $H^{(1)}_\nu$ truncates at half-integer $\nu$ — see Asymptotic forms of Bessel functions. Non-integer indices are not useful here anyway: the index comes from the angular separation constant $n(n+1)$ , and single-valuedness on the sphere forces it to be a non-negative integer.

$x > 0$ . The functions $y_n$ and $h^{(j)}_n$ are singular at the origin (item 6), and $j_n$ is defined there only by its series. In a problem posed on a ball including the centre, item 6 is what excludes $y_n$ ; in a spherical shell it does not, and both must be kept.

Where the plane-wave proof uses what. Step 1 uses completeness of the separated solutions of the Helmholtz equation and regularity of the exponential at the origin; drop regularity and $y_l$ terms reappear with arbitrary coefficients. Step 2 uses Legendre orthogonality on $[-1,1]$ , which is where the weight $du = -\sin\theta\,d\theta$ of the sphere enters. Step 3 uses term-by-term integration, justified by uniform convergence of the exponential series on the compact interval. Step 4 uses that both sides start at order $(kr)^{n}$ — if $j_n$ began at a lower power the matching would be between different powers and would give nothing. Nothing in the argument shows the expansion converges for large $kr$ ; what it shows is that if the expansion holds, the coefficients are these. Convergence is the content of the completeness statement invoked in Step 1.

The $x \gg n(n+1)/2$ in item 7. The asymptotic forms are stated with this condition (Arfken, Weber & Harris (2013), §14.7, the line after Eq. (14.182), p. 703), and it is not a formality: at fixed $x$ the functions $j_n$ switch, as $n$ grows past about $x$ , from oscillating with amplitude $1/x$ to decaying like $x^{n}/(2n+1)!!$ . Using item 7 in the regime $n \gtrsim x$ gives answers that are wrong by many orders of magnitude. The crossover is why partial-wave sums in scattering can be truncated at $l \approx kR$ .

A numerical trap in item 4. The recurrence $f_{n+1} = \frac{2n+1}{x}f_n - f_{n-1}$ run upwards is unstable for $j_n$ when $n \gtrsim x$ , because then $j_{n+1}$ is much smaller than either term on the right and is computed as the difference of two nearly equal numbers. In the worked example below, $j_3(1)$ comes out as $0.31017525 - 0.30116867$ , already a loss of one and a half digits, and $j_2(1)$ loses two before that. For $n \gtrsim x$ one runs the recurrence downwards instead, starting from a guessed pair at high $n$ and normalising at the end against the exact $j_0$ ; downward recurrence is stable because it amplifies the wanted solution. This is our observation, not a statement taken from the sources.

Conventions to watch. Three. (i) $y_n$ versus $n_\ell$ , recorded in the definition; the functions are the same. (ii) The scale factor: $j_n$ , $y_n$ , $h^{(1)}_n$ and $h^{(2)}_n$ all carry $\sqrt{\pi/2x}$ , but the modified spherical functions do not match each other — Arfken defines $i_n(x) = \sqrt{\pi/2x}\,I_{n+1/2}(x)$ and $k_n(x) = \sqrt{2/\pi x}\,K_{n+1/2}(x)$ , and notes explicitly that the factor in $k_n$ differs from the rest of the family (Arfken, Weber & Harris (2013), §14.7, Eqs. (14.193)-(14.194) and the sentence after them, p. 706); see Modified Bessel functions for $I_\nu$ and $K_\nu$ themselves. (iii) In the plane-wave expansion the sign convention of the exponent matters: with $e^{-i\mathbf{k}\cdot\mathbf{r}}$ the coefficients become $(2l+1)(-i)^{l}$ , as conjugating the expansion shows.

Method(Solving a radial Helmholtz problem with spherical Bessel functions)Original

Given $\nabla^{2}\psi + k^{2}\psi = 0$ on a spherical domain with data on the boundary:

Separate. Write $\psi = R(r)Y_{lm}(\theta,\varphi)$ ; the angular factor is a spherical harmonic and the radial factor satisfies the spherical Bessel equation with $n = l$ .
Write the general radial solution. $R(r) = A\,j_l(kr) + B\,y_l(kr)$ on a shell; equivalently $A'h^{(1)}_l(kr) + B'h^{(2)}_l(kr)$ if the natural conditions are at infinity.
Kill one solution with the condition at the inner end. If the origin is in the domain and $\psi$ must be finite, item 6 forces $B = 0$ .
Kill one solution with the condition at the outer end. If the domain is unbounded and only outgoing radiation is allowed, item 7 forces the combination $h^{(1)}_l$ (with the time convention $e^{-i\omega t}$ ).
Impose the remaining boundary condition and read off the eigenvalues. A Dirichlet condition $R(a) = 0$ with $B = 0$ requires $j_l(ka) = 0$ , so $k = \alpha_{lp}/a$ where $\alpha_{lp}$ is the $p$ -th positive zero of $j_l$ ; these are tabulated (Arfken, Weber & Harris (2013), §14.7, Table 14.2, p. 704). For $l = 0$ they are exactly $p\pi$ , since $j_0(x) = \sin x/x$ .
Expand the data. The eigenfunctions $j_l(\alpha_{lp}r/a)$ are orthogonal on $(0,a)$ with weight $r^{2}$ — the volume element of the ball rather than the area element of the disc, which is why the weight is $r^{2}$ and not $r$ as in Orthogonality and Fourier-Bessel series (Arfken, Weber & Harris (2013), §14.7, Eq. (14.183), p. 704).
Evaluate with closed forms or the Rayleigh formula, not with the series. For small $l$ use item 3 directly; for larger $l$ use item 4 upwards if $x \gtrsim l$ and downwards otherwise (see the remark).

Example(Checking a closed form, then truncating the plane-wave expansion)Original

Part (a): $j_2(3)$ two ways.

From the Rayleigh formula with $n = 2$ : applying $\tfrac1x\tfrac{d}{dx}$ to $j_0 = \sin x/x$ gives

\frac1x\frac{d}{dx}\frac{\sin x}{x} = \frac{1}{x}\Bigl(\frac{\cos x}{x} - \frac{\sin x}{x^{2}}\Bigr) = \frac{\cos x}{x^{2}} - \frac{\sin x}{x^{3}} ,

and applying it once more,

\Bigl(\frac1x\frac{d}{dx}\Bigr)^{2}\frac{\sin x}{x} = \frac{1}{x}\Bigl(-\frac{\sin x}{x^{2}} - \frac{2\cos x}{x^{3}} - \frac{\cos x}{x^{3}} + \frac{3\sin x}{x^{4}}\Bigr) = \Bigl(\frac{3}{x^{5}} - \frac{1}{x^{3}}\Bigr)\sin x - \frac{3}{x^{4}}\cos x .

Multiplying by $(-1)^{2}x^{2} = x^{2}$ reproduces the closed form of item 3,

j_2(x) = \Bigl(\frac{3}{x^{3}} - \frac{1}{x}\Bigr)\sin x - \frac{3}{x^{2}}\cos x .

At $x = 3$ , with $\sin 3 = 0.14112001$ and $\cos 3 = -0.98999250$ :

j_2(3) = \Bigl(\frac{3}{27} - \frac13\Bigr)(0.14112001) - \frac{3}{9}(-0.98999250) = (-0.22222222)(0.14112001) + 0.32999750 = 0.29863750 .

Independently, from the recurrence of item 4 with $n = 1$ , that is $j_2 = \frac{3}{x}j_1 - j_0$ : at $x = 3$ one has $j_0(3) = \sin 3/3 = 0.04704000$ and $j_1(3) = \sin3/9 - \cos3/3 = 0.01568000 + 0.32999750 = 0.34567750$ , so

j_2(3) = \frac{3}{3}(0.34567750) - 0.04704000 = 0.29863750 .

The two routes agree to all eight digits shown, and the exact value is $j_2(3) = 0.2986375$ .

Part (b): how fast the Rayleigh expansion converges.

Take $kr = 1$ and $\cos\gamma = \tfrac12$ , so that the exact answer is $e^{i/2} = \cos\tfrac12 + i\sin\tfrac12 = 0.8775826 + 0.4794255\,i$ .

The needed function values come from item 3 and then item 4 upwards, with $\sin1 = 0.84147098$ and $\cos1 = 0.54030231$ :

j_0(1) = \sin 1 = 0.84147098 , \qquad j_1(1) = \sin1 - \cos1 = 0.30116867 ,

j_2(1) = 3j_1(1) - j_0(1) = 0.90350601 - 0.84147098 = 0.06203503 ,

j_3(1) = 5j_2(1) - j_1(1) = 0.31017525 - 0.30116867 = 0.00900658 .

(The true values are $j_2(1) = 0.06203505$ and $j_3(1) = 0.00900658$ : the eight digits carried through the subtractions have already become six in $j_2$ , which is the instability noted in the remark.)

The Legendre values at $u = \tfrac12$ are $P_0 = 1$ , $P_1 = u = 0.5$ , $P_2 = \tfrac{3u^{2}-1}{2} = -0.125$ and $P_3 = \tfrac{5u^{3}-3u}{2} = -0.4375$ . The terms $(2l+1)i^{\,l}j_l(1)P_l(\tfrac12)$ are therefore

$l = 0$ : $1\cdot1\cdot0.84147098\cdot1 = +0.8414710$ ;
$l = 1$ : $3\cdot i\cdot0.30116867\cdot0.5 = +0.4517530\,i$ ;
$l = 2$ : $5\cdot(-1)\cdot0.06203505\cdot(-0.125) = +0.0387719$ ;
$l = 3$ : $7\cdot(-i)\cdot0.00900658\cdot(-0.4375) = +0.0275827\,i$ .

Summing,

\sum_{l=0}^{3} \;=\; 0.8802429 + 0.4793357\,i ,

against the exact $0.8775826 + 0.4794255\,i$ : an error of modulus $2.7\times10^{-3}$ . Continuing with $j_4(1) = 0.00101102$ and $j_5(1) = 0.00009256$ (computed from the series of item 2, which is stable where the upward recurrence is not), together with $P_4(\tfrac12) = -0.2890625$ and $P_5(\tfrac12) = 0.0898438$ , the next two terms are

9\cdot1\cdot0.00101102\cdot(-0.2890625) = -0.0026302 , \qquad 11\cdot i\cdot0.00009256\cdot0.0898438 = +0.0000915\,i ,

giving

\sum_{l=0}^{5} \;=\; 0.8776127 + 0.4794272\,i ,

an error of modulus $3.0\times10^{-5}$ , a hundredfold improvement for two more terms. The reason is item 6: at $kr = 1$ the factor $j_l(1) \approx 1/(2l+1)!!$ shrinks by roughly an order of magnitude at each step, as the list $0.841$ , $0.301$ , $0.0620$ , $0.00901$ , $0.00101$ , $0.0000926$ shows. All function values are rounded from a 25-digit computation.

Application(Scattering, radiation, and bound states in a sphere)AdaptedAdapted from Arfken, Weber & Harris (2013), §14.7, Example 14.7.1, pp. 705-706, and Eqs. (14.179)-(14.182), p. 703; Hassani (2013), §19.7.3, the hard-sphere energy levels, pp. 593-594, and Problem 22.35, p. 698

Partial-wave analysis. The reason the Rayleigh expansion matters in physics is that it converts an incoming plane wave into a sum over angular momenta, each term of which a spherically symmetric scatterer treats independently. Writing the incident wave as $\sum_l(2l+1)i^{l}j_l(kr)P_l(\cos\theta)$ and the scattered field as a sum of outgoing $h^{(1)}_l(kr)$ , the effect of the target on the $l$ -th channel reduces, by item 7, to a single phase shift: $j_l(kr) \sim \sin(kr - l\pi/2)/kr$ becomes $\sin(kr - l\pi/2 + \delta_l)/kr$ . Because $j_l(kr)$ is tiny for $l \gtrsim kR$ (item 6, and the crossover discussed in the remark), only a finite number of channels contribute to scattering from a target of radius $R$ , and that is what makes the method practical.

Multipole radiation. The outgoing solutions $h^{(1)}_l(kr)Y_{lm}$ are the electric and magnetic multipole fields. Item 7 gives their far field, $h^{(1)}_l(kr) \sim -i\,e^{i(kr - l\pi/2)}/kr$ , so every multipole radiates as $1/r$ with an $l$ -dependent phase; item 6 gives their behaviour near the source, $h^{(1)}_l \approx i\,y_l \approx -i(2l-1)!!/(kr)^{l+1}$ , which is the static $2^{l}$ -pole field. The single function $h^{(1)}_l$ interpolates between the static multipole and the radiation field.

The Helmholtz Green's function. The free Green's function of $\nabla^{2}+k^{2}$ in three dimensions decomposes as a sum over $l$ and $m$ of $j_l(kr_<)h^{(1)}_l(kr_>)$ times spherical harmonics of the two directions, the regular solution being used for the inner point and the outgoing one for the outer Hassani (Hassani (2013), Problem 22.35, p. 698). Letting the outer point go to infinity in that identity reproduces the plane-wave expansion, now in spherical-harmonic form; the construction itself is in Green's functions for the Helmholtz equation.

Bound states in a sphere. For a particle confined to a hard sphere of radius $a$ , the Schrödinger equation is the Helmholtz equation with $k^{2} = 2mE/\hbar^{2}$ , and the method above applies verbatim: regularity at the centre leaves $j_l$ , the wall condition requires $j_l(ka) = 0$ , and the energies are $E_{lp} = \hbar^{2}\alpha_{lp}^{2}/(2ma^{2})$ with $\alpha_{lp}$ the zeros of $j_l$ (Hassani (2013), §19.7.3, pp. 593-594); Arfken (Arfken, Weber & Harris (2013), Example 14.7.1, pp. 705-706) works the ground state. The lowest level is $l = 0$ , where $\alpha_{01} = \pi$ exactly because $j_0 = \sin x/x$ , giving $E_{\min} = \pi^{2}\hbar^{2}/(2ma^{2})$ . The pattern — a zero-point energy set by the first zero of $j_0$ , a discrete spectrum, and an $l$ -dependence that pushes higher angular momenta up in energy — is the generic shape of a three-dimensional bound-state problem. The related Schrödinger problems with a finite well use the modified spherical functions outside the well, where the solution must decay.

Theorem

Asymptotic forms of Bessel functions

Adapted

Add to PDF

At small argument every Bessel function is a power or a logarithm fixed by its leading series term; at large argument all of them are $x^{-1/2}$ times an exponential, with the whole family generated by two divergent-but-asymptotic series $P_\nu$ and $Q_\nu$ in $1/x$ — so that $J_\nu(x) \sim \sqrt{2/\pi x}\cos(x - \nu\pi/2 - \pi/4)$ and $K_\nu(x)\sim\sqrt{\pi/2x}\,e^{-x}$ .

Prerequisites: Neumann and Hankel functions, Modified Bessel functions, Integral representations of Bessel functions, The method of steepest descents, Asymptotic series

Theorem(Small argument)AdaptedAdapted from Arfken, Weber & Harris (2013), §14.3, Eq. (14.58) and Eq. (14.62), p. 668, and Eq. (14.66), p. 670; Arfken, Weber & Harris (2013), §14.4, Eqs. (14.79)-(14.82), p. 675; Arfken, Weber & Harris (2013), §14.5, Eq. (14.100), p. 681, and Eqs. (14.110)-(14.111), p. 683

As $x \to 0^{+}$ , with $\gamma$ the Euler-Mascheroni constant, the leading behaviour of each function is the leading term of its series (or, where the series is logarithmic, the logarithm):

J_\nu(x) \to \frac{1}{\Gamma(1+\nu)}\Bigl(\frac{x}{2}\Bigr)^{\nu} , \qquad J_\nu'(x) \to \frac{\nu}{2\Gamma(1+\nu)}\Bigl(\frac{x}{2}\Bigr)^{\nu-1} ,

Y_\nu(x) = -\frac{\Gamma(\nu)}{\pi}\Bigl(\frac{2}{x}\Bigr)^{\nu} + \cdots \quad(\nu > 0) , \qquad Y_0(x) = \frac{2}{\pi}\Bigl[\gamma + \ln\frac{x}{2}\Bigr] + \cdots ,

H^{(1)}_0(x) \approx 1 + i\frac{2}{\pi}\bigl(\ln x + \gamma - \ln 2\bigr) , \qquad H^{(1)}_\nu(x) \approx -i\,\frac{\Gamma(\nu)}{\pi}\Bigl(\frac{2}{x}\Bigr)^{\nu} \quad(\nu > 0) ,

with $H^{(2)}_\nu$ the complex conjugate of these for real $x$ , and

I_\nu(x) = \frac{x^{\nu}}{2^{\nu}\Gamma(\nu+1)} + \cdots , \qquad K_0(x) = -\ln x - \gamma + \ln 2 + \cdots , \qquad K_\nu(x) = 2^{\nu-1}\Gamma(\nu)\,x^{-\nu} + \cdots \quad(\nu > 0) .

The form of $Y_\nu$ given above is derived for non-integer $\nu$ ; for integer $n > 0$ the series contains a term $\tfrac{2}{\pi}J_n(x)\ln(x/2)$ as well, but the most divergent term is the same (Arfken, Weber & Harris (2013), §14.3, after Eq. (14.61), p. 668).

The pattern is uniform: in each pair exactly one member is regular at the origin ( $J_\nu$ and $I_\nu$ , vanishing like $x^{\nu}$ for $\nu>0$ and tending to $1$ for $\nu = 0$ ) and one is singular ( $Y_\nu$ and $K_\nu$ , diverging like $x^{-\nu}$ for $\nu>0$ and like a logarithm at $\nu = 0$ ). This is what makes "bounded at the origin" a usable boundary condition.

Theorem(Large argument)AdaptedAdapted from Arfken, Weber & Harris (2013), §14.6, Eqs. (14.133)-(14.143), pp. 692-693; Strauss (2008), §10.5, Eq. (5), p. 283, and Eqs. (14)-(15), p. 285

Write $\mu = 4\nu^{2}$ and define the two formal series

P_\nu(z) \sim 1 - \frac{(\mu-1)(\mu-9)}{2!\,(8z)^{2}} + \frac{(\mu-1)(\mu-9)(\mu-25)(\mu-49)}{4!\,(8z)^{4}} - \cdots ,

Q_\nu(z) \sim \frac{\mu-1}{1!\,(8z)} - \frac{(\mu-1)(\mu-9)(\mu-25)}{3!\,(8z)^{3}} + \cdots ,

and abbreviate the phase $\chi = z - \bigl(\nu+\tfrac12\bigr)\tfrac{\pi}{2} = z - \tfrac{\nu\pi}{2} - \tfrac{\pi}{4}$ . Then, as $|z| \to \infty$ :

H^{(1)}_\nu(z) = \sqrt{\frac{2}{\pi z}}\;e^{i\chi}\bigl[P_\nu(z) + iQ_\nu(z)\bigr] , \qquad -\pi < \arg z < 2\pi ,

H^{(2)}_\nu(z) = \sqrt{\frac{2}{\pi z}}\;e^{-i\chi}\bigl[P_\nu(z) - iQ_\nu(z)\bigr] , \qquad -2\pi < \arg z < \pi ,

J_\nu(z) = \sqrt{\frac{2}{\pi z}}\Bigl[P_\nu(z)\cos\chi - Q_\nu(z)\sin\chi\Bigr] , \qquad Y_\nu(z) = \sqrt{\frac{2}{\pi z}}\Bigl[P_\nu(z)\sin\chi + Q_\nu(z)\cos\chi\Bigr] ,

both for $-\pi < \arg z < \pi$ , and for the modified functions

K_\nu(z) = \sqrt{\frac{\pi}{2z}}\;e^{-z}\bigl[P_\nu(iz) + iQ_\nu(iz)\bigr] , \qquad I_\nu(z) = \frac{e^{z}}{\sqrt{2\pi z}}\bigl[P_\nu(iz) - iQ_\nu(iz)\bigr] ,

the latter for $-\pi/2 < \arg z < \pi/2$ . The $K_\nu$ expansion holds for $-3\pi/2 < \arg z < 3\pi/2$ and reads, written out,

K_\nu(z) \sim \sqrt{\frac{\pi}{2z}}\;e^{-z}\left[1 + \frac{4\nu^{2}-1^{2}}{1!\,8z} + \frac{(4\nu^{2}-1^{2})(4\nu^{2}-3^{2})}{2!\,(8z)^{2}} + \cdots\right] ,

in which every sign is $+$ : the alternating signs written into $P_\nu$ and $Q_\nu$ are cancelled by the rotation of the argument from $z$ to $iz$ , so the series for $P_\nu(iz)$ and $iQ_\nu(iz)$ carry all positive signs (Arfken, Weber & Harris (2013), §14.6, the remark after Eq. (14.136), p. 692). (The terms are positive only when every factor $4\nu^{2}-(2k-1)^{2}$ is; for $\nu$ of moderate size the early terms alternate.) Equivalently, in closed form,

K_\nu(z) \sim \sqrt{\frac{\pi}{2z}}\;e^{-z}\sum_{r=0}^{\infty}\frac{\Gamma\bigl(\nu+r+\tfrac12\bigr)}{r!\,\Gamma\bigl(\nu-r+\tfrac12\bigr)}\,(2z)^{-r} .

These series diverge. Each is an asymptotic series in Poincaré's sense: truncating after $N$ terms gives an error of the order of the first omitted term, so for large enough $z$ a few terms give high accuracy, but the sum over all terms does not converge for any $z$ . The symbol $\sim$ is used accordingly, and the equalities above are to be read term by term in $1/z$ .

Keeping only the leading term, and recording the error order, gives the statement most often used:

J_\nu(z) = \sqrt{\frac{2}{\pi z}}\cos\Bigl(z - \frac{\nu\pi}{2} - \frac{\pi}{4}\Bigr) + O\bigl(z^{-3/2}\bigr) , \qquad Y_\nu(z) = \sqrt{\frac{2}{\pi z}}\sin\Bigl(z - \frac{\nu\pi}{2} - \frac{\pi}{4}\Bigr) + O\bigl(z^{-3/2}\bigr) ,

H^{(1),(2)}_\nu(z) = \sqrt{\frac{2}{\pi z}}\;e^{\pm i(z - \nu\pi/2 - \pi/4)} + O\bigl(z^{-3/2}\bigr) , \qquad K_\nu(z) \sim \sqrt{\frac{\pi}{2z}}e^{-z} , \qquad I_\nu(z) \sim \frac{e^{z}}{\sqrt{2\pi z}} ,

the first three as stated by Strauss (Strauss (2008), §10.5, Eqs. (5), (14) and (15), pp. 283, 285).

Idea of the proofOriginal

Why $x^{-1/2}$ times an exponential, before any computation. Substituting $y = x^{-1/2}v$ into Bessel's equation removes the first-derivative term and leaves

v'' + \Bigl(1 - \frac{4\nu^{2}-1}{4x^{2}}\Bigr)v = 0

Hassani (Hassani (2013), Example 14.4.13, p. 432). At large $x$ the bracket is $1 + O(x^{-2})$ , so $v$ behaves like $\cos x$ or $\sin x$ , and $y = x^{-1/2}v$ like $x^{-1/2}$ times a sinusoid. The same substitution in the modified equation gives $v'' - (1 + \cdots)v = 0$ and hence $x^{-1/2}e^{\pm x}$ . That is the whole shape of the answer. What this argument cannot supply is the connection: the constant in front and the phase $-\nu\pi/2-\pi/4$ depend on which solution one started from at the origin, and a local analysis at infinity cannot see the origin.

How the connection is made. This is the reason the derivation goes through an integral representation rather than through the differential equation. A contour integral representing $H^{(1)}_\nu$ is a single formula valid at both ends: it defines the function unambiguously (so the normalisation is not lost), and its large- $x$ behaviour can be extracted by steepest descent. The exponent $\tfrac{x}{2}(z - z^{-1})$ has saddle points at $z = \pm i$ ; the contour for $H^{(1)}$ can be pushed over the one at $z = +i$ , where the exponent takes the value $ix$ — and there is the $e^{ix}$ . The phase $-\nu\pi/2 - \pi/4$ is assembled from three small factors: the $z^{-\nu-1}$ in the integrand evaluated at $z = i$ , the $1/\pi i$ in front, and the direction $3\pi/4$ in which the path crosses the saddle.

Why a divergent series rather than a convergent one. The point $z = \infty$ is an irregular singular point of Bessel's equation, so there is no theorem guaranteeing a convergent series in $1/z$ there, and in fact there is none. The mechanism is visible in the derivation of the $K_\nu$ expansion: one writes $K_\nu$ as $\int_0^\infty e^{-t}t^{\nu-1/2}(1+t/2z)^{\nu-1/2}dt$ times an elementary factor, expands the last bracket by the binomial theorem, and integrates term by term. The binomial series converges only for $t < 2|z|$ , but the integral runs to $t = \infty$ ; the exponential $e^{-t}$ keeps each individual term finite, but the interchange is illegitimate and the price is a series that is asymptotic and not convergent (Arfken, Weber & Harris (2013), §14.6, the paragraph after Eq. (14.133), p. 692). This is the same mechanism as in Watson's lemma.

Why one computation suffices for six functions. $K_\nu$ is defined from $H^{(1)}_\nu$ at a rotated argument, $J_\nu$ and $Y_\nu$ are the real and imaginary parts of $H^{(1)}_\nu$ on the real axis, $H^{(2)}_\nu$ is its conjugate, and $I_\nu$ is $i^{-\nu}J_\nu(ix)$ . So one expansion — the one for $K_\nu$ , which is the easiest because its integral representation has no oscillation in it — produces all the others by substitution.

Proof(The leading Hankel asymptotics by steepest descent)AdaptedAdapted from Arfken, Weber & Harris (2013), §14.6, Eqs. (14.119)-(14.126), pp. 688-689

Step 1 (the representation). Take as given the contour representation of Integral representations of Bessel functions,

H^{(1)}_\nu(t) = \frac{1}{\pi i}\int_{C_1}e^{(t/2)(z - 1/z)}\,z^{-\nu-1}\,dz ,

with $C_1$ a contour running from $z = 0$ out to $\infty e^{i\pi}$ ; the second Hankel function has the same integrand on the reflected contour $C_2$ (Arfken, Weber & Harris (2013), §14.6, Eqs. (14.119)-(14.120), p. 688). Take $t > 0$ large.

Step 2 (locate the saddle). Split the integrand into a slowly varying factor $g(z) = z^{-\nu-1}$ and an exponential $e^{w}$ with

w(z) = \frac{t}{2}\Bigl(z - \frac1z\Bigr) , \qquad w'(z) = \frac{t}{2}\Bigl(1 + \frac{1}{z^{2}}\Bigr) ,

so $w'(z) = 0$ at $z = \pm i$ . The contour $C_1$ can be deformed through $z_0 = +i$ and not through $-i$ . At that saddle,

w(i) = \frac{t}{2}\Bigl(i - \frac1i\Bigr) = \frac{t}{2}(i + i) = it , \qquad w''(z) = -\frac{t}{z^{3}} , \qquad w''(i) = -\frac{t}{i^{3}} = -it .

Step 3 (choose the direction of descent). The saddle-point formula of The method of steepest descents gives

\int_{C}g(z)e^{w(z)}\,dz \;\approx\; g(z_0)\,e^{w(z_0)}\,e^{i\theta}\sqrt{\frac{2\pi}{|w''(z_0)|}} , \qquad \theta = -\frac{\arg w''(z_0)}{2} + \frac{\pi}{2}\ \text{ or }\ \frac{3\pi}{2} .

Here $\arg w''(i) = \arg(-it) = -\pi/2$ , so $\theta = \pi/4 + \pi/2 = 3\pi/4$ or $\theta = \pi/4 + 3\pi/2 = 7\pi/4$ . The value $7\pi/4 \equiv -\pi/4$ is not available: the contour cannot be brought to cross the saddle in that direction without first passing through a region where $|e^{w}|$ exceeds its value at the saddle, which would make the saddle not dominant. So $\theta = 3\pi/4$ .

Step 4 (assemble). With $g(i) = i^{-\nu-1} = e^{i(\pi/2)(-\nu-1)}$ and $|w''(i)| = t$ ,

H^{(1)}_\nu(t) \approx \frac{1}{\pi i}\;e^{i(\pi/2)(-\nu-1)}\;e^{it}\;e^{3i\pi/4}\;\sqrt{\frac{2\pi}{t}} .

Collect the modulus: $\frac{1}{\pi}\sqrt{2\pi/t} = \sqrt{2/\pi t}$ . Collect the phases, using $1/i = e^{-i\pi/2}$ :

-\frac{\pi}{2} \;-\; \frac{(\nu+1)\pi}{2} \;+\; t \;+\; \frac{3\pi}{4} \;=\; t - \frac{\nu\pi}{2} - \frac{\pi}{2} - \frac{\pi}{2} + \frac{3\pi}{4} \;=\; t - \frac{\nu\pi}{2} - \frac{\pi}{4} .

Hence

H^{(1)}_\nu(t) \approx \sqrt{\frac{2}{\pi t}}\;e^{i(t - \nu\pi/2 - \pi/4)} ,

and the saddle at $z = -i$ gives, by the same computation with the conjugate phases,

H^{(2)}_\nu(t) \approx \sqrt{\frac{2}{\pi t}}\;e^{-i(t - \nu\pi/2 - \pi/4)} .

Taking real and imaginary parts for real $t$ and real $\nu$ gives the leading forms of $J_\nu$ and $Y_\nu$ . $\square$

Proof(The full expansion, from an integral representation for

K_\nu

)AdaptedAdapted from Arfken, Weber & Harris (2013), §14.6, Eqs. (14.128)-(14.143), pp. 690-693

Step 1 (the representation). For $\nu > -\tfrac12$ and $\operatorname{Re}z > 0$ ,

K_\nu(z) = \frac{\pi^{1/2}}{\Gamma\bigl(\nu+\tfrac12\bigr)}\Bigl(\frac{z}{2}\Bigr)^{\nu}\int_{1}^{\infty}e^{-zx}\bigl(x^{2}-1\bigr)^{\nu-1/2}\,dx

(Arfken, Weber & Harris (2013), §14.6, Eq. (14.128), p. 690). (That this integral solves the modified Bessel equation, has the right behaviour at the origin and decays at infinity are the three checks that identify it with $K_\nu$ ; Arfken carries them out on pp. 690-692, the third of them being a by-product of this proof.)

Step 2 (move the exponential outside). Substitute $x = 1 + t/z$ , so $dx = dt/z$ and $x$ runs from $1$ to $\infty$ as $t$ runs from $0$ to $\infty$ . Then $e^{-zx} = e^{-z}e^{-t}$ and

x^{2}-1 = \frac{t^{2}}{z^{2}} + \frac{2t}{z} = \frac{2t}{z}\Bigl(1 + \frac{t}{2z}\Bigr) ,

so $(x^{2}-1)^{\nu-1/2} = (2t/z)^{\nu-1/2}(1+t/2z)^{\nu-1/2}$ . Collecting the powers of $z$ — the prefactor supplies $z^{\nu}2^{-\nu}$ , the bracket $2^{\nu-1/2}z^{-\nu+1/2}$ , the measure $z^{-1}$ — everything combines to $2^{-1/2}z^{-1/2}$ , and

K_\nu(z) = \sqrt{\frac{\pi}{2z}}\;\frac{e^{-z}}{\Gamma\bigl(\nu+\tfrac12\bigr)}\int_{0}^{\infty}e^{-t}\,t^{\nu-1/2}\Bigl(1 + \frac{t}{2z}\Bigr)^{\nu-1/2}dt .

The exponential decay $e^{-z}$ is now explicit and outside the integral, which is the point of the substitution.

Step 3 (expand and integrate term by term). Apply the binomial series to the last factor,

\Bigl(1+\frac{t}{2z}\Bigr)^{\nu-1/2} = \sum_{r=0}^{\infty}\binom{\nu-\tfrac12}{r}\Bigl(\frac{t}{2z}\Bigr)^{r} , \qquad \binom{\nu-\tfrac12}{r} = \frac{\Gamma\bigl(\nu+\tfrac12\bigr)}{r!\,\Gamma\bigl(\nu-r+\tfrac12\bigr)} ,

and integrate each term against $e^{-t}t^{\nu-1/2}$ using $\int_0^{\infty}e^{-t}t^{\nu+r-1/2}\,dt = \Gamma\bigl(\nu+r+\tfrac12\bigr)$ , the Euler integral for the Gamma function. The factor $\Gamma(\nu+\tfrac12)$ cancels and

K_\nu(z) \sim \sqrt{\frac{\pi}{2z}}\;e^{-z}\sum_{r=0}^{\infty}\frac{\Gamma\bigl(\nu+r+\tfrac12\bigr)}{r!\,\Gamma\bigl(\nu-r+\tfrac12\bigr)}(2z)^{-r} .

Step 4 (why $\sim$ and not $=$ ). The binomial series converges only for $t < 2|z|$ , but the integration runs to $t = \infty$ . The interchange of sum and integral is therefore not justified, and the resulting series diverges for every $z$ . What survives is the asymptotic property: the $e^{-t}$ weight concentrates the integral near $t = O(1)$ , where the expansion is valid for large $|z|$ , so truncating after $N$ terms leaves an error of the size of the next term.

Step 5 (write it with $P_\nu$ and $Q_\nu$ ). Expanding the ratio of Gamma functions, $\Gamma(\nu+r+\tfrac12)/\Gamma(\nu-r+\tfrac12)$ is the product of the $2r$ factors $(\nu+r-\tfrac12)(\nu+r-\tfrac32)\cdots(\nu-r+\tfrac12)$ , which pairs off as $\prod_{k=1}^{r}\bigl(\nu^{2}-(k-\tfrac12)^{2}\bigr) = 4^{-r}\prod_{k=1}^{r}\bigl(\mu - (2k-1)^{2}\bigr)$ with $\mu = 4\nu^{2}$ . Hence the $r$ -th term is $\prod_{k=1}^{r}\bigl(\mu-(2k-1)^{2}\bigr)\big/\bigl(r!\,(8z)^{r}\bigr)$ , which is the written-out form in the statement. Splitting it into even $r$ (real when $z$ is replaced by $iz$ ) and odd $r$ (imaginary) is exactly the split into $P_\nu(iz)$ and $iQ_\nu(iz)$ defined in the statement.

Step 6 (the other five functions). Rearranging the definition $K_\nu(x) = \tfrac{\pi}{2}i^{\nu+1}H^{(1)}_\nu(ix)$ of Modified Bessel functions gives

H^{(1)}_\nu(z) = \frac{2}{\pi}\,e^{-(i\pi/2)(\nu+1)}\,K_\nu(-iz) ,

and substituting the expansion of Step 5, with $-i$ for the argument, produces the stated $H^{(1)}_\nu$ expansion; analytic continuation widens the sector in which it holds to $-\pi < \arg z < 2\pi$ . Then $H^{(2)}_\nu$ is the complex conjugate for real argument, $J_\nu$ and $Y_\nu$ are the real and imaginary parts of $H^{(1)}_\nu$ for real $z$ and real $\nu$ , which gives the $\cos\chi$ , $\sin\chi$ combinations, and $I_\nu(z) = i^{-\nu}J_\nu(iz)$ gives the last one. $\square$

Remark(Where the hypotheses are used, and what the sources differ about)AdaptedAdapted from Arfken, Weber & Harris (2013), §14.6, the sector statements with Eqs. (14.138)-(14.143), pp. 692-693, and the properties of the asymptotic forms, including the accuracy criterion Eq. (14.145) and the termination at half-integer order, pp. 693-694; Hassani (2013), §16.2.1, pp. 504-505

The sectors matter. Each expansion is quoted with a range of $\arg z$ , and they are not the same range. $H^{(1)}_\nu$ is valid on $-\pi < \arg z < 2\pi$ and $H^{(2)}_\nu$ on $-2\pi < \arg z < \pi$ ; $J_\nu$ and $Y_\nu$ only on $-\pi < \arg z < \pi$ ; $I_\nu$ only on the right half-plane $|\arg z| < \pi/2$ ; $K_\nu$ on $-3\pi/2 < \arg z < 3\pi/2$ , wider than the strip $|\arg z| < \pi/2$ on which the integral of Step 1 converges, because the finished series continues analytically. Using one outside its sector is the Stokes phenomenon in practice: a subdominant exponential that was invisible becomes dominant, and the formula silently loses its leading term. The narrow sector for $I_\nu$ is the clearest case — $I_\nu$ grows like $e^{z}$ on the positive axis but oscillates on the imaginary axis, where it is a Bessel function in disguise.

$\nu > -\tfrac12$ in the integral representation. Step 1 needs it for convergence of $\int_1^{\infty}(x^{2}-1)^{\nu-1/2}dx$ near $x = 1$ . The conclusion, however, holds for all $\nu$ , because $P_\nu$ and $Q_\nu$ depend on $\nu$ only through $\mu = 4\nu^{2}$ and $K_{-\nu} = K_\nu$ .

Where the divergence comes from. Step 4 above, and nowhere else: every other step is exact. This is worth saying because it tells you what the series can and cannot do. It cannot be summed; it can be truncated. The best accuracy available at a given $z$ is reached by stopping at the smallest term, and is of that term's order — the optimal truncation of Asymptotic series. The worked example below exhibits this.

Half-integer order: the series terminates. For $\nu = \pm\tfrac12, \pm\tfrac32, \dots$ one of the factors $\mu - (2k-1)^{2} = 4\nu^{2} - (2k-1)^{2}$ vanishes, so $P_\nu$ and $Q_\nu$ become polynomials in $1/z$ and the asymptotic forms become exact solutions (Arfken, Weber & Harris (2013), §14.6, p. 694). At $\nu = \tfrac12$ , $\mu = 1$ , so $P_{1/2} = 1$ and $Q_{1/2} = 0$ , and the expansion collapses to $J_{1/2}(z) = \sqrt{2/\pi z}\cos(z - \pi/2) = \sqrt{2/\pi z}\sin z$ — which is the exact closed form (Riley, Hobson & Bence (2006), §18.5.1, the worked example, p. 605). This is the cross-check that fixes every convention in the statement, and it is also the source of the elementary closed forms of Spherical Bessel functions.

How large is "large". The leading term of $J_n$ is accurate when the $Q_\nu\sin\chi$ term is negligible, that is when $8x \gg |4n^{2}-1|$ (Arfken, Weber & Harris (2013), §14.6, Eq. (14.145), p. 694). For $n = 0$ this is satisfied already around $x = 5$ ; for $n$ of order $10$ the asymptotic region begins near $x = 50$ . The condition is about $\nu^{2}/x$ , not about $x$ alone, which is why high-order Bessel functions need a different treatment (an expansion in which $\nu$ and $x$ are comparable) rather than more terms of this one.

Two sources, two regimes. We follow Arfken's $P_\nu$ , $Q_\nu$ form throughout. Hassani (Hassani (2013), §16.2.1, pp. 504-505) derives the same leading forms by a saddle-point calculation performed at fixed $\nu/x$ , writing $\nu = x\cos\beta$ , and then specialises to $x \gg \nu$ ; in that specialisation the first correction appears as a $\nu$ -independent $1/8ix$ , and his two-term results read $J_\nu(x)\approx\sqrt{2/\pi x}\bigl[\cos\chi + \tfrac{1}{8x}\sin\chi\bigr]$ and $Y_\nu(x)\approx\sqrt{2/\pi x}\bigl[\sin\chi - \tfrac{1}{8x}\cos\chi\bigr]$ . These agree with Arfken term by term at $\nu = 0$ , where $Q_0 = -1/8z$ . They are not in conflict for general $\nu$ either: Hassani's expansion parameter is $1/(x\sin\beta) = 1/\sqrt{x^{2}-\nu^{2}}$ and his phase is $x\sin\beta - \nu\beta$ , and expanding both for fixed $\nu$ and large $x$ restores the $\nu$ -dependence that Arfken's $Q_\nu = (4\nu^{2}-1)/8z$ carries explicitly. The practical point: if you use Hassani's two-term formulas with the simplified phase $\chi$ , they are reliable only for $\nu^{2} \ll x$ ; Arfken's are uniform in $\nu$ to the order stated.

What the small-argument statement does not claim. Each line of the small-argument statement is a leading term, not an expansion with an error bound; the symbol used in Arfken is $\to$ or " $+\cdots$ ", and we have kept it. For $Y_\nu$ at integer order the omitted terms include a $J_n(x)\ln(x/2)$ piece which is not negligible relative to the next power, only relative to the leading $x^{-n}$ .

Method(Using the expansions)Original

Decide which end you are at. Compare $x$ with $\nu$ . If $x \lesssim 1$ , use the small-argument statement. If $8x \gg |4\nu^{2}-1|$ , use the large-argument one. If $x$ and $\nu$ are comparable, neither applies: this is the transition region, where the function changes character, and it needs uniform (Airy-type) asymptotics instead — see Airy functions.
Check that your argument is in the sector. For real positive arguments every formula in the statement applies. For complex arguments, read the sector off the statement before using the formula, especially for $I_\nu$ .
Compute the phase first. Almost all errors in using these formulas are errors in $\chi = x - \nu\pi/2 - \pi/4$ . Note that it depends on $\nu$ , so $J_0$ and $J_1$ are a quarter-period out of step: this is why the zeros of successive Bessel functions interlace at large argument.
Add terms until they stop shrinking, then stop. Compute the terms $t_r = \prod_{k=1}^{r}\bigl(\mu-(2k-1)^{2}\bigr)/\bigl(r!\,(8z)^{r}\bigr)$ one at a time by the ratio $t_r/t_{r-1} = \bigl(\mu-(2r-1)^{2}\bigr)/(8zr)$ . Stop at the smallest $|t_r|$ ; the error is then of that order. Adding further terms makes the answer worse.
Estimate the error as the first omitted term. This is the working rule for an asymptotic series and is what the example below confirms numerically.
Check against a special value. At $\nu = \pm\tfrac12$ the series terminates and the result must reproduce $\sqrt{2/\pi x}\sin x$ or $\sqrt{2/\pi x}\cos x$ exactly. Any sign or phase error in a hand-derived formula shows up there immediately.

Example(Two numerical checks:

J_0(10)

and

K_0(5)

)Original

Part (a): the oscillatory case. Take $\nu = 0$ , $x = 10$ , so $\mu = 0$ and $\chi = 10 - \pi/4 = 9.2146018$ , with

\cos\chi = -0.9779942 , \qquad \sin\chi = 0.2086322 , \qquad \sqrt{\frac{2}{\pi\cdot10}} = 0.2523133 .

The leading term alone gives

J_0(10) \approx (0.2523133)(-0.9779942) = -0.2467609 ,

against the true value $J_0(10) = -0.2459358$ : an error of $-8.3\times10^{-4}$ .

Now include the next order. With $\mu = 0$ ,

P_0(10) = 1 - \frac{(-1)(-9)}{2!\,(80)^{2}} = 1 - \frac{9}{12800} = 0.9992969 , \qquad Q_0(10) = \frac{-1}{80} = -0.0125 ,

J_0(10) \approx 0.2523133\bigl[(0.9992969)(-0.9779942) - (-0.0125)(0.2086322)\bigr] = 0.2523133\,(-0.9746987) = -0.2459294 .

The error has dropped to $6.4\times10^{-6}$ , a factor of $130$ for one extra order in $1/x$ — consistent with the accuracy criterion $8x \gg |4\nu^{2}-1|$ , which here reads $80 \gg 1$ .

Part (b): the exponential case, and optimal truncation. Take $\nu = 0$ , $z = 5$ , so $8z = 40$ and the prefactor is

\sqrt{\frac{\pi}{10}}\,e^{-5} = (0.5604991)(0.0067379) = 0.0037766134 .

The terms $t_r$ are generated by $t_0 = 1$ and $t_r = t_{r-1}\cdot\bigl(-(2r-1)^{2}\bigr)/(40r)$ , giving

t_1 = -0.0250000 , \quad t_2 = 0.0028125 , \quad t_3 = -0.0005859 , \quad t_4 = 0.0001794 , \quad t_5 = -0.0000727 .

The partial sums of the bracket, and the resulting approximations to $K_0(5) = 0.00369110$ , are

$1$ term: bracket $1.0000000$ , value $0.00377661$ , error $+8.55\times10^{-5}$ ;
$2$ terms: bracket $0.9750000$ , value $0.00368220$ , error $-8.90\times10^{-6}$ ;
$3$ terms: bracket $0.9778125$ , value $0.00369282$ , error $+1.72\times10^{-6}$ ;
$4$ terms: bracket $0.9772266$ , value $0.00369061$ , error $-4.91\times10^{-7}$ ;
$5$ terms: bracket $0.9774060$ , value $0.00369128$ , error $+1.86\times10^{-7}$ .

Each error is smaller than, and has the sign of, the first omitted term multiplied by the prefactor — step 5 of the method. For instance after $4$ terms the next term is $t_4 = 1.794\times10^{-4}$ , giving $0.0037766\times1.794\times10^{-4} = 6.8\times10^{-7}$ as the predicted error bound, and the actual error is $4.9\times10^{-7}$ .

The series does not keep improving. Continuing the recursion, $|t_r|$ decreases to a minimum at $r = 10$ , where $t_{10} = 1.127\times10^{-5}$ , and increases from $r = 11$ onwards ( $|t_{11}| = 1.129\times10^{-5}$ , $|t_{12}| = 1.244\times10^{-5}$ , $|t_{13}| = 1.496\times10^{-5}$ ). Truncating at the smallest term gives $0.0036911191$ against the exact $0.0036910983$ , an error of $2.1\times10^{-8}$ ; truncating at $r = 13$ instead gives an error of $3.1\times10^{-8}$ , worse. So at $z = 5$ the series can deliver between five and six correct significant figures and no more, however many terms are kept. That ceiling is the defining feature of an asymptotic series. Where it sits can be estimated: Stirling's approximation applied to $|t_r| = \bigl((2r-1)!!\bigr)^{2}/\bigl(r!\,(8z)^{r}\bigr)$ gives $|t_r| \approx C\,(r/2ez)^{r}$ , which is smallest at $r \approx 2z$ with value of order $e^{-2z}$ . At $z = 5$ that predicts a minimum near $r = 10$ of relative size $e^{-10} = 4.5\times10^{-5}$ , and the computation above found it at $r = 10$ with relative size $1.1\times10^{-5}$ .

All values are rounded from a 20-digit computation; the quoted digits are correct as shown.

Application(Radiation conditions, phase shifts, and the zeros of

J_n

)AdaptedAdapted from Arfken, Weber & Harris (2013), §14.6, Example 14.6.1, p. 695, the discussion of the zeros and of the use of asymptotics in Wronskian formulas, pp. 693-694; Arfken, Weber & Harris (2013), §14.4, the role of the Hankel functions as travelling waves, p. 678

Picking the solution that radiates. This is the reason the asymptotic forms exist in a physics text at all. In a two-dimensional wave problem with a source at the origin and nothing coming in from infinity, the physical requirement is $U \sim e^{i(kr - \omega t)}$ at large $r$ . Of the solutions of Bessel's equation, only $H^{(1)}_0(kr)$ has that behaviour, by the large-argument statement; so the solution is $H^{(1)}_0(kr)e^{-i\omega t}$ , and its divergence as $r \to 0$ — logarithmic, by the small-argument statement — is exactly what a point source should look like (Arfken, Weber & Harris (2013), §14.6, Example 14.6.1, p. 695). The same argument in three dimensions, applied to $h^{(1)}_l$ , is what makes Spherical Bessel functions the right basis for radiation problems, and it is the content of the Sommerfeld radiation condition used in Green's functions for the Helmholtz equation.

Phase shifts in scattering. The entire phase-shift formalism rests on the leading large-argument form. A free radial wave behaves as $\sin(kr - l\pi/2)/kr$ at large $r$ ; a scattered one as $\sin(kr - l\pi/2 + \delta_l)/kr$ . The potential is encoded in the single number $\delta_l$ per channel precisely because the asymptotic form has one free phase and one free amplitude, and conservation of probability fixes the amplitude. The $-l\pi/2$ that $\delta_l$ is measured against comes from the $\nu\pi/2$ in $\chi$ with $\nu = l+\tfrac12$ .

Where the zeros are. Setting the leading term to zero predicts that the large zeros of $J_n$ are near the zeros of $\cos(x - n\pi/2 - \pi/4)$ , that is near $x = (m + \tfrac{n}{2} + \tfrac34)\pi$ , and that consecutive zeros are spaced by $\pi$ . Comparing with a table shows the agreement is already good for the first few zeros when $n$ is small (Arfken, Weber & Harris (2013), §14.6, p. 693). The same observation explains why the zeros of $J_n$ and $Y_n$ interlace at large argument: their leading terms are a cosine and a sine of the same phase. These estimates are the usual starting values for a Newton iteration when the zeros are needed accurately, as they are for the Fourier-Bessel expansions of Orthogonality and Fourier-Bessel series.

Fixing Wronskian constants. A Wronskian of two solutions of Bessel's equation is known to be $\mathrm{const}/x$ by Abel's identity, and the constant can be read off from the asymptotic forms as easily as from the behaviour at the origin — more easily, when one of the functions is singular there (Arfken, Weber & Harris (2013), §14.6, p. 694). This is a standard use of the expansions that has nothing to do with large arguments as such.

Exponential selection. For the modified functions the asymptotic statement does the same job as the regularity condition does at the origin: in a problem on an unbounded domain, $I_\nu(x) \sim e^{x}/\sqrt{2\pi x}$ disqualifies $I_\nu$ immediately and leaves $K_\nu$ . Screened potentials, cylindrical Green's functions and boundary layers all select their solution this way; see Modified Bessel functions.

Theorem · advanced

Integral representations of Bessel functions

Adapted

Add to PDF

The loop integral $J_\nu(z) = (2\pi i)^{-1}\int_C e^{(z/2)(t-1/t)}t^{-\nu-1}dt$ over a contour wrapped round the cut on the negative axis reproduces $J_\nu$ at every order when $\operatorname{Re}z>0$ , collapses to Bessel's integral $\tfrac1\pi\int_0^\pi\cos(n\theta - z\sin\theta)\,d\theta$ when $\nu$ is an integer, and splits along two halves of its contour into the Hankel functions.

Prerequisites: Bessel functions of the first kind, Neumann and Hankel functions, The Gamma function, Branch points and branch cuts, Integrals around branch cuts, Cauchy's theorem, Analytic continuation, Integrating and differentiating series term by term, Fourier series, Pointwise convergence of Fourier series

Definition(The Schläfli contour and the loop integral)AdaptedAdapted from Arfken, Weber & Harris (2013), §14.1, Eq. (14.37) and Fig. 14.5, pp. 653-654; Hassani (2013), §16.2, Eqs. (16.12)-(16.13) and Fig. 16.1, pp. 498-499

Fix $\nu \in \mathbb{C}$ and cut the $t$ -plane along the negative real axis, so that $t^{-\nu-1} = e^{(-\nu-1)\log t}$ is single valued once we fix $\arg t \in [-\pi,\pi]$ , reading $\arg t = -\pi$ on the lower edge of the cut and $\arg t = +\pi$ on the upper edge.

The Schläfli contour $C$ is the loop that comes in from $t = \infty e^{-i\pi}$ along the underside of the cut, encircles the origin once anticlockwise on a circle of some radius $\rho > 0$ , and returns to $t = \infty e^{+i\pi}$ along the upper side. For $\operatorname{Re}z > 0$ define

F_\nu(z) \;=\; \frac{1}{2\pi i}\int_{C} \exp\!\left[\frac{z}{2}\Bigl(t - \frac1t\Bigr)\right] t^{-\nu-1}\,dt .

Two remarks on the set-up.

Why the loop is open. At integer order the integrand is single valued and $C$ may be closed into a circle; the integral is then the Laurent coefficient of the generating function and returns $J_n(z)$ by the residue theorem (Bessel functions of the first kind). At non-integer order $t = 0$ is a branch point, not a pole (Branch points and branch cuts), so no closed contour round the origin is available and the circle must be opened out along a cut, exactly as in a keyhole evaluation (Integrals around branch cuts).
Why the cut is taken on the negative axis. The exponent $\tfrac{z}{2}(t - 1/t)$ has real part $\approx -\tfrac12 s\operatorname{Re}z$ on the ray $t = -s$ , $s \to +\infty$ . So for $\operatorname{Re}z > 0$ the integrand decays like $e^{-s\operatorname{Re}z/2}$ at both ends of the loop and the integral converges absolutely; on the positive axis it would instead blow up. The condition $\operatorname{Re}z > 0$ is part of the definition and cannot be dropped.

Because $\rho$ may be changed and the rays displaced without crossing a singularity, $F_\nu$ does not depend on the shape of $C$ , only on its class (Cauchy's theorem). We therefore fix $\rho = 1$ once and for all, so that $|t| \ge 1$ everywhere on the two rays; this normalisation is what the bound of Step 2 and the substitution $t = e^{w}$ of Step 4 use. Statements about "a loop of this kind" below mean a loop of this shape with any radius $\rho>0$ .

Theorem(Schläfli's representations)AdaptedAdapted from Arfken, Weber & Harris (2013), §14.1, Eqs. (14.37)-(14.40), pp. 653-654, and Eqs. (14.18)-(14.20), p. 647; Hassani (2013), §16.2, Eqs. (16.13), (16.17), (16.18), pp. 499-501; Riley, Hobson & Bence (2006), §18.5.3, Eq. (18.99), p. 613

Let $\operatorname{Re}z > 0$ and let $\nu \in \mathbb{C}$ . Then $F_\nu(z) = J_\nu(z)$ ; that is,

J_\nu(z) \;=\; \frac{1}{2\pi i}\int_{C}\exp\!\left[\frac{z}{2}\Bigl(t-\frac1t\Bigr)\right]t^{-\nu-1}\,dt . \tag{S1}

Substituting $t = e^{w}$ unrolls the loop into a strip and gives the exponential form

J_\nu(z) \;=\; \frac{1}{2\pi i}\int_{C'} e^{\,z\sinh w \,-\, \nu w}\,dw , \tag{S2}

where $C'$ runs from $w = +\infty - i\pi$ to $-i\pi$ , up the imaginary axis to $+i\pi$ , and out to $+\infty + i\pi$ ; this is the image of $C$ with its circle at the radius $\rho = 1$ fixed in the definition, the circle becoming the segment of the imaginary axis.

Evaluating (S2) leg by leg gives the real form, valid for $\operatorname{Re}z > 0$ and every $\nu \in \mathbb{C}$ ,

J_\nu(z) \;=\; \frac{1}{\pi}\int_0^{\pi}\cos\bigl(\nu\theta - z\sin\theta\bigr)\,d\theta \;-\; \frac{\sin \nu\pi}{\pi}\int_0^{\infty} e^{-\nu t \,-\, z\sinh t}\,dt . \tag{S3}

For real $\nu$ and real $z > 0$ both integrands are real, which is what makes (S3) the form to compute with.

At integer order $\nu = n$ the second term vanishes and (S3) reduces to Bessel's integral

J_n(z) \;=\; \frac{1}{\pi}\int_0^{\pi}\cos\bigl(n\theta - z\sin\theta\bigr)\,d\theta ,

which holds then for all $z$ , not only in a half-plane (Bessel functions of the first kind; (Riley, Hobson & Bence (2006), §18.5.3, Eq. (18.99), p. 613)).

Idea of the proofOriginal

Everything rests on one algebraic accident: when Bessel's operator in $z$ is applied to the integrand $e^{(z/2)(t-1/t)}t^{-\nu-1}$ , the result is an exact derivative in $t$ . So the integral solves Bessel's equation for any contour whose endpoint term vanishes, and the whole question becomes which contour, and which multiple of which solution. The Schläfli loop is the one that answers both at once: it is the circle of the generating function, opened out so that it still works when $t^{-\nu-1}$ has a branch point, and rescaling $t$ by $z$ turns it into Hankel's loop integral for $1/\Gamma$ , which delivers the series for $J_\nu$ coefficient by coefficient — normalisation included. Opening the loop costs exactly one thing, the $\sin\nu\pi$ term in (S3), which is the cut making itself felt and which disappears at integer order.

Proof(Steps 1-3: the loop integral is

J_\nu

)AdaptedAdapted from Arfken, Weber & Harris (2013), §14.1, Eqs. (14.37)-(14.40), pp. 653-654; Arfken, Weber & Harris (2013), §13.1, Eqs. (13.30)-(13.31), pp. 604-606; Hassani (2013), §16.2, pp. 498-499

Write $\Phi(z,t) = e^{(z/2)(t-1/t)}\,t^{-\nu-1}$ and let $\mathcal{B}_z = z^{2}\partial_z^{2} + z\partial_z + (z^{2}-\nu^{2})$ be Bessel's operator.

Step 1 (the integrand is annihilated up to a total derivative). Differentiating $\Phi$ in $z$ brings down a factor $\tfrac12(t-1/t)$ each time, so

\mathcal{B}_z\Phi = \Phi\left[\frac{z^{2}}{4}\Bigl(t-\frac1t\Bigr)^{2} + \frac{z}{2}\Bigl(t-\frac1t\Bigr) + z^{2} - \nu^{2}\right] .

Since $(t-1/t)^{2} + 4 = (t+1/t)^{2}$ , the two terms carrying $z^{2}$ combine into $\tfrac{z^{2}}{4}(t+1/t)^{2}$ , leaving

\mathcal{B}_z\Phi = \Phi\left[\frac{z^{2}}{4}\Bigl(t+\frac1t\Bigr)^{2} - \nu^{2} + \frac{z}{2}\Bigl(t-\frac1t\Bigr)\right] . \tag{$*$}

Now set $P(z,t) = e^{(z/2)(t-1/t)}\,t^{-\nu}\bigl[\nu + \tfrac{z}{2}(t+1/t)\bigr]$ and differentiate in $t$ by the product rule. Differentiating the prefactor,

\frac{\partial}{\partial t}\Bigl(e^{(z/2)(t-1/t)}t^{-\nu}\Bigr) = e^{(z/2)(t-1/t)}\Bigl[\frac{z}{2}\Bigl(1+\frac{1}{t^{2}}\Bigr)t^{-\nu} - \nu\,t^{-\nu-1}\Bigr] = \Phi\Bigl[\frac{z}{2}\Bigl(t+\frac1t\Bigr) - \nu\Bigr] ,

while differentiating the bracket gives $e^{(z/2)(t-1/t)}t^{-\nu}\cdot\tfrac{z}{2}(1-1/t^{2}) = \Phi\cdot\tfrac{z}{2}(t-1/t)$ . Hence

\frac{\partial P}{\partial t} = \Phi\left[\Bigl(\frac{z}{2}\Bigl(t+\frac1t\Bigr)-\nu\Bigr)\Bigl(\nu + \frac{z}{2}\Bigl(t+\frac1t\Bigr)\Bigr) + \frac{z}{2}\Bigl(t-\frac1t\Bigr)\right] ,

and the product of the two brackets is $\tfrac{z^{2}}{4}(t+1/t)^{2}-\nu^{2}$ . Comparing with $(*)$ ,

\mathcal{B}_z\Phi(z,t) \;=\; \frac{\partial P}{\partial t}(z,t) .

Step 2 (so $F_\nu$ solves Bessel's equation). Write $s = |t|$ ; on the two rays $s \ge 1$ , since the circle of $C$ has radius $1$ . Each $z$ -derivative of $\Phi$ brings down a factor $\tfrac12(t-1/t)$ , which on the rays has modulus at most $\tfrac12(s + 1/s) \le s$ , so there the integrand and its first two $z$ -derivatives are bounded in modulus by $\mathrm{const}\cdot(1+s^{2})\,e^{-s\delta/2}\,s^{-\operatorname{Re}\nu-1}$ , uniformly for $z$ in a compact subset of $\{\operatorname{Re}z \ge \delta\}$ with $\delta>0$ . That bound is integrable, so differentiation under the integral sign is legitimate and

\mathcal{B}_zF_\nu(z) = \frac{1}{2\pi i}\int_{C}\frac{\partial P}{\partial t}\,dt = \frac{1}{2\pi i}\Bigl[P(z,t)\Bigr]_{\text{start}}^{\text{end}} .

Both ends sit at $t = \infty e^{\mp i\pi}$ , where $|e^{(z/2)(t-1/t)}| = e^{-(s-1/s)\operatorname{Re}z/2} \to 0$ faster than the algebraic factors $t^{-\nu}$ and $\tfrac{z}{2}(t+1/t)$ in $P$ can grow. So the bracket vanishes and $F_\nu$ solves Bessel's equation.

Step 3 (it is $J_\nu$ , with the right scale). Take first $z = x$ real and positive, and substitute $t = 2u/x$ . A positive scaling maps the cut plane to itself and carries $C$ to a loop $C_u$ of the same kind, which by Cauchy's theorem we may take with unit circle. Since $\tfrac{x}{2}(t-1/t) = u - x^{2}/(4u)$ and $t^{-\nu-1}dt = (x/2)^{\nu}u^{-\nu-1}du$ ,

F_\nu(x) = \Bigl(\frac{x}{2}\Bigr)^{\nu}\,\frac{1}{2\pi i}\int_{C_u} e^{u}\,e^{-x^{2}/(4u)}\,u^{-\nu-1}\,du .

Expand the middle factor, $e^{-x^{2}/(4u)} = \sum_{m\ge0}(-1)^{m}(x^{2}/4)^{m}u^{-m}/m!$ . On $C_u$ we have $|u| \ge 1$ , so the $m$ th term is bounded by $(x^{2}/4)^{m}/m!$ , a convergent series of constants; by the Weierstrass $M$ -test the expansion converges uniformly on $C_u$ , and since $e^{u}u^{-\nu-1}$ is absolutely integrable along $C_u$ , the sum and the contour integral may be interchanged (Integrating and differentiating series term by term). Therefore

F_\nu(x) = \Bigl(\frac{x}{2}\Bigr)^{\nu}\sum_{m=0}^{\infty}\frac{(-1)^{m}}{m!}\Bigl(\frac{x^{2}}{4}\Bigr)^{m}\cdot\frac{1}{2\pi i}\int_{C_u}e^{u}u^{-\nu-m-1}\,du .

By the Lemma below each contour integral equals $1/\Gamma(\nu+m+1)$ , so

F_\nu(x) = \sum_{m=0}^{\infty}\frac{(-1)^{m}}{m!\,\Gamma(\nu+m+1)}\Bigl(\frac{x}{2}\Bigr)^{\nu+2m} = J_\nu(x) \qquad (x>0) .

Step 3a (from the ray to the half-plane). The restriction to real $x$ was only a convenience: the scaling $t \mapsto 2u/z$ by a complex $z$ would tilt the two rays of $C$ and move the cut with them. It costs nothing, because the bound of Step 2 also makes $F_\nu$ analytic on $\operatorname{Re}z>0$ , and $J_\nu$ is analytic there as well. Two functions analytic on a connected open set and agreeing on a subset with a limit point in it — here the ray $z>0$ — agree throughout (Analytic continuation), so $F_\nu = J_\nu$ throughout $\operatorname{Re}z>0$ , which is (S1). $\square$

Lemma(Lemma: Hankel's loop integral for

1/\Gamma

)AdaptedAdapted from Arfken, Weber & Harris (2013), §13.1, Eqs. (13.30)-(13.31) and Fig. 13.2, pp. 604-606

For every $s \in \mathbb{C}$ and every loop $C_u$ of the kind described above,

\frac{1}{2\pi i}\int_{C_u} e^{u}\,u^{-s}\,du \;=\; \frac{1}{\Gamma(s)} .

(Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §13.1, Eq. (13.31), p. 606) state the same contour integral with the cut placed on the positive axis, in the form $\int e^{-t}t^{\nu}dt = 2ie^{i\nu\pi}\Gamma(\nu+1)\sin\nu\pi$ ; rotating $t \mapsto -u$ and using the reflection formula turns one into the other. The short direct proof is as follows.)

By Cauchy's theorem the value does not depend on the radius of the circle, so take $\operatorname{Re}s < 1$ and shrink that radius to $r$ . On it $|e^{u}u^{-s}|\,|du| = O(r^{1-\operatorname{Re}s})$ , which tends to $0$ as $r \to 0$ — this is the step that needs $\operatorname{Re}s<1$ . On the lower edge $u = -\sigma$ with $\arg u = -\pi$ we have $u^{-s} = \sigma^{-s}e^{i\pi s}$ and $du = -d\sigma$ with $\sigma$ running from $\infty$ down to $r$ ; on the upper edge $\arg u = +\pi$ , so $u^{-s} = \sigma^{-s}e^{-i\pi s}$ , and $\sigma$ runs from $r$ out to $\infty$ . Adding the two edges,

\int_{C_u}e^{u}u^{-s}\,du \;\longrightarrow\; \bigl(e^{i\pi s} - e^{-i\pi s}\bigr)\int_0^{\infty}e^{-\sigma}\sigma^{-s}\,d\sigma = 2i\sin(\pi s)\,\Gamma(1-s) ,

the Euler integral converging because $\operatorname{Re}(1-s) > 0$ . Dividing by $2\pi i$ and applying the reflection formula $\Gamma(s)\Gamma(1-s) = \pi/\sin\pi s$ (The Gamma function) gives $\sin(\pi s)\Gamma(1-s)/\pi = 1/\Gamma(s)$ .

Finally, with $r > 0$ held fixed the contour integral is an entire function of $s$ (the integrand is entire in $s$ and dominated locally uniformly), and $1/\Gamma$ is entire (The Gamma function). Two entire functions agreeing on the half-plane $\operatorname{Re}s<1$ agree everywhere. $\square$

Proof(Step 4: the exponential and real forms)AdaptedAdapted from Hassani (2013), §16.2, Eqs. (16.17)-(16.18) and Fig. 16.2, p. 501

The circle of $C$ has radius $1$ , by the normalisation fixed in the definition; put $t = e^{w}$ . Then $t - 1/t = 2\sinh w$ and $t^{-\nu-1}dt = e^{-\nu w}dw$ , which is (S2) once the image contour is identified: the lower ray $t = \sigma e^{-i\pi}$ , $\sigma \ge 1$ , becomes $w = x - i\pi$ with $x = \ln\sigma \ge 0$ traversed inwards; the unit circle becomes the segment $w = i\theta$ , $\theta$ from $-\pi$ to $\pi$ ; the upper ray becomes $w = x + i\pi$ traversed outwards. That is the contour $C'$ .

Now evaluate the three legs.

The segment. With $w = i\theta$ , $\sinh(i\theta) = i\sin\theta$ and $dw = i\,d\theta$ , so this leg contributes $\frac{1}{2\pi i}\int_{-\pi}^{\pi}e^{\,i(z\sin\theta - \nu\theta)}\,i\,d\theta = \frac{1}{2\pi}\int_{-\pi}^{\pi}e^{\,i(z\sin\theta-\nu\theta)}\,d\theta .$ Write $A(\theta) = z\sin\theta - \nu\theta$ , so the integrand is $\cos A + i\sin A$ . Replacing $\theta$ by $-\theta$ replaces $A$ by $-A$ ; since $\cos$ is even and $\sin$ odd as entire functions — no reality of $z$ or $\nu$ is needed — the first part is even in $\theta$ and the second odd. The odd part integrates to zero and the even part doubles, leaving $\tfrac1\pi\int_0^{\pi}\cos(\nu\theta - z\sin\theta)\,d\theta$ : the first term of (S3).
The two rays. Since $\cosh(\mp i\pi) = -1$ and $\sinh(\mp i\pi) = 0$ , we get $\sinh(x \mp i\pi) = -\sinh x$ on both rays, while $e^{-\nu w} = e^{-\nu x}e^{\pm i\nu\pi}$ . The lower ray runs inwards ( $x$ from $\infty$ to $0$ ) and the upper outwards, so with $dw = dx$ their sum is $\frac{1}{2\pi i}\Bigl(e^{-i\nu\pi} - e^{\,i\nu\pi}\Bigr)\int_0^{\infty}e^{-\nu x - z\sinh x}\,dx = -\frac{\sin\nu\pi}{\pi}\int_0^{\infty}e^{-\nu x - z\sinh x}\,dx ,$ using $e^{-i\nu\pi}-e^{i\nu\pi} = -2i\sin\nu\pi$ . That is the second term of (S3).

Both integrals converge for $\operatorname{Re}z>0$ , the second because $e^{-z\sinh x}$ decays doubly exponentially. Adding the three legs gives (S3), and at $\nu = n \in \mathbb{Z}$ the factor $\sin n\pi$ kills the second term. $\square$

Corollary(Splitting the contour: Hankel and Neumann functions)AdaptedAdapted from Arfken, Weber & Harris (2013), §14.4, Eqs. (14.88)-(14.96) and Fig. 14.10, pp. 676-677; Hassani (2013), §16.2, pp. 501-502

By Step 1 of the proof, any contour on which $P(z,t) = e^{(z/2)(t-1/t)}t^{-\nu}\bigl[\nu+\tfrac{z}{2}(t+1/t)\bigr]$ vanishes at both ends produces a solution of Bessel's equation. Besides $t = \infty e^{\pm i\pi}$ , the point $t = 0$ approached along the positive real axis is such an end: there $|e^{(z/2)(t-1/t)}| = e^{-\operatorname{Re}(z)/(2t)\,+\,O(t)}$ vanishes faster than any power of $t$ , so it beats both the $t^{-\nu}$ and the $z/(2t)$ inside $P$ . So the Schläfli loop may be severed at the origin into

$C_1$ : from $t = 0^{+}$ through the upper half-plane out to $\infty e^{+i\pi}$ , and
$C_2$ : from $\infty e^{-i\pi}$ through the lower half-plane back to $t = 0^{+}$ ,

each of which carries a solution in its own right. These solutions are the Hankel functions: for $\operatorname{Re}z > 0$ ,

H^{(1)}_\nu(z) = \frac{1}{\pi i}\int_{C_1}e^{(z/2)(t-1/t)}t^{-\nu-1}\,dt , \qquad H^{(2)}_\nu(z) = \frac{1}{\pi i}\int_{C_2}e^{(z/2)(t-1/t)}t^{-\nu-1}\,dt .

Since $C_1$ and $C_2$ reassemble $C$ while the prefactor $1/\pi i$ is twice the $1/2\pi i$ of (S1), the two displayed integrals add to $2J_\nu$ ; equivalently, half their sum is $J_\nu$ , which is the identity $J_\nu = \tfrac12\bigl(H^{(1)}_\nu+H^{(2)}_\nu\bigr)$ of Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §14.4, Eq. (14.92), p. 677). That check pins down only the sum of the two halves, so it does not by itself establish the labelling; what completes it is a second computation, carried out in Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §14.4, Eqs. (14.94)-(14.96), p. 677), in which the substitution $t = e^{\pm i\pi}/s$ is applied to each half separately. It returns each half to its own contour and yields $H^{(1)}_\nu = e^{-i\nu\pi}H^{(1)}_{-\nu}$ , $H^{(2)}_\nu = e^{i\nu\pi}H^{(2)}_{-\nu}$ , hence a second linear combination giving $J_{-\nu}$ ; the two combinations together force the identification above and reproduce $Y_\nu = \tfrac{1}{2i}\bigl(H^{(1)}_\nu - H^{(2)}_\nu\bigr)$ . The matching real form of $Y_\nu$ is then

Y_\nu(z) = \frac{1}{\pi}\int_0^{\pi}\sin\bigl(z\sin\theta - \nu\theta\bigr)\,d\theta \;-\; \frac{1}{\pi}\int_0^{\infty}\Bigl(e^{\nu t} + e^{-\nu t}\cos\nu\pi\Bigr)e^{-z\sinh t}\,dt

for $\operatorname{Re}z>0$ . (Page 501 of Hassani (Hassani (2013), §16.2, p. 501) displays this as four separate integrals, in the uncombined form $Y_\nu = (\cot\nu\pi)J_\nu - (\sin\nu\pi)^{-1}J_{-\nu}$ , together with the substitution of $\pi-\theta$ for $\theta$ in the third of them; carrying out that substitution and collecting terms gives the two integrals displayed here. The half-plane $\operatorname{Re}z>0$ is the one Hassani states for the corresponding Hankel representations on p. 502; Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §14.4, Eqs. (14.90)-(14.91), p. 676) state the contour forms for real $x>0$ .) The split is the practical payoff of the whole construction: the exponent $\tfrac{z}{2}(t-1/t)$ has saddle points at $t = \pm i$ , one in each half-plane, so $C_1$ can be pushed over $t=+i$ and $C_2$ over $t=-i$ and each Hankel function gets its asymptotics from a single saddle (Asymptotic forms of Bessel functions, The method of steepest descents). On the undivided contour $C$ the two saddles interfere, which is why $J_\nu$ has an oscillatory and not a one-term large- $z$ form.

Remark(Where the hypotheses are used)Original

$\operatorname{Re}z > 0$ . Used twice: to make the rays of $C$ converge at all, and to kill the endpoint term $P$ in Step 2. It cannot be weakened for general $\nu$ . At integer order it is not needed, because the loop closes into a circle on which the integrand is entire in $z$ ; this is exactly why Bessel's integral holds for every complex $z$ while (S3) at non-integer order does not.

$\nu \notin \mathbb{Z}$ versus $\nu \in \mathbb{Z}$ . The branch point at $t=0$ is what forces the cut and the open loop, and the second term of (S3) is its only trace. Dropping that term is not a harmless approximation: at $\nu = \tfrac12$ , $z = 2$ the first integral alone is $0.6304911$ , while $J_{1/2}(2) = 0.5130161$ , so the correction $0.1174750$ carries $23\%$ of the answer.

The vanishing endpoint term, not the contour. Step 1 is an identity, not an estimate; it is the vanishing of $P$ at the ends that is the real hypothesis. This is why the same integrand on three different contours gives three different functions ( $J_\nu$ , $H^{(1)}_\nu$ , $H^{(2)}_\nu$ ), and why a contour with one end at $t = 0^{-}$ would give nothing: there $|e^{-z/2t}| = e^{+\operatorname{Re}z/2|t|}$ blows up.

$\operatorname{Re}s<1$ in the Lemma. Needed only to discard the small circle, and removed afterwards by analytic continuation in $s$ . Without that continuation step the Lemma would not cover the terms of Step 3 with large $m$ .

Real $z$ in Step 3. The series computation is done on the positive real axis only, because a complex scaling $t \mapsto 2u/z$ would drag the branch cut with it. Nothing is lost: the bound of Step 2 makes $F_\nu$ analytic on the open half-plane, and the identity theorem carries the result off the axis. Drop the analyticity step and the proof would establish (S1) for real positive $z$ alone.

Method(Getting an answer out of a representation)Original

Match the form to the question. The loop form (S1) is for questions about analyticity, normalisation and continuation in $\nu$ ; the exponential form (S2) is for order shifts and saddle points, because the order enters only through the linear factor $e^{-\nu w}$ ; the real form (S3) is for numerics and for bounds.
For a relation between neighbouring orders, write the shift as multiplication of the integrand of (S2) by an elementary function of $w$ : $J_{\nu-1}+J_{\nu+1}$ pulls out $2\cosh w$ , and $J_{\nu-1}-J_{\nu+1}$ pulls out $2\sinh w$ .
Remove the extra factor with a total derivative. Use $\partial_w e^{z\sinh w-\nu w} = (z\cosh w - \nu)e^{z\sinh w-\nu w}$ to trade $\cosh w$ for a boundary term plus a multiple of $J_\nu$ . Check that the boundary term vanishes: on $C'$ it does whenever $\operatorname{Re}z>0$ .
For a large-argument form, deform onto a steepest-descent path through one saddle of $\tfrac{z}{2}(t-1/t)$ — which first requires splitting $C$ so that only one saddle is in play (The method of steepest descents).
For a bound, put absolute values inside (S3). At integer order this gives $|J_n(x)| \le 1$ for real $x$ immediately, with no series manipulation.
To expand a periodic function in Bessel functions, write its Fourier coefficient as an integral and change variables until the integrand is $\cos(n\theta - x\sin\theta)$ ; the coefficient is then $J_n(x)$ by inspection. The application below does this.
Always re-check the scale. An integral that solves Bessel's equation determines the function only up to a constant. Fix it by letting $z \to 0$ and comparing with the leading term $(z/2)^{\nu}/\Gamma(\nu+1)$ .

Example(A recurrence at general order, and a numerical check at

\nu = \tfrac12

)Original

(a) The recurrence, from (S2). Let $\operatorname{Re}z>0$ and write $G(w) = e^{z\sinh w-\nu w}$ . Adding the representations of $J_{\nu-1}$ and $J_{\nu+1}$ and factoring out $G$ ,

J_{\nu-1}(z) + J_{\nu+1}(z) = \frac{1}{2\pi i}\int_{C'}G(w)\bigl(e^{w}+e^{-w}\bigr)\,dw = \frac{1}{2\pi i}\int_{C'}2G(w)\cosh w\,dw .

Now $G'(w) = (z\cosh w - \nu)G(w)$ , so $z\,G\cosh w = G' + \nu G$ and

\frac{z}{2}\bigl[J_{\nu-1}(z)+J_{\nu+1}(z)\bigr] = \frac{1}{2\pi i}\int_{C'}G'(w)\,dw \;+\; \nu\,\frac{1}{2\pi i}\int_{C'}G(w)\,dw .

The first integral is the difference of $G$ at the two ends $w = \infty \pm i\pi$ of $C'$ . There $\sinh w = -\sinh x$ with $x \to +\infty$ , while $e^{-\nu w} = e^{-\nu x}e^{\mp i\nu\pi}$ at $w = x \pm i\pi$ , so

|G| = e^{\pm\pi\operatorname{Im}\nu}\;e^{-\sinh x\,\operatorname{Re}z\,-\,\operatorname{Re}(\nu) x} \longrightarrow 0 ,

because $\operatorname{Re}z > 0$ and the leading factor $e^{\pm\pi\operatorname{Im}\nu}$ is a constant, independent of $x$ (it is $1$ for real $\nu$ ). The first integral is therefore zero. The second is $\nu J_\nu(z)$ by (S2). Dividing by $z/2$ ,

J_{\nu-1}(z) + J_{\nu+1}(z) = \frac{2\nu}{z}J_\nu(z) ,

the recurrence of Bessel functions of the first kind, obtained here for every $\nu$ at once rather than from a generating function that exists only at integer order.

(b) The numbers. Take $\nu = \tfrac12$ , $z = 2$ .

The target. From the half-integer closed form $J_{1/2}(x) = \sqrt{2/\pi x}\,\sin x$ (Bessel functions of the first kind), $J_{1/2}(2) = \sqrt{1/\pi}\,\sin 2 = 0.5641896 \times 0.9092974 = 0.5130161 .$
First term of (S3). $\displaystyle \frac1\pi\int_0^{\pi}\cos\Bigl(\frac{\theta}{2}-2\sin\theta\Bigr)d\theta = 0.6304911$ .
Second term. Here $\sin\nu\pi = \sin(\pi/2) = 1$ and $\displaystyle\int_0^{\infty}e^{-t/2-2\sinh t}\,dt = 0.3690585$ , so the term to subtract is $0.3690585/\pi = 0.1174750$ .
Combine. $0.6304911 - 0.1174750 = 0.5130161$ , agreeing with step 1 to all seven digits.
Cross-check of (a) at the same point. The neighbouring half-integer orders are $J_{-1/2}(x) = \sqrt{2/\pi x}\,\cos x$ and $J_{3/2}(x) = \sqrt{2/\pi x}\,(x^{-1}\sin x - \cos x)$ , so $J_{-1/2}(2) = \sqrt{1/\pi}\cos 2 = -0.2347857$ and $J_{3/2}(2) = \sqrt{1/\pi}\bigl(\tfrac12\sin 2 - \cos 2\bigr) = 0.4912938$ . Their sum is $0.2565081$ , and $(2\nu/z)J_\nu(z) = \tfrac12 \times 0.5130161 = 0.2565081$ .

Application(Kepler's equation and the Bessel series for the eccentric anomaly)Original

A body on an elliptical orbit of eccentricity $e \in [0,1)$ has its position fixed by the eccentric anomaly $E$ , which is related to the uniformly growing mean anomaly $M$ (proportional to elapsed time) by Kepler's equation

M = E - e\sin E .

This is transcendental in $E$ , and the classical way to invert it is with Bessel's integral.

The expansion. Replacing $M$ by $-M$ sends $E$ to $-E$ , and replacing $M$ by $M+2\pi$ sends $E$ to $E+2\pi$ , so $E - M$ is an odd, $2\pi$ -periodic, continuously differentiable function of $M$ (differentiable because $dM/dE = 1-e\cos E \ge 1-e > 0$ ). Its Fourier series is therefore a pure sine series converging to it (Fourier series, Pointwise convergence of Fourier series), $E - M = \sum_{n\ge1}b_n\sin nM$ , with

b_n = \frac{2}{\pi}\int_0^{\pi}(E-M)\sin nM\,dM .

Integrate by parts. The boundary term $\bigl[-(E-M)\cos nM / n\bigr]_0^{\pi}$ vanishes, because $E = M$ at both $M=0$ and $M=\pi$ , and $\int_0^{\pi}\cos nM\,dM = 0$ for $n \ge 1$ , so

b_n = \frac{2}{n\pi}\int_0^{\pi}\cos(nM)\,\frac{dE}{dM}\,dM = \frac{2}{n\pi}\int_0^{\pi}\cos\bigl(nE - ne\sin E\bigr)\,dE ,

where the second equality just changes the variable of integration from $M$ to $E$ using Kepler's equation, the limits $M = 0,\pi$ corresponding to $E = 0,\pi$ . The last integral is Bessel's integral at order $n$ and argument $ne$ . Hence

E = M + 2\sum_{n=1}^{\infty}\frac{J_n(ne)}{n}\,\sin nM .

Each Fourier coefficient of the solution of a transcendental equation is a Bessel function evaluated at a point that moves with the index — a fact that no series expansion of $J_n$ makes visible, and that follows in three lines from the integral representation.

A worked value. Take $e = 0.2$ and $M = 1$ radian. The first terms $2J_n(0.2n)\sin(n)/n$ are

$n$	$J_n(0.2n)$	term	partial sum for $E$
1	$0.0995008$	$+0.1674541$	$1.1674541$
2	$0.0197347$	$+0.0179447$	$1.1853988$
3	$0.0043997$	$+0.0004139$	$1.1858127$
4	$0.0010330$	$-0.0003909$	$1.1854218$
5	$0.0002498$	$-0.0000958$	$1.1853260$

The sum converges to $E = 1.1853242$ , and the check is immediate: $E - 0.2\sin E = 1.1853242 - 0.2\times0.9266210 = 1.0000000 = M$ .

Where the contour forms are used instead. In wave problems the real form is the wrong tool and the split contours of the statement above are the right one. Scattering of a plane wave by a cylinder, radiation from a line source, and guided modes in a cylindrical waveguide all require the outgoing solution at large $kr$ , and the one-saddle contour integrals for $H^{(1)}_\nu$ and $H^{(2)}_\nu$ are what deliver its amplitude and phase (Asymptotic forms of Bessel functions). The same integrals, read as functions of $\nu$ rather than $z$ , are the starting point for the uniform large-order expansions used when the order and the argument are comparable, as for a whispering-gallery mode or a grazing ray.