I asked a similar question before, but now I can formulate it more concretely. I am trying to perform an expansion of the function $$f(x) = \sum_{n=1}^{\infty} \frac{K_2(nx)}{n^2 x^2},$$ for $x \ll 1$. Here, $K_2(x)$ is the modified Bessel function of the second kind. This series is a result of solving the integral $$f(x) = \frac{1}{3}\int_1^\infty \frac{(t^2-1)^{3/2}}{\mathrm{e}^{xt}-1}\mathrm{d}t.$$ The stated result is $$f(x) \approx \frac{\pi^4}{45 x^4} - \frac{\pi^2}{12 x^2}+\frac{\pi}{6x}-\frac{1}{32}\left( \frac{3}{2}-2\gamma+2\ln4\pi-\ln x^2\right)+\mathcal{O}(x^2),$$ where $\gamma$ is the Euler-Mascheroni constant. It agrees numerically with $f(x)$ for small $x$. However, by using the series expansion of the Bessel function $$K_2(nx) = \frac{2}{n^2x^2}-\frac{1}{2}+\frac{1}{2}\sum_{k=0}^\infty \left[\psi(k+1)+\psi(k+3)-\ln\frac{n^2x^2}{4}\right]\frac{\left(\frac{n^2 x^2}{4}\right)^{k+1}}{k!(k+2)!},$$ with $\psi(x)$ being the digamma function and using the zeta regularization for summation over $n$, I am able to reproduce all the terms except $\frac{\pi}{6x}$. I.e., my result is $$f(x) = \frac{2\zeta(4)}{x^4} - \frac{\zeta(2)}{2x^2} + \frac{1}{8}\sum_{k=0}^\infty \left[\left(\psi(k+1)+\psi(k+3)-\ln\frac{x^2}{4}\right)\zeta(-2k) + 2 \zeta'(-2k)\right]\frac{\left(\frac{x^2}{4}\right)^{k}}{k!(k+2)!}.$$ It seems very strange that the $\frac{\pi}{6x}$ term should appear in the expansion since only even powers of $x$ appear in $K_2(nx)$. But, numerically, it is certainly there. How did I miss it?
Edit #1: I just got an idea where the $\frac{\pi}{6x}$ term might come from! Approximating the integral representation of $f(x)$ for $xt \ll 1$ and using the UV regulator $\Lambda$, we have $$f(x) \approx \frac{1}{3} \int_1^\Lambda \frac{(t^2-1)^{3/2}}{xt}\mathrm{d}t \approx \frac{\Lambda^3}{9x} - \frac{\Lambda}{2x} + \frac{\pi}{6x}.$$ OK, so now the question is why doesn't the Bessel series see this term and what to do about it?
Edit #2: The missing term might indicate that the zeta regularization isn't used properly. The term $\frac{\pi}{6x}$ appears right in the middle, separating the convergent sums $\zeta(4)$ and $\zeta(2)$ from the (regularized) divergent sums $\zeta(-2k)$ and $\zeta'(-2k)$. So, the missing term may be the price to pay for using the zeta regularization. Unfortunately, I don't know enough math to come to any decisive conclusion.
Edit #3: In edit #1 I argued that the missing term $\frac{\pi}{6x}$ comes from expansion of the exponential in the denominator. When I expand the numerator in the binomial series $$f(x) = \frac{1}{3} \int_1^\infty \frac{t^3}{\mathrm{e}^{xt}-1}\sum_{k=0}^\infty \binom{3/2}{k}(-t^{-2})^k\mathrm{d}t$$ and perform the integration, I again obtain my original result without $\frac{\pi}{6x}$. So, depending on what I choose to expand, I obtain different (incomplete) results??