V(M_f,T,\mu )&\!\! =\!\!& -\frac{1}{2}\int_0^{M_f}\!\!\sum_{f=u,d,s}h_f{\mathrm{d}}M_f-\frac{N_c}{4\pi^2}\int_0^{M_f}\!\!\sum_{f=u,d,s}M_fJ_0(M_f^2,T,\mu ){\mathrm{d}}M_f + C(T,\mu )\nonumber \\
     &\!\! = \!\!&\frac{1}{16}\left.\left(4Gh_f^2+\kappa h_uh_dh_s 
     +\frac{3g_1}{2}\left(h_f^2\right)^2+3g_2h_f^4\right)\right|_0^{M_f} 
     \!\! +\frac{N_c}{8\pi^2}\!\sum_{f=u,d,s}\!\! J_{-1}(M_f^2,T,\mu ) 
     + C(T,\mu ),\end{aligned}$$ where the function $C(T,\mu )$ does not depend on $M$, and therefore cannot be determined from the gap equation; it will be found from other arguments in the end of this section, but obviously $C(0,0)=0$. The integral $J_{-1}(M^2,T,\mu )$ is the immediate generalization of the $T=\mu =0$ case (\[J\_1\]) $$\label{J-1}
     J_{-1}(M^2,T,\mu )=-\int_0^{M^2}\!\!\! J_0(M^2,T,\mu ){\mathrm{d}}M^2
     =J_{-1}(M^2)+J^{\mathrm{med}}_{-1}(M^2,T,\mu ),$$ where the medium contribution to $J_{-1}(M^2,T,\mu )$ is $$\begin{aligned}
     J_{-1}^{\mathrm{med}}(M^2,T,\mu )&\!\! =\!\!& 4\int_0^{M^2}\!\!\!{\mathrm{d}}M^2 
     \int^\infty_0\!{\mathrm{d}}|\vec{p}_E||\vec{p}_E|^2 \hat{\rho}_{\Lambda}
     \frac{n_{q}+n_{\overline{q}}}{E_p} =  
     4\int^\infty_0\!{\mathrm{d}}|\vec{p}_E||\vec{p}_E|^2\hat{\rho}_{\Lambda\vec{p}_E}
     \int^{M^2}_0\!\!\!{\mathrm{d}}M^2\,\frac{n_{q}+n_{\overline{q}}}{E_p}
     \nonumber\\
     &\!\! =\!\!&8\int^\infty_0\!{\mathrm{d}}|\vec{p}_E||\vec{p}_E|^2
     \hat{\rho}_{\Lambda\vec{p}_E}\left(2\left(E_p(M)-E_p(0)\right) 
     +T\ln\frac{n_{q M}n_{\overline{q} M}}{n_{q 0}n_{\overline{q} 0}}\right)
     \nonumber\\
     &\!\! =\!\!&8T\int^\infty_0\!{\mathrm{d}}|\vec{p}_E||\vec{p}_E|^2
     \hat{\rho}_{\Lambda\vec{p}_E}\ln\frac{\left(1+e^{-\frac{E_p(0)-\mu}{T}}\right)
     \left(1+e^{-\frac{E_p(0)+\mu}{T}}\right)}{\left(1+e^{-\frac{E_p(M)-\mu}{T}}
     \right)\left(1+e^{-\frac{E_p(M)+\mu}{T}}\right)}\end{aligned}$$ with $n_{q0}, n_{\overline{q}0}$ and $n_{qM}, n_{\overline{q}M}$ referring to the occupation numbers for massless and massive particles correspondingly, $E_p(M)
=\sqrt{M^2+\vec{p}^{\, 2}_E}$, $E_p(0)=|\vec{p}_E|$. In spite of the fact that the integral $J_{-1}^{\mathrm{med}}(M^2,T,\mu )$ is convergent, we still keep the regularization $\hat{\rho}_\Lambda$ to be consistent. Note that the action of the operator $\hat{\rho}_\Lambda$ (see eq. (\[reg-2\])) on any smooth function, depending on $M^2$ through the energy, $f(E_p(M))$, can also be expressed in terms of momentum as $\hat{\rho}_\Lambda f(E_p)=
\hat{\rho}_{\Lambda\vec{p}_E}f(E_p)$, where $$\hat{\rho}_{\Lambda\vec{p}_E}
     =1-\left(1-\Lambda^2 \frac{\partial}{\partial \vec{p}_E^{\, 2}}\right)
     \exp\left(\Lambda^2\frac{\partial}{\partial\vec{p}_E^{\, 2}}\right).$$

Then, noting that, for instance, 