I'm working on a problem where the following integral showed up $$I\left( t \right)=\int\limits_{0}^{t}{d\tau {{e}^{i{{\omega }_{0}}\tau }}\int\limits_{0}^{\infty }{d\omega \sqrt{f\left( \omega \right)}{{e}^{-i\omega \left( \tau -\gamma \right)}}}}, $$ where $$f\left( \omega \right)=\frac{{{\omega }^{3}}}{{{e}^{\omega }}-1}.$$
I solved it numerically without much problem. For my surprise, however, I noticed that for large values of $t$ and $\gamma$, the numerical integral converges to $$\underset{t\to \infty }{\mathop{\lim }}\,I\left( t \right)=Cf\left( {{\omega }_{0}} \right),$$ where $C$ is a constant.
This is a beautiful result and I spent the last couple of hours trying to figure out analytically why it does happen. I played with the stationary phase method, with integral representations of Dirac $\delta$, and many other things, but I had no success so far.
I'm wondering whether anyone has seen such type of result or may have any insights on where it comes from.