Solution (source code)

= Solution

Apply the <law of total variance> to each epoch. The <Poisson distribution> has equal conditional mean and <variance>, while $\operatorname{Var}(T_j)=\binom j2^{-2}$. Hence
$$
\operatorname{Var}(Y_j)=\mathbb E\!\left[\frac{j\theta T_j}{2}\right]+\operatorname{Var}\!\left(\frac{j\theta T_j}{2}\right)=\frac{\theta}{j-1}+\frac{\theta^2}{(j-1)^2}.
$$
The <independent mutation counts in coalescent epochs> permit <variance additivity for independent random variables>. Thus
$$
\boxed{\operatorname{Var}(S)=\theta H_{n-1}+\theta^2H_{n-1}^{(2)},\qquad H_{n-1}^{(2)}=\sum_{i=1}^{n-1}\frac1{i^2}.}
$$
The first term is conditional <Poisson distribution> variation; the second comes from the random lengths of the ancestral tree. These are the <moments of the mutation count on a neutral coalescent tree>.