$\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.
We make use of the identity
$$ \sum_{n=-\infty}^{\infty} \frac{1}{a^{2} - (x + n\pi)^{2}} = \frac{\cot(x+a) - \cot(x-a)}{2a}, \quad a > 0 \text{ and } x \in \Bbb{R}. $$
Then for $\alpha, \beta > 0$ it follows that
\begin{align*}
I := \mathrm{PV}\int_{0}^{\infty} \frac{\log\cos^{2}(\alpha x)}{\beta^{2} - x^{2}}
&= \frac{1}{2} \mathrm{PV} \int_{-\infty}^{\infty} \frac{\log\cos^{2}(\alpha x)}{\beta^{2} - x^{2}} \, dx \\
&= \frac{\alpha}{2} \mathrm{PV} \int_{-\infty}^{\infty} \frac{\log\cos^{2}x}{(\alpha\beta)^{2} - x^{2}} \, dx \\
&= \frac{\alpha}{2} \sum_{n=-\infty}^{\infty} \mathrm{PV} \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} \frac{\log\cos^{2}x}{(\alpha\beta)^{2} - (x+n\pi)^{2}} \, dx \\
&= \frac{\alpha}{2} \mathrm{PV} \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} \left( \sum_{n=-\infty}^{\infty} \frac{1}{(\alpha\beta)^{2} - (x+n\pi)^{2}} \right) \log\cos^{2}x \, dx \\
&= \frac{1}{4\beta} \mathrm{PV} \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} (\cot(x+\alpha\beta) - \cot(x-\alpha\beta)) \log\cos^{2}x \, dx,
\end{align*}
where interchanging the order of integration and summation is justified by Tonelli's theorem applied to the summation over large indices $n$. Then
\begin{align*}
I
&= \frac{1}{4\beta} \mathrm{PV} \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} (\cot(x+\alpha\beta) - \cot(x-\alpha\beta)) \log\cos^{2}x \, dx \\
&= \frac{1}{2\beta} \mathrm{PV} \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} (\cot(x+\alpha\beta) - \cot(x-\alpha\beta)) \log\left|2\cos x\right| \, dx \tag{1}
\end{align*}
Here, we exploited the following identity to derive (1).
$$ \mathrm{PV} \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} \cot(x+a) \, dx = 0 \quad \forall a \in \Bbb{R}. $$
Now with the substitution $z = e^{2ix}$ and $\omega = e^{2i\alpha\beta}$, it follows that
\begin{align*}
I
&= \frac{1}{2\beta} \Re \mathrm{PV} \int_{|z|=1} \left( \frac{\bar{\omega}}{z - \bar{\omega}} - \frac{\omega}{z - \omega} \right) \log(1 + z) \, \frac{dz}{z}. \tag{2}
\end{align*}
Now consider the following unit circular contour $C$ with two $\epsilon$-indents $\gamma_{\omega,\epsilon}$ and $\gamma_{\bar{\omega},\epsilon}$.
Then the integrand of (2)
$$ f(z) = \left( \frac{\bar{\omega}}{z - \bar{\omega}} - \frac{\omega}{z - \omega} \right) \frac{\log(1 + z)}{z} $$
is holomorphic inside $C$ (since the only possible singularity at $z = 0$ is removable) and has only logarithmic singularity at $z = -1$. So we have
$$ \oint_{C} f(z) \, dz = 0. $$
This shows that
\begin{align*}
I
&= \frac{1}{2\beta} \Re \lim_{\epsilon \downarrow 0} \left( \int_{-\gamma_{\omega,\epsilon}} f(z) \, dz + \int_{-\gamma_{\bar{\omega},\epsilon}} f(z) \, dz \right) \\
&= \frac{1}{2\beta} \Re \left( \pi i \mathrm{Res}_{z=\omega} f(z) + \pi i \mathrm{Res}_{z=\bar{\omega}} f(z) \right) \\
&= \frac{1}{2\beta} \Re \left( - \pi i \log(1 + \omega) + \pi i \log(1 + \bar{\omega}) \right) \\
&= \frac{\pi}{\beta} \arg(1 + \omega)
= \frac{\pi}{\beta} \arctan(\tan (\alpha \beta)).
\end{align*}
In particular, if $\alpha\beta < \frac{\pi}{2}$ then we have
$$ I = \pi \alpha. $$
But due to the periodicity of $\arg$ function, this function draws a scaled saw-tooth function for $\alpha > 0$. Of course, $I$ is an even function of both $\alpha$ and $\beta$, so the final result is obtained by even extension of this saw-tooth function.
Best Answer
The similar integral
$$ \int_0^1\frac{\log^2(1-x)}x\mathrm dx=2\zeta(3) $$
is evaluated in this blog post using the substitution $u=-\log(1-x)$:
$$ \begin{align} \int_0^1\frac{\log^2(1-x)}x\mathrm dx &= \int_0^\infty\frac{u^2}{1-\mathrm e^{-u}}\mathrm e^{-u}\,\mathrm du \\ &= \int_0^\infty u^2\sum_{n=1}^\infty\mathrm e^{-nu}\mathrm du \\ &= \sum_{n=1}^\infty\int_0^\infty u^2\mathrm e^{-nu}\mathrm du \\ &= \sum_{n=1}^\infty\frac2{n^3} \\ &= 2\zeta(3)\;. \end{align} $$
Analogously substituting $u=\log(1+x)$ in the present integral leads to an integral up to $\log2$ that can be evaluated in terms of polylogarithms evaluated at $\frac12$:
$$ \begin{align} &\int_0^{\log2}\frac{\log^2(1+x)}x\mathrm dx \\ =& \int_0^{\log2}\frac{u^2}{\mathrm e^u-1}\mathrm e^u\,\mathrm du \\ =& \int_0^{\log2}\frac{u^2}{1-\mathrm e^{-u}}\mathrm du \\ =& \int_0^{\log2} u^2\sum_{n=0}^\infty\mathrm e^{-nu}\mathrm du \\ =& \sum_{n=0}^\infty\int_0^{\log2} u^2\mathrm e^{-nu}\mathrm du \\ =& \sum_{n=0}^\infty\int_0^{\log2} u^2\mathrm e^{-nu}\mathrm du \\ =& \frac13\log^32+\sum_{n=1}^\infty\frac1n\left(-2^{-n}\log^22+2\int_0^{\log2} u\mathrm e^{-nu}\mathrm du\right) \\ =& \frac13\log^32+\sum_{n=1}^\infty\left(-\frac1n2^{-n}\log^22+\frac2{n^2}\left(-2^{-n}\log2+\int_0^{\log2}\mathrm e^{-nu}\mathrm du\right)\right) \\ =& \frac13\log^32+\sum_{n=1}^\infty\left(-\frac1n2^{-n}\log^22-\frac2{n^2}2^{-n}\log2-\frac2{n^3}\left(2^{-n}-1\right)\right) \\ =& \def\Li{\operatorname{Li}} \frac13\log^32-\Li_1\left(\frac12\right)\log^22-2\Li_2\left(\frac12\right)\log2-2\Li_3\left(\frac12\right)+2\zeta(3) \\ =& \frac13\log^32-\log2\log^22-2\left(\frac{\pi^2}{12}-\frac12\log^22\right)\log2-2\left(\frac16\log^32-\frac{\pi^2}{12}\log2+\frac78\zeta(3)\right)+2\zeta(3) \\ =& \frac{\zeta(3)}4\;. \end{align} $$
Not only is this a rather complicated derivation of a much simpler result; it also looks as if the polylogarithm values may have been obtained using the present integral in the first place.