This may sound as a mathematical question, but it should be very familiar to physicists. 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 result should be $$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?
I think that the missing term indicates 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 formal math to, so to speak, "repair the damage". Any help on the subject is more than welcome.