In doing a computation regarding the electroweak phase transition in the early universe, I came across the following integral: $$I=\int_0 ^ {\infty} x^2 \text{ln}(1-\text{exp}(-\sqrt{x^2+u^2}))\text{d}x. $$
I am interested in the behaviour of this parametric integral near $u=0$. Specifically, my reference asks me to show that
$$I = -\frac{\pi^4}{45} + \frac{\pi^2}{12}u^2 - \frac{\pi}{6}u^3 - \frac{1}{32}\bigg(\text{ln}(u^2)-(\frac{3}{2}+2 \text{ln}(4\pi)-2\gamma)\bigg) u^4+\mathcal{O}(u^6) $$
This expression is rather mysterious to me. I'm pretty sure I can obtain the first two terms by putting $u=0$ and differentiating with respect $u^2$, respectively, and using some values for $\zeta(n)$. The cubic term and the logarithm seem to be beyond my abilities.
Any help is appreciated.

I think that you can do what you need, if you carefully use the inequality (due to Yves Coudene (2018) A Strange Inequality Concerning Alternating Series, The American Mathematical Monthly, 125:6, 554-557.): $$ \ln(1+x)\le x - {x^2\over 2}+{x^3\over 3}+x^4 \left ( \ln(2)-1+{1\over2}-{1\over3} \right). $$