As I show in my
Perl::PDQ book,
the residence time at an M/M/1 queue is trivial to derive and (unlike most queueing theory texts) does not require any probability theory arguments.
Great for Guerrillas!
However, by simply adding
another server (i.e., M/M/2), that same Guerrilla approach falls apart. This situation has always bothered me profoundly and on several occasions I thought I saw how to get to the exact formula—the
Erlang C formula—Guerrilla style.
But, on later review, I always found something wrong.
Although I've certainly had correct pieces of the puzzle, at various times, I could never get everything to fit in a completely consistent way.
No matter how creative I got,
I always found a fly in the ointment.
The best I had been able to come up with is what I call the "morphing model" approximation where you start out with $m$ parallel queues at low loads and it morphs into a single $m$times faster M/M/1 queue at high loads.
That model is also exact for $m = 2$ servers—which is some kind of progress, but not much. Consequently, despite a few misplaced enthusiastic announcements in the past, I've never been able to publish the fully corrected morphing model.
Previous misfires have included:
 Falsely claimed in this 2008 blog post.
There, I show a Table of how complicated the Erlang B and C functions are when expressed as rational functions of
$\rho$ (the perserver utilization) for $m = 1, 2, \ldots, 6$ service facilities.
 As a side note, it's precisely those impenetrable polynomials with unfathomable coefficients that completely put me off the approach I am about to describe here.
 In the Comments section of that same post, Boris Solovyov asked in 2013: "Did you ever write this out in full?"
I sheepishly had to report that I hadn't had time to pursue it (which is mostly true). I hope he reads this post.
 A more intuitive explanation of the motivation for the morphing model was given in this 2011 blog post, but no advance on the long soughtafter correction terms.
 Falsely claimed again in this 2012 blog post. There's a photo of one of my whiteboards, deliberately kept inscrutably small—a good choice, as it turned out.
I think I know how Kepler must've felt. The difference between an M/M/1 queue and an M/M/m queue is like the difference
between a circle and an ellipse. Just by changing the width slightly, the calculation of the perimeter becomes enormously
complicated—in fact, it doesn't even have an analytic solution!
Now, however, I believe I finally have the correct approach for sure. Really!... Yep... Uh huh... No, seriously. Well, see for yourself
The Starting Point
Start with the
morphing approximation for the M/M/m waiting time in my
Perl::PDQ book—Chap. 2 in the original 2004 edition or Chap. 4 in the 2nd 2011 edition.
\begin{equation}
\Phi_W^{approx} (m, \rho) \, = \, \frac{1}{\rho^{m} \, \sum_{k=0}^{m1} \rho^k} \label{eqn:morphfun}
\end{equation}
This definition has a denominator involving a truncated geometric series in $\rho$.
It is used to derive the Guerrilla approximation for the residence time at an M/M/m queue with mean service time $S$
\begin{equation*}
R^{approx} \, = \, \frac{S}{1  \rho^m}
\end{equation*}
which is exact for $m = 1, 2$.
The exact general formula for the residence time is given by
\begin{equation*}
R^{exact} \, = \, S \; + \; \frac{C(m, \rho)}{m( 1  \rho)} \, S
\end{equation*}
where $C(m, \rho)$ is the Erlang C function defined in (\ref{sec:realEC}).
Five Easy Pieces
The following 5 steps provide the corrections to the approximation in (\ref{eqn:morphfun}) and result in the exact Erlang C function without resorting to the usual probability arguments found in almost all queueing theory textbooks. Along the way,
and quite unexpectedly (for me), I derive the Erlang B function as an intermediate result. Originally, I thought I would have to introduce that function as a kind of axiomatic probability. I wasn't thinking of it as a major waypoint.
As far as I'm aware, this derivation has never been presented before because you have to be sufficiently perverse to have thought up the morphing approximation in the first place.

Replace $\rho$ in (\ref{eqn:morphfun}) by the traffic intensity $a = m \rho$
\begin{equation}
\Phi_W^{approx} (m, a) \, = \, \frac{1}{a^{m} \, \sum_{k=0}^{m1} ~a^k} \label{eqn:morpha}
\end{equation}

Convert $\Phi_W^{approx}$ to a truncated exponential series by applying $a^n \mapsto a^n / n!$
\begin{equation}
\frac{1}{m! a^{m} \, \sum_{k=0}^{m1}~\frac{a^k}{k!}}
\end{equation}

Extend the summation over all $m$ servers to yield the famous
Erlang B function for an M/M/m/m queue
\begin{equation}
B(m, a) \, = \, \frac{\frac{a^m}{m!}}{\sum_{k=0}^{m}~\frac{a^k}{k!}} \label{eqn:EB}
\end{equation}
Historically, $B(m, a)$ has been associated with the probability that an incoming telephone call is blocked and lost from the system, e.g., gets an engaged signal. I can't remember the last time I heard an engaged signal. I think they've been replaced by voicemail and Muzak.

Using the more compact notation: $A_m = a^m / m!$ and $\Sigma_k$ for the reduced sum over $(m  1)$ terms, rewrite (\ref{eqn:EB}) as
\begin{equation}
B(m, a) \, = \, \frac{A_m}{A_m + \Sigma_k}
\end{equation}

Scale $A_m$ by $(1  \rho)^{1}$ to introduce the infinite possible wait states and arrive at
\begin{equation}
\Phi_W^{exact} (m, a) \, = \, \frac{1}{1 \, + \, (1  \rho) \, m! \, a^{m} \, \sum_{k=0}^{m1}~\frac{a^k}{k!}} \label{eqn:myEC}
\end{equation}
which is the fully corrected version of (\ref{eqn:morphfun}). Q.E.D.
Furthermore, it can be shown that (\ref{eqn:myEC}) is identical to the famous
Erlang C function
\begin{equation}
C(m, a) \, = \, \frac{ \frac{a^m}{m!} \, \big(\frac{m}{ma}\big) }{1 \, + \, \frac{a}{1!} \, + \, \frac{a^2}{2!} \, + \, \cdots \, + \, \frac{a^{m1}}{(m1)!} \, + \, \frac{a^m}{m!} \, \big(\frac{m}{ma}\big) } \label{sec:realEC}
\end{equation}
Historically, $C(m, a)$ has been associated with the probability that an incoming telephone call must wait to get a connection (aka "call waiting"), rather than being dropped, as it is in B(m, a). However, since $C(m, a)$ determines the mean waiting time in any M/M/m queue, it applies to any multiserver system, e.g., modern multithreaded applications.
Equation (\ref{sec:realEC}) was first published by A. K. Erlang in 1917 (he of the Copenhagen Telephone Company), along with (\ref{eqn:EB}). Thus, next year will be its centennial. Nice timing on my part (although that was never the plan—there never being any plan).
It's worth emphasizing that the morphing approximation (\ref{eqn:morphfun}) accounts for about 90% of what is going on with $R^{exact}$ in the M/M/m queue. The remaining 10% contains the minutiae regarding how the waiting line actually forms. But, as you can see from the above transformations, it's a rather subtle 10%.
Next Steps
This is just a sketch proof. I've suppressed a lot of details because there are many and I have a 100 pages of typeset notes to back that up.
I know a lot about what
doesn't work. (Sigh!)
Now that I have the correct mathematical logic sketched out, I'm quite confident that it can also be supplemented with a more visual representation of how the
corrections to the morphing function (\ref{eqn:morphfun}) arise.