A First (Magical) Solution by Cornel Ioan Valean
As I also mentioned in this post, the present series, which plays a key part there, is a beast in every way, and finding an elegant approach is an extremely difficult task, often looking like an impossible task. Experienced users in the realm of harmonic series will immediately figure out this monster won't fall easily to its knees. However, it will fall down during the work below!
Let's start with a special particular case of one of the Cornel's generalizations presented during another solution at this address
$$\int_0^1 \frac{x\operatorname{arctanh}^n(x)}{1+a^2 x^2}\textrm{d}x \bigg|_{n=1}$$
$$=\frac{n!}{a^2 2^n}\left((1-2^{-n})\zeta(n+1)+ \Re\biggr\{\operatorname{Li}_{n+1}\left(\frac{a^2-1}{a^2+1}+i\frac{2a}{1+a^2}\right)\biggr\}\right)\bigg|_{n=1}$$
$$=\frac{1}{2}\frac{\arctan^2(a)}{a^2}, \ a \in \mathbb{R}\setminus\{0\}. \tag1$$
Proving the second equality: We want to exploit the Euler's formula, and then write
\begin{equation*}
\begin{aligned}
\int_0^1 \frac{x\operatorname{arctanh}(x)}{1+a^2 x^2}\textrm{d}x\overset{\text{use} \ a=\tan(\theta/2)}{=}&\frac{1}{2\tan^2(\theta/2) }\left(\frac{\pi^2}{12}+ \Re\bigr\{\operatorname{Li}_2\left(-\cos(\theta)+i\sin(\theta)\right)\bigr\}\right)\\
&=\frac{1}{2\tan^2(\theta/2) }\left(\frac{\pi^2}{12}+ \Re\biggr\{\operatorname{Li}_2\left(-e^{-i \theta}\right)\biggr\}\right)\\
&=\frac{1}{2\tan^2(\theta/2) }\left(\frac{\pi^2}{12}- \sum_{n=1}^{\infty} (-1)^{n-1}\frac{\Re\{e^{-i n\theta}\}}{n^2}\right)\\
&=\frac{1}{2\tan^2(\theta/2) }\left(\frac{\pi^2}{12}- \sum_{n=1}^{\infty} (-1)^{n-1}\frac{\cos(n\theta)}{n^2}\right)\\
&=\frac{\theta^2}{8\tan^2(\theta/2)}\overset{\text{use} \ \tan(\theta/2)=a}{=}\frac{1}{2}\frac{\arctan^2(a)}{a^2},
\end{aligned}
\end{equation*}
where during the calculations I used the well-known Fourier series, $\displaystyle \sum_{n=1}^{\infty}(-1)^{n-1}\frac{\cos(nx)}{n^2}=\frac{\pi^2}{12}-\frac{x^2}{4}$,$-\pi\le x \le\pi$, and the calculations are complete.
Another solution: In this solution I'll only consider the restriction of the type $|x|\le1, \ x\neq 0$,
\begin{equation*}
\begin{aligned}
\int_0^1 \frac{x\operatorname{arctanh}(x)}{1+a^2 x^2}\textrm{d}x&=\int_0^1 \sum_{n=1}^{\infty}(-1)^{n-1} x(a x)^{2n-2}\operatorname{arctanh}(x)\textrm{d}x\\
&=\sum_{n=1}^{\infty}(-1)^{n-1}a^{2n-2} \int_0^1 x^{2n-1}\operatorname{arctanh}(x)\textrm{d}x\\
&=\frac{1}{4}\sum_{n=1}^{\infty}(-1)^{n-1}a^{2n-2} \frac{2H_{2n}-H_n}{n}\\
&=\frac{1}{2}\frac{\arctan^2(a)}{a^2},
\end{aligned}
\end{equation*}
where in the calculations I used the logarithmic integrals $\displaystyle \int_0^1 x^{2n-1}\log(1-x)\textrm{d}x=-\frac{H_{2n}}{2n}$, then $\displaystyle \int_0^1 x^{2n-1}\log(1+x)\textrm{d}x=\frac{H_{2n}-H_n}{2n}$, (A simple proof may be found in this RG article), and finally the classical Cauchy product, $\displaystyle \arctan^2(x)=\frac{1}{2}\sum_{n=1}^{\infty}(-1)^{n-1}x^{2n}\frac{2H_{2n}-H_n}{n}, \ |x|\le1$. Details on how to start from this variant and extend it to $a \in \mathbb{R}\setminus\{0\}$ will be found in the sequel of (Almost) Impossible Integrals, Sums, and Series.
Yet, one more solution: This way is also going to be extremely elegant, and multiplying the auxiliary integral by $2a^2$ and cleverly integrating by parts, we get
\begin{equation*}
\begin{aligned}
\int_0^1 2a^2\frac{x\operatorname{arctanh}(x)}{1+a^2 x^2}\textrm{d}x&=-\int_0^1 \log\left(\frac{1+a^2x^2}{1+a^2}\right)\frac{1}{1-x^2}\textrm{d}x\\
&=-\int_0^1 \int_0^a \frac{\partial}{\partial t}\biggr\{\log\left(\frac{1+t^2x^2}{1+t^2}\right)\biggr\}\frac{1}{1-x^2}\textrm{d}t \textrm{d}x\\
&=\int_0^1 \int_0^a\frac{2 t}{(1+t^2)(1+t^2 x^2)}\textrm{d}t \textrm{d}x\\
&=\int_0^a\int_0^1 \frac{2 t}{(1+t^2)(1+t^2 x^2)}\textrm{d}x \textrm{d}t\\
&=2\int_0^a \frac{\arctan(t)}{1+t^2}\textrm{d}t\\
&=\arctan^2(a),
\end{aligned}
\end{equation*}
which brings an end to the present solution.
The main calculations: We may start from the function $\displaystyle f(x)=\frac{(\pi/2-\arctan(t))^2\operatorname{arctanh}(t)}{t}$ which we integrate from $0$ to $\infty$ (or I could start directly from the blue integral), and then using the simple fact that based on $(1)$ we have $\displaystyle \frac{(\pi/2-\arctan(t))^2}{t}=2\int_0^1 \frac{x\operatorname{arctanh}(x)}{t(t^2+x^2)}\textrm{d}x$, we write that
$$\Re\biggr \{\int_0^{\infty} \frac{(\pi/2-\arctan(t))^2\operatorname{arctanh}(t)}{t}\textrm{d}t\biggr\}=\Re\biggr \{2\int_0^{\infty}\operatorname{arctanh}(t)\left( \int_0^1 \frac{x\operatorname{arctanh}(x)}{t(t^2+x^2)}\textrm{d}x\right)\textrm{d}t\biggr\}$$
$$=\Re\biggr \{2\int_0^1x\operatorname{arctanh}(x)\left(\int_0^{\infty} \frac{\operatorname{arctanh}(t)}{t(x^2+t^2)}\textrm{d}t\right)\textrm{d}x\biggr\}$$
$$=2\int_0^1x\operatorname{arctanh}(x)\left(\int_0^{\infty}\left(\operatorname{PV} \int_0^1 \frac{1}{(x^2+t^2)(1-t^2u^2)}\textrm{d}u\right)\textrm{d}t\right)\textrm{d}x$$
$$=2\int_0^1x\operatorname{arctanh}(x)\left(\int_0^1\left(\operatorname{PV} \int_0^{\infty} \frac{1}{(x^2+t^2)(1-u^2t^2)}\textrm{d}t\right)\textrm{d}u\right)\textrm{d}x$$
$$=\pi\int_0^1\operatorname{arctanh}(x)\left(\int_0^1\frac{1}{1+x^2 u^2}\textrm{d}u\right)\textrm{d}x$$
$$=\pi\int_0^1\frac{\arctan(x)\operatorname{arctanh}(x)}{x}\textrm{d}x.\tag2$$
Splitting the first integral in $(2)$ at $t=1$, using trivial identities with $\arctan(x)$ and $\operatorname{arctanh}(x)$, employing elementary variable changes to reduce the integration interval to the space from $t=0$ to $t=1$, we write
$$\Re\biggr \{\int_0^{\infty} \frac{(\pi/2-\arctan(t))^2\operatorname{arctanh}(t)}{t}\textrm{d}t\biggr\}$$
$$\int_0^1 \frac{(\pi/2-\arctan(t))^2\operatorname{arctanh}(t)}{t}\textrm{d}t+\Re\biggr \{\int_1^{\infty} \frac{(\pi/2-\arctan(t))^2\operatorname{arctanh}(t)}{t}\textrm{d}t\biggr\}$$
$$=\frac{\pi^2}{4}\int_0^1\frac{\operatorname{arctanh}(t)}{t}\textrm{d}t-\pi\int_0^1\frac{\arctan(t)\operatorname{arctanh}(t)}{t}\textrm{d}t+2\color{blue}{\int_0^1\frac{\arctan^2(t)\operatorname{arctanh}(t)}{t}\textrm{d}t}$$
$$=\frac{45}{16}\zeta(4)-\pi\int_0^1\frac{\arctan(t)\operatorname{arctanh}(t)}{t}\textrm{d}t+2\color{blue}{\int_0^1\frac{\arctan^2(t)\operatorname{arctanh}(t)}{t}\textrm{d}t}\tag3.$$
Upon combining $(2)$ and $(3)$ we get a first useful and powerful relation which for the art's sake I'll put as a representation of $\zeta(4)$,
$$\zeta(4)$$
$$=\frac{32}{45} \pi \int_0^1 \frac{\arctan(t)\operatorname{arctanh(t)}}{t}\textrm{d}t-\frac{32}{45} \color{blue}{\int_0^1 \frac{\arctan^2(t)\operatorname{arctanh}(t)}{t} \textrm{d}t}.\tag4$$
The first integral in $(4)$ is a known one, immediately reducible to the work with classical results, and we may put it in the form
$$\int_0^1 \frac{\arctan(t)\operatorname{arctanh(t)}}{t}\textrm{d}t=\frac{\pi ^3}{32}+\frac{1}{2} \int_0^1 \frac{\arctan(t)\log(t)}{1+t}\textrm{d}t+\frac{1}{2} \Im\biggr\{\int_0^1 \frac{\log (1+i t) \log (t)}{1-t} \textrm{d}t\biggr\}$$
$$=\frac{1}{2}\log(2)G+\frac{\pi}{16}\log^2(2)+\frac{3}{64}\pi^3-\Im\{\operatorname{Li}_3(1+i)\}.\tag5$$
The first integral in $(5)$ appears in this post with more solutions, and it also appears in (Almost) Impossible Integrals, Sums, and Series (see page $14$, Sect. $1.24$). The second integral is very elegantly derived if we consider $\displaystyle I(a)=\int_0^1 \frac{\log (1+a t) \log (t)}{1-t} \textrm{d}t$ and differentiate under the integral sign with respect to $a$.
Using the logarithmic integrals and the Cauchy product employed during the calculation of the second solution to the auxiliary result in $(1)$, we get the following integral transformation to series
$$\color{blue}{\int_0^1 \frac{\arctan^2(t)\operatorname{arctanh}(t)}{t} \textrm{d}t}$$
$$=\frac{1}{8}\sum _{n=1}^{\infty}(-1)^{n-1} \frac{ H_n^2}{n^2}+\frac{1}{2}\sum_{n=1}^{\infty} (-1)^{n-1}\frac{H_{2 n}^2}{n^2}-\frac{1}{2}\color{red}{\sum _{n=1}^{\infty} (-1)^{n-1}\frac{ H_n H_{2 n}}{n^2}}.\tag6$$
Combining $(4)$, $(5)$, and $(6)$ and observing the first two series are known, we arrive at the desired value of the series,
$$\color{red}{\sum _{n=1}^{\infty} (-1)^{n-1}\frac{ H_n H_{2 n}}{n^2}}$$
$$\color{red}{=2 G^2-2\log(2)\pi G-\frac{1}{8}\log^4(2)-\frac{21}{8}\log(2)\zeta(3)+\frac{1}{4}\log^2(2)\pi ^2+\frac{773}{5760}\pi ^4}$$
$$\color{red}{-4 \pi \Im\biggr\{\operatorname{Li}_3\left(\frac{1+i}{2}\right)\biggr\}-3 \operatorname{Li}_4\left(\frac{1}{2}\right)},$$
where I also used polylogarithmic relations as the ones found here.
Note that for the first two resulting harmonic series in $(6)$, we have
$$\sum_{n=1}^{\infty} (-1)^{n-1} \frac{H_n^2}{n^2}$$
$$=\frac{41}{16}\zeta (4)-\frac{7}{4}\log (2) \zeta (3)+\frac{1}{2}\log ^2(2)\zeta (2)-\frac{1}{12} \log ^4(2) -2 \text{Li}_4\left(\frac{1}{2}\right),$$
which is a classical result, and it may also be found in (Almost) Impossible Integrals, Sums, and Series, page $310$,
and
$$\sum_{n=1}^{\infty} (-1)^{n-1} \frac{H_{2n}^2}{n^2}$$
$$=2 G^2-\log(2)\pi G+\frac{231}{32}\zeta(4)-\frac{35}{16}\log(2)\zeta(3)+\log^2(2)\zeta(2)-\frac{5}{48}\log^4(2)$$
$$-2 \pi \Im\biggr\{\operatorname{Li}_3\left(\frac{1+i}{2}\right)\biggr\}-\frac{5}{2}\operatorname{Li}_4\left(\frac{1}{2}\right),$$
which is given and easily derived in this post.
End of story
More series from the atypical alternating weight 4 class:
$$\sum_{n=1}^{\infty} (-1)^{n-1} \frac{H_{2n}^3}{n}$$
$$=2G^2-\frac{3}{4}\log(2)\pi G+\frac{1055}{256}\zeta(4)-\frac{93}{64}\log(2)\zeta(3)+\frac{21}{32}\log^2(2)\zeta(2)-\frac{1}{32}\log^4(2)$$
$$-\frac{3}{2}\pi \Im\biggr\{\operatorname{Li}_3\left(\frac{1+i}{2}\right)\biggr\};$$
$$\sum_{n=1}^{\infty} (-1)^{n-1} \frac{H_{2n}H_{2n}^{(2)}}{n}$$
$$
=\frac{1}{4}\log(2)\pi G-\frac{137}{128}\zeta(4)+\frac{35}{64}\log(2)\zeta(3)-\frac{3}{8}\log^2(2)\zeta(2)+\frac{5}{96}\log^4(2)$$
$$+\frac{\pi}{2}\Im\biggr\{\operatorname{Li}_3\left(\frac{1+i}{2}\right)\biggr\}+\frac{5}{4}\operatorname{Li}_4\left(\frac{1}{2}\right).$$
They are straightforward by exploiting their related generating functions that are immediately built based on the generating functions in (Almost) Impossible Integrals, Sums, and Series (see page $284$, Sect. $4.10$). More details and results are included in the sequel of the mentioned book.
And more series ... from the same class:
$$\sum_{n=1}^{\infty} (-1)^{n-1}\frac{H_{2n}^{(3)}}{n}=\frac{199}{128}\zeta (4)-\frac{3}{32} \log (2)\zeta (3)-G^2,$$
where this series is part of the generalization found in Cornel's answer here. It may also be found in the article A simple strategy of calculating two alternating harmonic series generalizations by C.I. Valean at this link.
Best Answer
A slightly different approach where I will make use of a particular Euler sum.
Let $$I = \int_0^1 \frac{\ln (1 - x) \ln (1 - x^3)}{x} \, dx.$$ Expanding the $\ln (1 - x^3)$ term gives \begin{align} I &= - \sum_{n = 1}^\infty \frac{1}{n} \int_0^1 x^{3n - 1} \ln (1 - x) \, dx\tag1 \end{align} Making use of the result (for a proof of this, see here) $$\int_0^1 x^{n - 1} \ln (1 - x) \, dx = -\frac{H_n}{n}.$$ Re-indexing, namely $n \mapsto 3n$ gives $$\int_0^1 x^{3n - 1} \ln (1 - x) \, dx = -\frac{H_{3n}}{3n}.$$ Substitution of this result into (1) reduces our integral $I$ to the following Euler sum $$I = \frac{1}{3} \sum_{n = 1}^\infty \frac{H_{3n}}{n^2} = 3 \sum_{n = 1}^\infty \frac{H_{3n}}{(3n)^2}.$$
For the Euler sum, since the series converges absolutely, terms in the sum may be rearranged. Doing so we have \begin{align} \sum_{n =1}^\infty \frac{H_{3n}}{(3n)^2} &= \frac{H_3}{3^2} + \frac{H_6}{6^2} + \frac{H_9}{9^2} + \cdots\\ &= \frac{2}{3} \left [\frac{3}{2} \frac{H_3}{3^2} + \frac{3}{2} \frac{H_6}{6^2} + \frac{3}{2} \frac{H_9}{9^2} + \cdots \right ]\\ &= \frac{2}{3} \left [\left (\frac{H_3}{3^2} + \frac{H_6}{6^2} + \frac{H_9}{9^2} + \cdots \right ) + \frac{1}{2} \left (\frac{H_3}{3^2} + \frac{H_6}{6^2} + \frac{H_9}{9^2} + \cdots \right ) \right ]\\ &= \frac{2}{3} \left [\left (-\frac{1}{2} \frac{H_1}{1^2} - \frac{1}{2} \frac{H_2}{2^2} + \frac{H_3}{3^2} - \frac{1}{2} \frac{H_4}{4^2} - \frac{1}{2} \frac{H_5}{5^2} + \frac{H_6}{6^2} - \cdots \right )\right.\\ & \qquad + \left. \frac{1}{2} \left (\frac{H_1}{1^2} + \frac{H_2}{2^2} + \frac{H_3}{3^2} + \frac{H_4}{4^2} + \cdots \right ) \right ]\\ &= \frac{2}{3} \sum_{n = 1}^\infty \frac{H_n}{n^2} \cos \left (\frac{2 \pi n}{3} \right ) + \frac{1}{3} \sum_{n = 1}^\infty \frac{H_n}{n^2} \tag2\\ &= \frac{2}{3} \zeta (3) + \frac{2}{3} \operatorname{Re} \sum_{n = 1}^\infty \frac{H_n}{n^2} \left (e^{\frac{2 \pi i}{3}} \right )^n \end{align} Note in (2) the well-known result of $\sum_{n = 1}^\infty \frac{H_n}{n^2} = 2 \zeta (3)$ has been used.
The sum can now be found by making use of the following generating function (for a simple proof of this result, see here) $$\sum_{n=1}^\infty\frac{H_n}{n^2}x^n=\operatorname{Li}_3(x)-\operatorname{Li}_3(1-x)+\ln(1-x)\operatorname{Li}_2(1-x)+\frac12\ln x\ln^2(1-x)+\zeta(3).$$ Setting $x = e^{\frac{2 \pi i}{3}}$ gives \begin{align} \sum_{n = 1}^\infty \frac{H_{3n}}{(3n)^2} &= \frac{2}{3} \zeta (3) + \frac{2}{3} \operatorname{Re} \left [\operatorname{Li}_3 (e^{\frac{2 \pi i}{3}}) - \operatorname{Li}_3 (1 - e^{\frac{2 \pi i}{3}}) + \ln (1 - e^{\frac{2 \pi i}{3}}) \operatorname{Li}_2 (1 - e^{\frac{2 \pi i}{3}}) \right.\\ & \qquad \left. + \frac{1}{2} \ln (e^{\frac{2 \pi i}{3}}) \ln^2 (1 - e^{\frac{2 \pi i}{3}}) + \zeta (3) \right ]\\ &= \frac{5}{3} \zeta (3) + \frac{2}{3} \operatorname{Re} \left [\operatorname{Li}_3 (e^{\frac{2 \pi i}{3}}) - \operatorname{Li}_3 (1 - e^{\frac{2 \pi i}{3}}) + \ln (1 - e^{\frac{2 \pi i}{3}}) \operatorname{Li}_2 (1 - e^{\frac{2 \pi i}{3}}) \right.\\ & \qquad \left. + \frac{\pi i}{3} \ln^2 (1 - e^{\frac{2 \pi i}{3}}) \right ]\tag3 \end{align} Now, since (this part is tedious, but readily doable) \begin{align} \operatorname{Re} \operatorname{Li}_3 (e^{\frac{2 \pi i}{3}}) &= -\frac{4}{9} \zeta (3)\\ \operatorname{Re} \operatorname{Li}_3 (1 - e^{\frac{2 \pi i}{3}}) &= \frac{\pi^2}{18} \ln 3 + \frac{13}{18} \zeta (3)\\ \operatorname{Re} \left [i \ln^2 (1 - e^{\frac{2 \pi i}{3}}) \right ] &= \frac{\pi}{6} \ln 3\\ \operatorname{Re} \left [\ln (1 - e^{\frac{2 \pi i}{3}}) \operatorname{Li}_2 (1 - e^{\frac{2 \pi i}{3}}) \right ] &= \frac{\pi^3}{27 \sqrt{3}} - \frac{\pi}{18 \sqrt{3}} \psi^{(1)} \left (\frac{1}{3} \right ) \end{align} The Euler sum in (3) thus becomes $$\sum_{n = 1}^\infty \frac{H_{3n}}{(3n)^2} = \frac{5}{9} \zeta (3) + \frac{2 \pi^3}{81 \sqrt{3}} - \frac{\pi}{27 \sqrt{3}} \psi^{(1)} \left (\frac{1}{3} \right ),$$ so that we finally arrive at the following value for the integral (and hence your double sum) of