# Missing term in series expansion

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??

Using Mellin transforms we can find the asymptotics in an intuitive and straightforward manner. We have that $$\mathfrak{M}(K_2(x)/x^2; s) = 2^{s-4} \Gamma(s/2) \Gamma(s/2-2),$$ so that viewing $f(x)$ as a harmonic sum, we get $$\mathfrak{M}(f(x); s) = f^*(s) = 2^{s-4} \Gamma(s/2) \Gamma(s/2-2) \zeta(s).$$ Now invert to get the expansion of $f(x).$ I will give a table of the contributions from the poles down to the pole at $s=-6.$ $$\begin{array} \operatorname{Res}(f^*(s) \, x^{-s}; s=4) & = & 1/45\,{\frac {{\pi }^{4}}{{x}^{4}}} \\ \operatorname{Res}(f^*(s) \, x^{-s}; s=2) & = & -1/12\,{\frac {{\pi }^{2}}{{x}^{2}}} \\ \operatorname{Res}(f^*(s) \, x^{-s}; s=1) & = & 1/6\,{\frac {\pi }{x}} \\ \operatorname{Res}(f^*(s) \, x^{-s}; s=0) & = & -1/16\,\ln \left( 4\,\pi \right) +1/16\,\gamma+1/16\, \ln \left( x \right) -{\frac {3}{64}} \\ \operatorname{Res}(f^*(s) \, x^{-s}; s=-2) & = & {\frac {1}{96}}\,\zeta \left( 1,-2 \right) {x}^{2} \\ \operatorname{Res}(f^*(s) \, x^{-s}; s=-4) & = & {\frac {1}{3072}}\,\zeta \left( 1,-4 \right) {x}^{4} \\ \operatorname{Res}(f^*(s) \, x^{-s}; s=-6) & = & {\frac {1}{184320}}\,\zeta \left( 1,-6 \right) {x}^{6}. \end{array}$$ As to how the Mellin transform of $K_2(x)$ is calculated, I can offer some ideas. Start with the known integral representation $$K_\alpha(x) = \int_0^\infty e^{-x \cosh t} \cosh (\alpha t) \; dt,$$ so that $$\mathfrak{M}(K_2(x); s) = \int_0^\infty \int_0^\infty e^{-x \cosh t} \cosh (2t) \; dt \; x^{s-1} \; dx.$$ This becomes $$\int_0^\infty \cosh (2t) \int_0^\infty e^{-x \cosh t} x^{s-1} \; dx \; dt = \Gamma(s) \int_0^\infty \frac{\cosh (2t)}{(\cosh t)^s} dt.$$ Now using $$\cosh(2t) = 2\cosh(t)^2 - 1,$$ we obtain $$\Gamma(s) \left(2 \int_0^\infty \frac{1}{(\cosh t)^{s-2}} dt - \int_0^\infty \frac{1}{(\cosh t)^s} dt \right).$$ Furthermore, $$\int_0^\infty \frac{1}{(\cosh t)^s} dt = 2^s \int_0^\infty \frac{1}{(e^t + e^{-t})^s} dt = 2^s \int_0^\infty \frac{1}{e^{ts}} \frac{1}{(1 + e^{-2t})^s} dt$$ which is $$2^s \int_0^\infty \frac{1}{e^{ts}} \sum_{q\ge 0} (-1)^q \binom{q+s-1}{q} e^{-2qt} dt = 2^s \int_0^\infty \sum_{q\ge 0} (-1)^q \binom{q+s-1}{q} e^{-(2q+s)t} dt$$ or $$2^s \sum_{q\ge 0} (-1)^q \binom{q+s-1}{q} \int_0^\infty e^{-(2q+s)t} dt = 2^s \sum_{q\ge 0} (-1)^q \binom{q+s-1}{q} \frac{1}{2q+s}.$$ Now this last term is $$2^s \frac{1}{s} \; _2F_1(s/2, s; 1+s/2; -1).$$ But we have $$_2F_1(a, 2a; a+1; -1) = \frac{1}{2a} \frac{\Gamma(a+1)^2}{\Gamma(2a)}$$ as can be seen e.g. here.
Hence $$\int_0^\infty \frac{1}{(\cosh t)^s} dt = 2^s \frac{1}{s^2} \frac{\Gamma(s/2+1)^2}{\Gamma(s)}.$$ Concluding, we have shown that $$\mathfrak{M}(K_2(x);s) = 2 \Gamma(s) 2^{s-2} \frac{1}{(s-2)^2} \frac{\Gamma(s/2)^2}{\Gamma(s-2)} - \Gamma(s) 2^s \frac{1}{s^2} \frac{\Gamma(s/2+1)^2}{\Gamma(s)} \\ = 2^{s-1} \frac{s-1}{s-2}\Gamma(s/2)^2 - 2^s \frac{1}{s^2} \Gamma(s/2+1)^2 = 2^{s-1} \frac{s-1}{s-2}\Gamma(s/2)^2 - 2^{s-2} \Gamma(s/2)^2 = 2^{s-2} \Gamma(s/2)^2 \left( 2 \frac{s-1}{s-2} - 1\right) = 2^{s-2} \Gamma(s/2)^2 \frac{s}{s-2} = 2^{s-2} \Gamma(s/2)^2 \frac{s/2}{s/2-1} = 2^{s-2} \Gamma(s/2+1) \Gamma(s/2-1).$$