Zeta Functions:
From Riemann to Feynman

Tigran Ishkhanyan, Wolfram Research, Inc.
Jayanta Phadikar, Wolfram Research, Inc.

#WolframTechConf

ABSTRACT

Zeta functions (such as the Riemann zeta and Multiple Zeta Values, MZVs), together with harmonic numbers and (multiple/generalized) polylogarithms are the native language of many difficult problems - from Euler sums in analytic number theory and nested sums in combinatorics to differential-equation methods for multi-scale Feynman integrals (iterated integrals and Fuchsian systems). We introduce a coherent family of  symbols - AlternatingHarmonicNumber, HyperHarmonicNumber, MultipleHarmonicNumber, MultipleZeta, HarmonicPolyLog, MultiplePolyLog, and GeneralizedPolyLog (GPL) - with exact evaluation and high-precision numerics, shuffle/stuffle-aware simplification, branch-aware transformations, and deep integration with Sum, RSolve, Integrate, and DSolve. These capabilities enable users to (i) collapse high-weight Euler sums to canonical MZV combinations, (ii) normalize mixed products of harmonic numbers/logs/polylogs for downstream Sum/Integrate, (iii) solve regular and perturbative Fuchsian systems directly in terms of GPLs, and (iv) verify results via one-line high-precision checks. The talk will open with an informal historical overview of zeta functions and then present a unified interface for Euler sums, harmonic number/polylog transforms, iterated integrals and Fuchsian systems.

Harmonic Numbers

Harmonic numbers are a bridge between discrete sums and analytic functions.

Classical definition: sum of the reciprocal of the first n natural numbers

(Subscript[H, n])=(∑_(k=1)^n)(1/k)

Sum[(1/k),{k,1,n}]

Table[HarmonicNumber[n],{n,1,10}]

Harmonic number (Subscript[H, n]) can be defined via integral:

(∫_0^1)((1-(x^n))/(1-x))dx

Integrate[((1-(x^n))/(1-x)),{x,0,1}]

Recurrence identity:

HarmonicNumber[z]==HarmonicNumber[z+1]-1/(z+1)//FullSimplify

The generating function for the harmonic numbers (Subscript[H, n]):

(∑_(n=1)^∞)(z^n)(Subscript[H, n])=((-ln(1-z))/(1-z))

GeneratingFunction[HarmonicNumber[n],n,z]

Harmonic numbers appear in several applications from number theory, special functions, etc.

Limit[HarmonicNumber[n]-Log[n],n->∞]

PolyGamma[n]==HarmonicNumber[n-1]-EulerGamma//FullSimplify

Generalized Harmonic Numbers

The n-th harmonic number (HarmonicNumber2[n, r]) of order r is given by:

(Subscript[H, n]^((r)))=(∑_(k=1)^n)(1/(k^r)) ,  (Subscript[H, n]^((1)))=(Subscript[H, n])

HarmonicNumber[4,r]

ComplexPlot3D[HarmonicNumber[n,1/2],{n,-1-I,1+I},PlotLegends->Automatic]

The n-th harmonic number (HarmonicNumber2[n, r])(s) of order r and decoration value s is given by:

(Subscript[H, n]^((r)))(s)=(∑_(k=1)^n)((s^k)/(k^r)) ,  (Subscript[H, n]^((r)))(1)=(Subscript[H, n]^((r)))

HarmonicNumber[4,r,s]

(∑_(k=1)^n)(s^k)(k^(-r))==HarmonicNumber[n,r,s]

ReImPlot[HarmonicNumber[3,0.3,s],{s,-4,4}]

Alternating Harmonic Numbers

Alternating harmonic numbers are also generalizations of the harmonic numbers:

(Subscript[H, n])=(∑_(k=1)^n)(1/k)

(Subscript[H, n]^((r)))=(∑_(k=1)^n)(1/(k^r))

(Subscript[(OverBar[H]), n])=(∑_(k=1)^n)((((-1))^(k+1))/k)

(Subscript[(OverBar[H]), n]^((r)))=(∑_(k=1)^n)((((-1))^(k+1))/(k^r))

Table[HarmonicNumber[n],{n,10}]
Table[AlternatingHarmonicNumber[n],{n,10}]

ReImPlot[AlternatingHarmonicNumber[t],{t,0,5}]

AlternatingHarmonicNumber[4,1/2]

(∑_(k=1)^n)(s^k)(k^(-r))==HarmonicNumber[n,r,s]

ReImPlot[HarmonicNumber[3,0.3,s],{s,-4,4}]

(AlternatingHarmonicNumber[z]==AlternatingHarmonicNumber[z+1]-((((-1))^z)/(Sqrt[(((1+z))^2)])))//FunctionExpand//FullSimplify

Hyper Harmonic Numbers

Hyperharmonic numbers (Subscript[H, n]^(<p>)) of level p are calculated recursively:

(Subscript[H, n]^(<p>))=(∑_(k=1)^n) (Subscript[H, k]^(<p-1>)), where  (Subscript[H, n]^(<0>))=1/n

For p=1 the hyperarmonic number is just the harmonic number:

(Subscript[H, n]^(<1>))=(Subscript[H, n])

Sum[HarmonicNumber[k],{k,n}]

Table[HarmonicNumber[n],{n,10}]
Table[HyperHarmonicNumber[2,n],{n,10}]

ReImPlot[HyperHarmonicNumber[p,n]/.p->2/3,{n,-∞,∞}]

By definition, the hyperharmonic numbers satisfy the recurrence relation (Subscript[H, n]^(<p>))= (Subscript[H, n-1]^(<p>))+ (Subscript[H, n]^(<p-1>))

HyperHarmonicNumber[p,n]==HyperHarmonicNumber[-1+p,n]+HyperHarmonicNumber[p,-1+n]//FullSimplify

Hyperharmonic numbers have generalizations:

(Subscript[H, n]^(<p>(r)))=(∑_(k=1)^n) (Subscript[H, k]^(<p-1>(r))),(Subscript[H, n]^(<1>(r)))=(Subscript[H, n]^((r))) gives the hyperharmonic number of level p with weight r

(Subscript[H, n]^(<p>(r)))(s)=(∑_(k=1)^n) (Subscript[H, k]^(<p-1>(r)))(s),(Subscript[H, n]^(<1>(r)))(s)=(Subscript[H, n]^((r)))(s) gives the hyperharmonic number of level p with weight r with decoration value s

Analytic Continuations

Harmonic numbers and their various variants were given above for positive integer n. To analytically continue these functions one need other special functions:

HarmonicNumber[n]//FunctionExpand
HarmonicNumber[n,r]//FunctionExpand
HarmonicNumber[n,r,s]//FunctionExpand

AlternatingHarmonicNumber[n]//FunctionExpand
AlternatingHarmonicNumber[n,r]//FunctionExpand

HyperHarmonicNumber[p,n]//FunctionExpand

Using the digamma, zeta, Hurwitz zeta, Lerch phi, polylog functions harmonic numbers can be further extended from initial discrete sums.

Also, these functions give the limits of harmonic numbers for n->∞:

HarmonicNumber[∞,r,s]
AlternatingHarmonicNumber[∞]
AlternatingHarmonicNumber[∞,r]

Multiple Harmonic Numbers

Another generalization of harmonic numbers is the multiple harmonic number MultipleHarmonicNumber[n,{(Subscript[r, 1]),(Subscript[r, 2]),…,(Subscript[r, k])}] of order n with weights (Subscript[r, i]):

(∑_(n >=(Subscript[n, 1])>...>(Subscript[n, k])>=1))((((Subscript[n, i])))^(-(Subscript[r, i])))

Table[MultipleHarmonicNumber[n,{2,1}],{n,1,10}]

DiscretePlot[MultipleHarmonicNumber[5,{11,4,t}],{t,0,15}]

MultipleHarmonicNumber[6,{4,3,5,2}]//N

NSum[(1/((n1^4) (n2^3) (n3^5)(n4^2) )),{n1,1,6},{n2,1,n1-1},{n3,1,n2-1},{n4,1,n3-1}]

MultipleHarmonicNumber[400,{3}]//N

HarmonicNumber[400,3]//N

For n->∞, multiple harmonic number is multiple zeta:

MultipleHarmonicNumber[400,{3,4,5}]//N

NSum[(1/((n1^3) (n2^4) (n3^5) )),{n1,1,400},{n2,1,n1-1},{n3,1,n2-1}]

MultipleZeta[{3,4,5}]//N

MultipleZeta

Finite Sum → HarmonicNumber[n, r]

(∑_(i=1)^n)(1/(i^r))

Infinite sum gives us Riemann Zeta (series definition of Zeta):

(∑_(i=1)^∞)(1/(i^r))

Zeta[4]==HarmonicNumber[∞,4]

Extend to multiple indices - we get MultipleZeta:

(∑_(i1=1)^∞)(∑_(i2=1)^(i1-1))(∑_(i3=1)^(i2-1))(∑_(i4=1)^(i3-1))(1/((i1^4) (i2^2) (i3^8) (i4^3)))

This is actually a MultipleHarmonicNumber at n = ∞ :

%==MultipleHarmonicNumber[∞,{4,2,8,3}]

For simple weights MultipleZeta will auto-evaluate:

MultipleZeta[{2,4,3}]

Numerical evaluation:

N[MultipleZeta[{4,2,8}]]

Integral definition of Zeta:

(∫_0^1)(∫_0^x1)(∫_0^x2)(∫_0^x3)(∫_0^x4)(1/(x1 x2 x3 x4 (1-x5)))\[DifferentialD]x5\[DifferentialD]x4\[DifferentialD]x3\[DifferentialD]x2\[DifferentialD]x1

Zeta[n]: n-1 x’s followed by (1-x) in Chen’s iterated integral. Upper limit of x1 is 1.

Insert more of (1-x)’s - we get MultipleZeta:

(∫_0^1)(∫_0^x1)(∫_0^x2)(∫_0^x3)(∫_0^x4)(∫_0^x5)(∫_0^x6)(∫_0^x7)(1/(x1 x2 x3 x4 (1-x5) x6 x7 (1-x8)))\[DifferentialD]x8\[DifferentialD]x7\[DifferentialD]x6\[DifferentialD]x5\[DifferentialD]x4\[DifferentialD]x3\[DifferentialD]x2\[DifferentialD]x1

Applications: Euler sum, Mordell-Tornheim sum, conical sum, Feynman integrals

MultiplePolyLog

Finite Sum → HarmonicNumber[n, r, s]

(∑_(i=1)^n)((s^i)/(i^r))

Infinite sum gives us PolyLog (series definition):

(∑_(i=1)^∞)((s^i)/(i^r))

HarmonicNumber[∞,2,1/3]==PolyLog[2,1/3]

Extend to multiple indices - we get MultiplePolyLog:

(∑_(i1=1)^∞)(∑_(i2=1)^(i1-1))(∑_(i3=1)^(i2-1))(∑_(i4=1)^(i3-1))((((((1/3)))^i1) ((((1/5)))^i2) ((((1/7)))^i3) ((((1/9)))^i4))/((i1^2) (i2^4) (i3^6) (i4^8)))

This is actually a decorated MultipleHarmonicNumber at n = ∞ :

%==MultipleHarmonicNumber[∞,{2,4,6,8},{1/3,1/5,1/7,1/9}]

For simple parameters MultiplePolyLog will auto-evaluate:

MultiplePolyLog[{2,1,2},{1,-1,1}]

Numerical evaluation:

N[MultiplePolyLog[{2,1,2},{(1/3),(1/4),(1/5)}]]

Applications: Euler sum, Feynman integrals

HarmonicPolyLog & GeneralizedPolyLog

Generalization of the series definition of PolyLog → MultiplePolyLog

Generalization of the integral definition of PolyLog → HarmonicPolyLog, GeneralizedPolyLog

Integral definition of PolyLog:

(∫_0^I)(∫_0^x1)(∫_0^x2)(∫_0^x3)(∫_0^x4)(1/(x1 x2 x3 x4 (x5-1)))\[DifferentialD]x5\[DifferentialD]x4\[DifferentialD]x3\[DifferentialD]x2\[DifferentialD]x1

PolyLog[n, k] → poles of the integrand are of the form: {0, 0, ..., 0, 1}; n-1 zeroes followed by 1

Integral definition of 3 argument PolyLog (Nielsen PolyLog):

(∫_0^I)(∫_0^x1)(∫_0^x2)(∫_0^x3)(∫_0^x4)(∫_0^x5)(1/(x1 x2 x3 x4 (x5-1) (x6-1)))\[DifferentialD]x6\[DifferentialD]x5\[DifferentialD]x4\[DifferentialD]x3\[DifferentialD]x2\[DifferentialD]x1

PolyLog[n1, n2, k] → poles are of the form: {0, 0, ..., 0, 1, 1, ..., 1}; n1 zeroes followed by n2 ones

Generalize to arbitrary poles from the set: {-1, 0, 1} → HarmonicPolyLog

(∫_0^I)(∫_0^x1)(∫_0^x2)(∫_0^x3)(∫_0^x4)(1/((x1+0) (x2+1) (x3-1) (x4+0) (x5-1)))\[DifferentialD]x5\[DifferentialD]x4\[DifferentialD]x3\[DifferentialD]x2\[DifferentialD]x1

%==HarmonicPolyLog[{0,-1,1,0,1},I]

Generalize to the set of arbitrary complex numbers → GeneralizedPolyLog (Goncharov PolyLog)

(∫_0^1)(∫_0^x1)(∫_0^x2)(∫_0^x3)(1/((x1-(2+6 I)) (x2-(3-I)) (x3-(1+2 I)) (x4-I)))\[DifferentialD]x4\[DifferentialD]x3\[DifferentialD]x2\[DifferentialD]x1

Hierarchy:

Log → PolyLog[n, z] → PolyLog[n, p, z] → HarmonicPolyLog → GeneralizedPolyLog

{GeneralizedPolyLog[{0},x],GeneralizedPolyLog[{0,0,1},x],GeneralizedPolyLog[{0,0,0,1,1},x],GeneralizedPolyLog[{-1,0,1,1},x]}

Reduction to simpler functions:

HarmonicPolyLog[{-1,1},1/2]

GeneralizedPolyLog[{I,2},1/3]

(* poles are drawn from {0, 1} and evaluated at endpoint 1 → MultipleZeta *)
GeneralizedPolyLog[{0,0,0,0,1,0,0,1},1]==HarmonicPolyLog[{0,0,0,0,1,0,0,1},1]==MultipleZeta[{5,3}]

Numerical evaluation:

N[HarmonicPolyLog[{-1,1,2},1/2]]

N[GeneralizedPolyLog[{I,2,3+I,1},1/2]]

Applications: Generalized Zeta, Feynman integrals (Iterated integrals, Fuchsian ODE system, Rational ODE system)

Applications

Series related applications

Solvers: Sum, RSolve, GeneratingFunction, ZTransform ...

Functions used: MultipleHarmonicNumber, MultipleZeta, MultiplePolyLog

Application area: Euler sums, Mordell-Tornheim sums. These are used in Number Theory and Feynman Integrals

Integral related applications

Solvers: Integrate, DSolve, AsymptoticDSolve ...

Functions used: MultipleZeta, HarmonicPolyLog, GeneralizedPolyLog

Application area: Iterated integrals, Fuchsian ODE system, Rational ODE system. These are used in Feynman Integrals.

Series Related Applications: Euler Sum

Definition of Euler Sum or “sum of harmonic numbers”:

(∑_(n=1)^k)(n^p) (m^n) (((Subscript[ℋ, 1])[n])^(Subscript[p, 1])) (((Subscript[ℋ, 2])[n])^(Subscript[p, 2]))...(((Subscript[ℋ, k])[n])^(Subscript[p, k]))

k: positive integer or ∞

p: integer

(Subscript[ℋ, 1]), (Subscript[ℋ, 2]), ... (Subscript[ℋ, k]) : harmonic numbers (HarmonicNumber, AlternatingHarmonicNumber, HyperHarmonicNumber, MultipleHarmonicNumber)

(Subscript[p, 1]), (Subscript[p, 2]), ... (Subscript[p, k]) : non-negative integers

Examples:

(∑_(n=1)^k)HarmonicNumber[n] AlternatingHarmonicNumber[n]

(∑_(n=1)^∞) (1/(n^3))HarmonicNumber[n]((HarmonicNumber[n,2])^3)

(∑_(n=1)^∞) (1/(n^2))((((1/2)))^n)HarmonicNumber[n]MultipleHarmonicNumber[n,{2,3},{1/4,1/5}]

Series Related Applications (contd..)

Certain solvers use Euler sum machinery in the background e.g. RSolve, GeneratingFunction, ExponentialGeneratingFunction, ZTransform and DirichletTransform

RSolveValue[{x[k]-2x[k-1]==HarmonicNumber[k,2]^2*k,x[1]==1},x[k],k]

GeneratingFunction[AlternatingHarmonicNumber[n,2]^2,n,x]

ExponentialGeneratingFunction[StirlingS1[n,3]HarmonicNumber[n],n,x]

ZTransform[MultipleHarmonicNumber[n,{2,3}],n,z]

DirichletTransform[HarmonicNumber[n,2],n,s]

Integral Related Applications: Iterated Integrals

In general, rational iterated integrals can be expressed as a sum of GeneralizedPolyLogs.

General form of rational iterated integral:

(∫_0^z)(∫_0^x1)...(∫_0^(Subscript[x, k-2]))(∫_0^(Subscript[x, k-1]))((Subscript[R, 1])[(Subscript[x, 1])](Subscript[R, 2])[(Subscript[x, 2])]...(Subscript[R, k-1])[(Subscript[x, k-1])](Subscript[R, k])[(Subscript[x, k])])\[DifferentialD](Subscript[x, k])\[DifferentialD](Subscript[x, k-1])...\[DifferentialD](Subscript[x, 2])\[DifferentialD](Subscript[x, 1])

(Subscript[R, 1])[(Subscript[x, 1])],(Subscript[R, 2])[(Subscript[x, 2])],...  are rational functions of  (Subscript[x, 1]),(Subscript[x, 2]),...

GeneralizedPolyLog may auto evaluate to HarmonicPolyLog, PolyLog, Log etc.

Examples:

(∫_0^1)(∫_0^x1)(∫_0^x2)(∫_0^x3)(∫_0^x4)(∫_0^x5)(∫_0^x6)(∫_0^x7)(1/(x1 x2 x3 x4 (1-x5) x6 x7 (1-x8)))\[DifferentialD]x8\[DifferentialD]x7\[DifferentialD]x6\[DifferentialD]x5\[DifferentialD]x4\[DifferentialD]x3\[DifferentialD]x2\[DifferentialD]x1

(∫_0^(I/2))(∫_0^x1)(∫_0^x2)(∫_0^x3)(∫_0^x4)(1/(x1 x2 x3 x4 (x5-1)))\[DifferentialD]x5\[DifferentialD]x4\[DifferentialD]x3\[DifferentialD]x2\[DifferentialD]x1

(∫_0^(I/2))(∫_0^x1)(∫_0^x2)(∫_0^x3)(1/((x1-1) x2 (x3-1) (x4+1)))\[DifferentialD]x4\[DifferentialD]x3\[DifferentialD]x2\[DifferentialD]x1

(∫_0^(I/2))(∫_0^x1)(∫_0^x2)(∫_0^x3)(1/((x1-(Sqrt[2])) (x2-(1/6)) (x3-(1+I)) (x4-(1-I))))\[DifferentialD]x4\[DifferentialD]x3\[DifferentialD]x2\[DifferentialD]x1

(∫_0^I)(∫_0^x1)(∫_0^x2)(((x1^2) x3)/(((x1^2)-1) ((x2^2)+1) (2 x3-6)))\[DifferentialD]x3\[DifferentialD]x2\[DifferentialD]x1

It has applications in multiple Hurwitz Zeta, multiple affine Zeta, solution of rational ODE systems etc.

Integral Related Applications: ODE system

Computation of iterated integrals lets us compute solution of linear Fuchsian/rational perturbative ODE system:

Y'[x]==ϵ M[x]Y[x]

Rational ODE system → elements of M[x] are rational functions of x

Fuchsian ODE system → elements of M[x] are rational functions of x with no polynomial part, and all poles are of multiplicity 1

(M={{1/(x-2)^2,-2/(x+1)},{1/(x+2),3/(x+I^3)}})//MatrixForm

AsymptoticDSolveValue[Join[Thread[{y1'[x],y2'[x]}==ϵ M.{y1[x],y2[x]}],{y1[0]==1,y2[0]==1}],{y1[x],y2[x]},x,{ϵ,0,3}]

In V14.3 such computations took more time and often came in terms of unevaluated integrals - as maximum order of ϵ increased.

Relationship with Feynman integrals:

Differential equation method for solving Feynman integrals - solution steps:

Step 1 → IBP reduction to set of master integrals (external package: LiteRed, FIRE ...)

Master integrals constitute a rational or a Fuchsian perturbative ODE system

Step 2 → reduction to a Fuchsian perturbative system in ϵ  form (Libra or Fuchsia)

Step 3 → perturbative solution (FeynCalc)

We will now natively support Step 3 (and to some extent Step 2) in Mathematica.

Concluding Remarks

New harmonic numbers:

HarmonicNumber[n] ( existing )
HarmonicNumber[n, r] ( existing )
HarmonicNumber[n, r, s] (new)

AlternatingHarmonicNumber[n]
AlternatingHarmonicNumber[n, r]
AlternatingHarmonicNumber[n, r, s]

HyperHarmonicNumber[p, n]
HyperHarmonicNumber[p, n, r]
HyperHarmonicNumber[p, n, r, s]

MultipleHarmonicNumber[n]
MultipleHarmonicNumber[n, {r1, r2, ..., rk}]
MultipleHarmonicNumber[n, {r1, r2, ..., rk}, {s1, s2, ..., sk}]

Applications: finite Euler sums, generating functions etc.

New Zetas and poly logs:

MultipleZeta
MultiplePolyLog
HarmonicPolyLog
GeneralizedPolyLog

Applications: Euler sums, iterated integrals, rational ODE system etc