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.
