A solution in large steps for:
$$I=\int _0^{\frac{\pi }{2}}x\ln \left(\sin \left(x\right)\right)\ln ^2\left(\cos \left(x\right)\right)\:dx=\frac{155}{128}\zeta \left(5\right)+\frac{13}{32}\zeta \left(2\right)\zeta \left(3\right)-\operatorname{Li}_5\left(\frac{1}{2}\right)$$
$$-\frac{49}{32}\ln \left(2\right)\zeta \left(4\right)-\frac{5}{6}\ln ^3\left(2\right)\zeta \left(2\right)+\frac{1}{120}\ln ^5\left(2\right),$$
$$J=\int _0^{\frac{\pi }{2}}x\ln ^2\left(\sin \left(x\right)\right)\ln \left(\cos \left(x\right)\right)\:dx=-\frac{155}{128}\zeta \left(5\right)-\frac{1}{32}\zeta \left(2\right)\zeta \left(3\right)+\operatorname{Li}_5\left(\frac{1}{2}\right)$$
$$+\frac{49}{32}\ln \left(2\right)\zeta \left(4\right)-\frac{2}{3}\ln ^3\left(2\right)\zeta \left(2\right)-\frac{1}{120}\ln ^5\left(2\right).$$
Let's consider calculating $I-J$ first.
$$I-J=\int _0^{\frac{\pi }{2}}x\ln \left(\sin \left(x\right)\right)\ln ^2\left(\cos \left(x\right)\right)\:dx-\int _0^{\frac{\pi }{2}}x\ln ^2\left(\sin \left(x\right)\right)\ln \left(\cos \left(x\right)\right)\:dx$$
$$=\frac{1}{3}\underbrace{\int _0^{\infty }\frac{\ln ^3\left(x\right)\arctan \left(x\right)}{1+x^2}\:dx}_{\mathcal{A}}-\frac{2}{3}\underbrace{\int _0^{\frac{\pi }{2}}x\ln ^3\left(\sin \left(x\right)\right)\:dx}_{\mathcal{B}}+\frac{\pi }{6}\underbrace{\int _0^{\frac{\pi }{2}}\ln ^3\left(\sin \left(x\right)\right)\:dx}_{\mathcal{C}}.$$
Here $\mathcal{A}$ and $\mathcal{C}$ are very straightforward, the former can be easily calculated through differentiation under the integral sign considering the following parameter:
$$\mathcal{A}\left(a\right)=\int _0^{\infty }\frac{\ln ^3\left(x\right)\arctan \left(ax\right)}{1+x^2}\:dx.$$
After the usual procedure it reduces to trivial/well-known integrals, in the other hand $\mathcal{C}$ may be calculated employing the Beta function, so:
$$\mathcal{A}=\int _0^{\infty }\frac{\ln ^3\left(x\right)\arctan \left(x\right)}{1+x^2}\:dx=\frac{93}{16}\zeta \left(5\right)+\frac{21}{16}\zeta \left(2\right)\zeta \left(3\right),$$
$$\mathcal{C}=\int _0^{\frac{\pi }{2}}\ln ^3\left(\sin \left(x\right)\right)\:dx=-\frac{3\pi }{4}\zeta \left(3\right)-\frac{\pi ^3}{8}\ln \left(2\right)-\frac{\pi }{2}\ln ^3\left(2\right).$$
This only leaves us with the toughest one which is:
$$\mathcal{B}=\int _0^{\frac{\pi }{2}}x\ln ^3\left(\sin \left(x\right)\right)\:dx=-\frac{\pi }{16}\underbrace{\int _0^{\infty }\frac{\ln ^3\left(1+x^2\right)}{1+x^2}\:dx}_{\operatorname{Beta function}}-\frac{1}{8}\int _0^{\infty }\frac{\arctan \left(\frac{1}{x}\right)\ln ^3\left(\frac{x^2}{1+x^2}\right)}{1+x^2}\:dx$$
$$=-\frac{9}{4}\zeta \left(2\right)\zeta \left(3\right)-\frac{45}{8}\ln \left(2\right)\zeta \left(4\right)-\frac{3}{2}\ln ^3\left(2\right)\zeta \left(2\right)-\frac{1}{4}\operatorname{\mathfrak{I}} \left\{\int _0^{\infty }\frac{\ln ^4\left(\frac{x}{x-i}\right)}{1+x^2}\:dx\right\}$$
$$+\frac{1}{2}\int _0^{\infty }\frac{\arctan ^3\left(x\right)\ln \left(1+x^2\right)}{1+x^2}\:dx$$
$$=-\frac{9}{4}\zeta \left(2\right)\zeta \left(3\right)-\frac{45}{8}\ln \left(2\right)\zeta \left(4\right)-\frac{3}{2}\ln ^3\left(2\right)\zeta \left(2\right)$$
$$-\frac{1}{4}\operatorname{\mathfrak{I}} \left\{-12i\operatorname{Li}_5\left(2\right)+\pi \ln ^4\left(2\right)\right\}-\int _0^{\frac{\pi }{2}}x^3\ln \left(\cos \left(x\right)\right)\:dx.$$
Where the fourier series of $\ln \left(\cos \left(x\right)\right)$ may be employed for that last integral for a fast calculation.
In the end we obtain that:
$$\mathcal{B}=\int _0^{\frac{\pi }{2}}x\ln ^3\left(\sin \left(x\right)\right)\:dx=-\frac{93}{128}\zeta \left(5\right)-\frac{9}{8}\zeta \left(2\right)\zeta \left(3\right)+3\operatorname{Li}_5\left(\frac{1}{2}\right)$$
$$+\frac{57}{32}\ln \left(2\right)\zeta \left(4\right)-\frac{1}{2}\ln ^3\left(2\right)\zeta \left(2\right)-\frac{1}{40}\ln ^5\left(2\right).$$
Collecting all the previous results leads us to the closed-form of $I-J$ which is:
$$I-J=\int _0^{\frac{\pi }{2}}x\ln \left(\sin \left(x\right)\right)\ln ^2\left(\cos \left(x\right)\right)\:dx-\int _0^{\frac{\pi }{2}}x\ln ^2\left(\sin \left(x\right)\right)\ln \left(\cos \left(x\right)\right)\:dx$$
$$=\frac{155}{64}\zeta \left(5\right)+\frac{7}{16}\zeta \left(2\right)\zeta \left(3\right)-2\operatorname{Li}_5\left(\frac{1}{2}\right)-\frac{49}{16}\ln \left(2\right)\zeta \left(4\right)-\frac{1}{6}\ln ^3\left(2\right)\zeta \left(2\right)+\frac{1}{60}\ln ^5\left(2\right).$$
Note that if we apply the substitution $x=\frac{\pi }{2}-t$ to either $I$ or $J$ we can easily find each of them separately.
And that completes the solution.
The following answer shows that $$\small \left[ K_{ib}(a) \right]^{2} + \frac{\pi}{2} \, \cosh(\pi b) \int_{0}^{\infty} Y_{0}\left(2 a \sinh \frac{t}{2} \right) \cos(bt) \, \mathrm dt + \frac{\pi}{2} \, \sinh(\pi b) \int_{0}^{\infty} J_{0}\left(2 a \sinh \frac{t}{2} \right) \sin(bt) \, \mathrm dt =0. $$
The two integral identities satisfy this equation since $- \cosh^{2}(\pi b)+ \sinh^{2}(\pi b) = -1$.
This equation can be used with Gary's answer to prove the second integral identity.
My starting point is the more well known integral representation $$\left[ K_{\nu}(x) \right]^{2} = 2 \int_{0}^{\infty}K_{0}(2x \cosh t) \cosh(2 \nu t) \, \mathrm dt , \, \quad x>0. $$
Using this representation, we have $$\left[ K_{ib}(a) \right]^{2} = 2 \int_{0}^{\infty} K_{0} (2a \cosh x) \cos(2bx) \, \mathrm dx = \int_{0}^{\infty} K_{0} \left(2a \cosh \frac{u}{2} \right) \cos(bu) \, \mathrm du.$$
The important thing to notice here is that $\left[ K_{ib}(a) \right]^{2}$ is a real value if $a$ and $b$ are positive.
The function $$ f(z) = K_{0} \left(2a \cosh \frac{z}{2} \right) \cos(bz) $$ has branch cuts were $\cosh \left(\frac{x}{2} \right)$ is real and negative.
None of these branch cuts fall inside or on a rectangular contour with vertices at $z=0$, $z= R$, $z= R +\pi i $, and $z= \pi i $.
There is, however, a branch point on the contour at $z= i \pi$.
So let's integrate $f(z)$ counterclockwise around the above contour with the addition of an quarter circle indentation about $z = i \pi$.
However, since $ \lim_{z \to i \pi} (z- i \pi)K_{0} \left(2a \cosh \frac{z}{2} \right) \cos(bz) =0$, there is no contribution from letting the radius of the indentation go to zero. (The function $K_{0}\left( 2a \cosh \frac{z}{2} \right)$ behaves like $-\log(z-i \pi)$ near $z= i \pi$.)
As $R \to \infty$, the integral vanishes on the right side of contour because the magnitude of $K_{0}(z)$ decays exponentially to zero as $\Re(z) \to + \infty$.
On the left side of the contour, we have $$-i \int_{0}^{\pi} K_{0} \left(2a \cos \frac{t}{2} \right) \cosh (bt) \, \mathrm dt. $$
And on the top side of the contour, we have $$ \begin{align} &-\int_{0}^{\infty} K_{0} \left (2a \cosh \frac{t+ i \pi}{2} \right) \cos \left(b(t + i \pi) \right) \, \mathrm dt \\ &= - \int_{0}^{\infty} K_{0} \left(2ai \sinh \frac{t}{2} \right) \left(\cosh(\pi b)\cos(bt) - i \sinh(\pi b) \sin(bt) \right) \, \mathrm dt \\ & \overset{(1)}{=} \frac{\pi}{2} \int_{0}^{\infty}\left(Y_{0} \left(2a \sinh \frac{t}{2} \right) + i J_{0} \left(2a \sinh \frac{t}{2} \right) \right)\left(\cosh(\pi b)\cos(bt) - i \sinh(\pi b) \sin(bt) \right) \, \mathrm dt. \end{align}$$
Therefore, since there are no singularities inside the contour, we have $$\left[ K_{ib}(a) \right]^{2} + \frac{\pi}{2} \int_{0}^{\infty}\left(Y_{0} \left(2a \sinh \frac{t}{2} \right) + i J_{0} \left(2a \sinh \frac{t}{2} \right) \right)\left(\cosh(\pi b)\cos(bt) - i \sinh(\pi b) \sin(bt) \right) \, \mathrm dt$$
$$-i \int_{0}^{\pi} K_{0} \left(2a \cos \frac{t}{2} \right) \cosh (bt) \, \mathrm dt =0. $$
And by equating the real parts on both sides of the above equation, we get $$ \small \left[ K_{ib}(a) \right]^{2} + \frac{\pi}{2} \, \cosh(\pi b) \int_{0}^{\infty} Y_{0}\left(2 a \sinh \frac{t}{2} \right) \cos(bt) \, \mathrm dt + \frac{\pi}{2} \, \sinh(\pi b) \int_{0}^{\infty} J_{0}\left(2 a \sinh \frac{t}{2} \right) \sin(bt) \, \mathrm dt =0. $$
$(1)$ See the answer here.
Best Answer
I believe that your conjecture is indeed true, by utilizing similar technique in the linked integral.
To begin with, we have $$ \mathrm{arctan} \left( \cot \left( \theta x \right) \right) =\pi \left( \left\{ -\frac{\theta x}{\pi} \right\} -\frac{1}{2} \right) $$ Now, by using the fourier expansion of the fractional part from $\left[-\pi,\pi\right]$. We may exclude the points that would be problematic for this expansion ($\left\{ x\in \mathbb{R} \mid \theta x\in \pi \mathbb{Z} \right\} $), since it's a measure-zero set and will not affect the evaluation of the integral. So $$ \pi \left( \left\{ -\frac{\theta x}{\pi} \right\} -\frac{1}{2} \right) =\pi \left( \frac{1}{2}+\frac{1}{\pi}\sum_{k=1}^{\infty}{\frac{\sin \left( 2\theta kx \right)}{k}}-\frac{1}{2} \right) =\sum_{k=1}^{\infty}{\frac{\sin \left( 2\theta kx \right)}{k}} $$ Pulgging this into the intergal, we have \begin{align*} &\int_0^{\pi}{\mathrm{arctan} \left( \cot \left( mx \right) \right) \mathrm{arctan} \left( \cot \left( nx \right) \right)}\mathrm{d}x \\ =&\int_0^{\pi}{\sum_{i=1}^{\infty}{\frac{\sin \left( 2mix \right)}{i}}\sum_{j=1}^{\infty}{\frac{\sin \left( 2njx \right)}{j}}}\mathrm{d}x \\ =&\sum_{i=1}^{\infty}{\sum_{j=1}^{\infty}{\frac{1}{ij}}}\int_0^{\pi}{\sin \left( 2mix \right) \sin \left( 2njx \right)}\mathrm{d}x \\ =&\frac{\pi}{2}\sum_{i=1}^{\infty}{\sum_{j=1}^{\infty}{\frac{1}{ij}}}\delta _{mi,nj} \end{align*} You would encounter the exact same sum when you try the same method on that linked integral. Let's evaluate it. first, we have $$ \forall i,j\in \mathbb{N} , a=ni=mj\Rightarrow \left[ m,n \right] |a\Rightarrow \frac{1}{i}=\frac{n}{a}, \frac{1}{j}=\frac{m}{a} $$ So, we can rewrite the sum in terms of $\left[ m,n \right]$. $$ \frac{\pi}{2}\sum_{i=1}^{\infty}{\sum_{j=1}^{\infty}{\begin{array}{c} \frac{1}{ij}\delta _{ni,mj}\\ \end{array}}}=\frac{\pi}{2}\sum_{k=1}^{\infty}{\frac{nm}{k^2\left[ m,n \right] ^2}}=\frac{\pi \left( m,n \right)}{2\left[ m,n \right]}\sum_{k=1}^{\infty}{\frac{1}{k^2}}=\frac{\pi ^3\left( m,n \right)}{12\left[ m,n \right]} $$ Which is your claim. At the end, we have $$ I\left( n,m \right) =\int_0^{\pi}{\mathrm{arctan} \left( \cot \left( mx \right) \right) \mathrm{arctan} \left( \cot \left( nx \right) \right)}\mathrm{d}x=\frac{\pi ^3\left( m,n \right)}{12\left[ m,n \right]} $$