Moments of the mutation count on a neutral coalescent tree ($\mathbb{E}S=\theta H_{n-1},\mathrm{Var}S=\theta H_{n-1}+{\theta }^2{H}_{n-1}^{(2)}$)

The independent mutation counts in coalescent epochs have the geometric mutation count in a coalescent epoch laws. Summing their expectations and variances gives the displayed formulas, where $H_m={\sum }_{i=1}^{m}1/i$ and $H_m^{(2)}={\sum }_{i=1}^{m}1/i^2$. The extra quadratic term records randomness of the ancestral tree, beyond the conditional Poisson distribution variance.