The Gamma function and related functions

Gamma, Beta, digamma, zeta and error functions.

In this chapter 5

Tick a unit to add it to your PDF basket.

Gamma and related functions

Gamma and related functions

Definition

The Gamma function

Adapted

Add to PDF

Defines $\Gamma$ by Euler's integral on the right half-plane, collects the functional equation $\Gamma(z+1) = z\Gamma(z)$ with the factorial and half-integer values, Gauss's limit formula, and the reflection and duplication formulas, and shows how to reduce an integral or a half-integral argument to a tabulated value.

Prerequisites: Differentiation under the integral sign, Analytic continuation, Multiple integrals and change of variables, Gaussian integrals, Integrals around branch cuts

Definition(Euler's integral)AdaptedAdapted from Arfken, Weber & Harris (2013), §13.1, eq. (13.5), p. 600; Riley, Hobson & Bence (2006), eq. (18.153), p. 636; Hassani (2013), eq. (12.9), p. 378; Folland (1999), §2.3, p. 58; Strauss (2008), §A.5, p. 425

For $z \in \mathbb{C}$ with $\operatorname{Re} z > 0$ , the Gamma function is

\Gamma(z) \;=\; \int_0^{\infty} t^{\,z-1} e^{-t}\,dt , \qquad t^{\,z-1} := e^{(z-1)\ln t} \ \ (t > 0).

The integral converges absolutely exactly on this half-plane. Write $x = \operatorname{Re} z$ , so that $\lvert t^{\,z-1}e^{-t}\rvert = t^{\,x-1}e^{-t}$ . Near $t = 0$ the exponential is bounded between two positive constants and $\int_0^1 t^{\,x-1}dt$ is finite if and only if $x > 0$ ; near $t = \infty$ the factor $e^{-t}$ beats every power, so the tail converges for every $z$ . The same estimate shows that the convergence is uniform on each strip $0 < a \le x \le b < \infty$ .

Two changes of variable give forms that occur as often as the original. Substituting $t = y^2$ ,

\Gamma(z) \;=\; 2\int_0^{\infty} y^{\,2z-1} e^{-y^2}\,dy ,

which connects $\Gamma$ to the Gaussian integrals (Gaussian integrals); and substituting $t = \ln(1/\tau)$ ,

\Gamma(z) \;=\; \int_0^{1} \Bigl(\ln\tfrac1\tau\Bigr)^{z-1} d\tau .

Both are listed by Arfken, Weber and Harris (Arfken, Weber & Harris (2013), eqs. (13.6)-(13.7), p. 600), the first also by Riley, Hobson and Bence (Riley, Hobson & Bence (2006), eq. (18.157), p. 636).

Definition(Properties of

\Gamma

)AdaptedAdapted from Arfken, Weber & Harris (2013), §13.1, eqs. (13.1), (13.2), (13.4), (13.8), (13.23), (13.27)-(13.29), pp. 599-604; Riley, Hobson & Bence (2006), eqs. (18.154)-(18.158), pp. 636-637; Hassani (2013), Box 12.4.1 and eqs. (12.11), (12.13), (12.20), (12.22), pp. 379-381; Folland (1999), §2.3 and Prop. 2.56, pp. 58-59, 80; Strauss (2008), §A.5, eqs. (1)-(4), p. 425

Let $\Gamma$ be given by Euler's integral on $H = \{\operatorname{Re} z > 0\}$ .

Analyticity and continuation. $\Gamma$ is analytic on $H$ and has a unique extension to a function analytic on $\mathbb{C}\setminus\{0,-1,-2,\dots\}$ with a simple pole at each non-positive integer, the residue at $z = -n$ being $(-1)^n/n!$ . (Carried out in Analytic continuation; see the remark below.) From here on $\Gamma$ means that extension.
Functional equation. $\Gamma(z+1) = z\,\Gamma(z)$ at every $z$ where both sides are finite. In particular $\Gamma(1) = 1$ and
$\Gamma(n+1) \;=\; n! \qquad (n = 0, 1, 2, \dots).$
Half-integer values. $\Gamma(\tfrac12) = \sqrt\pi$ , and for $n = 0, 1, 2, \dots$
$\Gamma\!\left(n+\tfrac12\right) \;=\; \frac{(2n)!}{4^{n}\,n!}\,\sqrt\pi .$
Gauss's limit formula. For every $z \notin \{0,-1,-2,\dots\}$ ,
$\Gamma(z) \;=\; \lim_{n\to\infty} \frac{n!\; n^{z}}{z(z+1)(z+2)\cdots(z+n)} .$
Arfken, Weber and Harris (Arfken, Weber & Harris (2013), eq. (13.1), p. 599) take this as their first definition of $\Gamma$ and derive the Euler integral from it; we go the other way.
Reflection formula. For $z \notin \mathbb{Z}$ ,
$\Gamma(z)\,\Gamma(1-z) \;=\; \frac{\pi}{\sin \pi z} , \qquad\text{equivalently}\qquad \Gamma(z)\,\Gamma(-z) \;=\; -\,\frac{\pi}{z \sin \pi z} .$
Legendre's duplication formula. For every $z$ at which both sides are finite,
$\Gamma(z)\,\Gamma\!\left(z+\tfrac12\right) \;=\; 2^{\,1-2z}\sqrt\pi\;\Gamma(2z) .$
No zeros. $\Gamma(z) \ne 0$ for every $z$ in its domain; equivalently $1/\Gamma$ is entire.

Idea of the proofOriginal

Everything follows from two computations. The functional equation is one integration by parts, and iterating it both extends $\Gamma$ to the left and turns it into the factorial on the integers. The second computation is the one that produces $\pi$ : multiply two Euler integrals together and change variables in the plane. Taking the product in the form $\Gamma(z)\Gamma(1-z)$ and separating the total $s+t$ from the ratio $s/t$ makes the double integral factor, and the ratio integral is a standard keyhole-contour integral worth $\pi/\sin\pi z$ ; that is the reflection formula, and $\Gamma(1/2)^2 = \pi$ is its value at $z = 1/2$ .

The duplication formula is not an integral computation at all. Both $\Gamma(z)\Gamma(z+\frac12)$ and $\Gamma(2z)$ are limits of explicit products of $2n+2$ linear factors, and those two products differ only by a factor free of $z$ together with a power $(n/(2n+1))^{2z}$ . So the ratio of the two sides is a constant, which one reads off at $z = 1/2$ . Gauss's limit formula, which supplies those products, is itself the Euler integral with $e^{-t}$ replaced by $(1-t/n)^n$ and then integrated by parts $n$ times.

ProofAdaptedAdapted from Arfken, Weber & Harris (2013), §13.1, eqs. (13.9)-(13.14) and (13.24)-(13.29), pp. 601-604; Riley, Hobson & Bence (2006), §18.12.1, p. 636; Hassani (2013), eqs. (12.11)-(12.13) and (12.20)-(12.22), pp. 379-381; Folland (1999), §2.3, pp. 58-59; Strauss (2008), §A.5, p. 425

(2) The functional equation. Fix $\operatorname{Re} z > 0$ and integrate by parts on $[\varepsilon, T] \subset (0,\infty)$ , differentiating $t^{\,z}$ and integrating $e^{-t}$ :

\int_{\varepsilon}^{T} t^{\,z} e^{-t}\,dt \;=\; \varepsilon^{\,z}e^{-\varepsilon} - T^{\,z}e^{-T} + z\int_{\varepsilon}^{T} t^{\,z-1}e^{-t}\,dt .

The first boundary term is $O(\varepsilon^{\operatorname{Re}z})$ and vanishes as $\varepsilon \to 0$ because $\operatorname{Re} z > 0$ ; the second vanishes as $T \to \infty$ because $e^{-T}$ beats $T^{\operatorname{Re}z}$ . Hence $\Gamma(z+1) = z\Gamma(z)$ on $H$ , and then on all of $\mathbb{C}\setminus\{0,-1,\dots\}$ because two functions analytic there and equal on $H$ are equal (Analytic continuation). Also $\Gamma(1) = \int_0^\infty e^{-t}dt = 1$ , so induction gives $\Gamma(n+1) = n\Gamma(n) = \dots = n!$ .

(3) Half-integer values. In the squared form of the definition put $z = \tfrac12$ :

\Gamma\!\left(\tfrac12\right) \;=\; 2\int_0^{\infty} e^{-y^2}\,dy \;=\; \int_{-\infty}^{\infty} e^{-y^2}\,dy \;=\; \sqrt\pi ,

the last step by the Gaussian integral (Gaussian integrals). Applying (2) $n$ times,

\Gamma\!\left(n+\tfrac12\right) \;=\; \left(n-\tfrac12\right)\left(n-\tfrac32\right)\cdots\tfrac12 \cdot \sqrt\pi \;=\; \frac{(2n-1)(2n-3)\cdots 1}{2^{n}}\,\sqrt\pi .

Multiplying numerator and denominator by $(2n)(2n-2)\cdots 2 = 2^{n} n!$ turns the odd product into $(2n)!$ and gives the stated closed form. (Arfken, Weber and Harris (Arfken, Weber & Harris (2013), eq. (13.28), p. 604) write the same thing with double factorials, $(2n-1)!!\,\sqrt\pi/2^{n}$ .)

(4) Gauss's limit formula. For $\operatorname{Re} z > 0$ and $n \in \mathbb{N}$ set

F_n(z) \;=\; \int_0^{n} \Bigl(1 - \frac{t}{n}\Bigr)^{n} t^{\,z-1}\,dt .

Step 1: $F_n(z) \to \Gamma(z)$ . Since $\ln(1-u) \le -u$ on $[0,1)$ , we have $0 \le (1-t/n)^n \le e^{-t}$ for $0 \le t \le n$ , so the integrands are dominated by $t^{\,x-1}e^{-t}$ , which is integrable on $(0,\infty)$ for $x = \operatorname{Re}z > 0$ . Pointwise $(1-t/n)^n \to e^{-t}$ for each fixed $t > 0$ . Dominated convergence (Dominated convergence) applied to the integrands extended by $0$ beyond $t = n$ gives $F_n(z) \to \Gamma(z)$ .

Step 2: evaluate $F_n(z)$ . Substituting $t = nu$ ,

F_n(z) \;=\; n^{z}\int_0^1 (1-u)^{n} u^{\,z-1}\,du \;=\; n^{z}\,J_n(z), \qquad J_m(w) := \int_0^1 (1-u)^{m} u^{\,w-1}\,du .

Integrate $J_m(w)$ by parts for $m \ge 1$ , differentiating $(1-u)^m$ and integrating $u^{\,w-1}$ to $u^{\,w}/w$ . The boundary term vanishes at $u = 1$ because $m \ge 1$ and at $u = 0$ because $\operatorname{Re} w > 0$ , leaving

J_m(w) \;=\; \frac{m}{w}\,J_{m-1}(w+1) .

This is the one place where $\operatorname{Re} w > 0$ is needed. Iterating $n$ times and finishing with $J_0(w) = \int_0^1 u^{\,w-1}du = 1/w$ ,

J_n(z) \;=\; \frac{n}{z}\cdot\frac{n-1}{z+1}\cdots\frac{1}{z+n-1}\cdot\frac{1}{z+n} \;=\; \frac{n!}{z(z+1)\cdots(z+n)} .

Combining the two steps proves (4) on $H$ . Step 3: all other $z$ . Let $z \notin\{0,-1,-2,\dots\}$ and pick an integer $m \ge 1$ with $\operatorname{Re}(z+m) > 0$ . Writing $P_n(w) = w(w+1)\cdots(w+n)$ , the exact algebraic identity

\frac{n!\,n^{z}}{P_n(z)} \;=\; \frac{n!\,n^{\,z+m}}{P_n(z+m)}\cdot\frac{\prod_{j=n+1}^{\,n+m}(z+j)}{n^{m}\prod_{j=0}^{\,m-1}(z+j)}

holds for every $n$ , because $P_n(z+m)/P_n(z)$ is the quotient of those two products of $m$ factors. As $n \to \infty$ the first factor tends to $\Gamma(z+m)$ by steps 1-2, and $\prod_{j=n+1}^{n+m}(z+j)/n^{m}\to 1$ since each of the $m$ factors is $n(1+O(1/n))$ . So the limit exists and equals $\Gamma(z+m)/\bigl(z(z+1)\cdots(z+m-1)\bigr)$ , which is $\Gamma(z)$ by part (2). (The integral $J_n$ is a Beta integral, and this calculation is the special case $B(z, n+1)$ of The Beta function.)

(5) Reflection. Take $0 < \operatorname{Re} z < 1$ and multiply two Euler integrals:

\Gamma(z)\,\Gamma(1-z) \;=\; \int_0^{\infty}\!\!\int_0^{\infty} s^{\,z-1} t^{-z} e^{-(s+t)}\,ds\,dt .

The double integral converges absolutely, since the modulus of the integrand is $s^{\,x-1}t^{-x}e^{-(s+t)}$ with $x = \operatorname{Re}z$ and its integral is the finite number $\Gamma(x)\Gamma(1-x)$ ; so the iterated integral may be treated as a double integral and the variables changed (Multiple integrals and change of variables). Replace $(s,t)$ by the total and the ratio,

u \;=\; s+t, \qquad v \;=\; s/t, \qquad\text{so}\qquad s \;=\; \frac{uv}{1+v}, \quad t \;=\; \frac{u}{1+v},

a $C^1$ bijection of $(0,\infty)^2$ onto itself with Jacobian determinant $-u/(1+v)^2$ . The two powers collapse,

s^{\,z-1}t^{-z} \;=\; \left(\frac{uv}{1+v}\right)^{z-1}\left(\frac{u}{1+v}\right)^{-z} \;=\; \frac{(1+v)\,v^{\,z-1}}{u} ,

and multiplying by the Jacobian factor $u/(1+v)^2$ leaves $v^{\,z-1}/(1+v)$ , free of $u$ . The integral therefore splits:

\Gamma(z)\,\Gamma(1-z) \;=\; \left(\int_0^{\infty} e^{-u}\,du\right)\left(\int_0^{\infty} \frac{v^{\,z-1}}{1+v}\,dv\right) \;=\; \int_0^{\infty} \frac{v^{\,z-1}}{1+v}\,dv .

For $0 < \operatorname{Re} z < 1$ the remaining integral is the keyhole-contour integral $\int_0^\infty v^{p-1}(v+a)^{-1}dv = \pi a^{\,p-1}/\sin p\pi$ of Integrals around branch cuts at $a = 1$ , $p = z$ , with value $\pi/\sin\pi z$ . Both sides of the resulting identity are analytic on $\mathbb{C}\setminus\mathbb{Z}$ , so it holds there. Finally, replacing $z$ by $-z$ and using $\Gamma(1+z) = z\Gamma(z)$ on the left gives $\Gamma(z)\Gamma(-z) = -\pi/(z\sin\pi z)$ .

(6) Duplication. With $P_n$ as in step 3 above, (4) reads $\Gamma(w) = \lim_n n!\,n^{w}/P_n(w)$ . Pair the factors of $P_n(z)$ and $P_n(z+\tfrac12)$ :

P_n(z)\,P_n\!\left(z+\tfrac12\right) \;=\; \prod_{k=0}^{n}\left(z+k\right)\left(z+k+\tfrac12\right) \;=\; \prod_{k=0}^{n}\frac{(2z+2k)(2z+2k+1)}{4} \;=\; \frac{P_{2n+1}(2z)}{4^{\,n+1}} ,

because the factors $2z+j$ for $j = 0,1,\dots,2n+1$ each occur exactly once. Hence, with $m = 2n+1$ ,

\frac{\bigl(n!\,n^{z}/P_n(z)\bigr)\bigl(n!\,n^{z+1/2}/P_n(z+\tfrac12)\bigr)}{m!\,m^{2z}/P_m(2z)} \;=\; \frac{(n!)^2\,\sqrt n\;4^{\,n+1}}{(2n+1)!}\cdot\left(\frac{n}{2n+1}\right)^{2z} .

Let $n \to \infty$ . The left side tends to $\Gamma(z)\Gamma(z+\tfrac12)/\Gamma(2z)$ by (4), a finite non-zero number when $\operatorname{Re}z > 0$ . On the right, $(n/(2n+1))^{2z} \to 2^{-2z} \ne 0$ . A quotient of two convergent sequences with non-zero limits converges, so the remaining factor $c_n = (n!)^2\sqrt n\,4^{\,n+1}/(2n+1)!$ converges as well, say to $c$ ; and $c_n$ does not involve $z$ , so $c$ is one constant for all $z$ . Hence

\frac{\Gamma(z)\,\Gamma(z+\tfrac12)}{\Gamma(2z)} \;=\; c\,2^{-2z} \qquad (\operatorname{Re} z > 0).

Evaluate at $z = \tfrac12$ : the left side is $\Gamma(\tfrac12)\Gamma(1)/\Gamma(1) = \sqrt\pi$ and the right side is $c/2$ , so $c = 2\sqrt\pi$ and the displayed identity is $\Gamma(z)\Gamma(z+\tfrac12) = 2^{\,1-2z}\sqrt\pi\,\Gamma(2z)$ . Both sides are meromorphic on $\mathbb{C}$ , so this holds wherever they are finite. (Numerically, $c_n = 3.41855$ , $3.53168$ , $3.54477$ at $n = 10$ , $10^2$ , $10^4$ , against $2\sqrt\pi = 3.5449077$ .)

(7) No zeros. If $\Gamma(z_0) = 0$ for some $z_0$ in the domain, then $z_0 \notin \mathbb{Z}$ would force the left side of the reflection formula to vanish while $\pi/\sin\pi z_0$ is finite and non-zero, which is impossible; and at $z_0 = n \ge 1$ we have $\Gamma(n) = (n-1)! \ne 0$ by (2). So $\Gamma$ never vanishes, and $1/\Gamma$ , which is analytic off the poles and takes the value $0$ at each of them, is entire. $\square$

Remark(Where the hypotheses are used)Original

$\operatorname{Re} z > 0$ in the definition is exactly the integrability condition at $t = 0$ , as the estimate in the definition block shows; for $\operatorname{Re} z \le 0$ the integral diverges outright rather than converging conditionally, so the continuation cannot be obtained by reinterpreting it. Folland (Folland (1999), §2.3, p. 58) makes this the point of departure: the integral defines $\Gamma$ on the half-plane, and the functional equation then does the continuing.
The same condition is what makes the integration by parts work twice over: it kills the boundary term $\varepsilon^{\,z}e^{-\varepsilon}$ in part (2) and the term at $u = 0$ in the recursion for $J_m$ in part (4).
Dominated convergence in part (4) needs the inequality $(1-t/n)^n \le e^{-t}$ , not merely the pointwise limit; without a dominating function the interchange of limit and integral is not automatic, because the interval of integration also grows.
$0 < \operatorname{Re} z < 1$ in part (5) is needed twice: for absolute convergence of the double integral (at $x \le 0$ the $s$ -integral diverges, at $x \ge 1$ the $t$ -integral does), and for the keyhole evaluation of $\int_0^\infty v^{\,z-1}/(1+v)\,dv$ . Only after both sides are known to be analytic does continuation extend the identity to all non-integral $z$ . Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §18.12.2, p. 639) are explicit that the elementary argument gives the formula on this strip and that more is needed elsewhere.
At integers the reflection formula must be read as an identity of meromorphic functions. At $z = n \ge 1$ , $\Gamma(1-z)$ has a pole and $\sin \pi z = 0$ , and the two singularities match; the useful content at those points is part (7) together with the pole structure of part (1).
Duplication at $2z \in \{0,-1,-2,\dots\}$ likewise compares two functions with poles at the same places, which is why we proved it first on $\operatorname{Re} z > 0$ and then continued.

Remark(Conventions, and what is proved where)Original

The analyticity of $\Gamma$ on the right half-plane, the continuation by $\Gamma(z) = \Gamma(z+1)/z$ , the simple poles at $0,-1,-2,\dots$ and the residues $(-1)^n/n!$ are proved in detail in Analytic continuation, where $\Gamma$ serves as the worked example of continuation by a functional relation; part (1) above only records the result, and this unit's own work is the functional equation, the limit formula, and the reflection and duplication formulas. Worth noting about that proof: the Leibniz rule for differentiating under the integral sign (Differentiation under the integral sign) does not apply to Euler's integral, because it is stated for a bounded interval and a real parameter while here the interval is unbounded and the integrand unbounded at $t=0$ ; the analyticity is obtained instead by truncating the integral and applying Morera's theorem to the limit. The Leibniz rule does apply to each truncation $\int_{1/n}^{n}$ , which is how the differentiated forms $\Gamma'(x) = \int_0^\infty t^{\,x-1}e^{-t}\ln t\,dt$ and its higher analogues are justified for real $x > 0$ .

Three conventions worth keeping straight.

Argument shifts. Many texts, Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §13.1-13.2, pp. 599-611) among them, state results for $\Gamma(z+1)$ rather than $\Gamma(z)$ , since $\Gamma(z+1) = z!$ when $z$ is a non-negative integer. Their duplication formula, $\Gamma(1+z)\Gamma(z+\frac12) = 2^{-2z}\sqrt\pi\,\Gamma(2z+1)$ , becomes ours on applying the functional equation to $\Gamma(1+z)$ and $\Gamma(2z+1)$ .
The factorial notation. $z!$ for $\Gamma(z+1)$ at non-integral $z$ is used by Riley, Hobson and Bence (Riley, Hobson & Bence (2006), eqs. (18.155)-(18.156), p. 636) and avoided by Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §13.1, p. 606). We write $\Gamma$ throughout and reserve $n!$ for integers.
Which property is the definition. Euler's integral (here, and in Riley, Hobson and Bence (Riley, Hobson & Bence (2006), p. 636), Hassani (Hassani (2013), p. 378), Folland (Folland (1999), p. 58) and Strauss (Strauss (2008), p. 425)), Gauss's limit (Arfken, Weber and Harris (Arfken, Weber & Harris (2013), eq. (13.1), p. 599)), and the Weierstrass product $1/\Gamma(z) = z e^{\gamma z}\prod_{k\ge1}(1+z/k)e^{-z/k}$ (Arfken, Weber and Harris (Arfken, Weber & Harris (2013), eq. (13.15), p. 602), Hassani (Hassani (2013), eq. (12.14), p. 379)) all define the same function; part (4) is the bridge between the first two, and taking reciprocals in part (4) and inserting $n^{-z} = e^{-z\ln n}$ gives the third, with the Euler-Mascheroni constant $\gamma$ appearing as $\lim_n(\sum_{k\le n}k^{-1} - \ln n)$ . The product form is the natural starting point for Digamma and polygamma functions.

Method(Reducing a quantity to

\Gamma

)Original

Powers times a stretched exponential. To evaluate $\int_0^\infty x^{b} e^{-c x^{k}}\,dx$ with $c > 0$ , $k > 0$ and $b > -1$ , substitute $u = c x^{k}$ . Then $x = (u/c)^{1/k}$ , $dx = (u/c)^{1/k}k^{-1}u^{-1}du$ and
$\int_0^{\infty} x^{b} e^{-c x^{k}}\,dx \;=\; \frac{1}{k\,c^{(b+1)/k}}\;\Gamma\!\left(\frac{b+1}{k}\right).$
Check the exponent of $c$ by dimensions: $x$ scales as $c^{-1/k}$ , so the integral scales as $c^{-(b+1)/k}$ .
Reduce the argument to $[1,2]$ . Use $\Gamma(z) = \Gamma(z+1)/z$ to step right until the argument lies in $[1,2]$ , where $\Gamma$ is between $0.8856$ and $1$ ; or $\Gamma(z+1) = z\Gamma(z)$ to step left. This is both the numerically stable way to evaluate $\Gamma$ and the way to read a half-integer or negative value off $\Gamma(\tfrac12) = \sqrt\pi$ .
Use reflection when two arguments add to an integer. A product $\Gamma(p)\Gamma(q)$ with $p + q = 1$ is $\pi/\sin\pi p$ with no further work; with $p+q = m \in \mathbb{Z}$ , first shift by the functional equation until the sum is $1$ .
Use duplication when two arguments differ by $\tfrac12$ . It converts $\Gamma(z)\Gamma(z+\frac12)$ into a single $\Gamma(2z)$ , and read backwards it splits $\Gamma$ of an even multiple into two factors — the move that generates half-integer values and that underlies the $\Gamma$ -quotients in the coefficients of Bessel and hypergeometric series.
For a large or complex argument, go asymptotic. Use Stirling's series (Stirling's formula) on $\ln\Gamma$ , after shifting the argument into the right half-plane with reflection if necessary; summing $\ln\Gamma$ is also how large factorials are handled in statistical mechanics.

Example(Half-integer values and a Gaussian moment)Original

Goal: evaluate $\displaystyle\int_0^{\infty} x^{6} e^{-x^{2}}\,dx$ and the two values $\Gamma(\tfrac72)$ , $\Gamma(-\tfrac52)$ .

Reduce the integral. Step 1 of the method with $b = 6$ , $c = 1$ , $k = 2$ gives
$\int_0^{\infty} x^{6}e^{-x^{2}}\,dx \;=\; \tfrac12\,\Gamma\!\left(\tfrac72\right).$
Step down to $\Gamma(\tfrac12)$ . By the functional equation three times,
$\Gamma\!\left(\tfrac72\right) \;=\; \tfrac52\cdot\tfrac32\cdot\tfrac12\,\Gamma\!\left(\tfrac12\right) \;=\; \tfrac{15}{8}\sqrt\pi \;=\; 3.3233510 .$
The closed form of part (3) agrees: $n = 3$ gives $(2n)!\sqrt\pi/(4^n n!) = 720\sqrt\pi/(64\cdot 6) = \tfrac{15}{8}\sqrt\pi$ .
The integral. Therefore $\int_0^\infty x^6 e^{-x^2}dx = \tfrac{15}{16}\sqrt\pi = 1.6616755$ . Numerical quadrature of the integral returns $1.6616755$ .
A negative argument. Step the functional equation leftwards from $\Gamma(\tfrac12)$ , i.e. use $\Gamma(z) = \Gamma(z+1)/z$ at $z = -\tfrac12, -\tfrac32, -\tfrac52$ :
$\Gamma\!\left(-\tfrac12\right) = \frac{\sqrt\pi}{-\tfrac12} = -2\sqrt\pi, \qquad \Gamma\!\left(-\tfrac32\right) = \frac{-2\sqrt\pi}{-\tfrac32} = \tfrac43\sqrt\pi, \qquad \Gamma\!\left(-\tfrac52\right) = \frac{\tfrac43\sqrt\pi}{-\tfrac52} = -\tfrac{8}{15}\sqrt\pi .$
So $\Gamma(-\tfrac52) = -0.9453087$ , and the sign alternates with each step, as it must: between consecutive poles $\Gamma$ keeps one sign, and the residues $(-1)^n/n!$ alternate.
Cross-check by reflection. With $z = \tfrac72$ , part (5) gives $\Gamma(\tfrac72)\Gamma(-\tfrac52) = \pi/\sin(7\pi/2) = -\pi$ . Indeed $3.3233510 \times (-0.9453087) = -3.1415927$ .

Example(An algebraic integral by reflection)Original

Goal: evaluate $\displaystyle I = \int_0^{\infty} \frac{dx}{x^{1/3}(1+x)}$ in closed form.

Recognise the shape. The proof of part (5) showed that for $0 < \operatorname{Re}z < 1$
$\int_0^{\infty}\frac{v^{\,z-1}}{1+v}\,dv \;=\; \Gamma(z)\,\Gamma(1-z) .$
Here $z - 1 = -\tfrac13$ , so $z = \tfrac23$ , which lies in the strip.
Apply reflection.
$I \;=\; \Gamma\!\left(\tfrac23\right)\Gamma\!\left(\tfrac13\right) \;=\; \frac{\pi}{\sin(2\pi/3)} \;=\; \frac{\pi}{\sqrt3/2} \;=\; \frac{2\pi}{\sqrt3} \;=\; 3.6275987 .$
Two independent checks. The tabulated values $\Gamma(\tfrac13) = 2.6789385$ and $\Gamma(\tfrac23) = 1.3541179$ multiply to $3.6275987$ ; and numerical quadrature of $I$ gives $3.6275987$ .
A duplication check at a different point. Reflection at $z = \tfrac14$ gives $\Gamma(\tfrac14)\Gamma(\tfrac34) = \pi/\sin(\pi/4) = \pi\sqrt2 = 4.4428829$ . Duplication at $z = \tfrac14$ gives the same product by a different route, since $\Gamma(\tfrac14)\Gamma(\tfrac34) = 2^{1/2}\sqrt\pi\,\Gamma(\tfrac12) = \sqrt2\,\pi$ . The tabulated values $\Gamma(\tfrac14) = 3.6256099$ and $\Gamma(\tfrac34) = 1.2254167$ confirm it.

Example(Mean and spread of a Weibull lifetime)Original

A batch of pump bearings has time to failure $T$ with density $f(t) = (\beta/\tau)(t/\tau)^{\beta-1}e^{-(t/\tau)^{\beta}}$ on $t > 0$ , with shape $\beta = \tfrac32$ and scale $\tau = 10$ h. Find the mean life and the standard deviation.

Moments. With $u = (t/\tau)^{\beta}$ , so $t = \tau u^{1/\beta}$ and $f(t)\,dt = e^{-u}\,du$ ,
$\mathbb{E}[T^{\,r}] \;=\; \int_0^{\infty} t^{\,r} f(t)\,dt \;=\; \tau^{\,r}\int_0^{\infty} u^{\,r/\beta}e^{-u}\,du \;=\; \tau^{\,r}\,\Gamma\!\left(1 + \frac{r}{\beta}\right).$
The substitution is step 1 of the method in disguise; note that $\int_0^\infty u^{r/\beta}e^{-u}du$ is $\Gamma(1+r/\beta)$ , not $\Gamma(r/\beta)$ , because the exponent of $u$ is $r/\beta$ and Euler's integral carries $z-1$ .
Mean. $r = 1$ , $r/\beta = \tfrac23$ :
$\mathbb{E}[T] \;=\; 10\,\Gamma\!\left(\tfrac53\right) \;=\; 10\cdot\tfrac23\,\Gamma\!\left(\tfrac23\right) \;=\; \tfrac{20}{3}(1.3541179) \;=\; 9.0274529\ \mathrm{h},$
using the functional equation to express $\Gamma(\tfrac53)$ through the value of $\Gamma$ at $\tfrac23$ already computed above.
Second moment. $r = 2$ , $r/\beta = \tfrac43$ , and two applications of the functional equation give $\Gamma(\tfrac73) = \tfrac43\cdot\tfrac13\,\Gamma(\tfrac13) = \tfrac49(2.6789385) = 1.1906393$ . So $\mathbb{E}[T^2] = 100(1.1906393) = 119.06393\ \mathrm{h^2}$ .
Spread. $\operatorname{Var}T = \mathbb{E}[T^{2}] - (\mathbb{E}[T])^{2}$ . Carrying one more figure than the rounded values quoted above, $\mathbb{E}[T] = 9.02745293$ h and $\mathbb{E}[T^{2}] = 119.0639349\ \mathrm{h^2}$ , so
$\operatorname{Var}T \;=\; 119.0639349 - 81.4949064 \;=\; 37.5690285\ \mathrm{h^2},$
and the standard deviation is $6.1293579$ h, about $68\%$ of the mean. Numerical integration of $\int t f(t)dt$ and $\int t^2 f(t)dt$ reproduces steps 2 and 3 to all digits shown.
Sanity check. At $\beta = 1$ the density is exponential with mean $\tau$ and standard deviation $\tau$ , and the formulas give $\tau\Gamma(2) = \tau$ and $\tau\sqrt{\Gamma(3)-\Gamma(2)^2} = \tau$ , as they must.

ApplicationAdaptedAdapted from Folland (1999), §2.7, Props. 2.54-2.56 and Cor. 2.55, pp. 79-80; Strauss (2008), §10.5, eq. (4), p. 283; Arfken, Weber & Harris (2013), Example 13.1.1, pp. 606-607

Volumes in many dimensions. The cleanest appearance of $\Gamma$ in analysis is the measure of the unit sphere and ball in $\mathbb{R}^n$ . Integrating $e^{-\lvert x\rvert^2}$ over $\mathbb{R}^n$ in two ways — as the $n$ -th power of the one-dimensional Gaussian integral, and in polar coordinates, where the radial part is $\int_0^\infty r^{\,n-1}e^{-r^2}dr = \tfrac12\Gamma(n/2)$ by the squared form of the definition — gives

\sigma(S^{n-1}) \;=\; \frac{2\pi^{n/2}}{\Gamma(n/2)} , \qquad \operatorname{vol}(B^{n}) \;=\; \frac{\sigma(S^{n-1})}{n} \;=\; \frac{\pi^{n/2}}{\Gamma\!\left(\tfrac n2 + 1\right)} ,

the second by the functional equation, as Folland (Folland (1999), Prop. 2.54 and Cor. 2.55, pp. 79-80) sets it out. At $n = 3$ these are $4\pi = 12.566371$ and $\tfrac43\pi = 4.1887902$ . Because $\Gamma$ grows faster than any power, the unit ball's volume collapses: it peaks at $n = 5$ with $\pi^{5/2}/\Gamma(7/2) = 5.2637890$ and has fallen to $3.2985089$ by $n = 9$ . That collapse is the reason a Monte Carlo sampler that encloses a high-dimensional ball in its circumscribed cube almost never lands inside, and the reason phase-space volumes in statistical mechanics are computed through $\Gamma$ rather than geometrically.

Thermal averages. Classical statistical mechanics weights a state of energy $E$ by $e^{-\beta E}$ with $\beta = 1/k_{\mathrm B}T$ . For a structureless ideal gas the density of states goes as $E^{1/2}$ , so normalisation and the mean energy are both integrals of the form of step 1 of the method, and their ratio is a ratio of $\Gamma$ 's:

\langle E\rangle \;=\; \frac{\int_0^\infty E^{3/2}e^{-\beta E}\,dE}{\int_0^\infty E^{1/2}e^{-\beta E}\,dE} \;=\; \frac{\Gamma(5/2)}{\beta\,\Gamma(3/2)} \;=\; \frac{3}{2}k_{\mathrm B}T ,

since $\Gamma(5/2)/\Gamma(3/2) = 3/2$ by the functional equation; this is Arfken, Weber and Harris (Arfken, Weber & Harris (2013), Example 13.1.1, pp. 606-607)'s route to the equipartition result. The same reduction gives every moment of the Maxwell speed distribution $f(v) \propto v^{2}e^{-av^{2}}$ , $a = m/(2k_{\mathrm B}T)$ : applying step 1 of the method to numerator and denominator,

\langle v^{n}\rangle \;=\; \frac{\int_0^\infty v^{\,n+2}e^{-av^{2}}dv}{\int_0^\infty v^{2}e^{-av^{2}}dv} \;=\; \left(\frac{2k_{\mathrm B}T}{m}\right)^{n/2}\frac{\Gamma\!\left(\frac{n+3}{2}\right)}{\Gamma\!\left(\frac32\right)} .

The power of $k_{\mathrm B}T/m$ carries the dimensions of $v^{n}$ and the $\Gamma$ -quotient is the pure number; at $n = 1$ it returns the mean speed $\sqrt{8k_{\mathrm B}T/\pi m}$ .

Normalising special functions. $1/\Gamma$ , entire by part (7), is the standard normaliser in series solutions of differential equations: the Bessel function of order $s$ is

J_s(z) \;=\; \sum_{j=0}^{\infty}\frac{(-1)^{j}}{\Gamma(j+1)\,\Gamma(j+s+1)}\left(\frac z2\right)^{2j+s} ,

in which $\Gamma$ replaces the factorials of the integer-order case (Strauss (Strauss (2008), §10.5, eq. (4), p. 283), whose series this is).

Two readings of that series should be kept apart. Strauss (Strauss (2008), §10.5, p. 283) treats negative integral $s$ as excluded: he states that the series solves Bessel's equation for every real $s$ except a negative integer, on the ground that $\Gamma$ is infinite at the negative integers, so that the series "is not defined" there. Our own reading, written as above with $1/\Gamma$ rather than with $\Gamma$ in a denominator, is that the coefficients remain defined at $s = -m$ : the first $m$ of them carry $1/\Gamma$ of a non-positive integer, which is $0$ by part (7), so those terms drop out and what is left sums to $(-1)^{m}J_{m}$ — a solution, but not a new one. The restriction Strauss places on the order is therefore a restriction on getting a second solution, which is the same conclusion Bessel functions of the first kind reaches; nothing above is claimed on his authority beyond the displayed series.

The same role is played in the hypergeometric series (The hypergeometric function) and in the Gamma and Beta probability densities (The Beta function).

Theorem

The Beta function

Adapted

Add to PDF

The Beta integral $B(p,q) = \int_0^1 t^{p-1}(1-t)^{q-1}dt$ equals $\Gamma(p)\Gamma(q)/\Gamma(p+q)$ , proved by writing the product of two Gamma integrals in polar coordinates; the trigonometric and $(1+u)^{-p-q}$ forms that follow give closed values for powers of sine and cosine and for a large family of algebraic integrals.

Prerequisites: The Gamma function, Multiple integrals and change of variables

Definition(Euler's integral of the first kind)AdaptedAdapted from Riley, Hobson & Bence (2006), eq. (18.162), p. 638; Hassani (2013), eq. (12.16), p. 380; Arfken, Weber & Harris (2013), eq. (13.45), p. 617; Grimmett & Stirzaker (2020), §4.4 (4.4.8), p. 108

For complex $p$ and $q$ with $\operatorname{Re} p > 0$ and $\operatorname{Re} q > 0$ , the Beta function is

B(p,q) \;=\; \int_0^1 t^{\,p-1}(1-t)^{\,q-1}\,dt .

The integral converges absolutely under exactly those two conditions. Writing $a = \operatorname{Re}p$ and $b = \operatorname{Re}q$ , the modulus of the integrand is $t^{\,a-1}(1-t)^{\,b-1}$ ; near $t = 0$ the second factor is bounded away from $0$ and $\infty$ , so integrability there needs $a > 0$ , and near $t = 1$ the same argument on the first factor needs $b > 0$ .

Hassani (Hassani (2013), §12.4, p. 380) calls this Euler's integral of the first kind, the Gamma integral being the second. Arfken, Weber and Harris (Arfken, Weber & Harris (2013), eq. (13.45), p. 617) go the other way and define $B(p,q)$ as $\Gamma(p)\Gamma(q)/\Gamma(p+q)$ , then prove that the integral above computes it; the content is the same, and part 2 of the statement is the bridge.

Theorem(The Beta function and its forms)AdaptedAdapted from Riley, Hobson & Bence (2006), eqs. (18.163)-(18.165), pp. 638-639; Arfken, Weber & Harris (2013), eqs. (13.45)-(13.51) and Exercise 13.3.1, pp. 617-619; Hassani (2013), eqs. (12.17)-(12.19), p. 380

Let $\operatorname{Re} p > 0$ and $\operatorname{Re} q > 0$ .

Symmetry. $B(p,q) = B(q,p)$ .
Relation to $\Gamma$ .
$B(p,q) \;=\; \frac{\Gamma(p)\,\Gamma(q)}{\Gamma(p+q)} .$
Trigonometric form.
$B(p,q) \;=\; 2\int_0^{\pi/2} \sin^{2p-1}\!\theta\;\cos^{2q-1}\!\theta\;d\theta .$
Form on the half-line.
$B(p,q) \;=\; \int_0^{\infty}\frac{u^{\,p-1}}{(1+u)^{\,p+q}}\,du \;=\; \int_1^{\infty} s^{-p-q}\,(s-1)^{\,q-1}\,ds .$
Recurrences.
$B(p,q) \;=\; B(p+1,q) + B(p,q+1), \qquad B(p,q+1) \;=\; \frac{q}{p+q}\,B(p,q) .$
Special values and extension. $B(p,1) = 1/p$ and $B(\tfrac12,\tfrac12) = \pi$ . Since the right-hand side of part 2 is a quotient of functions meromorphic on $\mathbb{C}$ , it extends $B$ to all $(p,q)$ with $p, q \notin \{0,-1,-2,\dots\}$ ; in particular, for $0 < \operatorname{Re}p < 1$ the reflection formula of The Gamma function gives
$B(p, 1-p) \;=\; \Gamma(p)\,\Gamma(1-p) \;=\; \frac{\pi}{\sin \pi p} .$

Idea of the proofOriginal

Part 2 is the whole unit, and it is proved by raising the dimension rather than by evaluating anything. A single Gamma integral, written with $e^{-y^2}$ instead of $e^{-t}$ , is an integral over a half-line; the product of two such integrals is an integral over a quadrant of the plane. In polar coordinates that quadrant splits: the radial variable carries all the decay and reassembles into one more Gamma integral, of argument $p+q$ because the powers add, while the angular variable carries all the ratio information and leaves behind an integral in $\sin$ and $\cos$ that depends on $p$ and $q$ separately. Recognising that angular integral as $B(p,q)$ — which is what the substitution $t = \sin^2\theta$ does — finishes the proof.

The other forms are one substitution each, chosen by what the integrand must look like: $t = \sin^2\theta$ sends $[0,1]$ to a quarter period and turns $1-t$ into $\cos^2\theta$ ; $t = u/(1+u)$ sends $[0,1]$ to $[0,\infty)$ and turns $1-t$ into $1/(1+u)$ . The recurrences need no integration at all: the first is the algebraic identity $t^{p-1}(1-t)^{q-1} = t^{p}(1-t)^{q-1} + t^{p-1}(1-t)^{q}$ under the integral sign, and the second is the functional equation of $\Gamma$ applied to part 2.

ProofAdaptedAdapted from Riley, Hobson & Bence (2006), §18.12.2, pp. 638-639; Arfken, Weber & Harris (2013), §13.3, eqs. (13.46)-(13.51), pp. 617-618; Hassani (2013), eqs. (12.17)-(12.19), p. 380

(1) Symmetry. The substitution $t = 1-\tau$ maps $[0,1]$ to itself, reverses orientation and exchanges the two factors, so $B(p,q) = \int_0^1 \tau^{\,q-1}(1-\tau)^{\,p-1}d\tau = B(q,p)$ .

(3) Trigonometric form. Put $t = \sin^2\theta$ for $\theta \in [0,\pi/2]$ , a $C^1$ increasing bijection onto $[0,1]$ with $dt = 2\sin\theta\cos\theta\,d\theta$ and $1 - t = \cos^2\theta$ . Then $t^{\,p-1}(1-t)^{\,q-1}dt = 2\sin^{2p-1}\theta\cos^{2q-1}\theta\,d\theta$ , which is part 3.

(4) Form on the half-line. Put $t = u/(1+u)$ for $u \in [0,\infty)$ , again an increasing bijection onto $[0,1)$ , with $1-t = 1/(1+u)$ and $dt = (1+u)^{-2}du$ . Then

t^{\,p-1}(1-t)^{\,q-1}\,dt \;=\; \frac{u^{\,p-1}}{(1+u)^{\,p-1}}\cdot\frac{1}{(1+u)^{\,q-1}}\cdot\frac{du}{(1+u)^{2}} \;=\; \frac{u^{\,p-1}\,du}{(1+u)^{\,p+q}} .

For the second expression in part 4, substitute $s = 1+v$ with $v \in [0,\infty)$ : it becomes $\int_0^{\infty} v^{\,q-1}(1+v)^{-p-q}\,dv$ , which is the first expression with $p$ and $q$ interchanged, hence equal to it by part 1. That is the form in which Hassani (Hassani (2013), eq. (12.17), p. 380) records it, reached there by $t \mapsto 1/t$ in the definition.

(2) Relation to $\Gamma$ . Use the squared form $\Gamma(w) = 2\int_0^\infty y^{\,2w-1}e^{-y^2}dy$ of Euler's integral (The Gamma function) for both factors:

\Gamma(p)\,\Gamma(q) \;=\; 4\int_0^{\infty}\!\!\int_0^{\infty} x^{\,2p-1} y^{\,2q-1} e^{-(x^2+y^2)}\,dx\,dy .

The double integral converges absolutely: replacing $p$ and $q$ by $a = \operatorname{Re}p$ and $b = \operatorname{Re}q$ turns the integrand into its modulus, and the same display then evaluates to the finite number $\Gamma(a)\Gamma(b)$ . So the iterated integral may be treated as a double integral and the variables changed (Multiple integrals and change of variables). Pass to plane polar coordinates in the open first quadrant,

x \;=\; \rho\sin\varphi, \qquad y \;=\; \rho\cos\varphi, \qquad (\rho,\varphi) \in (0,\infty)\times(0,\tfrac\pi2), \qquad dx\,dy \;=\; \rho\,d\rho\,d\varphi .

Then $x^2 + y^2 = \rho^2$ , and the powers separate as $x^{\,2p-1}y^{\,2q-1} = \rho^{\,2p+2q-2}\sin^{2p-1}\varphi\,\cos^{2q-1}\varphi$ . With the extra $\rho$ from the area element the radial exponent is $2(p+q)-1$ , so

\Gamma(p)\Gamma(q) \;=\; \underbrace{\left(2\int_0^{\pi/2}\sin^{2p-1}\varphi\,\cos^{2q-1}\varphi\;d\varphi\right)}_{=\;B(p,q)\ \text{by part 3}} \;\cdot\; \underbrace{\left(2\int_0^{\infty}\rho^{\,2(p+q)-1}e^{-\rho^2}\,d\rho\right)}_{=\;\Gamma(p+q)}

-- the first factor by part 3, the second by the squared form of Euler's integral at argument $p+q$ , which converges because $\operatorname{Re}(p+q) > 0$ . Since $\Gamma(p+q) \ne 0$ (The Gamma function, part 7), dividing gives part 2.

(5) Recurrences. Inside the integral, $t^{\,p-1}(1-t)^{\,q-1} = \bigl(t + (1-t)\bigr)t^{\,p-1}(1-t)^{\,q-1} = t^{\,p}(1-t)^{\,q-1} + t^{\,p-1}(1-t)^{\,q}$ ; integrating gives the first recurrence, with both new integrals convergent because the exponents only increased. For the second, part 2 and the functional equation $\Gamma(q+1) = q\Gamma(q)$ , $\Gamma(p+q+1) = (p+q)\Gamma(p+q)$ give

B(p,q+1) \;=\; \frac{\Gamma(p)\,\Gamma(q+1)}{\Gamma(p+q+1)} \;=\; \frac{q}{p+q}\cdot\frac{\Gamma(p)\Gamma(q)}{\Gamma(p+q)} \;=\; \frac{q}{p+q}\,B(p,q) .

(6) Special values. $B(p,1) = \int_0^1 t^{\,p-1}dt = 1/p$ , which also follows from part 2 since $\Gamma(1) = 1$ and $\Gamma(p+1) = p\Gamma(p)$ . For $p = q = \tfrac12$ , part 3 gives $B(\tfrac12,\tfrac12) = 2\int_0^{\pi/2}d\theta = \pi$ , in agreement with part 2 and $\Gamma(\tfrac12)^2 = \pi$ . The extension and the value of $B(p,1-p)$ are immediate from part 2 and the reflection formula. $\square$

Remark(Where the hypotheses are used)Original

$\operatorname{Re}p > 0$ and $\operatorname{Re}q > 0$ are two separate endpoint conditions, one per end of $[0,1]$ , and each is needed on its own. Dropping only the first leaves an integral that diverges at $t = 0$ and converges at $t = 1$ . In the trigonometric form of part 3 the same two conditions appear as integrability at $\theta = 0$ and at $\theta = \pi/2$ ; in the half-line form of part 4, as integrability at $u = 0$ and at $u = \infty$ — there the decay at infinity is $u^{-q-1}$ , so it is $\operatorname{Re}q > 0$ that controls the tail.
Absolute convergence is what licenses the polar change of variables in part 2. The integrand is complex-valued whenever $p$ or $q$ is, so the interchange is not covered by any monotone argument; the bound by $x^{\,a-1}y^{\,b-1}e^{-(x^2+y^2)}$ with $a,b > 0$ is.
The identity $\Gamma(p)\Gamma(q) = B(p,q)\Gamma(p+q)$ is only divided by $\Gamma(p+q)$ at the end, and that step uses the fact that $\Gamma$ has no zeros. Written as a product it also makes sense in the limiting cases, which is how part 2 extends to arguments where $\Gamma$ has poles: $B$ inherits those poles from the numerator and zeros from the denominator.
The extension in part 6 is by definition, not by continuation of the integral. For $\operatorname{Re}p \le 0$ the Beta integral diverges; what is meromorphic in $p$ and $q$ is the quotient of Gamma functions, and that quotient is what $B$ means outside the quadrant of convergence.

Remark(A difference of convention between our sources)Original

The two forms

B(p,q) \;=\; 2\int_0^{\pi/2}\sin^{2p-1}\!\theta\,\cos^{2q-1}\!\theta\,d\theta \qquad\text{and}\qquad B(p,q) \;=\; 2\int_0^{\pi/2}\cos^{2p-1}\!\theta\,\sin^{2q-1}\!\theta\,d\theta

are both correct and both in print: the first in Riley, Hobson and Bence (Riley, Hobson & Bence (2006), eq. (18.164), p. 638) and Hassani (Hassani (2013), eq. (12.18), p. 380), the second in Arfken, Weber and Harris (Arfken, Weber & Harris (2013), eq. (13.47), p. 617). They differ by which of $\sin$ and $\cos$ carries which parameter, and they agree because $B$ is symmetric (part 1) — equivalently, because $\theta \mapsto \pi/2 - \theta$ swaps $\sin$ and $\cos$ . We use the first throughout, so that $t = \sin^2\theta$ matches the order of the arguments in the definition. Only the pairing of exponent with parameter is at stake, and a formula taken from one source and used with the other's convention still gives the right number, because the two results are equal; but mixing them inside one calculation makes the bookkeeping unreadable, so we fix one.

Method(Evaluating an integral as a Beta function)Original

Powers of sine and cosine. For $\int_0^{\pi/2}\sin^{m}\theta\cos^{n}\theta\,d\theta$ with $m, n > -1$ , match exponents: $2p-1 = m$ and $2q-1 = n$ , so
$\int_0^{\pi/2}\sin^{m}\!\theta\,\cos^{n}\!\theta\,d\theta \;=\; \tfrac12 B\!\left(\frac{m+1}{2},\frac{n+1}{2}\right) \;=\; \frac{\Gamma\!\left(\frac{m+1}{2}\right)\Gamma\!\left(\frac{n+1}{2}\right)}{2\,\Gamma\!\left(\frac{m+n}{2}+1\right)} .$
Integrals over $[0,\pi]$ or $[0,2\pi]$ reduce to this one by symmetry, with attention to signs when an exponent is odd.
An algebraic integrand on $[0,1]$ . Bring $\int_0^1 x^{\alpha}(1-x^{k})^{\beta}dx$ to the standard shape with $t = x^{k}$ : it becomes $k^{-1}B\bigl(\frac{\alpha+1}{k},\,\beta+1\bigr)$ . The conditions for convergence are $\alpha > -1$ and $\beta > -1$ .
An algebraic integrand on $[0,\infty)$ . Use part 4 of the statement. The signature to look for is a numerator power and a denominator power of $1 + u$ (or, after scaling $u \to u/c$ , of $c + u$ ): the numerator exponent fixes $p$ and the denominator exponent fixes $p+q$ .
Then evaluate the $\Gamma$ 's. Use part 2 and the method of The Gamma function: step the arguments into $[1,2]$ , use $\Gamma(\tfrac12) = \sqrt\pi$ for half-integers, and reflection when two arguments add to an integer. If $p+q$ is an integer the answer is elementary; if $p + q$ is a half-integer it carries a $\sqrt\pi$ ; if $p$ and $q$ have denominator $3$ or $4$ it carries a $\Gamma(\tfrac13)$ or $\Gamma(\tfrac14)$ , which is not elementary.
Check the easy cases. Set one argument to $1$ : $B(p,1) = 1/p$ should reproduce the integral of a pure power. This catches factor-of- $k$ errors from step 2 and swapped parameters from step 1.

Example(A trigonometric integral)Original

Goal: evaluate $\displaystyle I = \int_0^{\pi/2}\sin^{4}\!\theta\,\cos^{6}\!\theta\;d\theta$ .

Match exponents. Step 1 of the method with $m = 4$ , $n = 6$ gives $p = \tfrac52$ , $q = \tfrac72$ , so $I = \tfrac12 B(\tfrac52,\tfrac72)$ .
Convert to $\Gamma$ .
$B\!\left(\tfrac52,\tfrac72\right) \;=\; \frac{\Gamma(\tfrac52)\,\Gamma(\tfrac72)}{\Gamma(6)} .$
Evaluate the Gammas. By the functional equation from $\Gamma(\tfrac12) = \sqrt\pi$ : $\Gamma(\tfrac52) = \tfrac34\sqrt\pi$ and $\Gamma(\tfrac72) = \tfrac{15}{8}\sqrt\pi$ , while $\Gamma(6) = 5! = 120$ . Hence
$B\!\left(\tfrac52,\tfrac72\right) \;=\; \frac{\tfrac34\cdot\tfrac{15}{8}\,\pi}{120} \;=\; \frac{45\pi}{3840} \;=\; \frac{3\pi}{256}, \qquad I \;=\; \frac{3\pi}{512} \;=\; 0.018407769 .$
Note that the two half-integer arguments supplied $\sqrt\pi$ each, so the answer is a rational multiple of $\pi$ , as step 4 of the method predicts for $p+q = 6$ an integer.
Checks. Numerical quadrature of $I$ returns $0.018407769$ . The symmetry check also holds: swapping the exponents, $\int_0^{\pi/2}\sin^6\theta\cos^4\theta\,d\theta = \tfrac12 B(\tfrac72,\tfrac52)$ , the same number, as it must be by $\theta \mapsto \pi/2-\theta$ .
The same formula for even powers of a single function. Setting $n = 0$ in step 1 gives the Wallis values, for instance $\int_0^{\pi/2}\sin^{6}\theta\,d\theta = \tfrac12 B(\tfrac72,\tfrac12) = \Gamma(\tfrac72)\Gamma(\tfrac12)/(2\Gamma(4)) = \tfrac{15}{96}\pi = 0.4908739$ .

Example(An algebraic integral on the half-line)Original

Goal: evaluate $\displaystyle J = \int_0^{\infty}\frac{\sqrt{x}}{(1+x)^{3}}\,dx$ and $\displaystyle K = \int_{-1}^{1}\left(1-x^{2}\right)^{3/2}dx$ .

$J$ by part 4. The numerator is $x^{\,p-1}$ with $p - 1 = \tfrac12$ , so $p = \tfrac32$ ; the denominator is $(1+x)^{\,p+q}$ with $p + q = 3$ , so $q = \tfrac32$ . Both have positive real part, so
$J \;=\; B\!\left(\tfrac32,\tfrac32\right) \;=\; \frac{\Gamma(\tfrac32)^{2}}{\Gamma(3)} \;=\; \frac{(\tfrac12\sqrt\pi)^{2}}{2} \;=\; \frac{\pi}{8} \;=\; 0.3926991 .$
Numerical quadrature gives $0.3926991$ .
$K$ by step 2 of the method. The integrand is even, so $K = 2\int_0^1 (1-x^2)^{3/2}dx$ . Substituting $t = x^{2}$ , that is $\int_0^1 t^{-1/2}(1-t)^{3/2}dt = B(\tfrac12,\tfrac52)$ , so
$K \;=\; B\!\left(\tfrac12,\tfrac52\right) \;=\; \frac{\Gamma(\tfrac12)\Gamma(\tfrac52)}{\Gamma(3)} \;=\; \frac{\sqrt\pi\cdot\tfrac34\sqrt\pi}{2} \;=\; \frac{3\pi}{8} \;=\; 1.1780972 .$
Check $K$ geometrically. $K$ is the area under $y = (1-x^2)^{3/2}$ , and with $x = \sin\theta$ it is $\int_{-\pi/2}^{\pi/2}\cos^{4}\theta\,d\theta = 2\cdot\tfrac12 B(\tfrac12,\tfrac52)$ by step 1 with $m = 0$ , $n = 4$ — the same number by a different route. Numerical quadrature gives $1.1780972$ .
Check against the easy case. Step 5 of the method: with the exponent $3/2$ replaced by $0$ , step 2 gives $B(\tfrac12,1) = 2$ , which is indeed $\int_{-1}^1 dx$ .

ApplicationAdaptedAdapted from Arfken, Weber & Harris (2013), Exercises 13.3.13 and 13.3.16, pp. 621-622; Grimmett & Stirzaker (2020), §4.4 (4.4.8) and Exercise 4.4.2, pp. 108-109

The period of an anharmonic oscillator. A particle of mass $m$ in the symmetric potential $V(x) = A\lvert x\rvert^{n}$ ( $A > 0$ , $n > 0$ ) with energy $E$ turns at $x_{\max} = (E/A)^{1/n}$ . From $\tfrac12 m\dot x^2 + V = E$ , a quarter period takes the particle from $0$ to $x_{\max}$ , so

\tau \;=\; 4\int_0^{x_{\max}}\frac{dx}{\sqrt{2(E - A x^{n})/m}} \;=\; 2\sqrt{2m}\int_0^{x_{\max}}\frac{dx}{\sqrt{E - Ax^{n}}} .

Substitute $x = x_{\max}u^{1/n}$ , so that $Ax^{n} = Eu$ and $dx = (x_{\max}/n)u^{1/n - 1}du$ . The integral becomes a Beta integral with $p = 1/n$ and $q = \tfrac12$ :

\tau \;=\; \frac{2\sqrt{2m}\;x_{\max}}{n\sqrt E}\int_0^1 u^{\frac1n - 1}(1-u)^{-1/2}\,du \;=\; \frac{2}{n}\sqrt{\frac{2\pi m}{E}}\left(\frac EA\right)^{1/n}\frac{\Gamma(1/n)}{\Gamma\!\left(\frac1n + \frac12\right)} ,

using $B(1/n,\tfrac12) = \sqrt\pi\,\Gamma(1/n)/\Gamma(\frac1n+\frac12)$ . This is the form given as Arfken, Weber and Harris (Arfken, Weber & Harris (2013), Exercise 13.3.13, p. 621). Two features are physics, not algebra. At $n = 2$ the Gamma quotient is $\Gamma(\tfrac12)/\Gamma(1) = \sqrt\pi$ and the energy dependence cancels, leaving $\tau = 2\pi/\omega$ with $A = m\omega^2/2$ : the harmonic case is the only one whose period is amplitude-independent. For $n > 2$ the factor $E^{1/n - 1/2}$ makes a stiffer well oscillate faster at higher energy. Numerically, for $m = 1\ \mathrm{kg}$ , $A = 1\ \mathrm{J\,m^{-4}}$ , $n = 4$ and $E = 1\ \mathrm{J}$ , the quotient is $\Gamma(\tfrac14)/\Gamma(\tfrac34) = 3.6256099/1.2254167 = 2.9586751$ and

\tau \;=\; \tfrac12\sqrt{2\pi}\,(2.9586751) \;=\; 3.7081494\ \mathrm{s},

which direct numerical integration of the quarter-period integral confirms to seven digits.

The Beta distribution. $B(a,b)$ is by construction the normaliser of the density $f(x) = x^{\,a-1}(1-x)^{\,b-1}/B(a,b)$ on $[0,1]$ , for $a, b > 0$ (Grimmett and Stirzaker (Grimmett & Stirzaker (2020), §4.4 (4.4.8), p. 108)). Its moments are ratios of Beta functions, and part 2 collapses them: $\mathbb{E}[X] = B(a+1,b)/B(a,b) = a/(a+b)$ , and the same computation with $a+2$ gives

\operatorname{Var}X \;=\; \frac{ab}{(a+b)^{2}(a+b+1)} .

At $a = 3$ , $b = 2$ this is a mean of $0.6$ and a variance of $6/150 = 0.04$ , a standard deviation of $0.2$ . The family contains the uniform distribution at $a = b = 1$ , and because its density has the same shape as the likelihood of a sequence of independent trials, it is the conjugate prior for a success probability — the reason Beta functions appear throughout Bayesian inference (Bayesian inference).

Weakly singular kernels. The kernel $(z-x)^{\alpha-1}$ of the Abel integral equation is inverted by convolving with $(y-z)^{-\alpha}$ , and the computation that makes the inversion work is exactly $\int_0^1 s^{-\alpha}(1-s)^{\alpha-1}ds = B(1-\alpha,\alpha) = \pi/\sin\pi\alpha$ , part 6 above (Volterra equations). The same Beta integral is what gives the fractional integration operators their composition rule.

Definition

Digamma and polygamma functions

Adapted

Add to PDF

Defines $\psi = \Gamma'/\Gamma$ and its derivatives, derives the series $\psi(z) = -\gamma - 1/z + \sum_{k\ge1}(1/k - 1/(z+k))$ and $\psi^{(m)}(z) = (-1)^{m+1}m!\sum_{k\ge0}(z+k)^{-m-1}$ from the Weierstrass product, collects the values at integers and half-integers, and turns a convergent series of rational terms into a finite combination of digamma values.

Prerequisites: The Gamma function, Analytic continuation, Integrating and differentiating series term by term, Uniform convergence

Definition(Digamma and polygamma)AdaptedAdapted from Arfken, Weber & Harris (2013), §13.2, eqs. (13.36)-(13.37) and (13.41), pp. 611-612; Hassani (2013), eq. (12.15), p. 380

On $\Omega = \mathbb{C}\setminus\{0,-1,-2,\dots\}$ , where $\Gamma$ is analytic and non-vanishing (The Gamma function), the digamma function is the logarithmic derivative of $\Gamma$ ,

\psi(z) \;=\; \frac{\Gamma'(z)}{\Gamma(z)} \;=\; \frac{d}{dz}\ln\Gamma(z) ,

and for $m \in \mathbb{N}$ the polygamma functions are its derivatives,

\psi^{(m)}(z) \;=\; \frac{d^{m}\psi}{dz^{m}}(z) \;=\; \frac{d^{\,m+1}}{dz^{\,m+1}}\ln\Gamma(z) .

All of them are analytic on $\Omega$ . The expression $\ln\Gamma$ is many-valued in general, but $\Gamma'/\Gamma$ is not, and it is the quotient that is meant: on the interval $(0,\infty)$ , where $\Gamma > 0$ , the two readings agree with $\ln\Gamma$ taken real.

Also needed is the Euler-Mascheroni constant

\gamma \;=\; \lim_{n\to\infty}\bigl(H_n - \ln n\bigr) \;=\; 0.5772156649\dots, \qquad H_n \;=\; \sum_{k=1}^{n}\frac1k ,

the limit existing because $H_n - \ln n$ is decreasing and bounded below.

A warning about the argument. Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §13.2, eqs. (13.36)-(13.41), pp. 611-612) write every formula of this section for $\psi(z+1)$ and $\psi^{(m)}(z+1)$ , because $\Gamma(z+1)$ is the quantity they treat as the factorial. Their $\psi(z+1)$ is our $\psi$ evaluated at $z+1$ , so their formulas are ours shifted; part 2 of the statement is the dictionary. Mixing the two conventions is the commonest source of an answer wrong by $1/z$ .

Definition(Series, values and recurrences)AdaptedAdapted from Arfken, Weber & Harris (2013), §13.2, eqs. (13.38)-(13.44) and Exercises 13.2.4, 13.2.8, 13.2.10, pp. 611-616; Hassani (2013), eq. (12.15), p. 380

Series for $\psi$ . For every $z \in \Omega$ ,
$\psi(z) \;=\; -\gamma - \frac1z + \sum_{k=1}^{\infty}\left(\frac1k - \frac{1}{z+k}\right) \;=\; -\gamma - \frac1z + \sum_{k=1}^{\infty}\frac{z}{k(z+k)} ,$
the series converging absolutely, and uniformly on each compact subset of $\Omega$ . Equivalently $\psi(1+z) = -\gamma + \sum_{k\ge1} z/\bigl(k(k+z)\bigr)$ , which is the form Arfken, Weber and Harris (Arfken, Weber & Harris (2013), eq. (13.38), p. 611) give.
Recurrence and integer values. $\psi(z+1) = \psi(z) + 1/z$ on $\Omega$ , and
$\psi(1) \;=\; -\gamma , \qquad \psi(n+1) \;=\; -\gamma + H_n \quad (n = 1, 2, \dots).$
Half-integer values, reflection, duplication.
$\psi\!\left(\tfrac12\right) = -\gamma - 2\ln 2 , \qquad \psi(1-z) - \psi(z) = \pi\cot\pi z , \qquad \psi(z) + \psi\!\left(z+\tfrac12\right) = 2\psi(2z) - 2\ln 2 .$
Consequently $\psi(n+\tfrac12) = -\gamma - 2\ln2 + 2\sum_{j=1}^{n}(2j-1)^{-1}$ .
Polygamma series. For $m \ge 1$ and $z \in \Omega$ ,
$\psi^{(m)}(z) \;=\; (-1)^{m+1}\,m!\sum_{k=0}^{\infty}\frac{1}{(z+k)^{\,m+1}} ,$
with recurrence $\psi^{(m)}(z+1) = \psi^{(m)}(z) + (-1)^{m}m!\,z^{-m-1}$ , and at $z = 1$
$\psi^{(m)}(1) \;=\; (-1)^{m+1}\,m!\;\zeta(m+1) ,$
where $\zeta$ is the Riemann zeta function (Riemann zeta function and Bernoulli numbers). In particular $\psi^{(1)}(1) = \zeta(2) = \pi^2/6$ and $\psi^{(2)}(1) = -2\zeta(3)$ .
Maclaurin series of $\ln\Gamma$ . For $\lvert z\rvert < 1$ ,
$\ln\Gamma(1+z) \;=\; -\gamma z + \sum_{n=2}^{\infty}(-1)^{n}\,\frac{\zeta(n)}{n}\,z^{n} , \qquad \psi(1+z) \;=\; -\gamma + \sum_{n=2}^{\infty}(-1)^{n}\zeta(n)\,z^{\,n-1} .$
On the real axis the first series converges for $-1 < z \le 1$ .
Monotonicity on the positive axis. For $x > 0$ , $\psi^{(1)}(x) > 1/x > 0$ and $\psi^{(2)}(x) < 0$ ; so $\psi$ is strictly increasing and concave on $(0,\infty)$ , $\ln\Gamma$ is strictly convex there, and $x \mapsto \ln x - \psi(x)$ is strictly decreasing. Moreover $\psi(x) \to -\infty$ as $x \downarrow 0$ and $\psi(x) \to +\infty$ as $x \to \infty$ , so $\psi$ maps $(0,\infty)$ bijectively onto $\mathbb{R}$ .

Idea of the proofOriginal

A product becomes a sum under a logarithm, and a sum can be differentiated term by term. That single observation is the whole theory: $\Gamma$ itself is awkward to differentiate, but $1/\Gamma$ is an infinite product of linear factors, so $\ln\Gamma$ is a series of logarithms of linear factors and $\psi$ is a series of reciprocals of linear factors. Each further derivative raises the power and makes the series converge faster, which is why the polygamma series are the clean ones and the digamma series is the one that needs the Euler-Mascheroni constant to repair a divergence: $\sum_k 1/(z+k)$ diverges on its own, and $\gamma$ is precisely the bookkeeping left over when one subtracts the harmonic series from it.

Everything else is read off that series. Differences of $\psi$ at arguments shifted by integers telescope, giving the recurrence and the harmonic numbers. The functional equations of $\Gamma$ — reflection and duplication — become additive identities for $\psi$ on taking logarithms and differentiating, and those supply the half-integer values. Setting $z = 1$ in the polygamma series turns it into $\zeta(m+1)$ , which is what makes $\psi$ and its derivatives the natural language for summing series of rational terms: partial fractions expresses such a series as a combination of the very sums the polygamma functions tabulate.

ProofAdaptedAdapted from Arfken, Weber & Harris (2013), §13.2, eqs. (13.34)-(13.44), pp. 610-613, with the Weierstrass product, eq. (13.15), p. 602; Hassani (2013), eqs. (12.14)-(12.15), pp. 379-380

(1) The series. Step 1: the Weierstrass product. Gauss's limit formula (The Gamma function, part 4) says $\Gamma(z) = \lim_n n!\,n^{z}/\bigl(z(z+1)\cdots(z+n)\bigr)$ . Since $(z+1)\cdots(z+n)/n! = \prod_{k=1}^n(1+z/k)$ , taking reciprocals gives

\frac{1}{\Gamma(z)} \;=\; \lim_{n\to\infty} z\,n^{-z}\prod_{k=1}^{n}\left(1+\frac zk\right) .

Write $n^{-z} = e^{-z\ln n} = e^{\gamma z}\,e^{-zH_n}\,e^{\,z(H_n - \ln n - \gamma)}$ . The last factor tends to $1$ because $H_n - \ln n \to \gamma$ , and $e^{-zH_n} = \prod_{k=1}^n e^{-z/k}$ , so

\frac{1}{\Gamma(z)} \;=\; z\,e^{\gamma z}\prod_{k=1}^{\infty}\left(1+\frac zk\right)e^{-z/k} ,

the product of Arfken, Weber and Harris (Arfken, Weber & Harris (2013), eq. (13.15), p. 602) and Hassani (Hassani (2013), eq. (12.14), p. 379). Its factors tend to $1$ fast: $\ln\bigl[(1+u)e^{-u}\bigr] = -u^2/2 + O(u^3)$ , so with $u = z/k$ the logarithms are $O(\lvert z\rvert^2/k^2)$ .

Step 2: logarithms and differentiation, for real argument. Fix $0 < a < b$ and let $x \in [a,b]$ . Every factor above is positive, so logarithms may be taken:

-\ln\Gamma(x) \;=\; \ln x + \gamma x + \sum_{k=1}^{\infty}\left[\ln\left(1+\frac xk\right) - \frac xk\right] .

The series converges by the $O(x^2/k^2)$ bound of step 1. Differentiating the $k$ -th term gives $1/(k+x) - 1/k = -x/\bigl(k(k+x)\bigr)$ , of modulus at most $b/k^2$ on $[a,b]$ ; so the differentiated series converges uniformly there by the Weierstrass M-test (Uniform convergence), and term-by-term differentiation is legitimate (Integrating and differentiating series term by term, part 2, whose hypotheses — continuously differentiable terms, convergence at a point, uniform convergence of the differentiated series — are exactly what we have checked). Hence

-\psi(x) \;=\; \frac1x + \gamma + \sum_{k=1}^{\infty}\left(\frac{1}{k+x} - \frac1k\right),

which rearranges to part 1 for $x > 0$ .

Step 3: extension to $\Omega$ . Let $S_N(z) = -\gamma - 1/z + \sum_{k=1}^{N}\bigl(1/k - 1/(z+k)\bigr)$ , analytic on $\Omega$ . On a compact $K \subset \Omega$ take $R$ with $\lvert z\rvert \le R$ on $K$ ; for $k > 2R$ we have $\lvert z+k\rvert \ge k/2$ and so

\left\lvert \frac1k - \frac{1}{z+k}\right\rvert \;=\; \frac{\lvert z\rvert}{k\,\lvert z+k\rvert} \;\le\; \frac{2R}{k^{2}} ,

whence $S_N$ converges uniformly on $K$ by the M-test. A uniform limit of continuous functions is continuous, and for a closed contour $\gamma \subset \Omega$ uniform convergence lets the limit pass through the contour integral, so $\oint_\gamma \lim_N S_N = \lim_N \oint_\gamma S_N = 0$ by Cauchy's theorem; by Morera's theorem the limit $S$ is analytic on $\Omega$ (both steps exactly as in Analytic continuation). Now $S$ and $\psi$ are analytic on $\Omega$ and agree on $(0,\infty)$ by step 2, a set with accumulation points in $\Omega$ , so they agree on all of $\Omega$ by the identity theorem. The second form of part 1 is the algebraic identity $1/k - 1/(z+k) = z/\bigl(k(z+k)\bigr)$ , and adding $1/z$ to both sides of part 1 gives the form in $\psi(1+z)$ once part 2 is known.

(2) Recurrence and integer values. Differentiating the logarithm of $\Gamma(z+1) = z\Gamma(z)$ gives $\psi(z+1) = 1/z + \psi(z)$ on $(0,\infty)$ and hence on $\Omega$ by analyticity. At $z = 1$ , part 1 reads

\psi(1) \;=\; -\gamma - 1 + \sum_{k=1}^{\infty}\left(\frac1k - \frac{1}{1+k}\right) \;=\; -\gamma - 1 + 1 \;=\; -\gamma ,

the series telescoping to $1$ . Iterating the recurrence from $\psi(1)$ gives $\psi(n+1) = -\gamma + 1 + \tfrac12 + \dots + \tfrac1n$ .

(3) Reflection, duplication, half-integers. Take logarithmic derivatives of the two functional equations of The Gamma function. From $\Gamma(z)\Gamma(1-z) = \pi/\sin\pi z$ , differentiating the logarithm in $z$ and noting that the chain rule contributes $-1$ from the second factor,

\psi(z) - \psi(1-z) \;=\; -\pi\cot\pi z ,

which is the stated reflection relation. From $\Gamma(z)\Gamma(z+\tfrac12) = 2^{\,1-2z}\sqrt\pi\;\Gamma(2z)$ , taking logarithms gives $\ln\Gamma(z) + \ln\Gamma(z+\tfrac12) = (1-2z)\ln2 + \tfrac12\ln\pi + \ln\Gamma(2z)$ , and differentiating,

\psi(z) + \psi\!\left(z+\tfrac12\right) \;=\; -2\ln2 + 2\,\psi(2z) .

Put $z = \tfrac12$ : the left side is $\psi(\tfrac12) + \psi(1)$ and the right side is $-2\ln2 + 2\psi(1)$ , so $\psi(\tfrac12) = \psi(1) - 2\ln2 = -\gamma - 2\ln 2 = -1.9635100$ . The values at $n+\tfrac12$ follow from the recurrence of part 2, since $\sum_{j=1}^{n}(j-\tfrac12)^{-1} = 2\sum_{j=1}^{n}(2j-1)^{-1}$ .

(4) Polygamma. Differentiating part 1 once gives $\psi^{(1)}(z) = 1/z^2 + \sum_{k\ge1}(z+k)^{-2} = \sum_{k\ge0}(z+k)^{-2}$ ; the differentiated series converges uniformly on compact subsets of $\Omega$ by the M-test (the terms are $O(k^{-2})$ ), so the differentiation is termwise by the argument of step 3. Differentiating $m-1$ further times, each time legitimately because the series only converges faster, gives

\psi^{(m)}(z) \;=\; \sum_{k=0}^{\infty}\frac{d^{\,m-1}}{dz^{\,m-1}}\frac{1}{(z+k)^{2}} \;=\; (-1)^{m-1}m!\sum_{k=0}^{\infty}\frac{1}{(z+k)^{\,m+1}} ,

and $(-1)^{m-1} = (-1)^{m+1}$ . The recurrence is part 2 differentiated $m$ times, using $d^{m}(1/z)/dz^{m} = (-1)^{m}m!\,z^{-m-1}$ . At $z = 1$ the sum is $\sum_{k\ge0}(1+k)^{-m-1} = \zeta(m+1)$ , convergent because $m + 1 \ge 2$ .

(5) Maclaurin series. $\ln\Gamma(1+z)$ is analytic on the unit disc — $\Gamma(1+z)$ is analytic and non-vanishing for $\lvert z\rvert < 1$ , and $\ln\Gamma(1+0) = 0$ fixes the branch — so it equals its Taylor series there, with radius exactly $1$ because $\Gamma(1+z)$ has a pole at $z = -1$ . The coefficients are the derivatives at $0$ : the first is $\psi(1) = -\gamma$ by part 2, and for $n \ge 2$ the $n$ -th derivative is $\psi^{(n-1)}(1) = (-1)^{n}(n-1)!\,\zeta(n)$ by part 4, so the $n$ -th coefficient is $(-1)^{n}(n-1)!\zeta(n)/n! = (-1)^{n}\zeta(n)/n$ . Differentiating gives the series for $\psi(1+z)$ .

(6) Monotonicity. Part 4 with $m = 1$ gives $\psi^{(1)}(x) = \sum_{k\ge0}(x+k)^{-2}$ for $x > 0$ , and with $m = 2$ gives $\psi^{(2)}(x) = -2\sum_{k\ge0}(x+k)^{-3} < 0$ . The first sum is a left-endpoint Riemann sum, over the unit intervals $[k,k+1]$ , of the strictly decreasing function $t \mapsto (x+t)^{-2}$ , so it strictly exceeds the integral:

\psi^{(1)}(x) \;=\; \sum_{k=0}^{\infty}\frac{1}{(x+k)^{2}} \;>\; \int_0^{\infty}\frac{dt}{(x+t)^{2}} \;=\; \frac1x \;>\; 0 .

Hence $\psi$ is strictly increasing and concave; $(\ln\Gamma)'' = \psi^{(1)} > 0$ makes $\ln\Gamma$ strictly convex; and the derivative of $\ln x - \psi(x)$ is $1/x - \psi^{(1)}(x) < 0$ , so that function is strictly decreasing. For the limits: the single term $-1/x$ in part 1 drives $\psi(x) \to -\infty$ as $x \downarrow 0$ , the remaining terms staying bounded; and $\psi(n+1) = -\gamma + H_n \to \infty$ with monotonicity in between. A continuous strictly increasing function with those limits is a bijection onto $\mathbb{R}$ . $\square$

Remark(Where the hypotheses are used)Original

The constant $\gamma$ is not decoration. In part 1 neither $\sum_k 1/k$ nor $\sum_k 1/(z+k)$ converges; only the difference does, and $\gamma$ is the price of having split a convergent series into two divergent halves in the Weierstrass product. This is why the digamma series, alone among the family, contains a transcendental constant: for $m \ge 1$ the polygamma series in part 4 converges on its own and no constant appears.
Term-by-term differentiation needs the differentiated series to converge uniformly, not merely the original one; that is the hypothesis of Integrating and differentiating series term by term and it is what the $b/k^2$ bound in step 2 supplies. The bound degrades as $b \to \infty$ , which is why the argument is run on a compact interval and then extended, rather than on $(0,\infty)$ at once.
Step 2 is carried out on the positive axis for a reason. Taking the logarithm of the product requires choosing a branch, and for $x > 0$ every factor is a positive real so the choice is forced. The identity theorem then carries the result to all of $\Omega$ without any branch question, because $\psi = \Gamma'/\Gamma$ is single-valued even though $\ln\Gamma$ is not.
The reflection relation of part 3 holds only off the integers, where $\cot \pi z$ is finite and $\Gamma(z)\Gamma(1-z)$ makes sense; the poles of $\pi\cot\pi z$ at the integers match the poles of $\psi$ there. Each pole of $\psi$ is simple with residue $-1$ , as the term $-1/(z+k)$ of part 1 shows.
Part 4 excludes $m = 0$ . The polygamma series as written would be $\sum_k 1/(z+k)$ at $m = 0$ , which diverges; part 1 is the correct $m = 0$ statement, with the compensating $1/k$ .
The convergence condition in the method below is not optional. A partial-fraction decomposition $\sum_j c_j/(n+a_j)$ sums to a combination of digamma values only when $\sum_j c_j = 0$ , which is the same as requiring the numerator degree to be at most the denominator degree minus two. Otherwise the series diverges and no rearrangement can rescue it.

Method(Summing a series of rational terms)AdaptedAdapted from Arfken, Weber & Harris (2013), §13.2, Series Summation, p. 613, and Exercises 13.2.14-13.2.15, p. 616

To evaluate $\sum_{n\ge1} R(n)$ where $R = P/Q$ is a rational function with no pole at a positive integer:

Check convergence. Require $\deg P \le \deg Q - 2$ . If $\deg P = \deg Q - 1$ the series diverges logarithmically; split off the offending part first, or pair terms so that the leading behaviour cancels.
Factor and decompose. Write $Q(n) = \prod_j (n + a_j)^{\,m_j+1}$ and expand $R$ in partial fractions. With simple factors only, $R(n) = \sum_j c_j/(n+a_j)$ , and step 1 guarantees $\sum_j c_j = 0$ .
Replace each simple group by digamma values. From part 1, $\psi(1+a) = -\gamma + \sum_{k\ge1}\bigl(1/k - 1/(k+a)\bigr)$ . Summing over $j$ with weights $c_j$ , the $\gamma$ terms and the $1/k$ terms both cancel because $\sum_j c_j = 0$ , leaving
$\sum_{n=1}^{\infty}\sum_j \frac{c_j}{n+a_j} \;=\; -\sum_j c_j\,\psi(1+a_j) .$
The two-factor case is worth memorising: for $a \ne b$ ,
$\sum_{n=1}^{\infty}\frac{1}{(n+a)(n+b)} \;=\; \frac{\psi(1+b) - \psi(1+a)}{b-a} .$
Use polygamma for repeated factors. A term $(n+a)^{-m-1}$ with $m \ge 1$ needs no cancellation: by part 4,
$\sum_{n=0}^{\infty}\frac{1}{(n+a)^{\,m+1}} \;=\; \frac{(-1)^{m+1}}{m!}\,\psi^{(m)}(a) ,$
and a sum starting at $n = 1$ is this minus $a^{-m-1}$ .
Evaluate the digamma values. Use part 2 for arguments differing from $1$ by an integer (the answer is then rational, built from harmonic numbers), part 3 for half-integers (expect $\ln 2$ ), and the reflection relation for arguments summing to an integer (expect $\pi\cot$ ). Arguments with denominator $3$ , $4$ or $6$ give closed forms in $\pi$ and logarithms; other arguments must be left as $\psi$ or evaluated numerically.
Check against a partial sum. Add the first few hundred terms numerically. Convergence is only like $1/N$ for a two-factor denominator, so compare with the tail estimate rather than expecting agreement to full precision.

Example(A telescoping sum, and what the method automates)Original

Goal: evaluate $\displaystyle \Sigma_1 = \sum_{n=1}^{\infty}\frac{1}{n(n+3)}$ .

Convergence. The numerator has degree $0$ and the denominator degree $2$ , so step 1 is satisfied.
Partial fractions. $\dfrac{1}{n(n+3)} = \dfrac13\left(\dfrac1n - \dfrac{1}{n+3}\right)$ , so $c_1 = \tfrac13$ at $a_1 = 0$ and $c_2 = -\tfrac13$ at $a_2 = 3$ , and indeed $c_1 + c_2 = 0$ .
Apply step 3. With $a = 0$ , $b = 3$ ,
$\Sigma_1 \;=\; \frac{\psi(4) - \psi(1)}{3} .$
Evaluate. By part 2, $\psi(4) = -\gamma + H_3 = -\gamma + \tfrac{11}{6}$ and $\psi(1) = -\gamma$ , so the $\gamma$ cancels and
$\Sigma_1 \;=\; \frac{11}{18} \;=\; 0.6111111 .$
Checks. Direct summation of $N = 10^{6}$ terms gives $0.6111101$ , short of the answer by $1.0\times10^{-6}$ . That is exactly the omitted tail: $\tfrac13\sum_{n>N}\bigl(1/n - 1/(n+3)\bigr) = \tfrac13\bigl(1/(N{+}1)+1/(N{+}2)+1/(N{+}3)\bigr) \approx 1/N = 10^{-6}$ , the slow convergence step 6 warns about. The telescoping is visible here, so the answer can be had by hand; the point of the method is that it survives when the roots are not integers. The same three lines with $a = 0$ , $b = \tfrac12$ give
$\sum_{n=1}^{\infty}\frac{1}{n\bigl(n+\tfrac12\bigr)} \;=\; 2\bigl(\psi(\tfrac32)-\psi(1)\bigr) \;=\; 4 - 4\ln 2 \;=\; 1.2274113 ,$
where $\psi(\tfrac32) = \psi(\tfrac12) + 2 = 2-\gamma-2\ln2 = 0.0364900$ by parts 2 and 3; summing $10^{6}$ terms gives $1.2274103$ .

Example(An odd-denominator sum by polygamma)Original

Goal: evaluate $\displaystyle \Sigma_2 = \sum_{n=1}^{\infty}\frac{1}{(2n-1)^{2}} = 1 + \tfrac19 + \tfrac1{25} + \cdots$ without assuming the value of $\zeta(2)$ at the odd integers separately.

Shift to the polygamma shape. $(2n-1)^2 = 4\bigl(n - \tfrac12\bigr)^2$ , so with $k = n-1$
$\Sigma_2 \;=\; \frac14\sum_{k=0}^{\infty}\frac{1}{\bigl(k+\tfrac12\bigr)^{2}} \;=\; \frac14\,\psi^{(1)}\!\left(\tfrac12\right)$
by step 4 of the method with $m = 1$ , $a = \tfrac12$ .
Get $\psi^{(1)}(\tfrac12)$ from duplication. Differentiate the duplication relation of part 3 once in $z$ :
$\psi^{(1)}(z) + \psi^{(1)}\!\left(z+\tfrac12\right) \;=\; 4\,\psi^{(1)}(2z) .$
At $z = \tfrac12$ this reads $\psi^{(1)}(\tfrac12) + \psi^{(1)}(1) = 4\psi^{(1)}(1)$ , so $\psi^{(1)}(\tfrac12) = 3\psi^{(1)}(1) = 3\zeta(2)$ by part 4.
Conclude.
$\Sigma_2 \;=\; \tfrac34\,\zeta(2) \;=\; \frac{3}{4}\cdot\frac{\pi^{2}}{6} \;=\; \frac{\pi^{2}}{8} \;=\; 1.2337006 .$
Note that $\psi^{(1)}(\tfrac12) = \pi^2/2 = 4.9348022$ .
Checks. Summing $N = 10^{6}$ terms of $\Sigma_2$ gives $1.2337003$ , the tail being about $1/(4N) = 2.5\times10^{-7}$ . Independently, $\Sigma_2$ should be $\zeta(2)$ minus the even part $\sum_{n\ge1}(2n)^{-2} = \tfrac14\zeta(2)$ , which is $\tfrac34\zeta(2)$ — the same answer, and a reminder that step 2 did the work of that elementary split without needing it.

ApplicationAdaptedAdapted from Arfken, Weber & Harris (2013), Exercises 13.1.21 and 13.2.12-13.2.13, pp. 610 and 616, and §13.2, Series Summation, p. 613

The log-mean of a Gamma variable, and fitting a shape parameter. The Gamma density with rate $\lambda$ and shape $t$ is $f(x) = \lambda^{t}x^{\,t-1}e^{-\lambda x}/\Gamma(t)$ on $x > 0$ (Arfken, Weber and Harris (Arfken, Weber & Harris (2013), Exercise 13.1.21, p. 610)). Its fractional moments are

\mathbb{E}[X^{s}] \;=\; \frac{\lambda^{t}}{\Gamma(t)}\int_0^{\infty}x^{\,t+s-1}e^{-\lambda x}\,dx \;=\; \frac{\Gamma(t+s)}{\Gamma(t)\,\lambda^{s}} ,

and differentiating in $s$ at $s = 0$ — the derivative of $\Gamma(t+s)$ at $s=0$ being $\Gamma(t)\psi(t)$ — gives

\mathbb{E}[\ln X] \;=\; \psi(t) - \ln\lambda .

For $t = 3$ , $\lambda = 2$ this is $\psi(3) - \ln 2 = 0.9227843 - 0.6931472 = 0.2296372$ , so the geometric mean is $e^{0.2296372} = 1.2581434$ , well below the arithmetic mean $t/\lambda = 1.5$ ; the gap is a measure of the skewness that makes the Gamma family useful for rainfall totals and for waiting times. The same formula is the whole content of maximum-likelihood fitting of a Gamma shape: the likelihood equations reduce to $\ln\hat t - \psi(\hat t) = \ln\bar x - \overline{\ln x}$ , a scalar equation with at most one root, because $\ln x - \psi(x)$ is strictly decreasing on $(0,\infty)$ by part 6 (Maximum likelihood).

Phases and small expansions in quantum scattering. A charged particle in a Coulomb field has a wave function whose normalisation and phase involve $\Gamma(1 + i\eta)$ with $\eta$ the Sommerfeld parameter (Arfken, Weber and Harris (Arfken, Weber & Harris (2013), Exercise 13.2.13, p. 616)). Part 5 gives the small- $\eta$ behaviour directly: $\ln\Gamma(1+i\eta) = -\gamma\,i\eta - \zeta(2)\eta^{2}/2 + O(\eta^{3})$ , so $\arg\Gamma(1+i\eta) = -\gamma\eta + O(\eta^{3})$ and $\lvert\Gamma(1+i\eta)\rvert^{2} = 1 - \zeta(2)\eta^{2} + O(\eta^{4})$ . The same series, with $z$ a small algebraic quantity rather than imaginary, is how relativistic hydrogen wave functions are expanded in powers of $\alpha^{2}Z^{2}$ (Arfken, Weber and Harris (Arfken, Weber & Harris (2013), Exercise 13.2.12, p. 616)).

Lattice and loop sums. Any sum whose general term is a rational function of the index — a one-dimensional lattice Green's function, a dimensionally regularised loop integral expanded about an integer dimension, a mode sum in a cavity — reduces by the method above to a finite combination of $\psi$ and $\psi^{(m)}$ at rational arguments, and hence to $\gamma$ , $\ln 2$ , $\pi$ and the zeta values. That reduction is why $\gamma$ and $\zeta(2), \zeta(3)$ appear in the coefficients of perturbative expansions, and why tables of $\psi^{(m)}$ were worth computing: the alternating sum $\sum_{k\ge0}(-1)^{k}(2k+1)^{-2}$ , Catalan's constant, becomes $\bigl[\psi^{(1)}(\tfrac14) - \psi^{(1)}(\tfrac34)\bigr]/16 = 0.9159656$ , a value with no known closed form in elementary constants.

Definition

Riemann zeta function and Bernoulli numbers

Adapted

Add to PDF

Defines $\zeta(s)=\sum_{n\ge1}n^{-s}$ for real $s>1$ and the Bernoulli numbers $B_n$ by the generating function $t/(e^t-1)$ , and proves Euler's identity $\zeta(2k)=(-1)^{k+1}(2\pi)^{2k}B_{2k}/\bigl(2\,(2k)!\bigr)$ , which turns a rational recursion into the closed values $\pi^2/6$ , $\pi^4/90$ , $\pi^6/945$ and gives every $\zeta$ at an even integer.

Prerequisites: Euler-Maclaurin summation, Convergence tests for series, Uniform convergence, Power series, Integrating and differentiating series term by term, Laurent series, The residue theorem, Taylor series of analytic functions, Analytic continuation, Pole expansions and infinite products, The Gamma function

Definition(The Riemann zeta function on the real half line)AdaptedAdapted from Arfken, Weber & Harris (2013), §12.4, Eq. (12.60), p. 571, and §12.2, p. 564; Stein & Shakarchi (2003), Ch. 5, §3.1, p. 155; Ahlfors (1979), Ch. 5, §4, Eq. (55), p. 212

For a real number $s > 1$ the Riemann zeta function is the sum of the convergent series

\zeta(s) \;=\; \sum_{n=1}^{\infty}\frac{1}{n^{s}} \;=\; 1 + \frac{1}{2^{s}} + \frac{1}{3^{s}} + \cdots .

Convergence for $s>1$ and divergence for $s\le1$ are decided by the integral test (Convergence tests for series). The test needs a positive decreasing function, so take $s>0$ first; then $x\mapsto x^{-s}$ is positive and decreasing on $[1,\infty)$ , and for $s \neq 1$

\int_{1}^{N}x^{-s}\,dx \;=\; \frac{N^{1-s}-1}{1-s} ,

which has a finite limit as $N\to\infty$ exactly when $s>1$ , while for $s=1$ the integral is $\ln N \to \infty$ . Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §4.3.2, pp. 128-130) record both halves of this, and also give a second proof of the divergence at $s=1$ by grouping the terms $2^{-k}+\cdots+(2^{k+1}-1)^{-1} \ge \tfrac12$ . For $s\le0$ the function $x^{-s}$ is not decreasing and the test says nothing, but there it is not needed: the terms $n^{-s}\ge1$ do not tend to $0$ . In fact a single comparison covers the whole of $s<1$ at once, since $n^{-s} \ge n^{-1}$ there.

Two matters of notation. The argument is written $s$ by Ahlfors (Ahlfors (1979), Ch. 5, §4, p. 212) and Stein and Shakarchi (Stein & Shakarchi (2003), Ch. 5, §3.1, p. 155). Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §12.4, Eq. (12.60), p. 571) write $\zeta(s)$ where they introduce the Dirichlet series, but $\zeta(z)$ where they treat the function as a function of a complex variable — at §12.2, p. 564, and throughout §13.5, pp. 629-631. We use $s$ everywhere. And the same series defines an analytic function of a complex argument on the half plane $\operatorname{Re}s > 1$ , because it converges uniformly on each closed half plane $\operatorname{Re}s\ge\sigma_0>1$ , where it is dominated by $\sum n^{-\sigma_0}$ (Ahlfors, p. 212, with the Weierstrass $M$ -test of Uniform convergence). This unit needs only real $s>1$ ; see the note below for what lies beyond.

Definition(Bernoulli numbers)AdaptedAdapted from Arfken, Weber & Harris (2013), §12.2, Eqs. (12.28)-(12.30) and Table 12.2, pp. 560-562; Stein & Shakarchi (2003), Ch. 3, Problem 4(a), p. 97

The Bernoulli numbers $B_0, B_1, B_2,\dots$ are the Taylor coefficients, scaled by factorials, of the function

G(t) \;=\; \frac{t}{e^{t}-1} \;=\; \sum_{n=0}^{\infty}\frac{B_{n}}{n!}\,t^{n} .

The quotient is defined and analytic on the punctured disc $0<\lvert t\rvert<2\pi$ and extends analytically across $t=0$ : since $e^t-1 = t\bigl(1+\tfrac{t}{2}+\tfrac{t^2}{6}+\cdots\bigr)$ has a simple zero at the origin with the bracket equal to $1$ there, $G$ has a removable singularity at $t=0$ with $G(0)=1$ . Its nearest singularities are the simple poles of $1/(e^t-1)$ at $t=\pm2\pi i$ , so the radius of convergence of the series above is exactly $2\pi$ (Taylor series of analytic functions).

Reading off the first two coefficients gives $B_0 = 1$ and, by one differentiation at the origin,

B_{1} \;=\; G'(0) \;=\; -\tfrac12 .

Multiplying the series for $G$ by $(e^t-1)/t = \sum_{m\ge0}t^{m}/(m+1)!$ and requiring the product to be $1$ gives, on comparing coefficients of $t^{N}$ , a rational recursion determining every $B_n$ from its predecessors. In the form Arfken, Weber and Harris (Arfken, Weber & Harris (2013), Eq. (12.34), p. 561) give for the even-index numbers,

N-1 \;=\; \sum_{n=1}^{N-1}\binom{2N}{2n}B_{2n} , \qquad N = 2,3,4,\dots

Stein and Shakarchi (Stein & Shakarchi (2003), Ch. 3, Problem 4(b), p. 97) state the equivalent recursion over all indices, $B_n = -\tfrac{1}{n+1}\sum_{k=0}^{n-1}\binom{n+1}{k}B_k$ for $n\ge1$ . Either gives

B_{0}=1,\quad B_{1}=-\tfrac12,\quad B_{2}=\tfrac16,\quad B_{4}=-\tfrac1{30},\quad B_{6}=\tfrac1{42},\quad B_{8}=-\tfrac1{30},\quad B_{10}=\tfrac{5}{66},\quad B_{12}=-\tfrac{691}{2730},

with $B_{n}=0$ for every odd $n\ge3$ (proved below). Arfken, Weber and Harris tabulate these as far as $B_{10}$ (Table 12.2, p. 562); the value of $B_{12}$ quoted here is the next one either recursion returns, and is not taken from the table. The Bernoulli polynomials $B_n(x)$ , generated by $te^{xt}/(e^t-1)$ , and the summation formula in which they are used belong to Euler-Maclaurin summation.

Definition(Euler's formula for

\zeta

at the even integers)AdaptedAdapted from Arfken, Weber & Harris (2013), §12.2, Eqs. (12.37)-(12.38), pp. 563-564; Stein & Shakarchi (2003), Ch. 3, Problem 4(c),(e), pp. 97-98

For every integer $n\ge2$ ,

\frac{B_{n}}{n!} \;=\; -\sum_{\substack{m\in\mathbb{Z}\\ m\neq0}}\frac{1}{(2\pi m i)^{n}} .

Consequently:

Odd indices vanish. $B_{n} = 0$ for every odd $n \ge 3$ .
Even indices are zeta values. For every integer $k\ge1$ ,
$B_{2k} \;=\; (-1)^{k+1}\,\frac{2\,(2k)!}{(2\pi)^{2k}}\,\zeta(2k) ,$
which is the form in Arfken, Weber and Harris (Arfken, Weber & Harris (2013), Eq. (12.38), p. 564); equivalently, solving for the zeta value,
$\zeta(2k) \;=\; (-1)^{k+1}\,\frac{(2\pi)^{2k}}{2\,(2k)!}\,B_{2k} .$
In particular $\zeta(2)=\pi^{2}/6$ , $\zeta(4)=\pi^{4}/90$ and $\zeta(6)=\pi^{6}/945$ .
Signs and growth. Since $\zeta(2k)>0$ , the sign of $B_{2k}$ is $(-1)^{k+1}$ : the even-index Bernoulli numbers alternate in sign. And since $\zeta(2k)\to1$ , the formula forces $\lvert B_{2k}\rvert \sim 2\,(2k)!/(2\pi)^{2k}$ , which grows without bound — Arfken, Weber and Harris (Arfken, Weber & Harris (2013), p. 564) quote $B_{20} = -5.291\times10^{2}$ and $B_{200} = -3.647\times10^{215}$ .

The identity says nothing about $\zeta$ at odd arguments; no closed form for $\zeta(3), \zeta(5),\dots$ is known, and Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §12.4, p. 572) go no further than saying that these values seem to resist evaluation in closed form.

Idea of the proofOriginal

Everything follows from one observation: the function being expanded, $t/(e^t-1)$ , has poles, and the Bernoulli numbers are its Taylor coefficients, so the coefficients must remember where the poles are. The poles of $1/(e^t-1)$ sit at the points $2\pi m i$ , evenly spaced up the imaginary axis, all simple, all with residue $1$ . A sum over an arithmetic progression of poles is a sum of the form $\sum m^{-n}$ — a zeta value.

The mechanism that converts "Taylor coefficient" into "sum over all the poles" is the residue theorem. Write $B_n/n!$ as the residue at the origin of $t^{-n}/(e^t-1)$ , by Cauchy's coefficient formula. Then integrate that same function around a large circle. The circle encloses the origin and the $2N$ nearest poles on the imaginary axis, so the integral is $2\pi i$ times $B_n/n!$ plus the residues $(2\pi m i)^{-n}$ . On the other hand the integrand decays like $\lvert t\rvert^{-n}$ on the circle while the perimeter grows only linearly, so the integral itself tends to $0$ once $n\ge2$ . The two statements together say that the Taylor coefficient is minus the sum of all the other residues.

The last step is bookkeeping in the sum $\sum_{m\ne0}(2\pi m i)^{-n}$ : pairing $m$ with $-m$ makes it vanish for odd $n$ and doubles it for even $n$ , and for $n=2k$ the factor $i^{-2k}=(-1)^{k}$ comes out, leaving exactly $\zeta(2k)$ .

The one quantitative point needing care is the circle. The radius must be chosen to keep the contour away from the poles — the odd multiples $(2N+1)\pi$ do that, sitting half a period from every pole — and then $1/(e^t-1)$ is bounded on every such circle by a constant that does not depend on $N$ .

ProofAdaptedAdapted from Arfken, Weber & Harris (2013), §12.2, Eqs. (12.37)-(12.38) and Fig. 12.1, pp. 563-564

Fix $n\ge2$ and set

g(t) \;=\; \frac{t^{-n}}{e^{t}-1} .

Step 1: the poles of $g$ and the residues away from the origin. The zeros of $e^{t}-1$ are the points $t = 2\pi m i$ , $m \in \mathbb{Z}$ , and each is simple because $\frac{d}{dt}(e^t-1) = e^{t} = 1 \ne 0$ there. For $m \ne 0$ the factor $t^{-n}$ is analytic at $2\pi m i$ , so $g$ has a simple pole there with

\operatorname*{Res}_{t=2\pi m i} g \;=\; \frac{(2\pi m i)^{-n}}{e^{2\pi m i}} \;=\; \frac{1}{(2\pi m i)^{n}} ,

by the quotient rule for simple poles (The residue theorem). Apart from these, $g$ is analytic except at $t=0$ .

Step 2: the residue at the origin is the Bernoulli number. On $0<\lvert t\rvert<2\pi$ the definition of the $B_k$ gives the Laurent development

g(t) \;=\; t^{-n-1}\cdot\frac{t}{e^{t}-1} \;=\; \sum_{k=0}^{\infty}\frac{B_{k}}{k!}\,t^{\,k-n-1} ,

convergent on that punctured disc. The coefficient of $t^{-1}$ is the term $k=n$ , so by uniqueness of the Laurent expansion (Laurent series)

\operatorname*{Res}_{t=0} g \;=\; \frac{B_{n}}{n!} .

Step 3: a lower bound for $\lvert e^{t}-1\rvert$ on the circles $\lvert t\rvert = R_{N}$ . Put $R_{N} = (2N+1)\pi$ for $N\ge1$ . We claim

\lvert e^{t}-1\rvert \;\ge\; \tfrac12 \qquad\text{for all } \lvert t\rvert = R_{N},\ N\ge1 .

Write $t = x+iy$ and split on the size of $x$ .

If $\lvert x\rvert\ge1$ : the reverse triangle inequality gives $\lvert e^{t}-1\rvert \ge \bigl\lvert\,\lvert e^{t}\rvert-1\,\bigr\rvert = \lvert e^{x}-1\rvert$ , which is at least $e-1$ when $x\ge1$ and at least $1-e^{-1} = 0.6321$ when $x\le-1$ .
If $\lvert x\rvert\le1$ : then $\lvert y\rvert = \sqrt{R_{N}^{2}-x^{2}} \ge \sqrt{R_{N}^{2}-1}$ , so $0 \;\le\; R_{N}-\lvert y\rvert \;\le\; R_{N}-\sqrt{R_{N}^{2}-1} \;=\; \frac{1}{R_{N}+\sqrt{R_{N}^{2}-1}} \;\le\; \frac{1}{R_{N}} \;\le\; \frac{1}{3\pi} \;<\; 0.107 .$ Since $R_{N} = 2\pi N+\pi$ , writing $\eta = \lvert y\rvert - 2\pi N$ puts $\eta$ in the interval $[\pi-0.107,\ \pi]$ . The function $e^{t}-1$ has period $2\pi i$ , so $\lvert e^{t}-1\rvert = \lvert e^{x+i\eta}-1\rvert$ and $\lvert e^{x+i\eta}-1\rvert^{2} \;=\; e^{2x}-2e^{x}\cos\eta+1 \;\ge\; e^{2x}+2e^{x}\cos(0.107)+1 \;\ge\; e^{-2}+1.98\,e^{-1}+1 \;>\; 1.8 ,$ using $\cos\eta \le \cos(\pi-0.107) = -\cos(0.107) < -0.994$ in the middle step and $x \ge -1$ with monotonicity in $x$ at the end. Hence $\lvert e^{t}-1\rvert > 1.3$ here.

Both cases give at least $\tfrac12$ , and in particular no pole of $g$ lies on any circle $\lvert t\rvert = R_{N}$ .

Step 4: the contour integral vanishes in the limit. The poles of $g$ inside $\lvert t\rvert<R_{N}$ are $t=0$ and $t = 2\pi m i$ with $1 \le \lvert m\rvert \le N$ , since $\lvert 2\pi m\rvert < (2N+1)\pi$ is exactly $\lvert m\rvert\le N$ . The residue theorem on the positively oriented circle $C_{N} = \{\lvert t\rvert = R_{N}\}$ therefore gives

\frac{1}{2\pi i}\oint_{C_{N}} g(t)\,dt \;=\; \frac{B_{n}}{n!} \;+\; \sum_{1\le\lvert m\rvert\le N}\frac{1}{(2\pi m i)^{n}} .

By Step 3 the integrand obeys $\lvert g(t)\rvert \le 2R_{N}^{-n}$ on $C_{N}$ , so the estimate for contour integrals gives

\left\lvert\frac{1}{2\pi i}\oint_{C_{N}}g\,dt\right\rvert \;\le\; \frac{1}{2\pi}\cdot 2\pi R_{N}\cdot\frac{2}{R_{N}^{\,n}} \;=\; \frac{2}{R_{N}^{\,n-1}} \;\xrightarrow[N\to\infty]{}\;0 ,

because $n\ge2$ makes the exponent $n-1$ positive. The remaining sum converges absolutely, its terms having modulus $(2\pi\lvert m\rvert)^{-n}$ with $n\ge2$ , so letting $N\to\infty$ yields the displayed identity

\frac{B_{n}}{n!} \;=\; -\sum_{m\ne0}\frac{1}{(2\pi m i)^{n}} .

Step 5: separate the parities. Pair the terms $m$ and $-m$ . For odd $n$ , $(-m)^{-n} = -m^{-n}$ , so each pair cancels and the sum is $0$ ; hence $B_{n}=0$ for odd $n\ge3$ , which is claim 1. For $n = 2k$ the two terms of a pair are equal, and $i^{-2k} = (i^{2})^{-k} = (-1)^{k}$ , so

\frac{B_{2k}}{(2k)!} \;=\; -\,2\sum_{m=1}^{\infty}\frac{1}{(2\pi m)^{2k}\,i^{2k}} \;=\; \frac{(-1)^{k+1}\,2}{(2\pi)^{2k}}\sum_{m=1}^{\infty}\frac{1}{m^{2k}} \;=\; (-1)^{k+1}\,\frac{2\,\zeta(2k)}{(2\pi)^{2k}} ,

the last step by the definition of $\zeta$ at the real argument $2k>1$ . Multiplying by $(2k)!$ gives claim 2 in Arfken's form, and dividing instead by the prefactor gives the inverted form. Claim 3 is immediate: every term of $\zeta(2k)$ is positive, so $\operatorname{sign}B_{2k} = (-1)^{k+1}$ , and $1 < \zeta(2k) \le \zeta(2) = \pi^{2}/6$ for $k\ge1$ pins $\lvert B_{2k}\rvert$ between $2(2k)!/(2\pi)^{2k}$ and $\tfrac{\pi^{2}}{6}$ times that.

The sample values follow from the table of $B_{2k}$ : with $k=1$ , $(2\pi)^{2}/(2\cdot2!) = \pi^{2}$ and $B_2 = 1/6$ give $\zeta(2)=\pi^{2}/6$ ; with $k=2$ , $(2\pi)^{4}/(2\cdot4!) = \pi^{4}/3$ and $B_4 = -1/30$ with the sign $(-1)^{3}$ give $\zeta(4) = \pi^{4}/90$ ; $k=3$ is done in the example below. $\square$

Remark(Where the hypotheses are used)AdaptedAdapted from Arfken, Weber & Harris (2013), §12.2, pp. 560-564; Stein & Shakarchi (2003), Ch. 3, Problem 4(c), p. 98

$s>1$ in the definition. This is what makes the defining series converge, and it is used at the very end of the proof, where $\sum m^{-2k}$ is identified as $\zeta(2k)$ : the identification is legitimate only because $2k\ge2>1$ . The proof produces no information at odd arguments and none at $s\le1$ .
$n\ge2$ in the main identity. It enters twice: to make the large-circle bound $2R_N^{1-n}$ tend to zero, and to make the residue sum converge absolutely. For $n=1$ the identity is false — its right-hand side is a divergent symmetric sum, while $B_1 = -\tfrac12$ — and this is exactly why $B_1$ is the one odd-index Bernoulli number that does not vanish. It has to be computed separately from the generating function, as in the definition above.
The radius $R_N = (2N+1)\pi$ . An odd multiple of $\pi$ sits half a period away from every pole $2\pi m i$ , which is what Step 3 converts into a bound independent of $N$ . A circle of radius $2N\pi$ would pass through two poles and the argument would not even be stated. Arfken, Weber and Harris (Arfken, Weber & Harris (2013), Fig. 12.1, p. 563) use a different contour for the same purpose — one that detours around the origin and runs along the real axis — and must then argue that the two real-axis segments cancel and that the large arc is negligible; the single growing circle avoids both discussions.
Simplicity of the poles of $1/(e^t-1)$ . Used in Step 1 to get the residue $1$ . At $t=0$ the pole of $g$ has order $n+1$ , so Step 2 reads the residue off the Laurent series instead of using a residue formula — this is the step that makes the Bernoulli numbers appear at all.
The radius of convergence $2\pi$ is a real restriction. The series $\sum B_nt^n/n!$ is useless beyond $\lvert t\rvert = 2\pi$ , which is why the expansions of $\cot$ , $\tan$ and $\coth$ in Bernoulli numbers all carry a condition like $\lvert t\rvert<2\pi$ or $\lvert t\rvert<\pi$ . Claim 3 is the same fact seen from the coefficients: $\lvert B_{2k}\rvert^{1/2k}/(2k)!^{1/2k} \to 1/2\pi$ .
An alternative route. Stein and Shakarchi (Stein & Shakarchi (2003), Ch. 3, Problem 4(c)-(e), pp. 97-98) reach the same formula without contours: the generating function gives $z\cot z = 1+\sum_{k\ge1}(-1)^{k}2^{2k}B_{2k}z^{2k}/(2k)!$ , the partial-fraction expansion of the cotangent (Pole expansions and infinite products) gives $x\cot x = 1-2\sum_{m\ge1}\zeta(2m)x^{2m}/\pi^{2m}$ , and comparing coefficients of $x^{2m}$ does the rest. The two routes are the same computation in different dress, since the partial-fraction expansion of the cotangent is itself a statement about the poles at the multiples of $\pi$ . We follow the contour route because it needs only the residue theorem, which The residue theorem already provides.

Note(Three conventions for the Bernoulli numbers, and what lies outside

s>1

)AdaptedAdapted from Arfken, Weber & Harris (2013), §12.2, Eq. (12.39), p. 564; Ahlfors (1979), Ch. 5, §1.3, Exercise 4, p. 186; Arfken, Weber & Harris (2013), §13.5, Eqs. (13.67)-(13.72) and Exercise 13.5.10, pp. 629-633

Conventions. Arfken, Weber and Harris (Arfken, Weber & Harris (2013), p. 560) warn that formulas from different authors must not be combined without checking the definition, and the library bears them out. Three indexings occur in it:

ours, $t/(e^t-1)=\sum B_nt^n/n!$ , which is the convention of the Handbook of Mathematical Functions that Arfken, Weber and Harris adopt (p. 560);
the variant they record at Eq. (12.39), $\widetilde B_{k} = 2(2k)!\,\zeta(2k)/(2\pi)^{2k}$ , with the subscript halved and all signs positive, so that $\widetilde B_{k} = \lvert B_{2k}\rvert$ ;
Ahlfors (Ahlfors (1979), Ch. 5, §1.3, Exercise 4, p. 186), who writes the Laurent development of $(e^{z}-1)^{-1}$ at the origin as $z^{-1}-\tfrac12+\sum_{k\ge1}(-1)^{k-1}B_{k}^{\mathrm{A}}z^{2k-1}/(2k)!$ and notes that these $B^{\mathrm{A}}_{k}$ are all positive. Comparing with our Step 2 gives $B^{\mathrm{A}}_{k} = (-1)^{k-1}B_{2k} = \lvert B_{2k}\rvert$ , so Ahlfors' indexing agrees with the second variant.

So the only differences in this library are the halved subscript and the suppressed sign; with ours, $B_{2k}$ and $\zeta(2k)$ are related with the factor $(-1)^{k+1}$ , and with the other two the sign disappears. Our numerical values $B_2=\tfrac16$ , $B_4=-\tfrac1{30}$ , $B_6=\tfrac1{42}$ are on our convention throughout.

Outside the half line. Everything above is confined to real $s>1$ , which is this unit's scope, but the function does not stop there. Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §13.5, pp. 629-633) continue $\zeta$ to the whole plane by the reflection formula $\Gamma(\tfrac{z}{2})\pi^{-z/2}\zeta(z) = \Gamma(\tfrac{1-z}{2})\pi^{-(1-z)/2}\zeta(1-z)$ (Eq. (13.67)) together with the alternating series $\zeta(z) = (1-2^{1-z})^{-1}\sum_{n\ge1}(-1)^{n-1}n^{-z}$ , which converges on $\operatorname{Re}z>0$ and so covers the critical strip; the continued function is analytic except for a simple pole at $z=1$ of residue $1$ (Eq. (13.69) and Exercise 13.5.10). They also derive Euler's prime product $\zeta(s) = \prod_{p\text{ prime}}(1-p^{-s})^{-1}$ (Eq. (13.72)), which the application below meets again in probabilistic clothing. The machinery is that of Analytic continuation; none of it is needed here.

Method(Values of

\zeta

at integer arguments)Original

Even argument $2k$ , closed form. Get $B_{2k}$ from the recursion $N-1 = \sum_{n=1}^{N-1}\binom{2N}{2n}B_{2n}$ with $N=k+1$ , which involves only the earlier even-index numbers; then $\zeta(2k) \;=\; (-1)^{k+1}\frac{(2\pi)^{2k}}{2\,(2k)!}B_{2k} .$ The arithmetic is exact and rational apart from the power of $\pi$ .
Odd argument, or non-integer $s$ . No closed form is available, and summing the series directly converges only algebraically: its terms decay like $n^{-s}$ $n^{- s}$ and the tail after $N$ $N$ terms is about $N^{1-s}/(s-1)$ $N^{1 - s} / (s - 1)$ . At $s=3$ $s = 3$ that tail is about $1/(2N^{2})$ $1/ (2 N^{2})$ , so a few dozen terms are enough for four figures — $N=32$ $N = 32$ puts the partial sum within $5\times10^{-4}$ $5 \times 1 0^{- 4}$ of $\zeta(3)=1.2020569\ldots$ $ζ (3) = 1.2020569 \dots$ — but every further decimal digit costs a factor $\sqrt{10}\approx3.2$ $10 \approx 3.2$ in $N$ $N$ : $N=707$ $N = 707$ for $10^{-6}$ $1 0^{- 6}$ and $N\approx2.2\times10^{4}$ $N \approx 2.2 \times 1 0^{4}$ for $10^{-9}$ $1 0^{- 9}$ . The cost worsens without limit as $s\downarrow1$ $s ↓ 1$ , where the exponent $s-1$ $s - 1$ collapses. For many digits, or for $s$ $s$ near $1$ $1$ , either
- use Euler-Maclaurin summation (Euler-Maclaurin summation), which replaces the tail $\sum_{n\ge N}n^{-s}$ by $\frac{N^{1-s}}{s-1}+\frac{N^{-s}}{2}+\cdots$ , the further corrections carrying the same $B_{2k}$ as above, and gives many digits from a handful of terms; or
- use the alternating series $\zeta(s) = (1-2^{1-s})^{-1}\sum_{n\ge1}(-1)^{n-1}n^{-s}$ of Arfken, Weber and Harris (Arfken, Weber & Harris (2013), Eq. (13.68), p. 629), whose partial sums bracket the answer.
A consistency check worth doing. Any claimed $\zeta(2k)$ should satisfy $\zeta(2k)>1$ , decrease in $k$ , and sit above the two-term lower bound $1+2^{-2k}$ . The closed forms do: $\pi^{2}/6 = 1.644934$ , $\pi^{4}/90 = 1.082323$ , $\pi^{6}/945 = 1.017343$ .

Example(

\zeta(6)=\pi^{6}/945

, from the recursion to the number)Original

Step 1: $B_{6}$ from the recursion. Take $N=4$ in $N-1=\sum_{n=1}^{N-1}\binom{2N}{2n}B_{2n}$ :

3 \;=\; \binom{8}{2}B_{2}+\binom{8}{4}B_{4}+\binom{8}{6}B_{6} \;=\; 28\cdot\tfrac16 \;+\; 70\cdot\left(-\tfrac1{30}\right) \;+\; 28\,B_{6} .

The first two terms are $\tfrac{14}{3}$ and $-\tfrac73$ , whose sum is $\tfrac73$ , so $28B_{6} = 3-\tfrac73 = \tfrac23$ and

B_{6} \;=\; \frac{2}{3\cdot28} \;=\; \frac{1}{42} .

(The same recursion at $N=2$ and $N=3$ returns $B_2=\tfrac16$ and $B_4=-\tfrac1{30}$ , and at $N=5,6$ it returns $B_8=-\tfrac1{30}$ , $B_{10}=\tfrac5{66}$ , matching Table 12.2 of Arfken, Weber and Harris, p. 562.)

Step 2: insert it. With $k=3$ the sign is $(-1)^{4}=+1$ , the factorial is $6! = 720$ , and $(2\pi)^{6} = 64\pi^{6}$ :

\zeta(6) \;=\; \frac{64\pi^{6}}{2\cdot720}\cdot\frac{1}{42} \;=\; \frac{64\pi^{6}}{60480} \;=\; \frac{\pi^{6}}{945} .

The reduction is $60480/64 = 945$ . Numerically $\pi^{6} = 961.389194$ , so

\zeta(6) \;=\; 1.0173430620 .

Step 3: check it against the series. Summing $\sum_{n\le N}n^{-6}$ directly:

$N$	$\sum_{n\le N}n^{-6}$	$\zeta(6)-\sum$	tail bound $\tfrac1{5N^{5}}$
$10$	$1.01734151244$	$1.5495\times10^{-6}$	$2.0\times10^{-6}$
$100$	$1.01734306196$	$1.95\times10^{-11}$	$2.0\times10^{-11}$
$1000$	$1.01734306198$	$2.0\times10^{-16}$	$2.0\times10^{-16}$

The tail bound is the integral comparison $\sum_{n>N}n^{-6} \le \int_{N}^{\infty}x^{-6}dx = 1/(5N^{5})$ , the same comparison that proved convergence in the definition; the observed errors sit just below it, as they must. Note the cost: $1000$ terms of the series buy sixteen digits, while Step 2 gives them all exactly. That contrast is the whole point of the closed form, and it is also why the method above sends odd arguments to Euler-Maclaurin rather than to the series.

Application(Blackbody radiation and the Stefan-Boltzmann constant)AdaptedAdapted from Arfken, Weber & Harris (2013), §13.5, Exercises 13.5.3, 13.5.5, 13.5.6, 13.5.7 and 13.5.8, pp. 632-633

The energy density of thermal radiation in a cavity, integrated over frequency with the Planck distribution, is

u \;=\; \frac{8\pi k^{4}T^{4}}{c^{3}h^{3}}\int_{0}^{\infty}\frac{x^{3}}{e^{x}-1}\,dx ,

the set-up of Arfken, Weber and Harris (Arfken, Weber & Harris (2013), Exercise 13.5.5, p. 632), where $x = h\nu/kT$ . The integral is a zeta value. Expand the integrand as a geometric series in $e^{-x}$ , which is legitimate term by term because every term is positive and the partial sums increase ( $\sum_{n\le M}$ of positive measurable terms, so monotone convergence applies, Integrating and differentiating series term by term):

\int_{0}^{\infty}\frac{x^{3}}{e^{x}-1}dx \;=\; \int_{0}^{\infty}x^{3}\sum_{n=1}^{\infty}e^{-nx}\,dx \;=\; \sum_{n=1}^{\infty}\int_{0}^{\infty}x^{3}e^{-nx}dx \;=\; \sum_{n=1}^{\infty}\frac{3!}{n^{4}} \;=\; 3!\,\zeta(4) ,

the inner integral being $\Gamma(4)/n^{4}$ after $x\mapsto x/n$ (The Gamma function). This is the general identity $\int_{0}^{\infty}x^{s}(e^{x}-1)^{-1}dx = \Gamma(s+1)\zeta(s+1)$ of Arfken, Weber and Harris, Exercise 13.5.6. Now Euler's formula supplies the closed value: $3!\,\zeta(4) = 6\cdot\pi^{4}/90 = \pi^{4}/15 = 6.49393940$ .

So $u = 8\pi^{5}k^{4}T^{4}/(15c^{3}h^{3})$ , and with the radiated flux $\sigma T^{4} = cu/4$ ,

\sigma \;=\; \frac{2\pi^{5}k^{4}}{15\,c^{2}h^{3}} .

Putting in $k = 1.380649\times10^{-23}\,\mathrm{J\,K^{-1}}$ , $h = 6.62607015\times10^{-34}\,\mathrm{J\,s}$ and $c = 2.99792458\times10^{8}\,\mathrm{m\,s^{-1}}$ gives

\sigma \;=\; 5.6703744\times10^{-8}\ \mathrm{W\,m^{-2}\,K^{-4}} ,

the measured Stefan-Boltzmann constant. The $\pi^{5}$ in that number is $\pi^{4}$ from $\zeta(4)$ — that is, from $B_{4} = -1/30$ — times one $\pi$ from the solid-angle factor. The same pattern recurs: Arfken, Weber and Harris's Exercise 13.5.3 puts $\zeta(5)$ into the low-temperature resistivity of a monovalent metal, and Exercise 13.5.8 puts $7\pi^{5}/30$ into the neutrino energy density of the early universe. That last factor is the Fermi counterpart of the computation above: the sign in the denominator changes to $e^{x}+1$ , the geometric series becomes alternating, and Arfken, Weber and Harris's Exercise 13.5.7 gives $\int_{0}^{\infty}x^{s}(e^{x}+1)^{-1}dx = \Gamma(s+1)(1-2^{-s})\zeta(s+1)$ , so that $\int_{0}^{\infty}x^{3}(e^{x}+1)^{-1}dx = 6\cdot\tfrac78\cdot\tfrac{\pi^{4}}{90} = \tfrac{7\pi^{4}}{120}$ , which is $\tfrac78$ of the Bose value. Where a Bose or Fermi distribution is integrated against a power, a zeta value appears, and for an even argument this unit evaluates it.

Application(A random integer: Euler's product as independence)AdaptedAdapted from Grimmett & Stirzaker (2020), Problems 3.11.31 and 3.11.52, pp. 93, 97

Because $\zeta(\beta)<\infty$ for $\beta>1$ , the weights $m^{-\beta}$ can be normalised into a probability distribution on $\mathbb{N}$ :

\mathbb{P}(M=m) \;=\; \frac{m^{-\beta}}{\zeta(\beta)} , \qquad m = 1,2,3,\dots

Grimmett and Stirzaker (Grimmett & Stirzaker (2020), Problem 3.11.31, p. 93) construct this $M$ as $\prod_{i}p_{i}^{N(i)}$ over the primes $p_1,p_2,\dots$ , with the exponents $N(i)$ independent and geometric with parameter $1-p_{i}^{-\beta}$ ; the normalising constant then has the two expressions $\zeta(\beta)^{-1} = \prod_{i}(1-p_{i}^{-\beta})$ , which is Euler's product. Read the other way (their Problem 3.11.52, p. 97), the events $E_{p} = \{p \text{ divides } M\}$ are independent across primes with $\mathbb{P}(E_{p}) = p^{-\beta}$ , and Euler's product is just the statement that $\mathbb{P}(M=1) = \prod_{p}(1-p^{-\beta})$ .

Independence makes a divisibility question into a product. Grimmett and Stirzaker's part (c): $M$ is divisible by no perfect square above $1$ exactly when no $p^{2}$ divides it, and these events are independent with probability $p^{-2\beta}$ each, so

\mathbb{P}(M \text{ square-free}) \;=\; \prod_{p}\bigl(1-p^{-2\beta}\bigr) \;=\; \frac{1}{\zeta(2\beta)} .

At $\beta=2$ the answer is a number this unit can produce in closed form, because $2\beta = 4$ is even:

\mathbb{P}(M\text{ square-free}) \;=\; \frac{1}{\zeta(4)} \;=\; \frac{90}{\pi^{4}} \;=\; 0.9239384 .

The classical companion, with the uniform distribution on $\{1,\dots,n\}$ and $n\to\infty$ , is $1/\zeta(2) = 6/\pi^{2} = 0.6079271$ for the probability that two integers are coprime. Both numbers are closed-form only because the argument is even.

Definition

Error function and incomplete Gamma function

Adapted

Add to PDF

Defines the incomplete Gamma functions $\gamma(a,x)$ and $\Gamma(a,x)$ and the error functions $\operatorname{erf}$ and $\operatorname{erfc}$ as the case $a=\tfrac12$ , then gives both of the expansions a computation needs — a convergent power series for small argument and an asymptotic series with a rigorous remainder bound for large argument.

Prerequisites: The Gamma function, Asymptotic series, Power series, Integrating and differentiating series term by term, Uniform convergence, Gaussian integrals

Definition(The two incomplete Gamma functions)AdaptedAdapted from Arfken, Weber & Harris (2013), §13.6, Eqs. (13.73)-(13.74), p. 633; Riley, Hobson & Bence (2006), §18.12.3, Eq. (18.166), p. 639

Cutting Euler's integral for $\Gamma(a)$ (The Gamma function) at a point $x>0$ defines the lower and upper incomplete Gamma functions

\gamma(a,x) \;=\; \int_{0}^{x}e^{-t}t^{\,a-1}\,dt , \qquad \Gamma(a,x) \;=\; \int_{x}^{\infty}e^{-t}t^{\,a-1}\,dt .

The first needs $\operatorname{Re}a>0$ for convergence at the origin, where the integrand has modulus $t^{\operatorname{Re}a-1}e^{-t}$ ; the second converges for every $a$ once $x>0$ , because the only issue there is the decay at infinity, which $e^{-t}$ supplies. Adding them restores the whole integral:

\gamma(a,x)+\Gamma(a,x) \;=\; \Gamma(a) \qquad (\operatorname{Re}a>0) .

Because of that identity the two carry the same information, so which one a calculation is written in terms of is settled by the problem at hand and not by any difference in content; both Arfken, Weber and Harris (Arfken, Weber & Harris (2013), p. 633) and Riley, Hobson and Bence (Riley, Hobson & Bence (2006), p. 639) make the point. Riley, Hobson and Bence also record the normalised pair

P(a,x) \;=\; \frac{\gamma(a,x)}{\Gamma(a)} , \qquad Q(a,x) \;=\; \frac{\Gamma(a,x)}{\Gamma(a)} \;=\; 1-P(a,x) ,

which are the ones a statistician meets, since $P(a,x)$ is a cumulative distribution function in $x$ .

Definition(Error function and complementary error function)AdaptedAdapted from Arfken, Weber & Harris (2013), §13.6, Eqs. (13.90) and (13.93), pp. 637-638; Riley, Hobson & Bence (2006), §18.12.4, Eqs. (18.167)-(18.168), p. 640; Strauss (2008), §2.4, Eq. (10), p. 51

The error function and the complementary error function are

\operatorname{erf}(x) \;=\; \frac{2}{\sqrt{\pi}}\int_{0}^{x}e^{-u^{2}}\,du , \qquad \operatorname{erfc}(x) \;=\; 1-\operatorname{erf}(x) \;=\; \frac{2}{\sqrt{\pi}}\int_{x}^{\infty}e^{-u^{2}}\,du .

The factor $2/\sqrt{\pi}$ is chosen so that $\operatorname{erf}(\infty)=1$ , which is the Gaussian integral $\int_{0}^{\infty}e^{-u^{2}}du = \tfrac12\sqrt{\pi}$ (Gaussian integrals); this is the normalisation of Arfken, Weber and Harris (Arfken, Weber & Harris (2013), Table 1.2, p. 66, and Eq. (13.90), p. 637), Riley, Hobson and Bence (Riley, Hobson & Bence (2006), Eq. (18.167), p. 640) and Strauss (Strauss (2008), Eq. (10), p. 51), and the one we use throughout. Immediately from the definition: $\operatorname{erf}(0)=0$ , $\operatorname{erf}(-x) = -\operatorname{erf}(x)$ since the integrand is even, $\operatorname{erfc}(0)=1$ , and

\frac{d}{dx}\operatorname{erf}(x) \;=\; \frac{2}{\sqrt{\pi}}e^{-x^{2}} ,

so $\operatorname{erf}$ is strictly increasing, extends to an entire function of a complex argument, and satisfies the differential equation $y''+2xy'=0$ noted by Hassani (Hassani (2013), §14.6, p. 437) — on our normalisation with the data $y(0)=0$ , $y'(0)=2/\sqrt{\pi}$ .

The error function is the case $a=\tfrac12$ . Substituting $t=u^{2}$ in the lower incomplete Gamma integral, so that $dt = 2u\,du$ and $t^{-1/2}=u^{-1}$ ,

\gamma\bigl(\tfrac12,x^{2}\bigr) \;=\; \int_{0}^{x^{2}}e^{-t}t^{-1/2}dt \;=\; 2\int_{0}^{x}e^{-u^{2}}du \;=\; \sqrt{\pi}\,\operatorname{erf}(x) ,

and likewise $\Gamma(\tfrac12,x^{2}) = \sqrt{\pi}\,\operatorname{erfc}(x)$ . These are the relations

\operatorname{erf}(x) = \pi^{-1/2}\gamma\bigl(\tfrac12,x^{2}\bigr), \qquad \operatorname{erfc}(x) = \pi^{-1/2}\Gamma\bigl(\tfrac12,x^{2}\bigr)

of Arfken, Weber and Harris (Arfken, Weber & Harris (2013), Eq. (13.93), p. 638) and Riley, Hobson and Bence (Riley, Hobson & Bence (2006), Eqs. (18.167)-(18.168), p. 640); the identity $\operatorname{erf}+\operatorname{erfc}=1$ is the splitting of the previous definition divided by $\Gamma(\tfrac12)=\sqrt{\pi}$ . Everything proved below about $\gamma(a,x)$ and $\Gamma(a,x)$ therefore specialises to the error functions, and that is how we obtain their expansions.

The Gaussian distribution function. Riley, Hobson and Bence (p. 640) record the form in which the error function is met in probability: for the standard normal distribution,

\Phi(x) \;=\; \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{x}e^{-y^{2}/2}\,dy \;=\; \frac12+\frac12\operatorname{erf}\!\left(\frac{x}{\sqrt2}\right) .

The $\sqrt2$ is the whole content of the difference between the two conventions: one integrates $e^{-u^{2}}$ , the other $e^{-y^{2}/2}$ . As a check, $\Phi(1) = \tfrac12+\tfrac12\operatorname{erf}(0.7071068) = \tfrac12+\tfrac12(0.6826895) = 0.8413447$ , the familiar one-sigma figure.

Definition(Series, recurrences and large-argument expansions)AdaptedAdapted from Arfken, Weber & Harris (2013), §13.6, Eqs. (13.75)-(13.76), pp. 633-634, Eqs. (13.91)-(13.92), p. 637, and Exercises 13.6.3-13.6.4, p. 638; Riley, Hobson & Bence (2006), §18.12.3, the integration-by-parts recurrence and integer closed form, p. 639

Let $x>0$ and $a$ be real with $a>0$ where $\gamma$ appears.

Recurrences. $\displaystyle \Gamma(a+1,x) = a\,\Gamma(a,x)+x^{a}e^{-x}$ for every real $a$ , and $\displaystyle \gamma(a+1,x) = a\,\gamma(a,x)-x^{a}e^{-x}$ for $a>0$ .
Convergent series for the lower function. For every $x>0$ , $\gamma(a,x) \;=\; x^{a}\sum_{n=0}^{\infty}\frac{(-1)^{n}x^{n}}{n!\,(a+n)} ,$ the series converging for all $x$ (it is the term-by-term integral of the exponential series). In particular $\operatorname{erf}(x) \;=\; \frac{2}{\sqrt{\pi}}\sum_{n=0}^{\infty}\frac{(-1)^{n}x^{2n+1}}{(2n+1)\,n!} \;=\; \frac{2}{\sqrt{\pi}}\left(x-\frac{x^{3}}{3}+\frac{x^{5}}{10}-\frac{x^{7}}{42}+\cdots\right) ,$ valid for every real (indeed complex) $x$ .
Exact truncation with remainder. Put $c_{0}=1$ and $c_{m} = (a-1)(a-2)\cdots(a-m)$ for $m\ge1$ . Then for every $N\ge1$ $\Gamma(a,x) \;=\; e^{-x}x^{\,a-1}\sum_{m=0}^{N-1}\frac{c_{m}}{x^{m}} \;+\; c_{N}\,\Gamma(a-N,x) ,$ an identity, not an approximation. If moreover $N\ge a-1$ , the remainder obeys $\bigl\lvert c_{N}\Gamma(a-N,x)\bigr\rvert \;\le\; \lvert c_{N}\rvert\,e^{-x}x^{\,a-1-N} ,$ which is the modulus of the first omitted term.
Asymptotic expansion. Consequently, as $x\to\infty$ with $a$ fixed, $\Gamma(a,x) \;\sim\; e^{-x}x^{\,a-1}\sum_{m=0}^{\infty}\frac{(a-1)(a-2)\cdots(a-m)}{x^{m}}$ in the sense of Asymptotic series: each truncation has error of smaller order than its last retained term. This is the expansion of Arfken, Weber and Harris (Arfken, Weber & Harris (2013), Eq. (13.76), p. 634, and Exercise 13.6.4, p. 638), whose coefficient $\Gamma(a)/\Gamma(a-m)$ is the same product; they state it for non-integral $a$ . For such $a$ the series diverges for every $x$ . When $a$ is a positive integer it is not a series at all but a terminating sum, which is part 5.
Integer $a$ : the expansion terminates and is exact. If $a=n\in\mathbb{N}$ then $c_{n}=0$ , so part 3 with $N=n$ has no remainder and $\Gamma(n,x) \;=\; (n-1)!\,e^{-x}\sum_{s=0}^{n-1}\frac{x^{s}}{s!} , \qquad \gamma(n,x) \;=\; (n-1)!\left(1-e^{-x}\sum_{s=0}^{n-1}\frac{x^{s}}{s!}\right) ,$ the closed forms of Arfken, Weber and Harris (Arfken, Weber & Harris (2013), Eq. (13.75), p. 633) and Riley, Hobson and Bence (Riley, Hobson & Bence (2006), p. 639).
The error functions. Taking $a=\tfrac12$ and replacing $x$ by $x^{2}$ in parts 3 and 4: for $x>0$ and $N\ge1$ , $\operatorname{erfc}(x) \;=\; \frac{e^{-x^{2}}}{x\sqrt{\pi}}\left(\sum_{m=0}^{N-1}(-1)^{m}\frac{(2m-1)!!}{2^{m}x^{2m}}\right) \;+\; R_{N}(x) , \qquad \lvert R_{N}(x)\rvert \;\le\; \frac{e^{-x^{2}}}{x\sqrt{\pi}}\cdot\frac{(2N-1)!!}{2^{N}x^{2N}} ,$ with $(-1)!!=1$ and $(2m-1)!! = 1\cdot3\cdots(2m-1)$ . Hence $\operatorname{erfc}(x) \;\sim\; \frac{e^{-x^{2}}}{x\sqrt{\pi}}\left(1-\frac{1}{2x^{2}}+\frac{3}{4x^{4}}-\frac{15}{8x^{6}}+\cdots\right) \qquad (x\to\infty) ,$ which is the expansion Arfken, Weber and Harris (Arfken, Weber & Harris (2013), Eq. (13.92), p. 637) give in the equivalent form $\operatorname{erf}(x) \approx 1-\pi^{-1/2}x^{-1}e^{-x^{2}}\sum_{m}(-1)^{m}(2m-1)!!/(2^{m}x^{2m})$ . The terms alternate in sign, so the bound on $R_{N}$ says the truncation error is no larger in modulus than the first term dropped, and successive truncations bracket the answer.

Idea of the proofOriginal

Two expansions, one for each end of the range, and each is got by expanding a different factor of the integrand.

For small $x$ the exponential is the harmless factor: expand $e^{-t}$ in its everywhere-convergent power series and integrate term by term over the bounded interval $[0,x]$ . Each monomial integrates to a power of $x$ , and the result converges for every $x$ — it is an honest Taylor series. The reason this only works for the lower function is that the upper one integrates over an unbounded interval, where the exponential series converges nowhere uniformly.

For large $x$ the power $t^{a-1}$ is the harmless factor: on $[x,\infty)$ it varies slowly compared with the collapse of $e^{-t}$ , so almost all of the integral comes from a neighbourhood of the left endpoint. Integration by parts is the systematic way to exploit that: it peels off $x^{a-1}e^{-x}$ , the value at the endpoint, and leaves an integral of the same shape with the exponent lowered by one. Each repetition peels off a smaller term and lowers the exponent again. Nothing has been approximated so far — the identity in part 3 is exact for every $N$ .

The approximation enters only when we stop and discard the leftover integral, and the single fact that makes the discard controllable is that $t^{a-N-1}$ is decreasing once $N > a-1$ : then it is largest at the endpoint, so the leftover integral is at most $x^{a-N-1}\int_{x}^{\infty}e^{-t}dt$ . That is precisely the size of the next term, which is the classical statement that the error is no worse than the first term omitted.

Unless the peeling stops of its own accord — which happens exactly when $a$ is a positive integer, for then some factor $a-m$ is $0$ — the series cannot converge, because the exponents being lowered eventually make $\lvert c_{m}\rvert$ grow like $m!$ while $x^{-m}$ shrinks only geometrically. So for each $x$ there is a best place to stop — roughly where consecutive terms stop decreasing — and that best accuracy is all the expansion will ever give. The error functions inherit all of this at $a=\tfrac12$ , and the one accident of $a=\tfrac12$ is that the coefficients are the double factorials with strictly alternating signs.

ProofAdaptedAdapted from Arfken, Weber & Harris (2013), §13.6, Exercise 13.6.3 and Exercise 13.6.4, p. 638, and Eq. (13.76), p. 634; Riley, Hobson & Bence (2006), §18.12.3, the worked integration by parts, p. 639

Part 1 (recurrences). Integrate by parts with $u=t^{a}$ , $dv = e^{-t}dt$ . On $[x,\infty)$ ,

\Gamma(a+1,x) \;=\; \int_{x}^{\infty}t^{a}e^{-t}dt \;=\; \Bigl[-t^{a}e^{-t}\Bigr]_{x}^{\infty}+a\int_{x}^{\infty}t^{a-1}e^{-t}dt \;=\; x^{a}e^{-x}+a\,\Gamma(a,x) ,

the boundary term at infinity vanishing because $t^{a}e^{-t}\to0$ for every fixed $a$ . On $[0,x]$ the same computation gives the boundary term $-x^{a}e^{-x}$ at the upper end and $0$ at the lower end, the latter requiring $a>0$ so that $t^{a}\to0$ :

\gamma(a+1,x) \;=\; \Bigl[-t^{a}e^{-t}\Bigr]_{0}^{x}+a\int_{0}^{x}t^{a-1}e^{-t}dt \;=\; -x^{a}e^{-x}+a\,\gamma(a,x) .

Part 2 (convergent series). Fix $x>0$ . The exponential series $e^{-t} = \sum_{n\ge0}(-1)^{n}t^{n}/n!$ converges uniformly on the compact interval $[0,x]$ , by the Weierstrass $M$ -test with $M_{n}=x^{n}/n!$ (Uniform convergence). Multiplying a uniformly convergent series by the fixed function $t^{a-1}$ preserves uniform convergence on $[\varepsilon,x]$ for each $\varepsilon>0$ , and the resulting series of integrals converges absolutely because $\int_{0}^{x}t^{a+n-1}dt = x^{a+n}/(a+n)$ is summable; so we may integrate term by term (Integrating and differentiating series term by term) and let $\varepsilon\to0$ :

\gamma(a,x) \;=\; \sum_{n=0}^{\infty}\frac{(-1)^{n}}{n!}\int_{0}^{x}t^{\,a+n-1}dt \;=\; \sum_{n=0}^{\infty}\frac{(-1)^{n}}{n!}\cdot\frac{x^{\,a+n}}{a+n} \;=\; x^{a}\sum_{n=0}^{\infty}\frac{(-1)^{n}x^{n}}{n!\,(a+n)} .

For the error function put $a=\tfrac12$ and replace $x$ by $x^{2}$ : the general term becomes $(-1)^{n}x^{2n+1}\big/\bigl(n!(n+\tfrac12)\bigr) = 2(-1)^{n}x^{2n+1}/\bigl(n!(2n+1)\bigr)$ , so $\sqrt\pi\operatorname{erf}(x) = \gamma(\tfrac12,x^{2}) = 2\sum_{n}(-1)^{n}x^{2n+1}/\bigl((2n+1)n!\bigr)$ , which is part 2's second display. The first few coefficients are $1$ , $-\tfrac13$ , $\tfrac1{10}$ , $-\tfrac1{42}$ .

Part 3 (exact truncation). Apply part 1 in the direction that lowers the first argument: replacing $a$ by $a-1$ there gives $\Gamma(a,x) = x^{a-1}e^{-x}+(a-1)\Gamma(a-1,x)$ . We claim by induction on $N$ that

\Gamma(a,x) \;=\; e^{-x}x^{\,a-1}\sum_{m=0}^{N-1}\frac{c_{m}}{x^{m}}+c_{N}\,\Gamma(a-N,x) .

For $N=1$ this is exactly the displayed relation, since $c_{0}=1$ and $c_{1}=a-1$ . Assume it for $N$ and apply the same relation to $\Gamma(a-N,x)$ , whose first argument is $a-N$ :

\Gamma(a-N,x) \;=\; x^{\,a-N-1}e^{-x}+(a-N-1)\,\Gamma(a-N-1,x) .

Substituting, and using $c_{N}(a-N-1) = c_{N+1}$ together with $c_{N}x^{a-N-1} = x^{a-1}\cdot c_{N}x^{-N}$ , gives the statement for $N+1$ . Note that no convergence question arises: each $\Gamma(a-N,x)$ is a finite integral because $x>0$ .

For the bound, let $b = a-N$ and suppose $N\ge a-1$ , that is $b-1\le0$ . Then $t\mapsto t^{\,b-1}$ is non-increasing on $[x,\infty)$ , so $t^{\,b-1}\le x^{\,b-1}$ there and

0 \;<\; \Gamma(b,x) \;=\; \int_{x}^{\infty}t^{\,b-1}e^{-t}dt \;\le\; x^{\,b-1}\int_{x}^{\infty}e^{-t}dt \;=\; x^{\,a-N-1}e^{-x} .

Multiplying by $\lvert c_{N}\rvert$ gives part 3's inequality, and the right-hand side is $\lvert c_{N}\rvert x^{-N}$ times the prefactor $e^{-x}x^{a-1}$ , i.e. the modulus of the term with $m=N$ — the first one omitted.

Part 4 (asymptotic character). Fix $N$ with $N\ge a-1$ . By part 3 the error committed by truncating after the term $m=N-1$ is at most $\lvert c_{N}\rvert e^{-x}x^{a-1-N}$ , which is $O(x^{-1})$ times the last retained term $c_{N-1}e^{-x}x^{a-N}$ . So each truncation has error of smaller order than its final term as $x\to\infty$ , which is the definition of an asymptotic expansion (Asymptotic series). Divergence for every fixed $x$ , when $a$ is not a positive integer: no factor $a-m$ then vanishes, so no $c_{m}$ does, and the ratio of consecutive terms has modulus $\lvert a-m-1\rvert/x$ , which exceeds $1$ once $m>x+\lvert a\rvert$ ; the terms therefore eventually increase without bound. For $a$ a positive integer the terms are $0$ from $m=a$ on and the sum is finite, which is part 5.

Part 5 (integer $a$ ). If $a=n$ is a positive integer then $c_{n} = (n-1)(n-2)\cdots(n-n) = 0$ , because the last factor is $0$ . Part 3 with $N=n$ therefore has a vanishing remainder:

\Gamma(n,x) \;=\; e^{-x}x^{\,n-1}\sum_{m=0}^{n-1}\frac{c_{m}}{x^{m}} , \qquad c_{m} \;=\; (n-1)(n-2)\cdots(n-m) \;=\; \frac{(n-1)!}{(n-1-m)!} .

Re-index with $s = n-1-m$ , so $x^{n-1-m} = x^{s}$ and $c_{m} = (n-1)!/s!$ :

\Gamma(n,x) \;=\; (n-1)!\,e^{-x}\sum_{s=0}^{n-1}\frac{x^{s}}{s!} ,

and subtracting from $\Gamma(n)=(n-1)!$ gives the companion formula for $\gamma(n,x)$ .

Part 6 (error functions). Take $a=\tfrac12$ . Then

c_{m} \;=\; \left(-\tfrac12\right)\left(-\tfrac32\right)\cdots\left(\tfrac12-m\right) \;=\; (-1)^{m}\,\frac{1\cdot3\cdots(2m-1)}{2^{m}} \;=\; (-1)^{m}\frac{(2m-1)!!}{2^{m}} ,

there being $m$ factors, each negative, with moduli $\tfrac12,\tfrac32,\dots,\tfrac{2m-1}{2}$ . Since $N\ge1>a-1=-\tfrac12$ for every $N\ge1$ , the remainder bound of part 3 applies for all $N$ . Writing part 3 at the point $x^{2}$ and multiplying by $\pi^{-1/2}$ , using $\operatorname{erfc}(x)=\pi^{-1/2}\Gamma(\tfrac12,x^{2})$ and $e^{-x^{2}}(x^{2})^{-1/2} = e^{-x^{2}}/x$ for $x>0$ , gives exactly part 6, with the remainder bound as stated. The signs alternate because the $c_{m}$ do and $x^{-2m}>0$ , and for an alternating sequence the bound "error at most the first omitted term" also makes consecutive truncations bracket the value. $\square$

Remark(Where the hypotheses are used)Original

$\operatorname{Re}a>0$ for $\gamma$ , nothing for $\Gamma$ . The lower integral is the one that reaches the origin, where $t^{a-1}$ blows up; the condition is exactly integrability there. It is used again in part 1, where the boundary term $t^{a}e^{-t}$ at $t=0$ is discarded, and this is why only the $\gamma$ recurrence carries the restriction while the $\Gamma$ recurrence holds for every real $a$ — including the negative values $a-N$ that the repeated integration by parts generates.
$x>0$ throughout parts 3-6. Every remainder $\Gamma(a-N,x)$ is a convergent integral only because the lower limit is positive; at $x=0$ the upper function with a negative first argument diverges. Correspondingly the large-argument expansion says nothing at small $x$ , and the prefactor $e^{-x^{2}}/(x\sqrt\pi)$ in part 6 is singular at $x=0$ , where the true value $\operatorname{erfc}(0)=1$ is perfectly finite.
$N\ge a-1$ for the remainder bound. This is what makes $t^{\,a-N-1}$ non-increasing on $[x,\infty)$ , so that its maximum is at the endpoint. Without it the bound is false in the stated form: for $a-N-1>0$ the factor grows and the integral is not controlled by its endpoint value. For the error functions $a=\tfrac12$ and every $N\ge1$ qualifies, which is why part 6 needs no side condition.
Boundedness of the interval in part 2. Term-by-term integration of the exponential series is justified by uniform convergence on a compact interval. Attempting the same on $[x,\infty)$ fails at once — $\sum(-1)^nt^n/n!$ is not uniformly convergent there — and this is the structural reason the upper function has no convergent expansion of that kind and must be handled by integration by parts instead.
Why the asymptotic series is still useful although it diverges. Divergence is not a defect to be worked around; it bounds the accuracy obtainable. The modulus of the term $m$ is $(2m-1)!!/(2^{m}x^{2m})$ for $\operatorname{erfc}$ , and the ratio of consecutive terms is $(2m+1)/(2x^{2})$ , which is below $1$ exactly while $m < x^{2}-\tfrac12$ . So the terms decrease until about $m\approx x^{2}$ and increase thereafter: truncating near there gives the best this expansion can do, and the bound in part 6 says what that is. The example below carries this out at $x=3$ .
Alternating signs are special to $a=\tfrac12$ . The bracketing property in part 6 comes from the signs of the $c_{m}$ , all of which are negative factors for $a=\tfrac12$ . For $a>1$ some factors $a-m$ are positive and the signs are not alternating, so only the remainder bound survives, not the bracketing.

Note(Three normalisations of the error function in this library)AdaptedAdapted from Hassani (2013), §14.6, p. 436, and §15.6, Problem 15.23, p. 487; Arfken, Weber & Harris (2013), Table 1.2, p. 66; Riley, Hobson & Bence (2006), §18.12.4, p. 640

The symbol $\operatorname{erf}$ is not used consistently across the cited books, and a formula moved between them without checking will be wrong by a factor or an additive constant. In this library:

Ours, and that of Arfken, Weber and Harris (Arfken, Weber & Harris (2013), Table 1.2, p. 66), Riley, Hobson and Bence (Riley, Hobson & Bence (2006), Eq. (18.167), p. 640), Strauss (Strauss (2008), Eq. (10), p. 51) and Stakgold and Holst (Stakgold & Holst (2011), §7.3, p. 449): $\operatorname{erf}(x) = \tfrac{2}{\sqrt\pi}\int_{0}^{x}e^{-u^{2}}du$ , with $\operatorname{erf}(0)=0$ and $\operatorname{erf}(\infty)=1$ .
Hassani (Hassani (2013), §14.6, p. 436) writes $\operatorname{erf}(x) = \tfrac{1}{\sqrt\pi}\int_{-\infty}^{x}e^{-t^{2}}dt$ , integrating from $-\infty$ and with half our prefactor. That function equals $\tfrac12\bigl(1+\operatorname{erf}(x)\bigr)$ in our notation, so it runs from $0$ to $1$ as $x$ runs over $\mathbb{R}$ and takes the value $\tfrac12$ at the origin. His initial data for $y''+2xy'=0$ , namely $y(0)=\tfrac12$ and $y'(0)=1/\sqrt{\pi}$ (p. 437), are those of his function; ours satisfies the same equation with $y(0)=0$ , $y'(0)=2/\sqrt{\pi}$ .
Hassani (Hassani (2013), §15.6, Problem 15.23, p. 487) writes, in the same book, $\operatorname{erf}(z) = \int_{0}^{z}e^{-t^{2}}dt$ with no prefactor at all, and in that normalisation records the confluent hypergeometric form $\int_{0}^{z}e^{-t^{2}}dt \;=\; z\,\Phi\bigl(\tfrac12,\tfrac32;-z^{2}\bigr) ,$ with $\Phi$ the confluent hypergeometric function. Nothing in this unit rests on that identification, so Confluent hypergeometric functions is a pointer onward rather than a prerequisite; all that is used of $\Phi$ here is the ratio of Pochhammer symbols recorded at the end of this bullet. Translated into our normalisation Hassani's relation reads $\operatorname{erf}(z) = \tfrac{2}{\sqrt\pi}z\,\Phi(\tfrac12,\tfrac32;-z^{2})$ ; it is the series of part 2 in hypergeometric dress, since $(\tfrac12)_{n}/(\tfrac32)_{n} = 1/(2n+1)$ .

So the three differ only by a constant factor and, in the second case, an additive constant. Our statements and numbers use the first throughout, and the $P,Q$ pair of Riley, Hobson and Bence relates to it by $P(\tfrac12,x^{2}) = \operatorname{erf}(x)$ and $Q(\tfrac12,x^{2}) = \operatorname{erfc}(x)$ .

Method(Evaluating

\operatorname{erf}

\operatorname{erfc}

to a target accuracy)Original

Decide which function you want. For $x$ large, $\operatorname{erf}(x)$ is close to $1$ and all the information is in $\operatorname{erfc}(x)$ ; computing it as $1-\operatorname{erf}(x)$ in finite precision destroys that information. Always compute the small one directly.
Choose the expansion by the size of $x$ . Use the convergent series of part 2 for $x$ up to roughly $2$ , and the asymptotic series of part 6 for larger $x$ . The two regions overlap, which is what makes the pair usable for all $x$ — but the overlap is one of usability, not of accuracy. The crossover is exactly where the asymptotic series is at its weakest: at $x=2$ its best truncation ( $N=4$ ) is still $1.4\%$ off, so in that region the convergent series is the one to trust, and step 4 will say so.
Convergent series: sum until the terms are small. The terms $t_{n} = 2(-1)^{n}x^{2n+1}/\bigl(\sqrt\pi(2n+1)n!\bigr)$ decrease in modulus as soon as $x^{2}<(2n+3)(n+1)/(2n+1)$ , and after that the series alternates, so the sum of the discarded tail is at most $\lvert t_{n+1}\rvert$ . Stop when that is below the tolerance.
Asymptotic series: check the bound before trusting the value. Add the terms $(-1)^{m}(2m-1)!!/(2^{m}x^{2m})$ while they decrease, that is while $m<x^{2}-\tfrac12$ , and use part 6 to bound the error by the first omitted term. If that bound is above the tolerance, the expansion cannot reach it at this $x$ and no number of terms will help: go back to step 3, or use a continued fraction or quadrature.
Sanity checks that cost nothing. $0<\operatorname{erfc}(x)<1$ for $x>0$ ; $\operatorname{erfc}$ is decreasing; and the first two terms of part 6 bracket the answer, so the true value must lie between $e^{-x^{2}}/(x\sqrt\pi)\cdot(1-\tfrac1{2x^{2}})$ and $e^{-x^{2}}/(x\sqrt\pi)$ for $x>1$ .
For general $\gamma(a,x)$ , $\Gamma(a,x)$ the same split applies, with the crossover near $x\approx a+1$ : series from part 2 below it, the expansion of part 4 above it, and the exact closed form of part 5 whenever $a$ is a positive integer.

Example(Both expansions, at

x=0.6

and at

x=3

)Original

Small argument: $\operatorname{erf}(0.6)$ from the convergent series. The terms $t_{n} = \tfrac{2}{\sqrt\pi}(-1)^{n}x^{2n+1}/\bigl((2n+1)n!\bigr)$ at $x=0.6$ are

$n$	$t_{n}$	partial sum
$0$	$+0.67702750$	$0.67702750$
$1$	$-0.08124330$	$0.59578420$
$2$	$+0.00877428$	$0.60455848$
$3$	$-0.00075208$	$0.60380640$
$4$	$+0.00005265$	$0.60385904$
$5$	$-0.00000310$	$0.60385594$
$6$	$+1.5745\times10^{-7}$	$0.60385610$
$7$	$-7.0178\times10^{-9}$	$0.60385609$

Eight terms give $\operatorname{erf}(0.6) = 0.603856091$ , and the next term bounds the error by $2.8\times10^{-10}$ ; the true value is $0.6038560908479$ , so the actual error is $-2.7\times10^{-10}$ , inside the bound as the alternation guarantees. The terms shrink fast because $x^{2}=0.36$ is small, and a tolerance of $10^{-9}$ costs eight terms.

Large argument: $\operatorname{erfc}(3)$ from the asymptotic series. Here the prefactor is

\frac{e^{-9}}{3\sqrt{\pi}} \;=\; \frac{1.2340980\times10^{-4}}{5.3173616} \;=\; 2.3208842\times10^{-5} ,

and the bracket terms $(-1)^{m}(2m-1)!!/(2^{m}\cdot 9^{m})$ are $1$ , $-\tfrac1{18}$ , $+\tfrac{3}{324}$ , $-\tfrac{15}{5832}$ , $+\tfrac{105}{104976}$ , i.e. $1$ , $-0.0555556$ , $+0.0092593$ , $-0.0025720$ , $+0.0010002$ :

$N$	truncated sum	$\operatorname{erfc}(3)$ estimate	error	bound from part 6
$1$	$1.0000000$	$2.3208842\times10^{-5}$	$+1.12\times10^{-6}$	$1.29\times10^{-6}$
$2$	$0.9444444$	$2.1919462\times10^{-5}$	$-1.71\times10^{-7}$	$2.15\times10^{-7}$
$3$	$0.9537037$	$2.2134359\times10^{-5}$	$+4.39\times10^{-8}$	$5.97\times10^{-8}$
$4$	$0.9511317$	$2.2074665\times10^{-5}$	$-1.58\times10^{-8}$	$2.32\times10^{-8}$

The true value is $\operatorname{erfc}(3) = 2.20904970\times10^{-5}$ . Every error is below its bound and alternates in sign, so consecutive rows bracket the answer: the value lies between $2.2074665\times10^{-5}$ and $2.2134359\times10^{-5}$ after four terms. That bracket is $2.7\times10^{-3}$ wide in relative terms and pins $\operatorname{erfc}(3)$ to $2.21\times10^{-5}$ : three significant figures from four terms of arithmetic, with the $N=4$ entry itself in error by $7.2\times10^{-4}$ relative.

Where it stops improving. The terms decrease while $m<x^{2}-\tfrac12 = 8.5$ , so the best truncation is $N=9$ . It gives $2.2092512\times10^{-5}$ , an absolute error of $2.0\times10^{-9}$ and a relative error of $9.1\times10^{-5}$ . Adding further terms makes matters worse: $N=12$ gives error $-2.7\times10^{-9}$ . About four correct digits is all this expansion will deliver at $x=3$ , for any number of terms.

The contrast with the convergent series at the same point. The convergent series does converge at $x=3$ , but it must first build up and then cancel away a value near $1$ : its largest term is $214.17$ (at $n=7$ ), and reaching an absolute accuracy of $10^{-9}$ in $\operatorname{erf}(3) = 0.9999779095$ needs $37$ terms. To extract $\operatorname{erfc}(3)\approx2.2\times10^{-5}$ from that one then subtracts from $1$ , and the cancellation costs about eight digits: $\sum_{n}\lvert t_{n}\rvert = 1629.99$ against an answer of $2.2090497\times10^{-5}$ is a ratio of $7.4\times10^{7}$ , and the final subtraction $1-\operatorname{erf}(3)$ accounts for $\log_{10}(1/\operatorname{erfc}(3)) = 4.7$ of those digits on its own. Four terms of the asymptotic series are both cheaper and better conditioned — this is the whole case for a divergent expansion.

Application(Diffusion from a step: the error-function profile)AdaptedAdapted from Strauss (2008), §2.4, Steps 1-3 and Eq. (5), pp. 48-49, and Example 1, p. 51; Stakgold & Holst (2011), §7.3, Eq. (7.3.22) and Duhamel's formula Eqs. (7.3.23)-(7.3.25), p. 449

The error function is not merely a convenient notation for a Gaussian integral; it is forced on the diffusion equation by a scaling argument, and the route is instructive.

Seek a solution of $u_{t}=ku_{xx}$ on the whole line depending only on the dimensionless combination $p = x/\sqrt{4kt}$ — the only such combination available, as Strauss (Strauss (2008), §2.4, Step 1, p. 48) argues from the invariance of the equation under $x\mapsto a^{1/2}x$ , $t\mapsto at$ . Writing $u = g(p)$ and using the chain rule, $u_{t} = -\tfrac{p}{2t}g'(p)$ and $u_{xx} = \tfrac{1}{4kt}g''(p)$ , so the equation collapses to the ordinary differential equation

g''(p)+2p\,g'(p) \;=\; 0 ,

which is the equation satisfied by the error function. Its integrating factor $e^{p^{2}}$ gives $g'(p) = C_{1}e^{-p^{2}}$ and hence $g(p) = C_{1}\int_{0}^{p}e^{-q^{2}}dq+C_{2}$ : every scale-invariant solution is an affine function of $\operatorname{erf}(p)$ . Matching the step initial data $u(x,0)=1$ for $x>0$ and $0$ for $x<0$ fixes the constants (Strauss, Step 3, p. 49) and yields

u(x,t) \;=\; \frac12+\frac12\operatorname{erf}\!\left(\frac{x}{\sqrt{4kt}}\right) .

A number. Join two long copper bars at $t=0$ , the half-line $x<0$ at $0^{\circ}$ C and the half-line $x>0$ at $100^{\circ}$ C, and ask for the temperature one centimetre inside the hot bar after $10$ s. Copper has thermal diffusivity $k = 1.11\times10^{-4}\,\mathrm{m^{2}\,s^{-1}}$ , so

\sqrt{4kt} \;=\; \sqrt{4\cdot1.11\times10^{-4}\cdot10} \;=\; 0.0666333\ \mathrm{m} , \qquad p \;=\; \frac{0.01}{0.0666333} \;=\; 0.1500751 .

The convergent series of part 2 gives $\operatorname{erf}(0.1500751) = 0.1680788$ in four terms, so the temperature is

T \;=\; 100\left(\frac12+\frac12\cdot0.1680788\right)\ ^{\circ}\mathrm{C} \;=\; 58.40\ ^{\circ}\mathrm{C} ,

and by the oddness of $\operatorname{erf}$ the mirror point one centimetre into the cold bar sits at $41.60\,^{\circ}$ C, the two averaging to the $50\,^{\circ}$ C that the join itself holds for all time. The diffusion length $\sqrt{4kt}=6.7$ cm after ten seconds is why the answer is so close to the mean of $50\,^{\circ}$ C: one centimetre is well inside the smeared front.

Boundary data instead of initial data. Stakgold and Holst (Stakgold & Holst (2011), §7.3, Eq. (7.3.22), p. 449) solve the companion half-line problem $u_{t}=u_{xx}$ on $x>0$ with $u(x,0)=0$ and $u(0,t)=1$ , and obtain the complementary error function: $u(x,t) = \operatorname{erfc}\bigl(x/\sqrt{4t}\bigr)$ , which is again of the form $C_{1}\int_{0}^{p}e^{-q^{2}}dq + C_{2}$ and so solves the same reduced equation, now with $g(0)=1$ and $g(\infty)=0$ . On this solution they build Duhamel's formula (their Eqs. (7.3.23)-(7.3.25)), which superposes time-shifted copies to solve the problem with an arbitrary boundary temperature $h(t)$ . The large-argument expansion of part 6 is what tells one how far the heat has not yet reached: at $p=3$ the surface temperature has raised the interior by only $2.2\times10^{-5}$ of its value.