One way to get the closed form of the following arctanh integral,
$$\int_0^1 \frac{\operatorname{arctanh}^2\left(x^2\right) \log^2(x)}{1 +x^2} \mathrm dx$$ $$=\log (2)\frac{\pi ^2}{4}G-\frac{\pi }{2}G^2 +\frac{7\pi}{16} \log(2)\zeta(3)-\log^2(2)\frac{\pi ^3}{64} - \frac{17 \pi ^5}{768}+ \pi ^2 \Im\left\{\text{Li}_3\left(\frac{1+i}{2}\right)\right\},$$ is by exploiting the following result (I'll put the general form, although not needed),
$$\int_0^{\pi/4} \cos(2 n x) (\operatorname{arctanh}(\tan(x)))^2\mathrm dx=\frac{1}{4}\int_0^{\pi/4} \cos(2 n x) (\operatorname{arctanh}(\sin(2x)))^2\mathrm dx$$ $$=\frac{1}{4}\int_0^{\pi/4} \cos(2 n x) \log^2\left(\tan\left(\frac{\pi}{4}+x\right)\right)\mathrm dx=\frac{1}{4}\int_0^{\pi/4} \cos(2 n x) \log^2\left(\tan\left(\frac{\pi}{4}-x\right)\right)\mathrm dx$$ $$\small =\frac{1}{4}\int_0^{\pi/4} \cos(2 n x) \log^2\left(\sec(2x)+\tan(2x)\right)\mathrm dx=\frac{1}{4}\int_0^{\pi/4} \cos(2 n x) \log^2\left(\sec(2x)-\tan(2x)\right)\mathrm dx$$ \begin{equation*}= \begin{cases} \displaystyle \quad \frac{\pi^2}{96}\frac{(-1)^{m-1}}{2m-1}-\frac{1}{8}(-1)^{m-1} \frac{H_{m-1}^2}{2m-1}-\frac{1}{8}(-1)^{m-1}\frac{H_{m-1}^{(2)}}{2m-1} \\[10pt] \displaystyle \quad +\frac{1}{2}\frac{(-1)^{m-1}}{2m-1}\sum_{k=1}^{m-1}\frac{ H_{2k}}{k}, \qquad \qquad \quad n=2m-1, \ m \in \mathbb{N}; \\[10pt] \displaystyle \quad =\frac{\pi}{8}(-1)^m \frac{H_{2m}}{m}-\frac{\pi}{16}(-1)^m \frac{H_m}{m}, \ \ n=2m, \ m \in \mathbb{N}, \end{cases} \end{equation*} where $H_m^{(p)}=1+\frac{1}{2^p}+\cdots+\frac{1}{m^p}, \ p\ge1$, is the $m$th generalized harmonic number of order $p$, found in Another Valuable Elementary Generalized Trigonometric Integral with the Square of the Inverse Hyperbolic Tangent by C.I. Valean.
Also, we need the Fourier series of $\log^2(\sin(x))$ and the harmonic series,
$$\sum _{n=1}^{\infty} (-1)^{n-1} \frac{H_n H_{2n}}{n^2}$$
$$=2G^2-2\log(2)\pi G-\frac{1}{8}\log^4(2)-\frac{21}{8}\log(2)\zeta(3)+\frac{3}{2}\log^2(2)\zeta(2)+\frac{773}{64}\zeta(4)$$
$$-4\pi \Im\biggr \{\operatorname{Li}_3\left(\frac{1+i}{2}\right)\biggr \}-3\operatorname{Li}_4\left(\frac{1}{2}\right),$$
found in More (Almost) Impossible Integrals, Sums, and Series: A New Collection of Fiendish Problems and Surprising Solutions (2023), which is the sequel of (Almost) Impossible Integrals, Sums, and Series (2019).
Question: I'm not necessarily interested in full solutions (partial solutions are fine - e.g., showing the transformation to more convenient sums of key harmonic series), but in different views of approaching such integrals (in large steps).
The result in the mentioned paper seems to be of great utility, as seen in this other post Calculating $\int_0^1 \frac{\arctan^3(x) \log(x)} x\,\mathrm dx$.
A note: to get an idea of what may happen when using harmonic series (this is just one scenario), we may recall that $$\frac{\operatorname{arctanh}^2(x)}{1-x}$$ $$=\frac{1}{2}\sum_{n=1}^{\infty} x^n \left(\frac{1}{2}H_n^2+H_n \overline{H}_n-\frac{1}{2}\overline{H}_n^2+\overline{H}_n^{(2)}-2\sum_{k=1}^n (-1)^{k-1}\frac{H_k}{k}\right),$$ where $H_n^{(m)}=1+\frac{1}{2^m}+\cdots+\frac{1}{n^m}, \ m\ge1$, is the $n$th generalized harmonic number of order $m$ and $\overline{H}_n^{(m)}=1-\frac{1}{2^m}+\cdots+(-1)^{n-1}\frac{1}{n^m}, \ m\ge1$, denotes the $n$th generalized skew-harmonic number of order $m$, which is presented in Analogues of the established Landen-type identities in the form of series and some related Cauchy products by C.I. Valean.
BIG BONUS: Combining the initial integral and $\displaystyle \int_0^1 \frac{\log^2(1-x^4) \log^2(x)(1-x^2)}{1 -x^4} \mathrm dx$, which has a Beta function structure (so, easy to calculate), we immediately arrive at $$\int_0^1 \frac{\log(1-x^2)\log^2(x)\log(1+x^2)}{1 +x^2} \mathrm dx$$ $$=\pi G^2+\frac{1}{2}\log (2)\pi^2G-\frac{21}{2}\zeta (3)G+\frac{7}{8} \log (2)\pi \zeta (3) +\frac{5}{32} \log ^2(2)\pi ^3 +\frac{3}{16} \log (2)\pi ^4 +\frac{\pi ^5}{16}$$ $$-\frac{3}{128} \log (2) \psi ^{(3)}\left(\frac{1}{4}\right)-\pi ^2 \Im\left\{\operatorname{Li}_3\left(\frac{1+i}{2}\right)\right\}.$$
The simpler version, with $x$ instead of $x^2$, $$\int_0^1\frac{\log(1-x)\log^2(x)\log(1+x)}{1+x}\textrm{d}x$$ $$=\frac{213}{16}\zeta(5)-\frac{5}{2}\zeta(2)\zeta(3)-\frac{49}{8}\log(2)\zeta(4)+\frac{7}{4}\log^2(2)\zeta(3)+\frac{1}{3}\log^3(2)\zeta(2)-\frac{1}{10}\log^5(2)$$ $$-4\log(2)\operatorname{Li}_4\left(\frac{1}{2}\right)-8\operatorname{Li}_5\left(\frac{1}{2}\right),$$ is already given in the book above, page $56$.
More content: We also have the similar version (this time the argument of the arctanh function is $x$ instead of $x^2$) $$\int_0^1 \frac{\operatorname{arctanh}^2(x) \log^2(x)}{1 +x^2} \mathrm dx$$ $$= \frac{7 }{256}\pi ^5+\log ^2(2)\frac{\pi ^3}{32}-\log (2)\frac{\pi ^2}{2} G-\pi ^2 \Im\left\{\text{Li}_3\left(\frac{1+i}{2}\right)\right\}.$$
A first solution can be obtained by calculating and exploiting the version $\displaystyle \int_0^{\infty} \frac{\operatorname{arctanh}^2(1/x) \log^2(x)}{1 +x^2} \mathrm dx $, where it is very useful to consider an extension (in the sense of Cauchy principal value) of the auxiliary result, $\displaystyle \int_0^1 \frac{x\operatorname{arctanh}(x)}{1-a^2 x^2}\textrm{d}x=\frac{1}{2}\frac{\operatorname{arctanh}^2(a)}{a^2}, \ |a|<1$, found in the sequel, page $9$, which is a strategy that can also be exploited to get a second solution to the main integral, and for a second solution we can exploit again the auxiliary result in the mentioned paper.
