Approximation of Fermi-Dirac integral $\int \text{d}x \frac{f(x,\beta)}{1+e^{\beta x}}$

I don't know if this will be of use, but ...

If we assume that $f$ belongs to the space of test functions that are smooth and of compact support, then certainly we have

$$f(x)-f(-x)=2f'(0)x+O(x^3)$$

Hence,

$$\begin{align} \int_0^\infty \frac{f(x)-f(-x)}{1+e^{\beta x}}\,dx&=2f'(0)\int_0^\infty \frac{x}{1+e^{\beta x}}\,dx+O\left(\int_0^\infty \frac{x^3}{1+e^{\beta x}}\,dx\right)\\\\ &=\frac{2f'(0)}{\beta^2}\int_0^\infty \frac{x}{1+e^x}\,dx+O\left(\frac{1}{\beta^4}\int_0^\infty \frac{x^3}{1+e^{\beta x}}\,dx\right)\\\\ &=\frac{\pi^2}{6\beta^2}f'(0)+O\left(\frac{7\pi^4}{120\beta^4}\right)\\\\ &=\frac{\pi^2}{6\beta^2}f'(0)+O\left(\frac{1}{\beta^4}\right) \end{align}$$

Therefore, we can write

$$\int_{-\infty}^\infty \frac{f(x)}{1+e^{\beta x}}\,dx=\int_{-\infty}^0 f(x)\,dx+\frac{\pi^2}{6\beta^2}f'(0)+O\left(\frac{1}{\beta^4}\right)$$

Now note that $\frac{1}{1+e^{\beta x}}\sim \Theta(-x)-\frac{\pi^2}{6\beta^2}\delta'(x) +O(\beta^{-4})$ in the sense that

$$\langle f, \frac{1}{1+e^{\beta x}}\rangle =\langle f,\Theta(-x)-\frac{\pi^2}{6\beta^2}\delta'(x) )\rangle +O(\beta^{-4})$$

To justify the Sommerfeld expansion let assume that $f$ is a suitable testfunction (something like a Lorentzian should be good for now: integrable at $0$ and a sufficent decay at infinity). Furthermore let us denote the Fermi distribution by $n_{\beta}(x)$.

We want to explore the integral

$$ I(\beta)=\int_{\mathbb{R}} n_{\beta}(x)f(x)dx $$

in the limit of $\beta\gg1$.

Now, let us rescale $x\rightarrow \beta x$ and split the range of integration at the origin

$$ \beta I(\beta)=I_{+}(\beta)+I_{-}(\beta)=\int_0^{\infty}\frac{f(x/\beta)}{e^x+1}dx+\int_{-\infty}^0\frac{f(x/\beta)}{e^x+1}dx $$

Since $e^{x}<1$ on $\mathbb{R}_-$ we find

$$ I_-(\beta)=\int_{-\infty}^0f(x/\beta)dx+\int_{-\infty}^0dx(f(0)+\frac{x}{\beta}f'(0)+O{(\beta^{-2})})\sum_{n\geq 1}(-1)^ne^{n x}=\\\int_{-\infty}^0f(x/\beta)dx-\log(2)f(0)+\frac{\pi^2}{12\beta}f'(0)+O(\beta^{-2}) $$

Since on $\mathbb{R}_+$ the integral converges pretty fast we can expand $f$ directly and find

$$ I_+(\beta)=\int_0^{\infty}\frac{f(0)+\frac{x}{\beta}f'(0)+O{(\beta^{-2})}}{e^x+1}dx=\log(2)f(0)+\frac{\pi^2}{12\beta}f'(0)+O(\beta^{-2}) $$

which yields (after undoing the rescaling from the beginning)

$$ I(\beta)=\int_{-\infty}^{0}f(x)dx+\frac{\pi^2}{6\beta^2}f'(0)+O(\beta^{-3}) $$

which is in the sense of distributions can indeed be interpreted as

$$ n_{\beta}(x)\sim \theta(-x)-\frac{\pi^2}{6\beta^2}\delta'(x)+O(\beta^{-3})\quad \text{as}\,\, \beta >>1 $$

that $O(\beta^{-3})$ can be replaced by an $O(\beta^{-4})$ estimate is straightforward to include

---

Remark: If we allow a temperature dependence for $f$, so $f_{\beta}(x)$ Sommerfeld should be valid as long as we can assure that a Taylor expansion around the origin (in $x$) is valid, so $f_\beta(x/\beta)\sim f_{\beta}(0)+\frac{x}{\beta}f_{\beta}'(0)...$ where $f_{\beta}^{(n)}(0)$ are at least as such that $f_{\beta}^{(n)}(0)/\beta^n \rightarrow 0$ as $\beta\rightarrow \infty$