Zeta regularization gone bad 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.
 A: Doing this kind of regularization in a series of which each term is more divergent than the ones before is rather dangerous. For example, one might compute the following series:
\begin{equation}
\sum_{n=1}^\infty \frac1{n^2+\Lambda} = \frac1\Lambda \sum_{k=0}^\infty \left(-\frac1\Lambda\right)^k \sum_{n=1}^\infty n^{2k} \stackrel{??}= 0
\end{equation}
due to $\zeta(-2k) = 0 \ \forall k\in\mathbb N$. This is obviously wrong as the correct result is a hyperbolic tangent:
\begin{equation}
\sum_{n=1}^\infty \frac1{n^2+\Lambda} = -\frac1{2\Lambda} + \frac\pi{2\sqrt\Lambda} \coth\sqrt\Lambda\pi \not= 0 \;.
\end{equation}
To better understand, you could, for example, introduce a regulator $e^{-n\epsilon}$ in your sum to make it convergent, and then take the limit $\epsilon\to0$. You will find the result with naive zeta regularization, plus an infinite series of ever more divergent terms. In order to find the correct result, you will have to somehow deal with this series first. In applications such as quantum field theory we know (due to renormalizability) that these divergences cancel in physical quantities, but only if there's finitely many of them (so no infinite series of them).
