Solution (source code)

= Solution

With two offspring always, the backward equation becomes $\phi'=\mu(\phi^2-\phi)$. Separation and the initial value give
$$
\phi_t(s)=\frac{se^{-\mu t}}{1-s(1-e^{-\mu t})}=\sum_{k\ge1}e^{-\mu t}(1-e^{-\mu t})^{k-1}s^k.
$$
Thus $M_t$ is geometric on the positive integers, and
$$
\boxed{P(N_t=j)=e^{-\mu t}(1-e^{-\mu t})^j,\qquad j=0,1,\ldots.}
$$
This is a <Yule process> shifted down by one. It is \b[not an inhomogeneous <Poisson process>]. Given $N_t=j$, its instantaneous jump rate is $\mu(j+1)$, a state-dependent random rate rather than a deterministic function of time. Also $E[N_t]=e^{\mu t}-1$ whereas $\operatorname{Var}(N_t)=e^{\mu t}(e^{\mu t}-1)$, unequal to the mean for $t>0$. Its mean arrival intensity $\mu e^{\mu t}$ does not turn its dependent increments and geometric marginals into Poisson ones.