$\rho-3P$ for ultra-relativistic regime in early universe I am trying to read this paper, in which they try to get dark energy by modifying the EFEs to its unimodular form. As interesting as it may be for someone here, I'm struggling to understand the equation $(12)$ in the paper. Which reads as follows $$\rho-3P\approx\frac{m_t^2T^2}{2\hbar^3},$$ here $m_t$ is the top mass and $T$ is temperature of the universe. The paper claims the relationship should follow from standard thermodynamics. However, I am at a loss in understanding how to arrive at this expression. A derivation or even a hint for how to arrive at a relation like $\rho-3P\approx\frac{m^2T^2}{2\hbar^3}$ up to first order of approximation from standard thermodynamics would be appreciated.
 A: In statistical mechanics the energy density and pressure have the following definitions,
$$\rho=g\int\frac{d^3p}{(2\pi\hbar)^3}Ef(\vec{x},\vec{p},t)\\
P=g\int\frac{d^3p}{(2\pi\hbar)^3}\frac{p^2}{3E}f(\vec{x},\vec{p},t).$$
Here, $g$ is the degeneracy factor and $f(\vec{x},\vec{p},t)$ is the distribution function of particles in phase space and in thermal equilibrium will be equal to either Bose-Einstein or Fermi-Dirac distribution for bosons and fermions respectively,
$$f=\frac{1}{\exp\left[\frac{E}{T}\right]\pm1}.$$
For massless particles, $|\vec{p}|=E$ and therefore $\rho-3P=0$. For ultrarelativistic particles which are in a very high temperature (conditions resembling the early universe) $\frac{m}{T}\ll1$, the distribution fuction still behaves like that of a massless one,
$$f=\frac{1}{\exp\left[\frac{p}{T}\right]\pm1},$$ however, we can obtain the first order correction for massive particles by expanding $E$ for ultrarelativistic particles of small mass,
$$E=\sqrt{p^2+m^2}\approx p(1+\frac{m^2}{2p^2})=p+\frac{m^2}{2p}.$$
The first order correction in energy density is,
$$\rho^{(1)}=\frac{gm^2}{2(2\pi\hbar)^3}\int_0^\infty dp4\pi p^2\frac{\frac{1}{p}}{\exp\left[\frac{p}{T}\right]\pm1}=\frac{gm^2}{4\pi^2\hbar^3}\int_0^\infty dp\frac{p}{\exp\left[\frac{p}{T}\right]\pm1}=\frac{gm^2T^2}{4\pi^2\hbar^3}\int_0^\infty dy\frac{y}{e^y\pm1}\propto m^2T^2.$$
One must perform the zeta function integral to get the correct numerical values. For pressure expand as follows,
$$\frac{p^2}{3E}\approx \frac{p}{3}(1-\frac{m^2}{2p^2}).$$
And similarly one can find the first order correction in pressure which will be of the order of $m^2T^2$. Therefore, for ultrarelativistic massive particles,
$$\rho-3P\approx\mathcal{O}[m^2T^2].$$
I leave it as an exercise to the OP to calculate the correct numerical factors.
Mathematical aid
Let,
\begin{align}
    I_n^\pm=\int_0^\infty\frac{x^n}{e^x\pm1}dx.
\end{align}
Now,
\begin{align}
    I_n^--I_n^+=\int_0^\infty x^n\frac{2}{e^{2x}-1}dx=2^{-n}\int_0^\infty\frac{y^n}{e^y-1}dy=2^{-n}I_n^-.
\end{align}
Thus,
\begin{align}
    I_n^+=\left(1-\frac{1}{2^n}\right)I_n^-.
\end{align}
And,
\begin{align}
    I_n^-=\int_0^\infty\frac{x^n}{e^x-1}dx=\Gamma(n+1)\zeta(n+1).
\end{align}
