Let's state with $\psi^{(1)}$ the trigamma. Calculate the order:
$$\mathcal{S} = \sum_{n=1}^\infty (-1)^{n-1} (\psi^{(1)}(n))^2$$
(Cornel Ioan Valean)
I uploaded this question here https://math.stackexchange.com/q/2412715/1295548 and https://math.stackexchange.com/q/2411733/1295548 from my old account
$$
\psi_1(z) = -\int\limits_0^1 \frac{x^{z-1} \ln x}{1-x} \, dx
$$
\begin{align}
(\psi_1(z))^2 & =\left( -\int\limits_0^1 \frac{x^{z-1}\ln x}{1-x} \, dx \right) \left( -\int\limits_0^1 \frac{y^{z-1}\ln y}{1-y} \, dy \right) = \int\limits_0^1 \int\limits_0^1 \frac{(xy)^{z-1} \ln x\ln y}{(1-x)(1-y)} \, dx \, dy \\[8pt]
& \sum_{n=1}^{+\infty} (-1)^{n-1} (\psi_1(n))^2 = \sum_{n=1}^{+\infty} ( (\psi_1 (2n-1))^2-(\psi_1(2n))^2)
\end{align}
$$=\sum_{n=1}^{+\infty} \int\limits_0^1 \int\limits_0^1 \frac{( (xy)^{2n-2}-(xy)^{2n-1} ) \ln x\ln y}{(1-x) (1-y)} \, dx \, dy = \int\limits_0^1 \int\limits_0^1 \frac{(1-xy)\ln x\ln y}{xy(1-x)(1-y)}\cdot \frac{(xy)^2}{1-(xy)^2} \, dx \, dy $$
$$=\int\limits_0^1 \int\limits_0^1 \frac{\ln x\ln y}{(1-x) (1-y)} \cdot \frac{xy}{1+xy} \, dx \, dy = -\frac{1}{6} \int\limits_0^1 \frac{6 \operatorname{Li}_2 (-y)+\pi^2 y}{1-y^2} \ln y \, dy $$
$$=-\frac{\pi ^{2}}{6}\int\limits_0^1 \frac{y}{1-y^2}\ln y \, dy – \int\limits_0^1 \frac{\ln y\operatorname{Li}_2 (-y)}{1-y^2} \, dy = \frac{\pi^2}{6} \cdot \frac{\pi^2}{24}-\int\limits_0^1 \frac{\ln y \operatorname{Li}_2 (-y)}{1-y^2} \, dy $$
Best Answer
Let $J$ be the integral to be computed. The following solution was written so that the parallel "clean" formula or a related integral $J_1$ (with no minus sign under the blue dilogarithm) is also mentioned: $$ \tag{$*$} J_1 = \int_0^1\log y\; \color{blue}{\operatorname{Li}_2(y)}\;\frac1{1-y^2}\; dy =\bbox[yellow]{\ -\frac 1{192}\pi^4\ }\ . $$ We record formulas for its "pieces": $$ \begin{aligned} J_{11} &= \int_0^1\log y\; \color{blue}{\operatorname{Li}_2(y)}\;\frac1{1-y}\; dy =\bbox[yellow]{\ -\frac 1{120}\pi^4\ }\ , \\ J_{12} &= \int_0^1\log y\; \color{blue}{\operatorname{Li}_2(y)}\;\frac1{1+y}\; dy =\bbox[yellow]{\ -\frac 1{480}\pi^4\ }\ , \end{aligned} $$
Up to some point, integrals are manipulated using the "usual" properties of integrals. The solution is given up to referenced results. However, afterwards, also a symbolic point of view is taken, key words are multiple zeta values (MZV) and generalized polylogarithms.
We will use first relations from the wolfram dilog page to split the wanted $J$ into pieces. For instance the relation $(6)$ is useful: $$ \tag{6} \begin{aligned} \operatorname{Li}_2(-y) - \operatorname{Li}_2(1-y) +\frac 12 \operatorname{Li}_2(1-y^2) &=-\frac {\pi^2}{12} - \log(+y)\log(1+y)\ , \\ \operatorname{Li}_2(+y) - \operatorname{Li}_2(1+y) +\frac 12 \operatorname{Li}_2(1-y^2) &=-\frac {\pi^2}{12} - \log(-y)\log(1-y)\ . \end{aligned} $$ Now subtract to get rid of the dilog computed in $(1-y^2)$, we obtain a way to express the dilog of $-y$ in terms of the dilog in $y$, $1\pm y$, and arguably simpler products of logarithms. Our $y$ runs in $[0,1]$, so we replace tacitly $\log(-y)$ by $\log y$, and pass to the real part of all expressions involved. (Equalities in the real part only, so modulo $i\Bbb R$ are denoted by $\equiv$.) From the above we obtain: $$ \tag{$\dagger$} \operatorname{Li}_2(-y) \equiv \color{blue}{\operatorname{Li}_2(y)} + \color{brown}{\operatorname{Li}_2(1-y) - \operatorname{Li}_2(1+y)} + \color{green}{\log(y)\log\frac{1-y}{1+y}}\ . $$
Note that we may have used $(5)$ from loc. cit. instead, which apparently should be simpler, since it involves only two dilogs, so we replace in integral involving $\operatorname{Li}_2(-y)$ with an other one involving $\operatorname{Li}_2(1+y)$, and products of $\log$'s. Using $(6)$ there are three dilogs appearing! Why? It turns out that the integral corresponding to $\operatorname{Li}_2(y)$ is "clean", and the difference $\operatorname{Li}_2(1+y)-\operatorname{Li}_2(1-y)$ also leads to a "clean" integral. Details follow.
We will use $(5)$ from loc. cit. to isolate the real part from $$ \operatorname{Li}_2(1+y) = \underbrace{ -\operatorname{Li}_2(-y) +\frac{\pi^2}6 -\log (-y)\log(1+y)}_{\Re \operatorname{Li}_2(1+y)} \pm i\pi\log(1+y)\ ,\qquad y\in[0,1]\ , $$ where the monodromy is hidden in the imaginary part, and we try to avoid powers of $\operatorname{Li}_2(1+y)$ in calculus, to mix different branches. Below, taking the real part is pointing tacitly to this dilog of $(1+y)$, $y\in[0,1]$, where $1+y$ leaves the disk of convergence of the dilog. In particular, $\Re \operatorname{Li}_2(2) =\Re \operatorname{Li}_2(1+1) =-\operatorname{Li}_2(-1)+\frac {\pi^2}6 =\frac{\pi^2}4$.
For these simpler products, for the integrals that are obtained from it is convenient to also introduce the notation: $$ \Bbb J_{abc,def}= \Bbb J\binom{a\ b\ c}{d\ e\ f} = \int_0^1 \frac {\log^a(1-y)\log^b(y)\log^c(1+y)}{(1-y)^d(y)^e(1+y)^f}\; dy\ . $$ A similar notation is introduced in the arXiv the article [Au1], page 1, of Kam Cheong Au, as mentioned in the other answer of . We use two references, to his articles, the information can be also extracted from the Book [Z], Chapter 13:
In the notations from [Au1] we have
Considering the results from [Au1], uur stand point is that for powers $d,e,f$ between $0,1$, and small values of the "weight" $w=a+b+c+1$ ($w\le 4$) we can compute / look up in a table for $\Bbb J_{abc,def}$, and take the value as a present. Note that a partial fraction decomposition reduces the cases for $(d,e,f)$ among $110$, $101$, $011$, $111$ to the simpler $100$, $010$, $001$. In fact, we will soon do a similar split for $\frac 1{1-y^2}=\frac 12\left(\frac 1{1-y}+\frac 1{1+y}\right)$.
Then we have from $(\dagger)$ a linear dependence between the real part of the following integrals $J, J_1,J_2,J_3$: $$ \begin{aligned} J &= \int_0^1\log y\; \operatorname{Li}_2(-y)\;\frac1{1-y^2}\; dy\ ,\text{ the integral of interest,} \\ J_1 &= \int_0^1\log y\; \color{blue}{\operatorname{Li}_2(y)}\;\frac1{1-y^2}\; dy =\bbox[yellow]{\ -\frac 1{192}\pi^4\ }\ ,\text{ a clean value to be shown below,} \\[3mm] J_2 &= \Re \int_0^1\log y\; \color{brown}{\Big( \operatorname{Li}_2(1+y) - \operatorname{Li}_2(1-y) \Big)} \;\frac1{1-y^2}\; dy \\ &= \frac 12 \Re\int_0^1\log y\; \left( \frac{\operatorname{Li}_2(1+y)}{1+y} + \frac{\operatorname{Li}_2(1+y)}{1-y} - \frac{\operatorname{Li}_2(1-y)}{1+y} - \frac{\operatorname{Li}_2(1-y)}{1-y} \right)\; dy \\ &\equiv -\frac 12 \int_0^1 \Big(\Re\operatorname{Li}_2(1+y)\Big)'\Re\operatorname{Li}_2(1+y)\; dy - \frac 12 \int_0^1 \Big(\operatorname{Li}_2(1-y)\Big)'\operatorname{Li}_2(1-y)\; dy \\ &\qquad\qquad + \frac 12 \int_0^1 \underbrace{\frac{\log y}{1-y}} _{=(\ +\operatorname{Li}_2(1-y)\ )'} \Re\operatorname{Li}_2(1+y)\; dy - \frac 12 \underbrace{\int_0^1\frac{\log y}{1+y}} _{\sim (\ -\Re\operatorname{Li}_2(1+y)\ )'} \operatorname{Li}_2(1-y)\; dy \\ &= -\frac 14\Bigg[\ (\Re\operatorname{Li}_2(1+y))^2\ \Bigg]_0^1 -\frac 14\Bigg[\ \operatorname{Li}_2(1-y)^2\ \Bigg]_0^1 \\ &\qquad\qquad +\frac 12\Bigg[\ \operatorname{Li}_2(1-y)\Re\operatorname{Li}_2(1+y)\ \Bigg]_0^1 \\ &= -\frac 14\Big[\ (\Re \operatorname{Li}_2(2))^2 - \operatorname{Li}_2(1)^2\ \Big] -\frac 14\Big[\ 0-\operatorname{Li}_2(1)^2 \ \Big] +\frac 12\Big[ \ 0-\operatorname{Li}_2(1))^2 \ \Big] \\ &=-\frac 14\cdot (\Re \operatorname{Li}_2(2))^2 =\bbox[yellow]{\ -\frac {\pi^4}{64}\ }\ , \\[3mm] J_3 &=\int_0^1\log y\cdot \color{green}{\log y\log\frac{1-y}{1+y}}\cdot \frac1{1-y^2}\; dy \\ &=\frac 12\Bbb J_{120-021,100+001}\qquad\text{(with a linear split w.r.t. the symbolic indices)} \\ &= \frac{53}{1440}\pi^4 +\frac 16\pi^2\log^2 2 -4\operatorname{Li}_{\color{red}4}\left(\frac 12\right) -\frac 72\log(2)\,\zeta(3) -\frac 16\log^42\ . \end{aligned} $$ Then we obtain a formula for $J$, from $J=J_1-J_2+J_3$, $$ \bbox[lightyellow]{\qquad J= \frac{17}{360}\pi^4 +\frac 16\pi^2\log^2 2 -4\operatorname{Li}_{\color{red}4}\left(\frac 12\right) -\frac 72\log(2)\,\zeta(3) -\frac 16\log^42\ . \qquad} $$
MZV and symbolic computations. I will try to collect some facts, see the above references. The problem with MZV's (Multiple Zeta *Values) is that all articles go in the direction of using the (Hopf-)algebraic part of them,
there are not so many examples - so that the many MSE integrals that appear still look like each one is a world for itself. Dictionary and formulas follow.
We fix a ("chromatic") level $N$. For our needs, $N=1$ and $N=2$ are enough. We deal with iterated integrals, computed for functions $f,g,\dots,h$,
identified in notation with the right objects to be taken, which are forms $f(t)\; dt$, $g(t)\; dt$, ... , $h(t)\; dt$ - which are defined on some interval (or manifold parametrized interval) $\gamma$. In our case (the manifold is $\Bbb R$ and) $\gamma=[0,1]$. So i will write the integral limits $0,1$ below. $$ \operatorname{It}(f,g,\dots,h) =\int_0^1 f(t)\; dt\int_0^t g(u)\; du\dots \int_0^? h(v)\; dv\ , $$ where the question mark stays for the previously used integration variable. So we integrate $f\otimes g\otimes\dots\otimes h$ on the simplex $1\ge t\ge u\ge\dots\ge v\ge 0$. Write $$ \operatorname{It}_\gamma(f,g,\dots,h) $$ if needed in order to integrate on an other path, e.g. $[0,t]$ with a dynamic $t$. We are mainly interested in functions $f,g,\dots,h$ in a very small world of functions. Here are they:
level $N=1$, $\frac 1t$ and $\frac 1{1-t}$ - written as forms $\frac{dt}t$, $\frac{dt}{1-t}$,
level $N=2$, $\frac 1t$ and $\frac 1{1-t}$ and $\frac 1{-1-t}$ - written as forms $\frac{dt}t$, $\frac{dt}{1-t}$, $\frac{dt}{-1-t}$,
and with a general $N$ we associate $\mu$ the primitive root of unity which is $\exp$ of $2\pi i/N$, and the alphabet contains now $\frac 1t$ and all fractions $\frac 1{\mu^{-k}-t}=\frac{\mu^k}{1-\mu^k t}$.
Iterated integrals look then like (non-commutative) words in this alphabet.
Since the alphabet is so small, and the the characters in it so complicated, we really use established "simple letters". For level $1$, in a first decade the convention was to take $0$ for $1/t$ and $1$ for $1/(1-t)$. Then further levels came in, and an other notation was with $a$ and $b_0$ instead, as in the references [Au1], [Au2], [Z]. The index zero of $b_0$ points to the power $k=0$ of the unit $\mu$ to be taken. So for a general $N$ we have $b_0$, $b_1$, ... till $b_{N-1}$. My convention below is to use
$$ \begin{array}{|c|c|} \hline \text{Words in few symbols} & \text{Tuples of few functions}\\\hline A &\frac{dt}t\\\hline B=B_0 &\frac{dt}{1-t}\\\hline B'=B_1 &-\frac{dt}{1+t}\\\hline X=\text{Alphabet }\{A,B,C\}^* & \Omega=\text{Alphabet in }\frac{dt}t,\ \frac{dt}{1-t},\ -\frac{dt}{1+t}\\\hline X\text{-word} & \Omega\text{-word}\\\hline X\text{-word }w & \Omega\text{-word often also written }w\\\hline \int_0^1 w\in\Bbb R & \operatorname{It}(w)=\operatorname{It}_{[0,1]}(w)\in\Bbb R\\\hline \operatorname{Li}_w(t)=\operatorname{Li}(w,t)= \int_0^t w & \operatorname{It}_{[0,t]}(w)\\\hline \int_\gamma w & \operatorname{It}_\gamma(w)\\\hline \bbox[lightyellow]{\qquad \int_\gamma w\cdot \int_\gamma w' = \int_\gamma w\ ш\ w'\qquad} & \operatorname{It}_\gamma(w)\operatorname{It}_\gamma(w') = \operatorname{It}_\gamma(w \ ш\ w')\\\hline ш\text{ is the shuffle product of words} &\\\hline \int_t^1 A^n & \frac 1{n!}(-1)^n\log^n t\\\hline \int_0^t B^n & \frac 1{n!}(-1)^n\log^n (1-t)=\frac 1{n!}\operatorname{Li}_1(t)^n\\\hline \int_0^t C^n & \frac 1{n!}(-1)^n\log^n (1+t)=\frac 1{n!}\operatorname{Li}_1(-t)^n\\\hline \operatorname{Li}(B,t) & \operatorname{Li}_1(t) =-\log(1-t)=\int_0^tB=\int_0^t\frac{dt}{1-t}\\\hline \operatorname{Li}(AB,t) & \operatorname{Li}_2(t) =\int_0^t\frac {\operatorname{Li}_1(t)}{t}\;dt \\ &=\int_0^t\frac {dt}t\int_0^t\frac {du}{1-u} = \operatorname{It}_{[0,t]}\left(\frac 1t,\frac1{1-t}\right) \\\hline \operatorname{Li}(AAB,t) & \operatorname{Li}_3(t) =\int_0^t\frac {\operatorname{Li}_2(t)}{t}\;dt = \operatorname{It}_{[0,t]}\left(\frac 1t,\frac 1t,\frac1{1-t}\right) \\\hline \operatorname{Li}(\underbrace{AA\dots AB}_{n\text{ letters}},t) & \operatorname{Li}_n(t) \\\hline Z(s) = Z_0(s)=\underbrace{AA\dots AB}_{s\text{ letters}} &\\ Z'(s) = Z_1(s) = \underbrace{AA\dots AB'}_{s\text{ letters}} &\\\hline Z_j(s)= \underbrace{AA\dots AB_j}_{s\text{ letters}} &\\\hline s =(s_1,s_2,\dots,s_d)\text{ poly-index} & s =(s_1,s_2,\dots,s_d)\in\Bbb N_{>0}^d\\\hline j =(j_1,j_2,\dots,j_d)\text{ poly-index} & j =(j_1,j_2,\dots,j_d)\in(\Bbb Z/N)^d\\\hline \Delta j=(j_1-0, j_2-j_1,\dots,j_d-j_{d-1}) & \\\hline x =(x_1,x_2,\dots,x_d)\text{ poly-argument} & x =(x_1,x_2,\dots,x_d)\\\hline \mu =(\mu,\mu,\dots,\mu)\text{ as poly-argument} & \mu =(\mu,\mu,\dots,\mu)\\\hline \prod Z(s,j) =\prod_j Z_j(s)\\ =Z_{j_1}(s_1)\dots Z_{j_d}(s_d) &\\\hline & \displaystyle \zeta(s) = L(s)=L(s,1)\\ & \displaystyle =\sum_{n_1>\dots >n_d\ge 1}\frac 1{n_1^{s_1}\dots n_d^{s_d}}\\\hline & \displaystyle L(s,x)=\sum_{n_1>\dots >n_d\ge 1}\frac {x_1^{n_1}\dots x_d^{n_d}}{n_1^{s_1}\dots n_d^{s_d}}\\\hline \int_0^1\prod Z(s,j) & L(s, \mu^{\Delta j})\\\hline \end{array} $$
We are computing integrals in this answer. Some part of the manipulations is purely algebraic, e.g. hidden in the shuffle algebra. Making the computations on the symbolic side is easier. Often, a MZV is obtained, and relations in the MZV-algebra are finally used to obtain the "usual" (even) powers of $\pi$, and zeta-values of odd arguments.
There are some convergence issues, for some words, but we neglect them here.
Let us illustrate the above with some computations. Words will occur, corresponding to functions/forms to be considered inside an iterated integral.
The iterated integrals and the $L$-values are extended to follow this rule.
The computation of $J_1$. The following part $J_{11}$ of $J_1$ is simpler, because only $y$ and $1-y$ appear, so the symbolic part involve only $A,B$: $$ \begin{aligned} J_{11} &= \int_0^1\log y\; \operatorname{Li}_2(y)\;\frac1{1-y}\; dy \\ &\overset{\text{up to }\pm} = \int_0^1\frac {dy}{1-y} \left(\int_y^1\frac{du}{u}\right) \left(\int_0^y\frac{du}{u}\int_u^1\frac {dv}{1-v}\right) \\ &= \int_0^1 B\int_0^t A\ ш\ AB \\ &= \int_0^1 B( \color{blue}{A}\ ш\ AB ) \\ &= \int_0^1 ( B\color{blue}{A}AB + BA\color{blue}{A}B + BAB\color{blue}{A} ) =2\int_0^1\color{brown}{B}AAB+\int_0^1 \color{brown}{B}ABA \\ &\qquad\text{ move $\color{brown}{B}$ from the first place behind}\\ &= -2\int_0^1 (A\color{brown}{B}AB + AA\color{brown}{B}B + AAB\color{brown}{B}) -\int_0^1 (A\color{brown}{B}BA + AB\color{brown}{B}A+ ABA\color{brown}{B}) \\ &= -3\int_0^1ABAB -4\int_0^1AABB -2\int_0^1ABB\color{blue}{A} \\ &\qquad\text{ move $\color{blue}{A}$ from the last place to the left}\\ \\ &=-3\int_0^1ABAB -4\int_0^1AABB +2\int_0^1(AB\color{blue}{A}B + A\color{blue}{A}BB + \color{blue}{A}ABB) \\ &=-\int_0^1ABAB=-\int_0^1Z(2)Z(2) \\ &\overset{\text{up to }\pm}=\int_0^1 Z(2,2)=\zeta(2,2) \\ &\qquad\text{ recall $\zeta(k)\zeta(n)=\zeta(k,n)+\zeta(n,k)+\zeta(n+k)$, use $k=n=2$} \\ &=\frac 12(\zeta(2)^2-\zeta(4)) = \frac 1{120}\pi^4\ . \end{aligned} $$ Numerical check, pari/gp:
The computation of the remained part $J_{12}$ from $J_1$ follows, we have algebraically a similar situation, but together with $y$ and $1-y$ there appears also $1+y$, so the symbolic part involve $A,B$ and also the $2$-chromatic piece $B'$: $$ \begin{aligned} J_{12} &= \int_0^1\log y\; \operatorname{Li}_2(y)\;\frac1{1+y}\; dy \\ &\overset{\text{up to }\pm} = \int_0^1\frac {dy}{1+y} \left(\int_y^1\frac{du}{u}\right) \left(\int_0^y\frac{du}{u}\int_u^1\frac {dv}{1-v}\right) \\ &= \int_0^1 B'\int_0^t A\ ш\ AB \\ &= \int_0^1 B'( \color{blue}{A}\ ш\ AB ) \\ &= \int_0^1 ( B'\color{blue}{A}AB + B'A\color{blue}{A}B + B'AB\color{blue}{A} ) =2\int_0^1\color{brown}{B'}AAB+\int_0^1 \color{brown}{B'}ABA \\ &\qquad\text{ move $\color{brown}{B'}$ from the first place behind}\\ &= -2\int_0^1 (A\color{brown}{B'}AB + AA\color{brown}{B'}B + AAB\color{brown}{B'}) -\int_0^1 (A\color{brown}{B'}BA + AB\color{brown}{B'}A+ ABA\color{brown}{B'}) \\ &= -2\int AB'AB-2\int_0^1AA(BB'+B'B) -\int_0^1(AB'B+ABB')\color{blue}{A} -\int_0^1ABAB' \\ &\qquad\text{ move $\color{blue}{A}$ from the last place to the left} \\ &= -2\int AB'AB-2\int_0^1AA(BB'+B'B) +\int_0^1\color{blue}{A}A(B'B+BB') \\ &\qquad +\int_0^1A\color{blue}{A}(B'B+BB') +\int_0^1AB'\color{blue}{A}B +\int_0^1AB\color{blue}{A}B' \\ &=-\int_0^1AB'AB=-\int_0^1Z_1(2)Z_0(2) \\ &\overset{\text{up to }\pm}=\int_0^1 Z_{(1,0)}(2,2)=L_{(1,1)}(2,2) \\ &=\zeta(\bar 2,\bar 2) \\ &\qquad\text{ use $\zeta(\bar k)\zeta(\bar n)=\zeta(\bar k,\bar n)+\zeta(\bar n,\bar k)+\zeta(n+k)$, $n+k$ even, for $n=k=2$} \\ &\qquad\text{ recall $\zeta(\bar 2)=-\frac 1{1^2}+\frac 1{2^2}-\frac 1{3^2}+\frac 1{4^2}\pm \dots=-\left(1-\frac 2{2^2}\right)\zeta(2)=-\frac1{12}\pi^2$, } \\ &=\frac 12(\zeta(\bar 2)^2-\zeta(4)) = \frac 12\pi^4\left(\frac 1{12^2}-\frac 1{90}\right) = -\frac 1{480}\pi^4\ . \end{aligned} $$ Numerical check, pari/gp:
I have to stop here, typing kills my energy, and the MO-web-interpreter hangs. It remains to show the formula for $J_3$, again the same symbolic computations apply. Similar examples are in [Au1], [Au2]. Initially, i considered the case to be simpler to put down on the paper, but - as the result shows, there must be a mess in weight four words that show up. If a particular aspect is of interest, i will fill in. The notation $\Bbb J$ was introduced to conquer better the situation, but...
Later EDIT: The formula for $J_3$ can be obtained in a similar symbolic framework, and then checked in pari/gp, i will do this now in order to answer some comments. The fomula for $J_3$ is $$ J_3 = \int_0^1\log^2 y\;\log\frac{1-y}{1+y}\cdot \frac1{1-y^2}\; dy $$ so is a linear combination of integrals of $\frac 1{1\pm y}\log^2y\frac 1{1\pm y}$. Working symbolically, $\frac 1{2!}\log^2y$ corresponds to $A$, and $\frac 1{1-y}$, $-\frac 1{1+y}$ respectively to $B=B_0$, $B'=B_1$. We use the indices $j,k$ in the $B$-symbols to cover all four cases. We use $\mp_j$ for $(-1)^{j-1}$, so $\mp_0=\mp$. Then: $$ \begin{aligned} J_{3jk}&:=\int_0^1\frac1{1-\mp_j y}\log^2 y\;\log(1\mp_k y)\; dy \\ &\overset{\text{up to }\pm}= \int_0^1 B_j(2A^2\ ш\ B_k) =2\int_0^1B_j(B_kAA +AB_kA+AAB_k)\qquad\text{(right-$A^2$-moves)} \\ &=-\int_0^1(AB_jB_kA+B_jAB_kA)+2\int_0^1B_jAB_kA+2\int_0^1B_jAAB_k \\ &=-\int_0^1AB_jB_kA+\int_0^1B_jAB_kA+2\int_0^1B_jAAB_k\qquad\text{(right-$A$-moves)} \\ &=+\int_0^1(2A^2B_jB_kA + AB_jAB_k) - \int_0^1(AB_jAB_k + 2B_jAAB_k) +2\int_0^1B_jAAB_k \\ &=+2\int_0^1A^2B_jB_k=2\int_0^1Z_j(3)Z_k(1) \\ &= 2L_{3, 1}(\ (-1)^j\ ,\ (-1)^{j+k}\ )\ , \\ &\qquad\qquad\text{ and we use the values:} \\ &\qquad L_{3,1}(1,1) = \zeta(3,1)=\frac 14\zeta(4)=\frac 1{360}\pi^4\ ,\\ &\qquad L_{3,1}(1,-1) = \zeta(3,\bar 1)\ ,\\ &\qquad L_{3,1}(-1,1) = \zeta(\bar 3,1)\ ,\\ &\qquad L_{3,1}(-1,-1) = \zeta(\bar 3,\bar 1)\ ,\\ &\qquad\qquad\text{ and we have:}\\ \frac 12J_{300} &=-\zeta(3,1)=-\frac 1{360}\pi^4\ , \\ \frac 12 J_{310} &= +\zeta(\bar 3,\bar 1) = \frac 1{180}\pi^4 +\frac 1{12}\pi^2\log^2 2 - 2\operatorname{Li}_4\left(\frac 12\right) - \frac 1{12}\log^4 2\ , \\ -\frac 12J_{301} &=+\zeta(3,\bar 1)= \frac{19}{1440}\pi^4 -\frac 74\log 2\;\zeta(3) \ , \\ -\frac 12 J_{311} &= -\zeta(\bar 3,1) = \frac 1{48}\pi^4 +\frac 1{12}\pi^2\log^2 2 - 2\operatorname{Li}_4\left(\frac 12\right) -\frac 74\log2\;\zeta(3) - \frac 1{12}\log^4 2\ , \end{aligned} $$ and this gives finally the claimed value for $J_3$, $$ J_3=\frac 12(J_{300} +J_{310}-J_{301}-J_{311}) =\frac {53}{1440}\pi^4 +\frac 16\pi^2\log^2 2 - 4\operatorname{Li}_4\left(\frac 12\right) -\frac 72\log 2\;\zeta(3) - \frac 16\log^4 2 \ . $$
Note: The sum $\zeta(3,\bar 1)+\sum(\bar 3,\bar 1)$ can be also extracted from [Au2], Example 3.2, where $-\zeta(3,\bar 1)-\sum(\bar 3,\bar 1)$ is given (with no further reference). The value $\zeta(3,1)=\frac 14\zeta(4)$ can be obtained from $\zeta(2)^2=\int_0^1AB\cdot\int_0^1 AB=\int_0^1AB\ ш\ AB =2\int_0^1ABAB + 4\int_0^1AABB=2\int_0^1Z(2)Z(2)+4\int_0^1Z(3)Z(1)=2\zeta(2,2)+\zeta(3,1)$. The book [Z] has in Appendix, Euler sums for lower weights, page 537, also formulas relating these MZVs. All weight four alternating MZV are listed.