Seeking for other methods to evaluate $\int_0^{\infty} \frac{\ln \left(x^n+1\right)}{x^n+1} dx$ for $n\geq 2$. Inspired by my post, I go further to investigate the general integral and find a formula for
$$
I_n=\int_0^{\infty} \frac{\ln \left(x^n+1\right)}{x^n+1} dx =-\frac{\pi}{n} \csc \left(\frac{\pi}{n}\right)\left[\gamma+\psi\left(1-\frac{1}{n}\right)\right] \tag*{}
$$
Let’s start with its partner integral
$$
I(a)=\int_0^{\infty}\left(x^n+1\right)^a d x
$$
and transform $I(a)$, by putting $y=\frac{1}{x^n+1}$, into a beta function
$$
\begin{aligned}
I(a) &=\frac{1}{n} \int_0^1 y^{-a-\frac{1}{n}-1}(1-y)^{-\frac{1}{n}-1} d y \\
&=\frac{1}{n} B\left(-a-\frac{1}{n}, \frac{1}{n}\right)
\end{aligned}
$$
Differentiating $I(a)$ w.r.t. $a$ yields
$$
I^{\prime}(a)=\frac{1}{n} B\left(-a-\frac{1}{n}, \frac{1}{n}\right)\left(\psi(-a)-\psi\left(-a-\frac{1}{n}\right)\right)
$$
Then putting $a=-1$ gives our integral$$
\begin{aligned}
I_n&=I^{\prime}(-1) \\&=\frac{1}{n} B\left(1-\frac{1}{n}, \frac{1}{n}\right)\left[\psi(1)-\psi\left(1-\frac{1}{n}\right)\right] \\
&=-\frac{\pi}{n} \csc \left(\frac{\pi}{n}\right)\left[\gamma+\psi\left(1-\frac{1}{n}\right)\right]\end{aligned}
$$

For examples,
$$
\begin{aligned}& I_2=-\frac{\pi}{2} \csc \frac{\pi}{2}\left[\gamma+\psi\left(1-\frac{1}{2}\right)\right]=\pi \ln 2,\\ & I_3=-\frac{\pi}{3} \csc \left(\frac{\pi}{3}\right)\left[\gamma+\psi\left(\frac{2}{3}\right)\right]=\frac{\pi \ln 3}{\sqrt{3}}-\frac{\pi^2}{9} ,\\ &I_4=-\frac{\pi}{4} \csc \left(\frac{\pi}{4}\right)\left[\gamma+\psi\left(\frac{3}{4}\right)\right]=\frac{3 \pi}{2 \sqrt{2}}\ln 2-\frac{\pi^2}{4 \sqrt{2}},\\
& I_5=-\frac{\pi}{5} \csc \left(\frac{\pi}{5}\right)\left[\gamma+\psi\left(\frac{4}{5}\right)\right]=-\frac{2 \sqrt{2} \pi}{5 \sqrt{5-\sqrt{5}}}\left[\gamma+\psi\left(\frac{4}{5}\right)\right], \\
& I_6=-\frac{\pi}{6} \csc \left(\frac{\pi}{6}\right)\left[\gamma+\psi\left(\frac{5}{6}\right)\right]=\frac{2 \pi}{3} \ln 2+\frac{\pi}{2} \ln 3-\frac{\pi^2}{2 \sqrt{3}},
\end{aligned}
$$

Furthermore, putting $a=-m$, gives
$$\boxed{I(m,n)=\int_0^{\infty} \frac{\ln \left(x^n+1\right)}{(x^n+1)^m} dx = \frac{1}{n} B\left(m-\frac{1}{n}, \frac{1}{n}\right)\left[\psi(m)-\psi\left(m-\frac{1}{n}\right)\right] }$$
For example,
$$
\begin{aligned}
\int_0^{\infty} \frac{\ln \left(x^6+1\right)}{(x^6+1)^5} dx  & =\frac{1}{6} B\left(\frac{29}{6}, \frac{1}{6}\right)\left[\psi(5)-\psi\left(\frac{29}{6}\right)\right] \\
& =\frac{1}{6} \cdot \frac{21505 \pi}{15552} \cdot\left(-\frac{71207}{258060}-\frac{\sqrt{3} \pi}{2}+\frac{3 \ln 3}{2}+2 \ln 2\right) \\
& =\frac{21505 \pi}{93312}\left(-\frac{71207}{258060}-\frac{\sqrt{3} \pi}{2}+\frac{3 \ln 3}{2}+2 \ln 2\right)
\end{aligned}
$$

Are there any other methods?
Your comments and alternative methods are highly appreciated.

For large $n$, a simple asymptotic behaviour of the integral can be deduced.
Let us first examine the structure of the integrand
$$y(x)=\frac{\ln \left(x^n+1\right)}{x^n+1} $$
Take the derivative of the function and set the result to zero to get the point $x=x_{m}$ where $y(x)$ reaches his maximum
$$1+(x_{m})^{n}=e$$
Putting this into $y(x)$ gives
$$y_{max}=y(x_{m})=\frac{\ln e}{e}=\frac{1}{e}$$
$y_{max}$ is independent of $n$
For integrals whose integrand, $y(x)$, has a sharp maximum a simple asymptotic formula exists.
$$I\approx \sqrt{2\pi\frac{y^{3}(x_{m})}{\left|y''(x_{m}) \right|}}$$
Here $y''(x_{m})$ is the second derivative of $y(x)$ at $x=x_{m}$.
I will skip elementary computations and write down final result
$$I_{n}\approx \frac{\sqrt{2\pi}}{n(e-1)^{1-\frac{1}{n}}}$$
Below is a few numerical examples to show approximation errors produced by the last formula
$n=10$,  the approximation error about $0.03$
$n=20$, the approximation error about $0.01$
$n=30$, the approximation error about $0.007$
The higher n is, the higher the accuracy of the last formula.
