$$\begin{aligned}
     &&\int^\infty_0\!{\mathrm{d}}|\vec{p}_E||\vec{p}_E|^2 
     \hat{\rho}_{\Lambda\vec{p}_E}\ln\left(1+e^{-\frac{E_p(M)-\mu}{T}}\right)
     \nonumber\\
     &\!\! =\!\!&\hat{\rho}_{\Lambda}\left(\left.\frac{|\vec{p}_E|^3}{3}\ln
     \left(1+e^{-\frac{E_p(M)-\mu}{T}}\right)\right|^{\infty}_0\right) 
     -\hat{\rho}_{\Lambda}\int^\infty_0\!{\mathrm{d}}|\vec{p}_E|\frac{|\vec{p}_E|^3}{3}
     \frac{\partial}{\partial |\vec{p}_E|}\ln
     \left(1+e^{-\frac{E_p(M)-\mu}{T}}\right)
     \nonumber\\
     &\!\! =\!\!&\hat{\rho}_\Lambda\int^\infty_0\!{\mathrm{d}}|\vec{p}_E|
     \frac{|\vec{p}_E|^4n_{qM}}{3TE_p(M)}=\frac{1}{3T}\int^\infty_0\!{\mathrm{d}}|\vec{p}_E||\vec{p}_E|^4\hat{\rho}_{\Lambda\vec{p}_E}\frac{n_{qM}}{E_p(M)}\end{aligned}$$ where we used the fact the surface term disappears, we get finally $$\label{J-1med}
     J_{-1}^{\mathrm{med}}(M^2,T,\mu )=-\frac{8}{3}\int^\infty_0{\mathrm{d}}|\vec{p}_E||\vec{p}_E|^4\hat{\rho}_{\Lambda\vec{p}_E}\left(
     \frac{n_{qM}+n_{\overline{q}M}}{E_p(M)}-
     \frac{n_{q0}+n_{\overline{q}0}}{E_p(0)}\right).$$

It is important to realize that the expansion of the integral (\[J-1med\]) for small values of $T$ and at $\mu =0$ starts from the term $$\label{limit}
     J^{\mathrm{med}}_{-1}(M^2,T,0)=\frac{16}{3}\int_0^\infty
     \frac{|\vec{p}_E|^3{\mathrm{d}}|\vec{p}_E|}{1+e^{{\frac{|\vec{p}_E|}{T}}}}
     +{\cal O}(T^{\frac{5}{2}}M^{\frac{3}{2}}e^{-\frac{M}{T}})=\frac{14}{45}
     \pi^4T^4+{\cal O}(T^{\frac{5}{2}}M^{\frac{3}{2}}e^{-\frac{M}{T}}).$$ This leading contribution in $T$ arises from the combination $\sim (n_{q0}+
n_{\overline{q}0})$ and does not depend on the cut-off, i.e. the Pauli-Villars regulator $\hat{\rho}_{\Lambda\vec{p}_E}$ does not affect the leading order of the low-temperature asymptotics of the integral. To make this clear let us consider the typical integral in eq. (\[J-1med\]) $$\label{int4}
     \int^\infty_0{\mathrm{d}}|\vec{p}_E||\vec{p}_E|^4\frac{n_{qM}}{E_p(M)}
     =\int^\infty_M\frac{|\vec{p}_E|^3{\mathrm{d}}E_p}{1+e^{\frac{E_p(M)}{T}}} 
     =T^4\int^\infty_0\!{\mathrm{d}}x\frac{(x^2+2x\frac{M}{T})^{\frac{3}{2}}}{1+
     e^{x+\frac{M}{T}}}.$$ If $M=0$, then the integral can be evaluated explicitly and is found to be $\frac{7\pi^4}{120}T^4$, leading us to the result (\[limit\]). If $M\neq 0$, we obtain at once the estimate at small $T$ $$\int^\infty_0\!{\mathrm{d}}x\frac{(x^2+2x\frac{M}{T})^{\frac{3}{2}}}{1+
     e^{x+\frac{M}{T}}}
     =\frac{3\sqrt{\pi}}{4}\left(2\frac{M}{T}\right)^{\frac{3}{2}}
     e^{-\frac{M}{T}} 
     \left(1+{\cal O}(T)\right).$$ We conclude that the integral vanishes exponentially for small $T$ and does not contribute to the leading order term in eq. (\[limit\]). It is then clear that the action of the Pauli-Villars regulator $\hat{\rho}_{\Lambda\vec{p}_E}$ on the integrand, which consists in subtracting the 