Approximation of Fermi-Dirac integral $\int \text{d}x \frac{f(x,\beta)}{1+e^{\beta x}}$ In physics is quite common to find integrals of the type
\begin{align}
I(\beta) = \int_{-\infty}^{\infty}\text{d}x \frac{f(x)}{1+e^{\beta x}} \tag{1}
\end{align}
where $f(x)$ is some quantity we want to average over the Fermi-Dirac distribution $n(x) = \left(1+e^{\beta x}\right)^{-1}$, and $\beta >0$ is a positive real parameter representing the the 'inverse of the temperature'. Since many times physicists are only interested in the 'low temperature' regime $\beta\gg 1$, it is common to consider the following Sommerfeld approximation:
\begin{align}
I(\beta)\underset{\beta \gg 1}{=} \int_{-\infty}^{0}\text{d} x~f(x)+\frac{\pi^2}{6\beta^2} f'(0)+O(\beta^{-4}). \tag{2}
\end{align}
Which sometimes appears in a mnemonic fashion as an expansion for $n(x,\beta)$ itself,
\begin{align}
\frac{1}{1+e^{\beta x}}\underset{\beta\gg1}{=} \theta(-x)-\frac{\pi^2} {6\beta^2}\delta'(x)+O(\beta^{-4}) \tag{3}
\end{align}
My question is
Can we understand Eq. (3) in a rigorous fashion? E.g. as the expansion of a distribution. If yes, can we understand the convergence of Eq.(2) as the condition for a series expansion under integral sign? E.g. dominated convergence.
My motivation
I need to study an integral similar to Eq. (1) but with the crucial difference that f(x) is also a function of the parameter $\beta$, and wanted to make a similar expansion as in Eq. (2). My idea was to understand Eq. (2) as a series expansion Eq.(3) under integrals sign, and also expand my $f$ in $\beta$. To check if this is safe, I wanted to use dominated convergence. But I am not confident this makes sense.
References for rigorous discussions of Eq. (2) or Eq. (3) are welcome.
 A: 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})$$
A: 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$
A: $\newcommand{\bbx}[1]{\,\bbox[15px,border:1px groove navy]{\displaystyle{#1}}\,}
 \newcommand{\braces}[1]{\left\lbrace\,{#1}\,\right\rbrace}
 \newcommand{\bracks}[1]{\left\lbrack\,{#1}\,\right\rbrack}
 \newcommand{\dd}{\mathrm{d}}
 \newcommand{\ds}[1]{\displaystyle{#1}}
 \newcommand{\expo}[1]{\,\mathrm{e}^{#1}\,}
 \newcommand{\ic}{\mathrm{i}}
 \newcommand{\mc}[1]{\mathcal{#1}}
 \newcommand{\mrm}[1]{\mathrm{#1}}
 \newcommand{\pars}[1]{\left(\,{#1}\,\right)}
 \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}}
 \newcommand{\root}[2][]{\,\sqrt[#1]{\,{#2}\,}\,}
 \newcommand{\totald}[3][]{\frac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}}
 \newcommand{\verts}[1]{\left\vert\,{#1}\,\right\vert}$

Note that
  $\ds{{1 \over \expo{\beta x} + 1} =
\left\{\begin{array}{lcl}
\ds{\Theta\pars{-x} + {\mrm{sgn}\pars{x} \over \expo{\beta\,\verts{x}} + 1}}
& \mbox{if} & \ds{x \not= 0}
\\[2mm]
\ds{1 \over 2} & \mbox{if} & \ds{x = 0} 
\end{array}\right.}$

  $\ds{\Theta}$ is the Heaviside Step Function.

Then,
\begin{align}
\left.\vphantom{\Large A}\mrm{I}\pars{\beta}\right\vert_{\ \beta\ >\ 0} & = \int_{-\infty}^{\infty}{\mrm{f}\pars{x} \over \expo{\beta x} + 1}\,\dd x =
\int_{-\infty}^{0}\mrm{f}\pars{x}\,\dd x +
\int_{-\infty}^{\infty}{\mrm{sgn}\pars{x} \over
\expo{\beta\,\verts{x}} + 1}\,\mrm{f}\pars{x}\,\dd x
\\[5mm] & =
\int_{-\infty}^{0}\mrm{f}\pars{x}\,\dd x + \int_{0}^{\infty}
{\mrm{f}\pars{x} - \mrm{f}\pars{-x} \over
\expo{\beta x} + 1}\,\dd x
\\[5mm] & =
\int_{-\infty}^{0}\mrm{f}\pars{x}\,\dd x +
2\sum_{n = 0}^{\infty}
{\mrm{f}^{\pars{2n + 1}}\pars{0} \over \pars{2n + 1}!}\int_{0}^{\infty}
{x^{2n + 1} \over \expo{\beta x} + 1}\,\dd x
\\[5mm] & =
\int_{-\infty}^{0}\mrm{f}\pars{x}\,\dd x +
2\sum_{n = 0}^{\infty}
{\mrm{f}^{\pars{2n + 1}}\pars{0} \over \pars{2n + 1}!}
\,\beta^{-\pars{2n + 2}}\int_{0}^{\infty}{x^{2n + 1} \over \expo{x} + 1}\,\dd x
\\[5mm] & =
\bbx{\int_{-\infty}^{0}\mrm{f}\pars{x}\,\dd x + \sum_{n = 0}^{\infty}
\bracks{2\,\mrm{f}^{\pars{2n + 1}}\pars{0}
\pars{1 - 2^{-2n - 1}}\zeta\pars{2n + 2}}\,\beta^{-\pars{2n + 2}}}
\end{align}
