I am tempted to think 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.
