Solution (source code)

= Solution

Conditioning on the epoch length gives the <probability generating function>
$$
\mathbb E[s^{Y_j}\mid T_j]=\exp\left[-\frac{j\theta}{2}(1-s)T_j\right].
$$
The <Laplace transform of an exponential distribution> then yields
$$
\mathbb E[s^{Y_j}]=\frac{\binom j2}{\binom j2+(j\theta/2)(1-s)}=\frac{j-1}{j-1+\theta(1-s)}=\frac{p_j}{1-(1-p_j)s},\qquad p_j=\frac{j-1}{j-1+\theta}.
$$
This is the <geometric mutation count in a coalescent epoch>. With the zero-based <geometric distribution> convention,
$$
\boxed{P(Y_j=k)=p_j(1-p_j)^k,\qquad k=0,1,\ldots.}
$$
At $\theta=0$, this reduces to $Y_j=0$ almost surely.