This may only be a long comment. Since you are going to tune $\tilde{w} \to 0^+$, why not set it to $0$ right away? Since $W(\tilde{w}^2)$ is continuous and zero and $W(\tilde{w}^2) = \tilde{w}^2 - \tilde{w}^4 + O(\tilde{w}^6)$, this can be done in the second form of the integral you gave. With the change of variables $x=zy$, it becomes
$$
\Psi(z) = z^{1/2} \int_0^\infty \frac{\ln(1+e^{-zy}/y)}{\sqrt{y}} \, dy
+ z^{3/2} \int_0^\infty \sqrt{y} \ln(1+e^{-zy}/y) \, dy .
$$
You didn't say what is it that you want to know about $\Psi(z)$ (barring an "exact" formula, which seems unlikely). Depending on what you want, the above integral can be transformed and approximated in various ways.
Edit: For small $z\ll 1$, you need two separate (uniform) approximations, for small $y\ll 1/z$ and for large $y\gg 1$:
$$\tag{$*$}
\ln(1+e^{-zy}/y) \sim
\begin{cases}
\ln(1+1/y) - \frac{zy}{y+1} + O(zy)^2 & (y \ll 1/z) \\
e^{-zy}/y - e^{-2zy}/(2y^2) + O(y)^2 & (1 \ll y)
\end{cases}
$$
These approximations have an overlap region $1 \ll y \ll 1/z$. Then, picking a point within the overlap region, like $y=z^{-1/2}$, you can integrate one approximation over $(0,z^{-1/2}]$ and the other over $[z^{-1/2},\infty)$. If you only integrate one of them over the full $(0,\infty)$ interval, then you get divergences either at one or at the other end. The two separate integrals will still have divergent terms coming from the $y=z^{-1/2}$ cut, but they will cancel when added together.
The integrals with the logs and square roots can be computed using elementary functions (a computer algebra system should do this easily). The integrals involving exponentials will be expressible in terms of incomplete $\Gamma$ functions, whose asymptotic expansions can be found here. Putting the results together, I get
\begin{multline*}
\Psi(z)
= z^{1/2} (2\pi -2\sqrt{\pi} z^{1/2} + \pi z + \cdots)
+ z^{3/2} (\sqrt{\pi} z^{-1/2} + \frac{2}{3}\pi + \cdots) \\
= 2\pi z^{1/2} -\sqrt{\pi} z + \frac{5}{3} \pi z^{3/2} + O(z^{3/2}) .
\end{multline*}
Note that my estimate for the error term is just a guess (it could be $O(z^2)$ instead). It would need to be checked by computing the leading correction coming from the next order expansions in $(*)$.
For large $z\gg 1$, the same method doesn't quite work, since the approximations $(*)$ will be valid in the regions $y\ll 1/z$ and $1 \ll y$, which no longer overlap. Perhaps a third approximation that covers the intermediate region and overlaps with the other two could be found. Then the same procedure would work again.
