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?
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}