Variance of the Watterson estimator ($\mathrm{Var}({\hat{\theta }}_W)=\frac{\theta }{H_{n-1}}+\frac{{\theta }^2{H}_{n-1}^{(2)}}{H_{n-1}^2}$)

The moments of the mutation count on a neutral coalescent tree give this variance for ${\hat{\theta }}_W=S/H_{n-1}$. Since $H_{n-1}\sim \log n$ and the sum of inverse squares is bounded, it tends to zero. Together with unbiasedness, this gives mean-square convergence to the mutation parameter, so it is a consistent estimator.