The Riemann zeta function is one of the most celebrated functions in mathematics, with deep connections to prime numbers, number theory, and mathematical physics. A fascinating generalization arises when we replace its familiar single sum with nested sums, leading to multiple zeta values. These objects appear naturally in number theory, combinatorics, and quantum field theory, particularly in the evaluation of Feynman integrals. In Version 15 of the Wolfram Language, we introduced MultipleZeta, bringing these remarkable mathematical objects into the language as a built-in function. They are also now part of the answers that Wolfram Language can produce automatically from computations involving Sum, Integrate, and other symbolic operations.
The world of multiple zeta functions is much larger than one might initially expect. There are several families of multivariate zeta functions, each arising from different mathematical constructions and carrying its own interesting properties. In this article, we focus mainly on the Euler–Zagier multiple zeta function (MultipleZeta in WL V15), one of the most widely studied generalizations of the Riemann zeta function. We will explore its mathematical structure, some beautiful identities, and how to work with it in Wolfram Language. Along the way, we will also briefly touch upon a few other families of multiple zeta functions and see how they fit into this broader mathematical landscape.

Table and figure appearance


1. Riemann zeta

The natural starting point for our exploration of multiple zeta functions is the familiar Riemann zeta function. It is one of the most remarkable functions in mathematics. Its definition involves a simple infinite sum, yet it connects several seemingly unrelated subjects, including prime numbers, complex analysis, and mathematical physics.
Much of its fascination comes from the way questions about positive integers can be studied through the properties of a function of a complex variable. Its special values give exact answers to infinite sums, while its complex zeros reveal information about the distribution of prime numbers.
For
Re(s)>1
, the Riemann zeta function is defined as the sum of reciprocal powers of the positive integers:
ζ(s)=
∞
∑
n=1
1
s
n
Wolfram Language evaluates this sum directly and also provides the condition for convergence:
In[]:=
Sum
1
s
n
,{n,1,∞},GenerateConditionsTrue
Out[]=
Zeta[s] if Re[s]>1

From positive integers to prime numbers

One of the most beautiful properties of the zeta function is its connection with prime numbers. Euler discovered that the sum over all positive integers can also be written as an infinite product over primes [6]:
ζ(s)=
∏
pprime
1
1-
-s
p
Why should a sum over all integers equal a product involving only primes? Each factor can be expanded as a geometric series:
1
1-
-s
p
=1+
-s
p
+
-2s
p
+⋯
Multiplying these series together amounts to choosing a nonnegative power of every prime. Since each positive integer has a unique prime factorization, every integer appears exactly once. The product therefore reconstructs the original zeta series. Absolute convergence justifies this argument when
Re(s)>1
.
We can see this connection numerically. At
s2
, the zeta function has the famous value
ζ(2)
2
π
6
, the solution to the classical Basel problem. Taking a finite product over successively more primes gives increasingly accurate approximations - the last column records the difference from the exact value:
Out[]=
Number of primes
Finite Euler product
at s2
Error from ζ(2)
10
1.633070490
0.011863576
100
1.644515222
4.1884512×
-4
10
1000
1.644913175
2.0892142×
-5
10
Euler’s product is a starting point of analytic number theory, where techniques from complex analysis are used to investigate the distribution of primes.

Special values and integral representations

Another fascinating feature of the zeta function is its collection of special values. At positive even integers, the answers are rational multiples of powers of
π
. At positive odd integers, the situation is much more mysterious.
In[]:=
Table[Zeta[k],{k,2,6}]
Out[]=

2
π
6
,Zeta[3],
4
π
90
,Zeta[5],
6
π
945

The value
ζ(3)
, known as Apéry’s constant, is irrational, but comparably simple formulas are not known for the positive odd values in general.
The zeta function also appears naturally in integration. A classical representation, valid for
Re(s)>1
, is
∞
∫
0
s-1
x
x

-1
dx=Γ(s)ζ(s)
This follows by expanding
1
x

-1
into a geometric series in
-x

and integrating term by term. For instance, Wolfram Language evaluates the following integral in terms of
ζ(3)
:
In[]:=
Integrate
2
x
Exp[x]-1
,{x,0,∞}
Out[]=
2Zeta[3]
Thus the same constants arise from both infinite sums and integrals, a connection that becomes even richer for multiple zeta functions.

Analytic continuation and the Riemann hypothesis

The defining series converges only for
Re(s)>1
, but analytic continuation extends the zeta function to the entire complex plane, apart from a simple pole at
s1
. It also satisfies a functional equation relating its values at
s
and
1-s
.
The extended function has infinitely many complex zeros. The negative even integers
-2,-4,-6,…
are called its trivial zeros. All the other zeros, known as nontrivial zeros, lie in the critical strip
0<Re(s)<1
[6].
In 1859, Bernhard Riemann proposed one of the most famous unsolved problems in mathematics: every nontrivial zero lies on the critical line:
Re(s)=
1
2
Why is this important? The prime number theorem says that the number of primes up to
x
is approximately
x
log(x)
. The logarithmic integral
Li(x)
provides a more accurate smooth approximation. The Riemann hypothesis would give particularly strong bounds on how far the actual prime count can deviate from this approximation [2]. Thus a statement about complex zeros has direct consequences for the arithmetic of prime numbers.
Wolfram Language provides ZetaZero for computing zeros on the critical line. The first five such zeros in the upper half-plane are:
In[]:=
N[Table[ZetaZero[k],{k,1,5}]]
Out[]=
{0.5+14.1347,0.5+21.022,0.5+25.0109,0.5+30.4249,0.5+32.9351}
All five have real part
1
2
, just as the hypothesis predicts. Numerical verification, even for extremely many zeros, does not establish the statement for every nontrivial zero.

From one variable to several

The Riemann zeta function involves one summation index and one complex variable. What happens if we introduce multiple indices and multiple variables? This apparently simple extension opens up a much richer world of zeta functions.
There is no unique multivariate generalization. Different summation regions and denominators lead to distinct families, including Euler–Zagier, Barnes, Epstein, and Mordell–Tornheim zeta functions. They have different structures, identities, and applications.
An especially important example is the Euler–Zagier multiple zeta function. At depth two, instead of summing over a single integer, we sum over an ordered pair of positive integers:
ζ(
s
1
,
s
2
)=
∑
n
1
>
n
2
≥1
1
s
1
n
1
s
2
n
2
Here the summation indices satisfy
n
1
>
n
2
≥1
, and the formula is understood wherever the series converges. Such nested sums satisfy rich algebraic identities and arise in number theory as well as quantum field theory.
In Version 15, Wolfram Language introduced the built-in function MultipleZeta for Euler–Zagier multiple zeta values. We will now explore their definitions, convergence, special values, integral representations, and algebraic relations. Later, we will briefly return to some of the other multivariate zeta families.

2. First extension: MultipleZeta

The first extension of the Riemann zeta function that we will explore is obtained by replacing its single summation index with a sequence of strictly ordered indices. This leads to the Euler–Zagier multiple zeta function. Its definition looks like a modest change to an ordinary series, but it introduces a rich structure of nested sums, special values, and algebraic relations. Wolfram Language Version 15 represents this family with the built-in function MultipleZeta.

Series definition

Start with two indices. Instead of summing over every positive integer independently, require the inner index
m
to be smaller than the outer index
n
. This gives the depth-two multiple zeta function:
ζ(a,b)=
∑
n>m≥1
1
a
n
b
m
=
∞
∑
n=1
n-1
∑
m=1
1
a
n
b
m
The order matters: the first exponent belongs to the largest summation index. Some books use ascending indices and consequently reverse the arguments. Wolfram Language follows the descending-index convention displayed above. The corresponding symbolic Sum returns MultipleZeta together with the convergence conditions:
In[]:=
Sum
1
a
n
b
m
,{n,1,∞},{m,1,n-1},GenerateConditionsTrue
Out[]=
MultipleZeta[{a,b}] if Re[a]>1&&Re[a+b]>2
The same construction extends to any number
k
of indices:
ζ(
s
1
,
s
2
,…,
s
k
)=
∑
n
1
>
n
2
>…>
n
k
≥1
1
s
1
n
1
s
2
n
2
⋯
s
k
n
k
In Wolfram Language the general expression is MultipleZeta[{s1, s2, ..., sk}].
The nested summation defining
ζ(2,1,1)
can be visualized as a three-dimensional lattice of points satisfying
s
1
>
s
2
>
s
3
≥1
- each sphere represents a term in the sum, with its size and color indicating the magnitude of its contribution:
In[]:=
MultipleZetaRegionPlot[10]​​
Out[]=
When there is only one index, we recover the ordinary Riemann zeta function:
In[]:=
MultipleZeta[{s}]
Out[]=
Zeta[s]

Multiple zeta functions, values, depth, and weight

There is a useful distinction between a multiple zeta function and a multiple zeta value. Allowing the arguments
s
1
,…,
s
k
to vary over complex numbers gives a function of
k
variables, wherever the defining series converges. A classical multiple zeta value (MZV) is obtained by taking positive integer arguments with
s
1
≥2
. Those arguments are often called an admissible index list.
The number
k
of indices is the depth. For a classical MZV, the weight is their sum,
s
1
+…+
s
k
. For example,
ζ(2,1)
has depth 2 and weight 3, whereas
ζ(5,3)
has depth 2 and weight 8. Increasing the depth does not necessarily produce a genuinely new constant. Euler discovered the striking identity
ζ(2,1)ζ(3)
:
In[]:=
MultipleZeta[{2,1}]
The same MultipleZeta symbol handles both classical MZVs and general complex arguments. Here are two examples, one with rational and one with complex indices:

Why multiple zeta values/functions matter

This is the simplest example of a stuffle product. Other identities come from iterated integrals, giving shuffle products. Together they impose strong relations among MZVs and raise interesting questions about which constants are genuinely independent at each weight.
Multiple zeta values also arise in perturbative quantum field theory. In dimensional regularization, coefficients in the expansions of many multiloop Feynman integrals involve nested sums and iterated integrals that reduce to MZVs. Their appearance in Feynman diagrams, and connections with knot theory, were emphasized by Broadhurst and Kreimer [9]. The algebra of these values is therefore useful in both pure mathematics and symbolic physics computations; more general integrals may require larger classes of special functions.

MultipleZeta in the Wolfram Language V15

The built-in function can arise automatically from Sum, rather than being entered explicitly. For example, a depth-three nested sum evaluates to MultipleZeta[{8, 2, 4}]:
The Wolfram Function Repository resource PolyLogToSum converts a MultipleZeta expression into an inactive nested Sum, making the summation bounds easy to inspect:
Activating this expression evaluates the sum and returns the original multiple zeta function:
For numerical work, N evaluates MultipleZeta directly at a requested precision. We can compare it with NSum applied to the corresponding nested series. The agreement in this example provides a useful independent check:

Absolute convergence: which indices are allowed?

The following helper implements the full criterion, and a few examples show how it distinguishes convergent from divergent index lists:
Convergence criteria for MultipleZeta[{a, b, c}]:
Some examples of convergent and non-convergent multiple zetas:

Beyond convergence: analytic continuation

The convergence conditions describe the defining series, not every value of the underlying mathematical function. Like the ordinary Riemann zeta function, the Euler–Zagier multiple zeta function admits meromorphic continuation beyond the original region, although poles occur along certain hyperplanes [14]. A divergent series need not mean that the continued function has a pole at the same arguments.
The current MultipleZeta implementation does not provide a general-purpose analytic continuation for all divergent input series. For example, the following direct evaluations return Infinity; the first case is a genuine pole, whereas some others have finite values under meromorphic continuation:
These finite answers come from the continued function, not from summing a divergent defining series. Within the original convergence region, the helper agrees numerically with MultipleZeta:
A more systematic continuation subtracts additional asymptotic terms, for example through Euler–Maclaurin expansion, and expresses the removed pieces in terms of ordinary zeta functions [14]. We will otherwise work mainly in the convergence region, where the nested-sum interpretation is direct and unambiguous.

4. Multiple quasi zeta

ResourceFunction[“MultipleZetaStar”] implements the star function [3]:
A multiple zeta-star function can always be written as a rational linear combination of MultipleZeta e.g.:
The helper mzStarExpand implements the above expansion:
Some special values for different indices:

5. Integral representations: indices become words

The nested series defining a multiple zeta value has a second, rather different description: it is an iterated integral. Instead of summing over ordered integers, we integrate over ordered points of an interval. This is more than a change of notation. The integral picture explains many algebraic identities among multiple zeta values and provides the starting point for the shuffle relations discussed in the next section.

From indices to binary words

For convergent positive-integer multiple zeta values, the word begins with 0 and ends with 1. This is the admissibility condition that makes the corresponding endpoint integral convergent.
A couple of short Wolfram Language helpers implement the conversion in both directions:

Two differential forms and an ordered integral

The ordering matters: the integral is not a product of independent one-dimensional integrals. The word specifies which kernel belongs to each variable, and the nested bounds encode the same hierarchy that appears in the original multiple sum. Such objects are called Chen iterated integrals [1].
In Wolfram Language:
Wolfram Language evaluates this one-dimensional integral exactly:

Further examples

The Wolfram Function Repository provides PolyLogToIntegral [3], which can construct an inactive integral representation automatically:
Activating the resulting integral verifies the identity:
The connection also extends to harmonic polylogarithms, whose iterated-integral structure is closely related. At the endpoint 1, appropriate harmonic polylogarithms give multiple zeta values; for example:
The key insight is that a list of integer indices can be translated into a word, and that word into a sequence of differential forms. With the series and integral descriptions in hand, we can now multiply multiple zeta values in two different ways. The resulting stuffle and shuffle products are the subject of the next section.

6. Shuffle and stuffle: two ways to multiply

The preceding section gave multiple zeta values two complementary descriptions: nested sums and iterated integrals. Each description leads to a different rule for multiplying them. From sums we obtain the stuffle product; from integrals we obtain the shuffle product. Both calculate the same ordinary product, but they organize its terms differently. This is the source of the powerful double-shuffle identities.
We work with convergent multiple zeta values having positive integer indices and first index at least 2. The Wolfram Function Repository provides ShuffleProductExpand and StuffleProductExpand, which can act directly on MultipleZeta expressions.

Stuffle products: multiplying nested sums

The following plot makes this partition visible - the diagonal is what distinguishes stuffle multiplication from ordinary interleaving:
For longer index lists, the same idea recursively interlaces their entries, with an additional possibility of merging entries by addition. This is also called the quasi-shuffle or harmonic product. Here the repository function expands a product of two zeta values directly, without any auxiliary word function:

Shuffle products: multiplying iterated integrals

When two iterated integrals are multiplied, their integration domains can be divided into regions corresponding to all order-preserving interleavings of the two words. Unlike the stuffle rule, no letters merge. Different interleavings may nevertheless give the same word, producing integer coefficients.
Do the shuffle product expansion using the above mentioned resource function:

Double-shuffle relations: comparing the two answers

This gives the identity:
The following image summarizes the derivation:
Thus a multiple zeta value of depth two reduces to an ordinary zeta value. Activating the expressions checks the identity, and built-in evaluation gives the exact values of the two depth-two terms:
This gives another nontrivial identity:
The distinction is visible here: the stuffle product may merge adjacent indices and therefore produce terms of different depths, while the shuffle product preserves all the 1s in the binary words and thus has depth equal to the sum of the input depths. Both products preserve the total weight, because merging adds exponents in one description and interleaving preserves the total number of letters in the other.
These two algebraic viewpoints are fundamental tools for discovering and proving relations among multiple zeta values. They also provide a natural starting point for the exact-value experiments in the next section.

8. Colored multiple zeta values and MultiplePolyLog

The Euler–Zagier multiple zeta function sums products of reciprocal powers. We obtain a wider family by attaching a numerator factor to each summation index. Those factors can be positive, negative, or complex. Wolfram Language represents the resulting series with MultiplePolyLog. In the literature, colored MZVs often means root-of-unity weights. Here we use “colored” more broadly, allowing arbitrary complex numerator weights.
Other roots of unity produce periodic phases. More general rational or complex numerator weights can create exponentially damped or oscillatory sums:

9. Numerical evaluation

Partial sums and the Cauchy criterion

Certified evaluation of the tail

11. Other multivariate zeta functions

The Euler–Zagier multiple zeta function is one member of a much wider landscape. Different generalizations use lattices, linear or quadratic forms, combinatorial shapes, Lie-group representations, or enlarged nested-sum alphabets. The constructions below are important for different reasons and highlight the range of connections to analytic number theory, combinatorics, and mathematical physics.

Barnes zeta functions

For the two-dimensional example with equal coefficients, the definition becomes
Thus one Barnes specialization reduces to familiar zeta values. The general Barnes function, however, is a distinct family and is not the Euler–Zagier MultipleZeta.

Shintani zeta functions

Shintani-type zeta functions sum products of powers of several positive linear forms over lattice points in a cone. Separate exponents can be assigned to the forms. Such functions arise in the study of algebraic number fields and multivariate Dirichlet series [12].
The finite-sum identity can be checked exactly, without numerically approximating the infinite double sum:
The example illustrates how changing the summation region leads to a genuinely different multivariate zeta construction.

Epstein zeta functions

This special case is tied to a Dirichlet L-function rather than to nested zeta sums:

Schur multiple zeta functions

If every box has the same exponent s, symmetric-function identities turn this three-index sum into ordinary zeta functions:
The formula is not just a limiting coincidence. Truncating all three indices at r, we can check the corresponding identity exactly:
This family is closely connected to Schur functions, determinant formulas, and the combinatorics of Young tableaux. It is particularly relevant here because it extends the shuffle, stuffle, duality, and zeta-star themes already encountered [11].

Root-system zeta functions

The theory of zeta functions of root systems unifies several types of multiple zeta sums associated with Lie algebras and their Weyl-group symmetries. In particular, it provides a common framework for Euler–Zagier, Mordell–Tornheim, and Witten-type zeta functions [13].
The corresponding Wolfram Language calculation reduces the double sum to a single Euler sum:
Another important specialization is the Witten zeta function, defined by reciprocal powers of the dimensions of irreducible representations of a compact Lie group.
The simplest nontrivial even case has a clean closed form, which can also be written using Euler–Zagier multiple zeta values:
The example links the representation theory of SU(3) directly with MultipleZeta. Root-system zeta functions extend this phenomenon to much broader families of Lie groups.

Multiple affine and poly base zeta functions

Multiple affine zeta values. An affine letter can contain several shifted factors at one level:
Polynomial-base zeta values. Replace affine factors by arbitrary nonvanishing polynomial factors:
Their strict nested limits are
Thus ordinary MZVs are contained in the affine family, which is itself contained in the polynomial-base family. The distinction matters particularly for noninteger polynomial powers, which need not reduce to a finite collection of affine letters. Denominators must be nonzero on the summation range; complex powers require consistent branch choices.

A common thread

These are overlapping perspectives, not disjoint categories. Barnes and Shintani use linear forms on lattices, Epstein uses quadratic forms, Schur uses Young tableaux, and root-system zeta functions use coupled denominators associated with Lie theory. The affine and poly base families instead enlarge the permitted denominator letters in strictly nested sums; they contain ordinary MZVs as special cases.
The common computational theme is that reindexing, algebraic identities, and finite-sum stuffle relations can connect apparently different constructions. Wolfram Language V15 introduced MultipleZeta for Euler–Zagier values; Sum and numerical summation also let us explore these other families.

12. Concluding remarks

Multiple zeta functions form a rich and fascinating family of mathematical objects, connecting number theory, combinatorics, algebra, and mathematical physics. Their nested sums reveal remarkable identities and structures that extend far beyond the classical Riemann zeta function. With the introduction of MultipleZeta in Wolfram Language Version 15, we can explore these functions computationally in a systematic way. Symbolic simplification, numerical evaluation, identities involving shuffle and stuffle products, relations with harmonic sums and multiple polylogarithms, and connections with other special functions can now be investigated within a unified computational framework.
Multiple zeta functions are also part of a much broader landscape of multivariate zeta functions. Euler–Zagier, Mordell–Tornheim, Witten, Shintani and other generalizations arise from different summation structures, algebraic settings, and geometric constructions. Each brings its own mathematical properties, identities, and applications. Together, these functions illustrate how the simple idea of summing reciprocal powers can develop into an extensive mathematical theory. The new capabilities in Wolfram Language provide useful tools for exploring this landscape and open up further opportunities for symbolic and numerical computation.

References

[1] Wolfram Research. MultipleZeta, MultipleHarmonicNumber, and MultiplePolyLog: definitions, conventions, and examples. MultipleZeta MultipleHarmonicNumber MultiplePolyLog
[2] Clay Mathematics Institute. Riemann Hypothesis — statement and status. Open reference
[3] Jayanta Phadikar. Wolfram Function Repository tools for multiple zeta values, iterated integrals, and product identities. MultipleZetaStar PolyLogToIntegral StuffleProductExpand ShuffleProductExpand
[4] M. E. Hoffman and Y. Ohno. Relations of multiple zeta values and their algebraic expression. Journal of Algebra 262 (2003), 332–347. Open reference
[5] J. Blümlein, D. J. Broadhurst, and J. A. M. Vermaseren. The Multiple Zeta Value Data Mine. Computer Physics Communications 181 (2010), 582–625. Open reference Browse the Data Mine
[6] NIST Digital Library of Mathematical Functions. Riemann zeta: definition, analytic continuation, Euler product, and zeros. Definition and Euler product Zeros
[7] M. E. Hoffman. Multiple harmonic series. Pacific Journal of Mathematics 152 (1992), 275-290. Includes the Euler-Goldbach historical background. Open reference
[8] D. Zagier. Values of zeta functions and their applications. First European Congress of Mathematics, Volume II (1994), 497-512; especially Sections 8-9. Open reference
[9] D. J. Broadhurst and D. Kreimer. Association of multiple zeta values with positive knots via Feynman diagrams up to 9 loops. Physics Letters B 393 (1997), 403-412; preprint 1996. Open reference
[10] Jayanta Phadikar. Certified Arbitrary-Precision Evaluation of a Family of Generalized Multiple Zeta Functions. arXiv:2606.20091v2 (19 June 2026), especially Sections 3-4. Paper and ancillary files
[11] Maki Nakasuji, Ouamporn Phuksuwan, and Yoshinori Yamasaki. On Schur multiple zeta functions: A combinatoric generalization of multiple zeta functions. arXiv:1704.08511 (2017; revised 2018). Open reference
[12] Kohji Matsumoto. Analytic properties of multiple zeta-functions in several variables. Section 2 discusses sums of products of linear forms and Shintani-type constructions. Author's paper
[13] Yasushi Komori, Kohji Matsumoto and Hirofumi Tsumura. On Witten multiple zeta-functions associated with semisimple Lie algebras III. Introduction: representation dimensions, root-system zeta functions and the A2/Mordell-Tornheim correspondence. Author's paper
[14] Shigeki Akiyama, Shigeki Egami and Yoshio Tanigawa. Analytic continuation of multiple zeta-functions and their values at non-positive integers. Acta Arithmetica 98 (2001), 107-116. The paper uses increasing indices, so reverse the argument order to compare with this notebook. Author's paper Published article
[15] Jayanta Phadikar. Monoidal Alphabets for Generalized Harmonic Sums. arXiv:2605.21525 (2026). Open reference

CITE THIS NOTEBOOK

An introduction to multivariate zeta functions​
by Jayanta Phadikar​
Wolfram Community, STAFF PICKS, October 12, 2026
​https://community.wolfram.com/t/28317