Trans series for an integral

$\color{brown}{\textbf{Solution.}}$

Subsitution $$x=u\sinh y\tag1$$ changes the issue integral to the form of $$I= u^3\int\limits_0^\infty\,\,\ln\left(1-e^{\Large^{\LARGE-u\cosh y}}\right)\sinh^2y\cosh y\,dy.\tag2$$

Integration by parts gives $$I= \dfrac13 u^3\int\limits_0^\infty\,\ln\left(1-e^{\Large^{\LARGE-u\cosh y}}\right)\,d\sinh^3y$$ $$ = \dfrac13 u^3 \ln\left(1-e^{\Large^{\LARGE-u\cosh y}}\right)\,\sinh^3y\bigg|_0^\infty\, -\dfrac13 u^4\int\limits_0^\infty\,\dfrac{\sinh^4y\,e^{\Large^{\LARGE-u\cosh y}}} {1-e^{\Large^{\LARGE-u\cosh y}}}\,dy$$ $$I=-\dfrac13 u^4\int\limits_0^\infty\,\dfrac{\sinh^4y}{e^{\Large^{\LARGE u\cosh y}}-1}\,dy.\tag3$$

On the other hand,

$$I=-\dfrac13 u^4\sum\limits_{k=1}^\infty I_k,\tag4$$ where $$I_k = \dfrac13 u^4\int\limits_0^\infty\,\sinh^4y\,e^{\Large^{\LARGE-ku\cosh y}}\,dy.\tag{5a}$$ If $u\to \infty,$ then the most part of the integral is situated in the area $y >> 1,$ where $$\sinh y \approx \cosh y.$$

So $$I_k \approx \dfrac13 u^4\int\limits_0^\infty\,\sinh^3y\cos y\,e^{\Large^{\LARGE-ku\cosh y}}\,dy\tag{5b},$$ $$I_k = -\dfrac{e^{\Large^{\LARGE-ku\cosh y}}}{12k^4} (24 + 2 k^2u^2 - ku (-24 + k^2 u^2) \cosh(y) + 6 u^2 \cosh(2 y) + k^3u^3 \cosh(3 y))\bigg|_0^\infty$$ $$ = \dfrac2{3k^4} e^{-ku} (k^2u^2 + 3 k u + 3),$$ and summation via polylogarithm changes $(3)$ to the expression $$I\approx -\dfrac23 u^2\operatorname{Li}_2\left(e^{-u}\right) - 2u\operatorname{Li}_3\left(e^{-u}\right) -2 \operatorname{Li}_4\left(e^{-u}\right)\tag6$$ (see also Wolfram Alpha summation), where $\operatorname{Li}$ is polylogarithm.

Similarly from $(3)$

$$I=-\dfrac13 u^4\sum\limits_{k=1}^\infty \dfrac{J_k}k,\tag7$$ where $$I_k = u^3\int\limits_0^\infty\,\sinh^2y\cosh y\,e^{\Large^{\LARGE-ku\cosh y}}\,dy.\tag{8a}$$

If $u\to \infty,$ then more accurate approximation can be used in the form of $$\sinh y\cosh y = \cosh^2y\sqrt{1-\cosh^{-2}y}\approx \cosh^2y\left(1-\dfrac12\cosh^{-2}y\right) = \cosh^2y-\dfrac12.$$

So $$J_k \approx \dfrac12u^3\int\limits_0^\infty\,(2\cosh^2 y -1)\sinh y \,e^{\Large^{\LARGE-ku\cosh y}}\,dy,\tag{8b}$$ $$J_k = \dfrac{k^2u^2 cos(2y)+4ku\cosh(y)+4}{2k^3} e^{\Large^{\LARGE-ku\cosh y}}\bigg|_0^\infty = \dfrac{(ku+2)^2}{2k^3} e^{-ku},$$ and summation via polylogarithm changes $(3)$ to the expression $$I\approx -\dfrac12 u^2 \operatorname{Li}_2\left(e^{-u}\right) - 2 u \operatorname{Li}_3\left(e^{-u}\right) - 2 \operatorname{Li}_4\left(e^{-u}\right).\tag9$$ (see also Wolfram Alpha summation), where $\operatorname{Li}$ is polylogarithm.

$\color{brown}{\textbf{Testing.}}$

\begin{vmatrix} \text{Formulas} & \text{OP (issue)} & (4)+(5b) & (7)+(8b) & (6) & (9)\\ \text{Range} & (0,150) & (0,12) &(0,12) & - & - \\ u=0.1 & -2.15691\,41055 & -2.15691\,41055 & -2.15691\,41055& -2.15938\,14465 & -2.15719\,44641\\ u=0.01 & -2.16456\,47397 & -2.16456\,47397 & -2.17457\,47397 & -2.16459\,18581 & -2.16456\,53771\\ \end{vmatrix}

1. Formula from the OP reference $\color{brown}{\textbf{has the error about 10x}}$ and looks wrong.
2. Approximation of $\sinh y$ in the integrals $(5b),(8b)$ works very accurate.
3. Both of the closed expressions $(6),(9)$ give a good accuracy.
4. The accuracy of the closed expressions improves when $u\to 0.$
5. Formula $(9)$ $$\color{green}{\boxed{\phantom{\Big|}\mathbf{I\approx -\dfrac12 u^2 \operatorname{Li}_2\left(e^{-u}\right) - 2 u \operatorname{Li}_3\left(e^{-u}\right) - 2 \operatorname{Li}_4\left(e^{-u}\right)}}\tag9}$$ gives the better approximation.

Partial answer. Regarding $I$ as a function of $u$, we note that

\begin{align*} I(0) &= \int_{0}^{\infty} x^2 \log(1-e^{-x})\, \mathrm{d}x = -\frac{\pi^4}{45}, \\ I'(0) &= 0, \\ I''(0) &= \int_{0}^{\infty} \frac{x}{e^x - 1} \, \mathrm{d}x = \frac{\pi^2}{6}. \end{align*}

To get the higher terms, write

\begin{align*} I(u) &= I(0^+) + \frac{I''(0)}{2}u^2 + \int_{0}^{\infty} \bigg( \underbrace{ x^2 \log\left(1 - e^{-\sqrt{x^2+u^2}}\right) - x^2 \log\left(1 - e^{-x}\right) - \frac{u^2 x}{2(e^x - 1)} }_{=\text{(*)}} \bigg) \, \mathrm{d}x. \end{align*}

Here, the integrand $\text{(*)}$ can be recast as

\begin{align*} \text{(*)} &= \int_{0}^{u} \left[ \frac{\partial}{\partial s} x^2 \log\left(1 - e^{-\sqrt{x^2+u^2}}\right) - \frac{sx}{(e^x - 1)} \right]\, \mathrm{d}s \\ &= \int_{0}^{u} \bigg[ \frac{sx^2}{\sqrt{x^2+s^2}(e^{\sqrt{x^2+s^2}}-1)} - \frac{sx}{(e^x - 1)} \bigg]\, \mathrm{d}s \\ &= \int_{0}^{u} \int_{0}^{s} s x^2 \frac{\partial}{\partial t} \bigg( \frac{1}{\sqrt{x^2+t^2}(e^{\sqrt{x^2+t^2}}-1)} \bigg) \, \mathrm{d}t \mathrm{d}s \\ &= -\int_{0}^{u} \int_{0}^{s} \frac{s x^2 t}{(t^2 + x^2)^2} f\left(\sqrt{x^2+t^2}\right) \, \mathrm{d}t\mathrm{d}s, \end{align*}

where $f(a) = \frac{a(e^a(1+a) - 1)}{(e^a - 1)^2}$. For the future use, we remark that $a = 0$ is a removable singularity of $f$ with $f(0) = 2$ and $2 \geq f(a) \geq f(b) \geq 0$ for all $0 \leq a \leq b$. Interchanging the order of integration,

\begin{align*} \text{(*)} &= -\int_{0}^{u} \frac{x^2 t (u^2-t^2)}{2(t^2 + x^2)^2} f\left(\sqrt{x^2+t^2}\right) \, \mathrm{d}t. \end{align*}

Plugging this back and applying the substitution $(x, t) \mapsto (utx, ut)$,

\begin{align*} I(u) &= I(0^+) + \frac{I''(0)}{2}u^2 - u^3 \int_{0}^{\infty} \int_{0}^{1} \frac{x^2 (1-t^2)}{2(1 + x^2)^2} f\left(u t \sqrt{x^2+1}\right) \, \mathrm{d}t\mathrm{d}x. \end{align*}

As $u \to 0^+$, monotone convergence tells that

\begin{align*} &\int_{0}^{\infty} \int_{0}^{1} \frac{x^2 (1-t^2)}{2(1 + x^2)^2} f\left(u t \sqrt{x^2+1}\right) \, \mathrm{d}t\mathrm{d}x \\ &\xrightarrow[u \to 0^+]{} \int_{0}^{\infty} \int_{0}^{1} \frac{x^2 (1-t^2)}{(1 + x^2)^2} \, \mathrm{d}t\mathrm{d}x = \frac{\pi}{6}. \end{align*}

So it follows that

$$ I(u) = -\frac{\pi^4}{45} + \frac{\pi^2}{12}u^2 - \frac{\pi}{6}u^3 + o(u^3). $$

Numerical analysis. If $0 < x < \frac{1}{2}$, then $\left| \log (1 - x) \right| = -\log(1-x) \leq \frac{x}{1-x} \leq 2x $. Using this, we note that, for $R > 1$,

$$ \left|\int_{R}^{\infty} x^2 \log\left( 1 - e^{-\sqrt{x^2 + u^2}} \right) \, \mathrm{d}x\right| \leq \int_{R}^{\infty} 2x^2 e^{-x} \, \mathrm{d}x = 2(R^2 + 2R + 2) e^{-R}. $$

For instance, this bound is $< 10^{-30}$ for $R \geq 80$, and so, we may compute

$$ \int_{0}^{\infty} x^2 \log\left( 1 - e^{-\sqrt{x^2 + u^2}} \right) \, \mathrm{d}x = \int_{0}^{80} x^2 \log\left( 1 - e^{-\sqrt{x^2 + u^2}} \right) \, \mathrm{d}x \pm 10^{-30}$$

uniformly in $u > 0$. Since the integrand is continuous on $[0, 80]$, we may use any CAS to easily estimate the truncated integral with a precision of $10^{-30}$. Then for the quantities

\begin{align*} \text{(A)} &= \text{numerical integration of } \int_{0}^{80} x^2 \log\left( 1 - e^{-\sqrt{x^2 + u^2}} \right) \, \mathrm{d}x, \\ \text{(B)} &= \text{OP's asymtotic form} \\ &= -\frac{\pi^4}{45} + \frac{\pi^2}{12}u^2 - \frac{\pi}{6}u^3 + \frac{\frac{3}{2} + 2\log(4\pi) - 2\gamma - 2\log u}{32} u^4, \\ \text{(C)} &= \text{@Yuri Negometyanov's asymptotic form} \\ &= -\frac{u^2}{2}\operatorname{Li}_2(e^{-u}) -2u \operatorname{Li}_3(e^{-u}) - 2\operatorname{Li}_4(e^{-u}), \end{align*}

we obtain the following comparisons of numerical values:

Coinciding leading decimals are highlighted as red.