Generalized integrals for Bessel Moments $\int_{0}^{\infty} x^4K_0(x)K_1(x)^3 \ln(xK_1(x))^2\text{d}x=\frac{1}{32}$ Let $I_\nu(x)$ be the modified Bessel functions of first kind with order$\text{ }\nu$,
$K_\nu(x)$ be the modified Bessel functions of second kind with order$\text{ }\nu$.

Prerequisite Information:
The integral
$$\int_{0}^{\infty} x^4K_0(x)K_1(x)^3
\ln(xK_1(x))^2\text{d}x=\frac{1}{32}$$
can be shown as follows:

Note that $$\frac{\text{d}}{\mathrm{d}x}\left ( -\frac{x^\alpha}{\alpha}
K_1(x)^\alpha  \right ) 
=x^\alpha K_0(x)K_1(x)^{\alpha-1}.$$
Therefore,$$\int_{0}^{\infty}x^\alpha K_0(x)K_1(x)^{\alpha-1}\text{d}x
=\frac{1}{\alpha},\qquad{\Re(\alpha)>0} .$$
And the equality immediately follows by differentiating the expression.

There are also some integral identities involving Bessel functions, but not (quite) trivial. These integrals had studied in arXiv:0801.0891. For example,
$$\begin{aligned}
&\int_{0}^{\infty}K_0(x)^3\text{d}x=\frac{3\Gamma\left ( \frac{1}{3}  \right )^6 }{32\pi\cdot2^{2/3}}  ,\\
&\int_{0}^{\infty}xK_0(x)^4\text{d}x=\frac{7}{8}\zeta(3) ,\\
&\int_{0}^{\infty}xI_0(x)K_0(x)^2\text{d}x= \frac{\pi}{3\sqrt{3} },\\
&\int_{0}^{\infty}xI_0(x)K_0(x)^3\text{d}x=\frac{\pi^2}{16} .
\end{aligned}$$
In this paper, the authors determine some relations among the moments. For example,
$$
\int_{0}^{\infty}K_0(x)^4\text{d}x
=\pi^2\int_{0}^{\infty}K_0(x)^2I_0(x)^2\mathrm{d}x.
$$
These relations can be generalized in many ways. Using contour integration, we conclude that
$$
\int_{0}^{\infty} 
x^3 K_0(x)^5I_0(x)\left ( \pi^2I_0(x)^2-K_0(x)^2 \right ) 
\text{d}x=\frac{\pi^4}{128}.
$$
(Only one example.)
Moreover,
$$
\int_{0}^{\infty} 
x^{2k+1} K_0(x)^5I_0(x)\left ( \pi^2I_0(x)^2-K_0(x)^2 \right ) 
\text{d}x=
\begin{cases}
  0 & k=0, \\
 a_k\cdot\pi^4 & k\in\mathbb{Z}^{+}.
\end{cases}
$$
Where $a_k$ is always a rational number.
And we are able to compute
$$
\int_{0}^{\infty} 
xI_0(\alpha x) K_0(x)^5I_0(x)\left ( \pi^2I_0(x)^2-K_0(x)^2 \right ) 
\text{d}x
$$
by expanding the $I_0(\alpha x)$ into Maclaurin series.
Another simple identity is given by
$$
\int_{0}^{\infty}x^7K_0(x)K_1(x)^2K_2(x)\text{d}x
=\frac{1}{3}.
$$

Problem:
I am trying to find more results but failed. Can we find the closed-forms of other moments such as $\int_{0}^{\infty}K_0(x)^5\text{d}x,
\int_{0}^{\infty}K_0(x)I_0(x)J_0(x)Y_0(x)\text{d}x$? Any idea would be much appreciated.

Maybe interests:
Two integrals (both are easy to check):
$$\begin{aligned}
&\int_{0}^{\infty} \frac{x^2}{\alpha^2+x^2}K_0(x)^2\text{d}x
=\frac{\pi^2}{4}-\frac{\pi^3}{8}\alpha \left ( J_0(\alpha)^2+Y_0(\alpha)^2 \right ),
\\ 
&\int_{0}^{\infty}K_0(x)^2\cos(\alpha x)\mathrm{d}x
=\frac{\pi}{\sqrt{4+\alpha^2} }K\left ( \frac{\alpha}{\sqrt{4+\alpha^2} }  \right ). 
\end{aligned}$$
Where $K(x)=\frac\pi2{}_2F_1\left(\frac12,\frac12;1;x^2\right)$ and ${}_2F_1$ is Gauss hypergeometric function.
 A: In addition to the identity $$\int_{0}^{\infty} K_{0}(x)^{4} \, \mathrm dx = \pi^{2} \int_{0}^{\infty} K_{0}(x)^{2} I_{0}(x)^{2} \, \mathrm dx,$$ there is also the  identity $$\int_{0}^{\infty} J_{\alpha}(x)^{4} \, \mathrm dx = \int_{0}^{\infty} Y_{\alpha}(x)^{2}J_{\alpha}(x)^{2} \, \mathrm dx, \quad \alpha \ge 0. \tag{1}$$
I assume $(1)$ is a known identity.
The following is a way to proof $(1)$ using contour integration.

In the right half-plane, the function $$f(z) = K_{\alpha}(z)^{2} I_{\alpha}(z)^{2} $$ is $\mathcal{O} \left( \frac{1}{z^{2}} \right)$ as $|z| \to \infty$. (See here.)
The function $f(z)$ is also analytic in right half-plane.
So by integrating $f(z)$ around a closed quarter-circle in the first quadrant (that is indented at the origin), we get
$$\int_{0}^{\infty} K_{\alpha}(x)^{2}I_{\alpha}(x)^{2} \, \mathrm dx - \int_{0}^{\infty}K_{\alpha}(ix)^{2} I_{\alpha}(ix)^{2} \, i \, \mathrm dx = 0. $$
But  $$ \begin{align} I_{\alpha}(ix) &= \sum_{m=0}^{\infty} \frac{1}{m!\Gamma(m+\alpha +1)} \left(\frac{ix}{2} \right)^{2m+ \alpha} \\ &= e^{i  \alpha \pi  /2 } \sum_{m=0}^{\infty} \frac{(-1)^{m}}{m!\Gamma(m+\alpha +1)} \left(\frac{x}{2} \right)^{2m+ \alpha} \\ &= e^{i \alpha \pi /2}J_{\alpha}(x), \end{align}$$
and $$ \begin{align} K_{\alpha}(ix) &= \frac{\pi}{2} \frac{I_{-\alpha}(ix) - I_{\alpha}(ix)}{\sin(\alpha \pi)} \\ &= \frac{\pi}{2} \frac{e^{-i \alpha \pi /2} J_{-\alpha}(x) -  e^{i \alpha \pi /2} J_{\alpha}(x)}{\sin(\alpha \pi)} \\ & =\frac{\pi e^{- i  \alpha \pi /2} }{2}   \frac{J_{-\alpha}(x) -  e^{i  \alpha \pi } J_{\alpha}(x)}{\sin(\alpha \pi)} \\ &=  \frac{\pi e^{- i \alpha \pi /2} }{2}   \frac{J_{-\alpha}(x) -\frac{1}{2} \left(e^{i \alpha \pi}+e^{- i \alpha \pi }\right)J_{\alpha}(x)- \frac{1}{2}\left( e^{i \alpha \pi } -e^{-i \alpha \pi } \right)J_{\alpha}(x)}{\sin(\alpha \pi) } \\ &= \frac{\pi e^{- i \alpha \pi /2} }{2}   \frac{J_{-\alpha}(x) -\cos(\alpha \pi) J_{\alpha}(x)- i \sin(\alpha \pi) J_{\alpha}(x)}{\sin(\alpha \pi) } \\ &= -\frac{\pi e^{- i \alpha \pi /2} }{2} \left( Y_{\alpha}(x)+ i J_{\alpha}(x) \right). \end{align}$$
Therefore, we have $$\int_{0}^{\infty} K_{\alpha}(x)^{2}I_{\alpha}(x)^{2} \, \mathrm dx -  \frac{i \pi^{2}}{4} \int_{0}^{\infty} \left(Y_{\alpha}(x)+iJ_{\alpha}(x) \right)^{2} J_{\alpha}(x)^{2} \, \mathrm dx = 0. $$
And equating the imaginary parts on both sides of the equation, we get $$\int_{0}^{\infty} \left( Y_{\alpha}(x)^{2}J_{\alpha}(x)^{2}- J_{\alpha}(x)^{4} \right) \, \mathrm dx =0. $$

Also, by equating the real parts on both sides of the equation, we get $$\int_{0}^{\infty} K_{\alpha}(x)^{2} I_{\alpha}(x)^{2} \, \mathrm dx + \frac{\pi^{2}}{2} \int_{0}^{\infty} Y_{\alpha}(x) J_{\alpha}(x)^{3} \, \mathrm dx =0. $$
