Asymptotic power series expansion of $\int_0^\infty\frac{x^\nu J_\nu(x\alpha)}{e^x-1}{\rm d}x$ around $\alpha=1$ and $\alpha<1$ What is the power series expansion $f(\alpha)$ as $\alpha\to 1$ and also in the limit $\alpha\ll 1$, where
$$f(\alpha)=\int_0^\infty\frac{x^\nu J_\nu(x\alpha)}{e^x-1}{\rm d}x$$
At $\alpha=1$ the integrand is analytic so we should have a convergent power series. The power series expansion for the case $\alpha\gg1$ is already answered here. How does this differ from the expansion at $\alpha=1$ and in the region $\alpha<1$.
I see that the numerical plot ($\nu=2$) of the integration is the following

Thus something very interesting is happening at the point $\alpha=1$.
Edit:
If it is difficult to find an asymptotic expansion at $\alpha=1$, then at least can we prove that $f(\alpha)$ has a maxima at that point?
 A: For $a$ near 1 the 'multiplication theorem' (see wiki on Bessel functions) can be used:
$$ J_\nu(a \ x) = a^\nu \sum_{n=0}^\infty \frac{1}{n!} \Big( \frac{(1-a^2)x}{2}\Big)^n J_{\nu+n}(x) $$
Keeping the $n=0$ and $n=1$ terms gives
$$\int_0^\infty \frac{x^\nu J_\nu(a \ x)}{e^x-1} dx\sim a^\nu \big(C_0^\nu+C_1^\nu(1-a^2)/2 \big)$$
$$C_0^\nu= \int_0^\infty \frac{x^\nu J_\nu( x)}{e^x-1}dx \ ,\quad 
C_1^\nu= \int_0^\infty \frac{x^{\nu+1} J_{\nu+1}( x)}{e^x-1}dx$$
This gives a very good approximation near $a=1.$ For instance, with $a=0.99,$ and
$\nu = 0.4,$ the two-term approximation is 0.0037% different from the actual and for $\nu = 0.8,$ 0.0096% different.  Of course, if you have to calculate numerically two integrals at exactly $a=1,$ there seems to be little practical advantage.  The exception to this is if you are doing some theoretical work and you need the analytic behavior as $a \to 1,$ and you really don't care so much for the exact values of
$C_0^\nu$ and $C_1^\nu,$ but that they exist.
A: Via expanding the series for the Bessel function and then integrating term-wise, one gets: $$\frac{(2a)^v}{\sqrt{\pi}} \sum_{m=0}^{\infty} (-1)^m \frac{\Gamma(m+v+1/2)}{\Gamma(m+1)} \zeta(2m+2v+1) \: a^{2m} $$
where $ \zeta(m)$ is the Riemann zeta function.
This, however, is only convergent for $a < 1 $ though.
A: Let us use the Abel-Plana formula:
$$\sum_0^\infty f(x)=\frac12 f(0)+\int_0^\infty f(x)dx+i\int_0^\infty \frac{f(-ix)-f(ix)}{e^{2\pi x}-1 }dx$$
which works for “weak bounds” as stated in the article. Let $$f(x)=\frac i2\csc(\pi v)(2\pi x)^v\text I_v(2\pi  xa)$$
of which the Modified Bessel Functions of the First Kind simplify as $$f(-ix)-f(ix)=(2\pi)^v x^v \text J_v(2\pi xa)$$.
Therefore:
$$\sum_{x=0}^\infty \frac i2\csc(\pi v)(2\pi x)^v\text I_v(2\pi  xa) =\frac12 \frac i2\csc(\pi v)(2\pi 0)^v\text I_v(2\pi  0a) +\int_0^\infty \frac i2\csc(\pi v)(2\pi x)^v\text I_v(2\pi  xa) dx+i\int_0^\infty \frac{\frac i2\csc(\pi v)(-2\pi ix)^v\text I_v(-2\pi i xa) -\frac i2\csc(\pi v)(2\pi ix)^v\text I_v(2\pi  ixa)}{e^{2\pi x}-1 }dx\implies i\int_0^\infty \frac{(2\pi)^v x^v \text J_v(2\pi xa)}{e^{2\pi x}-1 }dx+ \int_0^\infty \frac i2\csc(\pi v)(2\pi x)^v\text I_v(2\pi  xa) dx = \sum_{x=0}^\infty\frac i2\csc(\pi v)(2\pi x)^v\text I_v(2\pi  xa)  $$
Rearranging:
$$ i\int_0^\infty \frac{(2\pi)^v x^v \text J_v(2\pi xa)}{e^{2\pi x}-1 }dx+ \int_0^\infty \frac i2\csc(\pi v)(2\pi x)^v\text I_v(2\pi  xa) dx = \sum_{x=0}^\infty\frac i2\csc(\pi v)(2\pi x)^v\text I_v(2\pi  xa) \implies i\int_0^\infty \frac{(2\pi)^v x^v \text J_v(2\pi xa)}{e^{2\pi x}-1 }dx = \sum_{x=0}^\infty\frac i2\csc(\pi v)(2\pi x)^v\text I_v(2\pi  xa) -\int_0^\infty \frac i2\csc(\pi v)(2\pi x)^v\text I_v(2\pi  xa) dx $$
Factoring and substituting $2\pi x=t\implies \frac{dt}{2\pi}= dx,t_1=2\pi 0=0,t_2=2\pi\infty=\infty$:
$$i\int_0^\infty \frac{(2\pi)^v x^v \text J_v(2\pi xa)}{e^{2\pi x}-1 }dx = \sum_{x=0}^\infty\frac i2\csc(\pi v)(2\pi x)^v\text I_v(2\pi  xa) -\int_0^\infty \frac i2\csc(\pi v)(2\pi x)^v\text I_v(2\pi  xa) dx\mathop\implies ^{2\pi x=t}\frac i{2\pi}\int_0^\infty \frac{t^v \text J_v(ta)}{e^t-1 }dt = \sum_{\frac t{2\pi}=0}^\infty\frac i2\csc(\pi v)t^v\text I_v(ta) -\frac1{2\pi}\int_0^\infty \frac i2\csc(\pi v)t^v\text I_v(ta) dt $$
Factoring and some algebra:
$$\int_0^\infty \frac{t^v \text J_v(ta)}{e^t-1 }dt =- 2\pi i\sum_{\frac t{2\pi}=0}^\infty\frac i2\csc(\pi v)t^v\text I_v(ta) - -2\pi i\frac1{2\pi}\int_0^\infty \frac i2\csc(\pi v)t^v\text I_v(ta) dt=  2\pi i \frac i2\csc(\pi v)\sum_{\frac t{2\pi}=0}^\infty t^v\text I_v(ta) - -2\pi i\frac1{2\pi} \frac i2\csc(\pi v)\int_0^\infty t^v\text I_v(ta) dt =-\csc(\pi v)\left[\pi\sum_{\frac t{2\pi}=0}^\infty t^v\text I_v(ta)+\frac12 \int_0^\infty t^v\text I_v(ta) dt\right]$$
So our final result according to the Abel-Plana formula is:
$$f(a)= \sum_0^\infty f(x)=\frac12 f(0)+\int_0^\infty f(x)dx+i\int_0^\infty \frac{f(-ix)-f(ix)}{e^{2\pi x}-1 }dx=-\csc(\pi v)\left[\pi\sum_{\frac x{2\pi}=0}^\infty x^v\text I_v(xa)+\frac12 \int_0^\infty x^v\text I_v(xa) dx\right]= -\csc(\pi v)\left[2^v\pi^{v+1}\sum_{ x=0}^\infty x^v\text I_v(2\pi ax)+\frac12 \int_0^\infty x^v\text I_v(xa) dx\right] $$
Note that $$\int x^v\text I_v(ax)dx=a^v \frac{x^{2v+1}}{2^{v+1}}\Gamma\left(v+\frac12\right)\,_1\tilde{\text F}_2\left(v+\frac12;v+1,v+\frac32;\frac{a^2 x^2}4\right)+C$$
where appears the Regularized $\,_2\text F_1$ Hypergeometric function.
The result may be a difference of large numbers.I hope this helps. You can also try this special case geometric series expansion. Please correct me and give me feedback!
