Net particle number density for relativistic particles at finite chemical potential (tricky integral)

Question

Question: How does one show that the chemical potential of relativistic fermions is negligible at high energies? In particular, I would like to show that the difference between the particle density $n$ and the anti-particle density $\bar n$ is $$ n - \bar n = \frac{gT^3}{6\pi^2}\left[ \pi^2 \left( \frac{\mu}{T} \right) + \left(\frac{\mu}{T}\right)^3 \right] \tag{*}\label{main} $$ This relation is presented in Baumann's Cosmology textbook, Problem 3.1, which notes that it is exact in the $m/T \to 0$ limit. (That is: there is no expansion in $\mu/T$.) It also appears in Mukhanov's Physical Foundations of Cosmology in equation (3.54).

In the following, I give a bit of background to write the expression for $n(T)$ and establish conventions. Then I present three attempts, though in each of these attempts I have been unable to prove the exact form in \eqref{main}. I seek advice on how to fix these approaches so that they are consistent with \eqref{main} and one another.

Background

$n$ and $\bar n$ are the number densities of the particles and antiparticles, which may be written in terms of the phase space distributions as

$$ n(T) = g \int \frac{d^3\vec{p}}{(2\pi)^3} f(p, T) = \frac{g}{2\pi^2} \int_0^\infty \frac{dp\, p^2}{e^{(E(p) - \mu)/T}+1} $$ $g$ is the number of degrees of freedom for the particle, $p = |\vec{p}|$ is the magnitude of the 3-momentum, and in the relativistic limit we have $E(p) = p$. The antiparticle number density is nearly indentical, except that the antiparticle chemical potential is minus that of the particle, $\bar\mu = - \mu$.

We write the expression for $n-\bar n$ using dimensionless variables: $x = E/T$, $\beta = \mu/T$: $$ \displaystyle n-\bar n = \frac{g}{2\pi}T^3 \int_0^\infty dx \left[ \frac{x^2}{e^{x-\beta} + 1} - \frac{x^2}{e^{x+\beta} + 1} \right] \tag{a}\label{dimless} $$ In this post I present three approaches to show that \eqref{dimless} matches \eqref{main}. However, I am unable to get any approach to exactly match \eqref{main} and I would like to find the mistake(s) that makes all three approaches consistent with \eqref{main}.

Baumann approach

Baumann's problem 3.1 offers the following hint: $$ \int_0^\infty dy \frac{y}{e^y + 1} = \frac{\pi^2}{12} \tag{**}\label{hint} $$

Using this, I am able to reproduce the first term in $\eqref{main}$ as follows:

$$ n-\bar n = \frac{g}{2\pi}T^3 \left[ \int_0^\infty dx \frac{x^2}{e^{x-\beta} + 1} - \int_0^\infty dx \frac{x^2}{e^{x+\beta} + 1} \right] $$ $$ =\frac{g}{2\pi}T^3 \left[ \int_{-\beta}^\infty dy \frac{(y+\beta)^2}{e^{y} + 1} - \int_{\beta}^\infty dy \frac{(y-\beta)^2}{e^{y} + 1} \right] $$ $$ =\frac{g}{2\pi}T^3 \left[ \int_{-\beta}^\infty dy \frac{y^2 + 2y\beta + \beta^2}{e^{y} + 1} - \int_{\beta}^\infty dy \frac{y^2 - 2y\beta + \beta^2}{e^{y} + 1} \right] $$

At this stage one may regroup like terms in the brackets: $$ \left[\cdots \right] = \int_{-\beta}^\beta dy \frac{y^2 + \beta^2}{e^{y} + 1} + \int_{-\beta}^0 dy \frac{2y\beta}{e^{y} + 1} - \int_{0}^\beta dy \frac{2y\beta}{e^{y} + 1} + \int_{0}^\infty dy \frac{4y\beta}{e^{y} + 1} $$

The first three terms are all $\mathcal O(\beta^3)$, which one may see from the integrand and the region of integration. The last term is $\mathcal O(\beta)$ and we may apply the hint \eqref{hint} to find:

$$ n-\bar n = \frac{gT^3}{2\pi^2} \left[\frac{4\pi^2 \beta}{12} + \mathcal O(\beta^3)\right] $$

Recalling that $\beta = \mu/T$, the leading order term indeed matches \eqref{main}. However, the remaining terms seem to be cumbersome integrals. While we can at least see that the remaining terms are $\mathcal O(\beta^3)$ as requried in the main result \eqref{main}, it is not clear to me how to evaluate these cleanly.

Reif approach

An alternative approach is motivated by Appendix A.11 of Reif's Fundamentals of statistical and thermal physics. Here we write the integrals in \eqref{dimless} as series expansions:

$$ \frac{x^2}{e^{x-\beta}+1} = x^2 \frac{e^{\beta-x}}{1-(-e^{\beta - x})} = x^2 e^{\beta-x} \sum_{n=0}^\infty (-e^{\beta - x})^n = - x^2 \sum_{n=1}^\infty (-e^{\beta - x})^n $$

The second integrand in \eqref{dimless} is related by $\beta \to -\beta$ and an overall sign.

The integral over $dx$ simplifies on each term in the power series: $$ \int_0^\infty dx \frac{x^2}{e^{x-b}+1} = -\sum_{n=1}^\infty (-)^n e^{n\beta} \int_0^\infty dx\, x^2 e^{-nx} $$ $$= -\sum_{n=1}^\infty (-)^n e^{n\beta} \frac{1}{n^3} \int_0^\infty du\, u^2 e^{-u} $$ $$= -2\sum_{n=1}^\infty (-)^n \frac{e^{n\beta}}{n^3} $$ where in the last step we have used the definition of the $\Gamma$ function on integers, $\Gamma(n) = \int_0^\infty dx\, x^n e^{-x} = n!$.

The second integral in \eqref{dimless} is similarly---where we now include the minus sign from the difference $n-\bar n$: $$ - \int_0^\infty dx \frac{x^2}{e^{x+b}+1} = +2\sum_{n=1}^\infty (-)^n \frac{e^{-n\beta}}{n^3} = -2\sum_{n=-1}^{-\infty} (-)^n \frac{e^{n\beta}}{n^3} \ , $$ where in the last step we changed the sum from $n$ over positive integers to negative integers.

Together this lets us write the dimensionless integral in \eqref{dimless} as $$ \int_0^\infty dx \left[ \frac{x^2}{e^{x-\beta} + 1} - \frac{x^2}{e^{x+\beta} + 1} \right] = -2\sum_{n\neq 0} \frac{e^{n\beta}}{n^3}e^{i\pi n} $$ where $e^{i\pi n} = (-)^n$. Now there is a clever way to convert the sum over non-zero integers into a contour integral around the origin. The key observation is that $(\tan\pi z)^{-1}$ is a function with simple poles at every integer. Because the terms in the sum above have no singularities at the non-zero integers $n$, we may write the sum as a contour integral over an appropriate contour that encloses all non-zero real integers. All we have to do is multiply the summand/integrand by $(\tan\pi z)^{-1}$ because its poles will pick out the summand evaluated on each non-zero integer. This is shown in Reif Appendix A.11, whose key figure is shown below.

Writing this out:

$$ -2\sum_{n\neq 0} \frac{e^n\beta}{n^3}e^{i\pi n} = \frac{-1}{2i} \oint_C \frac{2e^{z\beta}e^{i\pi z}}{z^3\, \tan\pi z} $$

Next we recognize that we may deform the contour so that it encloses the origin. There's a minus sign coming from the orientation of the new contour $C'$:

$$ -2\sum_{n\neq 0} \frac{e^n\beta}{n^3}e^{i\pi n} = \frac{+1}{2i} \oint_{C'} \frac{2e^{z\beta}e^{i\pi z}}{z^3\, \tan\pi z} = \frac{+1}{2i} \text{Res}\left(\frac{2e^{z\beta}e^{i\pi z}}{z^3\, \tan\pi z},0\right) $$

I have most of the ingredients but numerical factors do not work out. We write out:

$$ \tan\pi z = \pi z \left(1 + \frac{1}{2}\pi^2z^2 + \cdots\right) $$ $$e^{z\bar\beta}\equiv e^{z(\beta+i\pi)} = 1+ z\bar\beta + \frac{1}{2}z^2 \bar\beta^2 + \frac{1}{6}z^3\bar\beta^3 + \cdots $$

This allows us to write $$ \frac{e^{z\bar\beta}}{z^3 \tan\pi z} = \frac{1+ z\bar\beta + \frac{1}{2}z^2 \bar\beta^2 + \frac{1}{6}z^3\bar\beta^3}{\pi z^4\left( 1 + \frac{1}{2}\pi^2z^2 + \cdots \right)} $$ From which we may identify the terms with simple poles: $$ \frac{e^{z\bar\beta}}{z^3 \tan\pi z} = \frac{1}{z}\left[ \frac{-\pi\bar\beta}{3} + \frac{\bar\beta^3}{6\pi}\right] $$ When plugging this into \eqref{dimless} this almost gives the correct expression except:

• I am assuming $\bar\beta \to \beta$, which ignores an $i\pi$ piece that game from an $(-)^{n} = e^{i\pi n} \to e^{i\pi z}$ factor, which in turn came from having to insert a minus sign to use the geometric series.
• Even then, the $\mathcal O(\beta)$ term is off by a factor of $-1/2$.

So this approach seems to be close. I may have been to sloppy finding the residue at the origin, but I'm not sure how to deal with the $\bar\beta = \beta + i\pi$ or the relative sign between the $\mathcal O(\beta)$ and $\mathcal O(\beta^3)$ terms.

Mukhanov approach

A third approach comes from Mukhanov's Principles of Physical Cosmology Sections 3.3.3 - 3.3.4 with the result from Baumann Problem 3.1 explicitly stated in equation (3.54). Mukhanov introduces a class of integrals

$$ J^{(\nu)}_\mp (\alpha,\beta) \equiv \int_\alpha^\infty dx\, (x^2 -\alpha^2)^{\nu/2} \left(\frac{1}{e^{x-\beta}\mp 1} + \frac{1}{e^{x+\beta}\mp 1}\right) $$

The parameter $\alpha = m/T$ can be set to 0 since we take the $m \ll T$ limit. In Mukhanov's equation (3.37) he proposes that the net particle density for fermions is

$$ n-\bar n = \frac{gT^3}{6\pi^2} \frac{\partial J^{(3)}_+}{\partial \beta} $$

This is posed as problem 3.8 in Mukhanov's text. I believe the easy way to validate this expression is to refer to an earlier observation that $n = \partial p/\partial \mu$, the derivative of the pressure with respect to the chemical potential. This is shown straightforwardly earlier in the book (it follows from equation3.32). It is then straightforward to show that the total pressure $p+\bar p$ is proportional to $J_+^{(3)}$ and that taking a $\partial/\partial \beta$ derivative introduces a relative sign between the particle and anti-particle terms so that one ends up with $n-\bar n$.

What is vexing to me---and mentioned on a related stackexchange question---is how to see that a $\partial/\partial \beta$ derivative of $J^{(3)}_+(0,\beta)$ could possibly give something that matches \eqref{dimless}. This seems to require that

$$ \int_0^\infty dx\, x^2 \left(\frac{1}{e^{x-\beta}+1} - \frac{1}{e^{x+\beta}+1}\right) = \int_0^\infty dx\, x^3 \frac{\partial}{\partial \beta} \left(\frac{1}{e^{x-\beta}+1} + \frac{1}{e^{x+\beta}+1}\right) $$

The right-hand integrand does not obviously simplify: $$ x^3 \frac{\partial}{\partial \beta} \left(\frac{1}{e^{x-\beta}+1} + \frac{1}{e^{x+\beta}+1}\right) = x^3 \left( \frac{e^{x-\beta}}{(e^{x-\beta}+1)^2} - \frac{e^{x+\beta}}{(e^{x+\beta}+1)^2} \right) $$ And somehow I would like this to become $$ \int_0^\infty dx\, x^3 \left( \frac{e^{x-\beta}}{(e^{x-\beta}+1)^2} - \frac{e^{x+\beta}}{(e^{x+\beta}+1)^2} \right) \stackrel{?}{=} \int_0^\infty dx\, x^2 \left( \frac{1}{e^{x-\beta}+1} - \frac{1}{e^{x+\beta}+1} \right) $$

Again, Mukhanov offers a plausible explanation using pressure as an intermediate step: first, net pressure is related to $J^{(3)}_+$, then the book separately shows that the number density is a $\mu$ derivative of the pressure. However, I would like to see how the derivative of $J^{(3)}_+$ explicitly can be messaged into the form of the net particle density ($n-\bar n$) integral in \eqref{dimless}. Once that is clear, I am happy with the limiting forms of the $J^{(3)}_+$ function used later in the text to show $n-\bar n$ in the $m\ll T$ limit.

Answer 1

Your question is purely mathematical. I will set $\mu=1$ with no loss of generality in the following. If you think in terms of fugacity: $$ z = e^\beta $$ you can use polylogarithms: $$ \text{Li}_s(x) = \sum_{n=1}^\infty\frac{x^n}{n^s} = \frac1{\Gamma(s)}\int_0^\infty \frac{xe^{-t}}{1-xe^{-t}}t^{s-1}dt $$

In $D$ dimensions, for ultra relativistic particles, your density of states is: $$ D(E) = \frac{2g}{(4\pi)^{D/2}\Gamma(D/2)}E^{D-1} $$ so let $s=D$ and: $$ C = \frac{2\Gamma(D)g}{(4\pi)^{D/2}\Gamma(D/2)}T^D $$ You have in general: $$ \begin{align} \frac nC &= -\text{Li}_s(-e^\beta) \\ \frac{\bar n}C &= -\text{Li}_s(-e^{-\beta}) \end{align} $$ Since you are interested in the difference, you just want the odd part of the first function (odd terms of the Taylor series). In your case, $s=D=3$.

A quick and dirty way is to use the power series expansion of $\text{Li}$: $$ \begin{align} f(\beta) &:= -\text{Li}_s(-e^\beta) \\ &= \sum_{n=1}^\infty \frac{(-1)^{n-1}}{n^s}e^{n\beta} \\ &= \sum_{n=1}^\infty\sum_{m=0}^\infty \frac{(-1)^{n-1}}{n^s}\frac{(n\beta)^m}{m!}\\ &= \sum_{m=0}^\infty\left(\sum_{n=1}^\infty\frac{(-1)^{k-1}}{n^{s-m}}\right)\frac{\beta^m}{m!}\\ &= \sum_{m=0}^\infty \eta(s-m)\frac{\beta^m}{m!} \end{align} $$ using the Dirichlet eta function. The manipulations are formal: you cannot use Fubini since the double series is not summable. You can simply tweak the calculation to make it rigorous by introducing $x$ such that $|x|<1$: $$ \begin{align} f(\beta,x) &:= -\text{Li}_s(-xe^\beta) \\ &= \sum_{n=1}^\infty \frac{(-1)^{n-1}}{n^s}x^ne^{n\beta} \\ &= \sum_{n=1}^\infty\sum_{m=0}^\infty \frac{(-1)^{n-1}}{n^s}x^n\frac{(n\beta)^m}{m!}\\ &= \sum_{m=0}^\infty\left(\sum_{n=1}^\infty\frac{(-1)^{k-1}}{n^{s-m}}x^n\right)\frac{\beta^m}{m!}\\ &= \sum_{m=0}^\infty (-\text{Li}_{s-m}(-x))\frac{\beta^m}{m!} \end{align} $$ This time, the series is summable and you can apply Fubini. By sending $x\to1$, you recover the definition of $\eta$ as the Abel summation: $$ \eta(s-m) = -\text{Li}_{s-m}(-1) $$ In your case, you just need to use: $$ \eta(2) =\frac{\pi^2}6 \quad \eta(0) = \frac12 $$ which gives you the correct leading order asymptotics (the relative $\pi^2$ factor).

You can also prove this using Reif's approach. You should not have introduced $\bar\beta$, but rather expand in $\beta$. This time, I will assume $s\in\mathbb Z$. It is a standard approach, you may have seen it to prove the famous relationship between Bernoulli's numbers and the zeta function on the positive even numbers: $$ \zeta(2n) = \frac{B_{2n}(2\pi)^{2n}}{2(2n)!} \tag1 $$ where the Bernoulli numbers are defined by the generating function: $$ \frac x{e^x-1} = \sum_{n=0}^\infty \frac{B_n}{n!}x^n $$ I find it clearer to think in terms of exponentials rather than $\tan$. You start with the function: $$ \frac{2\pi i}{e^{i2\pi z}-1} $$ which has simple poles at $\mathbb Z$ of residue $1$. You want to calculate the sum: $$ S = -2\sum_{n\in\mathbb Z^*}e^{n\beta}\frac{e^{in\pi}}{n^s} $$ so you naturally introduce the function: $$ \begin{align} f(z) &= -\frac1{z^s}\frac{2e^{i\pi z}}{e^{2\pi iz}-1}e^{\beta z} \\ &= \sum_{m=0}^\infty \frac i{z^{s-m}\sin(\pi z)}\frac{\beta^m}{m!} \\ &= \sum_{m=0}^\infty f_m(z)\frac{\beta^m}{m!} \end{align} $$ who still has poles only on $\mathbb Z$ and whose sum of resides of non zero poles gives $S$ (this is why I assume $s$ is an integer, otherwise a branch cut would be needed). Using the same contour: $$ \oint f_m(z)dz = 0 $$ by Jordan's lemma. On the other hand, the residues at $\mathbb Z^*$ gives $S$, so you just need to calculate the residue at $0$ of $f_m$. You just need to isolate the $1/z$ term of the Laurent expansion. Either you already know the Taylor expansion of $\tan$ / Laurent expansion of $\csc$ or you rederive it using the Bernoulli numbers and its defining generating function (this is where thinking in terms of exponentials is handy): $$ f_m(z) = 2iz^{m-s}\sum_{n=0}^\infty\frac{(-1)^{n-1}(2^{2n-1}-1)B_{2n}}{(2n)!}z^{2n-1} $$ so: $$ Res(f_m,0) = \frac{(-1)^{s-m-1}(2^{s-m-1}-1)B_{s-m}}{(s-m)!} $$ This gives you a new formula from the Taylor series, which is consistent with the previous one by using: $$ \eta(s) = (1-2^{s-1})\zeta(s) $$ and $(1)$.

On a final note, the Taylor series has a finite radius of convergence. It turns out to be $\pi$ by studying the global property of the polylogarithm (single branch point at $1$). Therefore, it is more than just an asymptotic expansion and will converge for $|\beta|<\pi$.