Okay, finally I was able to prove it.
Step 0. Observations. In view of the following identity
$$ \int_{0}^{\frac{\pi}{2}} \arctan (r \sin\theta) \, d\theta = 2 \chi_{2} \left( \frac{\sqrt{1+r^{2}} - 1}{r} \right), $$
Vladimir's result suggests that there may exists a general formula connecting
$$ I(r, s) = \int_{0}^{\frac{\pi}{2}} \arctan (r \sin\theta) \arctan (s \sin\theta) \, d\theta $$
and the Legendre chi function $\chi_{2}$. Indeed, inspired by Vladimir's result, I conjectured that
$$ I(r, s) = \pi \chi_{2} \left( \frac{\sqrt{1+r^{2}} - 1}{r} \cdot \frac{\sqrt{1+s^{2}} - 1}{s} \right). \tag{1} $$
I succeeded in proving this identity, so I post a solution here.
Step 1. Proof of the identity $\text{(1)}$. It is easy to check that the following identity holds:
$$ \arctan(ab) = \int_{1/b}^{\infty} \frac{a \, dx}{a^{2} + x^{2}}. $$
So it follows that
\begin{align*}
I(r, s)
&= \int_{1/r}^{\infty} \int_{1/s}^{\infty} \int_{0}^{\frac{\pi}{2}} \frac{\sin^{2}\theta}{(x^{2} + \sin^{2}\theta)(y^{2} + \sin^{2}\theta)} \, d\theta dy dx \\
&= \int_{1/r}^{\infty} \int_{1/s}^{\infty} \int_{0}^{\frac{\pi}{2}} \frac{1}{x^{2} - y^{2}} \left( \frac{x^{2}}{x^{2} + \sin^{2}\theta} - \frac{y^{2}}{y^{2} + \sin^{2}\theta} \right) \, d\theta dy dx \\
&= \frac{\pi}{2} \int_{1/r}^{\infty} \int_{1/s}^{\infty} \frac{1}{x^{2} - y^{2}} \left( \frac{x}{\sqrt{x^{2} + 1}} - \frac{y}{\sqrt{y^{2} + 1}} \right) \, dy dx.
\end{align*}
For the convenience of notation, we put
$$ \alpha = \frac{\sqrt{r^{2} + 1} - 1}{r} \quad \text{and} \quad \beta = \frac{\sqrt{s^{2} + 1} - 1}{s}. $$
Then it is easy to check that $\mathrm{arsinh}(1/r) = - \log \alpha$ and likewise for $s$ and $\beta$. Thus with the substitution $x \mapsto \sinh x$ and $y \mapsto \sinh y$, we have
\begin{align*}
I(r, s)
&= \frac{\pi}{2} \int_{-\log\alpha}^{\infty} \int_{-\log\beta}^{\infty} \frac{\sinh x \cosh y - \sinh y \cosh x}{\sinh^{2}x - \sinh^{2}y} \, dy dx.
\end{align*}
Applying the substitution $e^{-x} \mapsto x$ and $e^{-y} \mapsto y$, it follows that
\begin{align*}
I(r, s)
&= \pi \int_{0}^{\alpha} \int_{0}^{\beta} \frac{dydx}{1 - x^{2}y^{2}} \\
&= \pi \sum_{n=0}^{\infty} \left( \int_{0}^{\alpha} x^{2n} \, dx \right) \left( \int_{0}^{\beta} y^{2n} \, dx \right)
= \pi \sum_{0}^{\infty} \frac{(\alpha \beta)^{2n+1}}{(2n+1)^{2}} \\
&= \pi \chi_{2}(\alpha \beta)
\end{align*}
as desired, proving the identity $\text{(1)}$.
EDIT. I found a much simpler and intuitive proof of $\text{(1)}$. We first observe that $\text{(1)}$ is equivalent to the following identity
$$ \int_{0}^{\frac{\pi}{2}} \arctan\left( \frac{2r\sin\theta}{1-r^{2}} \right) \arctan\left( \frac{2s\sin\theta}{1-s^{2}} \right) \, d\theta = \pi \chi_{2}(rs). $$
Now we first observe that from the addition formula for the hyperbolic tangent, we obtain the following formula
$$
\operatorname{artanh}x - \operatorname{artanh} y = \operatorname{artanh} \left( \frac{x - y}{1 - xy} \right) $$
which holds for sufficiently small $x, y$. Thus
\begin{align*}
\arctan\left( \frac{2r\sin\theta}{1-r^{2}} \right)
&= \frac{1}{i} \operatorname{artanh}\left( \frac{2ir\sin\theta}{1-r^{2}} \right)
= \frac{\operatorname{artanh}(re^{i\theta}) - \operatorname{artanh}(re^{-i\theta})}{i} \\
&= 2 \Im \operatorname{artanh}(re^{i\theta})
= 2 \sum_{n=0}^{\infty} \frac{\sin(2n+1)\theta}{2n+1} r^{2n+1}.
\end{align*}
We readily check this holds for any $|r| < 1$. Therefore
\begin{align*}
&\int_{0}^{\frac{\pi}{2}} \arctan\left( \frac{2r\sin\theta}{1-r^{2}} \right) \arctan\left( \frac{2s\sin\theta}{1-s^{2}} \right) \, d\theta \\
&\quad = 4 \sum_{m=0}^{\infty} \sum_{n=0}^{\infty} \frac{r^{2m+1}s^{2n+1}}{(2m+1)(2n+1)} \int_{0}^{\frac{\pi}{2}} \sin(2m+1)\theta \sin(2n+1)\theta \, d\theta\\
&\quad = 2 \sum_{m=0}^{\infty} \sum_{n=0}^{\infty} \frac{r^{2m+1}s^{2n+1}}{(2m+1)(2n+1)} \int_{0}^{\frac{\pi}{2}} \{ \cos(2m-2n)\theta - \cos(2m+2n+2)\theta \} \, d\theta\\
&\quad = \pi \sum_{m=0}^{\infty} \sum_{n=0}^{\infty} \frac{r^{2m+1}s^{2n+1}}{(2m+1)(2n+1)} \delta_{m,n} \\
&\quad = \pi \chi_{2}(rs).
\end{align*}
The substitution $\sin(x) = \sqrt{t}$ leads to the expression
$$I = \frac{1}{2} \int \limits_0^1 \frac{t^{-1/4} \arcsin^2 (\sqrt{t})}{\sqrt{1-t}} \, \mathrm{d} t \, . $$
Now you can use the power series for $\arcsin^2$ (see for example this question) and integrate term by term (monotone convergence). Using the beta function you will find
\begin{align}
I &= \frac{1}{4} \sum \limits_{n=1}^\infty \frac{(2n)!!}{n^2 (2n-1)!!} \int \limits_0^1 t^{n-\frac{1}{4}} (1-t)^{-\frac{1}{2}} \, \mathrm{d} t \\
&= \frac{1}{4} \sum \limits_{n=1}^\infty \frac{(2n)!!}{n^2 (2n-1)!!} \operatorname{B}\left(n+\frac{3}{4},\frac{1}{2}\right) \\
&= \frac{\sqrt{\pi}}{4} \sum \limits_{n=1}^\infty \frac{(2n)!!}{n^2 (2n-1)!!} \frac{\Gamma\left(n+\frac{3}{4}\right)}{\Gamma\left(n+\frac{5}{4}\right)} \\
&= \frac{\sqrt{\pi} \, \Gamma\left(\frac{3}{4}\right)}{\Gamma\left(\frac{1}{4}\right)} \sum \limits_{n=1}^\infty \frac{(2n)!!}{n^2 (2n-1)!!} \frac{\prod_{k=1}^n (4k-1)}{\prod_{l=1}^{n+1} (4l-3)} \\
&= \frac{\pi \sqrt{2 \pi}}{\Gamma\left(\frac{1}{4}\right)^2} \sum \limits_{n=1}^\infty \frac{(2n)!!}{n^2 (4n+1) (2n-1)!!} \prod \limits_{k=1}^n \frac{4k-1}{4k-3} \, .
\end{align}
Mathematica gives the following expression in terms of a hypergeometric function:
$$ I = \frac{6 \pi \sqrt{2 \pi}}{5 \Gamma\left(\frac{1}{4}\right)^2} \, {}_4 \! \operatorname{F}_3 \left(1,1,1,\frac{7}{4};\frac{3}{2},2,\frac{9}{4};1\right) \approx 1.208656578687 \, .$$
Inverse symbolic calculators do not seem to give any expression for this number, so this might be as good as it gets.
Best Answer
$$I=\int_0^\frac{\pi}{2} x\cot x \ln(\cos x)dx\overset{IBP}=\int_0^\frac{\pi}{2}x\tan x\ln(\sin x)dx-\int_0^\frac{\pi}{2}\ln(\sin x)\ln(\cos x)dx$$ The substitution $\frac{\pi}{2}-x=x$ in the first integral gives: $$ I=\frac{\pi}{2} \int_0^\frac{\pi}{2}\cot x\ln(\cos x)dx-I-\int_0^\frac{\pi}{2}\ln(\sin x)\ln(\cos x)dx$$ $$I=\frac{\pi}{4} \int_0^\frac{\pi}{2}\cot x\ln(\cos x)dx -\frac\pi4 \ln^22+\frac{\pi^3}{96}$$ I won't focus on the second integral since I believe there is a way to avoid all the calculation and magically simplify it, but here is an approach. $$J=\int_0^\frac{\pi}{2}\cot x\ln(\cos x)dx\overset{\tan x=t}=-\frac12 \int_0^\infty \frac{\ln(1+x^2)}{x(1+x^2)}dx$$ Split the integral in the point $1$ then let $\frac{1}{x}\to x$ in the second part. $$J=-\frac12 \int_0^1 \frac{\ln(1+x^2)}{x(1+x^2)}dx-\frac12 \int_0^1 \frac{x\ln(1+x^2)-x\ln (x^2)}{1+x^2}dx$$ $$=-\frac12 \int_0^1 \frac{\ln(1+x^2)}{x}+\int_0^1 \frac{x\ln x}{1+x^2}dx=-\int_0^1 \frac{\ln(1+x^2)}{x}dx=-\frac{\pi^2}{24}$$ $$\Rightarrow I= 4J-\frac\pi4 \ln^22+\frac{\pi^3}{96}=-\frac{\pi}{4}\ln^2 2\Rightarrow P=\frac{\pi^3}{12}-\pi \ln^2 2$$