When deriving the angular distribution of energy for synchrotron radiation one has to evaluate two tricky improper integrals (see [1] below): $$ I_1 \equiv \int_{0}^{\infty} x^2 [K_{2/3}(x)]^2 \, \text{d}x, \quad I_2 \equiv \int_{0}^{\infty} x^2 [K_{1/3}(x)]^2 \, \text{d}x, $$ where $K_{2/3}(x)$ and $K_{1/3}(x)$ are modified Bessel functions of the second kind. Mathematica yields $I_1 = \frac{7\pi^2}{144}$ and $I_2 = \frac{5\pi^2}{144}$, but how can I derive these results? Starting from the general integral $I(s,\nu) = \int_{0}^{\infty} x^{s-1} [K_{\nu}(x)]^2 \,\text{d}x$, I tried using the following integral representation of $K_{\nu}(x)$: $$ K_{\nu}(x) = \int_{0}^{\infty} e^{-x \cosh t} \cosh(\nu t) \, \text{d}t. $$ In doing so, I found $$ I(s,\nu) = \Gamma(s) \int_{0}^{\infty} \int_{0}^{\infty} \frac{\cosh(\nu t)\cosh(\nu u)}{(\cosh t + \cosh u)^s} \, \text{d}t \, \text{d}u. $$ Jackson offers no hint in how one can evaluate these integrals, so any guidance will be valuable. Maybe some obscure identity/integral from Watson's book on Bessel functions should help, but I haven't found anything immediately helpful while skimming it.
For the last integral above, Mathematica returns $$ I(s,\nu) = \frac{(1-4\nu^2)\pi^2}{32} \, \sec(\pi \nu), \quad -\frac{3}{2} < \text{Re}(\nu) < \frac{3}{2}, $$ which hints at some Gamma/Beta function gymnastics if my intuition doesn't fail me.
[1] Jackson - Classical Electrodynamics, 3rd edition (1999), equations (14.79) and (14.80).