In Le Bellac's statistical physics book he derives the Sommerfeld expansion by a contour integral.
The idea is to expand integrals of the type $I(\beta)\equiv \int_{0}^{\infty}d\epsilon\, \frac{\varphi'(\epsilon)}{e^{\beta(\epsilon-\mu)}+1}$. Integrate by parts and expand $\varphi(\epsilon)$
$$\varphi(\epsilon)=\sum_{m=0}^{\infty}\frac{(\epsilon-\mu)^{m}}{m!}\left(\frac{d ^{m}\varphi}{d \epsilon}\Bigg |_{\epsilon=\mu}\right)$$
Make a change of variables to write $I(\beta)=\sum_{m=0}^{\infty}\frac{\beta^{-m}\varphi^{(m)(\mu)}}{m!}I_{m}$ where $I_{m}=\int_{-\infty}^{\infty}dx\, \frac{x^{m}e^{x}}{(e^{x}+1)^{2}}$.
The idea then to evaluate the $I_m$ is to consider the following integral
$$J(p)=\int_{-\infty}^{\infty}dx\, \frac{e^{ipx}}{(e^{x}+1)(e^{-x}+1)}=\sum_{m=0}^{\infty}\frac{(ip)^{m}}{m!}\int_{-\infty}^{\infty}\, dx \frac{x^{m}}{(e^{x}+1)(e^{-x}+1)}= \sum_{m=0}^{\infty}\frac{(ip)^{m}}{m!}I_{m}$$
Integrating $J(p)$ gives $J(p)=\frac{2\pi pe^{-\pi p}}{1-e^{-2\pi p}}=1-\frac{1}{6}(\pi p)^2+\ldots$ hence,
$$I(\beta)= \varphi(\mu) + \frac{\pi^{2}}{6}(k_{B}T)^{2}\varphi^{(2)}(\mu)+O(k_{B}T^{4})$$
So far so good, my problem comes when he writes the expansion of the Fermi-Dirac distribution itself as
$$\frac{1}{e^{\beta(\epsilon-\mu)}+1}\simeq \Theta_{H}(\mu-\epsilon)-\frac{\pi^{2}}{6}(k_{B}T)^{2}\delta'(\epsilon-\mu)+O(k_{B}T^{4})$$
which is meant to be valid under an integral. It appears to mimic the expansion of the generating function, at least where the signs are concerned. If we consider the definition of the Sommerfeld expansion from Wikipedia, we have a formula in the form
$$\int_{-\infty}^\infty \frac{H(\varepsilon)}{e^{\beta(\varepsilon - \mu)} + 1}\,\mathrm{d}\varepsilon = \int_{-\infty}^\mu H(\varepsilon)\,\mathrm{d}\varepsilon + \frac{\pi^2}{6}\left(\frac{1}{\beta}\right)^2H^\prime(\mu) + \mathrm{O} \left(\frac{1}{\beta\mu}\right)^4$$
If we give Le Bellac's form a try,
$$I(\beta)\simeq\int_{0}^{\infty}d\epsilon\, \varphi'(\epsilon)\left[\Theta_{H}(\mu-\epsilon)-\frac{\pi^{2}}{6}(k_{B}T)^{2}\delta'(\epsilon-\mu)+O(k_{B}T^{4}) \right] \\ =\varphi(\mu)-\frac{\pi^{2}}{6}(k_{B}T)^{2}\underbrace{\int_{0}^{\infty}d\epsilon\, \varphi'(\epsilon)\delta'(\epsilon-\mu)}_{K}$$
but integrating $K$ by parts gives
$$K=\left[\delta(\epsilon-\mu)\varphi'(\epsilon) \right]_{0}^{\infty}- \int_{0}^{\infty}d\epsilon\, \varphi''(\epsilon)\delta(\epsilon-\mu)=-\varphi''(\mu)$$
So the signs work out, and the negative sign is probably not a typo. On the other hand, the next term in the expansion would be $+\frac{7\pi^4}{360}(k_{B}T)^{4}\delta^{(3)}(\epsilon-\mu)$ which would take a negative sign when integrated.
The alternating behavior of the signs, as far as I can tell, is due to the imaginary unit that was included in the exponential for the contour integral. A very similar approach is done on p. 255 of the following notes except there is no imaginary unit in the exponential. Also, the expansion is meant to operate only on $\varphi$ instead of actually correcting the Fermi-Dirac distribution for small $T$.
My questions are as follows:
Are the corrections to the Fermi-Dirac distribution given by Le Bellac correct to all orders? They seem to work for the 0th and 1st order terms.
Assuming (1) is correct, how do we "know", or what's the reasoning, that the series expansion of the generating function would give the corrections to the Fermi-Dirac distribution? Because it seemed to me that the generating function was only to determine the $I_m$ coefficients.