I think the statement in the OP that $W_2(a)$ and $W_3(a)$ remain bounded when $a\rightarrow 0$ is mistaken. The three roots of $(w+1)^3+a^3=0$ are 
$$w_1= -a-1,\;\; w_2= \frac{1}{2} \left(-i \sqrt{3} a+a-2\right),\;\;w_3= \frac{1}{2} \left(i \sqrt{3} a+a-2\right).$$
Then the denominator $(w+1)^2$ evaluates to
$$(w_1+1)^2=a^2,\;\;(w_2+1)^2=-\frac{1}{4} \left(2+2 i \sqrt{3}\right) a^2,\;\;(w_3+1)^2=\frac{1}{2} i \left(\sqrt{3}+i\right) a^2,$$
so this vanishes for all three roots when $a\rightarrow 0$, while the numerator remains finite (equal to $-\gamma_{\rm Euler}$.

And indeed, a numerical check indicates that the Mathematica output is actually correct, and the erroneous numerical result for small $a$ is a numerical instability in the computation of the digamma function. See these two plots that compare the digamma expression (blue) with a numerical evaluation of the sum (gold), as a function of $a$. For $a\gtrsim 0.01$ the two answers are nearly indistinguishable (blue and gold coincide on the scale of the plot), and the deviations for smaller $a$ are widely oscillating, indicating a numerical instability.

<IMG SRC="https://ilorentz.org/beenakker/MO/digammacomparison.png" WIDTH="400"/>
<IMG SRC="https://ilorentz.org/beenakker/MO/digamma_comparison_2.png" WIDTH="400"/>