How far can we go with the integral $I_n=\int_0^1 \frac{\ln \left(1-x^n\right)}{1+x^2} d x$ Inspired by my post, I decided to investigate the integral in general
$$
I_n=\int_0^1 \frac{\ln \left(1-x^n\right)}{1+x^2} d x$$
by the powerful substitution $x=\frac{1-t}{1+t} .$
where $n$ is a natural number greater $1$.
Let’s start with easy one
\begin{aligned}
I_1 &=\int_0^1 \frac{\ln \left(\frac{2 t}{1+t}\right)}{1+t^2} d t \\ &=\int_0^1 \frac{\ln 2+\ln t-\ln (1+t)}{1+t^2} d t \\&=\frac{\pi}{4} \ln 2+\int_0^1 \frac{\ln t}{1+t^2}-\int_0^1 \frac{\ln (1+t)}{1+t^2} d t \\&=\frac{\pi}{4} \ln 2-G-\frac{\pi}{8} \ln 2 \\
&=\frac{\pi}{8} \ln 2-G\end{aligned}

By my post
$$I_2= \frac{\pi}{4} \ln 2-G $$
and
$$\begin{aligned}I_4 &=\frac{3 \pi}{4} \ln 2-2 G
\end{aligned}
$$
$$
\begin{aligned}
I_3=& \int_0^1 \frac{\ln (1-x)}{1+x^2} d x +\int_0^1 \frac{\ln \left(1+x+x^2\right)}{1+x^2} d x \\=& \frac{\pi}{8} \ln 2-G+\frac{1}{2} \int_0^{\infty} \frac{\ln \left(1+x+x^2\right)}{1+x^2} d x-G
\\ =& \frac{\pi}{8} \ln 2-\frac{4G}{3} +\frac{\pi}{6} \ln (2+\sqrt{3})
\end{aligned}
$$
where the last integral refers to my post.

Let’s skip $I_5$ now.
$$
I_6=\int_0^1 \frac{\ln \left(1-x^6\right)}{1+x^2} d x=\int_0^1 \frac{\ln \left(1+x^3\right)}{1+x^2} d x+I_3\\
$$
$$\int_0^1 \frac{\ln \left(1+x^3\right)}{1+x^2} d x = \int_0^1 \frac{\ln (1+x)}{1+x^2} d x +\int_0^1 \frac{\ln \left(1-x+x^2\right)}{1+x^2} d x\\=\frac{\pi}{8}\ln 2+ \frac{1}{2} \int_0^{\infty} \frac{\ln \left(1-x+x^2\right)}{1+x^2} d x-G \\= \frac{\pi}{8}\ln 2+ \frac{1}{2}\left( \frac{2 \pi}{3} \ln (2+\sqrt{3})-\frac{4}{3} G \right)- G \\= \frac{\pi}{8} \ln 2+\frac{\pi}{3} \ln (2+\sqrt{3})-\frac{5}{3} G $$
Hence $$I_6= \frac{\pi}{4} \ln 2+\frac{\pi}{2} \ln (2+\sqrt{3})-3 G$$

How far can we go with  the integral  $I_n=\int_0^1 \frac{\ln \left(1-x^n\right)}{1+x^2} d x?$

As @FShrike’s solution, I want to express the integral in terms of an infinite series of Diagamma functions.
Using $\ln \left(1-x^n\right)=-\sum_{k=1}^{\infty} \frac{x^{n k}}{k}$ for $|x|<1$, we have
$$
\int_0^1 \frac{\ln \left(1-x^n\right)}{1+x^2} d x=-\sum_{k=1}^{\infty} \frac{1}{k} \underbrace{\int_0^1 \frac{x^{n k}}{1+x^2} d x}_{I_k}
$$
$$
\begin{aligned}
I_k &=\int_0^1 \frac{y^{\frac{n k}{2}}}{1+y} \cdot \frac{1}{2 \sqrt{y}} d y \quad \textrm{, where }y=x^2\\
&=\frac{1}{2} \int_0^1 \frac{y \frac{n k-1}{2}}{1+y} d y
\end{aligned}
$$
As $\int_0^1 \frac{x^a}{1+x} d x=\frac{1}{2}\left[\psi\left(\frac{a}{2}+1\right)-\psi\left(\frac{a}{2}+\frac{1}{2}\right)\right]$, we can conclude that
$$
\boxed{
\int_0^1 \frac{\ln \left(1-x^n\right)}{1+x^2} d x =\frac{1}{4} \sum_{k=1}^{\infty} \frac{1}{k}\left[\psi\left(\frac{n k+1}{4}\right)-\psi\left(\frac{n k+3}{4}\right)\right]}
$$
