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.

$\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}
