変形Bessel関数の４乗のMellin変換

こちら（ http://dlmf.nist.gov/10.32.E17 ）に載っている変形Beesel関数の積の積分表示
$$K_v(x)K_\mu(x)=2\int_0^\infty K_{v\pm\mu}(2x\cosh(t))\cosh((v\mp\mu)t)dt$$
において$\nu=\mu$と置いた式です。

こちら（ http://dlmf.nist.gov/10.43.E22 ）に載っている式
$$\int_0^\infty x^{s-1}e^{-ax}K_v(x)dx=\sqrt{\frac\pi2}\Gamma(s-v)\Gamma(s+v)(a^2-1)^\frac{1-2s}{4}P_{v-\frac12}^{\frac12-s}(a)$$
のLegendre陪関数を超幾何関数に書き換えたものです。（$a$の範囲によってLegendre陪関数の分岐の取り方が変わりますがあまり気にしなくていいです）

こちら（ http://dlmf.nist.gov/10.32.E18 ）に載っている式
$$K_v(x)K_v(y)=\frac12\int_0^\infty e^{-\frac{t}2-\frac{x^2+y^2}{2t}}K_v\left(\frac{xy}{t}\right)\frac{dt}{t}$$
で$x,y$を$ax,bx$と置いたものです。

補題４の表示を用いてから積分の順序を変え、$abx^2/t\to x$と変数変換をすると次のように書けます。
$$\int_0^\infty x^{s-1}K_v(ax)K_v(bx)dx=\frac1{4(ab)^{s/2}}\int_0^\infty t^{s/2-1}e^{-t/2}dt\int_0^\infty x^{s/2-1}e^{-\frac{a^2+b^2}{2ab}x}K_v\left(x\right)dx $$
これに補題３を適応することで得られます。

こちら（ http://dlmf.nist.gov/15.8.E21 ）に載っているものです。

補題２の積分表示を用いてから積分の順序を交換すると次のように書けます。
$$\begin{align*} \int_0^\infty &x^{s-1}K_v(x)^4dx=4\int_0^\infty\int_0^\infty\int_0^\infty x^{s-1}K_{2v}(2x\cosh(t))K_{2v}(2x\cosh(y))dxdtdy\end{align*}$$
補題５のMellin変換を適応して$1/\cosh(t)\to t,1/\cosh(y)\to y$と変数変換すると次を得ます。
$$\begin{align*}\frac{\sqrt{\pi}}{2^s}\frac{\Gamma(\frac{s}2)\Gamma(\frac{s}2-2v)\Gamma(\frac{s}2+2v)}{\Gamma(\frac{s+1}2)}\int_0^1\int_0^1\frac{(ty)^{s/2-1}}{\sqrt{(1-t^2)(1-y^2)}}\hyF{2}{1}{\frac{s}2-2v,\frac{s}2+2v}{\frac{s+1}2}{-\frac{(t-y)^2}{4ty}}dtdy \end{align*}$$
$u=t/y$と変数変換して、積分範囲に注意して積分の順序を交換するとこうなります。
$$\begin{align*}&\frac{\sqrt{\pi}}{2^s}\frac{\Gamma(\frac{s}2)\Gamma(\frac{s}2-2v)\Gamma(\frac{s}2+2v)}{\Gamma(\frac{s+1}2)}\int_0^1\int_0^{1}\frac{y^{s-1}u^{s/2-1}}{\sqrt{(1-y^2)(1-y^2u^2)}}\hyF{2}{1}{\frac{s}2-2v,\frac{s}2+2v}{\frac{s+1}2}{-\frac{(1-u)^2}{4u}}dydu\\ &+\frac{\sqrt{\pi}}{2^s}\frac{\Gamma(\frac{s}2)\Gamma(\frac{s}2-2v)\Gamma(\frac{s}2+2v)}{\Gamma(\frac{s+1}2)}\int_1^\infty\int_0^{1/u}\frac{y^{s-1}u^{s/2-1}}{\sqrt{(1-y^2)(1-y^2u^2)}}\hyF{2}{1}{\frac{s}2-2v,\frac{s}2+2v}{\frac{s+1}2}{-\frac{(1-u)^2}{4u}}dydu \end{align*}$$
二項目の積分に$t=uy$という変数変換を施します。
$$\begin{align*}&\frac{\sqrt{\pi}}{2^s}\frac{\Gamma(\frac{s}2)\Gamma(\frac{s}2-2v)\Gamma(\frac{s}2+2v)}{\Gamma(\frac{s+1}2)}\int_0^1\int_0^{1}\frac{y^{s-1}u^{s/2-1}}{\sqrt{(1-y^2)(1-y^2u^2)}}\hyF{2}{1}{\frac{s}2-2v,\frac{s}2+2v}{\frac{s+1}2}{-\frac{(1-u)^2}{4u}}dydu\\ &+\frac{\sqrt{\pi}}{2^s}\frac{\Gamma(\frac{s}2)\Gamma(\frac{s}2-2v)\Gamma(\frac{s}2+2v)}{\Gamma(\frac{s+1}2)}\int_1^\infty\int_0^{1}\frac{t^{s-1}u^{-s/2-1}}{\sqrt{(1-(t/u)^2)(1-t^2)}}\hyF{2}{1}{\frac{s}2-2v,\frac{s}2+2v}{\frac{s+1}2}{-\frac{(1-u)^2}{4u}}dtdu \end{align*}$$
Eulerの積分表示を使って各項の積分を超幾何関数に書き直します。
$$\begin{align*}&\frac{\pi}{2^{s+1}}\frac{\Gamma(\frac{s}2)^2\Gamma(\frac{s}2-2v)\Gamma(\frac{s}2+2v)}{\Gamma(\frac{s+1}2)^2}\int_0^1 u^{s/2-1}\hyF{2}{1}{\frac12,\frac{s}2}{\frac{s+1}2}{u^2}\hyF{2}{1}{\frac{s}2-2v,\frac{s}2+2v}{\frac{s+1}2}{-\frac{(1-u)^2}{4u}}du\\ &+\frac{\pi}{2^{s+1}}\frac{\Gamma(\frac{s}2)^2\Gamma(\frac{s}2-2v)\Gamma(\frac{s}2+2v)}{\Gamma(\frac{s+1}2)^2}\int_1^\infty u^{-s/2-1}\hyF{2}{1}{\frac12,\frac{s}2}{\frac{s+1}2}{\frac1{u^2}}\hyF{2}{1}{\frac{s}2-2v,\frac{s}2+2v}{\frac{s+1}2}{-\frac{(1-u)^2}{4u}}du \end{align*}$$
二項目に$u\to1/u$という変換を施せば二項は等しくなり、次のように一つにまとめられます。
$$\begin{align*}&\frac{\pi}{2^{s}}\frac{\Gamma(\frac{s}2)^2\Gamma(\frac{s}2-2v)\Gamma(\frac{s}2+2v)}{\Gamma(\frac{s+1}2)^2}\int_0^1 u^{s/2-1}\hyF{2}{1}{\frac12,\frac{s}2}{\frac{s+1}2}{u^2}\hyF{2}{1}{\frac{s}2-2v,\frac{s}2+2v}{\frac{s+1}2}{-\frac{(1-u)^2}{4u}}du \end{align*}$$

定理７において、一つ目の超幾何関数に補題６の変換公式を適応すると左辺は次に等しいことが分かります。
$$\frac{\pi}{2^{s}}\frac{\Gamma(\frac{s}2)^2\Gamma(\frac{s}2-2v)\Gamma(\frac{s}2+2v)}{\Gamma(\frac{s+1}2)^2}\int_0^1 u^{s/2-1}(1+u)^{-s}\hyF{2}{1}{\frac{s}2,\frac{s}2}{s}{\frac{4u}{(1+u)^2}}\hyF{2}{1}{\frac{s}2-2v,\frac{s}2+2v}{\frac{s+1}2}{-\frac{(1-u)^2}{4u}}du $$
ここで$u=(1-x)/(1+x)$と変数変換すると、次のように積分はかなりきれいになります。
$$\frac{\pi}{2^{2s-1}}\frac{\Gamma(\frac{s}2)^2\Gamma(\frac{s}2-2v)\Gamma(\frac{s}2+2v)}{\Gamma(\frac{s+1}2)^2}\int_0^1 (1-x^2)^{s/2-1}\hyF{2}{1}{\frac{s}2,\frac{s}2}{s}{1-x^2}\hyF{2}{1}{\frac{s}2-2v,\frac{s}2+2v}{\frac{s+1}2}{\frac{x^2}{x^2-1}}dx $$
二つ目の超幾何関数にPfaffの変換公式を適応して$x^2\to x$と変数変換すればこうなります。
$$\frac{\pi}{2^{2s}}\frac{\Gamma(\frac{s}2)^2\Gamma(\frac{s}2-2v)\Gamma(\frac{s}2+2v)}{\Gamma(\frac{s+1}2)^2}\int_0^1 x^{1/2-1}(1-x)^{s-2v-1-1}\hyF{2}{1}{\frac{s}2,\frac{s}2}{s}{1-x}\hyF{2}{1}{\frac{s}2-2v,\frac{1}2-2v}{\frac{s+1}2}{x}dx $$
二つの超幾何関数を級数に書き直してベータ積分を使えば次のように級数だけで表すことができます。
$$\frac{\pi\sqrt{\pi}}{2^{2s}}\Gamma\left[\begin{matrix}\frac{s}2,\frac{s}2,s-2v,\frac{s}2-2v,\frac{s}2+2v\\\frac{s+1}2,\frac{s+1}2,\frac12+s-2v\end{matrix}\right]\sum_{n=0}^\infty\sum_{m=0}^\infty \frac{(\frac{s}2,\frac{s}2,s-2v)_n(\frac{s}2-2v,\frac12-2v,\frac12)_m}{(s)_n(\frac{s+1}{2})_m(\frac12+s-2v)_{n+m}n!m!} $$
これは求めていた式ですね。

定理８は$(x)_{n+m}=(x)_n(x+n)_{m}$を使って次のように書くこともできます。
$$\int_0^\infty x^{s-1}K_v(x)^4dx=\frac{\pi\sqrt{\pi}}{2^{2s}}\Gamma\left[\begin{matrix}\frac{s}2,\frac{s}2,s-2v,\frac{s}2-2v,\frac{s}2+2v\\\frac{s+1}2,\frac{s+1}2,\frac12+s-2v\end{matrix}\right]\sum_{n=0}^\infty\frac{(\frac{s}2-2v,\frac12-2v,\frac12)_n}{(\frac{s+1}{2},\frac12+s-2v)_{n}n!}\hyF{3}{2}{\frac{s}2,\frac{s}2,s-2v}{s,\frac12+s-2v+n}{1}$$
$v=0$とすれば級数の中の${}_3F_2$が${}_2F_1$になり、Gaussの超幾何定理が使えます。
$$\begin{align*}\int_0^\infty x^{s-1}K_0(x)^4dx&=\frac{\pi\sqrt{\pi}}{2^{2s}}\Gamma\left[\begin{matrix}\frac{s}2,\frac{s}2,s,\frac{s}2,\frac{s}2\\\frac{s+1}2,\frac{s+1}2,\frac12+s\end{matrix}\right]\sum_{n=0}^\infty\frac{(\frac{s}2,\frac12,\frac12)_n}{(\frac{s+1}{2},\frac12+s)_{n}n!}\hyF{2}{1}{\frac{s}2,\frac{s}2}{\frac12+s}{1}\\ &=\frac{\pi^2}{2^{2s}}\Gamma\left[\begin{matrix}\frac{s}2,\frac{s}2,s,\frac{s}2,\frac{s}2\\\frac{s+1}2,\frac{s+1}2,\frac{s+1}2,\frac{s+1}2\end{matrix}\right]\sum_{n=0}^\infty\frac{(\frac{s}2,\frac12,\frac12,\frac12)_n}{(\frac{s+1}{2},\frac{s+1}{2},\frac{s+1}{2})_{n}n!}\\ &=\frac{\pi^2}{2^{2s}}\frac{\Gamma(s)\Gamma(\frac{s}2)^4}{\Gamma(\frac{s+1}{2})^4}{}_4F_{3}\Biggl[\begin{matrix}\frac{s}2,\frac12,\frac12,\frac12\\\frac{s+1}2,\frac{s+1}2,\frac{s+1}2\end{matrix};1\Biggl] \end{align*}$$
$v=1/4$では$(0)_n$が現れるので$n=0$以外の項が全て消え、次のように${}_3F_2$で書くことができます。
$$\begin{align*}\int_0^\infty x^{s-1}K_\frac14(x)^4dx&=\frac{\pi\sqrt{\pi}}{2^{2s}}\Gamma\left[\begin{matrix}\frac{s}2,\frac{s}2,s-\frac12,\frac{s-1}2,\frac{s+1}2\\\frac{s+1}2,\frac{s+1}2,s\end{matrix}\right]\sum_{n=0}^\infty\frac{(\frac{s-1}2,0,\frac12)_n}{(\frac{s+1}{2},s)_{n}n!}\hyF{3}{2}{\frac{s}2,\frac{s}2,s-\frac12}{s,s+n}{1}\\ &=\frac{\pi\sqrt{\pi}}{2^{2s}}\Gamma\left[\begin{matrix}\frac{s}2,\frac{s}2,s-\frac12,\frac{s-1}2\\\frac{s+1}2,s\end{matrix}\right]\hyF{3}{2}{\frac{s}2,\frac{s}2,s-\frac12}{s,s}{1}\\ \end{align*}$$
この${}_3F_2$にThomaeの変換公式を適応すれば求めていた式が得られます。
\begin{align*} \int_0^\infty x^{s-1}K_\frac14(x)^4dx=\frac{2\pi^2}{s-1}\frac{\Gamma(s)\Gamma(\frac{s}2)^2}{2^{2s}\Gamma(\frac{s+1}{2})^2}{}_3F_{2}\Biggl[\begin{matrix}\frac12,\frac12,\frac12\\\frac{s+1}2,\frac{s+1}2\end{matrix};1\Biggl] \end{align*}
この手法から、$v$が整数の時は${}_4F_3$の有限和、$v=\frac14+\frac{n}2$の時は${}_3F_2$の有限和で書けることもわかりますね。

定理７より次のように書けます。
$$\int_0^\infty xK_v(x)^4dx=\frac{2\pi v}{\sin(2\pi v)}\int_0^1 \hyF{2}{1}{\frac12,1}{\frac{3}2}{u^2}\hyF{2}{1}{1-2v,1+2v}{\frac{3}2}{-\frac{(1-u)^2}{4u}}du$$
ここで、このパラメータに関する超幾何関数は次のように書けることが知られています。
$$\hyF{2}{1}{\frac12,1}{\frac{3}2}{u^2}=\frac{\tanh^{-1}(u)}{u}\:\:,\:\:\hyF{2}{1}{1-2v,1+2v}{\frac{3}2}{-\frac{(1-u)^2}{4u}}=\frac1v\frac{u}{u^2-1}\sinh(2v\ln(u))$$
よってこの積分はこのように書けます。
$$\begin{align*}\int_0^\infty xK_v(x)^4dx&=\frac{\pi}{\sin(2\pi v)}\int_0^1 \frac{\tanh^{-1}(u)}{u^2-1}(u^{2v}-u^{-2v})du\\ &=\frac{\pi}{2\sin(2\pi v)}\sum_{n=0}^\infty\frac{1}{2n+1}\int_0^1\frac{u^{n+v}-u^{n-v}}{u-1}du\\ &=\frac{\pi}{2\sin(2\pi v)}\sum_{n=0}^\infty\frac{\psi(n+1+v)-\psi(n+1-v)}{2n+1}\\ &=\frac{\pi^2}{16\sin(\pi v)^2}\left(2\psi\left(\frac12\right)-\psi\left(\frac12+v\right)-\psi\left(\frac12-v\right)\right) \end{align*}$$
私には最後の変形ができませんでした…