How to expand $b\int_0^\infty \operatorname{sech}^2\big(b\cdot f(x)\big)\,dx$ for large $b$? Suppose $f(x)$ has a single zero in $(0,\infty)$ at $x=c$ and has a Taylor expansion about this point with some nonzero radius of convergence $0<R\leq\infty$. For concreteness, I'm working with the case $f(x)=\sqrt{x^2+a^2}-M^2$ where $0<a<M$. There are probably some other assumptions about $f(x)$ that are necessary to make this a well-formed question (continuity, smoothness, etc.), but I'm not sure what they are and hopefully my example makes clear the type of situation I'm talking about.
Playing around with the function $b\cdot \operatorname{sech}^2(bx)$, one is immediately convinced that this behaves like a Dirac delta distribution: The area under the whole curve is always $2$, and the curve consists of a single "peak" that gets narrower and taller as $b\to\infty.$ Assuming that properties of the Dirac delta distribution are valid here, it is straightforward to show that
$$\lim_{b\to\infty}\left(b\int_0^\infty \operatorname{sech}^2\big(b\cdot f(x)\big) \,dx\right)=\frac{2}{|f'(c)|}.$$
I'm trying to do slightly better, expanding this result for large $b$, rather than taking the full limit $b\to\infty$. Specifically, I'm looking for something like
$$b\int_0^\infty \operatorname{sech}^2\big(b\cdot f(x)\big)\,dx =\frac{2}{|f'(c)|}+\frac{1}{b}\bigg(\dots\bigg)+\frac{1}{b²}\bigg(\dots\bigg)+\dots.$$
How can I find and rigorously justify such an expansion?

Full disclosure, I've already found a method that gives an answer, but it requires many shaky assumptions that I can't rigorously justify. I'll give my own method below, and would be happy for someone to provide justification for the questionable steps, or present a completely different method that is better justified than this.


*

*$\operatorname{sech}^2(x)$ is well approximated by $e^{-x^2}$ near its peak, but the integral of a Gaussian is $\sqrt{\pi}$ whereas the integral of $\operatorname{sech}^2$ is $2$, so we can approximate $b\int_0^\infty \operatorname{sech}^2\big(b\cdot f(x)\big)\,dx$ by $\frac{2}{\sqrt{\pi}}b\int_0^\infty \exp\big[-\big(b\cdot f(x)\big)^2\big]\,dx$.

*We can expand $f(x)$ with its Taylor series about $x=c$ and then manipulate the integrand as follows:


$
\begin{align}
\exp\big[-\big(b\cdot f(x)\big)^2\big]&\approx \exp\big[-\big(b\cdot (a_1 (x-c)+a_2 (x-c)^2+\dots)\big)^2\big]\\\\
&= \exp\big[-\big(b\cdot a_1 (x-c)\big)^2\big]\Big[1+\alpha_1(x-c)+\alpha_2(x-c)^2+\dots\Big]
\end{align}$


*We can change the region of integration to $(-\infty,+\infty)$ because the contribution from the negative axis becomes negligible for large $b$. Each piece of the expanded integrand then becomes exactly solvable and collecting terms with like powers of $1/b$ gives the desired expansion.


Of course, each one of these steps has shaky assumptions. In (1), how do we know the Gaussian approximation is valid, especially the part where we multiply by the correction factor $2/\sqrt{\pi}$? For (2), the manipulated integrand only converges to the original integrand in the radius of convergence of the Taylor series for $f(x)$. Outside this region, the final error is "small" due to the Gaussian prefactor, but it still isn't clear that this must produce the correct result. For (3), how do we know this step doesn't spoil the expansion?
My method appears to give me reasonable results, but I don't fully trust it and at any rate, I want to understand it in a more rigorous way.
 A: *

*Since ${\rm sech}=1/\cosh$ is an even function we can replace $f$ with $|f|$. Assume for starters 


*

*that $|f|$ is continuous; 

*that $|f(x)|$ goes fast enough to $\infty$ for $x\to \infty$ such that OP's integral 
$$I(b)~:=~b\int_{\mathbb{R}_+}\! \mathrm{d}x~{\rm sech}^2(b|f(x)|), \qquad b>0,\tag{1} $$
is convergent; 

*that $x=c>0$ is the only zero for $|f|$; 

*and that $f$ is differentiable at $x=c$,
$$\lim_{x\to c}\frac{f(x)}{x-c}~=~\lim_{x\to c}\frac{f(x)-f(c)}{x-c}~=~f^{\prime}(c)~\neq~0,\tag{2} $$
with non-zero derivative.
Next make the substitution 
$$y~=~b(x-c)\qquad\Leftrightarrow \qquad x~=~c+\frac{y}{b}\tag{3}$$
in OP's integral
$$I(b)~\stackrel{(1)+(3)+(5)}{=}~\int_{\mathbb{R}}\!\mathrm{d}y~1_{[-cb,\infty[}(y)~g_b(y), \tag{4}$$ 
where the integrand 
$$ g_b(y)~:=~{\rm sech}^2(b|f(c+\frac{y}{b})|)\quad\stackrel{(2)}{\longrightarrow}\quad {\rm sech}^2(|f^{\prime}(c)|y) \quad {\rm for}\quad b~\to~\infty\tag{5}$$
$y$-pointwise. 
Then use the Lebesgue's dominated convergence theorem (i.e. assume there exists a majorant function) to argue that
$$\lim_{b\to\infty} I(b)~=~\int_{\mathbb{R}}\!\mathrm{d}y~{\rm sech}^2(|f^{\prime}(c)|y)~=~\frac{ 1}{|f'(c)|}\left[\tanh(|f^{\prime}(c)|y) \right]_{y=-\infty}^{y=\infty} ~=~\frac{2}{|f'(c)|}.\tag{6}$$

*Now we impose further conditions on $f$ to discuss possible (one-sided) power series in $1/b>0$. Let us assume 


*

*that $f:\mathbb{R}\to \mathbb{R}$ is defined on the entire real axis; 

*that it is real analytic in some interval $[c\!-\!R, c\!+\!R]$; more precisely, that the Taylor series $\sum_{n\in\mathbb{N}_0}a_n(x-c)^n$ for $f$ at $x=c$ is convergent to $f$ in the entire interval $[c\!-\!R, c\!+\!R]$; 

*and that $|f|$ grows to $\infty$ at a least linear rate for $|x|\to \infty$. 
From now on we will be somewhat sketchy. It is easy to see that OP's integral (1) outside the interval $[c\!-\!R, c\!+\!R]$ will be exponentially suppressed (in fact non-perturbative in $1/b>0$) and cannot contribute to the (one-sided) power series in $1/b>0$, cf. comment by OP. So it's enough to consider the integral
$$ I^R(b)~:=~b\int_{\mathbb{R}_+}\! \mathrm{d}x~1_{[c-R, c+R]}(x)~{\rm sech}^2(bf(x))~=~\int_{\mathbb{R}}\!\mathrm{d}y~1_{[-bR,bR]}(y)~g_b^R(y) ,\tag{7} $$
where the integrand 
$$ g_b^R(y)~:=~{\rm sech}^2(yh(c+\frac{y}{b}))
~=~\sum_{n\in\mathbb{N}_0}g_n^R(y)b^{-n},\tag{8} $$
and we have defined
$$ h(x)~:= ~\frac{f(x)}{(x-c)}.\tag{9} $$
There will be a corresponding (one-sided) Taylor series
$$ \sum_{n\in\mathbb{N}_0}I_n^R b^{-n}, \qquad 1/b>0, \tag{10}$$
such that 
$$ I^R(b) -  \sum_{n=0}^NI_n^R b^{-n}~=~{o}(b^{-N}) .\tag{11}$$
