What's the value of $\int_0^1\frac{1}{2y} \ln(y) \ln^2(1-y) \, dy$?

We start by introducing the integral $$ I(a,b)=\frac{1}{2}\int_0^1y^{a-1}(1-y)^b\,dy=\frac{1}{2}B(a,1+b), $$ where $B$ denotes the beta function. Note that this integral is singular at $a=0$ and $b=-1$. Since $\partial_a y^a=y^a\ln y$ we are led to calculate $$ \partial_{a,b,b}I(a,b)=\frac{1}{2}\int_0^1 y^{a-1}(1-y)^b\ln y\bigl(\ln(1-y)\bigr)^2\,dy $$ as $a$ and $b$ tend to $0$. We will below inser the "non-dangerous" point $b=0$. In other words, we want to calculate $$ \partial_{a,b,b}B(a,1+b)\mid_{a\to 0^+,b\to 0}. $$ When differentiating the beta function, polygammas appear. Indeed, $$ \begin{aligned} \partial_bB(a,1+b)&=B(a,1+b)\bigl(\psi_0(1+b)-\psi_0(1+a+b)\bigr)\\ \partial_{b,b}B(a,1+b)&=B(a,1+b)\Bigl(\bigl(\psi_0(1+b)-\psi_0(1+a+b)\bigr)^2 +\psi_1(1+b)-\psi_1(1+a+b)\Bigr). \end{aligned} $$ Next, we can actually insert $b=0$ before we differentiate with respect to $a$ and take the limit $a\to 0$. We should differentiate the function (here we have used the facts that $\psi_0(1)=-\gamma$ (Euler's constant) and that $\psi_1(1)=\pi^2/6$) $$ f(a)=B(a,1)\Bigl(\bigl(\gamma+\psi_0(1+a)\bigr)^2+\frac{\pi^2}{6}-\psi_1(1+a)\Bigr) $$ and calculate $\lim_{a\to 0^+}f'(a)$. We get that $$ \begin{aligned} f'(a)&=B(a,1)\bigl(\psi_0(a)-\psi_0(1+a)\bigr)\Bigl(\bigl(\gamma+\psi_0(1+a)\bigr)^2+\frac{\pi^2}{6}-\psi_1(1+a)\Bigr)\\ &\quad+B(a,1)\Bigl(2\bigl(\gamma+\psi_0(1+a)\bigr)\psi_1(1+a)-\psi_2(1+a)\Bigr) \end{aligned} $$ Next, we use the (non-obvious) expansions around $a=0$ $$ \begin{aligned} B(a,1)&=\frac{1}{a}+O(1)\\ \psi_0(a)&=-\frac{1}{a}-\gamma+O(a)\\ \psi_0(1+a)&=-\gamma+\frac{\pi^2}{6}a+O(a^2)\\ \psi_1(1+a)&=\frac{\pi^2}{6}+\psi_2(1)a+\frac{\pi^4}{30}a^2+O(a^3)\\ \psi_2(1+a)&=\psi_2(1)+\frac{\pi^4}{15}a+O(a^2). \end{aligned} $$ to find that, as $a\to0^+$, $$ \begin{aligned} f'(a)&\approx -\frac{1}{a^2}\Bigl(\bigl(\frac{\pi^2}{6}a\bigr)^2-\psi_2(1)a-\frac{\pi^4}{30}a^2\Bigr)+\frac{1}{a}\Bigl(2\frac{\pi^2}{6}a\frac{\pi^2}{6}-\psi_2(1)-\frac{\pi^4}{15}a\Bigr)+O(a)\\ &=-\frac{\pi^4}{180}+O(a) \end{aligned} $$ as $a\to 0^+$. We conclude that $$ \partial_{a,b,b}B(a,1+b)\mid_{a\to 0^+,b\to 0}=-\frac{\pi^4}{180}. $$ Finally, dividing by $2$ (remember, we had a one-half in front of the beta function in the beginning), we get that $$ \int_0^1\frac{1}{2y}\ln y(\ln(1-y))^2\,dy=-\frac{1}{360}\pi^4. $$

There is a closed form antiderivative that can be found using repeated integration by parts: $$\begin{align}\int\frac{\ln y\cdot\ln^2(1-y)}{2y}dy=\frac{\ln^4(1-y)}8+\frac{\ln y\cdot\ln ^3(1-y)}6+\left(\frac{\pi^2}{12}-\frac{\ln ^2y}2\right)\cdot\ln ^2(1-y)\\ +\left[\left(\frac{\pi^2}6+\operatorname{Li}_2\left(\tfrac y{y-1}\right)\right)\cdot\ln y-\operatorname{Li}_3(1-y)-\operatorname{Li}_3\left(\tfrac y{y-1}\right)\right]\cdot\ln(1-y)\\ +\left[\vphantom{\Large|}\zeta(3)-\operatorname{Li}_3(1-y)\right]\cdot\ln y+\operatorname{Li}_4(1-y)-\operatorname{Li}_4(y)-\operatorname{Li}_4\left(\tfrac y{y-1}\right)\color{gray}{+C}\end{align}$$

Another one:

Write $\log^2(1-y)=\sum_{n,m=1}^{\infty}\frac{y^{m+n}}{mn}$

we obtain (exchanging summation and integration)

$$ 2I=\sum_{n,m=1}^{\infty}\frac{1}{mn}\int_0^1\log(y)y^{m+n-1}=\\ -\underbrace{\sum_{n,m=1}^{\infty}\frac{1}{mn}\frac{1}{(m+n)^2}}_{S} $$

the double sum can be tackeled by writing (i shamelessly benefit from this awesome answer)

$$ S=\sum_{n,m=1}^{\infty}\frac{1}{mn}\frac{1}{(m+n)^2}=\frac{1}{2}\sum_{n,m=1}^{\infty}\frac{1}{m^2n^2}-\frac{1}{2}\sum_{n,m=1}^{\infty}\frac{1}{(m+n)^2n^2}-\frac{1}{2}\sum_{n,m=1}^{\infty}\frac{1}{(m+n)^2m^2} $$

shifting arguments in the last two sums gives

$$ 2S=\sum_{n,m=1}^{\infty}\frac{1}{mn}\frac{1}{(m+n)^2}-\sum_{m=1,n<m}^{\infty}\frac{1}{n^2m^2}-\sum_{m=1,n>m}^{\infty}\frac{1}{n^2m^2} $$

which yields

$$ 2S=\sum_{n=m=1}^{\infty}\frac{1}{m^2n^2}=\zeta(4)=\frac{\pi^4}{90} $$

and therefore

$$ I=-S/2=-\frac{\pi^4}{360} $$