How does one regularize a divergent integral of the form,
$$
I = \int_0^\infty dx \, \cosh^4 x \hspace{0.5 in} ?
$$
Regularizing via simple subtraction of divergences (as we commonly do in physics) is not a viable option, since expanding $\cosh^4x$ in a power series reveals an infinite number of divergent contributions, with no finite piece to remain. However, I do know that there must be some regularization technique out there enabling us to evaluate this, since we have a similar result
$$
\int_0^\infty dx \, \cosh^4 x \sinh x = -{8 \pi \over 15}
$$
in common use. Any tips on how to obtain this sort of thing (or even reference to a source where this is calculated) would be very helpful!
[Math] Regularization of integral
definite integralsdivergence-operatorintegrationregularization
Related Solutions
Regularization is a way to make a divergent integral into a parametric family of finite integrals. Say, we need to evaluate $\int f(x) \mathrm{d} x$. We introduce (and this is a craft in itself) $f_\alpha(x)$, such that $\lim_{\alpha \to 0} f_\alpha(x) = f(x)$, such that $\int f_\alpha(x) \mathrm{d} x$ exists. We then perform, unjustified, interchanging of integration and taking the limit: $$ \int f(x) \mathrm{d} x = \int \lim_{\alpha \to 0} f_\alpha(x) \mathrm{d} x \stackrel{?}{=} \lim_{\alpha \to 0} \int f_\alpha(x) \mathrm{d} x $$
In your particular case, integral can be defined using Abel summations, since $$ \int_0^\infty \sin(x) \mathrm{d} x = \sum_{n=0}^\infty \int_{\pi n}^{\pi (n+1)} \sin(x) \mathrm{d} x = \sum_{n=0}^\infty 2 (-1)^n \stackrel{\text{Abel}}{=} 1 $$ Different regularization of divergent series do not need to give the same result, however, under certain conditions they do.
Coming back to your integral, its regularization is done, just like you did, by taking advantage of the hypergeometric nature of the integrand, and representing the integral as limit of Mellin convolutions (i.e. $\int_0^\infty h_\alpha(x) \sin(x) \mathrm{d} x$, where $h_\alpha(x)$ is another hypergeometric function, that makes the integral convergent, and such that $\lim_{\alpha \downarrow 0} h_\alpha(x) = 1$): $$ \lim_{\alpha \downarrow 0} \int_0^\infty \frac{\sin(x)}{x^\alpha}\mathrm{d} x = \lim_{\alpha \downarrow 0} \left( \Gamma(1-\alpha) \cos\left( \frac{\pi \alpha}{2} \right) \right) = 1 $$ $$ \lim_{\alpha \downarrow 0} \int_0^\infty \mathrm{e}^{-\alpha x} \sin(x) \mathrm{d} x = \lim_{\alpha \downarrow 0} \frac{1}{\alpha^2 + 1} = 1 $$ $$ \lim_{\alpha \downarrow 0} \int_0^\infty \frac{\sin(x)}{1+\alpha x} \mathrm{d} x = \lim_{\alpha \downarrow 0} \left( \frac{2 \sin \left(\frac{1}{\alpha }\right) \text{Ci}\left(\frac{1}{\alpha }\right)+\cos \left(\frac{1}{\alpha }\right) \left(\pi -2 \text{Si}\left(\frac{1}{\alpha }\right)\right)}{2 \alpha } \right) = 1 $$
Edit
The formula proposed by Tom may be obtained with the classical method : $$\int_0^{\infty} \frac{x^{s-1}}{e^x(e^x-1)} dx $$ $$ = \sum_{n=2}^{\infty} \int_0^{\infty} x^{s-1} e^{-nx} dx$$ $$ = \sum_{n=2}^{\infty} \int_0^{\infty} \left(\frac{t}{n}\right)^{s-1} e^{-t} dt/n $$ $$ = \left(\sum_{n=2}^{\infty} \frac{1}{n^s}\right) \int_0^{\infty} t^{s-1} e^{-t} dt $$ $$ = \left(\zeta(s)-1\right) \Gamma(s)$$
As $s\to 0$ we have $\displaystyle \left(\zeta(s)-1\right) \Gamma(s)\sim -\frac 3{2s} +\frac 32\gamma-\frac{\ln(2\pi)}2$ (with $\gamma=\lim_{s\to 0} \frac 1s-\Gamma(s)$ the Euler constant $0.5772156649\cdots$ and $\frac{-\ln(2\pi)}2=\zeta'(0)$) which is still divergent...
I am not sure that we can avoid this because of the equivalence $\int \frac 1{x^2} dx$ near $0$ but perhaps that your 'regularization' could consist in the limit : $\displaystyle \lim_{\epsilon\to 0} \frac{\left(\zeta(-\epsilon)-1\right) \Gamma(-\epsilon)+\left(\zeta(\epsilon)-1\right) \Gamma(\epsilon)}2=\frac 32\gamma-\frac{\ln(2\pi)}2$
or a symmetrical pondered integral centered at $0$ with the same result...
Corresponding series
Let's substitute this definition of the Mittag-Leffler function : $$\operatorname{E}_s(z)=\sum_{k=0}^\infty \frac{z^k}{\Gamma(sk+1)}$$ in your integral : $$I(s)=s\int_{0}^{\infty}\frac{E_{s}(x^{s})-1}{xe^{x}(e^{x}-1)}dx$$ $$I(s)=s\int_{0}^{\infty}\frac{\sum_{k=1}^\infty \frac{x^{sk}}{\Gamma(sk+1)}}{xe^{x}(e^{x}-1)}dx=s\sum_{k=1}^\infty \frac 1{\Gamma(sk+1)}\int_{0}^{\infty}\frac{ x^{sk-1}}{e^{x}(e^{x}-1)}dx$$
We may use our previous result to rewrite $I(s)$ as : $$I(s)=s\sum_{k=1}^\infty \frac {\left(\zeta(sk)-1\right) \Gamma(sk)}{\Gamma(sk+1)}$$ this becomes the simple (and convergent if $\Re(s)>0$ and $\frac 1s$ not integer at least...) : $$\boxed{I(s)=\sum_{k=1}^\infty \frac {\zeta(sk)-1}{k}}$$
Approximation and Visualization
We may remove the simple poles in this series at $s=\frac 1n$ for $n\in \mathbb{N}^*$ by subtracting $\frac 1{sk-1}$ at the numerator (Stieltjes expansion) with the hope of getting something simpler.
Indeed for $s$ small I found this rather accurate numerical approximation (relative error $<\frac 1{20000}$ in the rectangle $\Re(s)\in (10^{-5},2),\ |\Im(s)|<1$) :
$$\sum_{k=1}^\infty \frac {\zeta(sk)-1-\frac 1{sk-1}}{k}\approx\frac {\ln(s)}2-0.95973512990747991661376-0.04053073339766363s+0.00026485233s^2$$
I think that the $\frac {\ln(s)}2$ term is exact (from numerical 'scaling properties') as well as the provided digits of the last three constants.
Since $\displaystyle \sum_{k=1}^\infty \frac 1{k(sk-1)}=-\gamma -\psi\left(1-\frac 1s\right)$ with $\psi$ the digamma function we have the approximation : $$I(s)\approx -\gamma -\psi\left(1-\frac 1s\right) +\frac {\ln(s)}2-0.95973512990747991661376-0.04053073339766363s+0.00026485233s^2$$
that allows us to produce this picture of the imaginary part of $I(s)$ for $\Re(s)>0$ :
We observe a jump of amplitude $\pi$ (near the line $x=\Re(s)=0$ at the right around the line $y=\Im(s)=0$) coming from $\frac {\ln(s)}2$.
Let's observe that our expansion is analytic everywhere except at the simple poles $s=\frac 1n$ and the branch point $s=0$ allowing us to show the pictures of the imaginary and real part of the (possible) Riemann surface of the analytic extension of $I(s)$ on the whole plane (imagine an infinity of sheets superposed around the branch point $0$ generated by the $\frac {\ln(s)}2$ term) :
Asymptotic Expansion (direct method)
We want to rewrite following variant of $I(s)$ as a Maclaurin series $P(s)$ : $$J(s):=\sum_{k=1}^\infty \frac {\zeta(sk)-1-\frac 1{sk-1}}{k}-\dfrac{\ln(s)}2$$
Let's split the $-1$ at the numerator in two $-\frac 12$ terms to study the asymptotic behavior of this series (and obtain convergence of the individual 'power $\log$ integrals' we will see later...) : $$J(s)= \lim_{N\to\infty} J_N(s)\ \ \text{with}\ \ J_N(s)=\sum_{k=1}^N \frac {\zeta(sk)-\frac 12-\frac 1{sk-1}}{k}-\frac 12 \sum_{k=1}^N \frac 1k -\frac{\ln(s)}2$$
Let's first show that $\ \lim_{N\to\infty} J_N(s)$ exists for any $s>0$ such that $\frac 1s \notin \mathbb{N}$ (so that all the terms of the sum are defined).
Since $\ \sum_{k=1}^N \frac 1k= H_n=\gamma +\ln\bigl(N+\frac 12\bigr) + \operatorname{O}\bigl(\frac 1N\bigr)\ $ we have $\frac 12 \sum_{k=1}^N \frac 1k +\frac{\ln(s)}2=\frac {\gamma}2+\frac 12 \ln\bigl(s\bigl(N+\frac 12\bigr)\bigr)+\operatorname{O}\bigl(\frac 1N\bigr)=\frac {\gamma}2+\frac {\ln(sN)}2+\operatorname{O}\bigl(\frac 1N\bigr)$
For $sN \gg 1$ we have (setting $x:=sk$) :
$$s\sum_{k=N}^M \frac {\zeta(sk)-\frac 12-\frac 1{sk-1}}{sk}\sim \int_{sN}^{sM} \frac {\zeta(x)-\frac 12-\frac 1{x-1}}{x}\,dx$$
the error is of order $\frac 1N$ because of the $H_M-H_N$ terms and because $\zeta(x)-1 \sim 2^{-x}$ for $x \gg 1$.
A precise evaluation of the last integral is : $\left[\frac 32 \ln(x) -\ln(x-1)+\sum_{k>1} \operatorname{Ei}(-x\ln(k))\right]^{sM}_{sN}$ (with $\operatorname{Ei}$ the exponential integral and $\operatorname{Ei}(-u)\sim \frac {e^{-u}}{-u}$ for $u\gg 1$) but this isn't needed since the integral is clearly equivalent to $\frac 12 \left[\ln(x)\right]^{sM}_{sN}\ $ with again an error of order $\frac 1N$ for $sN\gg 1$.
Both $\frac {\ln(sN)}2$ terms cancel in $J_N(s)$ for large $sN$ so that the limit $J(s)$ exists for $s>0$ (for $\frac 1s \in \mathbb{N}$ $J(s)$ may be defined as a limit since we removed the poles).
Let's use the Euler Maclaurin asymptotic formula with $f(k)=\frac {\zeta(sk)-\frac 12-\frac 1{sk-1}}{k}$ to expand directly $J_N(s)$ : $$\sum_{k=0}^N f(k) \sim \int_0^N f(x) dx +\frac {f(0)+f(N)}2+\sum_{k=1}^\infty \frac{B_{2k}}{(2k)!} \bigl(f^{(2k-1)}(N)-f^{(2k-1)}(0)\bigr)$$ becomes : $$ s\sum_{k=0}^N \frac {\zeta(sk)-\frac 12-\frac 1{sk-1}}{sk}\sim \int_0^{sN} \frac {\zeta(x)-\frac 12-\frac 1{x-1}}{x}\,dx + \frac {f(0)+f(N)}2+\sum_{k=1}^\infty \frac{B_{2k}}{(2k)!} \left(f^{(2k-1)}(N)-f^{(2k-1)}(0)\right)$$
Note that we started at $k=0$ with $\displaystyle\ f(0)=\lim_{k\to 0}\ s\frac {\zeta(sk)-\frac 12-\frac 1{sk-1}}{sk}=(1+\zeta'(0))s\ $
and that, for any $n\in \mathbb{N}^*$, we have :
$f^{(n-1)}(0)=\lim_{N\to 0}\ s\left(\frac d{dN}\right)^{n-1}\ \left(\frac {\zeta(sN)-\frac 12-\frac 1{sN-1}}{sN}\right)=\lim_{N\to 0}\ s\left(\frac {s\cdot d}{d(sN)}\right)^{n-1} (\cdots)=\lim_{x\to 0}\ s^n\left(\frac {\zeta(x)-\frac 12-\frac 1{x-1}}x\right)^{(n-1)}=s^n\lim_{x\to 0} \frac{\left(\zeta(x)-\frac 1{x-1}\right)^{(n)}}n$
or simply $\ \displaystyle f^{(n-1)}(0)=\frac{n!+\zeta^{(n)}(0)}n s^n\ $ for $n>0$
For $n\ge 0$ and $N \gg 1$ we may show that $f^{(n)}(N)\sim \frac {(-1)^n n!}{2 N^{n+1}}\ $ so that the Euler Maclaurin asymptotic formula becomes (after removing $\frac{H_N+\ln(s)}2$ on both sides to get $J_N(s)$ at the left) :
$$J_N(s)= s\sum_{k=0}^N \frac {\zeta(sk)-1-\frac 1{sk-1}}{sk}-\ln(s)\sim \int_0^{sN} \frac {\zeta(x)-\frac 12-\frac 1{x-1}}x\,du - \frac {\gamma}2-\frac {\ln(sN)}2+\operatorname{O}\bigl(\frac 1N\bigr)+ \frac {f(0)}2-\sum_{k=1}^\infty \frac{B_{2k}}{(2k)!} f^{(2k-1)}(0)$$
At this point let's rewrite $\ln(sN)=\int_0^{sN} \frac 1{x+1}\,du+\operatorname{O}\bigl(\frac 1N\bigr)$ and go at the limit $N\to \infty$ : $$J(s)\sim - \frac {\gamma}2+\int_0^\infty \frac {\zeta(x)-\frac 12-\frac 1{x-1}}x-\frac 1{2(x+1)}\,dx + \frac {f(0)}2-\sum_{k=1}^\infty \frac{B_{2k}}{(2k)!} f^{(2k-1)}(0)$$
We recognize a Maclaurin series in $s$ : $$J(s)\sim K_0 -K_1 \frac s2-\sum_{k=1}^\infty K_{2k} \frac{B_{2k}}{(2k)!} s^{2k} $$
with the constant term given by $\displaystyle\ \boxed{K_0:= -\frac{\gamma}2 +\int_0^\infty \frac {\zeta(x)-\frac 1{x-1} -\frac 12 }x -\frac 1{2(1+x)}\,dx}$
(the numerical value is $K_0=-0.95973512990\cdots$ as wished ; I don't know a closed form...)
and the $K_n$ coefficients of $s^n$ by $\boxed{K_n:=\frac{n!+\zeta^{(n)}(0)}n}\ \text{for} \ n>0$
The wished analytic extension of $I(s)$ will be :
$$\boxed{I(s)\sim -\gamma -\psi\left(1-\frac 1s\right) +\frac {\ln(s)}2+K_0 -K_1 \frac s2-\sum_{k=1}^\infty K_{2k} \frac{B_{2k}}{(2k)!} s^{2k}}$$
Short table of $K_n$ values :
($\gamma_n$ is the nth Stieltjes constant : $\gamma_0= \gamma,\ \gamma_1= 0.07281584548367672486058\cdots$)
$\displaystyle K_1=1+\zeta'(0)=1-\frac{\ln(2\pi)}2=0.0810614667953272582\cdots$
$\displaystyle K_2=1+\frac{\zeta''(0)}2=1-\frac{\ln(2\pi)^2}4+\frac{\gamma^2}4-\frac{\pi^2}{48}-\frac{\gamma_1}2=
-0.0031782279542924256\cdots$
$K_3=\frac{3!+\zeta'''(0)}3=-0.00157038895408481592\cdots$
$K_4=\frac{4!+\zeta^{(4)}(0)}4=0.0007242029965730102531\cdots$
($-\frac {K_1}2$ and $-\frac {K_2}{12}$ may be compared to the coefficients for $s$ and $s^2$ obtained previously : $-0.040530733397\cdots$ and $0.00026485233\cdots$).
(it is interesting to observe that $\zeta^{(n)}(0)\approx -n!$ as well as $\zeta^{(n)}\left(\frac 12\right)\approx -2^{n+1}n!$)
Asymptotic Expansions (the same result after an excursion with Mittag-Leffler, another expression for $K_0$)
To get the exact terms of the expansion of $I(s)$ let's provide the Maclaurin expansion in $s$ of $s\left(E_s(x^s)-1\right)$ (the coefficients are functions of $x$ that we will have to integrate later over $\mathbb{R}^+$) :
$$s\left(E_s(x^s)-1\right)=C_0(x) e^x-\frac s2-\frac{s^2}{12}\left[\ln(x)+\gamma\right]+\frac{s^4}{720}\left[\ln(x)^3+3\gamma\ln(x)^2+3\left(\gamma^2-\zeta(2)\right)\ln(x)+\gamma^3-3\gamma\zeta(2)+2\zeta(3)\right]-\frac{s^6}{30240}\left[\ln(x)^5+5\gamma\ln(x)^4+10\left(\gamma^2-\zeta(2)\right)\ln(x)^3+10\left(\gamma^3-3\gamma \zeta(2)+2\zeta(3)\right)\ln(x)^2+5\left(\gamma^4-6\gamma^2\zeta(2)+8\gamma \zeta(3)+\pi^4/60\right)\ln(x)+\gamma^5-10\gamma^3\zeta(2)+20\gamma^2\zeta(3)+\gamma\pi^4/12-20\zeta(2)\zeta(3)+24\zeta(5) \right]+P_8(\ln(x))\operatorname{O}\bigl(s^8\bigr)$$
The term independent of $s$ : $C_0(x) e^x$ with $\displaystyle C_0(x)=\int_0^1 \frac{\gamma(t,x)}{\Gamma(t)}\,dt\ $ (and $\gamma(t,x)$ the lower 'incomplete gamma function') was obtained by considering the limit as $n\to \infty$ of formula (3.4) from Haubold's paper (with $s=\frac 1n$) : $$E_{\frac 1n}\left(x^{\frac 1n}\right)=e^x\left[1+\sum_{r=1}^{n-1}\frac{\gamma\left(1-\frac rn ,x\right)}{\Gamma\left(1-\frac rn\right)}\right]$$
A shorter expression for the Maclaurin expansion (that could again be obtained with Euler-Maclaurin) should be :
$$\boxed{s\left(E_s(x^s)-1\right)\sim C_0(x) e^x-\frac s2-\sum_{n=2}^\infty s^n\frac{B_n}{n!}\ \sum_{k=0}^{n-1} \binom{n-1}{k} S_k \ln(x)^{n-1-k}} $$
with $B_n$ the 'Bernoulli numbers' ($B_{2n+1}=0$ for $n>0$) and $S_k$ the coefficients appearing in the expansion of the reciprocal gamma function at $z=0$ (and at $z=1$) : $$\frac 1{\Gamma(s)}=\sum_{k=0}^\infty \frac {S_k}{k!} s^{k+1}\ \ \text{so that}\ \ S_k=\lim_{s\to 0} \left(\frac 1{\Gamma(s+1)}\right)^{(k)}$$
$S_0=1$
$S_1=\gamma$
$S_2=\gamma^2-\zeta(2)$
$S_3=\gamma^3-3\gamma \zeta(2)+2\zeta(3)$
$S_4=\gamma^4-6\gamma^2\zeta(2)+8\gamma \zeta(3)+\pi^4/60$
$S_5=\gamma^5-10\gamma^3\zeta(2)+20\gamma^2\zeta(3)+\gamma\pi^4/12-20\zeta(2)\zeta(3)+24\zeta(5)$
$C_0$ too may be rewritten with the $S_k$ terms as $\displaystyle C_0(x)=\sum_{k=0}^\infty \frac{S_k}{(-\ln(x))^{k+1}}$.
We will again start with : $\displaystyle J(s)= \lim_{N\to\infty} J_N(s)$ with $\ \displaystyle J_N(s)=\sum_{k=1}^N \frac {\zeta(sk)-\frac 12-\frac 1{sk-1}}{k}-\frac 12 \sum_{k=1}^N \frac 1k -\frac{\ln(s)}2$
but will use Mittag-Leffler after that : $$J_N(s)=s\sum_{k=1}^N \frac {\zeta(sk)\Gamma(sk)-\frac {\Gamma(sk)}2-\frac {(sk-1)\Gamma(sk-1)}{sk-1}}{(sk)\Gamma(sk)}-\frac 12 \sum_{k=1}^N \frac 1k -\frac{\ln(s)}2$$
$$J_N(s)=s\sum_{k=1}^N \frac {\zeta(sk)\Gamma(sk)-\frac {\Gamma(sk)}2-\Gamma(sk-1)}{\Gamma(sk+1)}-\frac 12 \sum_{k=1}^N \frac 1k -\frac{\ln(s)}2$$
Let's rewrite everything in integral form (using $\ \sum_{k=1}^N \frac 1k= \gamma +\int_0^N \frac 1{1+x} + \operatorname{O}\bigl(\frac 1N\bigr)\ $, $\Gamma(u)=\int_0^\infty t^{u-1} e^{-t} dt\ $ and considering the limit as $N\to \infty$) :
$$J_N(s)=\sum_{k=1}^N \left( \frac s{\Gamma(sk+1)}\int_{0}^{\infty}\frac{x^{sk-1}}{e^x-1}-\frac {x^{sk-1}}{2e^x}-\frac {x^{sk-1}}{x e^x}\,dx\right)-\frac{\gamma}2 -\frac{\ln(N s)}2+\operatorname{O}\bigl(\frac 1N\bigr)$$
$$J_N(s)=-\frac{\gamma}2+\int_{0}^{\infty}\frac{\sum_{k=1}^N \frac {sx^{sk}}{\Gamma(sk+1)}}{xe^{x}}\left(\frac {e^x}{e^{x}-1}-\frac 12-\frac 1x\right)\,dx -\frac{\ln(N s)}2+\operatorname{O}\bigl(\frac 1N\bigr)$$
$$J_N(s)=-\frac{\gamma}2 +\int_{0}^{\infty}\frac{\sum_{k=1}^N \frac {sx^{sk}}{\Gamma(sk+1)}}{xe^{x}}\left(\frac {e^x}{e^{x}-1}-\frac 12-\frac 1x\right)+\frac{e^{-N s x}-e^{-x}}{2x}\,dx +\operatorname{O}\bigl(\frac 1N\bigr)$$
that becomes at the limit $N \to \infty$ :
$$J(s)\sim -\frac{\gamma}2 +\int_{0}^{\infty}\frac{\sum_{k=1}^\infty\frac {sx^{sk}}{\Gamma(sk+1)}}{xe^{x}}\left(\frac {e^x}{e^{x}-1}-\frac 12-\frac 1x\right)-\frac 1{2(1+x)}\,dx$$
$$\boxed{J(s)\sim -\frac{\gamma}2 +\int_{0}^{\infty}\frac{s \left(E_{s}(x^{s})-1\right)}{xe^{x}}\left(\frac 1{e^{x}-1}-\frac 1x+\frac 12\right)-\frac 1{2(1+x)}\,dx}$$
From the MacLaurin expansion of $s \left(E_{s}(x^{s})-1\right)$ in powers of $s$, we get :
$$J(s)\sim K_0 -K_1 \frac s2-\sum_{n=2}^\infty K_n \frac{B_n}{n!} s^n $$
and the wished analytic extension of $I(s)$ :
$$\boxed{I(s)\sim -\gamma -\psi\left(1-\frac 1s\right) +\frac {\ln(s)}2+K_0 -K_1 \frac s2-\sum_{n=2}^\infty K_n \frac{B_n}{n!} s^n}$$
with the constant term given by : $$K_0:= -\frac{\gamma}2 +\int_0^\infty \frac{C_0(x)}x\left(\frac 1{e^{x}-1}-\frac 1x+\frac 12\right)-\frac 1{2(1+x)}\,dx$$
with the same numerical value and the $K_n$ coefficients of $s^n$ for $n>0$ by :
$$K_n:= \sum_{m=0}^{n-1} f_m \binom{n-1}{m} S_{n-1-m}$$
with :
$$f_m:=\int_{0}^{\infty}\frac{\ln(x)^m}{xe^{x}}\left(\frac 1{e^{x}-1}-\frac 1x+\frac 12\right)\,dx$$
Let's use (with a derivation similar to the previous one!) :
$$Z(s):=\left(\zeta(s)-\frac 1{s-1}-\frac 12\right) \Gamma(s)=\int_0^{\infty} \frac{x^s}{xe^{x}}\left(\frac 1{e^{x}-1}-\frac 1x+\frac 12\right)\,dx$$ (the idea of the additional term $-\frac 1x+\frac 12$ appeared in the Hermite-Stieltjes letters of 1885 allowing to extend the integral giving zeta near $x=0$ by removing the $\frac 1{x^2}$ and $\frac 1x$ singularities ; just some pages later we find the famous claim of Stieltjes of a proof of R.H. !)
since $x^s=e^{s\ln(x)}$ we get : $$f_m= Z^{(m)}(0)=\lim_{s\to 0}\left(\frac {d}{ds}\right)^m \left[\left(\zeta(s)-\frac 1{s-1}-\frac 12\right)\Gamma(s)\right]$$
From the binomial theorem it is clear that : $$K_{n+1}= \sum_{m=0}^n \binom{n}{m} f_m S_{n-m}= \lim_{s\to 0}\sum_{m=0}^n \binom{n}{m} Z^{(m)}\left(\frac 1{\Gamma(s+1)}\right)^{(n-m)}= \lim_{s\to 0} \left(\frac {Z(s)}{\Gamma(s+1)}\right)^{(n)}$$
allowing some simplifications : $$K_n=\lim_{s\to 0} \Bigl(\frac{\zeta(s)-\frac 1{s-1}-\frac 12}s\Bigr)^{(n-1)}=\lim_{s\to 0} \frac{\left(\zeta(s)-\frac 1{s-1}\right)^{(n)}}n\ \text{for} \ n>0$$ $$\text{i.e.}\ \ \boxed{K_n=\frac{n!+\zeta^{(n)}(0)}n}\ \text{for} \ n>0$$ i.e. the same result as previously
We could continue this very interesting subject and study the multiplicative partition (“Factorisatio Numerorum”) obtained by rewriting $\displaystyle I(s)=\sum_{k=1}^\infty \frac {\zeta(sk)-1}{k}$ as the logarithm of $\displaystyle \prod_{n>1} \frac 1{1-n^{-s}}$ getting this Dirichlet series or... but I'll have to make all this visible one day so...
Sorry if all this is not exactly what you hoped...
Best Answer
One may note that:
$$I=\int_0^\infty\cosh^4(x)~\mathrm dx=\frac1{16}\int_0^\infty e^{4x}+4e^{2x}+6+4e^{-2x}+e^{-4x}~\mathrm dx$$
Using $\displaystyle\int_0^\infty e^{-ax}~\mathrm dx=\frac1a$, we may regularize our integral to
$$I=\int_0^\infty\frac38~\mathrm dx$$
Which may be regularized to $\frac38\zeta(0)=-\frac3{16}$.
Same method for the second integral.