My problem is as follows. I have a reaction system $$X_1 \rightarrow X_2 \rightarrow X_3 \rightarrow ... \rightarrow X_N \rightarrow X_1$$ where each transition happens at an exponential rate with parameter $\lambda$. Obviously, it is straightforward to find the steady state distributions of particles in $X_1,X_2,...$. Further, the time it takes a particle to transit from $X_1$ back to $X_1$ is clearly $\tau \sim Erlang(\lambda, N)$. If we say, however, that a particle begins with age $0$ in stage $X_1$, and its age resets to $0$ every time it returns to stage $X_1$ - at steady state, what is the distribution of particle ages?
Any insight into this would be greatly appreciated. Obviously I could use a Monte-Carlo method to figure this out (and I have done), but I was wondering if there's a more analytic approach I could use.
In the steady state distribution, a random sample of the age of particles would follow a uniform distribution on the lifespan of each particle. Let $t$ denote the lifespan of particles, given by the Erlang distribution with $(\lambda,N)$ parameters and let $a$ be their age. Then the distribution of $a$ given $t$ should be Uniform on $[0,t]$. In other words it is a compound probability distribution. $$ p(a) = \int p(a|t)p(t)\,dt = \int_0^\infty \frac{1}{t}\mathbb I_{[0,t]}(a) \frac{\lambda^Nt^{N-1}e^{-\lambda t}}{(N-1)!} \, dt$$ Here $\mathbb I_{[a,b]}$ denotes the indicator function on $[a,b]$. Note that $\frac{1}{t}\mathbb I_{[0,t]}(a)$ can be rewritten as $\frac{1}{t}\mathbb I_{[a,\infty)}(t)$, so the bounds of integration only need to consider $[a,\infty)$. $$ p(a) = \int_{a}^{\infty} \frac{1}{t} \frac{\lambda^Nt^{N-1}e^{-\lambda t}}{(N-1)!} \,dt $$ $$ = \frac{\lambda}{N-2}\int_{a}^{\infty} \frac{\lambda^{N-1}t^{N-2}e^{-\lambda t}}{(N-2)!} \,dt $$ The integrand is now the density function of $Y\sim $ Erlang $(\lambda,N-1)$. $$ p(a) = \frac{\lambda}{N-2} \cdot P(Y\geq a) $$ $$ = \frac{\lambda e^{-\lambda a}}{N-2}\sum_{n=0}^{N-2}\frac{(\lambda a)^n}{n!}.$$