General formula for the variation of the chemical potential with temperature For small temperatures $T$, such that $k_\mathrm{B}T\ll \mu(T=0)\equiv \mu(0)$, the variation of chemical potential with temperature is given by $$\mu(T)=\mu(0)\left[1-\frac{\pi^2}{12}\left(\frac{k_\mathrm{B}T}{\mu(0)}\right)^2\right].\label{eqn:1}\tag{1}$$ Is there an exact analytical solution for the chemical potential at any temperature $T$ in a way that $\eqref{eqn:1}$ and its high temperature limit $(k_\mathrm{B}T\gg\mu(0))$ can be deduced from that? I am interested in such a situation because one often has to deal with Fermi gases at very high temperatures.
 A: The short answer is that there is no explicit formula for $\mu(T)$ valid at any $T$. In Peter Young's notes you can see his attempt to connect the low-$T$ and high-$T$ limits, and the basic equations can be found in most books on condensed matter physics. I'll just give the starting point here. One can equate the standard expression for $N$ at $T=0$ with the expression at higher $T$, given the 3D density of states for free electrons and the known occupation number formula at chemical potential $\mu(T)$. This gives
$$
\frac{2}{3}\epsilon_F^{3/2} = \int_0^\infty d\epsilon \, 
\frac{\epsilon^{1/2}}{\exp[\beta(\epsilon-\mu)]+1}
$$
where $\epsilon_F=\mu(0)$ and $\beta=1/k_BT$. The chemical potential as a function of $T$ may be obtained by solving this equation. 
As those notes point out, it is convenient to express it in dimensionless form by defining 
$x=\epsilon/\epsilon_F$,
$\tilde{\mu}=\mu/\epsilon_F$, and $\tilde{\beta}=\beta\epsilon_F$ or equivalently $\tilde{T}=k_BT/\epsilon_F$:
$$
\frac{2}{3} = \int_0^\infty dx \, 
\frac{x^{1/2}}{\exp[\tilde{\beta}(x-\tilde{\mu})]+1}
= \int_0^\infty dx \, 
\frac{x^{1/2}}{\exp[(x-\tilde{\mu})/\tilde{T}]+1} .
$$
For $\tilde{T}=0$ the Fermi function is zero for $x>\tilde{\mu}$ and the equation is satisfied for $\tilde{\mu}=1$. You are after the solution $\tilde{\mu}(\tilde{T})$ for nonzero $\tilde{T}$.
Evidently, this solution cannot be obtained analytically. It can be done numerically (see Fig 1 of the above-cited notes) and, of course, one can get low-$T$ and high-$T$ asymptotic forms for $\mu(T)$. The low-$T$ form you have already. The high-temperature limit follows from the fact that $\mu$ becomes large and negative at high $T$, and so the $+1$ in the denominator of the integrand may be dropped in this limit. Then the $\mu$-dependence may be taken out of the integral, and the integral itself can be done; the resulting equation can be rearranged to give $\mu$. See this answer for details.
