chemical potential, $\mu$, through the substitutions [@Kapusta] $$\label{subst}
      \int\!{\mathrm{d}}p_{0E}\to 2\pi T\!\!\!\sum_{n=-\infty}^\infty,
      \quad p_{0E}\to\omega_n-i\mu,\quad\omega_n= \pi T (2n+1).$$ Inserting (\[subst\]) into $J_0$ we obtain $$\begin{aligned}
    J_0(M^2)&\!\!\to\!\! &J_0(M^2,T,\mu )=16\pi^2 T\int
    \frac{{\mathrm{d}}^3\vec{p}_E}{(2\pi )^3}\hat{\rho}_\Lambda
    \sum_{n=-\infty}^{+\infty}\frac{1}{C_n-2i\mu\pi T(2n+1)}
    \nonumber\\
    C_n&\!\! =\!\!&E_p^2+\pi^2 T^2 \left(2n+1\right)^2 -\mu^2, \quad
    E_p = \sqrt{M^2+\vec{p}^{\, 2}_E}.\end{aligned}$$ The sum over $n$ is evaluated to give $$\sum_{n=-\infty}^{+\infty}\frac{1}{C_n-i2\mu\omega_n}
    =\frac{1}{(2\pi T)^2}\sum^{+\infty}_{n=-\infty}\frac{1}{(n+a)(n+b)}
    =\frac{1}{4TE_p}\left(\tanh\frac{\mu +E_p}{2T}-\tanh \frac{\mu 
    -E_p}{2T}\right)$$ where we use the abbreviations $a=\frac{1}{2}+\frac{i}{2\pi T}\left(E_p-\mu
\right)$ and $b=\frac{1}{2}-\frac{i}{2\pi T}\left(E_p+\mu\right)$. Thus we have $$\label{J0t}
     J_0(M^2,T,\mu )=4\int_0^\infty{\mathrm{d}}|\vec{p}_E| |\vec{p}_E|^2
     \hat{\rho}_\Lambda\frac{1}{E_p}\left(1-n_{q}-n_{\overline{q}}\right)
     =J_0(M^2)- 4\int_0^\infty{\mathrm{d}}|\vec{p}_E| |\vec{p}_E|^2
     \hat{\rho}_\Lambda\frac{n_{q}+n_{\overline{q}}}{E_p}$$ where the quark, anti-quark occupation numbers are given by $$\label{n-occup}
     n_{q}=\frac{1}{1+e^{\frac{E_p-\mu}{T}}}, \quad 
     n_{\overline{q}}=\frac{1}{1+e^{\frac{E_p+\mu}{T}}}.$$ Notice that $J_0(M^2,0,0)=J_0(M^2)$, therefore the vacuum piece is well isolated from the matter part. The remaining integral containing the quark number occupation densities $n_q, n_{\bar q}$ is strictly finite, the $\Lambda$ dependent terms being a remnant of the Pauli–Villars regularization scheme. At small $T$ and $M\neq 0$ we have $$\begin{aligned}
\label{J0-asymp}
     \int_0^\infty\!{\mathrm{d}}|\vec{p}_E| |\vec{p}_E|^2\frac{n_{q}}{E_p}
     &\!\!\simeq\!\!&
     \int_0^\infty\!{\mathrm{d}}|\vec{p}_E| \frac{|\vec{p}_E|^2}{E_p}e^{-\frac{E_p}{T}}
     =T\sqrt{2TM}e^{-\frac{M}{T}}\int_0^\infty\!{\mathrm{d}}x\sqrt{x}e^{-x}\sqrt{1
     +\frac{xT}{2M}} \nonumber \\
     &\!\! =\!\!&\sqrt{\frac{\pi M}{2}}T^{\frac{3}{2}}e^{-\frac{M}{T}}
     \left(1+\frac{3T}{8M}+\ldots \right).\end{aligned}$$ The special feature of this integral is that it vanishes exponentially with $T\to 0$.

Now we are ready to evaluate the thermodynamical potential. Indeed, the gap-equations at finite $T$ and $\mu$ are $$\label{gap-t}
            h_f + \frac{N_c}{2\pi^2} M_f J_0(M_f^2,T,\mu )=0.$$ Consequently, $$\begin{aligned}
\label{effpot-t}
     