Transient analysis of the Erlang A model

We consider the Erlang A model, or M/M/m+M\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M/M/m+M$$\end{document} queue, with Poisson arrivals, exponential service times, and m parallel servers, and the property that waiting customers abandon the queue after an exponential time. The queue length process is in this case a birth–death process, for which we obtain explicit expressions for the Laplace transforms of the time-dependent distribution and the first passage time. These two transient characteristics were generally presumed to be intractable. Solving for the Laplace transforms involves using Green’s functions and contour integrals related to hypergeometric functions. Our results are specialized to the M/M/∞\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M/M/\infty $$\end{document} queue, the M / M / m queue, and the M / M / m / m loss model. We also obtain some corresponding results for diffusion approximations to these models.


Introduction
In many real-world systems customers that are waiting for service may decide to abandon the system before entering service. In the process of designing systems, it is important to understand the effect of this abandonment phenomenon on the system's behavior. There has been a huge effort in developing models for systems that incorporate the effect of abandonments, also referred to as reneging or impatience (see, e.g., Dai et al. 2010;Garnett et al. 2002;Ward and Glynn 2005;Whitt 2004Whitt , 2006aKang and Ramanan 2010;Mandelbaum 2004, 2005). The simplest yet widely used model is the completely Markovian M/M/m + M model, also known as the Erlang A model. Its performance analysis has been an important subject of study in the literature (see for example Garnett et al. 2002;Whitt 2006c), not only because the Erlang A model is being used in practice (Mandelbaum and Zeltyn 2007), but also because it delivers valuable approximations for more general abandonment models (Whitt 2005).
The Erlang A model assumes Poisson arrivals with rate λ, exponential service times with mean 1/μ, m parallel servers, and most importantly, it incorporates the feature that waiting customers abandon the system after exponentially distributed times with mean 1/η. Let N (t) denote the queue length at time t. Assuming independence across the interarrival, service and reneging times, the queue length process is a birthdeath process (N (t)) t≥0 . The stationary distribution of this process, and associated performance measures like delay or abandonment probabilities, are easy to obtain (Garnett et al. 2002;Mandelbaum and Zeltyn 2007). In contrast, studying the timedependent behavior of (N (t)) t≥0 is generally judged to be prohibitively difficult (Fralix 2013;Ward 2012) because, among other things, the Kolmogorov forward equations do not seem to allow for a tractable solution. The main contributions of this paper are the exact solutions of both the forward and backward Kolmogorov equations, leading to exact expressions for the Laplace transforms of the time-dependent queue length distribution in Sect. 2 and first-passage times in Sect. 3.
The birth-death process describing the Erlang A model has birth rates, conditioned on N (t) = j, λ j = λ and death rates μ j = jμ for j ≤ m and μ j = mμ + ( j − m)η for j > m. There are available general results for the time-dependent behavior of birth-death processes. McGregor (1957a, b, 1958) have shown that the backward and forward Kolmogorov equations satisfied by the transition probabilities of a birth-death process can be solved via the introduction of a system of orthogonal polynomials and a spectral measure. For each set of birth and death rates (λ j , μ j ) there is an associated family of orthogonal polynomials. In some cases, when the set (λ j , μ j ) is assumed to have a special structure, these orthogonal polynomials can be identified. One such special case is the M/M/m queue, with λ j = λ and μ j = min{ j, m}μ. Notice that the Erlang A model incorporates the M/M/m queue as the special case η → 0 + . Karlin and McGregor (1958) have shown for the M/M/m queue that the relevant orthogonal polynomials are the Poisson-Charlier polynomials. Determining the spectral measure, though, is rather complicated, which is why van Doorn (1979) made a separate study of determining the spectral properties of the M/M/m queue, starting from the general expression for the spectral measure in Karlin and McGregor (1958) in terms of the Stieltjes transform. For the same M/M/m queue, Saaty (1960) derived the Laplace transform of Prob[N (t) = n] over time, in terms of hypergeometric functions. As in Saaty (1960), we do not resort to the approach in McGregor (1957a, b, 1958) for solving the Erlang A model, but instead opt to derive the explicit solution for the Laplace transform of Prob[N (t) = n] in a direct manner. The inverse transform then gives the desired solution for the time-dependent distribution, and we can also obtain the time-dependent moments. Mathematically, we shall use discrete Green's functions, contour integrals, and special functions related to hypergeometric functions. Having explicit expressions for the Laplace transforms is useful for ultimately obtaining various asymptotic formulas, which would likely be simpler than the full solution and yield insight into model behavior.
Due to the cumbersome expressions for some of the stationary characteristics, and the presumed intractability of the time-dependent distribution, simpler analytically tractable processes (D(t)) t≥0 have been constructed that have similar time-dependent and stationary behaviors as (N (t)) t≥0 . This can be done by imposing limiting regimes in which such approximating processes naturally arise as stochastic-process limits. Ward and Glynn (2002) make precise when the sample paths of the Erlang A model (and extensions using more general assumptions Ward and Glynn 2005) can be approximated by a diffusion process, where the type of diffusion process depends on the heavy-traffic regime. The diffusion process (D(t)) t≥0 is generally easier to study than the birth-death process (N (t)) t≥0 , and can thus be employed to obtain simple approximations for both the stationary and the time-dependent system behavior. In Ward and Glynn (2002, 2003 the limiting diffusion process is a reflected Ornstein-Uhlenbeck process, whose properties are well understood (Fricker et al. 1999;Linetsky 2005;Ward and Glynn 2003). Garnett et al. (2002) proved a diffusion limit for the Erlang A model in another heavy-traffic regime, known as the Halfin-Whitt or QED regime. In this regime, the diffusion process (D(t)) t≥0 is a combination of two Ornstein-Uhlenbeck processes with different restraining forces, depending on whether the process is below or above zero. Both the stationary behavior (Garnett et al. 2002) and the time-dependent behavior (Leeuwaarden and Knessl 2012) of this process are well understood. From our general result for the Laplace transform of Prob[N (t) = n] we show how the results obtained in Leeuwaarden and Knessl (2012) for the above diffusion processes can be recovered. See the survey paper Ward (2012) for a comprehensive overview of diffusion approximations for many-server systems with abandonments.
The paper is structured as follows. In Sect. 2 we work towards Theorems 1 and 2 that provide explicit expressions for the Laplace transform of the time-dependent distribution of (N (t)) t≥0 . In order to do so, we first reformulate the forward Kolmogorov equations in terms of difference equations in Sect. 2.1 and Laplace transforms in Sect. 2.2. In these steps we identify several key special functions of which some properties are proved in Sect. 2.3. Then, in Sect. 2.4 we prove Theorems 1 and 2 using all preliminary results in Sects. 2.1-2.3. In Sect. 2.5 we consider the limiting steady state case, and in Sect. 2.6 we treat the special cases η = μ = 1 (M/M/∞ queue), η → 0 + (M/M/m queue) and η → ∞ (the M/M/m/m loss model).
In Sect. 3 we follow a similar approach to obtain in Theorem 3 the Laplace transform of the distribution of the first time that (N (t)) t≥0 reaches some level n * > m. This result is established in Sect. 3.1. We then again specialize the general results in Theorem 3 to several limiting cases. In particular, we derive in Sect. 3.2 results for the Halfin-Whitt regime. Finally, we give in Sect. 3.3 an expression for the mean first passage time.

Transient distribution 2.1 Expressions in terms of difference equations
We let N (t) be the number of customers in the system and set (2.1) so that p n (t) depends parametrically on the initial condition n 0 , as well as the model parameters m, η and ρ = λ/μ. Since N (t) is a birth-death process with birth rate λ, death rate (setting μ = 1) N (t), for N (t) m, and death rate m and for n m + 1, with the initial condition p n (0) = δ(n, n 0 ), (2.6) with δ(n, n 0 ) = 1 for n = n 0 and δ(n, n 0 ) = 0 for n = n 0 . Setting and assuming that 0 < n 0 < m we obtain from (2.2) to (2.6) (2.10) If n 0 = 0 the right side of (2.8) must be replaced by −1, while if n 0 m the right side of (2.10) must be replaced by −δ(n, n 0 ), and then the right side of (2.9) is zero. We proceed to explicitly solve (2.8)-(2.10), distinguishing the cases 0 < n 0 < m and n 0 > m, and then we show that the results also apply for n 0 = m and n 0 = 0.

Special function solutions to the difference equations
Since the coefficients in the difference equations in (2.9) and (2.10) are linear functions of n, we can solve these explicitly with the help of contour integrals. First consider the integral where C 0 is a small circle in the z-plane, on which |z| < 1. The integrand in (2.11) is analytic inside the unit circle, if we define with | arg(1 − z)| < π, so that for z real and z < 1, arg(1 − z) = 0. By expanding as a binomial series, we obtain the alternate form where (·) is the Gamma function. It follows that F −1 (θ ) = 0, F 0 (θ ) = 1 and F 1 (θ ) = ρ + θ , and hence F n (θ ) satisfies Eq. (2.8). Furthermore, from (2.11) we have (2.14) as the contour C 0 is closed and the integrand in (2.14) is a perfect derivative. Thus F n (θ ) provides a solution to the homogeneous version of (2.9) (with the right side replaced by zero). We shall solve (2.8)-(2.10) using a discrete Green's function approach, and this will require a second, linearly independent, solution to (2.9). Such a solution may be obtained by using the same integrand as in (2.11) but integrating over a different contour. Thus we let where C 1 goes from −∞ − iε to −∞ + iε (ε > 0), encircling z = 1 in the counterclockwise sense (see Fig. 1). In (2.15) we use the branch (z −1) θ = |z −1| θ e iθ arg(z−1) , where | arg(z − 1)| < π, so the integrand is analytic in C − {Im(z) = 0, Re(z) 1}. By a calculation completely analogous to (2.14), and noting that C 1 begins and ends at z = −∞, where the integrand in (2.15) decays exponentially to zero, we see that G n (θ ) satisfies the homogeneous form of (2.9). However, G n does not satisfy the boundary equation in (2.8), and we now have (2.16) Now consider (2.10). We shall again construct two independent solutions to this difference equation. Let f n satisfy [ρ + θ + m + (n − m)η] f n = ρ f n−1 + [m + (n − m + 1)η] f n+1 and represent f n as a contour integral, with for some function F(·) and contour C. Then we have (2.18) We use integration by parts in (2.18) with and for now assume that C is such that there are no boundary contributions arising in (2.19), from endpoints of C. Using (2.19) in (2.18) we can rewrite (2.18) as a contour integral of z −n−1 times a function of z only, and if (2.18) is to hold for all n we argue that this function must vanish. We thus obtain the following differential equation for F(z): whose solution is, up to a multiplicative constant, Now we use (2.21) in (2.17) and make two different choices of C and different branches of (2.21), to obtain two independent solutions to (2.10). Note that now (2.21)

Fig. 1
A sketch of the branch cuts and the contours C 1 and C 2 has branch points both at z = 1 and z = 0. We let where C 1 is as in (2.15) (or Fig. 1) and the branch of (z − 1) θ/η is For the second solution we set where C 2 goes from −∞ − iε to −∞ + iε, encircling z = 0 in the counterclockwise sense, and (1−z) θ/η is defined to be analytic exterior to the branch cut where Im(z) = 0 and Re(z) 1, similarly to (2.11). Also, z m/η is defined as below (2.22), so the integrand in (2.24) is analytic exterior to the branch cuts where Im(z) = 0 and Re(z) 0 or Re(z) 1, and in particular on the contour C 2 (see again Fig. 1).
We have thus shown that the general solution to (2.10) is a linear combination of H n and I n , while that of (the homogeneous version of) (2.9) is a combination of F n and G n .

Properties of the key special functions
We now establish several useful properties of these functions. The integrals in (2.11), (2.15), (2.22) and (2.24) may all be expressed in terms of generalized hypergeometric functions, but we shall not use this fact. First, we note that if η = 1 then F n = I n and G n = H n . The latter is obvious since (2.15) and (2.22) have the same contour C 1 , while if η = 1 in (2.24) the branch point at z = 0 disappears and C 2 may be deformed to the loop C 0 in (2.11).
The functions H n and I n have very different asymptotic behaviors as n → ∞. For n large, standard singularity analysis shows that the asymptotics of I n are governed by the singularity at z = 1, and then setting z = 1 − ξ/n and letting n → ∞ in (2.24) yields Here C ξ goes from −∞ − iε to −∞ + iε, with ε > 0, and encircles ξ = 0. Thus I n has an algebraic dependence on n for n large. To expand H n we can simply dilate the contour C 1 in (2.22) to the range |z| 1 and then expand ( and hence H n decays roughly as 1/n!. Here we also used (n + x) ∼ (n)n x , which holds for n → ∞ and x fixed. Next we consider the discrete Wronskian Using the fact that H n and I n satisfy (2.10) we find that and thus W n [m + (n − m + 1)η] = ρW n−1 . Solving this simple recurrence leads to (2.29) To determine ω * (·) we let n → ∞ in (2.27), and use (2.25) and (2.26). Then H n+1 H n with I n+1 = O(I n ), so that Comparing (2.29) with (2.30) we determine ω * (·), and thus (2.31) A completely analogous calculation shows that which also follows by setting η = 1 in (2.31).

Main results: solutions of the difference equations
We solve (2.8)-(2.10) for 0 < n 0 < m, writing the solution as (2.33) Then P n will decay faster than exponentially as n → ∞, in view of (2.26), and satisfy the boundary equation in (2.8), since F n does but G n does not. It remains only to determine A, B, C, D; these functions will depend only on θ and the parameters m, η, ρ. By continuity at n = n 0 and n = m, we have (2.34) Then using (2.10) with n = m leads to and using (2.9) with n = n 0 (then δ(n, n 0 ) = 1) leads to If we introduce a and α by setting then both equations in (2.34) are satisfied. Using the fact that Then using the Wronskian identity in (2.32), with n = n 0 , and (2.38) we see that (2.42) Using (2.41) and (2.42) in (2.37)-(2.39), and then in (2.32) we have thus solved for P n (θ ), which we summarize below.
Theorem 1 For initial conditions 0 n 0 m, the Laplace transform P n (θ ) = ∞ 0 e −θt p n (t) dt of the time dependent distribution of N (t) is given by , n m; (2.43) (2.44) (2.45) Here F n G n , and H n are given by the contour integrals in (2.13), (2.15) and (2.22).
We note that if n 0 = m, (2.43)-(2.45) somewhat simplify, to Next we assume that N (0) = n 0 > m. Now we must solve the homogeneous form of (2.9) with the boundary condition in (2.8), and these imply that P n (θ ) must be proportional to F n for all 0 n m. Thus now G n will not enter the analysis. For n large, P n (θ ) must again be proportional to H n , which has the appropriate decay as n → ∞. Now we write Theorem 2 For initial conditions n 0 m, P n (θ ) is given by (2.55) Here F n , H n and I n are given by the contour integrals in (2.11), (2.22) and (2.24).
We note that if n 0 = m, expression (2.54) is not needed, and then (2.53) and (2.55) agree with the expression(s) in (2.47). Setting n 0 = m in (2.53) and using which follows from (2.31), we obtain (2.47) for n m.

Limiting case: steady state
We proceed to examine some limiting cases of Theorems 1 and 2, where the expressions simplify, sometimes considerably. First we consider the steady state limit of p n (t) as t → ∞, which corresponds to the limit of θ P n (θ ) as θ → 0. First, observe that as θ → 0, (2.11) and (2.15) yields (2.61) We can take |z| > 1 on C 1 , and then expand (z − 1) −1 as a Laurent series on C 1 . Using (2.58) and (2.11) we have where now on C 0 we can expand (1−z) −1 as ∞ =0 z , since |z| < 1. Combining (2.61) with (2.62), (2.60) then yields . Since I m = H m when θ = 0, the difference between 1 and 2 is  ). We have thus obtained the steady state limit from Theorem 2 as stated below (see e.g., Garnett et al. 2002;Ward 2012).

Corollary 1 The steady state distribution is
(2.69) Note that K and 1 are related by ρ 1 (1 + m/η)(ρ/η) −m/η K = 1. While we obtained Corollary 1 from Theorem 2, which applies for n 0 m, the result is independent of n 0 and Corollary 1 will also follow from Theorem 1 using very similar calculations to those in (2.60)-(2.66), which we omit. Of course, p n (∞) is more easily obtained by letting t → ∞ in (2.2)-(2.5) and solving the resulting elementary difference equations.

Limiting cases: extreme abandonments rates
Next we evaluate Theorems 1 and 2 for the special cases η = 1, η → 0 + (vanishing abandonment effects) and η → ∞. For η = 1 the model reduces to the standard infinite server M/M/∞ queue, and from Theorems 1 and 2 we obtain the following.

Corollary 2 When η = 1 the Laplace transform of p n (t) is given by
(2.70)

72)
and an alternate form is given by and then p n (∞) = e −ρ ρ n /n!.
We have already seen that when η = 1, H n = G n and I n = F n and we have the Wronskian identities in (2.31) and (2.32). Then both Theorems 1 and 2 reduce to (2.70), and we need not distinguish the cases n 0 ≷ m, as m disappears altogether from the expressions. Now, from (2.11) and (2.15) it is clear that F n (θ ) and G n (θ ) are entire functions of θ , for every n. Thus the only singularities of (2.70) are the poles of (θ), which occur at θ = −N , N = 0, 1, 2, . . . and the corresponding residues are (−1) N /N !. When θ = −N , G n and F n are no longer linearly independent, and in fact G n (−N ) = (−1) N F n (−N ), which follows by comparing (2.15) with (2.11). Thus evaluating the contour integral p n (t) = (2πi) −1 Br e θt P n (θ ) dθ (where Re(θ ) > 0 on the vertical Bromwich contour) as a residue series we obtain precisely (2.71), with (2.72). To obtain the expression in (2.73) we represent the F n (−k) in (2.71) as contour integrals, yielding Then expanding we −t + 1 − e −t n 0 using the binomial theorem leads to (2.73).
Note that as t → 0 1 − e −t n+n 0 −2 j → 0 unless n = n 0 and j = n, so that p n (0) = δ(n, n 0 ). As t → ∞ only the term with j = 0 in (2.73) remains and we obtain the steady state Poisson distribution. If η = 1 it is easier to solve (2.2)-(2.6) using the generating function G(t, u) = ∞ n=0 p n (t)u n which leads to the first order PDE Inverting the generating function then regains (2.73). Next we let η → 0 + , so that the model reduces to the m-server M/M/m queue. Then we obtain the following.
Corollary 3 When η = 0 the Laplace transform of p n (t) is given by, for 0 n 0 m, , n m, (2.77) (2.80) We also note that the transient distribution for the M/M/m model was previously obtained, in different forms, by Saaty (1960) andvan Doorn (1979). In van Doorn (1979) spectral methods are used, while in Saaty (1960) the Laplace transform is expressed in terms of hypergeometric functions.
To establish (2.77)-(2.84) we need to evaluate H n (θ ; m) and I n (θ ; m) for η → 0 + . We write H n in (2.22) as so that the integrand has saddle points where ∂ f /∂z = 0, and this occurs at (2.87) We can take |z| > 1 on C 1 and then the saddle at Z + determines the asymptotic behavior of H n as It follows that H n−1 /H n ∼ Z + in this limit, as A = 1/Z + and Z ± satisfy the quadratic equation Hence (2.43) reduces to (2.77) as η → 0 + . Also, (2.79) and (2.80) follow from (2.44) and (2.45), in view of (2.89) and the fact that Now consider n 0 m. We shall obtain (2.81)-(2.84) from Theorem 2. We must then expand I n (θ ; m) for η → 0 + . Using the saddle point method we find that as the expansion of (2.24), which involves the contour C 2 , is determined by the other saddle point in (2.87). Using (2.88) with n replaced by n 0 , (2.91), and Stirling's formula we obtain, after a lengthy calculation, the following limit (as η → 0 + ): After factoring out I n , the bracketed factor in (2.54) becomes and we also use (2.90) with n replaced by n 0 , and Then (2.55) goes to the limit in (2.81), since ρ Z + /m = B. This completes the proof of Corollary 3.
In the limit η → ∞, we expect our results to reduce to the Erlang loss model, or the M/M/m/m queue. We then obtain the following.

Corollary 4
As η → ∞ the Laplace transform of p n (t), for 0 n 0 m, approaches the limit P n (θ ) = n 0 ! (θ)e −ρ ρ n 0 +θ F n 0 [G n + ωF n ], n 0 n m (2.94) In particular the blocking probability p m (t) has the Laplace transform . (2.96) Note that (2.35) follows by setting n = m in (2.94) and using the Wronskian W m in (2.32). To establish Corollary 4, we note that by expanding the integrand in (2.22) for η → ∞ we obtain so with (2.98), (2.47) and (2.45) yields (2.94) and(2.95) in the limit η → ∞. The blocking probability in (2.96) may also be written as , which follows from (2.11) with n = m and an integration by parts. If N (0) = 0 (starting with an empty system) we obtain from (2.99) (2.100) Previously expressions for the Laplace transform of the blocking probability were obtained by Jagerman (2000), who showed that (if n 0 = 0) (2.101) and this can easily be shown to agree with both (2.100) and (2.99), as follows by expanding both integrands using the binomial theorem. For general n 0 ∈ [0, m] the blocking probability is given by Since the expressions in Theorems 1 and 2 and even Corollaries 3 and 4, are quite complicated, it is useful to expand these in various asymptotic limits. One such limit would have m → ∞, ρ → ∞ with m/ρ → 1 and m − ρ = O( √ m). This is a diffusion limit, sometimes referred to as the Halfin-Whitt regime. Here we would scale n, n 0 and ρ, for m → ∞, as and x, x 0 and β are O(1). In this limit we can approximate the contour integrals F n , G n , H n and I n by simpler special functions, namely parabolic cylinder functions. We discuss this limit in detail in van Leeuwaarden and Knessl (2011) for the M/M/m model with η = 0, and in Leeuwaarden and Knessl (2012) for the M/M/m + M model with η > 0. We can obtain then p n (t) ∼ m −1/2 P(x, t) where P will satisfy a parabolic PDE, which we explicitly solved in van Leeuwaarden and Knessl (2011), Leeuwaarden andKnessl (2012). An alternate approach is to evaluate Theorems 1 and 2, or Corollary 3 in the limit in (2.104), and thus identify P(x, t) directly. We shall discuss in more detail the limit in (2.104) for the first passage distributions. We also comment that the transient behavior of the M/M/m/m model was analyzed thoroughly in Knessl (1990) and Xie and Knessl (1993), for m → ∞ and various cases of ρ, including the scaling in (2.104). There we used mostly singular perturbation methods, but equivalent results could be obtained using Corollary 4 and methods for asymptotically expanding integrals.

Main result: Laplace transform of the first passage time
Here we compute the distribution of the time for the number N (t) of customers to reach some level n * , which may be viewed as a measure of congestion. We take n * > m, for otherwise the problem reduces to that of the M/M/∞ or M/M/m/m models. Thus we define the stopping time τ (n * ) = min{t : N (t) = n * }, (3.1) and its conditional distribution is (3.2) When n = n * we clearly have and for n < n * , Q n (t) satisfies the backward Kolmogorov equation(s) and, expecting that Q n (0) = 0 for n < n * , we obtain (3.10) The recurrences in (3.9) and (3.10) are similar to those in (2.9) and (2.10), and indeed we can convert the former to the latter by setting Then from (3.7) and (3.12) we have (3.13) and R n (θ ) will satisfy (ρ + n + θ)R n = (n + 1)R n+1 + ρ R n−1 for 0 < n < m, which is just the homogeneous version of (2.9), while for n > m, R n (θ ) will satisfy (2.10). Also, R 1 (θ ) = (ρ + θ)R 0 (θ ), so that R n (θ ) will satisfy the boundary equation in (2.8). We can thus write R n in terms of the special functions F n , G n , H n , I n that we introduced in Sect. 2, and since F n satisfies (2.8) we write and if both (3.14) and (3.15) apply for n = m we have the continuity equation (3.17) Finally, using (3.5) with n = m and noting that, in view of (3.11) and (3.12), we find that (m + η)R m+1 + ρ R m−1 = (θ + ρ + m)R m and thus Then (3.16), (3.17) and (3.19) yield three equations for the unknowns c 1 , c 2 , c 3 . After some algebra and use of (2.31) with n = m we obtain R n , and then Q n follows from (3.11) and (3.12). We summarize below the final results.

Theorem 3
The distribution of the first passage time to a level n * (> m) has the Laplace transform Q n (θ ) = E e −θτ (n * ) | N (0) = n : Note that actually (3.20) can be used even if n = m + 1 and it then agrees with (3.21). Similarly, (3.21) holds even if n = m − 1. If η = 1 we have F n = I n and then both (3.20) and (3.21) reduce to which is the result for the M/M/∞ model. We can again get results for the standard M/M/m model by letting η → 0 + in Theorem 3. Using the asymptotic results in (2.88) and (2.91), after some calculations that we omit we obtain the following.

Corollary 5 For the M/M/m model the first passage distribution to a level n * (> m) is given by
m n n * . (3.24) Here Z ± are as in (2.87).
Using the fact that F n (0) = ρ n /n! and Z ± (0) = [m + ρ ± |m − ρ|]/(2ρ) we can easily verify that Q n (0) = 1 for all n, so that the density is properly normalized. We shall discuss later the mean first passage time, which is equal to − Q n (0).

Halfin-Whitt regime
We next consider the limit in (2.104) in Corollary 5, also scaling the exit point n * as From (2.87) we obtain, using (2.104), (3.26) By scaling z = 1 − ξ/ √ m in (2.11) and noting that ρz − n log z = ρ + (x + β)ξ + 1 2 ξ 2 + o(1) with the Halfin-Whitt scaling in (2.104), the integral in (2.11) can be approximated by where D p (z) is the parabolic cylinder function of index p and argument z. In (3.27) the approximating contour Br + is a vertical contour in the ξ -plane, on which Re(ξ ) > 0, and ξ −θ is defined to be analytic for Re(ξ ) > 0 and real and positive for ξ real and positive. In view of (2.104), setting n = m corresponds to x = 0 and thus (3.28) A similar calculation shows that and we note that the difference F m+1 − F m is smaller than F m by a factor of m −1/2 .

Corollary 6
In the limit m → ∞, with the scaling in (2.104) and (3.25), the transform of the first passage distribution Q n (θ ) for the M/M/m model has the limit P(x, θ) where and We have previously obtained these results in Fralix et al. (2014), by directly solving the parabolic PDE satisfied by the diffusion approximation. Since 2D 1−θ (−β) + β D −θ (−β) = −2D −θ (−β), Corollary 6 agrees with Theorems 1 and 2 in Fralix et al. (2014). Now, we can also consider the Halfin-Whitt limit for the first passage distribution in the M/M/m + M model (with a fixed η > 0), and then Theorem 3 reduces to the following.