I am currently looking at the differential equations by Yule concerning the development of the number of species over time.
The (infinite) system of DEs for the probability $p_n(t)$ of $n$ species at time $t\in\mathbb R^+$, where $n\in \mathbb N$ and $r\in\mathbb R^+$, is given by Yule as: $$\frac{dp_1}{dt}=-rp_1,\\ \frac{dp_n}{dt}=(n-1)\,r\,p_{n-1}-n\,r\,p_n, \qquad n\geq2$$
If $p_1(0)=1$, the solution is given as a "geometric" probability: $$p_n(t)=e^{-rt}(1-e^{-rt})^{n-1}.$$
My question regards the interpretation. I do know, that $p_1(t)=exp(-rt)$ is the exponential probability that no new species has arrived until time $t$. In that sense, $p_1(t)$ is the survival function of an exponential distribution with parameter $r$.
But then, the case $n\geq 2$ becomes a bit more cumbersome, interpretationwise. It is kind of obvious that $(1-e^{-rt})^{n-1}$ is the cumulative distribution function of the maximum of $n-1$ independent exponentialy distributed random variables. But, $p_n(t)$ should be the probability that $n-1$ species existed somewhen at $s<t$. At that elusive time $s$ a new $n$-th species arrived, but no new species has arrived since (one would have to average over all possible times $s$). In other words, the next species will arrive somewhen after $t$.
I am trying to come up with a convincing explanation for the geometric nature of the probability, but I fail somewhat to put in into precise terms. Has anybody a clearer idea what's going on or can recommend a reference?