Proving that $\int_0^1 \frac{\log^2(x)\tanh^{-1}(x)}{1+x^2}dx=\beta(4)-\frac{\pi^2}{12}G$
I was able to solve this problem on my own.
Using integration by parts, $$\begin{align*} &\; \int_0^1 \frac{\log^2(x)\tanh^{-1}(x)}{1+x^2}dx \\ &= -2\int_0^1 \frac{\log(x)\tan^{-1}(x)\tanh^{-1}(x)}{x}dx-\int_0^1 \frac{\log^2(x)\tan^{-1}(x)}{1-x^2}dx \tag{1} \end{align*}$$
I posted the solution to both these integrals on another forum. Here are the links:
http://integralsandseries.prophpbb.com/topic711.html#p3975
http://integralsandseries.prophpbb.com/topic245.html#p1680
$$\begin{align*}\int_0^1\frac{\log(x)\tan^{-1}(x)\tanh^{-1}(x)}{x}dx &= \frac{\pi^2}{16}G-\frac{7\pi\zeta(3)}{32} \tag{2}\\ \int_0^1\frac{\log^2(x)\tan^{-1}(x)}{1-x^2}dx &= -\beta(4)-\frac{\pi^2}{24}G+\frac{7\pi}{16}\zeta(3)\tag{3} \end{align*}$$ $G$ denotes the Catalan's constant and $\beta(4)=\sum_{n=0}^\infty\frac{(-1)^n}{(2n+1)^4}$. Substituting these two results in equation (1) gives: $$\int_0^1 \frac{\log^2(x)\tanh^{-1}(x)}{1+x^2}dx=\beta(4)-\frac{\pi^2}{12}G \tag{4}$$
Proof sketch of integrals (2) and (3) : (Please see the above links for a more detailed answer)
The idea behind evaluating (2) and (3) is breaking them down into Euler Sums. Using the taylor series expansion $\tan^{-1}(x)=\sum_{n=0}^\infty\frac{(-1)^n x^{2n+1}}{2n+1}$ and integrating term-wise, we obtain the following relations:
\begin{align*} \int_0^1\frac{\log(x)\tan^{-1}(x)\tanh^{-1}(x)}{x}dx &= -\log(2)\frac{\pi^3}{32}-\frac{1}{2}\sum_{n=0}^\infty \frac{(-1)^n}{(2n+1)^3}\left( \gamma+\psi_0(n+1)\right) \\ &\;+\frac{1}{4}\sum_{n=0}^\infty \frac{(-1)^n \psi_1(n+1)}{(2n+1)^2} \tag{5}\\ \int_0^1 \frac{\log^2(x)\tan^{-1}(x)}{1-x^2}dx &=-\frac{1}{8}\sum_{n=0}^\infty\frac{(-1)^n\psi_2(n+1)}{2n+1}\tag{6} \end{align*}
These Euler Sums can be evaluated using the techniques shown in the paper "Euler Sums and Contour Integral Representations" by Philippe Flajolet and Bruno Salvy. Here is it's link. \begin{align*} \sum_{n=0}^\infty\frac{(-1)^n\psi_2(n+1)}{2n+1} &= 8\beta(4)+\frac{\pi^2}{3}G-\frac{7\pi}{2}\zeta(3) \\ \sum_{n=0}^\infty\frac{(-1)^n\psi_1(n+1)}{(2n+1)^2} &= 6\beta(4)+\frac{\pi^2}{4}G-\frac{7\pi}{4}\zeta(3) \\ \sum_{n=0}^\infty \frac{(-1)^n\left( \gamma+\psi_0(n+1)\right)}{(2n+1)^3} &= 3\beta(4)-\frac{7\pi}{16}\zeta(3)-\frac{\pi^3}{16}\log(2) \end{align*} Substituting these into equations (5) and (6) gives us the integrals (2) and (3).
A related integral
Using similar techniques, we can show that $$\displaystyle \int_0^1 \frac{\log^2(x)\tan^{-1}(x)}{x\left(1-x^2 \right)}dx=\beta(4)+\frac{7\pi \zeta(3)}{64}+\frac{\pi^3 \log(2)}{32}$$
UPDATE
Here is a beautiful generalization of the integral calculated below.
Let $n$ be a natural number. Then, we have $$\int_0^1 \frac{\log^{2n-1}(x)\operatorname{arctanh}(x)\arctan(x)}{x}\textrm{d}x$$ $$=\frac{\pi}{4}\left(2^{-2 n-1}-1\right) \zeta (2 n+1) (2n-1)!$$ $$+\frac{\pi}{16} \lim_{s\to0}\left(\frac{d^{2n-1}}{ds^{2n-1}}\left(\frac{1}{s}\tan\left(\frac{\pi}{4}s\right)\left(\pi+\psi\left(\frac{3}{4}-\frac{s}{4}\right)-\psi\left(\frac{1}{4}-\frac{s}{4}\right)\right)\right)\right),$$ where $\zeta$ represents the Riemann zeta function and $\psi$ denotes the Digamma function.
All details will appear soon in a new paper.
A solution by Cornel I. Valean to one of the resulting integrals in Shobhit Bhatnagar's post
Let's show that $$\int_0^1 \frac{\log(x)\operatorname{arctanh}(x)\arctan(x)}{x}\textrm{d}x=\frac{3}{8}\zeta(2)G-\frac{7}{32}\pi \zeta(3),$$ without using harmonic series.
We want to begin with the variable change $x \mapsto 1/x$ and use that $\arctan(x)+\arctan\left(\frac{1}{x}\right)=\frac{\pi}{2}, \ x>0$, and $\operatorname{arctanh}\left(\frac{1}{x}\right)-\operatorname{arctanh}(x)=\frac{\pi}{2}i, \ x>1$. Then, we write
$$I=\int_0^1 \frac{\log(x)\operatorname{arctanh}(x)\arctan(x)}{x}\textrm{d}x=-\int_1^{\infty} \frac{\log(x)\operatorname{arctanh}(1/x)\arctan(1/x)}{x}\textrm{d}x$$ $$=-\int_1^{\infty} \frac{\log(x)(\pi/2i+\operatorname{arctanh}(x))(\pi/2-\arctan(x))}{x}\textrm{d}x$$ $$=\Re\biggr \{\int_1^{\infty}\left(\frac{\log(x)\operatorname{arctanh}(x)\arctan(x)}{x}-\frac{\pi}{2}\frac{\operatorname{arctanh}(x)\log(x)}{x}\right)\textrm{d}x \biggr\}$$ $$=\Re\biggr \{\left(\int_0^{\infty}-\int_0^{1}\right)\left(\frac{\log(x)\operatorname{arctanh}(x)\arctan(x)}{x}-\frac{\pi}{2}\frac{\operatorname{arctanh}(x)\log(x)}{x}\right)\textrm{d}x \biggr\}$$ $$=\Re\biggr \{\int_0^{\infty}\left(\frac{\log(x)\operatorname{arctanh}(x)\arctan(x)}{x}-\frac{\pi}{2}\frac{\operatorname{arctanh}(x)\log(x)}{x}\right)\textrm{d}x \biggr\}$$ $$-\underbrace{\int_0^1 \frac{\log(x)\operatorname{arctanh}(x)\arctan(x)}{x}\textrm{d}x}_{\displaystyle I}+\frac{\pi}{2}\underbrace{\int_0^1\frac{\operatorname{arctanh}(x)\log(x)}{x}\textrm{d}x}_{\displaystyle -7/8\zeta(3)},$$ whence we obtain that $$I=\int_0^1 \frac{\log(x)\operatorname{arctanh}(x)\arctan(x)}{x}\textrm{d}x$$ $$=\frac{1}{2}\Re\biggr \{\int_0^{\infty}\left(\frac{\log(x)\operatorname{arctanh}(x)\arctan(x)}{x}-\frac{\pi}{2}\frac{\operatorname{arctanh}(x)\log(x)}{x}\right)\textrm{d}x \biggr\}-\frac{7}{32}\pi \zeta(3)$$ $$=\frac{1}{2}\Re\biggr \{\int_0^{\infty}\frac{\log(x)\operatorname{arctanh}(x)\arctan(x)}{x}\textrm{d}x \biggr\}-\frac{\pi}{4}\Re\biggr \{\int_0^{\infty}\frac{\operatorname{arctanh}(x)\log(x)}{x}\textrm{d}x \biggr\}$$ $$-\frac{7}{32}\pi \zeta(3),\tag1$$
where in the calculations I used the fact that $\int_0^1\frac{\operatorname{arctanh}(x)\log(x)}{x}\textrm{d}x=\int_0^1 \log(x) \sum_{n=1}^{\infty} \frac{x^{2n-2}}{2n-1}\textrm{d}x$$ $$=\sum_{n=1}^{\infty} \frac{1}{2n-1}\int_0^1 x^{2n-2}\log(x)\textrm{d}x=-\sum_{n=1}^{\infty} \frac{1}{(2n-1)^3}=-\frac{7}{8}\zeta(3)$.
For the second integral in $(1)$ , we have $$\Re\biggr \{\int_0^{\infty}\frac{\operatorname{arctanh}(x)\log(x)}{x}\textrm{d}x \biggr\}=\int_0^1\frac{\operatorname{arctanh}(x)\log(x)}{x}\textrm{d}x+\Re\biggr \{\int_1^{\infty}\frac{\operatorname{arctanh}(x)\log(x)}{x}\textrm{d}x \biggr\}$$ $$=\int_0^1\frac{\operatorname{arctanh}(x)\log(x)}{x}\textrm{d}x+\Re\biggr \{\int_1^{\infty}\frac{(\operatorname{arctanh}(1/x)-\pi/2 i)\log(x)}{x}\textrm{d}x \biggr\}$$ $$=\int_0^1\frac{\operatorname{arctanh}(x)\log(x)}{x}\textrm{d}x-\Re\biggr \{\int_0^1\frac{(\operatorname{arctanh}(x)-\pi/2 i)\log(x)}{x}\textrm{d}x \biggr\}$$ $$=0.\tag2$$
Combining the results in $(1)$ and $(2)$, we arrive at $$I=\int_0^1 \frac{\log(x)\operatorname{arctanh}(x)\arctan(x)}{x}\textrm{d}x$$ $$=\frac{1}{2}\Re\biggr \{\int_0^{\infty}\frac{\log(x)\operatorname{arctanh}(x)\arctan(x)}{x}\textrm{d}x \biggr\}-\frac{7}{32}\pi \zeta(3). \tag3$$
At this point, we consider the generalized integral result,
$$\displaystyle J(s)=\Re \biggr\{\int_0^{\infty}x^{s-1}\operatorname{arctanh}(x)\arctan(x)\textrm{d}x\biggr\}$$ $$=\frac{\pi}{8s}\tan\left(\frac{\pi}{4}s\right)\left(\pi+\psi\left(\frac{3}{4}-\frac{s}{4}\right)-\psi\left(\frac{1}{4}-\frac{s}{4}\right)\right),$$ $0>s>-2$ (which can be extended to $1>s>-2$), we want to prove. (Is this new in the mathematical literature?)
Using integral representations of $\arctan(x)$ and $\operatorname{arctanh}(x)$, we write $$J(s)=\Re \biggr\{\int_0^{\infty}x^{s-1}\operatorname{arctanh}(x)\arctan(x)\textrm{d}x\biggr\}$$ $$=\int_0^{\infty}\left(\int_0^1\left( PV\int_0^1\frac{x^{s+1}}{(1-y^2 x^2)(1+z^2 x^2)}\textrm{d}y\right) \textrm{d}z\right)\textrm{d}x$$ $$=\int_0^1\left( \int_0^1\left(PV\int_0^{\infty}\frac{x^{s+1}}{(1-y^2 x^2)(1+z^2 x^2)}\textrm{d}x\right) \textrm{d}z\right)\textrm{d}y$$ $$=\int_0^1\left( \int_0^1\frac{y^2}{y^2+z^2}\left(PV\int_0^{\infty}\frac{x^{s+1}}{1-y^2 x^2}\textrm{d}x\right) \textrm{d}z\right)\textrm{d}y$$ $$+\int_0^1\left( \int_0^1\frac{z^2}{y^2+z^2}\left(\int_0^{\infty}\frac{x^{s+1}}{1+z^2 x^2}\textrm{d}x\right) \textrm{d}z\right)\textrm{d}y$$ $$=\frac{1}{2}\int_0^1\left( \int_0^1\frac{y^{-s}}{y^2+z^2}\left(PV\int_0^{\infty}\frac{x^{s/2}}{1-x}\textrm{d}x\right) \textrm{d}y\right)\textrm{d}z$$ $$+\frac{1}{2}\int_0^1\left( \int_0^1\frac{z^{-s}}{y^2+z^2}\left(\int_0^{\infty}\frac{x^{s/2}}{1+x}\textrm{d}x\right) \textrm{d}y\right)\textrm{d}z$$ $$=\frac{\pi}{2}\cot\left(\frac{\pi}{2}s\right)\int_0^1\left( \int_0^1\frac{y^{-s}}{y^2+z^2} \textrm{d}z\right)\textrm{d}y-\frac{\pi}{2}\csc\left(\frac{\pi}{2}s\right)\int_0^1\left( \int_0^1\frac{z^{-s}}{y^2+z^2} \textrm{d}y\right)\textrm{d}z$$ $$=-\frac{\pi}{2}\tan\left(\frac{\pi}{4}s\right)\int_0^1\left( \int_0^1\frac{y^{-s}}{y^2+z^2} \textrm{d}z\right)\textrm{d}y=-\frac{\pi}{2}\tan\left(\frac{\pi}{4}s\right)\int_0^1 y^{-1-s}\left(\frac{\pi}{2}-\arctan(y)\right)\textrm{d}y$$ $$=-\frac{\pi}{2}\tan\left(\frac{\pi}{4}s\right)\left(-\frac{\pi}{4s}-\frac{1}{s}\int_0^1\frac{y^{-s}}{1+y^2}\textrm{d}y\right)$$ $$=\frac{\pi}{8s}\tan\left(\frac{\pi}{4}s\right)\left(\pi+\psi\left(\frac{3}{4}-\frac{s}{4}\right)-\psi\left(\frac{1}{4}-\frac{s}{4}\right)\right).$$
Now, based upon the previous result, it's easy to see that $$\lim_{s\to0}\frac{d}{ds}\{\Re\{J(s)\}\}=\lim_{s\to0}\frac{d}{ds}\biggr\{\Re\biggr\{\int_0^{\infty}x^{s-1}\operatorname{arctanh}(x)\arctan(x)\textrm{d}x\biggr \}\biggr\}$$ $$=\lim_{s\to0}\frac{d}{ds}\biggr\{\frac{\pi}{8s}\tan\left(\frac{\pi}{4}s\right)\left(\pi+\psi\left(\frac{3}{4}-\frac{s}{4}\right)-\psi\left(\frac{1}{4}-\frac{s}{4}\right)\right)\biggr\}$$ $$=\frac{3}{64}\zeta(2)\left(\psi^{(1)}\left(\frac{1}{4}\right)-\psi^{(1)}\left(\frac{3}{4}\right)\right)$$ $$=\frac{3}{4}\zeta(2)G,\tag4$$ which is immediately clear if we use the Trigamma series representation, $\displaystyle \psi^{(1)}(z)= \sum_{k=0}^{\infty} \frac{1}{(z+k)^2}$, and then recognize in the difference of the Trigamma special values the series representation of the Catalan's constant, $\displaystyle G=\sum_{k=0}^{\infty} \frac{(-1)^k}{(2k+1)^2}$.
Combining the results in $(3)$ and $(4)$, we conclude that
$$I=\int_0^1 \frac{\log(x)\operatorname{arctanh}(x)\arctan(x)}{x}\textrm{d}x=\frac{3}{8}\zeta(2)G-\frac{7}{32}\pi \zeta(3),$$ which is the desired result.
For example, using the same strategy we may get a generalization of the present integral. Another such exotic integral we may get is
$$\int_0^1 \frac{\log^3(x)\operatorname{arctanh}(x)\arctan(x)}{x}\textrm{d}x$$ $$=\frac{3}{1024}\zeta(2)\psi^{(3)}\left(\frac{1}{4}\right)-\frac{945}{256}\zeta(6)-\frac{93}{64}\pi\zeta(5)+\frac{45}{64}\zeta(4)G,$$ which looks really nice, isn't it?
A first note: The other integral in Shobhit Bhatnagar's post may be done in a similar style $\displaystyle \int_0^1 \frac{\log^2(x)\tanh^{-1}(x)}{1+x^2}\textrm{d}x,$ which can also be reduced to integrals already found in the book, (Almost) Impossible Integrals, Sums, and Series, like
$$\int_0^1 \frac{\arctan(x)\log^2(x)}{1+x} \textrm{d}x=\frac{21}{64}\pi \zeta(3)-\frac{\pi^3}{32}\log(2)-\frac{\pi^2}{24}G.$$
The other integral with $1-x$ in denominator may be calculated by a strategy similar to the one presented in this post.
A second note: Using the Cauchy product $\displaystyle \arctan(x)\operatorname{arctanh}(x)=\sum _{k=1}^{\infty} \sum _{n=1}^{2 k-1} \frac{(-1)^{n-1} x^{4 k-2}}{(2 k-1) (2 n-1)}$, is another way of attacking the integrals. For example, using the results with integrals previously given, we get other beautiful results with series
$$\sum_{k=1}^{\infty} \frac{1}{(2k-1)^3} \sum _{n=1}^{2 k-1}\frac{(-1)^{n-1}}{2 n-1}=\frac{7 }{8}\pi \zeta (3)-\frac{3 }{2}\zeta(2)G,$$
or
$$\sum_{k=1}^{\infty} \frac{1}{(2k-1)^5}\sum _{n=1}^{2 k-1}\frac{(-1)^{n-1}}{2 n-1}$$ $$=\frac{315}{32}\zeta (6)-\frac{15 }{8}\zeta(4)G+\frac{31 }{8}\pi \zeta (5)-\frac{1}{128} \zeta (2) \psi ^{(3)}\left(\frac{1}{4}\right).$$
The generalization of the main integral follows easily by employing the same ideas used in this post and the previous one.
Let $n$ be a natural number. Then, we have $$\int_0^1 \frac{\log^{2n}(x)\operatorname{arctanh}(x)}{1+x^2}\textrm{d}x$$ $$=\lim_{s\to0}\frac{d^{2n}}{ds^{2n}}\left(\frac{\pi}{16}\cot \left(\frac{\pi s}{2}\right) \left(\psi \left(\frac{3}{4}-\frac{s}{4}\right)-\psi\left(\frac{1}{4}-\frac{s}{4}\right)\right)-\frac{\pi ^2 }{16} \csc \left(\frac{\pi s}{2}\right)\right),$$ where $\psi$ represents the Digamma function.
Another similar generalization
Let $n$ be a natural number. Then, we get $$\int_0^1 \frac{\log^{2n}(x)\arctan(x)}{1-x^2}\textrm{d}x$$ $$=\frac{\pi}{4} \left(1-2^{-2 n-1}\right) \zeta (2 n+1)(2 n)!$$ $$-\lim_{s\to0}\frac{d^{2n}}{ds^{2n}}\left(\frac{\pi}{16} \csc \left(\frac{\pi s}{2}\right) \left(\pi \cos \left(\frac{\pi s}{2}\right)+\psi\left(\frac{s+1}{4}\right)-\psi\left(\frac{s+3}{4}\right)\right)\right),$$ where $\zeta$ represents the Riemann zeta function and $\psi$ denotes the Digamma function.
A solution in large steps by Cornel I. Valean to the main integral $$\int_0^1 \frac{\log^2(x)\operatorname{arctanh}(x)}{1+x^2}dx$$
We follow the strategy used to the auxiliary result from the previous post, and then we immediately arrive at
$$\int_0^1 \frac{\log^2(x)\operatorname{arctanh}(x)}{1+x^2}dx=\frac{1}{2}\Re\biggr\{ \int_0^{\infty } \frac{\log ^2(x) \operatorname{arctanh}(x)}{1+x^2} \textrm{d}x\biggr \}$$ $$=\frac{1}{2} \int_0^{\infty }\left(PV\int_0^1 \frac{x \log ^2(x)}{(1-y^2 x^2)(1+x^2)} \textrm{d}y\right)\textrm{d}x$$ $$=\frac{1}{2}\int_0^1\left(PV\int_0^{\infty} \frac{x \log ^2(x)}{(1-y^2 x^2)(1+x^2)} \textrm{d}x\right)\textrm{d}y$$ $$=\frac{\pi^2}{12}\int_0^1 \frac{\log(y)}{1+y^2}\textrm{d}y-\frac{1}{6}\int_0^1 \frac{\log^3(y)}{1+y^2}\textrm{d}y=\beta(4)-\frac{\pi^2}{12}G,$$ as desired.
End of story.
A note: Using the Cauchy product $\displaystyle \frac{\operatorname{arctanh}(x)}{1+x^2}=\sum _{n=1}^{\infty } \sum _{k=1}^n \frac{(-1)^{n+k} x^{2 n-1}}{2 k-1}$, an the value of the main integral, we immediately obtain the beautiful series
$$\sum _{n=1}^{\infty }\frac{(-1)^{n-1}}{n^3} \sum _{k=1}^n \frac{(-1)^{k-1}}{2 k-1}=4\beta(4)-\frac{\pi^2}{3}G.$$
Some kind of bonus: Using the integrals relation obained with integration by parts as shown in Shobhit Bhatnagar's post and combining it with the results obtained in this post and the previous one, we obtain the value of the other integral,
$$\int_0^1\frac{\log^2(x)\arctan(x)}{1-x^2}\textrm{d}x= -\beta(4)-\frac{\pi^2}{24}G+\frac{7\pi}{16}\zeta(3).$$
A note: It's clear the generalization $\displaystyle \int_0^1 \frac{\log^{2n}(x)\arctan(x)}{1-x^2}\textrm{d}x$ may be approached in the same way as $\displaystyle \int_0^1 \frac{\log^2(x)\operatorname{arctanh}(x)}{1+x^2}dx$.