Electron density in Semi Conductors I want to determine the relationship between the carrier density $N$ and the chemical potential $\mu_e$ in the conduction band of an intrinsic semi-conductor (GaN for example).
The band structure of the conduction band is assumed parabolic with effective mass $m_e^*=0.2m_e$. Using a $3D$ free electron carrier model the relationship is then given by : 
\begin{equation}
N = \int_{E_c}^{\infty}g(E)f(E)dE 
~~~~\mbox{with}~~~~
 f(E) = \frac{1}{1+e^{\frac{E-\mu_e}{k_BT}}}
~~~~\mbox{And}~~~~
g(E)=\frac{m_e^*\sqrt 2}{\pi^2 \hbar^3}\sqrt E
\end{equation}
Since the semi conductor is meant to work at very high carrier concentrations ($N\approx 10^{19} cm^{-3}$), we won't take advantage of the approximation $E_c-\mu_e\gg k_BT$. Using Matlab to compute these integrals, I obtain the following curves for the distributions $f(E)g(E)$ and $N(\mu_e)$ at $T=300K$:


We can clearly see that $N=2.10^{19}cm^{-3}$ is achieved with $\mu_e-E_c=11.87meV$ confirming that the approximation $E_c-\mu_e \gg k_B T$ can not be used. However, in a research paper, I have found that using the exact same values they end up with $\mu_e=130meV$. 
Paper is : 
Modal Gain in a Semiconductor Nanowire LaserWith Anisotropic Bandstructure
A. V. Maslov, Member, IEEE, and C. Z. Ning, Senior Member, IEEE
I really do not see where the discrepancy could come from after spending an afternoon on it.
 A: The carrier concentration of conduction band of a intrinsic semiconductor at a temperature $T$ is
$$N=\frac{1}{2\pi^2}\left(\frac{ 2m_e^* }{ \hbar^2}\right)^{3/2}\int_{E_c}^{\infty}\sqrt{E-E_c}f_{FD}(E,T)dE$$
where $m_e^*$ is electron effective mass and $f_{FD}$ is the Fermi-Dirac distribution function. Take $\mu$ to be the Fermi-Energy. We usually use the Boltzman approximation to evaluate the integral by takin $(E_c-\mu)/k_BT>>1$. But, if you do not use the approximation, you will get
$$N=-\frac{1}{2\pi^2}\left(\frac{ 2m_e^* k_BT}{ \hbar^2}\right)^{3/2}\Gamma\left(\frac32\right)\text{Li}_{3/2}\left(-e^{-(E_c-E_F)/k_BT}\right)$$
where $\text{Li}_{3/2}$ is the Polylogarithm. Since $\Gamma\left(\frac32\right)=\frac{\sqrt{\pi}}2$, we have
$$N=-2\left(\frac{ m_c^* k_BT}{2\pi \hbar^2}\right)^{3/2}\text{Li}_{3/2}\left(-e^{-(E_c-E_F)/k_BT}\right)$$
Then,
$$\lim_{x\to \infty}\frac{\text{Li}_{3/2}\left(-e^{-x}\right)}{e^{-x}}=\lim_{x\to 0+}\frac{\text{Li}_{3/2}\left(-x\right)}{x}= -1$$
Therefore $-\text{Li}_{3/2}\left(-e^{-x}\right)\sim e^{-x}$ when $x\to\infty$, which gives, $n_i=2\left(\frac{ m_c^* k_BT}{2\pi \hbar^2}\right)^{3/2}e^{-(E_c-E_F)/k_BT}$ when $(E_c-E_F)/k_BT>>1$ which is the same result we got by the Boltzman approximation.
But, using the Maclaurin series expansion of $\text{Li}_{3/2}$, when $x\to\infty$ we have
$$-\text{Li}_{3/2}\left(-e^{-x}\right)=e^{-x}-0.35355e^{-2x}+0.19245e^{-3x}-0.125e^{-4x}+\mathcal{O}(e^{-5x})$$
Hence, it is clear that
$$N=2\left(\frac{ m_c^* k_BT}{2\pi \hbar^2}\right)^{3/2}e^{-x}\left(1-0.35355e^{-x}+0.19245e^{-2x}-0.125e^{-3x}+\mathcal{O}(e^{-4x})\right)$$
where $x=(E_c-E_F)/k_BT$.

Here is the graph of $y=-a^\frac32\text{Li}_{3/2}\left(-e^{-x/a}\right)$ for $a=1,2,3$.
A: The error was due to the density of states function that was starting that should have been $\approx \sqrt{E-Ec}$ instead of $\sqrt{E}$. Indeed, at the bottom of the conduction band $E_c$, the density should be equal to $0$ hence the graphs for $f(E)g(E)$ should more look like the standard bell-like shape : 

