Proving $\int_{0}^{\pi/2} \ln \sin x ~ \ln \cos x~ dx=-\frac{\pi^3}{48}+\frac{\pi}{2}\ln^22$ Interestingly the integral $$I=\int_{0}^{\pi/2} \ln \sin x ~ \ln \cos x~ dx~~~~~~~~~(1)$$
is doable by hand by using Fourier series: $$ \ln \sin x=-\sum_{j=1}^{\infty} \frac{\cos 2j x}{j}-\ln 2,\quad \ln \cos x=\sum_{k=1}^{\infty} (-1)^{k+1} \cfrac{\cos 2kx}{k}-\ln 2, ~x\in [0,\pi/2].~~~~~~~~~~(2)$$
$$I=\sum_{j=1}^{\infty} \sum_{k=1}^{\infty} (-1)^{k}\int_{0}^{\pi/2} \frac{\cos 2jx \cos 2kx}{jk} dx-\ln 2 \int_{0}^{\pi/2}\left (\sum_{j=1}^{\infty} \frac{\cos 2jx}{j}+\sum_{k=1}^{\infty} (-1)^k \frac{\cos 2kx}{k}\right) dx+\frac{\pi}{2}\ln^22~~~~~~~~~~~(3)$$
The second and third integrals vanish, then
$$\implies I=\frac{1}{2} \sum_{j=1}^{\infty} \sum_{k=1}^{\infty} (-1)^k \int_{0}^{\pi/2}\frac{\cos2(j+k)x+\cos2(j-k)x}{jk} dx+\frac{\pi}{2}\ln^22~~~~~~~~~~~~~~~~(4)$$
$$\implies I=\sum_{j=1}^{\infty} \sum_{k=1}^{\infty}(-1)^k \frac{\sin2(j-k)\pi}{4jk(j-k)}+\frac{\pi}{2}\ln^22~~~~~~~~~~~~~~~~(5)$$
Taking limit $(j-k)\to 0$, we get
$$I=\frac{\pi}{4} \sum_{k=1}^{\infty} \frac{(-1)^k }{k^2}+\frac{\pi}{2}\ln^22=
-\frac{\pi^3}{48}+\frac{\pi}{2}\ln^22.~~~~~~~~~~~~~~~(6)$$
The question is: How else we can get this interesting integral (1)?
 A: To the nice solutions above - to add an another possible way to evaluate the integral
$$I=\int_{0}^{\pi/2} \ln (\sin x) ~ \ln (\cos x)~ dx$$
Let's consider $I_1=\int_{0}^{\pi/2} \Bigl(\ln (\sin x) - \ln (\cos x)\Bigr)^2~ dx$ and $I_2=\int_{0}^{\pi/2} \Bigl(\ln (\sin x) + \ln (\cos x)\Bigr)^2~ dx$
Using the symmetry of the integrand $I_1=\int_{0}^{\pi/2} \Bigl(2\ln^2 (\sin x) - 2\ln (\cos x)\ln (\cos x)\Bigr)~ dx$.
On the other hand
$$I_1=\int_{0}^{\pi/2} \ln^2(\tan x) \,dx=\int_{0}^{\pi/2} \ln^2(\tan x)\frac{d(\tan x)}{1+\tan^2 x}=\int_0^\infty\frac{\ln^2 y}{1+y^2} dy$$
$$=\frac{1}{8}\int_0^\infty\frac{\ln^2t \,dt}{\sqrt t(1+t)} dt=\frac{1}{8}\frac{\partial^2}{\partial s^2}|_{s=0}\int_0^\infty\frac{t^{s-1/2}}{1+t}dt$$
Making the change $\frac{1}{1+t}=x$
$$I_1=\frac{1}{8}\frac{\partial^2}{\partial s^2}|_{s=0}\int_0^1(1-x)^{s-1/2}x^{-s-1/2}dx=\frac{1}{8}\frac{\partial^2}{\partial s^2}|_{s=0}B\Bigl(s+\frac{1}{2};\frac{1}{2}-s\Bigr)$$
$$=\frac{1}{8}\frac{\partial^2}{\partial s^2}|_{s=0}\frac{\Gamma\Bigl(\frac{1}{2}-s\Bigr)\Gamma\Bigl(\frac{1}{2}-s\Bigr)}{\Gamma(1)}=\frac{1}{8}\frac{\partial^2}{\partial s^2}|_{s=0}\frac{\pi}{\sin\Bigl(\frac{\pi}{2}+\pi s\Bigr)}=\frac{\pi^3}{8}$$
Therefore
$$I_1=2\int_{0}^{\pi/2} \ln^2 (\sin x)\,dx - 2\int_{0}^{\pi/2}\ln (\cos x)\ln (\cos x)\,dx=\frac{\pi^3}{8}\quad \mathbf{(1)}$$
In the same fashion
$$I_2=2\int_{0}^{\pi/2} \ln^2 (\sin x)\,dx + 2\int_{0}^{\pi/2}\ln (\cos x)\ln (\cos x)\,dx=\int_{0}^{\pi/2} \ln^2 \Bigl(\frac{\sin2x}{2}\Bigr)dx$$
$$=\int_{0}^{\pi/2} \ln^2 \Bigl(\frac{\sin y}{2}\Bigr)dy=\int_{0}^{\pi/2} \ln^2 (\sin y)\,dy-2\ln2\int_{0}^{\pi/2} \ln (\sin y)\,dy+\frac{\pi}{2}\ln^22$$
Using the fact that $\int_{0}^{\pi/2} \ln (\sin y)\,dy=-\frac{\pi}{2}\ln2$ we get the second equation
$$\int_{0}^{\pi/2} \ln^2 (\sin x)\,dx + 2\int_{0}^{\pi/2}\ln (\cos x)\ln (\cos x)\,dx=\frac{3\pi}{2}\ln^22\quad \mathbf{(2)}$$
Combining (1) and (2)
$$I=\int_{0}^{\pi/2} \ln (\sin x) ~ \ln (\cos x)~ dx=\frac{\pi}{2}\ln^22-\frac{\pi^3}{48}$$
$$\int_{0}^{\pi/2} \ln^2 (\sin x)\,dx=\int_{0}^{\pi/2} \ln^2 (\cos x)\,dx=\frac{\pi}{2}\ln^22+\frac{\pi^3}{24}$$
A: This integral is evaluated in this answer. The approach in that answer seems to be similar to the method used in your question.
In this answer, a general method is developed to handle integrals of this type. We can use it to evaluate this integral:
$$
\begin{align}
\int_0^{\pi/2}\log(\sin(x))\log(\cos(x))\,\mathrm{d}x
&=\frac{\partial_1\partial_2\mathrm{B}\!\left(\frac12,\frac12\right)}8\tag1\\
&=\frac18\partial_1(A_{0,1}B)\tag2\\
&=\frac18(A_{1,1}B+A_{0,1}A_{1,0}B)\tag3\\[3pt]
&=\frac\pi8\left(-\zeta(2)+4\log(2)^2\right)\tag4\\
&=-\frac{\pi^3}{48}+\frac\pi2\log(2)^2\tag5
\end{align}
$$
Explanation:
$(1)$: apply $(6)$ from the cited answer
$(2)$: apply $\text{(4b)}$ from the cited answer
$(3)$: apply $\text{(4c)}$ and $\text{(4a)}$ from the cited answer
$(4)$: apply $(5)$ from the cited answer
$(5)$: simplify
A: Solution without using beta function:
Let $a=\ln(\sin x)$ and $b=\ln(\cos x)$ in the algebraic identity
$$ab=\frac12a^2+\frac12b^2-\frac12(a-b)^2,$$
we have
$$\ln(\sin x)\ln(\cos x)=\frac12\ln^2(\sin x)+\frac12\ln^2(\cos x)-\frac12\ln^2(\tan x).$$
Integrate both sides from $x=0$ to $\pi/2$,
\begin{gather*}
\int_0^{\frac{\pi}{2}}\ln(\sin x)\ln(\cos x)\mathrm{d}x\\
=\frac12\int_0^{\frac{\pi}{2}}\ln^2(\sin x)\mathrm{d}x+\frac12\int_0^{\frac{\pi}{2}}\ln^2(\cos x)\mathrm{d}x-\frac12\int_0^{\frac{\pi}{2}}\ln^2(\tan x)\mathrm{d}x.
\end{gather*}
The second integral is equivalent to the first one by using the rule $\int_a^b f(x)\mathrm{d}x=\int_a^b f(a+b-x)\mathrm{d}x$. For the third integral, let $\tan x=y$, we have
$$\int_0^{\frac{\pi}{2}}\ln(\sin x)\ln(\cos x)\mathrm{d}x=\int_0^{\frac{\pi}{2}}\ln^2(\sin x)\mathrm{d}x-\frac12\int_0^\infty\frac{\ln^2(y)}{1+y^2}\mathrm{d}y.$$
First integral:
We have
$$\ln^2(2\sin x)=\ln^2(2)+2\ln(2)\ln(\sin x)+\ln^2(\sin x)$$
or
$$\ln^2(\sin x)=\ln^2(2\sin x)-2\ln(2)\ln(\sin x)-\ln^2(2).$$
Integrate both sides from $x=0$ to $\pi/2$,
$$\int_0^{\frac{\pi}{2}} \ln^2(\sin x)\mathrm{d}x=\int_0^{\frac{\pi}{2}}\ln^2(2\sin x)\mathrm{d}x-2\ln(2)\int_0^{\frac{\pi}{2}}\ln(\sin x)\mathrm{d}x-\int_0^{\frac{\pi}{2}}\ln^2(2)\mathrm{d}x.$$
The third integral is $\frac{\pi}{2}\ln^2(2)$ and the second integral is $-\frac{\pi}{2}\ln(2)$. For the first one, integrate both sides of the identity:
$$\ln^2(2\sin x)=\left(\frac{\pi}{2}-x\right)^2+2\sum_{n=1}^\infty\frac{H_{n-1}}{n}\cos(2nx)$$
from $x=0$ to $\pi/2$ then change the order of integration and summation,
\begin{gather*}
\int_0^{\frac{\pi}{2}}\ln^2(2\sin x)\mathrm{d}x=\int_0^{\frac{\pi}{2}}\left(\frac{\pi}{2}-x\right)^2\mathrm{d}x+2\sum_{n=1}^\infty\frac{H_{n-1}}{n}\int_0^{\frac{\pi}{2}}\cos(2nx)\mathrm{d}x\\
=-\frac13\left(\frac{\pi}{2}-x\right)^3\bigg|_0^{\frac{\pi}{2}}+2\sum_{n=1}^\infty\frac{H_{n-1}}{n}\cdot\frac{\sin(2nx)}{2n}\bigg|_0^{\frac{\pi}{2}}\\
=\frac{\pi^3}{24}+2\sum_{n=1}^\infty\frac{H_{n-1}}{n}\cdot\frac{\sin(n\pi)}{2n}\\
\{\text{the sum evaluates to $0$, since $\sin(n\pi)=0$ for integer $n$}\}\\
=\frac{\pi^3}{24}.
\end{gather*}
Therefore,
$$\int_0^{\frac{\pi}{2}} \ln^2(\sin x)\mathrm{d}x=\frac{\pi}{2}\ln^2(2)-\frac{\pi^3}{16}.$$
Second integral:
By the generalization:
$$ \int_0^\infty\frac{\ln^{2a}(x)}{1+x^2}\mathrm{d}x=2^{-2a-1}\pi\lim_{s\to \frac12}\frac{\mathrm{d}^{2a}}{\mathrm{d} s^{2a}}\csc(\pi s),$$
which follows from differentiating Euler's reflection formula, we have
$$ \int_0^\infty\frac{\ln^2(x)}{1+x^2}\mathrm{d}x=\frac{\pi^3}{8}.$$
Combining the two integrals, we reach
$$\int_0^{\frac{\pi}{2}} \ln(\sin x)\ln(\cos x)\mathrm{d}x=\frac{\pi}{2}\ln^2(2)-\frac{\pi^3}{48}.$$
