What is the theoretical value of this phase space invariant? So I wanted know how to theoretically calculate this phase space invariant (equation $3.31a$ )$R$ in our universe (FLRW metric) during the cosmological nucleosysthesis:
$$R = \int_{p} \frac{\mathcal{N}(x,p,t)}{\mathcal{E}} dp_x dp_y dp_z$$
where $\mathcal{N}(x,p,t)$ is the phase space distribution of matter in our universe, $\mathcal{E}= p_0$ is the energy and $p_x$, $p_y$ and $p_z$ are the spatial momentum.
Feel free to assume extra assumptions (but assume the formation of the molecules has taken place)?
 A: I'm going to start by using some relations that appear explicitly in the text you linked. Specifically:
\begin{equation}
\mathcal{N} = \frac{g_{s}}{h^{3}} \eta
\end{equation}
where, for the era of nucleosynthesis, the Universe is matter-dominated ($\rho \propto a^{-3} \, , \, P=0$ with $a$ being the time-dependent scale factor), and thus the particle statistics is dominated by non-relativistic protons, electrons and Helium nuclei ($g_{s} = 2s +1$ with $s$ being typically $1/2$). Since matter is non-relativistic and at that energy level behave roughly classically, the particles follow Boltzmann statistics as such:
\begin{equation}
\eta \approx e^{-\beta (E - \mu)} = e^{\beta \mu} \, e^{-\beta p^{2}/2m}
\end{equation}
where we defined $\beta \equiv (k_{B}T)^{-1}$ and approximated the kinetic energy to the Newtonian form. The integral for the field $R$ now reads:
\begin{equation}
R = \frac{g_{s}}{h^{3}} e^{\beta \mu} \int \frac{e^{-\beta p^{2}/2m}}{\frac{p^{2}}{2m} + m} \, d\mathcal{V}_{p} = \frac{4\pi g_{s}}{h^{3}} e^{\beta \mu} \int _{0} ^{\Lambda} \frac{p^{2} e^{-\beta p^{2}/2m}}{\frac{p^{2}}{2m} + m} \, dp
\end{equation}
The parameter $\Lambda$ is a cutoff we introduce since the momentum of each particle has a form of cap due to the assumption they behave non-relativistically. Next we can use the trick:
\begin{equation}
p^{2} e^{-\beta p^{2}/2m} = -2m \frac{\partial }{\partial \beta} e^{-\beta p^{2}/2m}
\end{equation}
so the integral now reads:
\begin{equation}
R = -4m^{2} \frac{4\pi g_{s}}{h^{3}} e^{\beta \mu} \frac{\partial}{\partial \beta} \int _{0} ^{\Lambda} \frac{e^{-\beta p^{2}/2m}}{p^{2} + 2m^{2}} \, dp
\end{equation}
The integral now can be calculated using the Cauchy integral formula by choosing an appropriate contour that circumvents the two imaginary poles $p = \pm i\sqrt{2} m$:
\begin{equation}
\oint  \frac{e^{-\beta p^{2}/2m}}{(p+i\sqrt{2}m)(p-i\sqrt{2}m)} \, dp = 2\pi i \left [ \frac{e^{\beta m}}{i2\sqrt{2}m} - \frac{e^{\beta m}}{-i2\sqrt{2}m} \right ] = \sqrt{2} \pi \frac{e^{\beta m}}{m}
\end{equation}
The effect of the derivative with respect to $\beta$ only acts on the exponential which yields $m \, e^{\beta m}$, so the final result is:
\begin{equation}
R = -\frac{16\pi g_{s} m^{2}}{\sqrt{2}h^{3}} e^{\beta \tilde{\mu}} 
\end{equation}
where as defined in the text $\tilde{\mu} = \mu + m$. The result comes out negative, but since the quantity $R$ is a scalar field and not a positive-definite quantity like energy or the number of particles, I expect the result should be allowed.
EDIT:
As an addendum and in response to your comment about not taking the non-relativistic limit, let me provide some additional calculations. Let's assume that we are at some point during the radiation-dominated era of the Universe where most particles behave relativistically. This time they can be either photons/other massless bosons or hard fermions. The coupling between them is very weak due to being of high energy, so $\beta \mu \ll 1$ and the particle statistics can be both Bose-Einstein or Fermi-Dirac:
\begin{equation}
\eta = \frac{1}{e^{\beta(E - \mu)} \pm 1} \approx \frac{1}{e^{\beta p} \pm 1}
\end{equation}
where we further assumed $p^{2} \gg m^{2} \Rightarrow E \approx \mathcal{E} \approx p$. The scalar field $R$ will now be:
\begin{equation}
R = \frac{g_{s}}{h^{3}} \int \frac{1}{e^{\beta(E - \mu)} \pm 1} \frac{1}{p} \, d\mathcal{V}_{p} = \frac{4\pi g_{s}}{h^{3}} \int _{\lambda} ^{\Lambda} \frac{p}{e^{\beta(E - \mu)} \pm 1} \, dp
\end{equation}
As I mentioned in the comments, in addition to the upper cutoff $\Lambda$ (which you may freely now assume to be very large/approaching infinity), you impose a lower cutoff $\lambda$ which restricts the momentum of your particles to energies that would qualify as in the relativistic regime. The result of the integral had we taken no limits would be:
\begin{equation}
\int \frac{p}{e^{\beta p} \pm 1} \, dp = 
\begin{cases}
&\frac{p \, \log{(1 - e^{-\beta p})}}{\beta } - \frac{\mathrm{Li}_{2}(e^{-\beta p})}{\beta ^{2} } \, , \qquad \text{(EB)} \\
&\frac{\mathrm{Li} _{2}(e^{-\beta p})}{\beta ^{2} } - \frac{p \, \log{(1 + e^{-\beta p})}}{\beta } \, , \qquad \text{(FD)} \\
\end{cases}
\end{equation}
where $\mathrm{Li}_{n}$ is the polylogarithmic function. For the upper limit $\Lambda \rightarrow +\infty$:
\begin{equation}
\mathrm{Li}_{2}(0) = 0
\end{equation}
\begin{equation}
\lim _{p \rightarrow +\infty} {p \, \log{(1 \pm e^{-\beta p})}} = 0
\end{equation}
So all is left are the terms with the lower limit supplanted:
\begin{equation}
\int _{\lambda}^{\Lambda} \frac{p}{e^{\beta p} \pm 1} \, dp = 
\begin{cases}
&-\frac{\lambda \, \log{(1 - e^{-\beta \lambda})}}{\beta } + \frac{\mathrm{Li}_{2}(e^{-\beta \lambda})}{\beta ^{2} } \, , \qquad \text{(EB)} \\
&-\frac{\mathrm{Li} _{2}(e^{-\beta \lambda})}{\beta ^{2} } + \frac{\lambda \, \log{(1 + e^{-\beta \lambda})}}{\beta } \, , \qquad \text{(FD)} \\
\end{cases}
\end{equation}
For relatively large $\lambda$, Fermi-Dirac statistics leads to $R \approx 0$, but negative due to having a higher power of $\beta$ in the denominator in the first term (remember that $T$ in the radiation-dominated era is much higher than in the matter-dominated one). For Bose-Einstein statistics, had the lower limit been 0, the logarithm would tend towards $-\infty$, which is where I presume the divergence you found occurs. But for relatively large values $\lambda$, the logarithm is now regulated and is expected to be something close to 0. The first term will also be very small, so to determine whether $R$ is positive or negative you'd have to input more precise values of your parameters.
A: After pondering for a while, I thought it might be worth sharing a different angle to this problem. So let's start with the stress energy tensor $T^{\mu \nu}$ of a perfect fluid:
$$T^{\mu \nu} = \left(\rho + \frac{p}{c^2} \right) U^{\mu} U^\nu + p g^{\mu \nu}, $$
where $p$ is pressure, $\rho$ is density, $U^\mu$ is four velocity and $g^{\mu \nu}$ is the metric. Consider the trace of this quantity:
$$T \equiv T^{i}_i = (\rho + \frac{p}{c^2})   c^2+ 4p =3p + \rho c^2 $$
Now, we know for a particular regime this parameter is approximately constant. Hence, we have a conserved quantity or equation of state:
$$ 3p + \rho c^2 = \text{constant}$$
However, we also know from the linked text's equation $(3.31 c)$ in the question. The stress energy tensor is given by:
$$ T^{\mu \nu} = \int \mathcal{N}(x,p,t) p^\mu \otimes p^\nu \frac{d V_p}{E}$$
Taking the trace as mentioned before we get and using $m$ as mass:
$$ \text{constant}=3p + \rho c^2= m^2 c^4 \int \mathcal{N}(x,p,t)  \frac{d V_p}{E} $$
Thus,
$$  \int \mathcal{N}(x,p,t)  \frac{d V_p}{E} =\frac{1}{m^2c^4} (3p + \rho{c^2} ) = \frac{\text{constant}}{m^2c^4} $$
