How to find $ \int_0^\infty \dfrac x{1+e^x}\ dx$

$\newcommand{\+}{^{\dagger}} \newcommand{\angles}[1]{\left\langle\, #1 \,\right\rangle} \newcommand{\braces}[1]{\left\lbrace\, #1 \,\right\rbrace} \newcommand{\bracks}[1]{\left\lbrack\, #1 \,\right\rbrack} \newcommand{\ceil}[1]{\,\left\lceil\, #1 \,\right\rceil\,} \newcommand{\dd}{{\rm d}} \newcommand{\down}{\downarrow} \newcommand{\ds}[1]{\displaystyle{#1}} \newcommand{\expo}[1]{\,{\rm e}^{#1}\,} \newcommand{\fermi}{\,{\rm f}} \newcommand{\floor}[1]{\,\left\lfloor #1 \right\rfloor\,} \newcommand{\half}{{1 \over 2}} \newcommand{\ic}{{\rm i}} \newcommand{\iff}{\Longleftrightarrow} \newcommand{\imp}{\Longrightarrow} \newcommand{\isdiv}{\,\left.\right\vert\,} \newcommand{\ket}[1]{\left\vert #1\right\rangle} \newcommand{\ol}[1]{\overline{#1}} \newcommand{\pars}[1]{\left(\, #1 \,\right)} \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}} \newcommand{\pp}{{\cal P}} \newcommand{\root}[2][]{\,\sqrt[#1]{\vphantom{\large A}\,#2\,}\,} \newcommand{\sech}{\,{\rm sech}} \newcommand{\sgn}{\,{\rm sgn}} \newcommand{\totald}[3][]{\frac{{\rm d}^{#1} #2}{{\rm d} #3^{#1}}} \newcommand{\ul}[1]{\underline{#1}} \newcommand{\verts}[1]{\left\vert\, #1 \,\right\vert} \newcommand{\wt}[1]{\widetilde{#1}}$ \begin{align} &\color{#00f}{\large\int_{0}^{\infty}{x \over \expo{x} + 1}\,\dd x} =\int_{0}^{\infty}x\pars{{1 \over \expo{x} + 1} - {1 \over \expo{x} - 1}} \,\dd x + \int_{0}^{\infty}{x \over \expo{x} - 1}\,\dd x \\[3mm]&=-2\int_{0}^{\infty}{x \over \expo{2x} - 1}\,\dd x + \int_{0}^{\infty}{x \over \expo{x} - 1}\,\dd x =\half\int_{0}^{\infty}{x\expo{-x} \over 1 - \expo{-x}}\,\dd x \\[3mm]&=-\,\half\int_{0}^{\infty}\ln\pars{1 - \expo{-x}}\,\dd x =\half\int_{0}^{\infty}\sum_{n = 1}^{\infty}{\expo{-nx} \over n}\,\dd x =\half\sum_{n = 1}^{\infty}{1 \over n}\ \overbrace{\int_{0}^{\infty}\expo{-nx}\,\dd x}^{\ds{=\ {1 \over n}}} \\[3mm]&=\half\ \overbrace{\sum_{n = 1}^{\infty}{1 \over n^{2}}} ^{\ds{\zeta\pars{2} = {\pi^{2} \over 6}}} =\color{#00f}{\large{\pi^{2} \over 12}} \end{align}

$\ds{\zeta\pars{z}}$ is the Riemann Zeta Function.

One way to work this one out is to rewrite the integral as

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

then expand the denominator in a Taylor series. The representation of this integral becomes

$$\int_0^{\infty} dx \, x \sum_{k=0}^{\infty} (-1)^k e^{-(k+1) x} $$

Because the sum and integral converge, we can change the order of summation and integration:

$$\sum_{k=0}^{\infty} (-1)^k \int_0^{\infty} dx \, x \, e^{-(k+1) x} = \sum_{k=0}^{\infty} \frac{(-1)^k}{(k+1)^2}$$

Another way to work this, which you may or may not be aware of, is to use Cauchy's theorem. Consider the contour integral

$$\oint_C dz \frac{z^2}{1+e^z} $$

where $C$ is the rectangle with vertices $0$, $R$, $R+i 2 \pi$, and $i 2 \pi$, with a semicircular segment of radius $\epsilon$ centered at $i \pi$ into the rectangle. The contour integral then becomes

$$\int_0^R dx \frac{x^2}{1+e^x} + i \int_0^{2 \pi} dy \frac{(R+i y)^2}{1+e^R e^{i y}} +\int_R^0 dx \frac{(x+i 2 \pi)^2}{1+e^x} \\ +i \int_{2 \pi}^{\pi+\epsilon} dy \frac{(i y)^2}{1+e^{i y}}+ i \epsilon \int_{\pi/2}^{-\pi/2} d\phi \, e^{i \phi} \frac{(i \pi+\epsilon e^{i \phi})^2}{1-e^{\epsilon e^{i \phi}}} + i \int_{\pi-\epsilon}^0 dy \frac{(i y)^2}{1+e^{i y}}$$

As $R\to\infty$, the second integral vanishes. As $\epsilon \to 0$, the fifth integral becomes

$$i \epsilon \int_{\pi/2}^{-\pi/2} d\phi \, e^{i \phi} \frac{-\pi^2}{-\epsilon e^{i \phi}} = -i \pi^3$$

The contour integral is then

$$-i 4 \pi \int_0^{\infty} dx \frac{x}{1+e^x} + 4 \pi^2 \int_0^{\infty} \frac{dx}{1+e^x} + i PV \int_0^{2 \pi} dy \frac{y^2}{1+e^{i y}}-i \pi^3$$

Note that $PV$ denotes the Cauchy principal value. Also note that

$$\begin{align}i PV \int_0^{2 \pi} dy \frac{y^2}{1+e^{i y}} &= i \frac12 \int_0^{2 \pi} dy \, y^2 + \frac12 PV \int_0^{2 \pi} dy \, y^2 \tan{\left ( \frac{y}{2}\right )}\\ &= i \frac{4 \pi^3}{3} + \frac12 PV \int_0^{2 \pi} dy \, y^2 \tan{\left ( \frac{y}{2}\right )} \end{align}$$

By Cauchy's theorem, the contour integral is zero, which means that both the real and imaginary parts are zero. From the imaginary part of the contour integral being zero, we have

$$\int_0^{\infty} dx \frac{x}{1+e^x} = \frac{\pi^2}{12} $$

Another way to evaluate the integral. Rewrite $$ \int_0^{\infty} \frac{x}{1+e^{x}}\ dx=\int_0^{\infty} \frac{x \, e^{-x}}{1+e^{-x}}\ dx. $$ Using IBP by taking $u=x$ and $dv=-\dfrac{d(e^{-x})}{1+e^{-x}}$, you will obtain $$ \int_0^\infty\ln(1+e^{-x})\ dx. $$ Now, expand the integrand using Mclaurin's series for natural logarithm. After integrating the series, compare the result with Dirichlet eta function and you will obtain the integral is equal to $\eta(2)=\dfrac12\zeta(2)=\dfrac{\pi^2}{12}$.