# 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.

• I am in the same class. Would appreciate an update if you get it. – Lenol Dec 11 '17 at 23:57
• Sure. I'll post anything relevant I can find. – Sahand Tabatabaei Dec 11 '17 at 23:59
• I found $\mu(T)=\mu(0)+cT^2$, and that $\mu$ barely deviates from $\mu(0)$. I am not sure wether this is what he expects. – Lenol Dec 13 '17 at 1:37
• I don't know maybe? But that's just for low temperatures. He says prove $\mu=0$ for all T. I have no idea really. – Sahand Tabatabaei Dec 13 '17 at 2:14

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.

• Thank you so much for your answer. I fail to see where the specific dispersion relation $\epsilon_k^{\pm} =\pm \hbar c \mathbf k$ comes into your derivation, apart from the symmetry used in your density of states functions. Is this true for any form of (symmetric) single particle energies for the fermions ? – Sahand Tabatabaei Dec 13 '17 at 21:55
• @SahandTabatabaei The dispersion relation merely gives you the density of states and nothing more. This is true whenever you calculate averages in statistical mechanics. In this case the important property of the density of state is evenness, i.e. $D(\epsilon)=D(-\epsilon)$, which is the consequence of having symmetric energy level structure around $\epsilon=0$. So if you have a dispersion relation that is symmetric around $\epsilon=0$ you will get the same result that $\mu(T)=0$ for all $T$'s. This should not be a surprise to you. There is noting special about this specific case. – eranreches Dec 13 '17 at 22:03
• @SahandTabatabaei Also, if you know the dispersion relation you can even calculate explicitly, under several approximations, both $n_{\rm e}(T)$ and $n_{\rm h}(T)$. If you'll do this you will see that as the temperature rises, more electron-hole pairs are created. – eranreches Dec 13 '17 at 22:09
• @SahandTabatabaei Your observation about the Fermi-Dirac distribution at $\epsilon\rightarrow-\infty$ is correct. But that's exactly why we talk about holes instead of electrons in the valance band. You can easily see that $f_{\rm h}(\epsilon\rightarrow-\infty)\sim e^{\beta(\epsilon-\mu)}\ll 1$. The key here is that you can calculate all the quantities you want, but you have to treat the system as electrons in the conduction band and holes in the valance band. By doing so you can discard the infinite occupancy of electrons in the limit $\epsilon\rightarrow-\infty$ which causes problems. – eranreches Dec 14 '17 at 11:11
• @SahandTabatabaei All the physics happens between the conduction and the valance band. Energy levels that are too deep down (or too deep up above the conduction band) do not contribute to any measurable effects in reasonable temperatures. – eranreches Dec 14 '17 at 11:16

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.）

• If I to calculate the internal energy of the system, $U=\int_{-\infty}^{\infty}\epsilon f(\epsilon,\mu=0) D(\epsilon) d\epsilon$, the integral would be divergent. How can I make sense of this? – Lenol Dec 14 '17 at 4:01
• @Lenol You just don't do that this way. Instead, lets fix $U(T=0)=0$. We can do this because energy is defined up to a constant. Now imagine you excite an electron from level $-\epsilon$ into level $\epsilon$. It means that you create a hole at $-\epsilon$ with energy $-\left(-\epsilon\right)=\epsilon$ and an electron at $\epsilon$ with energy $\epsilon$. The total energy of the system has increased by $2\epsilon$. – eranreches Dec 14 '17 at 11:22
• @Lenol By the same reasoning you can construct a general expression for the energy of the system $$U=\int_{-\infty}^{0}(-\epsilon)f_{\rm h}(\epsilon)D(\epsilon){\rm d}\epsilon+\int_{0}^{\infty}\epsilon f_{\rm e}(\epsilon)D(\epsilon){\rm d}\epsilon$$ This expression converges because $$f_{\rm h}(\epsilon\rightarrow-\infty)\rightarrow 0$$ and $$f_{\rm e}(\epsilon\rightarrow\infty)\rightarrow 0$$ exponentially fast. – eranreches Dec 14 '17 at 11:23
• @Lenol Yup, eranreches's already answers your question. Or, as I mentioned above, if you really just wanted to stare at the question itself, and deem it completely as a "single electron energy level filling" problem without considering holes in solids or energy constant shifting, then the problem, strictly speaking, is not perfectly defined --- there should be always a lower bound of energy for any system. So to solve it, you must explain it properly first. – Kite.Y Dec 14 '17 at 14:42