Evaluate $\int_0^{\infty}\frac{\log x}{1+e^x}\,dx$

There's a way avoiding special functions and transforms, similar to the device for calculating Frullani integrals: Let's calculate $\displaystyle\int^\infty_\epsilon\frac{\log x}{1+e^x}\,dx$ up to terms $o(1)$ for $\epsilon\to0+$. With the identity $$\frac1{e^x+1}=\frac1{e^x-1}-\frac2{e^{2x}-1},$$ we have \begin{align}\int^\infty_\epsilon\frac{\log x}{1+e^x}\,dx&=\int^\infty_\epsilon\frac{\log x}{e^x-1}\,dx-\int^\infty_\epsilon\frac{2\log x}{e^{2x}-1}\,dx\\&=\int^\infty_\epsilon\frac{\log x}{e^x-1}\,dx-\int^\infty_{2\epsilon}\frac{\log x-\log2}{e^x-1}\,dx \\&=\int^{2\epsilon}_\epsilon\frac{\log x}{e^x-1}\,dx+\int^\infty_{2\epsilon}\frac{\log2}{e^x-1}\,dx \end{align} Using $\displaystyle\frac1{e^x-1}=\frac1x+O(1)$ and $\displaystyle\int^{2\epsilon}_\epsilon|\log x|\,dx=o(1)$, we see $$\int^{2\epsilon}_\epsilon\frac{\log x}{e^x-1}\,dx=\int^{2\epsilon}_\epsilon\frac{\log x}{x}\,dx+o(1)=\frac12\log^22+\log2\log\epsilon+o(1)$$ and $$\int^\infty_{2\epsilon}\frac{\log2}{e^x-1}\,dx=\log2\log\frac1{1-e^{-2\epsilon}}=-\log^22-\log2\log\epsilon+o(1),$$ so $$\int^\infty_\epsilon\frac{\log x}{1+e^x}\,dx=-\frac12\log^22+o(1),$$ and our integral is $$\int^\infty_0\frac{\log x}{1+e^x}\,dx=-\frac12\log^22.$$


By the inverse Laplace transform $$ \sum_{n\geq 1}\frac{(-1)^{n+1}}{n^s} = \frac{1}{\Gamma(s)}\int_{0}^{+\infty}\frac{x^{s-1}}{e^x+1}\,dx $$ and by differentiating both sides with respect to $s$ $$ \sum_{n\geq 1}\frac{(-1)^n \log n}{n^s} = -\frac{\Gamma'(s)}{\Gamma(s)^2}\int_{0}^{+\infty}\frac{x^{s-1}}{e^x+1}\,dx + \frac{1}{\Gamma(s)}\int_{0}^{+\infty}\frac{x^{s-1}\log(x)}{e^x+1}\,dx $$ so by evaluating at $s=1$ $$\int_{0}^{+\infty}\frac{\log x}{e^x+1}\,dx = \sum_{n\geq 1}\frac{(-1)^n\log n}{n}+\underbrace{\Gamma'(1)}_{-\gamma}\underbrace{\int_{0}^{+\infty}\frac{dx}{e^x+1}}_{\log 2} $$ and it just remains to crack the mysterious series $\sum_{n\geq 1}\frac{(-1)^n\log n}{n}$. On the other hand by Frullani's integral, the inverse Laplace transform or Feynman's trick we have $\log(n)=\int_{0}^{+\infty}\frac{e^{-x}-e^{-nx}}{x}\,dx$, so $$\sum_{n\geq 1}\frac{(-1)^n\log n}{n}=\int_{0}^{+\infty}\frac{\log(1+e^{-x})-e^{-x}\log 2}{x}\,dx=\gamma\log(2)-\frac{1}{2}\log^2(2)\tag{J}$$ where the last identity follows from the integral representation for the Euler-Mascheroni constant, got by applying the inverse Laplace transform to the series definition $\gamma=\sum_{n\geq 1}\left[\frac{1}{n}-\log\left(1+\frac{1}{n}\right)\right]$. Summarizing, we simply have $$ \int_{0}^{+\infty}\frac{\log(x)}{e^x+1}\,dx = \color{red}{-\frac{1}{2}\log^2(2)}.$$ It is possible to prove the equality between the LHS and the RHS of $(J)$ by summation by parts and Euler sums, too.


As pointed out in the comments, let $I(a)$ be the following integral: $$I(a)=\int\limits_0^{\infty} \frac{x^a}{1+e^x}\,dx$$ We are then looking for $I'(0)$.

$$\begin{align} I(a)&=\int\limits_0^{\infty}e^{-x}x^a\frac{1}{1+e^{-x}}\,dx\\ I(a)&= \int\limits_0^{\infty}e^{-x}x^a \sum_{n=0}^{+\infty}(-1)^n e^{-nx}\,dx\\ I(a)&=\sum_{n=0}^{+\infty}(-1)^n\int_\limits0^{\infty}e^{-x(n+1)}x^a\,dx\\ I(a)&=\sum_{n=1}^{+\infty}(-1)^{n-1}\int\limits_0^{\infty}e^{-nx}x^a\,dx &u=nx\\ I(a)&=\sum_{n=1}^{+\infty}(-1)^{n-1}\int\limits_0^{\infty} e^{-u}\frac{u^a}{n^a}\frac{du}{n}\\ I(a)&=\sum_{n=1}^{+\infty}\frac{(-1)^{n-1}}{n^{a+1}} \int\limits_0^{\infty}e^{-u}u^{a}\,du\\ I(a)&=\eta(a+1)\Gamma(a+1) \end{align}$$ Where $\eta(s)$ is the Dirichlet eta function and $\Gamma(s)$ is the Gamma function. Taking the derivative of both sides, $$\begin{align} I'(a)&=\eta'(a+1)\Gamma(a+1)+\eta(a+1)\psi(a+1)\Gamma(a+1)\\ I'(0)&=\eta'(1)\Gamma(1)+\eta(1)\psi(1)\Gamma(1)\\ I'(0)&=\log(2)\gamma-\frac{1}{2}\log^2(2)-\log(2)\gamma\\ I'(0)&=-\frac{1}{2}\log^2(2) \end{align}$$ Thus,

$$\int\limits_0^{\infty} \frac{\log(x)}{1+e^x}\,dx=-\log^2\left(\sqrt{2}^{\sqrt{2}}\right)$$