Chemical potential for ideal fermi gas with single particle energies $\epsilon_k^{\pm}=\pm \hbar c |\mathbf k|$ I came across the following question in some example problems for my graduate statistical mechanics course final exam:

Consider an ideal gas of fermions of density $n$ in three dimensions with the single-particle eigenstate energies given by $\epsilon_k^{\pm}=\pm \hbar c |\mathbf k|$.  Assume the chemical potential $\mu = 0$ at $T=0$. Prove that if the fermion density is constant, $\mu(T)=0$ for all $T$.

My first approach was to use the Fermi-Dirac distribution (using energy units for the temperature):
$$N=\sum_k\langle n_k\rangle = \sum_k \frac{1}{e^{\frac {\epsilon_k-\mu}T}+1}=\sum_k \frac{1}{e^{\frac {\hbar  c |\mathbf k|-\mu}T}+1}$$
Converting this into an integral in momentum space:
$$N=V\int_{\Bbb R^3}\frac {d^3k}{(2 \pi)^3}\frac{1}{e^{\frac {\hbar  c |\mathbf k|-\mu}T}+1}=\frac {V}{2 \pi^2}\int_0^\infty dk \ \frac{k^2}{z^{-1}e^{\frac {\hbar  c k}T}+1}$$
Where $z=e^{\frac {\mu}{T}}$ is the fugacity. Using change of variables, we get:
$$N=\frac {V}{2 \pi^2} \left(\frac{T}{\hbar c}\right)^3\int_0^\infty dx \ \frac{x^2}{z^{-1}e^{x}+1}$$
Which in terms of the Fermi-Dirac functions, is :
$$n=\frac {1}{ \pi^2} \left(\frac{T}{\hbar c}\right)^3f_3(z)$$
Where $n$ is the density. At low temperatures, using the Sommerfeld expansion $f_3\left(e^{\frac{\mu}{T}}\right) \simeq \frac 16\left(\frac {\mu}{T}\right)^3\left(1+\pi^2 \left(\frac T{\mu}\right)^2\right)$, we get:
$$n\simeq\frac {1}{ \pi^2} \left(\frac{T}{\hbar c}\right)^3\frac 16\left(\frac {\mu}{T}\right)^3\left(1+\pi^2 \left(\frac T{\mu}\right)^2\right)\simeq \frac {V}{ 6\pi^2} \left(\frac{\mu}{\hbar c}\right)^3 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ T\rightarrow0^+ $$ 
Which gives : 
$$\mu \simeq \hbar c\left(6 \pi^2 n\right)^{\frac{1}{3}}  \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ T\rightarrow0^+ $$
I can't see how we can assume the chemical potential to be zero at zero temperature. Clearly, the  above solution shows that that is impossible. How could the gas have zero fermi energy? I get the same result when calculating the grand canonical partition function from scratch and using the grand potential. A zero chemical potential doesn't make intuitive sense either, because it implies that adding or subtracting particles from the gas would have no energy cost whatsoever; which can't be true, because when we add a fermion it can't share an energy level with another particle (Pauli principle), and thus has to go to a higher energy level, which means that there is a nonzero energy cost.

Another problem is that I can't understand what the question means by the $\pm$ in the $\epsilon_k^{\pm}=\pm \hbar  c |\mathbf k|$ relation. Does it mean that each eigenstate of momentum corresponds to two energy eigenstates, one negative and one positive? For example, consider the - case $\epsilon_k=- \hbar  c |\mathbf k|$. The same reasoning as above gives:
$$n=\frac {1}{2 \pi^2} \left(\frac{T}{\hbar c}\right)^3\int_0^\infty dx \ \frac{x^2}{z^{-1}e^{-x}+1} \ \rightarrow +\infty$$
And this integral clearly diverges for all values of $z$.
I feel like I'm misinterpreting the question. Any help would be appreciated.
 A: Before delving into mathematics, lets try to get some intuition about the energy function you're dealing with. You are given
$$\epsilon_{\boldsymbol{k}}^{\pm}=\pm\hbar c\left|\boldsymbol{k}\right|$$
$$\epsilon_{\rm F}=\mu\left(T=0\right)=0$$
This means that at $T=0$ all the negative energy states are occupied and all the positive energy states are empty. This is typical for systems in solid-state physics. You basically have here a Band Structure: the negative energies form what is known as the Valence Band, and the positive energies from what is known as the Conduction Band. The valence band is fully occupied with electrons, which makes it hard to deal with. Therefore, physicists create the concept of Holes, or lack of electrons. While the valence band is full of electrons, it is equivalent to say that it is empty of holes - easy to treat. In summary, holes live in the valence band and electrons live in the conduction band.
Now comes the key part. When you excite an electron from the valence band to the conduction band you are essentially creating a pair of electron-hole: an hole (electron vacancy) appears in the valence band and an electron appears in the conduction band. It is very important now to note that because of this symmetry the number of electrons is equal to the number of holes. This is true only because the Fermi level is located exactly in the middle between the two bands.
Now we understand the intuition behind this problem. The next step is to calculate the number of electrons and the number of holes, require the two to be the same - and get $\mu=\mu\left(T\right)$. As you said, the probability of electron occupation at an energy $\epsilon$ is given by the Fermi-Dirac Distribution function
$$f_{\rm e}\left(\epsilon\right)=\frac{1}{e^{\beta\left(\epsilon-\mu\right)}+1}$$
with $\beta\equiv \frac{1}{K_{\rm B}T}$. On the other hand, the probability of occupation of hole in the valence band is equal to the probability of lack of electron
$$f_{\rm h}\left(\epsilon\right)=1-f_{\rm e}\left(\epsilon\right)=1-\frac{1}{e^{\beta\left(\epsilon-\mu\right)}+1}=\frac{1}{e^{-\beta\left(\epsilon-\mu\right)}+1}$$
Thus the number of electrons and holes respectively is given by
$$n_{\rm e}\left(T\right)=\int_{0}^{\infty}{\rm d}\epsilon D\left(\epsilon\right)f_{\rm e}\left(\epsilon\right)$$
$$n_{\rm h}\left(T\right)=\int_{-\infty}^{0}{\rm d}\epsilon D\left(\epsilon\right)f_{\rm h}\left(\epsilon\right)=\int_{0}^{\infty}{\rm d}\epsilon D\left(-\epsilon\right)f_{\rm h}\left(-\epsilon\right)$$
at all temperatures. Here $D\left(\epsilon\right)$ is the density of states, and you can argue that from symmetry $D\left(\epsilon\right)=D\left(-\epsilon\right)$. Plugging in the distributions, we get
$$n_{\rm e}\left(T\right)=\int_{0}^{\infty}{\rm d}\epsilon D\left(\epsilon\right)\frac{1}{e^{\beta\left(\epsilon-\mu\right)}+1}$$
$$n_{\rm h}\left(T\right)=\int_{0}^{\infty}{\rm d}\epsilon D\left(\epsilon\right)\frac{1}{e^{\beta\left(\epsilon+\mu\right)}+1}$$
You can now immediately see that if $n_{\rm e}\left(T\right)=n_{\rm h}\left(T\right)$ is true for all temperatures, then you must have
$$\mu\left(T\right)=0$$
for all temperatures, as wanted.
A: @eranreches's answer is complete and good enough. Just wanna add some tiny comments about the solution in your original question: 
Firstly, as you already noticed, there are two branches while in your calculation you only used one of them: that means you are trying to fill electrons into only positive energy states, then any finite number of filling would result into a positive chemical potential $\mu$ at $T=0$, which is just the highest energy level reached after you fill all electrons in.
Secondly, if you also take the negative energy into account, but still try to use a standard Fermi-Dirac integral for calculation, then you would notice that there are infinite number of negative energy levels can be filled -- no lower bound for energy levels, which results into a disaster since you don't even know where to start filling electrons in. And that's part of the reason @eranreches used density of states $D(\epsilon)$ in his calculation, which could avoid the problem in an explicit calculation. While in the real world, for example, in solid state systems, any energy band of course is restricted in a finite energy window, and also the momentum $k$ would take a finite value as well; then, more specifically, you may have $D(|\epsilon| > \epsilon_0) = 0$.
Thirdly, the linear dispersion is quite common in solid state, related to many interesting topics: e.g. Luttinger liquid, Dirac cones.... It is usually an approximation for only a finite number of states labled by $k$ in the full Brilioun zone, and is mostly seen for the band crossing (or touching) case: when two band cross each other, the region around the crossing point would hold a linear dispersion $\epsilon \sim |k|$. And the meaning of $\mu(T=0)=0$ in this case is just that electrons fill all of the states below this crossing point.
In general there would be other dispersions for band touching, e.g. $\epsilon\sim k^2$. (You could find something related here: arXiv 1603.03093.）
