Hypergeometric and other functions

Hypergeometric functions, elliptic integrals and Airy functions.

In this chapter 4

Tick a unit to add it to your PDF basket.

Hypergeometric and related functions

Hypergeometric and related functions

Definition

The hypergeometric function

Adapted

Add to PDF

The Gauss series ${}_2F_1(a,b;c;x)$ and the three-singular-point equation it solves, with Euler's integral representation, Gauss's closed form at $x=1$ , and the reduction of elementary and classical-polynomial functions to it.

Prerequisites: The Frobenius method, The Gamma function, The Beta function, Power series

Definition(Pochhammer symbols and the Gauss series)AdaptedAdapted from Arfken, Weber & Harris (2013), §18.5, eqs. (18.121) and the definition of $(a)_n$, p. 912; Riley, Hobson & Bence (2006), §18.10, eqs. (18.141)–(18.142) and the footnote on the notation ${}_2F_1$, pp. 629–630; Hassani (2013), §15.4, eq. (15.22) and Box 15.4.1, p. 473

For a complex number $a$ and an integer $n \ge 0$ the Pochhammer symbol (rising factorial) is

(a)_0 = 1, \qquad (a)_n = a(a+1)(a+2)\cdots(a+n-1) .

When $a$ is not a non-positive integer, $(a)_n = \Gamma(a+n)/\Gamma(a)$ (The Gamma function).

Let $a$ , $b$ be arbitrary and let $c$ be neither $0$ nor a negative integer. The Gauss series, or hypergeometric series, is

{}_2F_1(a,b;c;x) \;=\; \sum_{n=0}^{\infty} \frac{(a)_n (b)_n}{(c)_n}\,\frac{x^n}{n!} \;=\; 1 + \frac{ab}{c}\,\frac{x}{1!} + \frac{a(a+1)\,b(b+1)}{c(c+1)}\,\frac{x^2}{2!} + \cdots ,

and its sum is the hypergeometric function. Three features are immediate from the display.

Symmetry. $a$ and $b$ enter identically, so ${}_2F_1(a,b;c;x) = {}_2F_1(b,a;c;x)$ .
Termination. If $a$ or $b$ is $0$ or a negative integer $-m$ , then $(a)_n = 0$ for every $n > m$ and the series is a polynomial of degree at most $m$ . Hassani reaches the same conclusion from the Gamma-quotient form of the coefficients, where the poles of $\Gamma(a+k)$ and of $\Gamma(a)$ cancel while $k \le m$ and the pole of $\Gamma(a)$ alone survives afterwards.
It generalises the geometric series. Taking $b = c$ cancels the two middle factors and leaves $\sum_n (a)_n x^n/n!$ , and with $a = 1$ as well this is $\sum_n x^n$ — which is where the name comes from.

Notation. The subscripts in ${}_2F_1$ count Pochhammer symbols: two upstairs, one downstairs. They exist to distinguish the function from the generalised series ${}_pF_q$ , and in particular from the ${}_1F_1$ of Confluent hypergeometric functions. Where no confusion can arise, Riley, Hobson and Bence and Hassani write $F(a,b,c;x)$ and $F(\alpha,\beta;\gamma;z)$ for the same object; Arfken, Weber and Harris separate the numerator from the denominator parameters with a semicolon, and we follow them. Our variable is written $x$ when it is real and $z$ when the complex plane is in view; the function is the same.

Where the series converges. The ratio of consecutive terms is $(a+n)(b+n)x/[(n+1)(c+n)] \to x$ , so unless the series terminates its radius of convergence is $1$ . On the boundary, Riley, Hobson and Bence record for real parameters that the series converges at $x = 1$ when $c > a+b$ and at $x = -1$ when $c > a+b-1$ ; Arfken, Weber and Harris give the same two thresholds, and describe the series as divergent when $c \le a+b-1$ , which is a statement about the endpoints only, the radius of convergence being $1$ in every non-terminating case.

Definition(The hypergeometric equation, its solutions, and Gauss's summation formula)AdaptedAdapted from Riley, Hobson & Bence (2006), §18.10, eqs. (18.136)–(18.143), pp. 628–630, and eq. (18.144) with the value at $x=1$, pp. 631–632; Arfken, Weber & Harris (2013), §18.5, eqs. (18.120), (18.123), pp. 911–912, and Exercises 18.5.7–18.5.8, pp. 915–916; Hassani (2013), §15.4, eqs. (15.21)–(15.24), pp. 473–474; Ahlfors (1979), §4.4, eq. (24) and the following page, pp. 316–318; Coddington (1989), Ch. 4, §9, p. 182

Let $a$ , $b$ , $c$ be given numbers. The hypergeometric equation (Gauss's equation) is

x(1-x)\,y'' + \bigl[\,c - (a+b+1)x\,\bigr]y' - ab\,y = 0 . \tag{H}

Singular points. The points $x = 0$ , $x = 1$ and $x = \infty$ are regular singular points of (H), and there are no others. Conversely, every second-order linear equation whose only singularities in the extended plane are three regular singular points can be brought to the form (H) by a Möbius change of the independent variable followed by multiplication of the unknown by a power of $x$ and a power of $x-1$ .
Solutions at the origin. The exponents of (H) at $x = 0$ are $0$ and $1-c$ . If $c$ is neither $0$ nor a negative integer then $y_1(x) = {}_2F_1(a,b;c;x)$ solves (H) on $\lvert x\rvert < 1$ . If in addition $c$ is not a positive integer, then $y_2(x) = x^{1-c}\,{}_2F_1(a-c+1,\,b-c+1;\,2-c;\,x)$ is a second solution, and $y_1$ , $y_2$ are independent, so the general solution on $0 < \lvert x\rvert < 1$ is $Ay_1 + By_2$ . The solutions regular at the origin are the multiples of $y_1$ . When $c$ is an integer the two formulas coincide or one of them breaks down, and a logarithm appears in the second solution.
Euler's integral representation. If $c > b > 0$ and $\lvert x \rvert < 1$ , ${}_2F_1(a,b;c;x) = \frac{\Gamma(c)}{\Gamma(b)\,\Gamma(c-b)}\int_{0}^{1} t^{\,b-1}(1-t)^{\,c-b-1}(1-tx)^{-a}\,\mathrm{d}t . \tag{E}$
Gauss's summation formula. If $a$ , $b$ , $c$ are real, $c$ is not $0$ or a negative integer, and $c > a+b$ , then the series converges at $x = 1$ and ${}_2F_1(a,b;c;1) = \frac{\Gamma(c)\,\Gamma(c-a-b)}{\Gamma(c-a)\,\Gamma(c-b)} . \tag{G}$

Idea of the proofOriginal

Everything here is forced by the shape of (H). Dividing by $x(1-x)$ shows that the coefficient of $y'$ has simple poles only at $0$ and $1$ and the coefficient of $y$ a simple pole at each as well, which is exactly the Frobenius condition at both points; counting degrees at infinity gives the third singularity and no more. Because the singularities are at $0$ and $1$ and nowhere else in the finite plane, a power series about the origin can only fail at distance $1$ , so the radius of convergence is $1$ without any estimate being needed.

For the series itself, Frobenius turns (H) into a two-term recurrence — the price of the particular quadratic $x(1-x)$ multiplying $y''$ — and a two-term recurrence can simply be iterated. The coefficient that comes out, $(a)_n(b)_n/[(c)_n n!]$ , is the definition of ${}_2F_1$ ; the numerator parameters are those of the recurrence's numerator and the denominator parameter is $c$ , which is why the three parameters of the equation reappear unchanged in the function.

Gauss's formula (G) is the only statement here that is not bookkeeping, and it is Euler's integral in disguise: put $x = 1$ in (E) and the integrand collapses to $t^{b-1}(1-t)^{c-a-b-1}$ , a Beta integral. All the Gamma functions in (G) come from that one Beta function together with the normalising constant already present in (E).

ProofAdaptedAdapted from Riley, Hobson & Bence (2006), §18.10, eqs. (18.137)–(18.143), pp. 629–630, and the derivation of (18.144) and of the value at $x=1$, pp. 631–632; Hassani (2013), §15.4, the recurrence and eqs. (15.22)–(15.24), pp. 473–474; Ahlfors (1979), §4.4, pp. 316–318

1. The singular points. Write (H) as $y'' + p\,y' + q\,y = 0$ with

p(x) = \frac{c-(a+b+1)x}{x(1-x)}, \qquad q(x) = \frac{-ab}{x(1-x)} .

Both are analytic except at $0$ and $1$ . At $x=0$ , $xp(x)$ and $x^2q(x)$ are analytic (the second vanishes there), so $0$ is a regular singular point in the sense of The Frobenius method; the same computation at $x=1$ with $(x-1)p$ and $(x-1)^2q$ gives the second. For infinity substitute $x = 1/t$ : the transformed equation has $t=0$ as a regular singular point because $p(1/t)$ is $O(t)$ and $q(1/t)$ is $O(t^2)$ as $t \to 0$ , which is precisely the criterion. Coddington lists (H) as the example of an equation with these three regular singular points and no others.

For the converse we follow Ahlfors. Place the three singularities at $0$ , $1$ , $\infty$ by a Möbius map of the independent variable, which preserves both the order of the equation and the character of each singularity. Requiring regularity at $0$ and $1$ forces

p(x) = \frac{A}{x} + \frac{B}{x-1}, \qquad q(x) = \frac{C}{x^2} + \frac{D}{x(x-1)} + \frac{E}{(x-1)^2}

with no entire parts, since any polynomial added would make $\infty$ irregular, and regularity at infinity also forces one linear relation among the constants. The indicial equations at $0$ , $1$ and $\infty$ then express $A,\dots,E$ in terms of the six exponents, whose sum is $1$ . Finally, substituting $y = x^{\alpha_1}(x-1)^{\beta_1}w$ , where $\alpha_1$ is one exponent at $0$ and $\beta_1$ one exponent at $1$ , lowers those two exponents to $0$ and leaves an equation of exactly the form (H), with

c = 1+\alpha_1-\alpha_2, \qquad a = \alpha_1+\beta_1+\gamma_1, \qquad b = \alpha_1+\beta_1+\gamma_2 .

Ahlfors assumes here that no exponent difference is an integer, so that the two exponents at each singularity really are distinguishable and the substitution does not collapse.

2. The series. Put $y = \sum_{n\ge0} a_n x^{n+\sigma}$ into (H) and divide by $x^{\sigma-1}$ :

\sum_{n\ge0}\Bigl\{(1-x)(n+\sigma)(n+\sigma-1) + \bigl[c-(a+b+1)x\bigr](n+\sigma) - abx\Bigr\}a_n x^{n} = 0 .

Setting $x=0$ keeps only $n=0$ and gives the indicial equation $\sigma(\sigma-1)+c\sigma = 0$ , with roots $\sigma = 0$ and $\sigma = 1-c$ , as Hassani also records.

Take $\sigma = 0$ and collect the coefficient of $x^{n+1}$ . The terms with $(1-x)$ and with $c$ contribute from $a_{n+1}$ , the rest from $a_n$ :

(n+1)(n+c)\,a_{n+1} = \bigl[n(n-1) + (a+b+1)n + ab\bigr]a_n = (n+a)(n+b)\,a_n ,

the factorisation on the right being the identity $n^2+(a+b)n+ab = (n+a)(n+b)$ . Hence

a_{n+1} = \frac{(a+n)(b+n)}{(n+1)(c+n)}\,a_n ,

which is solvable in closed form because only two consecutive coefficients occur. With $a_0=1$ , iterating gives $a_n = (a)_n(b)_n/[(c)_n\,n!]$ , and no denominator vanishes precisely because $c \notin \{0,-1,-2,\dots\}$ . This is the Gauss series, so $y_1 = {}_2F_1(a,b;c;x)$ solves (H) wherever the series converges.

For the convergence radius: by the ratio test it is $1$ unless the series terminates. Ahlfors gives a reason that uses no estimate at all. The series continues analytically along every path avoiding $x=1$ ; if its radius exceeded $1$ it would be entire, and near infinity an entire solution would have to be a combination of the two solutions $x^{-a}g_1$ , $x^{-b}g_2$ that the theory provides there, which is single-valued only if $a$ or $b$ is an integer. Liouville's theorem then rules out a positive integer, leaving only the terminating case.

Second solution. Put $y = x^{1-c}u$ in (H). A direct substitution shows $u$ satisfies the hypergeometric equation with parameters

a_1 = a-c+1, \qquad b_1 = b-c+1, \qquad c_1 = 2-c ,

which is the same bookkeeping Riley, Hobson and Bence perform by taking $\sigma = 1-c$ in the Frobenius sum and comparing the resulting recurrence with the first one. Provided $c_1 = 2-c$ is not $0$ or a negative integer — that is, $c$ is not $2,3,4,\dots$ — the first part applies to $u$ , giving $y_2 = x^{1-c}{}_2F_1(a_1,b_1;c_1;x)$ . Independence is free: $y_1(0) = 1$ while $y_2$ behaves like $x^{1-c}$ at the origin, so for non-integer $c$ no constant multiple of one is the other.

3. Euler's integral. Assume $c>b>0$ and $\lvert x\rvert<1$ . Start from the series and write each coefficient as a Beta function (The Beta function): since $(b)_n/(c)_n = \Gamma(b+n)\Gamma(c)/[\Gamma(b)\Gamma(c+n)]$ and $B(b+n,c-b) = \Gamma(b+n)\Gamma(c-b)/\Gamma(c+n)$ ,

\frac{(b)_n}{(c)_n} = \frac{\Gamma(c)}{\Gamma(b)\Gamma(c-b)}\,B(b+n,\,c-b) = \frac{\Gamma(c)}{\Gamma(b)\Gamma(c-b)}\int_0^1 t^{\,b+n-1}(1-t)^{\,c-b-1}\,\mathrm{d}t ,

the Beta integral converging because $b+n>0$ and $c-b>0$ . Substituting this into the series,

{}_2F_1(a,b;c;x) = \frac{\Gamma(c)}{\Gamma(b)\Gamma(c-b)}\int_0^1 t^{\,b-1}(1-t)^{\,c-b-1}\sum_{n\ge0}\frac{(a)_n}{n!}(tx)^n\,\mathrm{d}t .

Interchanging sum and integral is legitimate: for $0\le t\le1$ and fixed $\lvert x\rvert = r < 1$ the inner series is dominated termwise by $\lvert(a)_n\rvert r^n/n!$ , whose sum is finite, so the convergence is uniform in $t$ and the integrable factor $t^{b-1}(1-t)^{c-b-1}$ may be carried through. Finally the inner sum is the binomial series $\sum_n (a)_n w^n/n! = (1-w)^{-a}$ at $w=tx$ , which gives (E). Riley, Hobson and Bence run this computation in the same direction; Arfken, Weber and Harris set it as Exercise 18.5.7 with the same hypothesis $c>b>0$ .

4. Gauss's formula. Suppose first $c>b>0$ and $c-a-b>0$ . Let $x \to 1^-$ in (E). The integrand increases to $t^{b-1}(1-t)^{c-a-b-1}$ , which is integrable on $(0,1)$ exactly because $b>0$ and $c-a-b>0$ ; monotone convergence (for $a>0$ ; for $a\le0$ the integrand is bounded by an integrable function uniformly in $x$ , and dominated convergence applies) therefore permits the limit under the integral sign, and

\lim_{x\to1^-}{}_2F_1(a,b;c;x) = \frac{\Gamma(c)}{\Gamma(b)\Gamma(c-b)}\int_0^1 t^{\,b-1}(1-t)^{\,c-a-b-1}\mathrm{d}t = \frac{\Gamma(c)}{\Gamma(b)\Gamma(c-b)}\cdot\frac{\Gamma(b)\Gamma(c-a-b)}{\Gamma(c-a)} ,

which is the right-hand side of (G). Since $c>a+b$ , the series converges at $x=1$ , and by Abel's theorem (Power series) its sum there is that same limit. The restriction $c>b$ is removable by the symmetry in $a$ and $b$ whenever $c>a$ , and the remaining cases follow by analytic continuation in the parameters, both sides of (G) being analytic where $c>a+b$ holds. $\square$

Remark(Where the hypotheses are used)Original

$c \notin \{0,-1,-2,\dots\}$ is used exactly once, at the step where the recurrence is iterated: a vanishing $(c)_n$ would make the coefficient undefined. If $c$ is a non-positive integer the exponent $0$ is still a root of the indicial equation, but the Frobenius recurrence breaks down at $n = -c$ and the series solution must be replaced by the other exponent's solution together with a logarithm.
$c$ not a positive integer is what makes $y_2$ exist and differ from $y_1$ . For $c=1$ the two formulas are literally the same function. For $c = 2,3,4,\dots$ the parameter $2-c$ of $y_2$ is a non-positive integer, so $y_2$ is not defined at all. Arfken, Weber and Harris note the escape hatch: for integer $c$ one of the two may still be rescued if $a$ or $b$ is an integer too, so that the offending denominator is cancelled by a vanishing numerator.
$c > b > 0$ is used twice in (E), and for different reasons. $b>0$ makes $t^{b-1}$ integrable at $t=0$ ; $c-b>0$ makes $(1-t)^{c-b-1}$ integrable at $t=1$ . Neither is needed for the series, so (E) is strictly weaker than the series as a definition — but it is the representation that continues ${}_2F_1$ analytically past $x=1$ into the plane cut along $[1,\infty)$ , which the series cannot do.
$c > a+b$ is needed twice as well: once for the convergence of the series at $x=1$ , and once for the convergence of the Beta integral at $t=1$ after $x$ has been set to $1$ . The two conditions are the same inequality, which is not a coincidence: both measure the same failure. If $c \le a+b$ the limit is $+\infty$ for positive parameters, in line with the pole that the right-hand side of (G) has at $c-a-b = 0$ . A concrete case: $a=b=\tfrac12$ , $c=1$ gives $c-a-b=0$ and ${}_2F_1(\tfrac12,\tfrac12;1;x) = 2K(x)/\pi$ diverges logarithmically as $x\to1$ (Elliptic integrals).
No exponent difference an integer is Ahlfors's hypothesis in the converse part of statement 1, not ours in the rest. It is what lets him pick one exponent at $0$ and one at $1$ and strip them off; if two exponents at a point coincide, the second solution there is logarithmic and the normal form is reached only after a limit.

Method(Recognising and reducing a hypergeometric equation)AdaptedAdapted from Riley, Hobson & Bence (2006), §18.10 and §18.10.1, the canonical-form discussion and the Chebyshev worked example, pp. 628–631; Arfken, Weber & Harris (2013), §18.5, the alternate forms (18.124)–(18.125) and the representations (18.128)–(18.135), pp. 913–914; Hassani (2013), §15.4, the Jacobi, Gegenbauer and Legendre reductions, eqs. (15.32)–(15.35), pp. 477–478; Coddington (1989), Ch. 4, §9, Exercise 5, p. 182

Given a second-order linear equation and asked whether its solutions are hypergeometric:

Count the singular points in the extended plane. Find every point where the equation, written as $y''+py'+qy=0$ , has $p$ or $q$ singular, and test infinity by $x=1/t$ . If there are exactly three and all are regular, the equation is hypergeometric up to the changes below; if there are more, or one is irregular, it is not (two singular points give only powers and logarithms, and an irregular point at infinity leads instead to Confluent hypergeometric functions).
Move the singular points to $0$ , $1$ , $\infty$ with the Möbius map determined by where they are. For an equation singular at $\pm1$ and $\infty$ — the usual case in physics, because that is where Legendre, Gegenbauer and Chebyshev live — the map is $x = (1-z)/2$ , which sends $z=1$ to $x=0$ and $z=-1$ to $x=1$ . Remember $\mathrm{d}/\mathrm{d}z = -\tfrac12\,\mathrm{d}/\mathrm{d}x$ .
Strip the exponents. Compute the indicial roots at $x=0$ and at $x=1$ ; if the smaller ones are not already $0$ , substitute $y = x^{\alpha}(1-x)^{\beta}w$ with those values, which shifts them to zero and leaves an equation in the canonical form.
Read off $a$ , $b$ , $c$ by matching coefficients with (H): $c$ is the constant term in the bracket, $a+b+1$ is minus the coefficient of $x$ in it, and $ab$ is minus the coefficient of $y$ . So $a$ and $b$ are the two roots of $\lambda^2 - (a+b)\lambda + ab = 0$ with these numbers, and which root is called $a$ does not matter.
Fix the normalisation. The solution you want is a multiple of ${}_2F_1$ ; find the multiple by comparing both sides at one convenient point, usually $x=0$ , where ${}_2F_1 = 1$ . Riley, Hobson and Bence do exactly this for the Chebyshev polynomials, noting that the two normalisations happen to agree there, and Arfken, Weber and Harris list the resulting prefactors for the ultraspherical, Legendre and Chebyshev families.
If you need the solution near $x=1$ or $x=\infty$ instead, do not re-derive it: Hassani's bases ${}_2F_1(a,b;a+b-c+1;1-x), \qquad x^{-a}\,{}_2F_1\!\left(a,\,1+a-c;\,a-b+1;\,\tfrac1x\right)$ and their partners are obtained from the one at the origin by the substitutions $x\mapsto1-x$ and $x\mapsto1/x$ . There are twenty-four such expressions in all — Kummer's solutions — of which only two are independent at any one point.

Shortcuts worth knowing. Many identifications are quicker from the series than from the equation: $(1-x)^{-a} = {}_2F_1(a,b;b;x)$ , $x^{-1}\ln(1+x) = {}_2F_1(1,1;2;-x)$ , $x^{-1}\arcsin x = {}_2F_1(\tfrac12,\tfrac12;\tfrac32;x^2)$ and $x^{-1}\arctan x = {}_2F_1(\tfrac12,1;\tfrac32;-x^2)$ are all read off by comparing coefficients. And the derivative relation

\frac{\mathrm{d}}{\mathrm{d}x}{}_2F_1(a,b;c;x) = \frac{ab}{c}\,{}_2F_1(a+1,b+1;c+1;x)

means that differentiating never leaves the family.

Example(Legendre's equation reduced, and a number from Gauss's formula)Original

Part 1: the reduction. Legendre's equation of degree $\lambda$ is

(1-z^2)\,u'' - 2z\,u' + \lambda(\lambda+1)u = 0 .

Its singular points are $z=\pm1$ and $\infty$ , all regular, so step 1 of the method succeeds. Apply step 2 with $z = 1-2x$ , so that $x = (1-z)/2$ and $\mathrm{d}/\mathrm{d}z = -\tfrac12\,\mathrm{d}/\mathrm{d}x$ . Then $1-z^2 = (1-z)(1+z) = 2x\cdot2(1-x) = 4x(1-x)$ and, writing $y(x) = u(z)$ ,

4x(1-x)\cdot\tfrac14 y'' - 2(1-2x)\cdot\bigl(-\tfrac12 y'\bigr) + \lambda(\lambda+1)y = 0 ,

that is

x(1-x)y'' + (1-2x)y' + \lambda(\lambda+1)y = 0 .

This is already in the form (H) — step 3 is vacuous, because the exponents at $x=0$ are $0$ and $1-c=0$ — and step 4 reads off

c = 1, \qquad a+b+1 = 2, \qquad ab = -\lambda(\lambda+1) .

So $a+b=1$ and $ab=-\lambda(\lambda+1)$ , whose roots are $a = -\lambda$ and $b=\lambda+1$ : indeed $a+b = 1$ and $ab = -\lambda(\lambda+1)$ . Step 5: at $x=0$ , i.e. $z=1$ , the Legendre polynomial satisfies $P_n(1)=1$ and ${}_2F_1=1$ , so the normalising constant is $1$ and

P_{\lambda}(z) = {}_2F_1\!\left(-\lambda,\;\lambda+1;\;1;\;\frac{1-z}{2}\right) .

(Coddington sets the same reduction as an exercise, in the form $\psi(t) = \varphi(2t-1)$ , and records the parameters $\alpha = p+1$ , $\beta = -p$ , $\gamma = 1$ — the same pair $\{-\lambda,\lambda+1\}$ with the labels exchanged, which the symmetry of ${}_2F_1$ makes irrelevant.)

A checkable number. Take $\lambda = 3$ and $z = \tfrac12$ , so the argument is $x = (1-\tfrac12)/2 = \tfrac14$ . Since $a=-3$ the series terminates after four terms:

{}_2F_1\!\left(-3,4;1;\tfrac14\right) = \sum_{n=0}^{3}\frac{(-3)_n(4)_n}{(1)_n}\frac{(1/4)^n}{n!} .

Term by term, with $(-3)_1=-3$ , $(-3)_2=6$ , $(-3)_3=-6$ and $(4)_1=4$ , $(4)_2=20$ , $(4)_3=120$ , $(1)_n = n!$ :

n=0:\;1, \qquad n=1:\;\frac{(-3)(4)}{1}\cdot\frac{1/4}{1} = -3, \qquad n=2:\;\frac{(6)(20)}{2}\cdot\frac{1/16}{2} = \frac{30}{16} = 1.875,

n=3:\;\frac{(-6)(120)}{6}\cdot\frac{1/64}{6} = -120\cdot\frac{1}{384} = -0.3125 .

The sum is $1-3+1.875-0.3125 = -0.4375$ . Against this, $P_3(z) = \tfrac12(5z^3-3z)$ gives

P_3\!\left(\tfrac12\right) = \tfrac12\left(\tfrac58-\tfrac32\right) = \tfrac12\left(-\tfrac78\right) = -\frac{7}{16} = -0.4375 ,

in exact agreement (Legendre polynomials).

Part 2: a value at $x=1$ . Take $a=\tfrac12$ , $b=\tfrac13$ , $c=2$ . Then $c-a-b = 2-\tfrac56 = \tfrac76 > 0$ , so (G) applies:

{}_2F_1\!\left(\tfrac12,\tfrac13;2;1\right) = \frac{\Gamma(2)\,\Gamma(7/6)}{\Gamma(3/2)\,\Gamma(5/3)} = \frac{1 \times 0.9277\ldots}{0.8862\ldots \times 0.9027\ldots} = 1.159595\ldots

The series itself converges to the same value, but slowly, as it must: its terms decay only like $n^{-(1+(c-a-b))} = n^{-13/6}$ . Summing the first $60$ terms gives $1.158077$ , and the first $2000$ give $1.159570$ — creeping towards $1.1595953$ from below, which is the practical argument for knowing (G) rather than summing. All the numbers in this example were checked in exact or high-precision arithmetic.

ApplicationOriginal

One function for the whole classical zoo. Nearly every named function that solves a second-order equation with three regular singular points is a ${}_2F_1$ with particular parameters: Legendre and associated Legendre functions, Gegenbauer (ultraspherical) and Jacobi polynomials, Chebyshev polynomials of all four kinds, and the Legendre functions of the second kind. The practical payoff is not elegance but transfer. A single library routine for ${}_2F_1$ evaluates all of them; a single contiguous relation among ${}_2F_1(a\pm1,b,c;x)$ specialises to the three-term recurrences that each family is usually introduced with; the derivative relation above specialises to each family's derivative formula; and Gauss's formula (G) supplies normalisation constants that would otherwise each need their own integral. Arfken, Weber and Harris make the point in reverse, by listing the prefactors needed to put each classical polynomial on its conventional scale — the mathematics is one function, the conventions are many.

Exactly solvable quantum mechanics. A one-dimensional Schrödinger problem is solvable in closed form essentially when the substitution that removes the asymptotic behaviour leaves a hypergeometric equation. The Pöschl–Teller well $V(x) = -V_0\operatorname{sech}^2(x/L)$ is the standard case: with $\xi = \tanh^2(x/L)$ or $\xi = \tfrac12(1+\tanh(x/L))$ the radial equation becomes (H), the bound states are the terminating cases where a numerator parameter is a negative integer — exactly the termination clause of the definition — and the reflection and transmission amplitudes above the barrier are ratios of Gamma functions, which is Gauss's formula (G) at work, since the connection between the solution that is purely outgoing at $x=+\infty$ and the incoming-plus-reflected solution at $x=-\infty$ is a linear relation among the hypergeometric bases at $0$ , $1$ and $\infty$ . The famous reflectionless potentials, where $T=1$ at every energy, are the parameter values at which one Gamma function in the denominator hits a pole.

Angular problems in three dimensions. Because $P_\lambda(\cos\theta) = {}_2F_1(-\lambda,\lambda+1;1;\sin^2(\theta/2))$ by the reduction worked above, the hypergeometric form is what makes non-integer degree $\lambda$ meaningful. That is needed whenever a boundary is a cone rather than a full sphere — the field near the tip of a conical conductor, or near the vertex of a wedge — where the separation constant is fixed by requiring $P_\lambda(\cos\theta_0)=0$ at the cone's half-angle and the admissible $\lambda$ come out non-integral. There is no Rodrigues formula or finite polynomial to fall back on in that situation; the series definition is the definition.

Numerical practice. The convergence behaviour in Part 2 of the example is the reason production code never evaluates ${}_2F_1$ from its series alone. Near $x=1$ one applies a transformation such as ${}_2F_1(a,b;c;x) = (1-x)^{c-a-b}{}_2F_1(c-a,c-b;c;x)$ or one of the connection formulas to move the argument into a region where the series converges geometrically; for $\lvert x\rvert>1$ one uses the basis at infinity. The twenty-four Kummer solutions are, from this point of view, a numerical toolkit rather than a curiosity.

Definition

Confluent hypergeometric functions

Adapted

Add to PDF

Kummer's equation, its entire solution $M(a,c,x)$ and the second solution $U(a,c,x)$ , obtained by merging two singular points of the hypergeometric equation, together with the reduction of Hermite, Laguerre, Bessel, error and incomplete-Gamma functions to them.

Prerequisites: The hypergeometric function, The Gamma function, The Frobenius method

Definition(Kummer's equation and the functions

M

and

U

)AdaptedAdapted from Arfken, Weber & Harris (2013), §18.6, eqs. (18.136)–(18.138), (18.141)–(18.142), p. 917; Riley, Hobson & Bence (2006), §18.11, eqs. (18.147)–(18.149) and the definition of $U$, pp. 633–634; Hassani (2013), §15.5, eq. (15.37), the series (15.38) and Box 15.5.1, p. 479

Kummer's equation, also called the confluent hypergeometric equation, is

x\,y'' + (c-x)\,y' - a\,y = 0 . \tag{K}

It has a regular singular point at $x=0$ and an irregular (essential) one at $x=\infty$ ; there are no others.

For $c$ neither $0$ nor a negative integer, Kummer's function is

M(a,c,x) \;=\; {}_1F_1(a;c;x) \;=\; \sum_{n=0}^{\infty}\frac{(a)_n}{(c)_n}\,\frac{x^n}{n!} \;=\; 1+\frac{a}{c}\frac{x}{1!}+\frac{a(a+1)}{c(c+1)}\frac{x^2}{2!}+\cdots

with $(a)_n$ the Pochhammer symbol of The hypergeometric function. The ratio of consecutive terms is $(a+n)x/[(c+n)(n+1)] \to 0$ , so the series converges for every finite $x$ (and every finite complex $z$ ): $M$ is an entire function of its argument, which is what the loss of the finite singularity at $x=1$ buys. As with ${}_2F_1$ , if $a$ is $0$ or a negative integer $-m$ the series terminates and $M$ is a polynomial of degree $m$ .

Arfken, Weber and Harris note that the commas in $M(a,c,x)$ and the semicolons in ${}_1F_1(a;c;x)$ are both standard and denote the same function; we use $M$ . Riley, Hobson and Bence write $M(a,c;x)$ and Hassani writes $\Phi(a;\gamma;z)$ .

The second solution is less tidy. For $c$ not an integer, $x^{1-c}M(a-c+1,2-c,x)$ solves (K) and is independent of $M$ ; but it is singular as $c \to 2,3,4,\dots$ while $M$ itself is singular as $c \to 0,-1,-2,\dots$ , so neither formula is usable across integer $c$ . The combination that is,

U(a,c,x) \;=\; \frac{\pi}{\sin \pi c}\left[\frac{M(a,c,x)}{\Gamma(a-c+1)\Gamma(c)} \;-\; x^{1-c}\,\frac{M(a-c+1,2-c,x)}{\Gamma(a)\Gamma(2-c)}\right] ,

is Tricomi's function. The two divergences cancel and $U$ has a finite limit as $c$ approaches an integer, which is exactly the device used for the Neumann function $Y_\nu$ in Neumann and Hankel functions; Arfken, Weber and Harris point out the parallel explicitly, including that the formula as written is indeterminate at integer $c$ and must be read as a limit.

A conflict of sign to be aware of. Riley, Hobson and Bence print the same combination with $\Gamma(a)\Gamma(2-c)$ in the second denominator, Arfken, Weber and Harris print $\Gamma(a)\Gamma(-c)$ . The two differ, since $\Gamma(2-c)\ne\Gamma(-c)$ in general. We follow Riley, Hobson and Bence, whose version is the one that reduces to a finite multiple of the standard $U$ ; the discrepancy is recorded rather than silently resolved.

Definition(Solutions, confluence, and the elementary reductions)AdaptedAdapted from Riley, Hobson & Bence (2006), §18.11, eqs. (18.147)–(18.150), pp. 633–634, and eqs. (18.151)–(18.152), p. 635; Arfken, Weber & Harris (2013), §18.6, eqs. (18.137)–(18.151), pp. 917–919; Hassani (2013), §15.5, eqs. (15.36)–(15.38) and Proposition 15.5.2, p. 479, and eqs. (15.42)–(15.43), pp. 482–483

Solutions at the origin. The exponents of (K) at $x=0$ are $0$ and $1-c$ . If $c\notin\{0,-1,-2,\dots\}$ then $y_1 = M(a,c,x)$ solves (K) on all of $\mathbb{C}$ ; if in addition $c$ is not a positive integer, $y_2 = x^{1-c}M(a-c+1,2-c,x)$ is a second, independent solution, and every solution is a combination of the two. For $c=1$ the two coincide.
Confluence. With $b$ as the parameter sent to infinity, $x(1-x)y''+[c-(a+b+1)x]y'-ab\,y=0 \quad\xrightarrow[\;b\to\infty\;]{\;x = z/b\;}\quad z\,u''+(c-z)u'-a\,u=0 ,$ and correspondingly $M(a,c,z) = \lim_{b\to\infty}{}_2F_1(a,b;c;z/b)$ . The two regular singular points of the hypergeometric equation at $x=1$ and $x=\infty$ merge into the single irregular point at $z=\infty$ ; this merging is what "confluent" names.
Integral representations. $M(a,c,x) = \frac{\Gamma(c)}{\Gamma(a)\Gamma(c-a)}\int_0^1 \mathrm{e}^{xt}\,t^{\,a-1}(1-t)^{\,c-a-1}\,\mathrm{d}t , \qquad c>a>0 , \tag{IM}$ $U(a,c,x) = \frac{1}{\Gamma(a)}\int_0^{\infty}\mathrm{e}^{-xt}\,t^{\,a-1}(1+t)^{\,c-a-1}\,\mathrm{d}t , \qquad \operatorname{Re}x>0,\ a>0 . \tag{IU}$
Kummer's transformation and the derivative. $M(a,c,x) = \mathrm{e}^{x}\,M(c-a,c,-x), \qquad \frac{\mathrm{d}}{\mathrm{d}x}M(a,c,x) = \frac{a}{c}\,M(a+1,c+1,x) .$
Reductions. With $n$ , $m$ non-negative integers and $\nu$ arbitrary, $\mathrm{e}^{x} = M(a,a,x), \qquad L_n(x) = M(-n,1,x), \qquad L_n^{m}(x) = \frac{(n+m)!}{n!\,m!}\,M(-n,m+1,x),$ $H_{2n}(x) = (-1)^n\frac{(2n)!}{n!}\,M\!\left(-n,\tfrac12,x^2\right), \qquad H_{2n+1}(x) = (-1)^n\frac{2\,(2n+1)!}{n!}\,x\,M\!\left(-n,\tfrac32,x^2\right),$ $J_\nu(x) = \frac{\mathrm{e}^{-\mathrm{i}x}}{\Gamma(\nu+1)}\left(\frac{x}{2}\right)^{\nu} M\!\left(\nu+\tfrac12,\,2\nu+1,\,2\mathrm{i}x\right), \qquad I_\nu(x) = \frac{\mathrm{e}^{-x}}{\Gamma(\nu+1)}\left(\frac{x}{2}\right)^{\nu} M\!\left(\nu+\tfrac12,\,2\nu+1,\,2x\right),$ $\operatorname{erf}(x) = \frac{2x}{\sqrt{\pi}}\,M\!\left(\tfrac12,\tfrac32,-x^2\right), \qquad \gamma(a,x) = \int_0^x \mathrm{e}^{-t}t^{\,a-1}\mathrm{d}t = \frac{x^a}{a}\,M(a,a+1,-x) \quad (\operatorname{Re}a>0) .$

Idea of the proofOriginal

Two things make (K) easier than its parent, and both come from the same source: the singularity that used to sit at $x=1$ has been pushed off to infinity and absorbed there. First, there is no finite obstruction left for a power series about the origin to run into, so the series converges everywhere instead of on a disc. Second, the recurrence that the Frobenius method produces loses one factor — it is $(c+n)(n+1)a_{n+1} = (a+n)a_n$ rather than $(c+n)(n+1)a_{n+1}=(a+n)(b+n)a_n$ — which is precisely the statement that one numerator parameter has been spent.

The price is paid at infinity, where the merged singularity is irregular. There is no convergent Frobenius series there, only an asymptotic one, and the two solutions behave quite differently: generically $M$ grows like $\mathrm{e}^{x}x^{a-c}$ while $U$ decays like $x^{-a}$ . This asymmetry is the whole physical content of the function. Almost every quantum bound-state problem that is exactly solvable is solvable because a boundary condition at infinity kills the growing solution, and the only way $M$ can fail to grow is by terminating — which happens exactly when $a$ is a non-positive integer. Quantization conditions of the form " $a = -N$ " are therefore not an accident of hydrogen; they are the generic shape of the answer.

The reductions in part 5 all have the same explanation. Bessel's equation, the Hermite equation and the Laguerre equation each have one regular singular point and one irregular point at infinity, so by the classification each of them must become (K) once the asymptotic factor $z^{\mu}\mathrm{e}^{-\nu z}$ is removed; the only work is finding $\mu$ and $\nu$ .

ProofAdaptedAdapted from Riley, Hobson & Bence (2006), §18.11, the worked confluence and eqs. (18.148)–(18.150), pp. 633–634, and the derivation of (18.150) with the relations on p. 635; Hassani (2013), §15.5, eqs. (15.36)–(15.38), p. 479, and the substitution scheme of §15.5.2, pp. 482–483; Arfken, Weber & Harris (2013), §18.6, eqs. (18.137), (18.141), (18.144)–(18.145), pp. 917–918

1. The series. Insert $y=\sum_{n\ge0}a_nx^{n+\sigma}$ into (K) and divide by $x^{\sigma-1}$ . The $n=0$ term at $x=0$ gives the indicial equation $\sigma(\sigma-1)+c\sigma=0$ with roots $0$ and $1-c$ , the same roots as for the hypergeometric equation because (K) has the same behaviour at the origin. With $\sigma=0$ , the coefficient of $x^{n}$ gives

(n+1)(n+c)\,a_{n+1} = (n+a)\,a_n ,

so with $a_0=1$ , $a_n = (a)_n/[(c)_n n!]$ , which is $M$ . Convergence for all $x$ was checked in the definition, and the series for $y'$ and $y''$ converge there too (differentiating a power series does not change its radius), so the termwise manipulation is valid on the whole plane. For the second solution substitute $y=x^{1-c}u$ : the equation for $u$ is again (K), with $a \mapsto a-c+1$ and $c \mapsto 2-c$ , so $u = M(a-c+1,2-c,x)$ whenever $2-c$ avoids the forbidden set, i.e. $c\ne2,3,4,\dots$ . Independence follows from the behaviour at $0$ : $y_1(0)=1$ while $y_2 \sim x^{1-c}$ , and for non-integer $c$ these cannot be proportional. This is Hassani's Proposition 15.5.2.

2. Confluence. Substitute $x = z/b$ in the hypergeometric equation, so $\mathrm{d}/\mathrm{d}x = b\,\mathrm{d}/\mathrm{d}z$ , and write $u(z)=y(x)$ :

b z\left(1-\frac{z}{b}\right)u'' + \bigl[bc-(a+b+1)z\bigr]u' - ab\,u = 0 .

Divide by $b$ :

z\left(1-\frac{z}{b}\right)u'' + \left[c-\frac{(a+1)z}{b}-z\right]u' - a\,u = 0 ,

and let $b\to\infty$ with $a$ , $c$ , $z$ fixed: the result is $zu''+(c-z)u'-au=0$ , which is (K). Hassani is careful about which parameter goes to infinity and why: the singular points of the intermediate equation sit at $z=0$ , $z=b$ and $z=\infty$ , and it is $z=b$ running out to meet $z=\infty$ that merges them; keeping $c$ finite is forced, since otherwise the coefficient of $u'$ blows up, and $a$ and $b$ may be exchanged, so sending $b$ rather than $a$ is a convention. The same limit applied termwise to the series works because $(b)_n/b^n \to 1$ for each fixed $n$ :

{}_2F_1\!\left(a,b;c;\frac{z}{b}\right) = \sum_{n\ge0}\frac{(a)_n}{(c)_n}\frac{z^n}{n!}\cdot\frac{(b)_n}{b^{n}} \longrightarrow \sum_{n\ge0}\frac{(a)_n}{(c)_n}\frac{z^n}{n!} = M(a,c,z) .

3. The integral representations. For (IM), take Euler's integral (E) of The hypergeometric function in the form obtained by exchanging $a$ and $b$ — legitimate by the symmetry of ${}_2F_1$ —

{}_2F_1(a,b;c;x) = \frac{\Gamma(c)}{\Gamma(a)\Gamma(c-a)}\int_0^1 t^{\,a-1}(1-t)^{\,c-a-1}(1-tx)^{-b}\,\mathrm{d}t ,

valid for $c>a>0$ , and set $x=z/b$ . Then $(1-tz/b)^{-b}\to\mathrm{e}^{tz}$ pointwise in $t$ as $b\to\infty$ , and for $b$ large the factor is bounded on $[0,1]$ uniformly, so the limit passes inside the integral and gives (IM). Riley, Hobson and Bence run exactly this argument; Arfken, Weber and Harris quote (IM) with the same hypothesis $c>a>0$ and offer the alternative route of verifying directly that the integral satisfies (K), excluding the other solution and checking the value at $x=0$ .

For (IU), the integral converges at $t=0$ because $a>0$ and at $t=\infty$ because $\operatorname{Re}x>0$ dominates the algebraic factor; Arfken, Weber and Harris obtain it by the verification route just described.

4. Kummer's transformation. In (IM) substitute $t \mapsto 1-t$ . The measure is unchanged, $t^{a-1}(1-t)^{c-a-1}$ becomes $(1-t)^{a-1}t^{c-a-1}$ , and $\mathrm{e}^{xt}$ becomes $\mathrm{e}^{x}\mathrm{e}^{-xt}$ . Pulling the constant $\mathrm{e}^x$ out,

M(a,c,x) = \mathrm{e}^{x}\,\frac{\Gamma(c)}{\Gamma(a)\Gamma(c-a)}\int_0^1 \mathrm{e}^{-xt}\,t^{\,(c-a)-1}(1-t)^{\,c-(c-a)-1}\,\mathrm{d}t = \mathrm{e}^{x}M(c-a,c,-x) ,

the last step being (IM) itself with $a$ replaced by $c-a$ and $x$ by $-x$ . The identity holds for all $x$ by analytic continuation, both sides being entire, and for all admissible $a$ , $c$ because both sides are analytic in the parameters. The derivative relation is immediate from the series: differentiating termwise and shifting the index,

\frac{\mathrm{d}}{\mathrm{d}x}\sum_{n\ge1}\frac{(a)_n}{(c)_n}\frac{x^n}{n!} = \sum_{n\ge1}\frac{(a)_n}{(c)_n}\frac{x^{n-1}}{(n-1)!} = \frac{a}{c}\sum_{m\ge0}\frac{(a+1)_m}{(c+1)_m}\frac{x^m}{m!} ,

using $(a)_{m+1} = a(a+1)_m$ and $(c)_{m+1}=c(c+1)_m$ .

5. The reductions. Each is the same two-line computation, and we give the general form, following Hassani's scheme. To reduce

w'' + \left(\frac{A}{z}+B\right)w' + \left(\frac{C}{z^2}+\frac{D}{z}+E\right)w = 0

to (K), put $w = z^{\mu}\mathrm{e}^{-\nu z}f(z)$ and choose $\nu$ and $\mu$ to kill the two terms that (K) does not have: $\nu$ is fixed by the requirement that the constant term in the coefficient of $f$ vanish, and $\mu$ by the requirement that the $z^{-2}$ term vanish. What is left is $zf''+(\alpha-\beta z)f'-\delta f=0$ , which after rescaling $z$ is (K). Two instances:

Bessel. For $w''+z^{-1}w'+(1-\nu^2z^{-2})w=0$ , the choices $\mu=\nu$ and the exponential factor $\mathrm{e}^{-\mathrm{i}z}$ give, after setting $t=2\mathrm{i}z$ , the equation $t f''+(2\nu+1-t)f' - (\nu+\tfrac12)f=0$ , which is (K) with $a=\nu+\tfrac12$ , $c=2\nu+1$ . Hence $J_\nu(z)$ is a multiple of $z^{\nu}\mathrm{e}^{-\mathrm{i}z}M(\nu+\tfrac12,2\nu+1,2\mathrm{i}z)$ , and comparing the leading coefficient with the series $J_\nu(z)=(z/2)^\nu/\Gamma(\nu+1)+\cdots$ of Bessel functions of the first kind fixes the multiple to be $2^{-\nu}/\Gamma(\nu+1)$ , giving the formula in part 5. The modified case is the same with $\mathrm{e}^{-z}$ in place of $\mathrm{e}^{-\mathrm{i}z}$ .
Laguerre. The Laguerre equation $xy''+(1-x)y'+ny=0$ is (K) with $a=-n$ , $c=1$ , with no substitution at all; since $L_n(0)=1=M(-n,1,0)$ , the two agree exactly.

The Hermite and error-function cases are the same calculation with $\mathrm{e}^{-x^2}$ -type factors and the change of variable $x\mapsto x^2$ , which Arfken, Weber and Harris record as the alternate form of (K) in the variable $x^2$ . The incomplete Gamma identity follows from (IM) with $c=a+1$ , where the weight $(1-t)^{c-a-1}$ is $1$ and the integral is elementary. $\square$

Remark(Where the hypotheses are used)Original

$c\notin\{0,-1,-2,\dots\}$ is used only in iterating the recurrence, exactly as for ${}_2F_1$ : $(c)_n$ must not vanish. $c$ not a positive integer is what keeps $x^{1-c}M(a-c+1,2-c,x)$ both defined ( $2-c$ must avoid the forbidden set) and distinct from $M$ . Both conditions are conditions on the exponent difference $1-c$ being a non-integer, which is the standard Frobenius proviso of The Frobenius method; when it fails, a logarithm enters, and $U$ is the bookkeeping device that absorbs it.
$c>a>0$ in (IM) plays the same two roles as $c>b>0$ did in Euler's integral: $a>0$ for integrability at $t=0$ , $c-a>0$ for integrability at $t=1$ . Note which parameter is which — the exponent in the first factor is $a-1$ , not $c-1$ — because the confluence is applied to the form of Euler's integral with $a$ and $b$ swapped.
$\operatorname{Re}x>0$ in (IU) is what makes $\mathrm{e}^{-xt}$ beat the algebraic growth $(1+t)^{c-a-1}$ at $t=\infty$ . It is a genuine restriction, not a technicality: $U$ has a branch cut along the negative real axis, and (IU) is the representation of its principal branch.
The confluence limit needs $a$ , $c$ , $z$ held fixed. If $c$ were allowed to grow with $b$ the limiting equation would be different (Hassani's point: a careless limit collapses the equation to the elementary $\ddot w + (\gamma/t)\dot w = 0$ ). Interchanging the limit with the sum in part 2 of the proof is justified termwise for each fixed $n$ , and the resulting series converges everywhere — but this is a statement about formal coefficients made rigorous by the entirety of the limit function, not a uniform statement about ${}_2F_1$ for finite $b$ .
Termination needs $a$ , not $c$ . $M$ is a polynomial exactly when $a \in\{0,-1,-2,\dots\}$ . This is one-sided in a way that matters physically: making $c$ a negative integer does not terminate the series, it destroys it.

Method(From a radial equation to Kummer's equation)AdaptedAdapted from Hassani (2013), §15.5.1, the reduction of the hydrogen-like radial equation, pp. 480–481, and §15.5.2, the Bessel reduction, pp. 482–483; Arfken, Weber & Harris (2013), §18.6, the alternate form (18.143) and the list of representations, pp. 918–919

Given a second-order equation on $(0,\infty)$ with a regular singular point at the origin and an irregular one at infinity — which covers nearly every radial Schrödinger equation:

Clear the first derivative if necessary and bring the equation to the normal form $u'' + \left(\lambda + \frac{A}{r} - \frac{B}{r^{2}}\right)u = 0 ,$ where $u = rR(r)$ for a three-dimensional radial function $R$ and $B$ is typically $\ell(\ell+1)$ .
Rescale so that the constant is $\pm\tfrac14$ . Put $r = kz$ and choose $k$ with $\lambda k^2=-\tfrac14$ for a bound state ( $\lambda<0$ ). The equation becomes $u''+(-\tfrac14+\alpha/z-B/z^2)u=0$ with $\alpha = A k$ .
Peel off both asymptotic behaviours. Write $u(z) = z^{\mu}\mathrm{e}^{-\nu z}f(z)$ and fix the two exponents by the two conditions $\nu^2 = \tfrac14 , \qquad \mu(\mu-1) = B ,$ choosing the signs by the physics: $\nu=+\tfrac12$ so that $u\to0$ at infinity, and $\mu = \ell+1$ (not $-\ell$ ) so that $u$ is finite at the origin.
Read off the parameters. What remains is $z f'' + \bigl(2\mu-z\bigr)f' - (\mu-\alpha)f = 0 ,$ which is (K) with $c = 2\mu$ and $a = \mu-\alpha$ , so $f = M(\mu-\alpha,\,2\mu,\,z)$ for the solution regular at the origin.
Impose the condition at infinity. $\mathrm{e}^{-z/2}M(a,c,z)$ blows up unless the series terminates, because $M(a,c,z)$ grows like $\mathrm{e}^{z}z^{a-c}$ . So demand $a = -N, \qquad N = 0,1,2,\dots ,$ which is the quantization condition. Translating it back through steps 2–4 gives the spectrum, and $f$ is then a polynomial — in the Coulomb case, an associated Laguerre polynomial.

Note on step 5. The growth of $M$ at infinity is the one fact in this recipe that the series does not show directly; it is read from the integral representation (IM), whose integrand is dominated near $t=1$ by $\mathrm{e}^{x}$ . Arfken, Weber and Harris and Riley, Hobson and Bence both state the behaviour without proof, and Hassani argues it by comparing the tail coefficients of $M$ with those of $\mathrm{e}^{z}$ . If instead the physical boundary condition selects decay at infinity without a condition at the origin, the right solution is $U$ , whose own integral representation (IU) makes the decay $x^{-a}$ visible at once.

Example(Two reductions checked to the last digit)Original

(a) A Laguerre polynomial. Take $a=-3$ , $c=1$ . The series terminates after four terms, and with $(1)_n = n!$ ,

M(-3,1,x) = \sum_{n=0}^{3}\frac{(-3)_n}{n!}\frac{x^n}{n!} = 1 + \frac{-3}{1}x + \frac{6}{2}\frac{x^2}{2} + \frac{-6}{6}\frac{x^3}{6} = 1-3x+\frac{3x^2}{2}-\frac{x^3}{6} ,

which is $L_3(x)$ . Check the claim that this solves Kummer's equation with these parameters. With $y=L_3$ , $y'=-3+3x-\tfrac{x^2}{2}$ and $y''=3-x$ ,

x(3-x)+(1-x)\!\left(-3+3x-\tfrac{x^2}{2}\right)+3\!\left(1-3x+\tfrac{3x^2}{2}-\tfrac{x^3}{6}\right)

= \left(3x-x^2\right)+\left(-3+6x-\tfrac{7x^2}{2}+\tfrac{x^3}{2}\right)+\left(3-9x+\tfrac{9x^2}{2}-\tfrac{x^3}{2}\right) = 0 ,

term by term: the constants give $-3+3=0$ , the $x$ terms $3+6-9=0$ , the $x^2$ terms $-1-\tfrac72+\tfrac92=0$ , the $x^3$ terms $\tfrac12-\tfrac12=0$ . A checkable number: at $x=2$ ,

L_3(2) = 1-6+6-\tfrac{8}{6} = -\tfrac13 = -0.333333\ldots

(b) A Bessel function in closed form. Take $\nu=\tfrac12$ in the representation of part 5, so $a=\nu+\tfrac12=1$ and $c=2\nu+1=2$ . First evaluate $M(1,2,w)$ exactly: since $(1)_n=n!$ and $(2)_n=(n+1)!$ ,

M(1,2,w) = \sum_{n\ge0}\frac{n!}{(n+1)!}\frac{w^n}{n!} = \sum_{n\ge0}\frac{w^n}{(n+1)!} = \frac{\mathrm{e}^{w}-1}{w} .

Now put $w=2\mathrm{i}x$ and use $\Gamma(\tfrac32)=\tfrac{\sqrt\pi}{2}$ :

J_{1/2}(x) = \frac{\mathrm{e}^{-\mathrm{i}x}}{\Gamma(3/2)}\left(\frac{x}{2}\right)^{1/2}\frac{\mathrm{e}^{2\mathrm{i}x}-1}{2\mathrm{i}x} = \frac{2}{\sqrt{\pi}}\sqrt{\frac{x}{2}}\cdot\frac{\mathrm{e}^{\mathrm{i}x}-\mathrm{e}^{-\mathrm{i}x}}{2\mathrm{i}x} = \frac{2}{\sqrt{\pi}}\cdot\frac{1}{\sqrt{2x}}\cdot\frac{\sin x}{1}\cdot\frac{\sqrt{x}\sqrt{x}}{x} ,

and simplifying the powers of $x$ ,

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

the familiar closed form. A checkable number at $x=1.3$ :

\sqrt{\frac{2}{\pi(1.3)}} = \sqrt{0.489707\ldots} = 0.699791\ldots, \qquad \sin(1.3)=0.963558\ldots,

J_{1/2}(1.3) = 0.699791\times0.963558 = 0.674289\ldots

and high-precision evaluation of the Bessel function gives $J_{1/2}(1.3) = 0.6742893968$ , agreeing to every digit shown. The identity $M(1,2,w)=(\mathrm{e}^w-1)/w$ was checked numerically as well, at $w=0.7$ : both sides equal $1.4482181535$ .

ApplicationOriginal

The hydrogen spectrum. The radial Schrödinger equation for a charge $Z$ with $V(r)=-Ze^2/r$ , written for $u=rR$ , is exactly the normal form of step 1 of the method, with $\lambda=2E<0$ , $A = 2Ze^2$ and $B=\ell(\ell+1)$ in the units $\hbar=m=1$ . Steps 2–4 then produce $f = M(\ell+1-\alpha,\,2\ell+2,\,z)$ with $\alpha = A/(2\sqrt{-\lambda})$ , and step 5 gives the quantization condition $\ell+1-\alpha = -N$ with $N\ge0$ an integer. Writing $n = N+\ell+1$ and restoring units, this is

E_n = -\frac{Z^2 m e^4}{2\hbar^2 n^2} ,

the Bohr spectrum, with $f$ an associated Laguerre polynomial and $R_{n\ell}(r)\propto r^{\ell}\mathrm{e}^{-Zr/na_0}M(-n+\ell+1,\,2\ell+2,\,2Zr/na_0)$ . This is Hassani's derivation, and what it shows is that the integer $n$ is not put in by hand anywhere: it is forced by the single requirement that an entire function of the form $\mathrm{e}^{-z/2}M$ stay bounded, i.e. that the Kummer series terminate. The accidental $\ell$ -independence of $E_n$ — the degeneracy peculiar to the Coulomb potential — appears as the fact that $\alpha$ , and hence the energy, depends on $N$ and $\ell$ only through their sum.

Coulomb scattering and the Gamow factor. Above threshold the same equation has $\lambda>0$ and no termination, and the physically relevant solutions are the regular and irregular Coulomb wave functions $F_\ell$ and $G_\ell$ , both built from $M$ and $U$ with imaginary arguments. The near-threshold penetration factor $\lvert M \rvert^2$ at the origin evaluates to the Gamow factor $2\pi\eta/(\mathrm{e}^{2\pi\eta}-1)$ , which sets the rate of every charged-particle fusion reaction in a star. It is the Gamma functions in $U$ — the same ones that appear in the definition above — that supply the exponential.

Three-dimensional harmonic trap, and why $U$ matters. For $V=\tfrac12m\omega^2r^2$ the same recipe gives $M(-N,\ell+\tfrac32,m\omega r^2/\hbar)$ , so the eigenfunctions are Laguerre polynomials in $r^2$ and the spectrum is $(2N+\ell+\tfrac32)\hbar\omega$ . But if the trap is cut off, or the problem is posed on the exterior of a sphere with only a decay condition at infinity, there is no condition at the origin to prefer $M$ , and the correct solution is $U$ . Its algebraic decay $x^{-a}$ , read straight off (IU), is then the statement that the wavefunction leaks out with a power-law tail rather than a Gaussian one.

Everyday special functions. The last two reductions in part 5 are why $\operatorname{erf}$ and the incomplete Gamma function (Error function and incomplete Gamma function) need no separate theory: their series, recurrences, asymptotics and continued fractions are the Kummer ones specialised. A numerical library that implements $M$ and $U$ well — typically by series for small $\lvert x\rvert$ , (IU) or an asymptotic expansion for large $\lvert x\rvert$ , and Kummer's transformation to avoid cancellation when $x<0$ — gets the error function, the incomplete Gamma function, the Coulomb functions, the Laguerre and Hermite polynomials and the modified Bessel functions in the same stroke. The cancellation point is worth stressing: evaluating $M(a,c,x)$ directly for $x$ large and negative subtracts nearly equal terms, and the standard cure is to compute $\mathrm{e}^{x}M(c-a,c,-x)$ instead, which is Kummer's transformation used as a numerical device.

Definition

Elliptic integrals

Adapted

Add to PDF

The incomplete and complete elliptic integrals of the first and second kinds, their binomial series and hypergeometric forms, and the exact period of a pendulum of arbitrary amplitude.

Prerequisites: Power series, The hypergeometric function

Definition(The four standard integrals, and the three conventions)AdaptedAdapted from Arfken, Weber & Harris (2013), §18.8, eqs. (18.169)–(18.174), pp. 928–929; Strogatz (1994), Exercise 6.7.4(c), the displayed definition of $K(m)$, p. 193; Ahlfors (1979), Ch. 6, §2.3, the displayed integrals for $K$ and $K'$, pp. 239–240

The incomplete elliptic integral of the first kind is

F(\varphi \,\backslash\, \alpha) = \int_0^{\varphi}\bigl(1-\sin^2\alpha\,\sin^2\theta\bigr)^{-1/2}\,\mathrm{d}\theta , \qquad\text{equivalently}\qquad F(x\,\vert\, m) = \int_0^{x}\Bigl[(1-t^2)(1-mt^2)\Bigr]^{-1/2}\mathrm{d}t ,

and the incomplete elliptic integral of the second kind is

E(\varphi\,\backslash\,\alpha) = \int_0^{\varphi}\bigl(1-\sin^2\alpha\,\sin^2\theta\bigr)^{1/2}\,\mathrm{d}\theta , \qquad E(x\,\vert\, m) = \int_0^{x}\left(\frac{1-mt^2}{1-t^2}\right)^{1/2}\mathrm{d}t .

The two forms of each are the same integral after $t=\sin\theta$ . Setting the upper limit to $\varphi=\pi/2$ , equivalently $x=1$ , gives the complete integrals

K(m) = \int_0^{\pi/2}\bigl(1-m\sin^2\theta\bigr)^{-1/2}\mathrm{d}\theta , \qquad E(m) = \int_0^{\pi/2}\bigl(1-m\sin^2\theta\bigr)^{1/2}\mathrm{d}\theta ,

with $0 \le m < 1$ for $K$ and $0 \le m \le 1$ for $E$ .

Three names for one number. The quantity $m$ in the displays above is the parameter; Arfken, Weber and Harris relate it to the modular angle $\alpha$ by $m = \sin^2\alpha$ , and it is also written $m = k^2$ , where $k$ is the modulus. All three are in use:

Written	Argument	Who
$K(m)$ , $E(m)$	parameter $m$	Arfken, Weber and Harris; Strogatz; most software, including the `ellipk`/`ellipe` of the standard numerical libraries
$K(k)$ , $E(k)$	modulus $k=\sqrt m$	much of the classical literature, and Arfken, Weber and Harris themselves in their exercises, where $K(k^2)$ is written
$F(\varphi\backslash\alpha)$	modular angle	AMS-55 tables

The backslash and the vertical bar in the AMS-55 notation are not decoration: they say which of the three the second argument is. Confusing $K(1/2)=1.8541$ with $K(k=1/2)=K(m=1/4)=1.6858$ is the commonest error in using these functions, and the only protection is to state the convention. In this unit the argument is always the parameter $m$ . Ahlfors, writing the integrals as $\int_0^1\mathrm{d}t/\sqrt{(1-t^2)(1-k^2t^2)}$ , uses the modulus; his $K$ is our $K(k^2)$ .

The integrals are called elliptic because the second kind first appeared as the arc length of an ellipse, not because their graphs are elliptical, and the adjective has since attached itself to the whole family — including the doubly periodic elliptic functions, which are the inverse functions of the first-kind integral.

Definition(Series, hypergeometric form, limits, and the pendulum period)AdaptedAdapted from Arfken, Weber & Harris (2013), §18.8, eqs. (18.175)–(18.179), pp. 929–930, and Example 18.8.1, eqs. (18.164)–(18.168), pp. 927–928; Strogatz (1994), Exercise 6.7.4, pp. 192–193

Series. For $0 \le m < 1$ , $K(m) = \frac{\pi}{2}\left\{1 + \sum_{n=1}^{\infty}\left[\frac{(2n-1)!!}{(2n)!!}\right]^2 m^n\right\} = \frac{\pi}{2}\left(1+\frac{m}{4}+\frac{9m^2}{64}+\frac{25m^3}{256}+\cdots\right),$ $E(m) = \frac{\pi}{2}\left\{1 - \sum_{n=1}^{\infty}\left[\frac{(2n-1)!!}{(2n)!!}\right]^2 \frac{m^n}{2n-1}\right\} = \frac{\pi}{2}\left(1-\frac{m}{4}-\frac{3m^2}{64}-\frac{5m^3}{256}-\cdots\right).$
Hypergeometric form. Both series are instances of the Gauss series of The hypergeometric function: $K(m) = \frac{\pi}{2}\,{}_2F_1\!\left(\tfrac12,\tfrac12;1;m\right), \qquad E(m) = \frac{\pi}{2}\,{}_2F_1\!\left(-\tfrac12,\tfrac12;1;m\right).$
Limits. $\lim_{m\to0}K(m) = \lim_{m\to0}E(m) = \frac{\pi}{2}, \qquad \lim_{m\to1}E(m) = 1, \qquad K(m)\to\infty \ \text{as}\ m\to1 ,$ the last divergence being logarithmic.
The pendulum. A rigid pendulum of length $l$ in gravity $g$ , released from rest at angle $\theta_M$ with $0<\theta_M<\pi$ , has exact period $T = 4\sqrt{\frac{l}{g}}\;K\!\left(\sin^2\frac{\theta_M}{2}\right) .$ In particular $T \to 2\pi\sqrt{l/g}$ as $\theta_M\to0$ and $T\to\infty$ as $\theta_M\to\pi$ .

Idea of the proofOriginal

Neither $K$ nor $E$ is elementary, and the reason is visible in the pendulum: the exact equation of motion is $\ddot\theta = -(g/l)\sin\theta$ , and the one first integral available — energy — turns the second-order equation into a quadrature $\mathrm{d}t \propto \mathrm{d}\theta/\sqrt{\cos\theta-\cos\theta_M}$ . That is as far as elementary methods go. The integral that remains is not elementary, so one names it, and the whole theory is the study of that name.

The series come from the only expansion available. In $K(m)$ the integrand is $(1-u)^{-1/2}$ with $u = m\sin^2\theta$ never exceeding $m<1$ , so the binomial series converges uniformly in $\theta$ and may be integrated term by term; each term is then $\int_0^{\pi/2}\sin^{2n}\theta\,\mathrm{d}\theta$ , a Wallis integral with the closed form $(\pi/2)(2n-1)!!/(2n)!!$ . Multiplying the two double-factorial ratios — one from the binomial coefficient, one from Wallis — is what produces the square in the answer, and that square is the signature of ${}_2F_1(\tfrac12,\tfrac12;1;m)$ .

The pendulum formula needs one further trick, and it is worth seeing why. After the half-angle identity the integral is $\int_0^{\theta_M}\mathrm{d}\theta/\sqrt{\sin^2(\theta_M/2)-\sin^2(\theta/2)}$ , which is singular at the upper limit and has the amplitude in both the integrand and the limit. The substitution $\sin(\theta/2) = \sin(\theta_M/2)\sin\varphi$ fixes both complaints at once: the upper limit becomes $\pi/2$ regardless of amplitude, and the square root that vanished at the endpoint cancels identically against the Jacobian. What is left is $K$ by definition. Why this amplitude ends up as $m=\sin^2(\theta_M/2)$ rather than anything else is the content of that single substitution.

ProofAdaptedAdapted from Arfken, Weber & Harris (2013), §18.8, the term-by-term integration leading to (18.175)–(18.176) and the comparison with §18.5 giving (18.177), p. 929, and Example 18.8.1, eqs. (18.164)–(18.168), pp. 927–928; Strogatz (1994), Exercise 6.7.4(a)–(c), pp. 192–193

1. The series for $K$ . Fix $0\le m<1$ . For $\lvert u\rvert<1$ the binomial series is

(1-u)^{-1/2} = \sum_{n=0}^{\infty}\frac{(2n-1)!!}{(2n)!!}\,u^{n},

with $(-1)!!=0!!=1$ ; the coefficients are $\binom{-1/2}{n}(-1)^n$ , and the double-factorial form follows by induction from the ratio $(2n-1)/(2n)$ between consecutive coefficients. Put $u = m\sin^2\theta$ . Then $0\le u\le m<1$ for every $\theta$ , so the series converges uniformly in $\theta$ by the Weierstrass test, dominated by the convergent numerical series $\sum_n \bigl[(2n-1)!!/(2n)!!\bigr]m^n$ . Uniform convergence permits term-by-term integration over the bounded interval $[0,\pi/2]$ (Power series):

K(m) = \sum_{n=0}^{\infty}\frac{(2n-1)!!}{(2n)!!}\,m^{n}\int_0^{\pi/2}\sin^{2n}\theta\,\mathrm{d}\theta .

The remaining integral is Wallis's: integrating by parts once gives the recursion $W_n = \frac{2n-1}{2n}W_{n-1}$ with $W_0=\pi/2$ , hence $W_n = \frac{\pi}{2}\frac{(2n-1)!!}{(2n)!!}$ . Substituting it produces the squared coefficient and the stated series. The same computation with the exponent $+\tfrac12$ , whose binomial coefficients are $-\bigl[(2n-1)!!/(2n)!!\bigr]/(2n-1)$ for $n\ge1$ , gives the series for $E$ ; Arfken, Weber and Harris leave this one as an exercise and we have reproduced it.

2. Hypergeometric identification. Compare coefficients. In ${}_2F_1(\tfrac12,\tfrac12;1;m) = \sum_n \bigl[(\tfrac12)_n\bigr]^2 m^n/\bigl[(1)_n n!\bigr]$ we have $(1)_n = n!$ and

\left(\tfrac12\right)_n = \tfrac12\cdot\tfrac32\cdots\tfrac{2n-1}{2} = \frac{(2n-1)!!}{2^n}, \qquad n! = \frac{(2n)!!}{2^n},

the second because $(2n)!! = 2^n n!$ . Hence $(\tfrac12)_n/n! = (2n-1)!!/(2n)!!$ and the $n$ th coefficient of the Gauss series is exactly $\bigl[(2n-1)!!/(2n)!!\bigr]^2$ , which is the coefficient in part 1. For $E$ the numerator parameter $-\tfrac12$ gives $(-\tfrac12)_n = -\frac{(2n-3)!!}{2^n}$ for $n\ge1$ , and $(2n-3)!! = (2n-1)!!/(2n-1)$ , which supplies the extra factor $-1/(2n-1)$ . Both identities hold wherever both sides converge, that is on $\lvert m\rvert<1$ .

3. Limits. At $m=0$ both integrands are $1$ and the integrals are $\pi/2$ . At $m=1$ the second integrand is $\lvert\cos\theta\rvert$ , so $E(1)=\int_0^{\pi/2}\cos\theta\,\mathrm{d}\theta=1$ . For $K$ near $m=1$ , the integrand $(1-m\sin^2\theta)^{-1/2}$ develops a non-integrable square-root singularity at $\theta=\pi/2$ in the limit; for $m<1$ write $1-m\sin^2\theta = (1-m)+m\cos^2\theta$ , so near $\theta=\pi/2$ the integrand behaves like $[(1-m)+m(\pi/2-\theta)^2]^{-1/2}$ , whose integral over a neighbourhood is $\sim m^{-1/2}\ln\bigl[1/\sqrt{1-m}\bigr]$ . The divergence is therefore logarithmic in $1-m$ , as Arfken, Weber and Harris state. Note that the series of part 1 is useless here — its coefficients tend to zero only like $1/n$ , so at $m=1$ it is a divergent harmonic-type series, which is precisely the divergence showing itself.

4. The pendulum period. Let $\theta$ be the angle from the downward vertical. Taking the potential energy to be $-mgl\cos\theta$ and the kinetic energy $\tfrac12 ml^2\dot\theta^2$ , and evaluating the total at the release point, where $\dot\theta=0$ and $\theta=\theta_M$ ,

\tfrac12 ml^2\dot\theta^2 - mgl\cos\theta = -mgl\cos\theta_M ,

so the mass cancels and

\dot\theta = \pm\sqrt{\frac{2g}{l}}\,\bigl(\cos\theta-\cos\theta_M\bigr)^{1/2} .

Take the quarter period from $\theta=0$ (moving in the positive direction) to $\theta=\theta_M$ , where the motion reverses. Separating variables,

\frac{T}{4} = \sqrt{\frac{l}{2g}}\int_0^{\theta_M}\frac{\mathrm{d}\theta}{\sqrt{\cos\theta-\cos\theta_M}} .

That this quarter arc really is a quarter of the period is the symmetry of the problem under $t\mapsto-t$ and under $\theta\mapsto-\theta$ , both of which the equation of motion respects.

Now the half-angle identity $\cos\theta = 1-2\sin^2(\theta/2)$ turns the radicand into $2\bigl[\sin^2(\theta_M/2)-\sin^2(\theta/2)\bigr]$ , so

\frac{T}{4} = \frac12\sqrt{\frac{l}{g}}\int_0^{\theta_M}\frac{\mathrm{d}\theta}{\sqrt{\sin^2(\theta_M/2)-\sin^2(\theta/2)}} .

Write $\sqrt m = \sin(\theta_M/2)$ and substitute

\sin\frac{\theta}{2} = \sqrt m\,\sin\varphi ,

which is a genuine change of variable because $\theta_M<\pi$ makes $\sin(\theta/2)$ strictly increasing on $[0,\theta_M]$ , and which carries $\theta=0$ to $\varphi=0$ and $\theta=\theta_M$ to $\varphi=\pi/2$ . Differentiating, $\tfrac12\cos(\theta/2)\,\mathrm{d}\theta = \sqrt m\cos\varphi\,\mathrm{d}\varphi$ , and $\cos(\theta/2) = \sqrt{1-m\sin^2\varphi}$ , so

\mathrm{d}\theta = \frac{2\sqrt m\cos\varphi}{\sqrt{1-m\sin^2\varphi}}\,\mathrm{d}\varphi ,

while the denominator becomes $\sqrt{m-m\sin^2\varphi} = \sqrt m\cos\varphi$ . The factors $\sqrt m\cos\varphi$ cancel exactly — this is the point of the substitution — leaving

\frac{T}{4} = \frac12\sqrt{\frac{l}{g}}\int_0^{\pi/2}\frac{2\,\mathrm{d}\varphi}{\sqrt{1-m\sin^2\varphi}} = \sqrt{\frac{l}{g}}\,K(m) ,

which is the assertion. Arfken, Weber and Harris carry out the same two steps, remarking that the half-angle substitution takes "a bit of clairvoyance"; Strogatz's Exercise 6.7.4 is this derivation split into the three parts (a), (b), (c). The small-amplitude limit is $m\to0$ , where $K\to\pi/2$ and $T\to2\pi\sqrt{l/g}$ ; the large-amplitude limit is $\theta_M\to\pi$ , where $m\to1$ and $K$ diverges, the physical statement that a pendulum balanced at the top takes forever to fall. $\square$

Remark(Where the hypotheses are used)Original

$m<1$ for the series. Uniform convergence of the binomial expansion needs $\sup_\theta m\sin^2\theta = m < 1$ . At $m=1$ the expansion still converges pointwise for $\theta<\pi/2$ but not uniformly, and term-by-term integration fails — correctly, since the sum diverges while every term is finite. This is not a gap in the proof but the logarithmic divergence of part 3 seen from the series side.
$m\le1$ for $E$ , $m<1$ for $K$ . The integrand of $E$ stays bounded as $m\to1$ (it becomes $\cos\theta$ ), so $E$ extends continuously to $m=1$ ; the integrand of $K$ does not.
$0<\theta_M<\pi$ in the pendulum. The lower bound is needed for the problem to be non-trivial; the upper bound is needed twice. First, $\sin(\theta_M/2)<1$ keeps $m<1$ , so that $K(m)$ is finite. Second, the substitution $\sin(\theta/2)=\sqrt m\sin\varphi$ must be invertible on the range of motion, which requires $\theta/2$ to stay in $[0,\pi/2)$ where the sine is injective. At $\theta_M=\pi$ both fail together, and the solution is no longer periodic at all: it is the separatrix, on which the bob approaches the inverted position asymptotically and never returns. The formula's divergence is therefore the right answer, not a breakdown.
Rigid rod, no damping, no driving. The energy integral is the whole method, and it exists only because the system is conservative. Add the smallest linear damping and $T$ is no longer defined; the motion spirals in. A pendulum on a string rather than a rod additionally requires the string to stay taut, which holds for $\theta_M\le\pi/2$ but not beyond.
Which $m$ ? Every statement above is in the parameter convention. In the modulus convention the pendulum period is $4\sqrt{l/g}\,K(k)$ with $k=\sin(\theta_M/2)$ , which looks different and is the same. Mixing the two silently is how a factor appears in a period and is then blamed on the physics.

Method(Reducing an integral to standard form)AdaptedAdapted from Arfken, Weber & Harris (2013), §18.8, the remark on integrals of a rational function of $t$ and the square root of a quartic, p. 930, and the definitions (18.169)–(18.174), pp. 928–929

Arfken, Weber and Harris record the general fact: any integral

\int_0^{x} R\!\left(t,\;\sqrt{a_4t^4+a_3t^3+a_2t^2+a_1t+a_0}\right)\mathrm{d}t ,

with $R$ rational in its two arguments, can be expressed in terms of elliptic integrals (and elementary functions). In practice one proceeds as follows.

Check the degree. A square root of a quadratic gives an elementary answer (inverse trigonometric or inverse hyperbolic); a cubic or a quartic gives elliptic integrals; a quintic or higher generally does not — those are hyperelliptic and have no classical closed form.
Factor the quartic into real linear and quadratic factors. The positions of the roots relative to the range of integration decide which standard form is reachable.
Normalise by a Möbius substitution $t\mapsto(\alpha u+\beta)/(\gamma u+\delta)$ chosen to send two of the roots to $\pm1$ and the other two to $\pm1/k$ . This is always possible and is what brings the radicand to the Legendre form $(1-u^2)(1-k^2u^2)$ .
Split the rational part. With the radicand in standard form, write $R$ as a sum of terms of the shapes $1/\sqrt{\cdot}$ , $u^2/\sqrt{\cdot}$ and $1/[(1-nu^2)\sqrt{\cdot}]$ . These integrate to the first, second and third kinds respectively.
Recognise a trigonometric problem directly. When the integral already has the shape $\int\mathrm{d}\theta/\sqrt{1-m\sin^2\theta}$ — as the pendulum does after the half-angle step — stop: it is $F(\varphi\backslash\alpha)$ with $m=\sin^2\alpha$ , complete if $\varphi=\pi/2$ .
Evaluate. Use the series of part 1 of the statement for $m\lesssim0.5$ , where it converges quickly; for $m$ near $1$ use the arithmetic–geometric mean, which converges quadratically for every $m<1$ and is what numerical libraries actually implement, or the logarithmic expansion about $m=1$ .

A caution on step 6. The series is a poor tool near $m=1$ for the reason given in the remark; conversely the arithmetic–geometric mean iteration $a_{j+1}=(a_j+b_j)/2$ , $b_{j+1}=\sqrt{a_jb_j}$ started from $a_0=1$ , $b_0=\sqrt{1-m}$ converges so fast that $K(m)=\pi/(2\,\mathrm{AGM})$ is accurate to machine precision in four or five steps. It degrades only as $m\to1$ , where $b_0\to0$ .

Example(A pendulum released from the horizontal)Original

Take a rigid pendulum of length $l=1.000\ \mathrm{m}$ in $g=9.80665\ \mathrm{m\,s^{-2}}$ , released from rest at $\theta_M = 90^\circ$ .

Step 1: the small-angle period, for comparison.

T_0 = 2\pi\sqrt{\frac{l}{g}} = 2\pi\sqrt{0.1019716} = 2\pi(0.3193299) = 2.006409\ \mathrm{s}.

Step 2: the parameter. $\theta_M/2 = 45^\circ$ , so

m = \sin^2 45^\circ = \left(\frac{\sqrt2}{2}\right)^2 = \frac12 .

Step 3: evaluate $K(1/2)$ from the series. With $m=\tfrac12$ the terms of part 1 of the statement are

1,\quad \frac{m}{4}=0.125,\quad \frac{9m^2}{64}=0.035156,\quad \frac{25m^3}{256}=0.012207,\quad \frac{1225 m^4}{16384}=0.004673,\ \dots

Partial sums times $\pi/2$ give $1.570796,\ 1.767146,\ 1.822369,\ 1.841544,\ 1.848884,\ \dots$ , and after twelve terms $1.854056$ . The limit, which the arithmetic–geometric mean returns at once, is

K(1/2) = 1.8540746773 .

(As a check that nothing has been dropped, direct numerical quadrature of $\int_0^{\pi/2}(1-\tfrac12\sin^2\theta)^{-1/2}\mathrm{d}\theta$ gives the same ten digits, as does $\tfrac{\pi}{2}\,{}_2F_1(\tfrac12,\tfrac12;1;\tfrac12)$ from part 2.)

Step 4: the period.

T = 4\sqrt{\frac{l}{g}}\,K(1/2) = 4(0.3193299)(1.8540747) = 2.368246\ \mathrm{s} ,

\frac{T}{T_0} = \frac{2}{\pi}K(1/2) = 1.1803406 .

A pendulum released from the horizontal runs $18.03\%$ slow compared with the small-angle formula. For a seconds pendulum that is more than half an hour a day.

Step 5: how good is the leading correction? Expanding $m = \sin^2(\theta_M/2) = \theta_M^2/4-\theta_M^4/48+\cdots$ and feeding it into the series of part 1 gives

T = T_0\left(1+\frac{\theta_M^2}{16}+\frac{11\,\theta_M^4}{3072}+\cdots\right) ,

the $\theta_M^4$ coefficient coming from $-\tfrac{1}{192}+\tfrac{9}{1024} = \tfrac{11}{3072}$ . At $\theta_M=\pi/2$ the first correction alone gives $1+\pi^2/64 = 1.15421$ , which captures only $85.5\%$ of the true correction; adding the quartic term gives $1.17601$ , and the exact answer is $1.18034$ . The moral is that the small-angle expansion is an expansion in $\theta_M^2/16$ , so it is excellent below $20^\circ$ (error under $10^{-4}$ ) and merely indicative at $90^\circ$ .

A second number, for the other kind. The quarter arc of the ellipse $x^2/a^2+y^2/b^2=1$ is $a\,E(m)$ with $m=(a^2-b^2)/a^2$ . For $a=5$ , $b=3$ , $m=16/25=0.64$ and $E(0.64)=1.27634994$ , so the full perimeter is

P = 4aE(m) = 20\times1.27634994 = 25.526999 ,

against the crude estimate $2\pi\sqrt{(a^2+b^2)/2} = 2\pi\sqrt{17} = 25.906237$ , which is $1.5\%$ too large. Direct quadrature of $4\int_0^{\pi/2}\sqrt{25\cos^2 t+9\sin^2 t}\,\mathrm{d}t$ returns $25.526999$ , confirming the identification. Every number in this example was computed to at least ten significant figures and checked against an independent evaluation.

ApplicationOriginal

Clocks and the circular error. The amplitude dependence computed above is the circular error of a pendulum clock, and it is the reason a clock's rate drifts as its driving force or its pivot friction changes the swing. Huygens's cycloidal cheeks were an attempt to remove it geometrically by forcing the bob onto a cycloid, on which the period is exactly amplitude-independent; in practice the cheeks introduced more friction than they removed error, and the surviving fix is the opposite one — keep the amplitude small and constant. The quantitative design rule follows straight from the expansion: to hold a clock to one second a day, that is $1.16\times10^{-5}$ in rate, the amplitude must be held to $\delta\theta_M$ with $\theta_M\,\delta\theta_M/8 < 1.2\times10^{-5}$ , so at a typical $\theta_M = 3^\circ = 0.0524$ rad the amplitude must be stable to about $0.0018$ rad, a tenth of a degree.

Orbits and arcs. The perimeter of an ellipse has no elementary formula, which is why $E$ exists. In celestial mechanics this means the arc length travelled along a Keplerian orbit, and the path length of a charged particle in a dipole field, are both second-kind integrals. The same function governs the length of a sagging cable and the shape of the elastica, the curve taken by a thin elastic rod bent by end loads — Euler's problem, whose solution is written with $F$ and $E$ and whose buckling thresholds are the values of $m$ at which a complete integral reaches a critical value.

Fields of circular sources. Let a ring of radius $a$ carry a steady current $I$ . Its azimuthal vector potential at the field point $(\rho,z)$ is

A_\varphi = \frac{\mu_0 I}{\pi k}\sqrt{\frac{a}{\rho}}\left[\left(1-\frac{k^2}{2}\right)K(k^2)-E(k^2)\right], \qquad k^2 = \frac{4a\rho}{(a+\rho)^2+z^2} ,

which Arfken, Weber and Harris set as Exercise 18.8.4 — note that the complete integrals there are written with the modulus squared as the argument, so in our convention the argument is $m=k^2$ throughout. The same combination appears for the mutual inductance of two coaxial loops and for the gravitational potential of a ring; the near-field limit $\rho\to a$ , $z\to0$ is $m\to1$ , and the logarithmic divergence of $K$ is exactly the logarithmic divergence of the field at a line current.

The gateway to elliptic functions. Ahlfors arrives at the same integrals from a completely different direction: the Schwarz–Christoffel map carrying the upper half-plane onto a rectangle is an elliptic integral, and the rectangle's two side lengths are $K$ and $K'$ , the complete integrals for the modulus $k$ and for the complementary modulus $\sqrt{1-k^2}$ . Reflecting the map across the sides of the rectangle extends its inverse to a meromorphic function on the whole plane with periods $2K$ and $2\mathrm{i}K'$ — a doubly periodic, or elliptic, function. That is the structural reason elliptic integrals keep appearing in integrable mechanics: the motion of a symmetric top, of a pendulum, and of a free rigid body are all linear flows on a torus, and a torus is a plane modulo a lattice. The lattice periods are $2K$ and $2\mathrm{i}K'$ , and the real one is the physical period computed above.

Definition

Airy functions

Adapted

Add to PDF

Airy's equation $y''=zy$ , its two standard solutions $\mathrm{Ai}$ and $\mathrm{Bi}$ as series, contour integrals and Bessel functions of order $1/3$ , their opposite behaviour on the two sides of the origin, and the role of $\mathrm{Ai}$ at a WKB turning point.

Prerequisites: Bessel functions of the first kind, The Gamma function, Power-series solutions at ordinary points

Definition(Airy's equation, and the functions

\mathrm{Ai}

and

\mathrm{Bi}

)AdaptedAdapted from Riley, Hobson & Bence (2006), §25.6, eq. (25.32) and the behaviour of the solutions near the origin, pp. 888–889, §25.6.2, the two series, p. 890, and §25.6.3, eqs. (25.35)–(25.36) and (25.39), pp. 891–892; Hassani (2013), Problem 14.24, p. 456

Airy's equation is

\frac{\mathrm{d}^2y}{\mathrm{d}z^2} = z\,y . \tag{A}

Riley, Hobson and Bence call it Stokes' equation and list the variants $y''+\lambda xy=0$ and $y''+xy=0$ , all of which reduce to (A) by a complex rescaling of the independent variable; we use (A) and write $\mathrm{Ai}(-x)$ when the oscillatory side is wanted.

Every finite point is an ordinary point of (A) — the only singularity is at $z=\infty$ , and it is irregular — so by Power-series solutions at ordinary points there are two power-series solutions about the origin, and both are entire. Writing a solution as $y=\sum_{n\ge0}c_nz^n$ and matching the coefficient of $z^{n+1}$ on the two sides of (A) gives the single recurrence

(n+3)(n+2)\,c_{n+3} = c_n \quad (n\ge0), \qquad\text{together with}\qquad 2c_2 = 0

from the coefficient of $z^0$ ; so $c_2=c_5=c_8=\cdots=0$ , and the two seeds $(c_0,c_1)=(1,0)$ and $(c_0,c_1)=(0,1)$ give

y_1(z) = 1 + \frac{z^3}{6} + \frac{z^6}{180} + \cdots = \sum_{k=0}^{\infty}\frac{z^{3k}}{\prod_{j=1}^{k}(3j-1)(3j)} , \qquad y_2(z) = z + \frac{z^4}{12} + \frac{z^7}{504} + \cdots = \sum_{k=0}^{\infty}\frac{z^{3k+1}}{\prod_{j=1}^{k}(3j)(3j+1)} .

Hassani sets the same computation for $y''+xy=0$ and asks for the radius of convergence, which is infinite by the ratio test: successive terms carry a factor $z^3/[(3k+2)(3k+3)]$ , which tends to zero for every $z$ .

These two series are a basis, but not a useful one, because both of them blow up as $z\to+\infty$ . The standard basis is chosen by behaviour at infinity instead:

$\mathrm{Ai}(z)$ is the solution that decays to zero as $z\to+\infty$ along the real axis, normalised as below;
$\mathrm{Bi}(z)$ is the solution that grows there, normalised to be in quadrature with $\mathrm{Ai}$ on the other side.

Riley, Hobson and Bence motivate the first by a curvature argument that is worth repeating because it shows the solution is unique up to scale. On $z>0$ the equation says $y''$ has the sign of $y$ , so wherever $y>0$ the graph is convex and bends away from the axis. If $y(0)>0$ and $y'(0)$ is positive or only slightly negative, the graph is carried up to $+\infty$ ; if $y'(0)$ is very negative the graph crosses the axis, after which the curvature reverses and carries it down to $-\infty$ . Exactly one intermediate slope gives a graph that approaches the axis from above and never crosses it. That solution is $\mathrm{Ai}$ . On $z<0$ the sign of $y''$ is opposite to that of $y$ , so every solution bends back towards the axis, crosses it, and oscillates — with frequency increasing and amplitude decreasing as $-z$ grows. Hassani's Problem 14.24 asks for the same dichotomy, proved by the Sturm comparison theorem: infinitely many zeros on the oscillatory side, at most one on the other.

Integral definitions. The integrand $\exp(-\tfrac13t^3+zt)$ decays at infinity exactly in the three sectors $\lvert\arg t-2\pi n/3\rvert<\pi/6$ , $n=0,1,2$ , where $\operatorname{Re}(t^3)>0$ . Let $C_1$ be a contour in the complex $t$ -plane coming in from infinity in the sector about $\arg t = 4\pi/3$ and going back out to infinity in the sector about $\arg t = 2\pi/3$ , and let $C_2$ , $C_3$ be the two contours obtained from it by rotating through $\pm2\pi/3$ . Then

\mathrm{Ai}(z) = \frac{1}{2\pi\mathrm{i}}\int_{C_1}\exp\!\left(-\tfrac13 t^3+zt\right)\mathrm{d}t , \qquad \mathrm{Bi}(z) = \frac{1}{2\pi}\int_{C_2}\exp\!\left(-\tfrac13 t^3+zt\right)\mathrm{d}t - \frac{1}{2\pi}\int_{C_3}\exp\!\left(-\tfrac13 t^3+zt\right)\mathrm{d}t .

Deforming $C_1$ onto the imaginary axis, $t=\mathrm{i}s$ , and discarding the odd part of the integrand gives the real form

\mathrm{Ai}(x) = \frac{1}{\pi}\int_0^{\infty}\cos\!\left(\tfrac13 s^3+xs\right)\mathrm{d}s \qquad (x\in\mathbb{R}),

which shows at once that $\mathrm{Ai}$ is real on the real axis. Strauss writes the same object, without the symmetrisation, as $A(\xi)=\frac{1}{2\pi}\int_{-\infty}^{\infty}\mathrm{e}^{\mathrm{i}k\xi+\mathrm{i}k^3/3}\mathrm{d}k$ .

A name used twice. Evans's catalogue of linear partial differential equations lists "Airy's equation" as $u_t+u_{xxx}=0$ . That is a different equation — third order, a PDE, and dispersive — and the connection is that its solution kernel is the function of this unit. Context disambiguates, but the clash is real.

Definition(Values, Wronskian, Bessel form and asymptotics)AdaptedAdapted from Riley, Hobson & Bence (2006), §25.6.3, the evaluation of Ai(0), pp. 893–894, the asymptotic forms (25.40)–(25.41), p. 893, the remark on Bessel functions of order 1/3, pp. 894–895, and §25.7, eqs. (25.52)–(25.53), p. 899; Hassani (2013), Problem 15.41, p. 490

Values at the origin, and the series. $\mathrm{Ai}(0) = \frac{3^{-2/3}}{\Gamma(2/3)} = \frac{3^{-1/6}\,\Gamma(1/3)}{2\pi} = 0.3550280539\ldots, \qquad \mathrm{Ai}'(0) = -\frac{3^{-1/3}}{\Gamma(1/3)} = -0.2588194038\ldots,$ $\mathrm{Bi}(0) = \frac{3^{-1/6}}{\Gamma(2/3)} = \sqrt3\,\mathrm{Ai}(0), \qquad \mathrm{Bi}'(0) = \frac{3^{1/6}}{\Gamma(1/3)} = -\sqrt3\,\mathrm{Ai}'(0) ,$ and consequently $\mathrm{Ai} = \mathrm{Ai}(0)\,y_1 + \mathrm{Ai}'(0)\,y_2$ and $\mathrm{Bi} = \mathrm{Bi}(0)\,y_1+\mathrm{Bi}'(0)\,y_2$ in terms of the two series above.
Wronskian. $\mathrm{Ai}(z)\mathrm{Bi}'(z)-\mathrm{Ai}'(z)\mathrm{Bi}(z) = 1/\pi$ , constant in $z$ ; in particular the pair is a basis.
Bessel form. Put $\zeta = \tfrac23 x^{3/2}$ for $x>0$ . Then $\sqrt x\,\mathcal{Z}_{1/3}(\zeta)$ solves Airy's equation for any solution $\mathcal{Z}_{1/3}$ of the modified Bessel equation of order $1/3$ , and $\sqrt x\,C_{1/3}(\zeta)$ solves $y''+xy=0$ for any solution $C_{1/3}$ of the ordinary Bessel equation of order $1/3$ . With the standard normalisations, $\mathrm{Ai}(x) = \frac{1}{\pi}\sqrt{\frac{x}{3}}\;K_{1/3}(\zeta) = \frac{\sqrt x}{3}\bigl[I_{-1/3}(\zeta)-I_{1/3}(\zeta)\bigr], \qquad \mathrm{Bi}(x) = \sqrt{\frac{x}{3}}\,\bigl[I_{-1/3}(\zeta)+I_{1/3}(\zeta)\bigr],$ $\mathrm{Ai}(-x) = \frac{\sqrt x}{3}\bigl[J_{1/3}(\zeta)+J_{-1/3}(\zeta)\bigr], \qquad \mathrm{Bi}(-x) = \sqrt{\frac{x}{3}}\,\bigl[J_{-1/3}(\zeta)-J_{1/3}(\zeta)\bigr].$
Asymptotics. As $x\to+\infty$ through real values, $\mathrm{Ai}(x) \sim \frac{1}{2\sqrt{\pi}\,x^{1/4}}\exp\!\left(-\tfrac23 x^{3/2}\right), \qquad \mathrm{Ai}(-x) \sim \frac{1}{\sqrt{\pi}\,x^{1/4}}\sin\!\left(\tfrac23 x^{3/2}+\frac{\pi}{4}\right),$ and $\mathrm{Bi}(-x) \sim \pi^{-1/2}x^{-1/4}\cos\bigl(\tfrac23x^{3/2}+\pi/4\bigr)$ , in quadrature with $\mathrm{Ai}(-x)$ as the choice of $\mathrm{Bi}$ was designed to make it.
Turning points. Let $k^2(x) = 2m[E-V(x)]/\hbar^2$ vanish simply at $x=a$ , with $k^2>0$ to the left and $k^2<0$ to the right. Near $a$ the Schrödinger equation $\psi''+k^2\psi=0$ becomes Airy's equation in the stretched variable $z = (x-a)/\mu$ with $\mu = \bigl[\hbar^2/2mV'(a)\bigr]^{1/3}$ , and the solution decaying on the forbidden side is $\mathrm{Ai}(z)$ . For a well with two simple turning points $\pm a$ , matching the $\mathrm{Ai}$ behaviour at each end to the WKB solution between them forces $\int_{-a}^{a}k(x)\,\mathrm{d}x = \left(n+\tfrac12\right)\pi , \qquad n=0,1,2,\dots$

Idea of the proofOriginal

Airy's equation is the simplest equation whose character changes across a point: $y''=zy$ is exponential where $z>0$ and oscillatory where $z<0$ , with the change happening at a single simple zero of the coefficient. Everything about the functions follows from that one fact, and so does their importance — any second-order equation $y''+k^2(x)y=0$ whose coefficient has a simple zero looks like Airy's equation when you zoom in on that zero and rescale. The scaling exponent $1/3$ in $\mu=[\hbar^2/2mV'(a)]^{1/3}$ is forced by dimensional balance between $y''$ and $zy$ , and it is the same $1/3$ that appears in the order of the Bessel functions, in the $t^{1/3}$ of the dispersive kernel, and in the $\tfrac23 x^{3/2}$ of the asymptotics (which is $\int\sqrt{x}\,\mathrm{d}x$ ).

The two constructions of the functions answer two different questions. The power series answers "what are they near zero", and is useless far out. The contour integral answers "what are they far out", and is the natural object because the equation has an irregular singular point at infinity: solutions there are Laplace-type integrals $\int f(t)\mathrm{e}^{zt}\mathrm{d}t$ , and substituting into (A) turns the differential equation for $y$ into a first-order equation for $f$ , namely $f'+t^2f=0$ , solved by $f=\mathrm{e}^{-t^3/3}$ . The only freedom left is the contour, and the requirement that the boundary terms vanish admits exactly three directions to infinity, in which $\operatorname{Re}t^3>0$ . Three endpoints give three contours, of which only two are independent because together the three form a closed loop. That is where the two named functions come from, and why they are naturally indexed by sectors of the complex plane.

The Bessel connection is the observation that $\int\sqrt{z}\,\mathrm{d}z = \tfrac23z^{3/2}$ is the phase, so the substitution $\zeta=\tfrac23z^{3/2}$ should straighten the equation out; it does, up to a factor $\sqrt z$ needed to remove the first-derivative term the substitution creates, and the residual order $\nu=1/3$ is $1/(2\cdot3/2)$ , the reciprocal of twice the exponent.

ProofAdaptedAdapted from Riley, Hobson & Bence (2006), §25.6.3, the Laplace-type construction (25.33)–(25.36), the deformation to (25.39) and the Gamma-integral evaluation of Ai(0), pp. 891–894, §25.7, the WKB solutions and the matching (25.50)–(25.53), pp. 898–899, and the turning-point worked example, pp. 900–901; Hassani (2013), Problem 15.41, p. 490

The contour solutions. Try $y(z)=\int_a^b f(t)\mathrm{e}^{zt}\,\mathrm{d}t$ . Then $y''=\int t^2f\mathrm{e}^{zt}$ , while

zy = \int_a^b f(t)\,z\mathrm{e}^{zt}\,\mathrm{d}t = \Bigl[f(t)\mathrm{e}^{zt}\Bigr]_a^b - \int_a^b f'(t)\mathrm{e}^{zt}\,\mathrm{d}t

after one integration by parts. So (A) holds provided the bracket vanishes and $f'+t^2f=0$ , that is $f(t)=A\mathrm{e}^{-t^3/3}$ . The bracket vanishes only if $\lvert t\rvert\to\infty$ with $\operatorname{Re}(t^3)>0$ , which holds in the three sectors $\lvert \arg t - 2\pi n/3\rvert<\pi/6$ , $n=0,1,2$ . A contour must begin in one sector and end in a different one, or the integral is zero; there are three such contours and their sum is a closed loop, so only two are independent. Defining $\mathrm{Ai}$ and $\mathrm{Bi}$ by the combinations in the definition fixes the basis.

The real form. Deform $C_1$ to run along the imaginary axis (permissible: the integrand is entire, and the ends can be swung inside the sectors without crossing a singularity). With $t=\mathrm{i}s$ , $\mathrm{d}t=\mathrm{i}\,\mathrm{d}s$ , and $-\tfrac13t^3+zt = \mathrm{i}(\tfrac13s^3+zs)$ , so

\mathrm{Ai}(z) = \frac{1}{2\pi}\int_{-\infty}^{\infty}\exp\!\left[\mathrm{i}\left(\tfrac13s^3+zs\right)\right]\mathrm{d}s .

For real $z$ the exponent is an odd function of $s$ , so the imaginary part of the integrand is odd and integrates to zero, leaving the cosine form with the range halved.

The value at the origin. Set $z=0$ and take $C_1$ as two rays, in along $\arg t = 4\pi/3$ and out along $\arg t=2\pi/3$ . On each ray put $\tfrac13t^3=s$ with $s>0$ real, so that $t=(3s)^{1/3}\mathrm{e}^{\mathrm{i}\vartheta}$ and $\mathrm{d}t = \mathrm{e}^{\mathrm{i}\vartheta}(3s)^{-2/3}\mathrm{d}s$ with $\vartheta=4\pi/3$ on the first ray (traversed inwards) and $\vartheta=2\pi/3$ on the second. Then

\mathrm{Ai}(0) = \frac{3^{-2/3}}{2\pi\mathrm{i}}\left(-\mathrm{e}^{4\pi\mathrm{i}/3}+\mathrm{e}^{2\pi\mathrm{i}/3}\right)\int_0^{\infty}\mathrm{e}^{-s}s^{-2/3}\,\mathrm{d}s .

The bracket is $-(-\tfrac12-\tfrac{\sqrt3}{2}\mathrm{i})+(-\tfrac12+\tfrac{\sqrt3}{2}\mathrm{i}) = \mathrm{i}\sqrt3$ , and the integral is $\Gamma(1/3)$ by the Euler integral of The Gamma function. Hence

\mathrm{Ai}(0) = \frac{3^{-2/3}\sqrt3}{2\pi}\,\Gamma(1/3) = \frac{3^{-1/6}\Gamma(1/3)}{2\pi} .

This is Riley, Hobson and Bence's worked evaluation. The alternative form $3^{-2/3}/\Gamma(2/3)$ follows from the reflection formula $\Gamma(1/3)\Gamma(2/3)=\pi/\sin(\pi/3)=2\pi/\sqrt3$ : substituting $\Gamma(1/3)=2\pi/[\sqrt3\,\Gamma(2/3)]$ into the display gives $3^{-1/6}\cdot2\pi/[\sqrt3\,\Gamma(2/3)\cdot2\pi] = 3^{-2/3}/\Gamma(2/3)$ . The value of $\mathrm{Ai}'(0)$ comes from the same contour with an extra factor $t$ in the integrand, which turns the exponent $-2/3$ into $-1/3$ and the Gamma function into $\Gamma(2/3)$ , and reverses the sign of the angular factor.

The Wronskian. $W=\mathrm{Ai}\,\mathrm{Bi}'-\mathrm{Ai}'\mathrm{Bi}$ satisfies $W' = \mathrm{Ai}\,\mathrm{Bi}''-\mathrm{Ai}''\mathrm{Bi} = z\,\mathrm{Ai}\,\mathrm{Bi}-z\,\mathrm{Ai}\,\mathrm{Bi} = 0$ , since both functions satisfy (A) and the equation has no first-derivative term to damp the Wronskian. So $W$ equals its value at the origin,

W = \frac{3^{-2/3}}{\Gamma(2/3)}\cdot\frac{3^{1/6}}{\Gamma(1/3)} + \frac{3^{-1/3}}{\Gamma(1/3)}\cdot\frac{3^{-1/6}}{\Gamma(2/3)} = \frac{3^{-1/2}+3^{-1/2}}{\Gamma(1/3)\Gamma(2/3)} = \frac{2/\sqrt3}{2\pi/\sqrt3} = \frac{1}{\pi},

the second term being positive because $\mathrm{Ai}'(0)$ and $\mathrm{Bi}(0)$ enter with a minus sign each. Equivalently $W = -2\sqrt3\,\mathrm{Ai}(0)\mathrm{Ai}'(0)$ , which is $-2\sqrt3\,(0.3550281)(-0.2588194) = 0.3183099 = 1/\pi$ .

The Bessel form. Follow Hassani's scheme. He records that the substitutions $z=\beta t^{\gamma}$ for the independent variable and $w=t^{\alpha}u$ for the dependent one — so that $w=\mathcal{Z}_{\nu}(z)$ is the Bessel function and $u=t^{-\alpha}w$ is the new unknown, a function of $t$ — turn Bessel's equation into

t^2\frac{\mathrm{d}^2u}{\mathrm{d}t^2}+(2\alpha+1)t\frac{\mathrm{d}u}{\mathrm{d}t}+\bigl(\beta^2\gamma^2t^{2\gamma}+\alpha^2-\nu^2\gamma^2\bigr)u = 0 .

We want this to be $t^2u''-t^3u=0$ , which is Airy's equation multiplied by $t^2$ . Matching term by term:

2\alpha+1 = 0 \Rightarrow \alpha=-\tfrac12; \qquad 2\gamma = 3 \Rightarrow \gamma = \tfrac32; \qquad \beta^2\gamma^2 = -1 \Rightarrow \beta = \pm\tfrac{2\mathrm{i}}{3}; \qquad \alpha^2=\nu^2\gamma^2 \Rightarrow \nu = \pm\tfrac13 .

(The last matching uses $\gamma=\tfrac32$ : $\alpha^2=\tfrac14$ and $\nu^2\gamma^2=\tfrac94\nu^2$ .) So the new unknown

u = t^{-\alpha}w = t^{1/2}J_{\pm1/3}\bigl(\tfrac23\mathrm{i}t^{3/2}\bigr)

solves Airy's equation. Since $J_\nu(\mathrm{i}\zeta) = \mathrm{i}^{\nu}I_\nu(\zeta)$ , these are the same as $t^{1/2}I_{\pm1/3}(\tfrac23t^{3/2})$ up to constants, which is the modified form in part 3. (Hassani's Problem 15.41 states the conclusion as " $J_{1/3}(\tfrac23\mathrm{i}t^{3/2})$ and $J_{-1/3}(\tfrac23\mathrm{i}t^{3/2})$ ", without the factor $t^{-\alpha}$ that his own substitution $w=t^{\alpha}u$ introduces; with $\alpha=-\tfrac12$ that factor is $t^{1/2}=\sqrt t$ , and without it the functions do not satisfy the equation. We follow the substitution rather than the printed conclusion; see the remark below.) Replacing $t$ by $-x$ throughout turns the modified Bessel functions into ordinary ones and gives the oscillatory pair. The particular combinations quoted for $\mathrm{Ai}$ and $\mathrm{Bi}$ are fixed by matching the values at $x=0$ : as $\zeta\to0$ , $I_{\mp1/3}(\zeta)\sim(\zeta/2)^{\mp1/3}/\Gamma(1\mp\tfrac13)$ , and feeding $\zeta=\tfrac23x^{3/2}$ into $\tfrac13\sqrt x\,[I_{-1/3}-I_{1/3}]$ gives $3^{-2/3}/\Gamma(2/3)$ in the limit, which is $\mathrm{Ai}(0)$ ; the $x^{1/2}$ factors cancel against $\zeta^{\mp1/3}\propto x^{\mp1/2}$ exactly, which is why a finite limit exists at all.

The asymptotics and the turning point. Write the WKB solutions of (A) (The WKB approximation): for $y''=zy$ they are $z^{-1/4}\exp(\mp\tfrac23z^{3/2})$ . Which fixed combination equals $\mathrm{Ai}$ depends on the sector, because $z^{3/2}$ has a branch point at the origin — this is the Stokes phenomenon, and for (A) the Stokes lines are $\arg z\in\{0,2\pi/3,4\pi/3\}$ . Riley, Hobson and Bence record the two matched forms:

\mathrm{Ai}(z)\ \sim\ \frac{1}{2\sqrt{\pi}z^{1/4}}\exp\!\left(-\tfrac23z^{3/2}\right)\ (z>0), \qquad \mathrm{Ai}(z)\ \sim\ \frac{1}{\sqrt{\pi}(-z)^{1/4}}\sin\!\left(\tfrac23(-z)^{3/2}+\tfrac{\pi}{4}\right)\ (z<0),

their equations (25.52) and (25.53). Note the factor of $2$ between the two amplitudes: a decaying exponential on one side corresponds to a full-amplitude sine on the other, and this factor is the content of the connection formula.

For the turning point, expand $E-V(x) = -V'(a)(x-a)+O\bigl((x-a)^2\bigr)$ about the simple zero $x=a$ , so that $\psi''+k^2\psi=0$ becomes

\psi'' - \frac{2mV'(a)}{\hbar^2}(x-a)\psi = 0 .

Setting $x-a=\mu z$ with $\mu^3 = \hbar^2/[2mV'(a)]$ reduces it to (A) exactly. The physically admissible solution is $\mathrm{Ai}(z)$ , decaying on the forbidden side $z>0$ . Its oscillatory form on the allowed side, continued outwards by WKB, is

\psi \propto \frac{1}{\sqrt{k(x)}}\sin\!\left(\int_x^{a}k(u)\,\mathrm{d}u+\frac{\pi}{4}\right) ,

and the mirror-image argument at the left turning point gives $\psi\propto k^{-1/2}\sin\bigl(\int_{-a}^{x}k+\pi/4\bigr)$ . The two expressions describe the same function, so their phases must agree modulo $\pi$ :

\int_{x}^{a}k+\frac{\pi}{4} + \int_{-a}^{x}k+\frac{\pi}{4} = (n+1)\pi \quad\Longrightarrow\quad \int_{-a}^{a}k(u)\,\mathrm{d}u = \left(n+\tfrac12\right)\pi ,

which is the Bohr–Sommerfeld condition with the half-integer that the two $\pi/4$ 's supply — one from each turning point. $\square$

Remark(Where the hypotheses are used, and two things the sources print differently)Original

Simple zero. The turning-point reduction uses $V'(a)\ne0$ twice: to make $\mu$ finite and non-zero, and to make the linear term the leading one in the expansion. At a quadratic turning point ( $V'(a)=0$ , $V''(a)\ne0$ ) the local equation is the parabolic-cylinder equation instead, and the connection formula changes; this is the situation at the top of a barrier, which is why the WKB tunnelling formula degrades near the barrier maximum.
Separated turning points. The matching argument treats the two turning points independently, which presumes the regions where the Airy approximation is valid do not overlap, i.e. that the classically allowed region is many de Broglie wavelengths wide. That is the usual $n\gg1$ condition; it is a pleasant accident, not a theorem, that the resulting formula is exact for the harmonic oscillator at every $n$ .
$x>0$ in the Bessel form. The identities of part 3 are stated for $x>0$ because $\zeta=\tfrac23x^{3/2}$ and $\sqrt x$ both need a branch choice otherwise. They extend to the cut plane but the oscillatory and exponential pairs are then different branches of the same analytic object.
Riley, Hobson and Bence print two different amplitudes for the same asymptotic. Their equation (25.40), on p. 893, reads $\mathrm{Ai}(x)\sim\bigl(\sqrt{2\pi}\,x^{1/4}\bigr)^{-1}\sin\bigl(\tfrac23\lvert x\rvert^{3/2}+\pi/4\bigr)$ as $x\to-\infty$ , while their equation (25.53), on p. 899, reads $\bigl(\sqrt{\pi}\,(-z)^{1/4}\bigr)^{-1}\sin\bigl(\tfrac23(-z)^{3/2}+\pi/4\bigr)$ . The two differ by $\sqrt2$ , and (25.40) additionally writes $x^{1/4}$ where the argument is negative. We follow (25.53). The choice is not arbitrary: at $x=-20$ , numerical evaluation gives $\mathrm{Ai}(-20)=-0.1764061$ , while (25.53) gives $-0.1766408$ and (25.40) gives $-0.1249039$ . Only the first is an asymptotic approximation to the function. The same $\pi^{-1/2}$ is what makes the connection factor $2$ quoted in the proof come out, and it is forced by the Wronskian in part 2, since the product of the two asymptotic amplitudes must reproduce $1/\pi$ .
Hassani's Problem 15.41 omits a factor. Its substitution is $w=t^{\alpha}u$ with $w$ a Bessel function and $u$ the solution of Airy's equation, and matching forces $\alpha=-\tfrac12$ , so that $u=t^{-\alpha}w=\sqrt t\,w$ ; but the problem's concluding sentence names $J_{\pm1/3}(\tfrac23\mathrm{i}t^{3/2})$ as the solutions of Airy's equation, with no $t^{1/2}$ . Substituting $J_{1/3}(\tfrac23\mathrm{i}t^{3/2})$ alone into $\ddot u = tu$ does not satisfy it. We have used his substitution scheme and kept the factor; the statement as printed is recorded here rather than repeated.

Method(Recognising and using an Airy regime)AdaptedAdapted from Riley, Hobson & Bence (2006), §25.7, the WKB solutions (25.49)–(25.51), the matched forms (25.52)–(25.53) and the worked bound-state problem, pp. 898–901

When a linear second-order problem $y''+Q(x)y=0$ has a coefficient that changes sign:

Locate the zeros of $Q$ and check they are simple ( $Q'(x_0)\ne0$ ). These are the turning points.
Away from them, use WKB. On $Q>0$ the solutions are $Q^{-1/4}\exp\bigl(\pm\mathrm{i}\int\sqrt Q\bigr)$ , oscillatory; on $Q<0$ they are $\lvert Q\rvert^{-1/4}\exp\bigl(\pm\int\sqrt{\lvert Q\rvert}\bigr)$ , one growing and one decaying. Both break down at a turning point, where $Q^{-1/4}$ blows up.
Near a turning point, rescale. With $Q(x)\approx Q'(x_0)(x-x_0)$ , set $x-x_0 = \mu z$ and $\mu = \lvert Q'(x_0)\rvert^{-1/3}$ ; the equation becomes $y''\pm zy=0$ , i.e. Airy's with the sign fixed by the direction of the sign change. The width of the Airy region is $\mu$ , which is the geometric mean of the local wavelength and the scale on which $Q$ varies.
Pick the Airy solution by the boundary condition at infinity. Decay on the forbidden side selects $\mathrm{Ai}$ ; a condition that permits growth, or a problem posed only on the allowed side, may need $\mathrm{Bi}$ or a combination.
Connect. Use the asymptotic pair of part 4 of the statement to match the Airy solution's tails onto the WKB solutions on both sides, and read off the relative amplitude and the extra phase $\pi/4$ .
Count. For a bound-state problem with two turning points, step 5 at each end gives the quantization condition $\oint p\,\mathrm{d}q = 2\pi\hbar(n+\tfrac12)$ . For a one-sided problem — a hard wall at $x=0$ and a turning point at $x=a$ — only one $\pi/4$ appears and the condition is $\int_0^a k = (n+\tfrac34)\pi$ with $n\ge0$ , because the wall contributes a phase $\pi/2$ rather than $\pi/4$ .

A worked consequence. Riley, Hobson and Bence apply steps 1–6 to the symmetric well $V(x)=V_0x^{2s}$ and obtain $E_n = c_s\bigl(n+\tfrac12\bigr)^{2s/(s+1)}$ with $c_s$ independent of $n$ : the scaling comes out of the quantization condition alone, before any integral is evaluated. For $s=1$ this is the exact harmonic spectrum $\propto(n+\tfrac12)$ ; for $s\to\infty$ it tends to $(n+\tfrac12)^2$ , the infinite square well.

Example(Energy levels of a neutron bouncing in gravity)Original

A slow neutron above a horizontal mirror feels $V(z)=mgz$ for $z>0$ and an impenetrable floor at $z=0$ . This is the one-sided problem of step 6, and it is solved exactly by an Airy function. The experiment is real — ultracold neutrons over a polished surface — so the numbers are checkable against measurement as well as against arithmetic.

Step 1: reduce to Airy's equation. The stationary Schrödinger equation is

-\frac{\hbar^2}{2m}\psi'' + mgz\,\psi = E\psi , \qquad \psi(0)=0,\quad \psi(\infty)=0 .

Rearranged, $\psi'' = \dfrac{2m^2g}{\hbar^2}\left(z-\dfrac{E}{mg}\right)\psi$ . Set

\ell = \left(\frac{\hbar^2}{2m^2g}\right)^{1/3}, \qquad \zeta = \frac{z - E/(mg)}{\ell} ,

and the equation is exactly $\mathrm{d}^2\psi/\mathrm{d}\zeta^2 = \zeta\psi$ .

Step 2: apply the boundary conditions. Decay as $z\to\infty$ forces $\psi\propto\mathrm{Ai}(\zeta)$ , discarding $\mathrm{Bi}$ . The floor requires $\psi=0$ at $z=0$ , that is at $\zeta = -E/(mg\ell)$ . So

\mathrm{Ai}\!\left(-\frac{E}{mg\ell}\right)=0 ,

and the spectrum is read straight off the zeros of $\mathrm{Ai}$ . Writing $a_n<0$ for the $n$ th zero,

E_n = -a_n\,mg\ell = \lvert a_n\rvert\left(\frac{m g^2\hbar^2}{2}\right)^{1/3} .

The quantization is therefore exact here, not a WKB approximation: the potential is linear, so the local Airy equation is the global one.

Step 3: numbers. With $m = 1.67492750\times10^{-27}\ \mathrm{kg}$ , $g = 9.80665\ \mathrm{m\,s^{-2}}$ , $\hbar = 1.054571817\times10^{-34}\ \mathrm{J\,s}$ ,

\left(\frac{mg^2\hbar^2}{2}\right)^{1/3} = \bigl(8.9569\times10^{-94}\bigr)^{1/3} = 9.6395\times10^{-32}\ \mathrm{J} = 0.601648\ \mathrm{peV} .

The first three zeros of $\mathrm{Ai}$ are $a_1=-2.3381074$ , $a_2=-4.0879494$ , $a_3=-5.5205598$ , so

E_1 = 2.3381074\times0.601648\ \mathrm{peV} = \mathbf{1.4067\ peV}, \qquad E_2 = 2.4595\ \mathrm{peV}, \qquad E_3 = 3.3214\ \mathrm{peV} .

The ground-state energy is $1.41$ picoelectronvolts, and the corresponding classical turning height is

z_1 = \frac{E_1}{mg} = \frac{2.2538\times10^{-31}}{1.67493\times10^{-27}\times9.80665}\ \mathrm{m} = 13.72\ \mu\mathrm{m} ,

while the Airy length is $\ell = 5.87\ \mu\mathrm{m}$ . These are macroscopic distances, which is why the levels can be resolved by scanning an absorber across a beam.

Step 4: check against WKB. The one-sided quantization rule of step 6 of the method, with the levels numbered from $1$ rather than from $0$ , reads $\int_0^{z_n}k\,\mathrm{d}z=(n-\tfrac14)\pi$ for $n=1,2,\dots$ . With $k=\sqrt{2m(E-mgz)}/\hbar$ the integral is $\tfrac23\lvert a_n\rvert^{3/2}$ in the scaled variable, so WKB predicts

\lvert a_n\rvert \approx \left[\frac{3\pi}{2}\left(n-\tfrac14\right)\right]^{2/3} .

This gives $2.3203$ , $4.0818$ , $5.5172$ against the exact $2.3381$ , $4.0879$ , $5.5206$ — errors of $0.76\%$ , $0.15\%$ and $0.06\%$ , shrinking as $n$ grows, exactly as the "many wavelengths" hypothesis of the remark predicts. All of these values were computed to ten significant figures and cross-checked against independent evaluations of the Airy zeros.

ApplicationOriginal

Light at a caustic. Airy introduced these functions in 1838 to describe the intensity near a rainbow's caustic, where geometrical optics predicts an infinite brightness. Near a fold caustic the two interfering ray paths coalesce, the phase difference has a cubic stationary point, and the diffraction integral is exactly $\int\exp[\mathrm{i}(s^3/3+xs)]\mathrm{d}s$ — the Airy integral (The method of stationary phase). The resulting pattern, bright fringes on the lit side decaying as $x^{-1/4}$ in envelope with phase $\tfrac23\lvert x\rvert^{3/2}$ , and a smooth exponential tail on the dark side, is the supernumerary bows visible inside a bright rainbow. The $\lvert x\rvert^{-1/4}$ and the $\tfrac23\lvert x\rvert^{3/2}$ in part 4 of the statement are the fringe amplitude and the fringe spacing as measured.

Every turning point in quantum mechanics. Wherever a classical particle stops and turns around, the wavefunction is locally an Airy function. The consequences are everywhere: the exponential tail under a barrier (and so the whole theory of tunnelling, including $\alpha$ -decay and the scanning tunnelling microscope); the $\tfrac12$ in the Bohr–Sommerfeld rule, which is the zero-point energy; the shift of highly excited atomic levels in an external field; and the "Airy pattern" of the wavefunction just inside a classical turning point in the semiclassical limit. The quantum bouncer of the example is the cleanest realisation, because there the Airy function is the exact answer rather than a local approximation, and the measured level spacings in ultracold-neutron experiments agree with it.

Dispersive wave tails. Strauss's exercise solves the linearized Korteweg–de Vries equation $u_t+u_{xxx}=0$ by Fourier transform and finds the kernel to be an Airy function with the self-similar argument $(x-y)/t^{1/3}$ : the solution is $u(x,t) = t^{-1/3}\int A\bigl((x-y)/t^{1/3}\bigr)\varphi(y)\,\mathrm{d}y$ up to a constant, and the $t^{1/3}$ is the scaling $x\to kx$ , $t\to k^3t$ that the equation respects. Physically this is the oscillatory wave train that runs ahead of a solitary water wave or an undular bore, and its envelope decays as $t^{-1/3}$ — a decay rate measured in tidal bores and in the wake behind a ship in shallow water. The same kernel governs the leading edge of an optical pulse in a fibre with third-order dispersion.

Numerics. The two representations divide the work cleanly. For $\lvert z\rvert\lesssim5$ the entire series of the definition converge fast and are the method of choice, with the caveat that for $z$ negative and moderately large they suffer catastrophic cancellation: $\mathrm{Ai}(-10)=0.04024$ emerges from terms as large as $8.8\times10^{7}$ , so double precision loses nine digits there. For larger $\lvert z\rvert$ the asymptotic forms of part 4, with their correction series, are both faster and more accurate — at $z=10$ the leading term alone already gives $\mathrm{Ai}(10)=1.1083\times10^{-10}$ against the true $1.1048\times10^{-10}$ , an error of $0.3\%$ . In between, the Bessel representation of part 3 lets a library reuse code it already has for $K_{1/3}$ and $J_{\pm1/3}$ (Modified Bessel functions), which is how several standard implementations are in fact organised.