$\def\ZZ{\mathbb{Z}}\def\RR{\mathbb{R}}$I will evaluate the sum at the end of Jim Belk's post.
There is clearly a lot more going on here; I hope a number theorist will come and explain what is really going on.
I will be sloppy about convergence issues.
Set
$$S:= \sum_{k,n=1}^{\infty} \frac{(-1)^{k+n}}{k^2+4kn+n^2}.$$
We set $m = k+2n$, in order to diagonalize the quadratic form:
$$S = \sum_{0 < 2n < m} \frac{(-1)^{m-n}}{m^2 - 3 n^2}.$$
It turns out to be more elegant to work with
$$T : = \sum_{0 \leq 2n < m} \frac{(-1)^{m-n}}{m^2 - 3 n^2} = S + \sum_{m=1}^{\infty} \frac{(-1)^m}{m^2} = S - \frac{\pi^2}{12}.$$
Let $R$ be the ring $\ZZ[\sqrt{3}]$. Let $D \subset R$ be $\{ m+n\sqrt{3} : 0 \leq 2n < m \}$. So, simply as a matter of formal rewriting
$$T = \sum_{m+n \sqrt{3} \in D} \frac{(-1)^{m+n}}{m^2-3 n^2}.$$
We will now try to interpret each of these quantities in terms of the ring $R$. Recall that, for $m+n \sqrt{3} \in R$, the norm $N(m+n \sqrt{3})$ is $m^2 - 3 n^2$.
Define $R_{+} = \{ a \in R : N(a) > 0\}$ and $R_{-} = \{ a \in R : N(a) < 0 \}$.
For $m+n \sqrt{3} \in R$, define $\sigma_2(m+n \sqrt{3}) = (-1)^{m+n}$.
Let $\Gamma$ denote the unit group of $R$. Explicitly, $\Gamma$ is $\{ \pm 1 \} \times (2+\sqrt{3})^{\ZZ}$.
Note that multiplication by $\Gamma$ takes $R_{+}$ and $R_{-}$ to themselves. Here is the first miracle: $D$ is a fundamental domain for the action of $\Gamma$ on $R_{+}$! For example, multiplication by $2+\sqrt{3}$ maps the ray $\RR_{\geq 0} (1+0 \sqrt{3})$ bounding one side of $D$ to the ray $\RR_{\geq 0} (2+\sqrt{3})$ on the other side. Moreover, multiplication by units leaves $N( \ )$ and $\sigma_2(\ )$ unchanged.
So we can view the sum as
$$T = \sum_{a \in R_{+}/\Gamma} \frac{\sigma_2(a)}{N(a)}$$
where the sum means to pick one representative for each orbit of the $\Gamma$ action.
I'd rather sum over $R$ than $R_{+}$. Define $\chi_{\infty}$ to be $\pm 1$ on $R_{\pm}$. So we can rewrite
$$ T = \frac{1}{2} \sum_{a \in R_{\neq 0}/\Gamma} \frac{\sigma_2(a) \left( 1+\chi_{\infty}(a) \right)}{|N(a)|} = \frac{1}{2}\left( \sum_{a \in R_{\neq 0}/\Gamma} \frac{\sigma_2(a)}{|N(a)|} + \sum_{a \in R_{\neq 0}/\Gamma} \frac{\sigma_2(a) \chi_{\infty}(a)}{|N(a)|} \right) = : \frac{1}{2} (U+V).$$
Now, $a$ and $b$ in $R$ are in the same $\Gamma$ orbit if and only if the ideals $(a)$ and $(b)$ are equal. So the above sums are running over all principal ideals of $R$.
Moreover, $R$ is a PID! And, finally, for a principal ideal $I = (a)$, we have $N(I) = |N(a)|$. So we can write:
$$U = \sum_{I \subset R \ \mbox{an ideal}} \frac{\sigma_2(I)}{N(I)} \quad V = \sum_{I \subset R \ \mbox{an ideal}} \frac{\sigma_2(I) \chi_{\infty}(I)}{N(I)} .$$
We set
$$U(s) = \sum_{I \subset R \ \mbox{an ideal}} \frac{\sigma_2(I)}{N(I)^s} \quad V(s) = \sum_{I \subset R \ \mbox{an ideal}} \frac{\sigma_2(I) \chi_{\infty}(I)}{N(I)^s} .$$
Note that $\sigma_2(a)$ is $1$ if and only if $1+\sqrt{3}$ divides $a$.
So we have the Euler factorization
$$U(s) = \left(-1 + 2^{-s} + 2^{-2s} + 2^{-3s} + \cdots \right) \prod_{\pi \neq (1+\sqrt{3})} \frac{1}{1-N(\pi)^{-s}} =$$
$$\frac{-1+2^{-s}+2^{-2s}+ 2^{-3s} + \cdots}{1+2^{-s}+2^{-2s}+2^{-3s}+\cdots} \prod_{\pi} \frac{1}{1-N(\pi)^{-s}} =(-1+2\cdot 2^{-s}) \zeta_R(s).$$
where $\pi$ runs over prime ideals of $R$.
So
$$\lim_{s \to 1^{+}} U(s) = \lim_{s \to 1^{+}} \frac{-1+2 \cdot 2^{-s}}{s-1} \lim_{s \to 1^{+}} (s-1) \zeta_R(s) = - \log 2 \lim_{s \to 1^{+}} (s-1) \zeta_R(s).$$
From the class number formula, this last limit is
$$\frac{2^2 \log (2+\sqrt{3})}{2 \cdot \sqrt{12}}.$$
So
$$U(1) = - \frac{\log 2 \log(2+\sqrt{3})}{\sqrt{3}}.$$
We now run the same trick with $V$. Again, we start with the Euler product:
$$V(s) = \left(-1 - 2^{-s} + 2^{-2s} - 2^{-3s} + 2^{-4s} - \cdots \right) \prod_{\pi \neq (1+\sqrt{3})} \frac{1}{1-\chi_{\infty}(\pi) N(\pi)^{-s}} =$$
$$\frac{-1 - 2^{-s} + 2^{-2s} - 2^{-3s} +2^{-4s} - \cdots}{1-2^{-s} + 2^{-2s} - 2^{-3s} + 2^{-4s} - \cdots} \prod_{\pi} \frac{1}{1-\chi_{\infty}(\pi) N(\pi)^{-s}} =
\left( -1 - 2^{1-s} \right) L(s, \chi_{\infty}).$$
So
$$V = V(1) = - 2 L(1, \chi_{\infty}).$$
We now must evaluate $L(1, \chi_{\infty})$. This $L$-function is defined by a sum over the ideals of $R$; we will rewrite it in terms of $L$-functions for $\ZZ$.
Let $p \geq 5$ be prime. I claim that $p$ splits in $R$ if and only if $p \equiv \pm 1 \bmod 12$. For such a $p$, if we write $p = \pi \bar{\pi}$, I claim that $N(\pi) = N(\bar{\pi}) = p$ if $p \equiv 1 \bmod 12$ and $N(\pi) = N(\bar{\pi}) = -p$ if $p \equiv -1 \bmod 12$. Proof Sketch: The prime $p$ splits in $R$ if and only if $\left( \frac{3}{p} \right) = 1$, which is easily computed to be equivalent to $p \equiv \pm 1 \bmod 12$. Since $R$ is a PID, we can write such a $p$ as $\pi \bar{\pi}$ for $\pi$ and $\bar{\pi} \in R$. Then $N(\pi) = N(\bar{\pi}) = \pm p$. But, for $a=m+n \sqrt{3} \in R$, we have $N(a) = m^2-3 n^2 \not \equiv -1 \mod 3$. So only one of the two possibilities for $\pm p$ can occur. $\square$
So $L(1, \chi_{\infty})$ is
$$\left(1+2^{-1} \right)^{-1} \left( 1+3^{-1} \right)^{-1} \prod_{p \equiv 1 \bmod 12} \left( 1- p^{-1} \right)^{-2} \prod_{p \equiv -1 \bmod 12} \left( 1+ p^{-1} \right)^{-2} \prod_{p \equiv \pm 5 \bmod 12} \left( 1- p^{-2} \right)^{-1}.$$
Define $\chi_4(n)$ to be $0$ if $n$ is even, $1$ if $n \equiv 1 \bmod 4$ and $-1$ is $n \equiv -1 \bmod 4$. Define $\chi_3(n)$ to be $1$, $-1$ or $0$ according to whether $n \equiv 1$, $2$ or $0 \bmod 3$. Then we have
$$L(1, \chi_{\infty}) = \prod_p \left( 1- \chi_4(p) p^{-1} \right)^{-1} \prod_p \left( 1-\chi_3(p) p^{-1} \right)^{-1} = \left( \sum_{n=1}^{\infty} \frac{\chi_4(n)}{n} \right) \left( \sum_{n=1}^{\infty} \frac{\chi_3(n)}{n} \right) .$$
The sum $\sum \chi_4(n) = 1-1/3+1/5-1/7 + \cdots$ is well known to be $\pi/4$. The sum $\sum \chi_3(n)/n$ is only slightly less well known; it is $\pi/(3 \sqrt{3})$.
So I get $L(1,\chi_{\infty}) = \frac{\pi}{4} \frac{\pi}{3 \sqrt{3}} = \frac{\pi^2}{12 \sqrt{3}}$ and $V = - \frac{\pi^2}{6 \sqrt{3}}$.
Putting it all together,
$$T = \frac{1}{2} \left( - \frac{\log 2 \log(2+\sqrt{3})}{\sqrt{3}} - \frac{\pi^2}{6 \sqrt{3}} \right) = - \frac{\log 2 \log(2+\sqrt{3})}{2 \sqrt{3}} - \frac{\pi^2}{12 \sqrt{3}} $$
$$S = \frac{\pi^2}{12} - \frac{\log 2 \log(2+\sqrt{3})}{2 \sqrt{3}} - \frac{\pi^2}{12 \sqrt{3}} .$$
The original integral is
$$\frac{1}{2} \left( \log(2)^2 -2 \sqrt{3} S \right) = \frac{\log(2)^2}{2} - \frac{\pi^2 \sqrt{3}}{12} + \frac{\log 2 \log(2+\sqrt{3})}{2} + \frac{\pi^2}{12}$$
$$=\frac{\pi^2(1-\sqrt{3})}{12} + \log(2) \frac{\log(4 + 2 \sqrt{3})}{2} = \frac{\pi^2(1-\sqrt{3})}{12} + \log(2) \log(1+\sqrt{3})$$
as desired.
Let the considered integral be denoted by $I$. Our starting point is to reduce the number of logarithms of different arguments in the integrand. Thus, using the fact
that $6ab^2=(a+b)^3-2a^3+(a-b)^3$ we obtain
\begin{align*}
6I&=\underbrace{\int_0^1\frac{\log x}{x}\log^3(1-x^2)dx}_{x\leftarrow \sqrt{u}}
-2\int_0^1\frac{\log x}{x}\log^3(1-x)dx+
\underbrace{\int_0^1\frac{\log x}{x}\log^3\left(\frac{1-x}{1+x}\right)dx}_{x\leftarrow\frac{1-u}{1+u}}
\\
&= \frac{1}{4}\int_0^1\frac{\log u}{u}\log^3(1-u)du
-2\int_0^1\frac{\log x}{x}\log^3(1-x)dx+
2\int_0^1\frac{\log^3 u}{1-u^2}\log\left(\frac{1-u}{1+u}\right)du \\
&=\frac{-7}{4}\int_0^1\frac{\log x}{x}\log^3(1-x)dx+
2\int_0^1\frac{\log^3 x}{1-x^2}\log\left(\frac{1-x}{1+x}\right)dx \\
\end{align*}
Thus,
$$
I=-\frac{7}{24}J+\frac{1}{3}K \tag{1}
$$
with
$$
J=\int_0^1\frac{\log(1- x)}{1-x}\log^3x\,dx,\quad
K=\int_0^1\frac{1}{1-x^2}\log\left(\frac{1-x}{1+x}\right)\log^3 x\,dx
\tag{2}
$$
Now, for $x\in(-1,1)$, we have
$$
\frac{\log(1-x)}{1-x}=-\left(\sum_{k=0}^\infty x^k\right)\left(\sum_{k=1}^\infty \frac{x^k}{k}\right)
=-\sum_{n=1}^\infty H_nx^n
$$
where $H_n=\sum_{k=1}^n1/k$ is the $n$-th Harmonic number. Similarly,
$$
\frac{1}{1-x^2}\log\left(\frac{1-x}{1+x}\right)
=-\left(\sum_{k=0}^\infty x^{2k}\right)\left(\sum_{k=1}^\infty \frac{-2x^{2k-1}}{2k-1}\right)
=-2\sum_{n=1}^\infty \widetilde{H}_nx^{2n-1}
$$
where
$$\widetilde{H}_n=\sum_{k=1}^n\frac{1}{2k-1}
=H_{2n}-\frac{1}{2}H_n$$
Thus, using the fact that $\int_0^1x^n(-\log x)^3\,dx=-6/(n+1)^4$ we conclude that
\begin{align*}
J&=\sum_{n=1}^\infty H_n\int_0^1x^n(-\log x)^3\,dx=6\sum_{n=1}^\infty
\frac{H_n}{(n+1)^4}\\
&=6\sum_{n=0}^\infty
\frac{1}{(n+1)^4}\left(H_{n+1}-\frac{1}{n+1}\right)=6A-6\zeta(5)\tag{3}
\end{align*}
with
\begin{equation*}
A=\sum_{n=1}^\infty\frac{H_n}{n^4}\tag{4}
\end{equation*}
Similarly,
\begin{align*}
K&=\sum_{n=1}^\infty (2H_{2n}-H_n)\int_0^1x^{2n-1}(-\log x)^3\,dx=6\sum_{n=1}^\infty\frac{2H_{2n}-H_n}{(2n)^4}\\
&=6\sum_{n=1}^\infty\frac{(1+(-1)^{2n})H_{2n}}{(2n)^4}
-\frac{3}{8}\sum_{n=1}^\infty\frac{H_n}{n^4}\\
&=6\sum_{n=1}^\infty\frac{(1+(-1)^{n})H_{n}}{n^4}
-\frac{3}{8}\sum_{n=1}^\infty\frac{H_n}{n^4}\\
&=\frac{45}{8}\sum_{n=1}^\infty\frac{H_{n}}{n^4}
+6\sum_{n=1}^\infty\frac{(-1)^{n}H_{n}}{n^4}=\frac{45}{8}A+6B\tag{5}
\end{align*}
with
\begin{equation*}
B=\sum_{n=1}^\infty\frac{(-1)^nH_n}{n^4}\tag{6}
\end{equation*}
Combining (1) with (3) and (5) we find that
\begin{equation*}
I=\frac{1}{8}A+2B+\frac{7}{4}\zeta(5)\tag{7}
\end{equation*}
Now, sums $A$ and $B$ are known (see here ) (in a more general setting), and
we have
$$
A=3\zeta(5)-\zeta(2)\zeta(3),\qquad
B=-\frac{59}{32}\zeta(5)+\frac{1}{2}\zeta(2)\zeta(3)
$$
Thus
$$
I=\frac{7}{8}\zeta(2)\zeta(3)-\frac{25}{16}\zeta(5).
$$
which is the announced result.$\qquad\square$
Best Answer
I tried substitutions and the differentiation w.r.t a paramater trick like the other posters. Another partial result, or a trail of breadcrumbs to follow, is the following. We try a series expansion, $$ \frac{\log\left(1-x^4\right)}{1+x^2} = \displaystyle \sum_{k=1}^{\infty} x^{4k}\left(x^{2} -1\right)H_k, $$ where $H_k$ are the Harmonic numbers. Then \begin{align} \int_0^1 \frac{\log x \log \left(1-x^4 \right)}{1+x^2}\ \mathrm{d}x &=\displaystyle \sum_{k=1}^{\infty}\, H_k\int_0^1 x^{4k}\left(x^{2} -1\right)\log x \ \mathrm{d}x \\ &=\displaystyle \sum_{k=1}^{\infty} \, \frac{H_k}{(4k+1)^2}-\displaystyle \sum_{k=1}^{\infty} \, \frac{H_k}{(4k+3)^2}. \end{align} These sums look very similar to the ones evaluated in this post, in which they are transformed into alternating sums. Using the same techniques, or perhaps working back from the answers, we can hopefully show that $$ \displaystyle \sum_{k=1}^{\infty} \, \frac{H_k}{(4k+1)^2} = -G\left(\frac{\pi}{4}+\frac{\log 8}{2} \right) +\frac{7}{4}\zeta(3) +\frac{\pi^3}{32} - \frac{\pi^2}{16}\log 8, $$ $$ \displaystyle \sum_{k=1}^{\infty} \, \frac{H_k}{(4k+3)^2} = -G\left(\frac{\pi}{4}-\frac{\log 8}{2} \right) +\frac{7}{4}\zeta(3) -\frac{\pi^3}{32} - \frac{\pi^2}{16}\log 8, $$ Subtracting the second from the first gives us $$ \frac{\pi^3}{16}-G\log 8. $$