The Riemann Zeta Function and Probability Distributions

Phi and Psi densities — Figure 1: Probability densities used to extend the zeta function

The famous Riemann zeta function was first introduced by Riemann in order to describe the distribution of the prime numbers. It is defined by the infinite sum

$\displaystyle \begin{aligned} \zeta(s) &=1+2^{-s}+3^{-s}+4^{-s}+\cdots\\ &=\sum_{n=1}^\infty n^{-s}, \end{aligned}$

(1)

which is absolutely convergent for all complex s with real part greater than one. One of the first properties of this is that, as shown by Riemann, it extends to an analytic function on the entire complex plane, other than a simple pole at ${s=1}$ . By the theory of analytic continuation this extension is necessarily unique, so the importance of the result lies in showing that an extension exists. One way of doing this is to find an alternative expression for the zeta function which is well defined everywhere. For example, it can be expressed as an absolutely convergent integral, as performed by Riemann himself in his original 1859 paper on the subject. This leads to an explicit expression for the zeta function, scaled by an analytic prefactor, as the integral of ${x^s}$ multiplied by a function of x over the range ${ x > 0}$ . In fact, this can be done in a way such that the function of x is a probability density function, and hence expresses the Riemann zeta function over the entire complex plane in terms of the generating function ${{\mathbb E}[X^s]}$ of a positive random variable X. The probability distributions involved here are not the standard ones taught to students of probability theory, so may be new to many people. Although these distributions are intimately related to the Riemann zeta function they also, intriguingly, turn up in seemingly unrelated contexts involving Brownian motion.

In this post, I derive two probability distributions related to the extension of the Riemann zeta function, and describe some of their properties. I also show how they can be constructed as the sum of a sequence of gamma distributed random variables. For motivation, some examples are given of where they show up in apparently unrelated areas of probability theory, although I do not give proofs of these statements here. For more information, see the 2001 paper Probability laws related to the Jacobi theta and Riemann zeta functions, and Brownian excursions by Biane, Pitman, and Yor.

The functional equation discovered by Riemann relates the values of ${\zeta(s)}$ to those of ${\zeta(1-s)}$ . Note that, as the series (1) is only defined on ${\Re(s) > 1}$ , there are no values of s for which it gives well-defined values of both ${\zeta(s)}$ and ${\zeta(1-s)}$ . So, extending the zeta function is a prerequisite. The functional equation is often expressed in terms of the xi function, defined as

$\displaystyle \xi(s)=\frac12s(s-1)\pi^{-s/2}\Gamma(s/2)\zeta(s).$

Here, ${\Gamma}$ is the gamma function, which is an analytic function defined by

$\displaystyle \Gamma(s)=\int_0^\infty x^{s-1}e^{-x}dx$

over ${\Re(s) > 0}$ . A simple application of integration by parts shows that it satisfies the functional equation ${\Gamma(s+1)=s\Gamma(s)}$ and, hence, extends to an analytic function on the complex plane, other than for simple poles at the non-positive integers. The xi function extends to an analytic function on the entire complex plane, and the Riemann functional equation is

$\displaystyle \xi(1-s)=\xi(s).$

It can be shown that ${\xi(s)}$ has no zeros in the half-plane ${\Re(s) > 1}$ , for example by using product expansions of the gamma and zeta functions which are valid on this range. So, the functional equation shows that all of its zeros lie in the critical strip ${0\le\Re(s)\le1}$ . In fact, there are infinitely many zeros in this range which, according to the long unsolved Riemann hypothesis, all lie on the line of symmetry ${\Re(s)=1/2}$ .

I will now outline a method of extending the zeta and xi functions and of deriving the functional equation, showing how it involves a probability distribution. We start with a simple rearrangement of a particular integral. Fixing a positive constant ${a}$ , we evaluate the following using the substitution ${y=ax^2}$ ,

$\displaystyle \begin{aligned} \int_0^\infty x^{s-1}e^{-ax^2}dx &=\int_0^\infty\left(a^{-1}y\right)^{(s-1)/2}e^{-y}\frac{dy}{2a^{1/2}y^{1/2}}\\ &=\frac12 a^{-s/2}\int_0^\infty y^{s/2-1}e^{-y}dy\\ &=\frac12 a^{-s/2}\Gamma(s/2) \end{aligned}$

(2)

Hence, the definition of the zeta function over ${\Re(s) > 1}$ gives,

$\displaystyle \begin{aligned} a^{-s/2}\Gamma(s/2)\zeta(s) &=\sum_{n=1}^\infty(an^2)^{-s/2}\Gamma(s/2)\\ &=2\int_0^\infty x^{s-1}\sum_{n=1}^\infty e^{-an^2x^2}dx. \end{aligned}$

Before taking this further, I mention the Jacobi identity, which applies to the summation inside the integral above. This is,

$\displaystyle \sum_ne^{-\pi n^2 x^2}=x^{-1}\sum_ne^{-\pi n^2x^{-2}}$

for all real ${x > 0}$ . For brevity, where I write a summation over n with no limits specified, it should be understood to be the sum over all integers, from ${-\infty}$ to ${\infty}$ . This is a special case of the following identity at ${y=0}$ ,

$\displaystyle \sum_ne^{-\pi (n+y)^2 x^2}=x^{-1}\sum_ne^{-\pi n^2x^{-2}+2\pi iny}.$

(3)

Note that the left hand side is periodic in y with period 1, and its Fourier series expansion can be computed, giving the right hand side. In particular, setting ${\psi(x)=\sum_ne^{-\pi n^2x^2}}$ , the Jacobi identity is ${\psi(1/x)=x\psi(x)}$ . Hence, we take ${a=\pi}$ above to obtain,

$\displaystyle \pi^{-s/2}\Gamma(s/2)\zeta(s)=\int_0^\infty x^{s-1}(\psi(x)-1)dx$

Now, multiply through by ${s(s-1)}$ and use the fact that ${s(s-1)x^{s-2}}$ is the second derivative of ${x^s}$ with respect to ${x}$ to obtain,

$\displaystyle \begin{aligned} 2\xi(s) &=s(s-1)\pi^{-s/2}\Gamma(s/2)\zeta(s)\\ &=\int_0^\infty\left(\frac{d^2}{dx^2}x^s\right)x(\psi(x)-1)dx\\ &=\int_0^\infty x^s\frac{d^2}{dx^2}\left(x\psi(x)\right)dx. \end{aligned}$

The third equality is just integration by parts, so we have obtained

$\displaystyle 2\xi(s)=\int_0^\infty x^s\Psi(x)dx$

(4)

where ${\Psi}$ is defined by

$\displaystyle \begin{aligned} \Psi(x) &=\frac{d^2}{dx^2}\left(x\psi(x)\right)\\ &=2\pi x\sum_nn^2\left(2\pi n^2x^2-3\right)e^{-\pi n^2x^2} \end{aligned}$

(5)

We note that this vanishes quickly as x goes to infinity.

Lemma 1 The following asymptotic limit holds as x tends to infinity,

$\displaystyle \Psi(x)\sim 8\pi^2x^3e^{-\pi x^2}.$ (6)

Proof: Definition (5) gives

$\displaystyle \frac{\Psi(x)}{8\pi^2x^3e^{-\pi x^2}} =\sum_{n=1}^\infty n^2(n^2-3/(2\pi x^2))e^{-\pi(n^2-1)x^2}$

and, by dominated convergence, we can take the limit ${x\rightarrow\infty}$ inside the summation. The ${n=1}$ term tends to 1 and the other terms tend to zero, giving the result. ⬜

In particular, ${\Psi(x)}$ vanishes faster than any power of x as x goes to infinity. Differentiating the Jacobi identity ${\psi(1/x)=x\psi(x)}$ twice gives the identity,

$\displaystyle \begin{aligned} \Psi(x)&=x^{-3}\Psi(1/x)\\ &=2\pi x^{-6}\sum_n n^2(2\pi n^2-3x^2)e^{-\pi n^2x^{-2}} \end{aligned}$

(7)

So, ${\Psi(x)}$ also vanishes faster than any power of x as x goes to zero. Hence, the integral (4) is well-defined and gives an analytic function of s over the complex plane. This proves the analytic extension of ${\xi}$ and, hence, of ${\zeta}$ . Using the substitution ${x=1/y}$ in (4) and applying identity (7) gives the functional equation,

$\displaystyle \begin{aligned} 2\xi(s) &=\int_0^\infty y^{-s}\Psi(1/y)y^{-2}dy\\ &=\int_0^\infty y^{1-s}\Psi(y)dy=2\xi(1-s). \end{aligned}$

Looking at the function ${\Psi}$ appearing in (4), we can show that it is positive and has integral equal to one, so is a probability density.

Theorem 2 The function ${\Psi}$ defines a probability density over the positive reals, and any random variable X with this density satisfies,

$\displaystyle {\mathbb E}[X^s]=2\xi(s)$ (8)

for all ${s\in{\mathbb C}}$ .

Proof: For all ${x\ge1}$ and nonzero integer n then,

$\displaystyle 2\pi n^2x^2-3\ge2\pi-3 > 0$

and, hence, (5) gives ${\Psi(x) > 0}$ . By the Jacobi identity (7), ${\Psi(x)=x^{-3}\Psi(1/x) > 0}$ for ${x\le1}$ , showing that ${\Psi}$ is positive everywhere. To show that it is a probability density, it must have integral equal to 1. Applying (4) and the functional equation,

$\displaystyle \begin{aligned} \int_0^\infty\Psi(x)dx &= 2\xi(0) = 2\xi(1)\\ &=\pi^{-1/2}\Gamma(1/2)\lim_{s\rightarrow1}(s-1)\zeta(s). \end{aligned}$

The identity ${\Gamma(1/2)=\sqrt\pi}$ is standard, as is the fact that ${(1-s)\zeta(s)\rightarrow1}$ as ${s\rightarrow1}$ . This shows that ${\Psi}$ has unit integral as required. Finally, (8) is just a different way of expressing (4). ⬜

Equation (8) describes the moments of a random variable X with probability density ${\Psi}$ . In particular, using the special values ${2\xi(1)=1}$ and ${\zeta(2)=\pi^2/6}$ ,

$\displaystyle \begin{aligned} &{\mathbb E}[X]=1,\\ &{\mathbb E}[X^2]=\frac\pi3. \end{aligned}$

As an example of an apparently unrelated situation where this distribution occurs, consider a standard Brownian bridge ${\{B_t\}_{t\in[0,1]}}$ . This is standard Brownian motion conditioned on hitting zero at time ${t=1}$ . Look at the difference between its maximum and minimum values.

$\displaystyle \sup_tB_t-\inf_tB_t=\sup_{0\le s\le t\le1}\lvert B_t-B_s\rvert$

Up to a constant scaling factor, this has probability density ${\Psi}$ . To be precise, if we scale it by ${\sqrt{2/\pi}}$ , then its probability density is ${\Psi}$ . As a little bit of trivia, suppose that now B is a standard Brownian motion. Then, ${B_t-tB_1}$ is a Brownian bridge over ${t\le1}$ . Applying the result just stated together with (8) for the moments gives

$\displaystyle \begin{aligned} &Z=\max_{0\le s\le t\le1}\lvert B_t-B_s+(s-t)B_1\rvert\\ &{\mathbb E}\left[Z^u\right]=2^{-\frac u2}u(u-1)\Gamma\left(\frac u2\right)\zeta(u). \end{aligned}$

At the time of writing, these are the equations on the notepad in the banner image of this site.

Another example is given by normalised Brownian excursions, which can be constructed as a Brownian bridge conditioned on being nonnegative. The maximum value of the excursion, again after scaling by ${\sqrt{2/\pi}}$ , has probability density ${\Psi}$ .

Riemann’s functional equation can alternatively be expressed as a symmetry of of the probability distribution introduced above.

Theorem 3 (The Functional Equation) Let X be a positive random variable with probability density ${\Psi}$ . Then,

$\displaystyle {\mathbb E}[Xf(1/X)]={\mathbb E}[f(X)]$ (9)

for all measurable functions ${f\colon{\mathbb R}_+\rightarrow{\mathbb R}_+}$ .

Proof: We can define a new measure ${{\mathbb Q}}$ on the underlying probability space by

$\displaystyle {\mathbb E}_{\mathbb Q}[f(X)]={\mathbb E}[Xf(1/X)]$

for all measurable functions ${f\colon{\mathbb R}_+\rightarrow{\mathbb R}_+}$ . This has generating function,

$\displaystyle {\mathbb E}_{\mathbb Q}[X^{s}]={\mathbb E}[X/X^{s}]=2\xi(1-s)=2\xi(s)={\mathbb E}[X^{s}]$

for real s. Taking ${s=0}$ gives ${{\mathbb E}_{\mathbb Q}[1]=1}$ , so that it is a probability measure. Then, as the distribution of a positive random variable is uniquely determined by its generating function, ${{\mathbb Q}={\mathbb P}}$ , giving (9). ⬜

Actually, identity (9) is equivalent to the functional equation. By linearity, it extends to any ${f\colon{\mathbb R}_+\rightarrow{\mathbb C}}$ such that ${f(X)}$ is integrable. Taking ${f(x)=x^s}$ gives,

$\displaystyle 2\xi(1-s)={\mathbb E}[X^{1-s}]={\mathbb E}[X/X^s]={\mathbb E}[X^s]=2\xi(s).$

Interestingly, we have seen the symmetry (9) previously in these notes for an entirely different distribution. Lemma 10 of the post on the normal distribution states that it holds for any lognormal random variable X of mean 1.

Another Distribution

Instead of starting with the usual series (1) defining the Riemann zeta function, we can instead consider a similar series which has alternating signs, giving the Dirichlet eta function

$\displaystyle \begin{aligned} \eta(s) &=1-2^{-s}+3^{-s}-4^{-s}+\cdots\\ &=\sum_{n=1}^\infty (-1)^{n+1}n^{-s}. \end{aligned}$

This has the benefit of converging for all ${\Re(s) > 0}$ which, importantly, includes the critical strip. The new function can be expressed in terms of the original zeta function,

$\displaystyle \begin{aligned} \eta(s) &= \sum_{n=1}^\infty n^{-s} -2\sum_{n > 1{\rm\ is\ even}}n^{-s}\\ &= \sum_{n=1}^\infty n^{-s} -2\sum_{n=1}^\infty (2n)^{-s}\\ &=(1-2^{1-s})\zeta(s). \end{aligned}$

In a similar fashion to the derivation above, we apply identity (2) with ${a=\pi n^2}$ ,

$\displaystyle \begin{aligned} \frac{1-2^{1-s}}{s-1}2\xi(s) &=s\pi^{-s/2}\Gamma(s/2)(1-2^{1-s})\zeta(s)\\ &=\sum_{n=1}^\infty s(-1)^{n+1}(\pi n^2)^{-s/2}\Gamma(s/2)\\ &=2\int_0^\infty sx^{s-1}\sum_{n=1}^\infty(-1)^{n+1}e^{-\pi n^2x^2}dx\\ &=\int_0^\infty\left(\frac{d}{dx}x^s\right)(1-\phi(x))dx \end{aligned}$

where we set

$\displaystyle \begin{aligned} \phi(x) &=\sum_n(-1)^ne^{-\pi n^2x^2}\\ &=2\psi(2x)-\psi(x). \end{aligned}$

(10)

Applying integration by parts, we have obtained

$\displaystyle \frac{1-2^{1-s}}{s-1}2\xi(s) =\int_0^\infty x^s\Phi(x)dx$

(11)

where I define

$\displaystyle \begin{aligned} \Phi(x)&=\phi^\prime(x)=4\psi^\prime(2x)-\psi^\prime(x)\\ &=2\pi x\sum_n(-1)^{n+1}n^2e^{-\pi n^2x^2}. \end{aligned}$

A Jacobi style identity is also satisfied by ${\phi}$ , as can be seen either by taking ${y=1/2}$ in (3), or by using (10) to express in terms ${\psi}$ and applying the Jacobi identity for ${\psi}$ ,

$\displaystyle \begin{aligned} \phi(x) &=x^{-1}\sum_ne^{-\pi (n+1/2)^2x^{-2}}\\ &=x^{-1}\left(\psi(1/(2x))-\psi(1/x)\right). \end{aligned}$

(12)

Hence, ${\Phi}$ satisfies the identity,

$\displaystyle \Phi(x)=x^{-4}\sum_n(2\pi(n+1/2)^2-x^2)e^{-\pi(n+1/2)^2 x^{-2}}$

(13)

These expressions show that ${\Phi(x)}$ is asymptotic to ${4\pi xe^{-\pi x^2}}$ as x goes to infinity, and ${\pi x^{-4}e^{-\pi x^{-2}/4}}$ as x goes to zero. In either case, ${\Phi(x)}$ vanishes faster than any power of x so the integral in (11) is defined for all ${s\in{\mathbb C}}$ , giving an alternative extension of ${\xi(s)}$ to the complex plane. Also, using the expressions above for ${\Phi}$ and ${\Psi}$ in terms of ${\psi}$ gives the identity,

$\displaystyle x\Phi^\prime(x)+2\Phi(x)=4\Psi(2x)-\Psi(x).$

This is a differential equation for ${\Phi}$ in terms of ${\Psi}$ , which can be integrated to give

$\displaystyle x^2\Phi(x)=\int_x^{2x}y\Psi(y)dy.$

(14)

This leads us to another probability density.

Theorem 4 The function ${\Phi}$ defines a probability density over the positive reals, and any random variable X with this density satisfies,

$\displaystyle \begin{aligned} E[X^s] &=s(1-2^{1-s})\pi^{-s/2}\Gamma(s/2)\zeta(s)\\ &= \frac{1-2^{1-s}}{s-1}2\xi(s) \end{aligned}$ (15)

for all ${s\in{\mathbb C}}$ .

Proof: Equation (14) together with positivity of ${\Psi}$ shows that ${\Phi}$ is positive. To be a probability density, it must integrate to 1. Applying (11),

$\displaystyle \int_0^\infty\Phi(x)dx=\frac{1-2^{1-0}}{0-1}2\xi(0)=2\xi(1)=1$

as required. Finally, (15) follows from (11). ⬜

As an example of where this second probability distribution appears in a context apparently unrelated to the Riemann zeta function, consider an IID sequence of random variables ${X_1,X_2,\ldots}$ with continuous distribution function ${F(x)={\mathbb P}(X_n < x)}$ . We also consider the sample approximations to this given by the proportion of the first n values of the ${X_k}$ which are less than x,

$\displaystyle F_n(x)=\frac1n\sum_{k=1}^n1_{\{X_k < x\}}.$

The maximum discrepancy between this sample distribution function and the true underlying one is,

$\displaystyle \sup_{x\in{\mathbb R}}\lvert F_n(x)-F(x)\rvert.$

This is a random quantity which, after scaling by the factor ${\sqrt{2n/\pi}}$ , converges in distribution as ${n\rightarrow\infty}$ to the probability measure with density ${\Phi}$ .

Another example is given by a standard Brownian bridge B. The scaled absolute maximum ${\sqrt{2/\pi}\sup_t\lvert B_t\rvert}$ can be shown to have probability density ${\Phi}$ . Additionally, the ‘sample standard deviation’

$\displaystyle \left(\int_0^1\left(B_t-\int_0^1B_sds\right)^2dt\right)^{1/2}$

after scaling by a factor ${\sqrt{2\pi}}$ , also has probability density ${\Phi}$ . A proof of this is given below, in lemma 12.

A further example is provided by a Brownian meander ${\{B_t\}_{t\in[0,1]}}$ . This is standard Brownian motion conditioned on being nonnegative over the interval ${[0,1]}$ , and it is known that ${\sup_tB_t/\sqrt{2\pi}}$ has probability density ${\Phi}$ .

Properties of the Distributions

I will now describe and prove some of the properties of the distributions, including computing their cumulative distribution and moment generating functions. I also show how they can be realized as a sum of gamma distributed random variables. First, the two probability densities are related in the following simple but surprising way.

Lemma 5 Let X have probability density ${\Psi}$ and, independently, let Y be uniformly distributed on the interval ${[1,2]}$ . Then, ${X/Y}$ has probability density ${\Phi}$ .

Proof: The generating function of ${Z=X/Y}$ can be computed,

$\displaystyle \begin{aligned} {\mathbb E}[Z^s] &={\mathbb E}[Y^{-s}]{\mathbb E}[X^s]\\ &=\int_1^2y^{-s}dy\,2\xi(s)\\ &=\frac{1-2^{1-s}}{s-1}2\xi(s) \end{aligned}$

which, according to (15), is the generating function corresponding to a variable with probability density ${\Phi}$ . ⬜

Less surprisingly, the cumulative distribution functions can be computed directly from the definitions above.

Lemma 6 If X has probability density ${\Phi}$ , then its distribution function is,

$\displaystyle \begin{aligned} {\mathbb P}(X < x) &=\phi(x)\\ &=\sum_n(-1)^ne^{-\pi n^2x^2}\\ &=x^{-1}\sum_ne^{-\pi (n+1/2)^2x^{-2}}. \end{aligned}$ (16)

If it has probability density ${\Psi}$ , then its distribution function is,

$\displaystyle \begin{aligned} {\mathbb P}(X < x) &=\frac{d}{dx}\left(x\psi(x)\right)\\ &=\sum_n(1-2\pi n^2x^2)e^{-\pi n^2x^2}\\ &=2\pi x^{-3}\sum_n n^2e^{-\pi n^2x^{-2}} \end{aligned}$ (17)

Proof: For the first equation, if X has probability density ${\Phi=\phi^\prime}$ then, as ${\phi}$ vanishes at zero, the cumulative distribution is given by

$\displaystyle \int_0^x\phi^\prime(y)dy=\phi(x).$

For equation (17), equality of each of the three expressions on the right hand side is given by differentiating the defining series of ${\psi(x)}$ and of ${x\psi(x)=\psi(1/x)}$ term by term. In particular this shows that ${(d/dx)x\psi(x)}$ vanishes at zero, and the distribution function is given by integrating,

$\displaystyle \int_0^x\frac{d^2}{dy^2}\left(y\psi(y)\right)dy =\frac d{dx}\left(x\phi(x)\right).$

⬜

The calculation of the moment generating functions will require the following integral identity due to Lévy.

Lemma 7 For any positive reals ${a,b}$ , the integral identity

$\displaystyle \int_0^\infty ae^{-a^2x^2-b^2x^{-2}}dx=\frac{\sqrt\pi}2e^{-2ab}$ (18)

holds.

Proof: Using the substitution ${x=b/(ay)}$ ,

$\displaystyle ae^{-a^2x^2-b^2x^{-2}}dx=-by^{-2}e^{-b^2y^{-2}-a^2y^2}dy.$

The integral in the lemma can be split into the ranges ${x >\sqrt{b/a}}$ and ${x < \sqrt{b/a}}$ . Applying the substitution above over ${x < \sqrt{b/a}}$ gives,

$\displaystyle \int\limits_{\sqrt{b/a}}^\infty(a+bx^{-2}) e^{-a^2x^2-b^2x^{-2}}dx.$

Next, apply the substitution ${z=ax-b/x}$ so that,

$\displaystyle dz=(a+b/x^2)dx,$

which transforms the integral to,

$\displaystyle \int_0^\infty e^{-z^2-2ab}dz=\frac{\sqrt\pi}{2}e^{-2ab}$

as required. ⬜

We now compute the moment generating functions for squares of random variables with the distributions above. The results are surprisingly simple, and it is not clear why this should be so.

Lemma 8 If X has probability density ${\Phi}$ then,

$\displaystyle {\mathbb E}\left[e^{-\pi^{-1}\lambda^2X^2}\right]=\frac{\lambda}{\sinh\lambda}$

for all real ${\lambda\not=0}$ . If X has probability density ${\Psi}$ then,

$\displaystyle {\mathbb E}\left[e^{-\pi^{-1}\lambda^2X^2}\right]=\left(\frac{\lambda}{\sinh\lambda}\right)^2.$

Proof: As the exponential term inside the expectations tends to 1 as X goes to zero, rather than vanishing, we use the form for the probability densities which explicitly vanish at zero, to ensure integrability of all the terms. If X has probability density ${\phi^\prime(x)}$ , we use expression (12) for ${\phi}$ , and integration by parts,

$\displaystyle \begin{aligned} {\mathbb E}\left[e^{\pi^{-1}\lambda^2X^2}\right] &=\int_0^\infty \phi^\prime(x)e^{-\pi^{-1}\lambda^2 x^2}dx\\ &=2\lambda\int_0^\infty \pi^{-1}\lambda x\phi(x)e^{-\pi^{-1}\lambda^2 x^2}dx\\ &=4\lambda\sum_{n=0}^\infty\int_0^\infty\pi^{-1}\lambda e^{-\pi^{-1}\lambda^2 x^2-\pi(n+1/2)^2x^{-2}}dx\\ &=2\lambda\sum_{n=0}^\infty e^{-\lambda(2n+1)} =2\lambda e^{-\lambda}/(1-e^{-2\lambda})\\ &=\lambda/\sinh\lambda \end{aligned}$

as required. Here, (18) was used to perform the integral.

Next, suppose that X has probability density ${\Psi(x)=(x\psi(x))^{\prime\prime}}$ . We use the form,

$\displaystyle x\psi(x)=\psi(1/x)=1+2\sum_{n=1}^\infty e^{-n^2\pi x^{-2}}.$

Using integration by parts, and the notation ${\partial_\lambda}$ for the derivative with respect to ${\lambda}$ ,

$\displaystyle \begin{aligned} {\mathbb E}\left[e^{\pi^{-1}\lambda^2X^2}\right] &=\int_0^\infty (x\psi(x))^{\prime\prime}e^{-\pi^{-1}\lambda^2 x^2}dx\\ &=\int_0^\infty(4\pi^{-2}\lambda^4x^2-2\pi^{-1}\lambda^2) (x\psi(x)-1) e^{-\pi^{-1}\lambda^2 x^2}dx\\ &=-2\lambda^2\partial_\lambda\int_0^\infty\pi^{-1}\lambda (x\psi(x)-1) e^{-\pi^{-1}\lambda^2 x^2}dx\\ &=-4\lambda^2\partial_\lambda\sum_{n=1}^\infty\int_0^\infty\pi^{-1}\lambda e^{-\pi^{-1}\lambda^2 x^2-\pi n^2 x^{-2}}dx\\ &=-2\lambda^2\partial_\lambda\sum_{n=1}^\infty e^{-2\lambda n} =-2\lambda^2\partial_\lambda(1/(e^{2\lambda}-1))\\ &=4\lambda^2e^{2\lambda}/(e^{2\lambda}-1)^2 =(\lambda/\sinh\lambda)^2 \end{aligned}$

as required. ⬜

Knowing the moment generating functions opens up a world of possibilities. For example, we find the following simple relation between the two distributions.

Lemma 9 Let X and Y be independent random variables with probability density ${\Phi}$ . Then, ${\sqrt{X^2+Y^2}}$ has probability density ${\Psi}$ .

Proof: The moment generating function of ${Z^2=X^2+Y^2}$ can be computed using lemma 8.

$\displaystyle \begin{aligned} {\mathbb E}[e^{-\pi^{-1}\lambda^2Z^2}] &= {\mathbb E}[e^{-\pi^{-1}\lambda^2X^2}] {\mathbb E}[e^{-\pi^{-1}\lambda^2Y^2}]\\ &=(\lambda/\sinh\lambda)^2. \end{aligned}$

Applying lemma 8 again and using the fact that the distribution of a nonnegative random variable is uniquely determined by its moment generating function shows that Z has probability density ${\Psi}$ . ⬜

We now have lemma 5 which give us a method of converting a random variable with density ${\Psi}$ to one with density ${\Phi}$ . In the other direction, lemma 9 goes from a pair of independent variables with density ${\Phi}$ to one with density ${\Psi}$ . Combining these gives a characterization of these distributions.

Corollary 10 Let X,Y,U,V be independent random variables with U,V uniform on ${[1,2]}$ .

If X,Y have probability density ${\Phi}$ , then so does ${U^{-1}\sqrt{X^2+Y^2}}$ .

If X,Y have probability density ${\Psi}$ , then so does ${\sqrt{U^{-2}X^2+V^{-2}Y^2}}$ .

Proof: For the first statement, ${\sqrt{X^2+Y^2}}$ has density ${\Psi}$ by lemma 9, so the result follows from lemma 5.

For the second statement, ${U^{-1}X}$ and ${V^{-1}Y}$ are independent and have density ${\Phi}$ by lemma 5, so the result follows from lemma 9. ⬜

In fact, it is not difficult to show that the properties given by corollary 10 uniquely determine the probability distributions, up to a constant scaling factor. This can be done by iteratively applying the statements to approximate the distributions in terms of a sum over products of inverse squares of uniform random variables. I do not go through this here though, and instead show how we can construct random variables with the given densities as a sum of gamma distributed random variables.

Recall that, for real ${a > 0}$ , the ${{\rm gamma}(a)}$ distribution on the positive reals is given by the probability density ${x^{a-1}e^{-x}/\Gamma(a)}$ .

Theorem 11 Let ${a > 0}$ be constant and ${Z_1,Z_2,\ldots}$ be an IID sequence of nonnegative random variables with the ${{\rm gamma}(a)}$ distribution, and set

$\displaystyle X=\frac1{\pi^2}\sum_{n=1}^\infty\frac{Z_n}{n^2}.$

This has moment generating function

$\displaystyle {\mathbb E}[e^{-\lambda^2 X}]=\left(\frac{\lambda}{\sinh\lambda}\right)^a$ (19)

for real ${\lambda}$ . In particular,

if ${a=1}$ then ${\sqrt{\pi X}}$ has probability density ${\Phi}$ .

if ${a=2}$ then ${\sqrt{\pi X}}$ has probability density ${\Psi}$ .

Proof: As it has the gamma distribution with parameter ${a}$ , ${Z_n}$ has moment generating function

$\displaystyle {\mathbb E}[e^{-u Z_n}]=\Gamma(a)^{-1}\int_0^\infty z^{a-1}e^{-(1+u)z}dz=(1+u)^{-a}.$

So, using independence of the ${Z_n}$ , we compute

$\displaystyle {\mathbb E}[e^{-\lambda^2X}] =\prod_{n=1}^\infty{\mathbb E}[e^{-\lambda^2\pi^{-2}n^{-2}Z_n}] =\prod_{n=1}^\infty(1+\lambda^2\pi^{-2}n^{-2})^{-a}.$

Substituting in the product expansion for ${\sinh}$

$\displaystyle \sinh x=x\prod_{n=1}^\infty\left(1+\frac{x^2}{n^2\pi^2}\right)$

gives (19) as required.

Finally, setting ${Y=\sqrt{\pi X}}$ then, by what we have shown above,

$\displaystyle {\mathbb E}[e^{-\pi^{-1}\lambda^2Y^2}]=(\lambda/\sinh\lambda)^a.$

Using the fact that the distribution is uniquely determined by the moment generating function, lemma 8 says that Y has density ${\Phi}$ if ${a=1}$ and ${\Psi}$ if ${a=2}$ . ⬜

Theorem 11 enables us to describe the distribution of the sample standard deviation of a Brownian bridge in terms of the probability density ${\Phi}$ , as was promised above.

Lemma 12 Let B be a standard Brownian bridge with sample mean ${\bar B = \int_0^1 B_tdt}$ and sample variance

$\displaystyle V=\int_0^1(B_t-\bar B)^2dt$

Then, ${\sqrt{2\pi V}}$ has probability density ${\Phi}$ .

Proof: By the Fourier expansion of the Brownian bridge,

$\displaystyle B_t -\bar B = \sum_{n=1}^\infty\frac1{\sqrt2\pi n}(Y_n\cos(2\pi nt)+Z_n\sin(2\pi nt))$

for IID standard normals ${Y_n,Z_n}$ . By Parseval’s theorem,

$\displaystyle V = \sum_{n=1}^\infty\frac1{4\pi^2 n^2}(Y_n^2+Z_n^2).$

As ${(Y_n^2+Z_n^2)/2}$ are IID ${{\rm gamma}(1)}$ random variables, the result is given by theorem 11. ⬜

A Family of Distributions

In light of theorem 11, we see that the distributions introduced above with densities ${\Phi}$ and ${\Psi}$ are just two out of an infinite family, one for each positive real number ${a}$ . I briefly look at this family, and show how they appear as integrals of Bessel bridges.

Definition 13 For real ${a > 0}$ , a nonnegative random variable X will be said to have distribution ${\mathcal S_a}$ if

$\displaystyle {\mathbb E}[e^{-\frac12\lambda^2 X}]=\left(\frac{\lambda}{\sinh\lambda}\right)^a$

for all real ${\lambda}$ .

See the 2001 paper, Infinitely divisible laws associated with hyperbolic functions by Pitman and Yor for a study of this and other related families of distributions. Using this definition, theorem 11 states the following.

Lemma 14 A nonnegative random variable X

has distribution ${\mathcal S_1}$ if and only if ${\sqrt{\pi X/2}}$ has density ${\Phi}$ .

has distribution ${\mathcal S_2}$ if and only if ${\sqrt{\pi X/2}}$ has density ${\Psi}$ .

Expressing this back in terms of the Riemann zeta function, (8) gives the moments

$\displaystyle {\mathbb E}[X^s]=2^{1+s}s(2s-1)\pi^{-2s}\Gamma(s)\zeta(2s)$

for ${s\in{\mathbb C}}$ and X with distribution ${\mathcal S_2}$ . Similarly, equation (15) gives

$\displaystyle {\mathbb E}[X^s]=2^{1+s}s(1-2^{1-2s})\pi^{-2s}\Gamma(s)\zeta(2s)$

when X has distribution ${\mathcal S_1}$ .

Theorem 11 also realizes the distributions ${\mathcal S_a}$ as sums of gamma distributed random variables.

Lemma 15 Let ${a > 0}$ be constant and ${Z_1,Z_2,\ldots}$ be an IID sequence of nonnegative random variables with the ${{\rm gamma}(a)}$ distribution. Then,

$\displaystyle \frac2{\pi^2}\sum_{n=1}^\infty\frac{Z_n}{n^2}.$

has distribution ${\mathcal S_a}$ .

It is standard that the moment generating function of the sum of a pair of nonnegative variables is equal to the product of their moment generating functions if and only if they are independent. So, definition 13 immediately gives the following result.

Lemma 16 Let X and Y be independent with distributions ${\mathcal S_a}$ and ${\mathcal S_b}$ respectively. Then ${X+Y}$ has distribution ${\mathcal S_{a+b}}$ .

An example of the occurrence of law ${\mathcal S_{1/2}}$ is provided by the integral of a squared Brownian bridge.

Lemma 17 Let B be a Brownian bridge. Then ${\int_0^1B_t^2dt}$ has distribution ${\mathcal S_{1/2}}$ .

Proof: We use the sine series representation of the Brownian bridge,

$\displaystyle B_t=\frac{\sqrt2}{\pi}\sum_{n=1}^\infty\frac{X_n}{n}\sin(\pi n t)$

where ${X_1,X_2,\ldots}$ is an IID sequence of standard normals. By Parseval’s theorem,

$\displaystyle \int_0^1B_t^2dt=\frac2{\pi^2}\sum_{n=1}^\infty\frac{X_n^2}{2n^2}.$

However, ${X_n^2/2}$ are IID ${{\rm gamma}(1/2)}$ random variables, so the result follows from lemma 15. ⬜

Finally, I show how all of the laws ${\mathcal S_a}$ occur as integrals of squared Bessel bridges. Recall that for positive integer a, the sum of the squares of a independent Brownian motions is a squared Bessel process of order a. This generalizes to noninteger ${a > 0}$ , and such processes are denoted as ${{\rm BES}^2_0(a)}$ . I will not be concerned with a detailed description here, but note that they have the properties:

If B is standard Brownian motion, then ${B^2}$ is a ${{\rm BES}^2_0(1)}$ process.
If independent processes X,Y are respectively ${{\rm BES}^2_0(a)}$ and ${{\rm BES}_0^2(b)}$ , then ${X+Y}$ is a ${{\rm BES}_0^2(a+b)}$ process.

Bessel bridges are Bessel processes conditioned on hitting zero at time ${t=1}$ . Again, I will not be concerned with a full description, and just require the following properties. Using the terminology ‘ ${{\rm BES}^2_0(a)}$ bridge’ for a squared Bessel bridge of order a:

If B is a Brownian bridge, then ${B^2}$ is a ${{\rm BES}^2_0(1)}$ bridge.
If independent processes X,Y are respectively ${{\rm BES}^2_0(a)}$ and ${{\rm BES}^2_0(b)}$ bridges, then ${X+Y}$ is a ${{\rm BES}^2_0(a+b)}$ bridge.

Integrals of squared Bessel bridges give random variables in the family ${\mathcal S_a}$ .

Theorem 18 For any ${a > 0}$ , let X be a ${{\rm BES}^2_0(a)}$ bridge. Then ${\int_0^1 X_tdt}$ has distribution ${\mathcal S_{a/2}}$ .

Proof: Fixing ${\lambda\in{\mathbb R}}$ , define ${f\colon{\mathbb R}_{ > 0}\rightarrow{\mathbb R}_{ > 0 }}$ by

$\displaystyle f(a) = {\mathbb E}\left[e^{-\frac12\lambda^2\int_0^1X_tdt}\right],$

where X is a ${{\rm BES}^2_0(a)}$ bridge, defined on some probability space. As the sum of independent ${{\rm BES}^2_0(a)}$ and ${{\rm BES}^2_0(b)}$ bridges is a ${{\rm BES}^2_0(a+b)}$ bridge, we have ${f(a+b)=f(a)f(b)}$ . The only solutions to this functional equation on the positive reals are of the form ${f(a)=C^a}$ for constant ${C > 0}$ .

As the square of a Brownian bridge is a ${{\rm BES}^2_0(1)}$ bridge, lemma 17 gives ${f(1)=(\lambda/\sinh\lambda)^{1/2} = C}$ . We have shown that, for a ${{\rm BES}^2_0(a)}$ bridge X,

$\displaystyle {\mathbb E}\left[e^{-\frac12\lambda^2\int_0^1X_tdt}\right]=f(a)=\left(\frac{\lambda}{\sinh\lambda}\right)^{\frac a2}$

as required. ⬜

One thought on “The Riemann Zeta Function and Probability Distributions”

jack morava says:

3 July 22 at 9:02 PM

Wu’s lemma asserts that \\half(\exp(2x – 1) has 2-adic integral coefficients and reduces mod two to \sum_{k geq 0} x^{2^k) , cf https://arxiv.org/abs/1608.04702 lemma 2.1.3. May I send you a short note related this? It seems to bw related to the Stefan-Boltzmann relation in statistical mechanics \dots

	Harry on Filtrations and Adapted P…
	Anonymous on Analytic Sets
	Daniel on Analytic Sets
	Sid on Martingale Convergence
	northwiz on Preservation of the Local Mart…
	George Lowther on Martingale Convergence
	Sid on Martingale Convergence
	George Lowther on The Burkholder-Davis-Gundy Ine…
	George Lowther on Preservation of the Local Mart…
	northwiz on Preservation of the Local Mart…

The Riemann Zeta Function and Probability Distributions

Another Distribution

Properties of the Distributions

A Family of Distributions

Published by George Lowther

One thought on “The Riemann Zeta Function and Probability Distributions”

Leave a Reply Cancel reply

Another Distribution

Properties of the Distributions

A Family of Distributions

Like this:

Related

Published by George Lowther

One thought on “The Riemann Zeta Function and Probability Distributions”

Leave a Reply Cancel reply