What does a customer see when it begins to be served in $M/M/1$ queue? In queueing theory, the PASTA (Poisson Arrivals See Time Averages) principle [wiki] justifies $a_n = P_n$ where 
$$a_n = \text{proportion of customers that find } n \text{ customers in the system when they arrive}$$ and,
$$P_n = \lim_{t \to \infty} P \{ X(t) = n \}$$ ($X(t) \text{ here denotes the number of customers in system at time } t$).
However, I am now more interested in the moment when a customer begins to be served, rather than when it arrives.     
Specifically, I focus on the case of $M/M/1$ queue (i.e., a single-server exponential queueing system with FCFS service discipline) and consider
$$s_n = \text{proportion of customers that find } n \text{ in the } M/M/1 \text{ queue when it begins to be served.}$$ 
My question is:

Question: What does a customer see when it begins to be served in $M/M/1$ queue with the FCFS service discipline?     
In other words, what is the relationship between $s_n$ and $P_n$ (or, the relationship between $s_n$ and $a_n$)?

 A: Let $X$ denote the number of customers in the queue when a new customer arrives, the new customer excepted and the customer being served included, if any. Let $T$ denote the time it takes to serve these $X$ customers. Let $N$ denote the number of customers arriving during a time interval of length $T$. One asks for the distribution of $N$.
For a stationary M/M/1 queue, $X$ is geometrically distributed with parameter $\rho=\lambda/\mu$, that is, $P(X=x)=\rho^x(1-\rho)$ for every integer $x\geqslant0$. Conditionally on $X=0$, $T=0$. Conditionally on $X=x$ with $x\geqslant1$, the distribution of $T$ is gamma with parameters $(x,\mu)$. Conditionally on $T=0$, $N=0$. Conditionally on $X=x$ and $T=t$ with $t\gt0$, the distribution of $N$ is Poisson with parameter $\lambda t$. 
This decomposition in terms of conditional distributions yields the distribution of $N$ through careful computations. Namely, for every positive integer $n$, $$P(N=n)=\sum_{x=1}^\infty\rho^x(1-\rho)\int_0^\infty\frac{\mu^x}{\Gamma(x)}t^{x-1}\mathrm e^{-\mu t}\mathrm e^{-\lambda t}\frac{(\lambda t)^n}{n!}\mathrm dt,$$ that is, $$P(N=n)=\frac{1-\rho}{n!}\sum_{x=1}^\infty\frac{\lambda^{x+n}}{\Gamma(x)}\int_0^\infty t^{x+n-1}\mathrm e^{-(\lambda+\mu)t}\mathrm dt,
$$
or, equivalently, $$P(N=n)=\frac{1-\rho}{n!}\sum_{x=1}^\infty\frac{\lambda^{x+n}}{\Gamma(x)}\frac{\Gamma(x+n)}{(\lambda+\mu)^{x+n}}=\frac{1-\rho}{n!}\sum_{x=1}^\infty\frac{\Gamma(x+n)}{\Gamma(x)}\frac1{\sigma^{x+n}},\qquad\sigma=1+\frac1\rho.
$$
Note that, for every $\sigma\gt1$, $$\sum_{x=1}^\infty\frac{\Gamma(x+n)}{\Gamma(x)}\frac1{\sigma^{x+n}}=(-1)^n\frac{\mathrm d^n}{\mathrm d\sigma^n}\sum_{x=1}^\infty\frac1{\sigma^{x}}=(-1)^n\frac{\mathrm d^n}{\mathrm d\sigma^n}\left(\frac1{\sigma-1}\right)=\frac{n!}{(\sigma-1)^{n+1}},
$$
hence, for every $n\geqslant1$, $$P(N=n)=(1-\rho)\rho^{n+1},
$$ and, considering the complement to $1$ of the sum, $$P(N=0)=(1-\rho)(1+\rho).$$
One can note that $N$ is distributed as $\max\{G-1,0\}$ where $G$ is geometrically distributed with parameter $\rho$ in the sense that $P(G=n)=(1-\rho)\rho^n$ for every $n\geqslant0$.
