How can I calculate moments of the time until the number of busy servers decreases to a certain level for the first time? Consider an $M/E_2/3/3$ queueing system with arrival rate $\lambda$ and average service time $1/\mu$. It is a loss system (no waiting room) with Erlang-2 service time distribution.  
Assume an arrival just occurred and increased the number of busy servers to 2. I am interested in moments of the time until the number of busy servers becomes 1 for the first time. 
Any ideas on how I can approach this problem?
 A: Let $\{X_n:n=0,1,2,\ldots\}$ be the embedded Markov chain, where the time points selected are immediately after the moment of departure. Note that after a departure the system cannot be full, so the state space is $\{0,1,2\}$. As in the general $M/G/1$ queue analysis, we compute the probability that there are $j$ arrivals during a service time:
$$
a_j = \int_0^\infty e^{-\lambda t}\frac{(\lambda t)^j}{j!} \mu(\mu t)e^{-\mu t}\ \mathsf dt = (j+1)\left(\frac\mu\lambda\right)^2\left(1+\frac\mu\lambda\right)^{-(j+2)},\ j=0,1,2.
$$
The transition matrix is given by
$$
P = \begin{bmatrix}a_0&a_1&1-(a_0+a_1)\\a_0&a_1&1-(a_0+a_1)\\0&a_0&1-a_0\end{bmatrix} = \begin{bmatrix}
\left(\frac\mu{\lambda+\mu}\right)^2 & \frac{2\lambda\mu^2}{(\lambda+\mu)^3} & \frac{\lambda^3+3\lambda^2\mu}{(\lambda+\mu)^3}\\
\left(\frac\mu{\lambda+\mu}\right)^2 & \frac{2\lambda\mu^2}{(\lambda+\mu)^3} & \frac{\lambda^3+3\lambda^2\mu}{(\lambda+\mu)^3}\\
0 &\left(\frac\mu{\lambda+\mu}\right)^2 & \frac{\lambda ^2+2 \lambda  \mu }{(\lambda +\mu )^2}
\end{bmatrix}.
$$
Let $X(t)$ be the number of customers in the system at time $t$. Suppose $X(0)=2$. Let $\tau = \inf\{t>0:X(t)=1\}$. Then $\tau = Z + S$ where $S$ is a service time and $Z$ is a geometric sum of arrivals and services, that is, $Z=\sum_{n=1}^N (A_n+S_n)$ where $\mathbb P(N=n) = \left(\frac\lambda{\lambda+\mu}\right)^n\left(\frac\mu{\lambda+\mu}\right)$. Now, the distribution of the sum of an arrival and a service is given by convolution:
$$
f_{A+S}(t) = f_A\star f_S(t) = \int_0^t \lambda e^{-\lambda} \tau\mu(\mu(t-\tau))e^{-\mu(t-\tau)}\ \mathsf d\tau = \frac{\lambda  \mu ^2 e^{-t (\lambda +\mu )} \left(e^{\lambda  t} (t (\lambda -\mu )-1)+e^{\mu  t}\right)}{(\lambda -\mu )^2}.
$$
This has characteristic function
\begin{align}
\varphi_{A+S}(\theta) &= \mathbb E[e^{i\theta(A+S)}]\\ &= \int_0^\infty e^{i\theta t}\int_0^t \lambda e^{-\lambda} \tau\mu(\mu(t-\tau))e^{-\mu(t-\tau)}\ \mathsf d\tau\\ &= \frac{\lambda  \mu ^2 e^{-t (\lambda +\mu )} \left(e^{\lambda  t} (t (\lambda -\mu )-1)+e^{\mu  t}\right)}{(\lambda -\mu )^2}\\ &= -\frac{i \lambda  \mu ^2}{(\theta +i \lambda ) (\theta +i \mu )^2}.
\end{align}
Therefore $Z$ has characteristic function
\begin{align}
\varphi_Z(\theta) &= \mathbb E[\varphi_{A+S}(\theta)^N]\\
&= \sum_{n=0}^\infty \mathbb E[\varphi_{A+S}(\theta)^N\mid N=n]\mathbb P(N=n)\\
&= \sum_{n=0}^\infty\left(-\frac{i \lambda  \mu ^2}{(\theta +i \lambda ) (\theta +i \mu )^2}\right)^n \left(\frac\lambda{\lambda+\mu}\right)^n\left(\frac\mu{\lambda+\mu}\right)\\
&= \frac{\mu  (\theta +i \lambda ) (\theta +i \mu )^2}{\theta ^3 (\lambda +\mu )+i \theta ^2 (\lambda +\mu ) (\lambda +2 \mu )-\theta  \mu  (\lambda +\mu ) (2 \lambda +\mu )-i \lambda  \mu ^3}.
\end{align}
The characteristic function of $\tau=Z+S$ is given by the product
\begin{align}
\varphi_\tau(\theta) &= \varphi_Z(\theta)\varphi_S(\theta)\\
&= \left(\frac{\mu  (\theta +i \lambda ) (\theta +i \mu )^2}{\theta ^3 (\lambda +\mu )+i \theta ^2 (\lambda +\mu ) (\lambda +2 \mu )-\theta  \mu  (\lambda +\mu ) (2 \lambda +\mu )-i \lambda  \mu ^3}\right)\left(\frac{\mu ^2}{(\mu -i \theta )^2}\right)\\
&=-\frac{\mu ^3 (\theta +i \lambda )}{\theta ^3 (\lambda +\mu )+i \theta ^2 (\lambda +\mu ) (\lambda +2 \mu )-\theta  \mu  (\lambda +\mu ) (2 \lambda +\mu )-i \lambda  \mu ^3}.
\end{align}
The moments of $\tau$ are given by
$$
\mathbb E[\tau^n] = i^{-n}\lim_{\theta\to 0}\varphi^{(n)}_\tau(s).
$$
For example, the mean is
\begin{align}
\mathbb E[\tau] &= -i\cdot\lim_{\theta\to 0}\frac{\mathrm d}{\mathrm d \theta}\left[-\frac{\mu ^3 (\theta +i \lambda )}{\theta ^3 (\lambda +\mu )+i \theta ^2 (\lambda +\mu ) (\lambda +2 \mu )-\theta  \mu  (\lambda +\mu ) (2 \lambda +\mu )-i \lambda  \mu ^3}\right]\\
&= -i\cdot\lim_{\theta\to 0}\frac{\mu ^3 \left(i \mu ^2 \left(2 \theta ^2+4 i \theta  \lambda -3 \lambda ^2\right)+2 \mu  (\theta +i \lambda )^3+2 \theta  \lambda  (\theta +i \lambda )^2\right)}{\left(\theta ^3 (-(\lambda +\mu ))-i \theta ^2 (\lambda +\mu ) (\lambda +2 \mu )+\theta  \mu  (\lambda +\mu ) (2 \lambda +\mu )+i \lambda  \mu ^3\right)^2}\\
&= -i\cdot\left(-\frac{-2 i \lambda ^3 \mu -3 i \lambda ^2 \mu ^2}{\lambda ^2 \mu ^3}\right)\\
&= \frac{2 \lambda +3 \mu }{\mu ^2}.
\end{align}
