Show that $\int\limits_0^{\frac{\pi}2}4\cos^2(x)\log^2(\cos x)~\mathrm dx=-\pi\log 2+\pi\log^2 2-\frac{\pi}2+\frac{\pi^3}{12}$ Within this collection of definite integrals, number $30.$ is given by

$$\int_0^{\frac{\pi}2}4\cos^2(x)\log^2(\cos x)\mathrm dx~=~-\pi\log 2+\pi\log^2 2-\frac{\pi}2+\frac{\pi^3}{12}\tag1$$

I figured out a way to evaluate the integral using a substitution and further the values of the derivatives of the Beta Function. My approach can be found on the bottom of the post. Another solution can be found here within the original post of the $(1)$.
I was just curious whether there is an easier attempt. First of all I thought about applying a Weierstraß Substitution to get rid of the trigonometry and instead solving an algebraic integral. Anyway it did not worked out quite well since the substitution made the new integrals even more complicated. Hence the logarithm is squared I guess series expansions$-$although the linked answer invoked the Fourier Series Expansion of $\log(\cos x)$$-$or IBP are not the right way to approach to this integral either.

I am asking for interesting or elegant ways to evaluate $(1)$. Please provide a different attempt than the two I suggested within this post if you are aware of one.

Thanks in advance!

Own evaluation

First of all rewrite the integral sligthly
$$\small\int_0^{\frac{\pi}2}4\cos^2(x)\log^2(\cos x)\mathrm dx=\int_0^{\frac{\pi}2}\cos^2(x)\log^2(\cos^2 x)\mathrm dx$$
To notice that the subsititution $\cos^2 x=y$ is suitable in this case. Computing $\mathrm dy$ and changing the borders of integration yields to
$$\begin{align}
\small\int_0^{\frac{\pi}2}\cos^2(x)\log^2(\cos^2 x)\mathrm dx&=\small\int_1^0 (y)\log^2(y)~\left(\frac12\frac{-\mathrm dy}{\sqrt{y(1-y)}}\right)\\
&=\small\int_0^1\frac{\sqrt{y}\log^2(y)}{\sqrt{1-y}}\mathrm dy
\end{align}$$
One may note the familiar structure similiar to one given by the second derivative of the Beta Function. Therefore this integral can be written as
$$\small\int_0^1\frac{\sqrt{y}\log^2(y)}{\sqrt{1-y}}\mathrm dy=\left.B_{xx}\left(x,\frac32\right)\right|_{x=\frac12}$$
The second derivative of the Beta Function is given in terms of the Beta Function connected with the Digamma and Trigamma Function as
$$\small B_{xx}(x,y)=B(x,y)[(\psi^{(0)}(x)+\psi^{(0)}(x+y))^2+\psi^{(1)}(x)-\psi^{(1)}(x+y)]$$
Plugging in the values $x=\frac12$ and $y=\frac32$ yields to
$$\begin{align}
\small\left.B_{xx}\left(x,\frac32\right)\right|_{x=\frac12}&=\small B\left(\frac12,\frac32\right)\left[\left(\psi^{(0)}\left(\frac12\right)+\psi^{(0)}\left(\frac32+\frac12\right)\right)^2+\psi^{(1)}\left(\frac12\right)+\psi^{(1)}\left(\frac12+\frac32\right)\right]\\
&=\small\frac14\pi\left[\left(\left(\psi^{(0)}\left(\frac12\right)+\psi^{(0)}\left(2\right)\right)^2\right)+\psi^{(1)}\left(\frac12\right)+\psi^{(1)}\left(2\right)\right]\\
&=\small\frac14\pi\left[(-\gamma-2\log 2-1+\gamma)^2+\frac{\pi^2}2-\frac{\pi^2}6+1\right]\\
&=\small\frac14\pi\left[4\log^2 2+4\log 2+\frac{\pi^2}3+2\right]\\
&=\small -\pi\log 2+\pi\log^2 2-\frac{\pi}2+\frac{\pi^3}{12}
\end{align}$$
Form hereon the desired equality follows
$$\small\therefore~\int_0^{\frac{\pi}2}4\cos^2(x)\log^2(\cos x)\mathrm dx~=~-\pi\log 2+\pi\log^2 2-\frac{\pi}2+\frac{\pi^3}{12}$$

 A: You can first substitute $x=\arctan u$:
\begin{align}
4\int_0^{\pi/2} \cos^2 (x) \log^2 \cos x\,dx&=\int_0^\infty \frac{\log^2(1+u^2)}{(1+u^2)^2} \,du \\ 
&=\frac{\partial^2}{\partial^2 \beta}\Biggr|_{\beta=0} \int_0^\infty \frac{1}{(1+u^2)^{2-\beta}}\,du \\ 
&=\frac{\partial^2}{\partial^2 \beta}\Biggr|_{\beta=0} \frac{1}{\Gamma(2-\beta)}\int_0^\infty \int_0^\infty \nu^{1-\beta} e^{-\nu(1+u^2)}\,d\nu\, du \\
&=\frac{\partial^2}{\partial^2 \beta}\Biggr|_{\beta=0} \frac{1}{\Gamma(2-\beta)}\int_0^\infty e^{-\nu}\nu^{1-\beta} \int_0^\infty  e^{-\nu u^2}\,du\, d\nu \\
&= \frac{\partial^2}{\partial^2 \beta}\Biggr|_{\beta=0} \frac{\sqrt{\pi}}{\Gamma(2-\beta)}\int_0^\infty e^{-\nu}\nu^{1/2 -\beta}  d\nu \\
&=\frac{\partial^2}{\partial^2 \beta}\Biggr|_{\beta=0} \frac{\sqrt{\pi}\,\Gamma\left(3/2-\beta\right)}{\Gamma(2-\beta)}
\end{align} 
After doing the computations, we can find that this is equal to 
\begin{align}
\frac{\partial^2}{\partial^2 \beta}\Biggr|_{\beta=0} \frac{\sqrt{\pi}\,\Gamma\left(3/2-\beta\right)}{\Gamma(2-\beta)} &= \sqrt{\pi}\frac{\Gamma(3/2)}{\Gamma(2)}\left(\psi^{(0)}(3/2)^2-2\psi^{(0)}(2)\psi^{(0)}(3/2)+\psi^{(1)}(3/2)+\psi^{(0)}(2)^2-\psi^{(1)}(2)\right) \\ &=\frac{\pi}{4}\left(\frac{\pi^2}{3}+4\log^2-4\log 2 -2\right) \\
&=\frac{\pi^3}{12}+\pi\log^2 2 -\pi\log 2 -\frac{\pi}{2}
\end{align}
A: Let’s consider the general case.$$\mathfrak{I}=\int\limits_0^{\pi/2}\mathrm dx\,\cos^mx\sin^nx$$Make the substitution $x\mapsto\sin x$ to get that$$\mathfrak{I}=\int\limits_0^1\mathrm dx\, \left(1-x^2\right)^{(m+1)/2}x^n$$And now make the transformation $x\mapsto x^2$ so that$$\mathfrak{I}=\frac 12\int\limits_0^1\mathrm dx\,\left(1-x\right)^{(m-1)/2}x^{(n-1)/2}=\frac 12\operatorname{B}\left(\frac {m+1}2,\frac {n+1}2\right)$$Now substitute $n=0$ and differentiate with respect to $m$.
