Fact-checked by Grok 2 weeks ago

Polygamma function

The polygamma functions \psi^{(n)}(z), for positive integers n \geq 1, are meromorphic functions on the defined as the successive derivatives of the \psi(z) = \frac{\Gamma'(z)}{\Gamma(z)}, where \Gamma(z) is the . In particular, \psi^{(1)}(z) is known as the , \psi^{(2)}(z) as the tetragamma function, and higher-order cases follow analogously. These functions arise naturally in the analytic continuation and properties of the gamma function, providing essential tools for studying factorial generalizations and harmonic series in complex analysis. One key integral representation for the trigamma function is \psi'(z) = \int_0^\infty \frac{t e^{-zt}}{1 - e^{-t}} \, dt for \Re z > 0, with series expansions such as \psi^{(n)}(z) = (-1)^{n+1} n! \sum_{k=0}^\infty \frac{1}{(k+z)^{n+1}} extending their utility to broader domains via analytic continuation. The polygamma functions exhibit reflection formulas, like \psi^{(n)}(1 - z) = (-1)^{n+1} \psi^{(n)}(z) + (-1)^n \pi \frac{d^n}{dz^n} \cot(\pi z), which connect values across the complex plane and highlight their periodicity modulo transformations. Special values at integers link polygamma functions to the Riemann zeta function, for instance, \psi^{(n)}(1) = (-1)^{n+1} n! \zeta(n+1), underscoring their role in number theory and the evaluation of infinite series. Asymptotic expansions, such as \psi'(z) \sim \frac{1}{z} + \frac{1}{2z^2} + \sum_{k=1}^m \frac{B_{2k}}{z^{2k+1}} + O\left(\frac{1}{z^{2m+2}}\right) for large |z| in |\arg z| < \pi, facilitate approximations in large-argument regimes. Beyond pure mathematics, polygamma functions appear in inequalities for special functions and have applications in statistics (such as moments of logarithmic gamma distributions) and physics (including statistical and quantum mechanics).

Fundamentals

Definition

The polygamma function of order m at a complex number z is defined as the (m+1)-th derivative of the natural logarithm of the gamma function, \psi^{(m)}(z) = \frac{d^{m+1}}{dz^{m+1}} \ln \Gamma(z), where m is a nonnegative integer. This definition generalizes the class of functions arising from successive differentiation of the logarithmic derivative of the gamma function. For the special case m = 0, the polygamma function reduces to the digamma function, \psi(z) = \psi^{(0)}(z) = \frac{\Gamma'(z)}{\Gamma(z)} = \frac{d}{dz} \ln \Gamma(z). Higher-order polygamma functions, such as the trigamma function \psi^{(1)}(z) for m = 1, extend this by incorporating additional derivatives. The polygamma functions are meromorphic on the complex plane, exhibiting poles of order m+1 at the non-positive integers z = 0, -1, -2, \dots. The nomenclature "polygamma" originates from the Greek prefix "poly-," signifying many, to denote the multiple orders of derivatives extending beyond the digamma function in relation to the gamma function.

Relation to logarithm of gamma function

The polygamma function of order m \geq 1 arises as the mth derivative of the digamma function \psi(z) = \frac{d}{dz} \ln \Gamma(z), or equivalently, as the (m+1)th derivative of the logarithm of the gamma function: \psi^{(m)}(z) = \frac{d^{m+1}}{dz^{m+1}} \ln \Gamma(z). This relation establishes the polygamma function as the higher-order logarithmic derivative of the gamma function \Gamma(z). For \operatorname{Re}(z) > 0, an explicit series representation follows from the : \psi^{(m)}(z) = (-1)^{m+1} m! \sum_{k=0}^{\infty} \frac{1}{(z + k)^{m+1}}. This series converges absolutely in the right half-plane and provides an initial representation tied directly to the gamma function's properties. The connection to \ln \Gamma(z) enables of \psi^{(m)}(z) to the entire as a single-valued . The poles of \Gamma(z) at non-positive integers induce singularities in \ln \Gamma(z), which upon differentiation yield s of order m+1 for \psi^{(m)}(z) at z = 0, -1, -2, \dots. Near each such z = -k (k = 0, 1, 2, \dots), the principal part of the Laurent expansion begins with the term \frac{(-1)^{m+1} m!}{(z + k)^{m+1}}, with lower-order terms completing the expansion; this leading coefficient is independent of k. For m ≥ 1 and positive real arguments x > 0, the polygamma functions exhibit alternating monotonicity: \psi^{(m)}(x) is strictly decreasing when m is odd and strictly increasing when m is even. This behavior stems from the signs of the higher derivatives, with \psi^{(m)}(x) having sign (-1)^{m+1}. The polygamma function also connects to the through \psi^{(m)}(z) = (-1)^{m+1} m! \, \zeta(m+1, z), facilitating further and computation via the function's known properties in the .

Representations

Integral representation

The polygamma function of order m \geq 1 has a real integral representation \psi^{(m)}(z) = (-1)^{m+1} \int_0^\infty \frac{t^m e^{-z t}}{1 - e^{-t}} \, dt, valid for \Re(z) > 0. This form arises from the integral representation of the , \zeta(s, z) = \frac{1}{\Gamma(s)} \int_0^\infty \frac{t^{s-1} e^{-z t}}{1 - e^{-t}} \, dt, for \Re(s) > 1 and \Re(z) > 0, via the relation \psi^{(m)}(z) = (-1)^{m+1} m! \, \zeta(m+1, z). The integral converges for \Re(z) > 0 because the integrand behaves as t^{m-1} near t = 0 (integrable for m \geq 1) and decays exponentially as t \to \infty. This representation is well-suited for numerical evaluation through techniques, such as Gauss-Laguerre quadrature adapted for the non-standard weight $1/(1 - e^{-t}), providing efficient computation for moderate to large \Re(z) where the main contribution comes from small t. For small z, or subtraction of singularities may be needed to enhance stability. A contour integral representation follows from the Hankel contour form for the , \zeta(s, z) = \frac{1}{2\pi i} \int_H \frac{(-t)^{s-1} e^{z t}}{e^t - 1} \, dt, where H is Hankel's contour starting at +\infty just above the positive real axis, encircling the origin counterclockwise, and returning to +\infty just below the axis, valid for all complex s with appropriate branch cuts. Substituting s = m+1 and the relation to the polygamma function gives \psi^{(m)}(z) = (-1)^{m+1} \frac{m!}{2\pi i} \int_H \frac{(-t)^m e^{z t}}{e^t - 1} \, dt, for \Re(z) > 0. This contour form enables analytic continuation across the complex plane, avoiding branch points, and is valuable for deriving asymptotic expansions by deforming the contour. For specific orders, such as the trigamma function (m=1), the integral relates to Fermi-Dirac integrals via the identity $1/(1 - e^{-t}) = 1 + 1/(e^t - 1), linking to the form \int_0^\infty t e^{-z t} / (e^t + 1) \, dt = \Gamma(2) [\zeta(2, z) - (1/2) \zeta(2, z + 1/2)], though the primary representation uses the Bose-Einstein denominator. Even orders (m even) yield integrands with positive symmetry in certain Fourier or Laplace transforms, facilitating bounds and monotonicity proofs, while odd orders introduce antisymmetric components that align with Fermi-Dirac statistics in physical applications like quantum gases.

Series representation

The polygamma function of order m \geq 1 has an infinite series representation given by \psi^{(m)}(z) = (-1)^{m+1} m! \sum_{k=0}^{\infty} \frac{1}{(z + k)^{m+1}}, valid for \operatorname{Re}(z) > 0. This representation arises from repeated differentiation of the series for the and provides a direct method for numerical evaluation when the real part of z is positive, though convergence slows for small m and large |z|. This series connects to the through the relation \psi^{(m)}(z) = (-1)^{m+1} m! \zeta(m+1, z), where the is defined as \zeta(s, a) = \sum_{k=0}^{\infty} \frac{1}{(a + k)^s} for \operatorname{Re}(s) > 1 and \operatorname{Re}(a) > 0. This link facilitates of the polygamma function to the via the known properties of the Hurwitz zeta, which extends meromorphically beyond the convergence region of the series. For the trigamma function (m=1), an alternative involves numbers and even zeta values, though it is typically presented in asymptotic form for large |z|: \psi^{(1)}(z) \sim \frac{1}{z} + \frac{1}{2z^2} + \sum_{k=1}^{\infty} \frac{B_{2k}}{z^{2k+1}}, valid as |z| \to \infty with |\arg z| < \pi. Similar expansions hold for higher orders using generalized terms. The basic series converges slowly for low orders, prompting acceleration techniques such as the Euler-Maclaurin formula, which approximates the tail of the sum via integrals and number corrections, improving efficiency for computational purposes. For odd-order polygamma functions, the reflection formula enables reformulation as alternating series, enhancing convergence in regions where the standard sum is inefficient.

Functional Relations

Recurrence relation

The polygamma function of order m \geq 0 satisfies the recurrence relation \psi^{(m)}(z+1) = \psi^{(m)}(z) + (-1)^m \frac{m!}{z^{m+1}}, valid for z \neq 0, -1, -2, \dots . This relation is derived from the functional equation of the gamma function, \Gamma(z+1) = z \Gamma(z). Taking the natural logarithm yields \ln \Gamma(z+1) = \ln \Gamma(z) + \ln z, and differentiating once gives the digamma recurrence \psi(z+1) = \psi(z) + 1/z. Further differentiation m times applies the Leibniz rule to the term \ln z, resulting in the additional summand (-1)^m m! / z^{m+1}, while the polygamma terms on the left and right remain unchanged. The recurrence enables efficient numerical evaluation by reducing the argument to the interval (0,1]. For \operatorname{Re}(z) > 0, let n = \lfloor \operatorname{Re}(z) \rfloor and f = z - n \in [0,1); then \psi^{(m)}(z) is computed as \psi^{(m)}(f) plus the accumulated shifts, with \psi^{(m)}(f) obtained via series or integral representations. This approach avoids direct evaluation at large arguments where other methods may be less efficient. Iterating the recurrence for integer shifts n \geq 1 yields the generalization \psi^{(m)}(z + n) = \psi^{(m)}(z) + (-1)^m m! \sum_{k=0}^{n-1} \frac{1}{(z + k)^{m+1}}. This finite sum provides an exact relation for multiple steps, facilitating computation across the positive real axis. In numerical implementations, the choice between forward (increasing z) and backward (decreasing z) application of the recurrence affects stability, particularly for large |z|. Forward recurrence from a base value in (0,1] to larger arguments accumulates the sum terms directly, which is generally stable for moderate m and \operatorname{Re}(z) > 0 since the differences decrease with z. However, for very large |z| or higher orders m, backward recurrence—expressing values at smaller arguments in terms of larger ones—can help avoid overflow in intermediate power computations, akin to strategies for the itself.

Reflection formula

The reflection formula for the , which is the zeroth-order polygamma function \psi^{(0)}(z) = \psi(z), relates the values at z and $1 - z as follows: \psi(1 - z) - \psi(z) = \pi \cot(\pi z), valid for z \notin \mathbb{Z}. This identity is a direct consequence of the reflection formula for the , \Gamma(z) \Gamma(1 - z) = \frac{\pi}{\sin(\pi z)}, by taking the with respect to z. For higher-order polygamma functions \psi^{(m)}(z) with m \geq 1, the reflection formula generalizes to \psi^{(m)}(1 - z) = (-1)^m \psi^{(m)}(z) + (-1)^m \pi \frac{d^m}{dz^m} \cot(\pi z), again valid for z \notin \mathbb{Z}. This relation is obtained by differentiating the digamma reflection formula m times. For the first-order case (trigamma function, m=1), the formula simplifies to \psi^{(1)}(1 - z) + \psi^{(1)}(z) = \frac{\pi^2}{\sin^2(\pi z)}, which follows from the general form upon substitution and recognizing that \frac{d}{dz} \cot(\pi z) = -\pi \csc^2(\pi z). These reflection formulas are particularly useful for evaluating polygamma functions at arguments near the poles (non-positive integers) or in the negative real half-plane, where direct computation may be unstable, by relating them to values in the positive half-plane. They also facilitate computations for fractional arguments by mapping to the unit interval (0,1), often in combination with the to shift arguments while avoiding branch cuts of the underlying .

Multiplication theorem

The multiplication theorem for the polygamma function provides a relation between the value at a scaled argument kz and a sum of values at shifted arguments, where k is a positive . This theorem generalizes Gauss's multiplication formula for the through successive logarithmic . For the , which is the polygamma function of order zero denoted \psi(z) = \psi^{(0)}(z), the multiplication theorem states that \psi(kz) = \frac{1}{k} \sum_{j=0}^{k-1} \psi\left(z + \frac{j}{k}\right) + \ln k, valid for positive k and z such that kz and z + j/k avoid the non-positive integers where poles occur. This formula arises by taking the logarithm of Gauss's multiplication formula for the , \Gamma(kz) = (2\pi)^{(k-1)/2} k^{kz - 1/2} \prod_{j=0}^{k-1} \Gamma\left(z + \frac{j}{k}\right), which holds under the same conditions on z, and differentiating both sides with respect to z: k \psi(kz) = k \ln k + \sum_{j=0}^{k-1} \psi\left(z + \frac{j}{k}\right), yielding the stated relation. For polygamma functions of positive order n \geq 1, denoted \psi^{(n)}(z), the theorem simplifies by omitting the logarithmic term, as higher derivatives of the constant factor vanish: \psi^{(n)}(kz) = \frac{1}{k^{n+1}} \sum_{j=0}^{k-1} \psi^{(n)}\left(z + \frac{j}{k}\right), again for positive integer k and avoiding poles. This follows from further differentiation of the relation, where the \ln k term contributes nothing beyond the first . A representative example is the case k=2, known as the duplication formula. For the , \psi(2z) = \frac{1}{2} \left[ \psi(z) + \psi\left(z + \frac{1}{2}\right) \right] + \ln 2; for the (n=1), \psi^{(1)}(2z) = \frac{1}{4} \left[ \psi^{(1)}(z) + \psi^{(1)}\left(z + \frac{1}{2}\right) \right]. These relations facilitate computation of polygamma values at rational multiples of arguments, such as half-integers, often in conjunction with the reflection formula for auxiliary evaluation. The theorem applies only for positive [integer](/page/Integer) scaling factors $k$, as the underlying gamma multiplication formula is established for integers; extensions to non-integer $k$ do not hold in general and require alternative approaches. Convergence is ensured in the meromorphic sense within the [complex plane](/page/Complex_plane) excluding poles, with the sum converging uniformly on compact sets avoiding those singularities.[](https://dlmf.nist.gov/5.15) ## Series Expansions ### Taylor series expansion The polygamma function $\psi^{(m)}(z)$ is meromorphic in the [complex plane](/page/Complex_plane), with poles of order $m+1$ at the non-positive integers $z = 0, -1, -2, \dots$. Consequently, it admits a [Taylor series](/page/Taylor_series) expansion around any point $z_0$ avoiding these poles: \[ \psi^{(m)}(z_0 + h) = \sum_{n=0}^{\infty} \frac{\psi^{(m+n)}(z_0)}{n!} h^n. The radius of convergence of this series is the distance from z_0 to the nearest non-positive integer. This expansion follows from the general Taylor theorem for analytic functions, and can also be viewed through the generating function approach using the shift operator e^{h \frac{d}{dz}}, since \psi^{(m)}(z + h) = e^{h \frac{d}{dz}} \psi^{(m)}(z). Expanding the exponential yields the same power series in h, with coefficients involving higher-order polygamma functions evaluated at z_0. A explicit form for the coefficients arises when expanding around positive integers, where the values \psi^{(m+n)}(k) for positive integer k can be expressed using the Riemann zeta function via \psi^{(m+n)}(k) = (-1)^{m+n+1} (m+n)! \zeta(m+n+1, k), with the Hurwitz zeta \zeta(s, k) = \sum_{j=0}^{\infty} (k + j)^{-s}. This leads to coefficients that involve higher zeta values scaled by factorials. For expansions around non-integer points, coefficients generally require evaluation via integral representations or numerical methods, lacking closed forms in terms of zeta functions. For the specific case of expansion around z = 1 with m \geq 1 and |x| < 1, \psi^{(m)}(1 + x) = (-1)^{m+1} \sum_{j=0}^{\infty} (-1)^j \frac{(m + j)!}{j!} \zeta(m + j + 1) x^j. Equivalently, reindexing with k = j + 1, \psi^{(m)}(1 + x) = \sum_{k=1}^{\infty} (-1)^{m + k} \frac{(m + k - 1)!}{(k - 1)!} \zeta(m + k) x^{k-1}. The constant term matches the known value \psi^{(m)}(1) = (-1)^{m+1} m! \zeta(m+1). These coefficients are computed directly from higher-order zeta values \zeta(m + k), which for even arguments reduce to rational multiples of \pi^{m+k} via known formulas, while odd arguments involve more complex evaluations. For low orders like the trigamma function (m=1), the first few terms are \zeta(2) - 2 \zeta(3) x + 3 \zeta(4) x^2 - \cdots. In contrast to asymptotic expansions valid for large |h|, the Taylor series provides exact local approximations within the convergence disk.

Asymptotic expansion

The asymptotic expansion of the polygamma function \psi^{(m)}(z) for m \geq 1 and large |z| in the sector |\arg z| \leq \pi - \delta with \delta > 0 is given by \psi^{(m)}(z) \sim (-1)^{m-1} \left( \frac{(m-1)!}{z^m} + \frac{m!}{2z^{m+1}} + \sum_{k=1}^\infty \frac{(m + 2k - 1)!}{(2k)!} \frac{[B_{2k}](/page/Bernoulli_number)}{z^{m + 2k}} \right), where B_{2k} are the . This expansion provides a semi-convergent approximation, with terms involving even powers of $1/z after the leading contributions. An equivalent form expresses the sum using : \sum_{k=1}^\infty (-1)^{m+1} \frac{[B_{2k}](/page/Bernoulli_number)}{2k} \binom{m + 2k - 1}{2k - 1} \frac{1}{z^{m + 2k}}, preceded by the leading term (-1)^{m+1} (m-1)! / z^m and the next-order term (-1)^{m+1} m! / (2 z^{m+1}). For the digamma function (m = 0), the expansion takes a distinct form due to the logarithmic leading behavior: [\psi(z)](/page/Psi) \sim \ln z - \frac{1}{2z} - \sum_{k=1}^\infty \frac{B_{2k}}{2k z^{2k}}. This series arises directly from differentiating the series for \ln \Gamma(z), where \ln \Gamma(z) \sim \left(z - \frac{1}{2}\right) \ln z - z + \frac{1}{2} \ln (2\pi) + \sum_{k=1}^\infty \frac{B_{2k}}{2k(2k-1) z^{2k-1}}, since [\psi(z)](/page/Psi) = \frac{d}{dz} \ln \Gamma(z). Higher-order polygamma functions follow by further differentiation of this Stirling series. The expansions are valid in the complex plane excluding the branch cut along the negative real axis, corresponding to the principal sector |\arg z| < \pi, though practical use requires staying away from the boundary by a fixed angle \delta > 0 to ensure uniform convergence properties. As divergent asymptotic series, they do not converge for finite z, but partial sums up to an optimal number of terms provide the best approximation. The optimal truncation point occurs when the terms begin to increase in magnitude, typically after N \approx |z| terms, yielding an error on the order of the first omitted term.

Inequalities and Bounds

Trigamma function bounds

The trigamma function \psi^{(1)}(z) for z > 0 is strictly decreasing from +\infty as z \to 0^+ to 0 as z \to \infty, with the value \psi^{(1)}(1) = \pi^2/6. This monotonicity follows from the series representation \psi^{(1)}(z) = \sum_{n=0}^\infty \frac{1}{(z+n)^2} and the fact that its derivative, the tetragamma function \psi^{(2)}(z) = -2 \sum_{n=0}^\infty \frac{1}{(z+n)^3} < 0 for z > 0. For $0 < z < 1, the trigamma function satisfies the sharp bounds \frac{\pi^2}{6} < \psi^{(1)}(z) < \frac{\pi^2}{\sin^2(\pi z)}. The lower bound arises from the decreasing property and the value at z=1, while the upper bound is obtained from the reflection formula \psi^{(1)}(z) + \psi^{(1)}(1-z) = \pi^2 / \sin^2(\pi z) combined with the positivity of \psi^{(1)}(1-z) > 0. A fundamental lower bound for z > 0 is \psi^{(1)}(z) > 1/z, with equality approached only as z \to \infty. This inequality can be derived from the series using the integral test for the decreasing function f(x) = 1/(z+x)^2, yielding \sum_{n=0}^\infty f(n) > \int_0^\infty f(x) \, dx = 1/z. More refined lower and upper bounds near the asymptotic regime include, for x > 0, \frac{1}{x} - \frac{1}{2x^2} + \frac{1}{6x^3} - \frac{1}{30x^5} < \psi^{(1)}(x+1) < \frac{1}{x} - \frac{1}{2x^2} + \frac{1}{6x^3} - \frac{1}{30x^5} + \frac{1}{42x^7} - \frac{1}{30x^9} + \frac{5}{66x^{11}}, which provide error control for computational approximations. The asymptotic expansion for large z is \psi^{(1)}(z) \sim \frac{1}{z} + \frac{1}{2z^2} + \sum_{k=1}^\infty \frac{B_{2k}}{z^{2k+1}}, where B_{2k} are the Bernoulli numbers; this follows from the Euler-Maclaurin formula applied to the series representation. An exact integral representation underlying such expansions and bounds is \psi^{(1)}(z) = \int_0^\infty \frac{t e^{-z t}}{1 - e^{-t}} \, dt, \quad \operatorname{Re} z > 0. Double inequalities bounding the completely monotonic degree of the remainder after truncating the first few terms of the asymptotic expansion (e.g., after $1/z + 1/(2z^2)) have been established for improved error estimates in numerical evaluations.

General polygamma inequalities

The polygamma function of order m \geq 1 satisfies the sign inequality (-1)^{m+1} \psi^{(m)}(z) > 0 for all z > 0, implying that \psi^{(m)}(z) alternates in sign with increasing order: positive for odd m and negative for even m. Furthermore, |\psi^{(m)}(z)| is strictly decreasing in z for z > 0, with \psi^{(m)}(z+1) < \psi^{(m)}(z) when m is odd (both positive, approaching zero) and \psi^{(m)}(z+1) > \psi^{(m)}(z) when m is even (both negative, becoming less negative). These properties extend the base case of the trigamma function (m=1), where positivity and monotonic decrease hold without sign alternation. A fundamental bound for higher orders is |\psi^{(m)}(z)| \leq m! \, \zeta(m+1) for z \geq 1, derived from the series representation \psi^{(m)}(z) = (-1)^{m+1} m! \sum_{k=0}^\infty \frac{1}{(z+k)^{m+1}}, where the sum is maximized at z=1 and equals \zeta(m+1). Sharper bounds refine this, such as |\psi^{(m)}(z)| < \frac{m!}{z^{m+1}} + \frac{(m+1)!}{2z^{m+2}} + \sum_{i=1}^p \frac{B_{2i} \prod_{j=1}^{m} (2i + j - 1)}{z^{2i + m + 1}} for finite p, with alternating partial sums providing enclosures around the exact value; these are due to Alzer, who established S_m(2n; z) < |\psi^{(m)}(z)| < S_m(2n+1; z) using Bernoulli numbers B_{2i}, where S_m(p; z) denotes the partial asymptotic sum. Unlike the trigamma case, where bounds often involve simple constants like \pi^2/6, higher-order polygamma inequalities for m > 1 typically require zeta values or infinite series without closed-form simplification. Convexity properties hold for certain transformations of \psi^{(m)}(z). For instance, \psi^{(m)}(e^x) is on \mathbb{R} when m = 2n-1 (odd) and when m = 2n-2 (even, n \geq 2); more generally, \psi^{(m)}(x^c) exhibits or depending on the exponent c and of m, such as for odd m and c > 0. Log-convexity is less direct but follows for specific cases where the function and its logarithm are both , often in the of completely monotonic extensions. Recent refinements post-2000 focus on tight integral-based bounds, particularly for odd orders, with further developments as of 2024 including mean-value inequalities for polygamma and functions. and provided sharp double inequalities: for k \geq 1, x > 0, (k-1)! \left\{ x + \left[ (k-1)! \, |\psi^{(k)}(1)| \right]^{1/k} \right\}^k + \frac{k!}{x^{k+1}} < |\psi^{(k)}(x)| < (k-1)! \left( x + \frac{1}{2} \right)^k + \frac{k!}{x^{k+1}}, with the constants best possible; these improve upon bounds by incorporating the value at without full series expansion. For odd k, representations yield further refinements, emphasizing the loss of elementary constants compared to m=1.