$\color{brown}{\textbf{Alternative expressions for the integral.}}$
Firstly,
$$I = \int\limits_0^{\Large^\pi\hspace{-1pt}/_2}\ln(1-x\cot x)\,\mathrm dx
= \int\limits_0^{\Large^\pi\hspace{-1pt}/_2}\ln(\sin x - x\cos x)\,\mathrm dx
- \int\limits_0^{\Large^\pi\hspace{-1pt}/_2}\ln(\sin x)\,\mathrm dx = \dfrac\pi2\ln2 +I_1,$$
wherein $I_1$ allows integration by parts:
$$I_1 = \int\limits_0^{\Large^\pi\hspace{-1pt}/_2}\ln(\sin x - x\cos x)\,\mathrm dx
= x\ln(\sin x-x\cos x)\bigg|_{\ 0}^{\Large^\pi\hspace{-1pt}/_2}
- \int\limits_0^{\Large^\pi\hspace{-1pt}/_2}\dfrac{x^2\sin x}{\sin x - x\cos x}\,\mathrm dx,$$
$$ I_1 =-\int\limits_0^{\Large^\pi\hspace{-1pt}/_2}\dfrac{x^2\sin x}{\sin x- x\cos x}\,\mathrm dx
= - \int\limits_0^{\Large^\pi\hspace{-1pt}/_2}\dfrac{x^2}{1-x\cot x}\,\mathrm dx = - J_{21},\tag1$$
where
$$J_{mn} = \int\limits_0^{\Large^\pi\hspace{-1pt}/_2}\dfrac{x^m}{(1-x\cot x)^n}\,\mathrm dx.\tag2$$
On the other hand,
$$J_{21} = \int\limits_0^{\Large^\pi\hspace{-1pt}/_2}\dfrac{x^2(1-x\cot x + x\cot x)}{1-x\cot x}\,\mathrm dx = \dfrac{\pi^3}{24} + I_2,$$
where
$$I_2 = \int\limits_0^{\Large^\pi\hspace{-1pt}/_2}\dfrac{x^3\cot x}{1 - x\cot x}\,\mathrm dx
= \int\limits_0^{\Large^\pi\hspace{-1pt}/_2}\dfrac{x^3}{\tan x - x}\,\mathrm dx.\tag3$$
Formulas $(3)$ are not suitable for numeric calculations.
But integration by parts is possible,
$$I_2 = \dfrac14\int\limits_0^{\Large^\pi\hspace{-1pt}/_2}\dfrac{1}{\tan x - x}\,\mathrm dx^4
= \dfrac14\dfrac{x^4}{\tan x-x}\bigg|_{\,0}^{\Large^\pi\hspace{-1pt}/_2}
+ \dfrac14\int\limits_0^{\Large^\pi\hspace{-1pt}/_2}\dfrac{x^4(1+\tan^2x -1)}{(\tan x - x)^2}\,\mathrm dx,$$
$$ I_2 = \dfrac14\int\limits_0^{\Large^\pi\hspace{-1pt}/_2}\dfrac{x^4}{(1 - x\cot x)^2}\,\mathrm dx = \dfrac14 J_{42},$$
$$I = \dfrac\pi2\ln2 - \dfrac{\pi^3}{24} - \dfrac14 J_{42}.\tag4$$
Formula $(4)$ provides both suitable numeric calculations via Wolfram Alpha by the expression
with the result
and the further building of the series in the elementary functions via the transformations in the form of
$$ J_{42} = \int\limits_0^{\Large^\pi\hspace{-1pt}/_2}\dfrac{x^4((1 - x\cot x)^2 + 2x\cot x(1 - x\cot x) + x^2\cot^2 x) }{(1 - x\cot x)^2}\,\mathrm dx\\
= \int\limits_0^{\Large^\pi\hspace{-1pt}/_2}\left(x^4 + 2\,\dfrac{x^5\cot x}{1-x\cot x} + \dfrac{x^6\cot^2x}{(1 - x\cot x)^2}\right)\,\mathrm dx\\
= \int\limits_0^{\Large^\pi\hspace{-1pt}/_2} x^4\,\mathrm dx
+ \int\limits_0^{\Large^\pi\hspace{-1pt}/_2}\left(\dfrac{2x^5\cot x}{1-x\cot x} + \dfrac{x^6\cot^2x}{(1 - x\cot x)^2}\right)\,\mathrm dx,$$
$$J_{42} = \dfrac{\pi^5}{160} + I_3 + I_4,\tag5$$
where
$$I_3 = 2\int\limits_0^{\Large^\pi\hspace{-1pt}/_2}\dfrac{x^5\cot x}{1-x\cot x} \,\mathrm dx
= 2\int\limits_0^{\Large^\pi\hspace{-1pt}/_2}\dfrac{x^5}{\tan x - x}\,\mathrm dx
= \dfrac13\int\limits_0^{\Large^\pi\hspace{-1pt}/_2}\dfrac{1}{\tan x - x}\,\mathrm dx^6\\
= \dfrac13\dfrac{x^6}{\tan x-x}\bigg|_{\,0}^{\Large^\pi\hspace{-1pt}/_2}
+ \dfrac13\int\limits_0^{\Large^\pi\hspace{-1pt}/_2}\dfrac{x^6(1+\tan^2x -1)}{(\tan x - x)^2}\,\mathrm dx = \dfrac13 J_{62},$$
$$I_4 = \int\limits_0^{\Large^\pi\hspace{-1pt}/_2}\dfrac{x^6\cot^2x}{(1 - x\cot x)^2}\,\mathrm dx
= \int\limits_0^{\Large^\pi\hspace{-1pt}/_2}\dfrac{x^6}{(\tan x - x)^2} \,\mathrm dx
= \dfrac17\int\limits_0^{\Large^\pi\hspace{-1pt}/_2}\dfrac{1}{(\tan x - x)^2} \,\mathrm dx^7\\
= \dfrac27\dfrac{x^7}{(\tan x-x)^3}\bigg|_{\,0}^{\Large^\pi\hspace{-1pt}/_2}
+ \dfrac27\int\limits_0^{\Large^\pi\hspace{-1pt}/_2}\dfrac{x^7(1+\tan^2x -1)}{(\tan x - x)^3}\,\mathrm dx
= \dfrac27\int\limits_0^{\Large^\pi\hspace{-1pt}/_2}\dfrac{x^7\cot x}{(1 - x\cot x)^3}\,\mathrm dx\\
= \dfrac27\int\limits_0^{\Large^\pi\hspace{-1pt}/_2}\dfrac{x^6(1 - (1 - x\cot x))}{(1 - x\cot x)^3}\,\mathrm dx =\dfrac27(J_{63}-J_{62}),$$
Therefore,
$$I = \dfrac\pi2\ln2 - \dfrac{\pi^3}{24} - \dfrac{\pi^5}{640} - \dfrac1{84}J_{62} - \dfrac1{14}J_{63}.\tag6$$
Numeric calculations via Mathcad Alpha by the formula $(6)$
leads to the same result, and this confirms the correctness of the approach.
$\color{brown}{\textbf{Recurrence relations.}}$
For the arbitrary $m,n$
$$ J_{mn} = \int\limits_0^{\Large^\pi\hspace{-1pt}/_2}\,(x\cot x + (1-x\cot x))^n \dfrac{x^m}{(1 - x\cot x)^n}\,\mathrm dx
= \int\limits_0^{\Large^\pi\hspace{-1pt}/_2}\sum\limits_{k=0}^n\binom nk\dfrac{x^{m+k}\cot^k x}{(1 - x\cot x)^k}\,\mathrm dx
= \dfrac{\pi^{m+1}}{(m+1)2^{m+1}} + \sum\limits_{k=1}^n\dfrac{\dbinom nk}{m+k+1} \int\limits_0^{\Large^\pi\hspace{-1pt}/_2}\dfrac{\mathrm dx^{m+k+1}}{(\tan x - x)^k}
= \dfrac{\pi^{m+1}}{(m+1)2^{m+1}}\\
+ \sum\limits_{k=1}^n\dfrac{\dbinom nk}{m+k+1} \left(\dfrac{x^{m+k+1}}{(\tan x - x)^{k}}\bigg|_{\,0}^{\Large^\pi\hspace{-1pt}/_2}
+ k\int\limits_0^{\Large^\pi\hspace{-1pt}/_2}\dfrac{x^{m+k+1}(1+\tan^2x-1)}{(\tan x-x)^{k+1}}\,\mathrm dx\right)\\
= \dfrac{\pi^{m+1}}{(m+1)2^{m+1}}
+ \sum\limits_{k=1}^n\dfrac{k}{m+k+1} \dbinom nk
\int\limits_0^{\Large^\pi\hspace{-1pt}/_2}\dfrac{x^{m+2}(x\cot x)^{k-1}}{(1 -x\cot x)^{k+1}}\,\mathrm dx\\
= \dfrac{\pi^{m+1}}{(m+1)2^{m+1}}
+ \sum\limits_{k=1}^n\dfrac{k}{m+k+1} \dbinom nk
\int\limits_0^{\Large^\pi\hspace{-1pt}/_2}\dfrac{x^{m+2}(1-(1-x\cot x))^{k-1}}{(1 -x\cot x)^{k+1}}\,\mathrm dx\\
= \dfrac{\pi^{m+1}}{(m+1)2^{m+1}}
+ \sum\limits_{k=1}^n\dfrac{k}{m+k+1} \dbinom nk \sum\limits_{j=0}^{k-1}(-1)^{k-1-j}\dbinom{k-1}j J_{m+2,\,j+2},$$
$$J_{mn} = \dfrac{\pi^{m+1}}{(m+1)2^{m+1}}
+ \sum\limits_{j=0}^{n-1} F_{j} J_{m+2,\,j+2},\tag7$$
where
$$F_{j} = \sum\limits_{k=j+1}^n (-1)^{k-1-j} \dfrac{k}{m+k+1}\dbinom nk \dbinom{k-1}j.\tag8 $$
If $(m,n)=(2,1),\ $ then
$$F_{0} = \sum\limits_{k=1}^1 (-1)^{k-1} \dfrac{k}{2+k+1}\dbinom1k \dbinom{k-1}0
=\dfrac14,$$
$$J_{21} = \dfrac{\pi^{3}}{3\cdot2^3} + \sum\limits_{j=0}^0 F_{j} J_{4,\,j+2}
= \dfrac{\pi^{3}}{24} + J_{42}.$$
If $(m,n)=(4,2),\ $ then
$$F_{0} = \sum\limits_{k=1}^2 (-1)^{k-1} \dfrac{k}{4+k+1}\dbinom2k \dbinom{k-1}0
=\dfrac13 - \dfrac27 = \dfrac{1}{21},$$
$$F_{1} = \sum\limits_{k=2}^2 (-1)^{k} \dfrac{k}{4+k+1}\dbinom2k \dbinom{k-1}1
=\dfrac27,$$
$$J_{42} = \dfrac{\pi^{5}}{5\cdot2^5} + \sum\limits_{j=0}^1 F_{j} J_{2,\,j+2}
= \dfrac{\pi^5}{160} + \dfrac1{21}J_{62} + \dfrac27J_{63}.$$
Similarly,
$$J_{62} = \dfrac{\pi^7}{896}+\dfrac1{36}J_{82}+\dfrac29J_{83}\tag9$$
(see also Wolfram Alpha test).
Besides,
$$J_{63} = \dfrac{\pi^7}{896}+\dfrac1{120}J_{82} + \dfrac1{15}J_{83} + \dfrac3{20}J_{84}.\tag{10}$$
$\color{brown}{\textbf{Simple series.}}$
Obtained results are not the best way to get required series of the arbitrary length.
$$\boxed{
\begin{matrix}
I & = & -3.35333726288947201778500718670823032009876022464933939598 \\
\frac\pi2\ln2 & = &1.088793045151801065250344449118806973669291850184643147162 \\
J_{21} & = & 4.442130308041273083035351635930890531086461245854584994170 \\
\frac{\pi^3}{24} & = & 1.291928195012492507311513127795891466759387023578546153922 \\
J_{42} & = & 12.60080845211512230289535403253999625730829688910415536099 \\
\frac{\pi^5}{160} & = & 1.912623029908009082892133187771472540501879416425468690959 \\
J_{62} & = & 9.357325953756236734147158157553707227832359838953032605558 \\
J_{63} & = & 35.84909465209885681432007993043088180418373451454989791084 \\
\frac{\pi^7}{896} & = & 3.370862977429455432493534032446475258836420173320761453966 \\
J_{82} & = & 13.21743446830609099759197972403428192140938899336281280188 \\
J_{83} & = & 25.28690408493225448274231109747825862030555487117486858192 \\
J_{84} & = & 102.2743092725712233044348622015074565154951081384648503713 \\
\end{matrix}}$$
On the other hand, using of the simple Laurent series for the function
$$g(y) = \dfrac{35}{1-y\sqrt{15}\cot y\sqrt{15}} = \dfrac7{y^2}-\sum\limits_{i=0}^\infty c_iy^{2i}\tag{11}$$
gives evidently convergent series
$$J_{21} = \dfrac1{35}\int\limits_0^{\Large^\pi\hspace{-1pt}/_2}
\left(7 - \sum\limits_{i=0}^\infty c_i\left(\dfrac{x^2}{15}\right)^{i+1}\right)\,\mathrm dx,$$
$$J_{21} = \dfrac32\pi - \dfrac3{14}\pi\sum\limits_{i=0}^\infty \dfrac{c_i}{2i+3}\left(\dfrac{\pi^2}{60}\right)^{i+1}\,\mathrm dx,\tag{12}$$
wherein the first $8$ terms provide the accuracy of $8$ decimal digits.
Define the function $\mathcal{I}:\mathbb{R}_{>0}^{2}\rightarrow\mathbb{R}$ via the definite integral
$$\mathcal{I}{\left(p,q\right)}:=\int_{0}^{1}\mathrm{d}t\,\frac{\arctan{\left(\sqrt{t^{2}+p^{2}}\right)}}{\left(t^{2}+q^{2}\right)\sqrt{t^{2}+p^{2}}}.\tag{1}$$
We seek a closed-form expression for $\mathcal{I}{\left(p,q\right)}$ for $\left(p,q\right)\in\mathbb{R}_{>0}^{2}$ such that $p>q$. It should be clear this problem is completely equivalent to the one posed by the OP.
We will make use of the following Euler substitution:
$$\sqrt{t^{2}+p^{2}}=t+x;~~~\small{p>0\land t\in\mathbb{R}\land x\in\mathbb{R}_{>0}}.\tag{2}$$
Solving for $t$, we have
$$t=\frac{p^{2}-x^{2}}{2x},$$
$$\implies\sqrt{t^{2}+p^{2}}=t+x=\frac{p^{2}+x^{2}}{2x},$$
$$\implies dt=dx\,\frac{(-1)\left(p^{2}+x^{2}\right)}{2x^{2}}.$$
Suppose $\left(p,q\right)\in\mathbb{R}_{>0}^{2}\land p>q$. We find
$$\begin{align}
\mathcal{I}{\left(p,q\right)}
&=\int_{0}^{1}\mathrm{d}t\,\frac{\arctan{\left(\sqrt{t^{2}+p^{2}}\right)}}{\left(t^{2}+q^{2}\right)\sqrt{t^{2}+p^{2}}}\\
&=\int_{0}^{1}\mathrm{d}t\,\frac{\left[\frac{\pi}{2}-\arctan{\left(\frac{1}{\sqrt{t^{2}+p^{2}}}\right)}\right]}{\left(t^{2}+q^{2}\right)\sqrt{t^{2}+p^{2}}}\\
&=\frac{\pi}{2}\int_{0}^{1}\mathrm{d}t\,\frac{1}{\left(t^{2}+q^{2}\right)\sqrt{t^{2}+p^{2}}}-\int_{0}^{1}\mathrm{d}t\,\frac{\arctan{\left(\frac{1}{\sqrt{t^{2}+p^{2}}}\right)}}{\left(t^{2}+q^{2}\right)\sqrt{t^{2}+p^{2}}}\\
&=\frac{\pi}{2}\int_{0}^{\frac{1}{\sqrt{1+p^{2}}}}\mathrm{d}u\,\frac{1}{\left(1-u^{2}\right)}\cdot\frac{1}{\left(\frac{pu}{\sqrt{1-u^{2}}}\right)^{2}+q^{2}};~~~\small{\left[t=\frac{pu}{\sqrt{1-u^{2}}}\right]}\\
&~~~~~-\int_{0}^{1}\mathrm{d}t\,\frac{\arctan{\left(\frac{1}{\sqrt{t^{2}+p^{2}}}\right)}}{\left(t^{2}+q^{2}\right)\sqrt{t^{2}+p^{2}}}\\
&=\frac{\pi}{2}\int_{0}^{\frac{1}{\sqrt{1+p^{2}}}}\mathrm{d}u\,\frac{1}{p^{2}u^{2}+q^{2}\left(1-u^{2}\right)}-\int_{0}^{1}\mathrm{d}t\,\frac{\arctan{\left(\frac{1}{\sqrt{t^{2}+p^{2}}}\right)}}{\left(t^{2}+q^{2}\right)\sqrt{t^{2}+p^{2}}}\\
&=\frac{\pi}{2}\int_{0}^{\frac{1}{\sqrt{1+p^{2}}}}\mathrm{d}u\,\frac{1}{q^{2}+\left(p^{2}-q^{2}\right)u^{2}}-\int_{0}^{1}\mathrm{d}t\,\frac{\arctan{\left(\frac{1}{\sqrt{t^{2}+p^{2}}}\right)}}{\left(t^{2}+q^{2}\right)\sqrt{t^{2}+p^{2}}}\\
&=\frac{\pi}{2}\int_{0}^{\frac{\sqrt{p^{2}-q^{2}}}{q\sqrt{1+p^{2}}}}\mathrm{d}v\,\frac{q}{\sqrt{p^{2}-q^{2}}}\cdot\frac{1}{q^{2}+q^{2}v^{2}};~~~\small{\left[u=\frac{qv}{\sqrt{p^{2}-q^{2}}}\right]}\\
&~~~~~-\int_{0}^{1}\mathrm{d}t\,\frac{\arctan{\left(\frac{1}{\sqrt{t^{2}+p^{2}}}\right)}}{\left(t^{2}+q^{2}\right)\sqrt{t^{2}+p^{2}}}\\
&=\frac{\pi}{2}\cdot\frac{1}{q\sqrt{p^{2}-q^{2}}}\int_{0}^{\frac{\sqrt{p^{2}-q^{2}}}{q\sqrt{1+p^{2}}}}\mathrm{d}v\,\frac{1}{1+v^{2}}-\int_{0}^{1}\mathrm{d}t\,\frac{\arctan{\left(\frac{1}{\sqrt{t^{2}+p^{2}}}\right)}}{\left(t^{2}+q^{2}\right)\sqrt{t^{2}+p^{2}}}\\
&=\frac{\pi\arctan{\left(\frac{\sqrt{p^{2}-q^{2}}}{q\sqrt{1+p^{2}}}\right)}}{2q\sqrt{p^{2}-q^{2}}}-\int_{0}^{1}\mathrm{d}t\,\frac{\arctan{\left(\frac{1}{\sqrt{t^{2}+p^{2}}}\right)}}{\left(t^{2}+q^{2}\right)\sqrt{t^{2}+p^{2}}}.\tag{3a}\\
\end{align}$$
Then $0<-1+\sqrt{1+p^{2}}<p$, and
$$\begin{align}
\mathcal{I}{\left(p,q\right)}
&=\frac{\pi\arctan{\left(\frac{\sqrt{p^{2}-q^{2}}}{q\sqrt{1+p^{2}}}\right)}}{2q\sqrt{p^{2}-q^{2}}}-\int_{0}^{1}\mathrm{d}t\,\frac{\arctan{\left(\frac{1}{\sqrt{t^{2}+p^{2}}}\right)}}{\left(t^{2}+q^{2}\right)\sqrt{t^{2}+p^{2}}}\\
&=\frac{\pi\arctan{\left(\frac{\sqrt{p^{2}-q^{2}}}{q\sqrt{1+p^{2}}}\right)}}{2q\sqrt{p^{2}-q^{2}}}\\
&~~~~~-\int_{p}^{-1+\sqrt{1+p^{2}}}\mathrm{d}x\,\frac{(-1)\left(p^{2}+x^{2}\right)}{2x^{2}}\cdot\frac{1}{\left(\frac{p^{2}-x^{2}}{2x}\right)^{2}+q^{2}}\cdot\frac{2x}{p^{2}+x^{2}}\\
&~~~~~\times\arctan{\left(\frac{2x}{p^{2}+x^{2}}\right)};~~~\small{\left[-t+\sqrt{t^{2}+p^{2}}=x\right]}\\
&=\frac{\pi\arctan{\left(\frac{\sqrt{p^{2}-q^{2}}}{q\sqrt{1+p^{2}}}\right)}}{2q\sqrt{p^{2}-q^{2}}}-\int_{-1+\sqrt{1+p^{2}}}^{p}\mathrm{d}x\,\frac{4x\arctan{\left(\frac{2x}{p^{2}+x^{2}}\right)}}{\left(p^{2}-x^{2}\right)^{2}+4q^{2}x^{2}}.\tag{3b}\\
\end{align}$$
It follows from the arctangent addition formula that
$$\arctan{\left(\frac{x}{-1+\sqrt{1+p^{2}}}\right)}-\arctan{\left(\frac{x}{1+\sqrt{1+p^{2}}}\right)}=\arctan{\left(\frac{2x}{p^{2}+x^{2}}\right)};~~~\small{p\in\mathbb{R}_{>0}\land x\in\mathbb{R}}.\tag{4}$$
Set $r:=\frac{p}{1+\sqrt{1+p^{2}}}\land\sigma:=\arcsin{\left(\frac{q}{p}\right)}$. Then, $0<r<1\land0<\sigma<\frac{\pi}{2}\land r^{-1}=\frac{1+\sqrt{1+p^{2}}}{p}=\frac{p}{-1+\sqrt{1+p^{2}}}$, and
$$\begin{align}
\mathcal{I}{\left(p,q\right)}
&=\frac{\pi\arctan{\left(\frac{\sqrt{p^{2}-q^{2}}}{q\sqrt{1+p^{2}}}\right)}}{2q\sqrt{p^{2}-q^{2}}}-\int_{-1+\sqrt{1+p^{2}}}^{p}\mathrm{d}x\,\frac{4x\arctan{\left(\frac{2x}{p^{2}+x^{2}}\right)}}{\left(p^{2}-x^{2}\right)^{2}+4q^{2}x^{2}}\\
&=\frac{\pi\arctan{\left(\frac{\sqrt{p^{2}-q^{2}}}{q\sqrt{1+p^{2}}}\right)}}{2q\sqrt{p^{2}-q^{2}}}-\int_{-1+\sqrt{1+p^{2}}}^{p}\mathrm{d}x\,\frac{4x}{\left(p^{2}-x^{2}\right)^{2}+4q^{2}x^{2}}\\
&~~~~~\times\left[\arctan{\left(\frac{x}{-1+\sqrt{1+p^{2}}}\right)}-\arctan{\left(\frac{x}{1+\sqrt{1+p^{2}}}\right)}\right]\\
&=\frac{\pi\arctan{\left(\frac{\sqrt{p^{2}-q^{2}}}{q\sqrt{1+p^{2}}}\right)}}{2q\sqrt{p^{2}-q^{2}}}-\int_{\frac{-1+\sqrt{1+p^{2}}}{p}}^{1}\mathrm{d}y\,\frac{4p^{2}y}{\left(p^{2}-p^{2}y^{2}\right)^{2}+4p^{2}q^{2}y^{2}}\\
&~~~~~\times\left[\arctan{\left(\frac{py}{-1+\sqrt{1+p^{2}}}\right)}-\arctan{\left(\frac{py}{1+\sqrt{1+p^{2}}}\right)}\right];~~~\small{\left[x=py\right]}\\
&=\frac{\pi\arctan{\left(\frac{\sqrt{p^{2}-q^{2}}}{q\sqrt{1+p^{2}}}\right)}}{2q\sqrt{p^{2}-q^{2}}}-\frac{1}{p^{2}}\int_{\frac{-1+\sqrt{1+p^{2}}}{p}}^{1}\mathrm{d}y\,\frac{4y}{\left(1-y^{2}\right)^{2}+4\left(\frac{q}{p}\right)^{2}y^{2}}\\
&~~~~~\times\left[\arctan{\left(\frac{py}{-1+\sqrt{1+p^{2}}}\right)}-\arctan{\left(\frac{py}{1+\sqrt{1+p^{2}}}\right)}\right]\\
&=\frac{\pi\arctan{\left(\frac{\sqrt{p^{2}-q^{2}}}{q\sqrt{1+p^{2}}}\right)}}{2q\sqrt{p^{2}-q^{2}}}-\frac{1}{p^{2}}\int_{r}^{1}\mathrm{d}y\,\frac{4y}{\left(1-y^{2}\right)^{2}+4y^{2}\sin^{2}{\left(\sigma\right)}}\\
&~~~~~\times\left[\arctan{\left(\frac{y}{r}\right)}-\arctan{\left(ry\right)}\right]\\
&=\frac{\pi\arctan{\left(\frac{\sqrt{p^{2}-q^{2}}}{q\sqrt{1+p^{2}}}\right)}}{2q\sqrt{p^{2}-q^{2}}}-\frac{1}{p^{2}}\int_{r}^{1}\mathrm{d}y\,\frac{4y\left[\arctan{\left(\frac{y}{r}\right)}-\arctan{\left(ry\right)}\right]}{1-2y^{2}\left[1-2\sin^{2}{\left(\sigma\right)}\right]+y^{4}}\\
&=\frac{\pi\arctan{\left(\frac{\sqrt{p^{2}-q^{2}}}{q\sqrt{1+p^{2}}}\right)}}{2q\sqrt{p^{2}-q^{2}}}-\frac{1}{p^{2}}\int_{r}^{1}\mathrm{d}y\,\frac{4y\left[\arctan{\left(\frac{y}{r}\right)}-\arctan{\left(ry\right)}\right]}{1-2y^{2}\cos{\left(2\sigma\right)}+y^{4}}\\
&=\frac{\pi\arctan{\left(\frac{\sqrt{p^{2}-q^{2}}}{q\sqrt{1+p^{2}}}\right)}}{2q\sqrt{p^{2}-q^{2}}}\\
&~~~~~-\frac{1}{p^{2}}\int_{r}^{1}\mathrm{d}y\,\frac{4y\left[\arctan{\left(\frac{y}{r}\right)}-\arctan{\left(ry\right)}\right]}{\left[1-2y\cos{\left(\sigma\right)}+y^{2}\right]\left[1+2y\cos{\left(\sigma\right)}+y^{2}\right]}\\
&=\frac{\pi\arctan{\left(\frac{\sqrt{p^{2}-q^{2}}}{q\sqrt{1+p^{2}}}\right)}}{2q\sqrt{p^{2}-q^{2}}}\\
&~~~~~-\frac{1}{p^{2}\cos{\left(\sigma\right)}}\int_{r}^{1}\mathrm{d}y\,\frac{4y\cos{\left(\sigma\right)}}{\left[1-2y\cos{\left(\sigma\right)}+y^{2}\right]\left[1+2y\cos{\left(\sigma\right)}+y^{2}\right]}\\
&~~~~~\times\left[\arctan{\left(\frac{y}{r}\right)}-\arctan{\left(ry\right)}\right]\\
&=\frac{\pi\arctan{\left(\frac{\sqrt{p^{2}-q^{2}}}{q\sqrt{1+p^{2}}}\right)}}{2q\sqrt{p^{2}-q^{2}}}\\
&~~~~~-\frac{1}{p^{2}\cos{\left(\sigma\right)}}\int_{r}^{1}\mathrm{d}y\,\left[\frac{1}{1-2y\cos{\left(\sigma\right)}+y^{2}}-\frac{1}{1+2y\cos{\left(\sigma\right)}+y^{2}}\right]\\
&~~~~~\times\left[\arctan{\left(\frac{y}{r}\right)}-\arctan{\left(ry\right)}\right]\\
&=\frac{\pi\arctan{\left(\frac{\sqrt{p^{2}-q^{2}}}{q\sqrt{1+p^{2}}}\right)}}{2q\sqrt{p^{2}-q^{2}}}\\
&~~~~~-\frac{1}{p^{2}\sin{\left(\sigma\right)}\cos{\left(\sigma\right)}}\int_{r}^{1}\mathrm{d}y\,\left[\frac{\sin{\left(\sigma\right)}}{1-2y\cos{\left(\sigma\right)}+y^{2}}-\frac{\sin{\left(\sigma\right)}}{1+2y\cos{\left(\sigma\right)}+y^{2}}\right]\\
&~~~~~\times\left[\arctan{\left(\frac{y}{r}\right)}-\arctan{\left(ry\right)}\right]\\
&=\frac{\pi\arctan{\left(\frac{\sqrt{p^{2}-q^{2}}}{q\sqrt{1+p^{2}}}\right)}}{2q\sqrt{p^{2}-q^{2}}}\\
&~~~~~-\frac{1}{q\sqrt{p^{2}-q^{2}}}\int_{r}^{1}\mathrm{d}y\,\left[\frac{\sin{\left(\sigma\right)}}{1-2y\cos{\left(\sigma\right)}+y^{2}}-\frac{\sin{\left(\sigma\right)}}{1+2y\cos{\left(\sigma\right)}+y^{2}}\right]\\
&~~~~~\times\left[\arctan{\left(\frac{y}{r}\right)}-\arctan{\left(ry\right)}\right].\tag{5}\\
\end{align}$$
As such, let's introduce another auxiliary function $\mathcal{J}:\left(0,1\right)\times\left(0,\frac{\pi}{2}\right)\rightarrow\mathbb{R}$ defined via the last definite integral above:
$$\mathcal{J}{\left(r,\sigma\right)}:=\int_{r}^{1}\mathrm{d}y\,\left[\frac{\sin{\left(\sigma\right)}}{1-2y\cos{\left(\sigma\right)}+y^{2}}-\frac{\sin{\left(\sigma\right)}}{1+2y\cos{\left(\sigma\right)}+y^{2}}\right]\left[\arctan{\left(\frac{y}{r}\right)}-\arctan{\left(ry\right)}\right].\tag{6}$$
Given $\left(r,\sigma\right)\in\left(0,1\right)\times\left(0,\frac{\pi}{2}\right)$, we find
$$\begin{align}
\mathcal{J}{\left(r,\sigma\right)}
&=\int_{r}^{1}\mathrm{d}y\,\left[\frac{\sin{\left(\sigma\right)}}{1-2y\cos{\left(\sigma\right)}+y^{2}}-\frac{\sin{\left(\sigma\right)}}{1+2y\cos{\left(\sigma\right)}+y^{2}}\right]\\
&~~~~~\times\left[\arctan{\left(\frac{y}{r}\right)}-\arctan{\left(ry\right)}\right]\\
&=\int_{r}^{1}\mathrm{d}y\,\left[\frac{\sin{\left(\sigma\right)}}{1-2y\cos{\left(\sigma\right)}+y^{2}}-\frac{\sin{\left(\sigma\right)}}{1+2y\cos{\left(\sigma\right)}+y^{2}}\right]\arctan{\left(\frac{y}{r}\right)}\\
&~~~~~-\int_{r}^{1}\mathrm{d}y\,\left[\frac{\sin{\left(\sigma\right)}}{1-2y\cos{\left(\sigma\right)}+y^{2}}-\frac{\sin{\left(\sigma\right)}}{1+2y\cos{\left(\sigma\right)}+y^{2}}\right]\arctan{\left(ry\right)}\\
&=\int_{1}^{\frac{1}{r}}\mathrm{d}t\,\left[\frac{r\sin{\left(\sigma\right)}}{1-2rt\cos{\left(\sigma\right)}+r^{2}t^{2}}-\frac{r\sin{\left(\sigma\right)}}{1+2rt\cos{\left(\sigma\right)}+r^{2}t^{2}}\right]\arctan{\left(t\right)};~~~\small{\left[y=rt\right]}\\
&~~~~~-\int_{r^{2}}^{r}\mathrm{d}u\,\left[\frac{r\sin{\left(\sigma\right)}}{r^{2}-2ru\cos{\left(\sigma\right)}+u^{2}}-\frac{r\sin{\left(\sigma\right)}}{r^{2}+2ru\cos{\left(\sigma\right)}+u^{2}}\right]\arctan{\left(u\right)};~~~\small{\left[y=\frac{u}{r}\right]}\\
&=\int_{r}^{1}\mathrm{d}u\,\left[\frac{r\sin{\left(\sigma\right)}}{r^{2}-2ru\cos{\left(\sigma\right)}+u^{2}}-\frac{r\sin{\left(\sigma\right)}}{r^{2}+2ru\cos{\left(\sigma\right)}+u^{2}}\right]\\
&~~~~~\times\arctan{\left(\frac{1}{u}\right)};~~~\small{\left[t=u^{-1}\right]}\\
&~~~~~-\int_{r^{2}}^{r}\mathrm{d}u\,\left[\frac{r\sin{\left(\sigma\right)}}{r^{2}-2ru\cos{\left(\sigma\right)}+u^{2}}-\frac{r\sin{\left(\sigma\right)}}{r^{2}+2ru\cos{\left(\sigma\right)}+u^{2}}\right]\arctan{\left(u\right)}\\
&=\int_{r}^{1}\mathrm{d}u\,\left[\frac{r\sin{\left(\sigma\right)}}{r^{2}-2ru\cos{\left(\sigma\right)}+u^{2}}-\frac{r\sin{\left(\sigma\right)}}{r^{2}+2ru\cos{\left(\sigma\right)}+u^{2}}\right]\left[\frac{\pi}{2}-\arctan{\left(u\right)}\right]\\
&~~~~~-\int_{r^{2}}^{r}\mathrm{d}u\,\left[\frac{r\sin{\left(\sigma\right)}}{r^{2}-2ru\cos{\left(\sigma\right)}+u^{2}}-\frac{r\sin{\left(\sigma\right)}}{r^{2}+2ru\cos{\left(\sigma\right)}+u^{2}}\right]\arctan{\left(u\right)}\\
&=\frac{\pi}{2}\int_{r}^{1}\mathrm{d}u\,\left[\frac{r\sin{\left(\sigma\right)}}{r^{2}-2ru\cos{\left(\sigma\right)}+u^{2}}-\frac{r\sin{\left(\sigma\right)}}{r^{2}+2ru\cos{\left(\sigma\right)}+u^{2}}\right]\\
&~~~~~-\int_{r}^{1}\mathrm{d}u\,\left[\frac{r\sin{\left(\sigma\right)}}{r^{2}-2ru\cos{\left(\sigma\right)}+u^{2}}-\frac{r\sin{\left(\sigma\right)}}{r^{2}+2ru\cos{\left(\sigma\right)}+u^{2}}\right]\arctan{\left(u\right)}\\
&~~~~~-\int_{r^{2}}^{r}\mathrm{d}u\,\left[\frac{r\sin{\left(\sigma\right)}}{r^{2}-2ru\cos{\left(\sigma\right)}+u^{2}}-\frac{r\sin{\left(\sigma\right)}}{r^{2}+2ru\cos{\left(\sigma\right)}+u^{2}}\right]\arctan{\left(u\right)}\\
&=\frac{\pi}{2}\int_{r}^{1}\mathrm{d}u\,\frac{d}{du}\left[\arctan{\left(\frac{u-r\cos{\left(\sigma\right)}}{r\sin{\left(\sigma\right)}}\right)}-\arctan{\left(\frac{u+r\cos{\left(\sigma\right)}}{r\sin{\left(\sigma\right)}}\right)}\right]\\
&~~~~~-\int_{r^{2}}^{1}\mathrm{d}u\,\left[\frac{r\sin{\left(\sigma\right)}}{r^{2}-2ru\cos{\left(\sigma\right)}+u^{2}}-\frac{r\sin{\left(\sigma\right)}}{r^{2}+2ru\cos{\left(\sigma\right)}+u^{2}}\right]\arctan{\left(u\right)}\\
&=\frac{\pi}{2}\bigg{[}\arctan{\left(\frac{1-r\cos{\left(\sigma\right)}}{r\sin{\left(\sigma\right)}}\right)}-\arctan{\left(\frac{1+r\cos{\left(\sigma\right)}}{r\sin{\left(\sigma\right)}}\right)}\\
&~~~~~-\arctan{\left(\frac{1-\cos{\left(\sigma\right)}}{\sin{\left(\sigma\right)}}\right)}+\arctan{\left(\frac{1+\cos{\left(\sigma\right)}}{\sin{\left(\sigma\right)}}\right)}\bigg{]}\\
&~~~~~+\int_{r^{2}}^{1}\mathrm{d}u\,\frac{r\sin{\left(\sigma\right)}}{r^{2}+2ru\cos{\left(\sigma\right)}+u^{2}}\arctan{\left(u\right)}\\
&~~~~~-\int_{r^{2}}^{1}\mathrm{d}u\,\frac{r\sin{\left(\sigma\right)}}{r^{2}-2ru\cos{\left(\sigma\right)}+u^{2}}\arctan{\left(u\right)}\\
&=\frac{\pi}{2}\left[\frac{\pi}{2}-\sigma-\arctan{\left(\frac{r^{2}\sin{\left(2\sigma\right)}}{1-r^{2}\cos{\left(2\sigma\right)}}\right)}\right]\\
&~~~~~+\int_{r^{2}}^{1}\mathrm{d}u\,\frac{r\sin{\left(\sigma\right)}}{r^{2}+2ru\cos{\left(\sigma\right)}+u^{2}}\arctan{\left(u\right)}\\
&~~~~~+\int_{-1}^{-r^{2}}\mathrm{d}u\,\frac{r\sin{\left(\sigma\right)}}{r^{2}+2ru\cos{\left(\sigma\right)}+u^{2}}\arctan{\left(u\right)};~~~\small{\left[u\mapsto-u\right]}\\
&=\frac{\pi}{2}\left[\frac{\pi}{2}-\sigma-\arctan{\left(\frac{r^{2}\sin{\left(2\sigma\right)}}{1-r^{2}\cos{\left(2\sigma\right)}}\right)}\right]\\
&~~~~~+\int_{-1}^{1}\mathrm{d}u\,\frac{r\sin{\left(\sigma\right)}}{r^{2}+2ru\cos{\left(\sigma\right)}+u^{2}}\arctan{\left(u\right)}\\
&~~~~~-\int_{-r^{2}}^{r^{2}}\mathrm{d}u\,\frac{r\sin{\left(\sigma\right)}}{r^{2}+2ru\cos{\left(\sigma\right)}+u^{2}}\arctan{\left(u\right)}.\tag{7}\\
\end{align}$$
To facilitate the evaluation of the remaining two integrals in the last line above, we introduce yet another auxiliary function $\mathcal{K}:\left(0,1\right)\times\left(0,\frac{\pi}{2}\right)\times\left(0,1\right]\rightarrow\mathbb{R}$ defined via the definite integral
$$\mathcal{K}{\left(r,\sigma,z\right)}:=\int_{-z}^{z}\mathrm{d}u\,\frac{r\sin{\left(\sigma\right)}}{r^{2}+2ru\cos{\left(\sigma\right)}+u^{2}}\arctan{\left(u\right)}.\tag{8}$$
Suppose $\left(r,\sigma,z\right)\in\left(0,1\right)\times\left(0,\frac{\pi}{2}\right)\times\left(0,1\right]$. We then have
$$\begin{align}
\mathcal{K}{\left(r,\sigma,z\right)}
&=\int_{-z}^{z}\mathrm{d}u\,\frac{r\sin{\left(\sigma\right)}}{r^{2}+2ru\cos{\left(\sigma\right)}+u^{2}}\arctan{\left(u\right)}\\
&=\int_{-z}^{z}\mathrm{d}u\,\arctan{\left(u\right)}\frac{d}{du}\arctan{\left(\frac{u+r\cos{\left(\sigma\right)}}{r\sin{\left(\sigma\right)}}\right)}\\
&=\arctan{\left(z\right)}\arctan{\left(\frac{z+r\cos{\left(\sigma\right)}}{r\sin{\left(\sigma\right)}}\right)}-\arctan{\left(-z\right)}\arctan{\left(\frac{-z+r\cos{\left(\sigma\right)}}{r\sin{\left(\sigma\right)}}\right)}\\
&~~~~~-\int_{-z}^{z}\mathrm{d}u\,\arctan{\left(\frac{u+r\cos{\left(\sigma\right)}}{r\sin{\left(\sigma\right)}}\right)}\frac{d}{du}\arctan{\left(u\right)};~~~\small{I.B.P.}\\
&=\arctan{\left(z\right)}\arctan{\left(\frac{z+r\cos{\left(\sigma\right)}}{r\sin{\left(\sigma\right)}}\right)}-\arctan{\left(z\right)}\arctan{\left(\frac{z-r\cos{\left(\sigma\right)}}{r\sin{\left(\sigma\right)}}\right)}\\
&~~~~~-\int_{-z}^{z}\mathrm{d}u\,\frac{1}{1+u^{2}}\arctan{\left(\frac{r\cos{\left(\sigma\right)}+u}{r\sin{\left(\sigma\right)}}\right)}\\
&=\arctan{\left(z\right)}\left[\arctan{\left(\frac{z+r\cos{\left(\sigma\right)}}{r\sin{\left(\sigma\right)}}\right)}-\arctan{\left(\frac{z-r\cos{\left(\sigma\right)}}{r\sin{\left(\sigma\right)}}\right)}\right]\\
&~~~~~-\int_{-z}^{z}\mathrm{d}u\,\frac{1}{1+u^{2}}\arctan{\left(\cot{\left(\sigma\right)}+r^{-1}u\csc{\left(\sigma\right)}\right)}\\
&=\arctan{\left(z\right)}\left[\arctan{\left(\frac{z+r\cos{\left(\sigma\right)}}{r\sin{\left(\sigma\right)}}\right)}-\arctan{\left(\frac{z-r\cos{\left(\sigma\right)}}{r\sin{\left(\sigma\right)}}\right)}\right]\\
&~~~~~-\int_{-\arctan{\left(z\right)}}^{\arctan{\left(z\right)}}\mathrm{d}\varphi\,\arctan{\left(\cot{\left(\sigma\right)}+r^{-1}\csc{\left(\sigma\right)}\tan{\left(\varphi\right)}\right)};~~~\small{\left[\arctan{\left(u\right)}=\varphi\right]}.\tag{9}\\
\end{align}$$
Suppose $\left(a,b,\psi,\omega\right)\in\mathbb{R}\times\mathbb{R}\times\left(-\frac{\pi}{2},\frac{\pi}{2}\right)\times\left(-\frac{\pi}{2},\frac{\pi}{2}\right)$, and assume
$$0<a\land\sqrt{1+a^{2}}<b\land-\frac{\pi}{4}\le\psi<\omega\le\frac{\pi}{4}.$$
Set
$$B:=\frac{\left[\sqrt{a^{2}+\left(b+1\right)^{2}}+\sqrt{a^{2}+\left(b-1\right)^{2}}\right]^{2}}{4b}\land\phi:=\frac12\arctan{\left(\frac{2ab}{b^{2}-a^{2}-1}\right)},$$
and note that $1<B\land0<\phi<\frac{\pi}{4}$. Then, it can be shown that
$$\begin{align}
\int_{\psi}^{\omega}\mathrm{d}\varphi\,\arctan{\left(a+b\tan{\left(\varphi\right)}\right)}
&=\frac12\left(\phi+\omega\right)^{2}-\frac12\left(\phi+\psi\right)^{2}\\
&~~~~~-\left(\omega-\psi\right)\arctan{\left(\frac{\tan{\left(\phi\right)}}{B}\right)}\\
&~~~~~-\frac12\operatorname{Li}_{2}{\left(\frac{1-B}{1+B},\pi-2\phi-2\omega\right)}\\
&~~~~~+\frac12\operatorname{Li}_{2}{\left(\frac{1-B}{1+B},\pi-2\phi-2\psi\right)}.\tag{10}\\
\end{align}$$
The integration formula above is derived in this question.
Integration formula $(10)$ is sufficient to provide us with a closed-form expression for $\mathcal{K}$, and in turn $\mathcal{J}$ and $\mathcal{I}$ as well. So our work is complete in principle, and all that remains is to substitute back the chain of results to obtain a final expression in terms of the original variables. (Forgive me if I don't bother to do that here since the end result is such a cumbersome expression.)
Best Answer
Denoting $\displaystyle\alpha=\frac1{2r}$, the answer is the following: $$\Phi_0(r) = 4\pi r \int_0^{\infty} \frac{\sin(t)}{t} \left( \frac{\sin(r t)}{t(\pi^2-r^2 t^2)}\right)^2 \mathrm{d} t$$ $$=\frac{r^3}{\pi^2}\left(2\alpha(2-\alpha)+\frac{1-\alpha}\pi\sin(2\pi\alpha)+\frac2{\pi^2}\big(1-\cos(2\pi\alpha)\big)\right), \alpha\in[0;1]\,\,\big(r\geqslant 1/2\big)\tag{a}$$ $$=\frac{2r^3}{\pi^2};\,\,\alpha\geqslant1\,\big(r\in[0;1/2]\big)\tag{b}$$ Making the substitution $t=\frac xr$ $$\Phi_0(r)=4\pi r^3\int_0^\infty\frac{\sin\frac xr}x\left(\frac{\sin x}{x(\pi^2-x^2)}\right)^2dx\tag{1}$$ $$\Phi_1(r)=4\pi r^2\int_0^\infty\cos\frac xr\left(\frac{\sin x}{x(\pi^2-x^2)}\right)^2dx=r^2\frac d{d\big(\frac1r\big)}\frac{\Phi_0(r)}{r^3}$$ $$\Phi_2(r)=4\pi r\int_0^\infty x\sin\frac xr\left(\frac{\sin x}{x(\pi^2-x^2)}\right)^2dx=-r\,\frac d{d\big(\frac1r\big)}\frac{\Phi_1(r)}{r^2}$$ The evaluation is straightforward. Using (1) and $\sin^2x=\frac{1-\cos(2x)}2$, we can present $\Phi_0(r)$ in the form $$\Phi_0(r)=128\pi r^3\int_0^\infty\frac{\sin\frac x{2r}(1-\cos x)}{x^3\big((2\pi)^2-x^2)^2}dx$$ Using for convenience the notation $\displaystyle \alpha=\frac1{2r}$ $$\Phi_0(r)=64\pi r^3\int_{-\infty}^\infty\frac{\sin\alpha x-\frac12\sin x(1+\alpha)-\frac12\sin x(\alpha-1)}{x^3\big((2\pi)^2-x^2\big)^2}dx$$ $$=64\pi r^3\Im\int_{-\infty}^\infty\frac{e^{i\alpha x}-\frac12e^{ix(1+\alpha)}-\frac12e^{ix(\alpha-1)}}{x^3\big((2\pi)^2-x^2\big)^2}dx$$ Now we go in the complex plane and create the closed contour $C$, adding a big half-circle (counter-clockwise) in the upper half-plane, and three small half-circles (around the points $-2\pi; 0; 2\pi$, clockwise). There are no poles inside the contour; the integral along a big arch $\to0$. Therefore, at $\alpha>1$ (or at $r\in[0;1/2]$) $$J(\alpha)=\oint_C\frac{e^{i\alpha z}-\frac12e^{iz(1+\alpha)}-\frac12e^{iz(\alpha-1)}}{z^3\big((2\pi)^2-z^2\big)^2}dz=\pi\,i \underset{\binom{z=\pm2\pi}{z=0}}{\operatorname{Res}}\frac{e^{i\alpha z}-\frac12e^{iz(1+\alpha)}-\frac12e^{iz(\alpha-1)}}{z^3\big((2\pi)^2-z^2\big)^2}$$ $$\Phi_0(r)=64\pi r^3\Im J(\alpha)=64\pi^2r^3\Re\underset{\binom{z=\pm2\pi}{z=0}}{\operatorname{Res}}\frac{e^{i\alpha z}-\frac12e^{iz(1+\alpha)}-\frac12e^{iz(\alpha-1)}}{z^3\big((2\pi)^2-z^2\big)^2}$$ At $\alpha<1$ we present $\Phi_0(r)$ in the form $$\Phi_0(r)=64\pi r^3\int_{-\infty}^\infty\frac{\sin\alpha x-\frac12\sin x(1+\alpha)+\frac12\sin x(1-\alpha)}{x^3\big((2\pi)^2-x^2\big)^2}dx$$ $$\Phi_0(r)=64\pi^2r^3\Re\underset{\binom{z=\pm2\pi}{z=0}}{\operatorname{Res}}\frac{e^{i\alpha z}-\frac12e^{iz(1+\alpha)}+\frac12e^{iz(1-\alpha)}}{z^3\big((2\pi)^2-z^2\big)^2}$$ The residues evaluation in both cases is straightforward and gives the answer presented above.
We expect that $\Phi_0(r)$ is a continuous function. Indeed, at $\alpha =1\,\,(a)=(b)$ (above), and at $\alpha\to 0 \,\,\Phi_0=0$.
Quick numeric check (I did it at $\alpha=\frac14\,\,(r=2)$ ) confirms the result. The answer also reproduces the asymptotics obtained by @Sal in his post.