Doubt about interpreting Fermi level in semiconductors Before I say anything, I know there's already a bunch of questions about the Fermi level in semiconductors on this website but I don't think my doubt in particular has been addressed before.
From what I have gathered, the Fermi level of a semiconductor is equal to the chemical potential of the electrons, i.e., it's the work required to add an electron to the system. This definition is totally fine, but I can't reconcile it with the Fermi level being in the gap between the valence and conduction bands. I know that in an intrinsic semiconductor at zero temperature the Fermi level lies exactly in the middle of the valence and conduction bands. Logically, by the definition that I have found, I would think the Fermi level would have to be the bottom energy of the conduction band, since it's the minimum energy that an electron must have to occupy a vacant state in the system. If anything I would expect the Fermi level to be that or higher if we consider energy expended in simply inserting the electron in the system. How is it possible for the Fermi level to be lower than the bottom of the conduction band? I feel like I'm misinterpreting how the Fermi level actually functions and I would appreciate any help.
Just to make this clear, I understand mathematically from the Fermi-Dirac distribution that if you have 2 energy levels in a system at zero temperature, and say one had an occupancy of 1 and the other an occupancy of 0, then the chemical potential would indeed have to lie between the two levels. I can mathematically convince myself that the Fermi level lies between the two bands, I just can't reconcile it with the definition of the Fermi level.
 A: 
From what I have gathered, the Fermi level of a semiconductor is equal to the chemical potential of the electrons, i.e., it's the work required to add an electron to the system.

We've got to be a bit careful here. The chemical potential is the change in internal energy of the system when we add one electron at constant entropy. More generally, we would have that
$$\mu = \Delta U - T\Delta S$$
Only at $T=0$ do we have that $\mu = \Delta U$ when a single particle is added. This suggests that $\mu(T\rightarrow 0) = E_c$, the energy at the bottom of the conduction band, which turns out to be true in contrast to the conventional wisdom which holds that at $T=0$, the Fermi level lies in the middle of the gap.

To understand how this misconception arose and why it persists, we need to re-examine the derivation of the Fermi-Dirac distribution function
$$f(E) = \frac{1}{e^{(E-\mu)/kT}+1}$$
which gives the probability that an energy level with energy $E$ is occupied. I will follow the derivation given in ref 1.
Adopting the canonical ensemble perspective, we fix the temperature $T$ and particle number $N$; a microstate then consists of a list $\{n\}\equiv\{n_1,n_2,\ldots\}$ of the occupation numbers of each single-particle energy level.
The partition function becomes
$$Z(T,N) = \sum_{\{n\}\in M_N} e^{-\beta E\big(\{n\}\big)} = \sum_{\{n\}} \exp\left[-\beta\sum_i n_i \epsilon_i\right]=\sum_{\{n\}\in M_N} \prod_i e^{-\beta n_i \epsilon_i}$$
where $M_N$ is the set of microstates such that $\sum_i n_i = N$, the summation over $i$ corresponds to the single-particle states, and $\epsilon_i$ is the energy of the $i^{th}$ state.  This summation cannot be easily performed for arbitrary $N$ because of the constraint $\{n\}\in M_N$. However, we can extract the distribution function as follows.
Consider in particular the $I^{th}$ state. The probability that the state is occupied is given by
$$f_I(T,N) = \sum_{\matrix{\{n\}\in M_N\\n_I=1}}e^{-\beta E\big(\{n\}\big)}/Z(T,N)$$
where the sum in the numerator is taken over all microstates $\{n\}\in M_N$ in which $n_I = 1$.  Since either a state is occupied or it is not, we may equivalently write this as $1$ minus the probability that the state $I$ is not occupied, i.e.
$$f_I(T,N) = 1-\sum_{\matrix{\{n\}\in M_N\\n_I = 0}} e^{-\beta E\big(\{n\}\big)}/Z(T,N)$$
Next, we note that each $N$-particle state with energy $E$ in which state $I$ is not occupied corresponds to a unique $N+1$-particle state with energy $E+\epsilon_I$ in which state $I$ is occupied and vice-versa.  Therefore, we may write
$$f_I(T,N) = 1- \sum_{\matrix{\{n\}\in M_{N+1} \\ n_I = 1}} e^{-\beta\bigg(E\big(\{n\}\big)-\epsilon_I\bigg)} / Z(T,N)= 1-e^{\beta \epsilon_I} f_I(T,N+1) \frac{Z(T,N+1)}{Z(T,N)}$$
Recalling that $Z(T,N) = e^{-\beta \mathcal F(T,N)}$ with $\mathcal F$ the Helmholtz potential, we observe that $$Z(T,N+1)/Z(T,N)=\exp\bigg(-\beta\big(\mathcal F(T,N+1)-\mathcal F(T,N)\big)\bigg)\equiv \exp\big(-\beta \mu(T,N)\big)$$ and finally obtain
$$f_I(T,N) = 1 - e^{\beta(\epsilon_I-\mu)}f_I(T,N+1)$$
where the chemical potential has been defined as $\mu(T,N)\equiv \mathcal F(T,N+1)-\mathcal F(T,N)$.

The Fermi-Dirac distribution is easily obtained under the assumption that $f_I(T,N)\approx f_I(T,N+1)$. The conventional wisdom is that this assumption is valid universally, because it is seemingly ridiculous to assume that the addition of another electron to a system cannot have any meaningful effect on the occupation probabilities. However, this is not true for a band insulator at $T=0$.
Let $N_0$ be the number of electrons which would completely fill the valance band and leave the conduction band empty, and let $I$ be the lowest energy state in the conduction band.  Clearly $f_I(0,N_0)= 0$, but the addition of a single additional electron yields $f_I(0,N_0+1)=1$, rendering our assumption invalid.
A correct treatment of this problem can be found in ref 2.  The result is that the low-temperature behavior of $\mu$ is given by
$$\mu= E_c -\frac{1}{2}\Delta+ kT \ln\left[\left(\frac{m_v}{m_c}\right)^{3/4} \alpha(T,V)\right] -kT N_i(T,V) \left[\alpha(T,V)+\alpha^{-1}(T,V)-2\right]$$
$$N_i(T,V) = \frac{1}{4}V \left(\frac{2kT}{\pi\hbar^2}\right)^{3/2} (m_vm_c)^{3/4} e^{-\Delta/2kT}, \qquad \alpha(T,V) = \frac{1+\sqrt{1+4N_i^2(T,V)}}{2N_i(T,V)}$$
where $E_c$ is the energy at the bottom of the conduction band, $\Delta$ is the band gap, and $m_v$ and $m_c$ are the effective masses at the top of the valence band and bottom of the conduction band, respectively.  This rather unpleasant expression yields the following behavior as $T\rightarrow 0$ (figure taken from ref 2):



*

*Ashcroft and Mermin, p.40-42


*M. R. A. Shegelski, Solid State Commun. 58, 351–354, 1986
A: A possible way to see this is as a consequence of the particle-hole symmetry.
Consider the following Hamiltonian as a simple example of a material with gapped spectrum:
$$
H = \sum_{k\sigma} \left( \varepsilon_{k,\text{val}} c^{\dagger}_{k\sigma,\text{val}} c_{k\sigma,\text{val}} + \varepsilon_{k,\text{cond}} c^{\dagger}_{k\sigma,\text{cond}} c_{k\sigma,\text{cond}} \right),
$$
where $k$ is the lattice momentum, $\sigma=\uparrow,\downarrow$ is the spin index and $\text{val.}$ ($\text{con.}$) represents the valence (conduction) band.
The particle-hole symmetry occurs when $H$ is invariant under the transformation of the fermionic operators
$$
c_{k\sigma,\text{val.}} \to c^{\dagger}_{k\sigma,\text{cond.}} 
\;\;\;\;
c_{k\sigma,\text{cond.}} \to c^{\dagger}_{k\sigma,\text{val.}},
$$
which manifestly occurs when $\varepsilon_{k,\text{val.}} = - \varepsilon_{k,\text{cond.}}$ (i.e. when the bands are symmetric with respect to the line $\varepsilon=0$, which therefore represents the middle of the gap).
In physical terms this symmetry implies that you can describe any eigenstate of the system either in terms of electrons occupying some single particle energy states, or in terms of equivalent fermions (holes) occupying states with opposite energies.
For instance, the ground state can be equivalently seen as electrons filling up the valence band, or as holes filling up the conduction band.
The same property holds for excited states as well: these states can be described by the electron Fermi function
$$
f(\varepsilon) = \frac{1}{e^{\beta(\varepsilon-\mu)}+1},
$$
which describes the probability that electrons occupy a state with energy $\varepsilon$ when the inverse temperature of the system is $\beta$.
The same state can be described in terms of the hole Fermi function $h(\varepsilon)$, which is $h(\varepsilon) = f(-\varepsilon)$ (holes are still fermions, but they populate states with opposite energies).
Finally, since holes are "lack of electrons", the relation $h(\varepsilon) = 1-f(\varepsilon)$ should hold, leading to the mathematical equation
$$
f(-\varepsilon) = 1 - f(\varepsilon) \;\;\; \to \;\;\; \mu=0,
$$
which states that the Fermi level should be exactly in the middle of the gap.
I am not sure this is the best way to answer, nevertheless this is the most natural explanation that comes to my mind.
A: The idea of the Fermi energy is probably more relevant to the theory of metals than to semiconductors, and in all honesty probably shouldn't be used outside of a somewhat narrowly defined scope in jellium-like models (treating electrons as a confined gas of non-interacting fermions.)
The particle-like excitations that can be observed in condensed matter systems and that are described by band theory are in reality only approximately analogous to a gas of free fermions.  A more apt theoretical description requires ideas from quantum field theory, the formalism of the standard model.  In short, the material in question (whether an insulator, conductor, or semi-conductor) is modeled as a special kind of vacuum (a stable or metastable arrangement of matter fields) which supports particle-like excitations particular to that vacuum (the details of which being deemed too complex to be treated beyond a small number of (small) non-linear terms intended to model small perturbations from the ground state.)  The parameters of the resulting model (effective mass, a kinetic term, interactions) usually depend on values of experimentally adjustable external fields, which easily break any symmetries that the system might possess "in vacuo" (or in a "meta-vacuum" w.r.t. 'measurement apparatuses' constructed from intrinsic/material-specific fields and operators.)
But by abuse of notation, the phrase "Fermi Energy" is often used in settings where it probably shouldn't, where it is defined on a per-model basis as the limit of the chemical potential operator for quasi-particles as temperature approaches zero, ceterus paribus (including field definitions.)  Hence, it is possible to have two Fermi Energies: one that describes the energy required to excite a quasi-particle with zero crystal-momentum, and the other that describes the corresponding energy for quasi-holes.  In general, there is a Fermi energy for every distinct type of particle-like excitation (which might not even be anything like electrons, especially localized/non-bulk theories of two-dimensional surfaces and interfaces and one-dimensional boundaries or cusps.)
A: First, as you note, in the context of semiconductors, the Fermi level is synonymous with the chemical potential, which is a meaningful quantity at all temperatures. (In contrast, the Fermi energy is something defined at absolute zero.)
Absolute zero is often a special case, so I think it's helpful to think about other situations first. Where is the Fermi level at nonzero temperatures? What happens as you take the limit as $T \to 0$? (I.e. what's $\lim_{T \to 0} E_F\left(T\right)$?) Does it end up in the middle of the gap or at its edge?
I'd argue that, in some sense, the Fermi level is undefined at zero temperature; strictly speaking, it could be anywhere in the gap --- from its lower bound to its upper bound. The math would work out the same. However, discontinuities are annoying, so I think that it makes sense for $E_F\left(0\right) \equiv \lim_{T\to0} E_F\left(T\right)$
