Exact Results for the Distribution of the Partial Busy Period for a Multi-Server Queue

Exact explicit results are derived for the distribution of the partial busy period of the M/M/c multi-server queue for a general number of servers. A rudimentary spectral method leads to a representation that is amenable to efficient numerical computation across the entire ergodic region. An alternative algebraic approach yields a representation as a finite sum of Marcum Q-functions depending on the roots of certain polynomials that are explicitly determined for an arbitrary number of servers. Asymptotic forms are derived in the limit of a large number of servers under two scaling regimes, and also for the large-time limit. Connections are made with previous work. The present work is the first to offer tangible exact results for the distribution when the number of servers is greater than two.

1. Introduction One of the fundamental concepts in queueing theory is the busy period (BP).Other basic quantities include the stationary queue-size and waiting-time distributions.For the archetypal M/M/c multi-server queueing model, the latter are extensively covered in all introductory textbooks, e.g.[12,28].But they are relevant only in the ergodic region where the total traffic intensity is less than unity.By contrast, the BP persists as a valid concept beyond this region.However, there is little available to be said about its calculation, even in the ergodic case.
There are two main characterizations of a BP.The full BP refers to the time interval during which all servers are occupied.The partial BP refers to the time interval during which at least one server is occupied.For a single-server queue, the two definitions trivially coincide.It is the distribution of the partial BP that will be studied here, and use of the term 'busy period' on its own will assume the partial variant.
We shall confine our attention to the M/M/c queue, and denote the number of servers by N .Thus c = N in our discussion.Within this class of models, the BP distribution for the single-server case is known.For a multi-server queue, the full BP is simply related to that of the single-server case with an adjusted service time [9].On the other hand, results on the partial BP for a multi-server queue are sparse.
For more than two servers, there exists in the published literature neither numerical algorithms nor analytic expressions for the distribution of the partial BP.The present work provides both of these.We first demonstrate a spectral method that serves as the basis for a simple and efficient numerical algorithm that covers at least the entire ergodic region and handles server numbers up to around N = 80.We have tested it on the interval of traffic intensities 0 ≤ r ≤ 1.We also develop an algebraic method that gives rise to an explicit representation of the BP distribution for any number of servers as a finite sum of Marcum Q-functions, dependent on the roots of a certain family of polynomials whose coefficients we determine.Karlin and McGregor [16] give an explicit integral representation for the case of two servers based on a spectral method involving families of orthogonal polynomials.But, while they claim that their method can in principle be extended to a larger number of servers, this must be done on a case by case basis, and the calculations become so unwieldy beyond N = 2 that this line of attack has not hitherto been pursued by anyone.By contrast, we show that a much more rudimentary spectral method yields more generally applicable results with greater ease.Arora [3] has also derived the distribution for the two-server problem as a infinite series of modified Bessel functions.Our algebraic approach can be viewed as a generalization of this to an arbitrary number of servers.
All other work appears to be centered around computing moments of the BP distribution.Natvig [21] obtains the first and second order moments of the length of the partial BP.Second moments of the partial BP distribution are also derived by Omahen and Marathe [22].Further progress has been made more recently by Artalejo and López-Herrero [4] who show how to compute arbitrary moments.As commented upon by these authors, prior to their paper in 2001, general results were only known for the one and two server cases, and moments up to second order had appeared for the general problem.We are not aware of anything more recent that alters this state of affairs.

Motivation
The BP distribution for the basic M/M/c queue is equally applicable to a large family of more complex, albeit memoryless, queueing systems that incorporate multiple classes of arrivals each comprised of multiple priority levels that are processed under a variety of queuing disciplines (e.g.preemptive or non-preemptive).A concrete example is furnished by the ambulance ramping problem encountered by hospital emergency departments (ED), where one question that may be studied is to what extent an ambulance offload zone alleviates the ambulance queue that builds up while patients wait to be admitted and prevents ambulances from being dispatched to other calls [18,15].Patients are differentiated via arrival source into the ED as either walk-ins and ambulance arrivals, with ambulance arrivals further sub-divided into those who enter the ED from the ambulance queue or the offload zone.Multiple priorities are assigned to patients in each arrival class depending on their acuity levels, and each arrival-class/priority-level combination can exhibit a different arrival rate.The main constraint is that of a common mean service time (i.e.treatment time) for all patients.It does not matter to the servers how the patients are arranged prior to entry, or on the mixture currently in service.Hence the BP distribution is just that of the basic model.Rastpour [25] has recently studied partial BP distributions in the context of emergency medical services.
Our interest in the partial BP emerged from its relevance to regenerative simulation of the steady-state limit of multi-server, multi-priority, multi-class queues as discussed above.The regenerative simulation technique was introduced by Crane and Iglehart [8].Empirical distributions for queue lengths and waiting times per class or priority level can be ascertained from steady-state discrete event simulations, and subsequently compared with theoretical predictions.In the ergodic region, sample means for summary statistics are easily estimated.However, confidence intervals are also required to judge whether the null hypothesis that the empirical and theoretical values are generated by the same candidate distribution is to be rejected at some given level of statistical significance.In steady state simulation, since one is producing a very long single run of a stochastic (e.g.Markov) process, the data underlying the empirical distributions are highly correlated.Difficulties with estimating confidence intervals arising from unwanted correlation can be avoided by recognizing that one is dealing with a regenerative process, and partitioning the time series into consecutive regeneration cycles.These are statistically independent and identically distributed as each cycle effectively restarts the process without memory of the past.The natural regeneration point in a queuing problem is the empty state, which occurs an infinite number of times in an ergodic system.Thus, each regeneration cycle comprises an initial idle (i.e.empty) period followed by a partial BP.The distribution of the idle period is simply the inter-arrival distribution, assumed here to be exponential.Therefore, the interesting quantity in the regeneration cycle problem is the partial BP.In carrying out a state-state simulation of a queueing system, it is a straightforward matter to collect data on the lengths of regeneration cycles and BPs.Comparison of these empirical results with theoretical expectations serves as a useful diagnostic tool for judging the soundness of the simulation.It was the paucity of theoretical results for this comparison that led to the present investigation.
The remainder of the paper is organized as follows.Our starting point in Section 3 will be the recursive system of equations considered by Artalejo and López-Herrero [4], which can be solved to give the moment generating function (MGF) of the partial BP.We begin by recasting this as a matrix inversion problem for a tridiagonal matrix.The matrix, which has the form of a generator for a birth-death process, is inverted by means of its spectral decomposition.This leads to an explicit representation of the distribution in terms of the eigenstates that is amenable to efficient numerical implementation.This is followed in Section 4 by a algebraic approach where the recurrence relations of Artalejo and López-Herrero [4] are solved directly for the MGF.The method leads to a finite continued fraction representation similar to previous findings in the literature, but which does not appear to be suitable for further analysis.An alternative method yields an explicit form for the MGF in terms of a family of polynomials, whose degree increases with the number of servers, and which depend parametrically on the traffic intensity.A simple two-dimensional system of recurrence relations allows the coefficients of these polynomials to be computed.In Section 5, asymptotic limits as the number of servers N approaches infinity are derived under two different scaling regimes.These are subsequently used for comparison in validating exact results for large but finite N , and determining how large N must be for one to judge that asymptotic behaviour has effectively set in.In Section 6, the polynomial-based representation of the MGF is combined with a complex contour integral that implements the inverse Laplace transform, whose structure is determined by poles that arise from zeros of the polynomials.The contour integration evaluates to explicit closed-form expressions for the partial BP distribution as a finite sum of Marcum Qfunctions.In Section 7, various summary statistics are computed from the posited theoretical distributions and compared, across a large range of model parameters, with known exact results and with the output of Monte Carlo simulation.Excellent agreement is demonstrated over the entire ergodic range of traffic intensity, and for server numbers up to well within the asymptotic regime.Concluding remarks are presented in Section 8.

Spectral Method
We consider the M/M/c queue with c = N servers, Poisson arrival rate λ and mean treatment (or service) time 1/µ.Following [4], the quantities ϕ n (s), n = 1, 2, . .., where n labels the number of busy servers, are defined to be MGFs of the first passage times from the state n = 1, 2, . . ., N to the empty state 0. Also, let T bp be the random variable (RV) describing the partial BP, in which case its MGF is given by ϕ 1 (s) ≡ ⟨e −sT bp ⟩ T bp .Then, if we set where r ≡ λ/(N µ) is the total traffic intensity, the MGF for the single-server problem (N = 1) reads ϕ 1 (s) = ψ − (s).We confine our attention to the ergodic region 0 ≤ r < 1, but we shall also study the boundary value r = 1 as a special case.We also introduce µ n ≡ nµ, for n = 1, 2, . . ., N , and without loss of generality, we may set µ = 1/N , so that r = λ and Finally, we write µ n ≡ nµ 1 with µ 1 = 1/N , so that 0 ≤ µ n ≤ 1 for all n.
The linear recurrence relations of Artalejo and López-Herrero [4] may be cast as follows: For and ϕ N = ψ − (s)ϕ N −1 , with ϕ 0 ≡ 1.As the original paper does not make the derivation of this result explicit, we supply the details in the Appendix.These equations have an (N − 1) × (N − 1) tridiagonal matrix structure that may be expressed as In Dirac notation [31], this matrix equation reads where the matrix M is the birth-death process generator It may be compared with the matrices introduced in (1.3) and (4.1) of Karlin and McGregor [16].Its trace and determinant, useful for checking the numerical solution of the eigenvalue problem, are If we introduce the standard orthonormal Cartesian basis on R N −1 as |e i ⟩, i = 1, 2, . . ., N − 1, then we have that Since M is non-symmetric, it has distinct right and left eigenvectors so that its eigenvalue problem is written as With the eigenvectors required to satisfy the orthonormality conditions ⟨v To determine ϕ 1 (s), we invert (5) to obtain and observe that it suffices to consider just the first and last elements of this equation, namely The matrix M may be symmetrized by means of a diagonal similarity transformation and the eigenvectors (assumed orthonormal) are related by Let us define It follows that where It is useful to note that the following consequence of the completeness relation for the eigenvector basis for m, n = 1, 2, where δ mn is the Kronecker delta.We can now express (12) as where we have introduced the resolvent functions We can eliminate ϕ N by applying the relation ϕ N = ψ − (s)ϕ N −1 to the second equation in (19), which yields upon noting that ψ + (s)ψ − (s) = 1/r.This expression is then inserted into the first equation in (19) to provide the explicit result for the MGF It should be noted that, due to cancellations, there are actually no singularities at the eigenvalues s = ξ k .The singularity structure in the complex s-plane of ϕ 1 (s) comprises a cut and potentially a finite number of poles.The cut is due to the presence of ψ + (s), and lies along the negative real axis corresponding to the finite interval for which the discriminant ∆(s) ≡ b 2 (s) − 4r is negative, and where b(s) ≡ r + 1 + s.Thus, the cut spans x − ≤ −s ≤ x + , where the cut limits are given by The poles arise from the vanishing of the denominator in (22), and only lie on the negative real axis between the origin and the near cut boundary, as discussed in [16].This is illustrated in Figure 1 for the case r = 0.5, N = 30, the details of which will be explained later on.
The BP distribution is recovered from the MGF ϕ 1 (s) by an inverse Laplace transform implemented as a Bromwich contour integral in the complex s-plane.The contour may be deformed onto the negative real axis, giving rise to discrete residue contributions from any relevant poles, plus a real-valued integral on the cut.An explicit representation of this integral has been given for the two-server case in equation (6.25) of [16].It is also of the same type as the integral arising in the waiting-time distribution for priority queues, as can be seen in [10].One should observe that, in the present formalism, the two-server problem is trivial as the matrix M is one dimensional, so that there is no eigenvalue problem to solve.
To evaluate the contribution to the BP distribution from the cut, we introduce the product function whose square will multiply both the numerator and denominator in (22) in order to eliminate spurious singularities at the eigenvalues.This results in the appearance of the non-singular functions Next, letting D ± (s) ≡ ψ ± (s) − G 22 (s) so that D + (s) is the denominator in (22), we see that By appealing to the identity 1/D + (s) = D − (s)/[D + (s)D − (s)], we are able to construct the cut contribution as where the cut function R cut (s) is given by This definition takes account of the fact that the cut is traversed in two opposite directions.For N = 1, the cut function is just the constant R cut (s) = 1/r while, for N = 2, we obtain Now, since Π(ξ k ) = 0 for all k = 1, 2, . . ., N − 1, it follows that But we also have that Thus, which provides a useful numerical check.By considering the point s = 0, we obtain the identities which imply that The eigenvalue problem for the symmetric matrix S is solved numerically using the LAPACK routine dstevd, that implements a divide-and-conquer method [2].Identical results are produced when using Matlab's eig function, which is probably an indication that eig invokes the same algorithm when applied to this problem.On making the change of integration variable x = −s to obtain positive values for the coordinates of the cut, we can express the cut contribution to the BP probability density function (PDF) as For the survival function1 (SF), r denote the length of the cut.Then a further change of variable, such that x = x − + ∆x•u, puts the integral over the unit interval where A further change of variable v = 2u − 1 casts the integral into a form that is amenable to efficient Gauss-Chebyshev quadrature: It is known that Gauss-Chebyshev quadrature of the second kind is equivalent to the trapezoidal rule on the unit circle [6].Thus, if we set v = cos(πτ ), so that 0 ≤ τ ≤ 1, and subsequently uniformly discretize the unit interval into L sub-intervals according to τ n = n/L, then we obtain the L-point quadrature rule for the SF given by with nodes and weights, respectively, Quadrature for the PDF is given by Equation ( 38) may also be expressed in a form suitable for Gauss-Chebyshev quadrature of the first kind, namely, It is known that Gauss-Chebyshev quadrature of the first kind is equivalent to the mid-point rule on the unit circle, which gives rise to the same quadrature scheme as described for the trapezoidal rule, but where the nodes are generated by τ n = (n − 1/2)/L for n = 1, 2, . . ., L, so that 0 < τ n < 1 which avoids any potential singularities.We find that the most numerically robust scheme over the entire range of traffic intensities r corresponds to the mid-point rule with the quadrature weights expressed in the form The trapezoidal rule is generally a poorly performing quadrature scheme.However, it has been shown to be exponentially convergent for some specific classes of integrand, such as those comprising smooth periodic functions in which the integration extends over a period [26,29].It should be noted that the integrand above is unit-periodic in τ .Because of this, and the fact that it is equivalent to Gauss-Chebyshev quadrature, the trapezoidal rule is an efficient way to compute the BP distribution.This is also true for the mid-point rule.
Alternatively, in order to attain a guaranteed level of accuracy, we may apply an iterative trapezoidal or mid-point rule to the integral  with and Gaussian quadrature schemes do not generally allow nodes and weights from a previous iteration to be re-used when sequentially moving to a finer grid in order to improve accuracy on the way to satisfying some predetermined convergence criterion.However, the trapezoidal rule can be iterated by doubling the number of grid points at each step.Only half the grid points need to be computed, as the grid points from the previous iteration supply the other half [24].Similar considerations apply to the iterative mid-point rule, but the grid size needs to be tripled at each step, with a third of the total number of points coming from the previous iteration.The large-t asymptotics are derived from (36) by We shall now look at the pole contribution to the BP distribution.Poles occur in (22) when ψ + (s) − G 22 (s) = 0. Let us suppose that these are found at s = ζ ℓ , ℓ = 1, 2, . ... Since the first term G 11 (s) will not contribute to the residue, the pole components of the MGF are given by the residues where the prime denotes differentiation.Since, for general s, ψ ′ + (s) = ψ + (s)/ ∆(s), we have at a pole Then, on setting we obtain Therefore, the general result for the pole contribution to the BP SF is with and The number of ζ-poles can be determined a priori according to The notation I[A] is for the indicator function such that I[A] = 1 if A is true, and is zero otherwise.The initial sum counts the contributing ξ-asymptotes, while the next two terms account for boundary effects.The maximum number of poles that can contribute is N − 1, and this occurs when r → 0 + , in which case x ± = 1, so that the cut disappears.Figure 1 displays the generation of the ζ-poles as the intersection of a sideways half-parabola defined by ψ + (s) with a tan-like function generated by the ξ-poles of G 22 (s).The tip of the parabola lies on the near edge of the cut.Thus, the parabola can only intercept the ξ-pole curves at points between the origin and the near cut boundary.The full BP distribution is given by the sum of the pole and cut contributions, i.e.F (t) = Fpol (t) + Fcut (t) in the case of the SF.By differentiating ϕ(s) at s = 0, we can obtain an expression for the mean BP given by This provides a useful numerical diagnostic when compared with the exact result [4] ⟨T bp ⟩ T bp = 1 r Here, Γ(n, x) denotes the upper incomplete gamma function with shape parameter n.
The upshot of the foregoing analysis is that the BP distribution is arbitrarily well represented as an exponential mixture.Thus, suppose that the RV X is an exponential mixture, so that for the SF.Then, we have for the mean For the log-mean, where γ e ≡ −ψ(1) ≃ 0.5772 denotes Euler's constant.For the differential entropy, which is easily computed via Gauss-Laguerre quadrature.In the present application to the BP, the nodes should account for the discrete poles as well as the cut.Thus, L = L cut + L pol .The discrete poles contribute nodes From the requirement that FX (0) = 1, we have that L ℓ=1 w ℓ = 1.Thus, the node/weight pairs (u ℓ , w ℓ ) may be thought of as describing the empirical SF for some RV 0 ≤ U < ∞ according to which we may refer to as the texture distribution by analogy with compound Gaussian clutter distributions in radar detection theory [27].This compound representation implies that the BP RV has the product form T bp = U V , where here V is a unit-mean exponentially distributed RV, which provides a very convenient scheme for generating random numbers drawn from the BP distribution.Figures 2 and 3 display the logarithm of the BP SF for the case of N = 30 servers and a variety of traffic intensities.In Figure 3, the horizontal axis records the (dimensionless) time as directly given by F (t) = Fpol (t) + Fcut (t), in which case one is observing the short-time behaviour of the distribution.Note that the MGF for the problem formulated in terms of physical units can be recovered from the dimensionless problem via ϕ phys (s) = ϕ 1 (s/s scl ) with frequency scale s scl = N µ.For the SF, this implies F (t) = Fphys (t scl •t) with timescale t scl = 1/s scl = 1/(N µ).In Figure 2, the horizontal axis measures time in units of the mean BP.Thus, the function F (m bp t) is being plotted, and this illustrates the quite distinct long-time behaviour of the distribution.
Since the regeneration cycle time is the sum of a BP and an idle period, and these are mutually independent, it follows that the regeneration cycle distribution is the convolution of the BP distribution and the inter-arrival distribution.Therefore, the SF for the regeneration cycle time can formally be expressed as This expression should be numerically unproblematic for r < 1/4.Above this value, the inter-arrival pole at s = −r lies within the cut {s : It should be noted, however, that for appreciable values of the traffic intensity, the idle period gives a negligible contribution to the regeneration cycle.3.1.Boundary Case: r → 0 Numerical instabilities arise for small values of the total traffic intensity r.However, in the r → 0 + limit, the v 12 become negligible, in which case we have It follows, for the BP PDF and SD, respectively, that recalling that ξ k < 0. A criterion can be developed to signal when to use this approximation.We have found that the test max |v

Boundary Case
When the smallest pole is very close to zero, numerical solution of this equation can present difficulties.In this case, we can linearize about s = 0 to obtain which leads to the approximate solution .
In order to catch bad numerical behaviour due to arithmetic underflow, when the condition is performed.If this check succeeds, then ζ 1 is replaced by the value obtained from the latter of (66).On the other hand, failure of the test indicates that the lowest pole has missed altogether, in which case ζ 1 obtained from (66) is appended to the collection of poles ζ k .A variety of diagnostic tests is available to track the validity of the numerical computations.For example, it can be checked (i) whether the expected number of ζ-poles is being generated, (ii) whether equation (56) for the mean holds to some desired accuracy, (iii) whether the relationship (33) holds to some desired accuracy and, among others, (iv) whether the MGF as given by (22) evaluates to unity at the origin.

Boundary
Case: r = 1 At the boundary of the ergodic region, where r = 1, there is no pole contribution.Hence, the BP SF comes entirely from the cut, and reads where we may recall that, for r = 1, We see that the BP distribution becomes power-law tailed in the non-ergodic region.We can write to obtain a form that is unity at the origin, rather than diverging.One can also adopt the Lomax form F (t) For N = 1, 2, the Lomax form gives a close fit across the whole domain.

Algebraic Approach
In this section, we exploit the tridiagonal structure of the problem to derive a continued-fraction representation of the BP MGF that can be used to yield explicit expressions on a case-by-case basis for small numbers of servers.We then proceed to analyse the underlying tridiagonal matrix in a different way that allows explicit representation of the cut function for any number of servers in terms of a family of polynomials (to which we refer as the cut polynomials).Finally, we derive an easily solvable two-dimensional system of recurrence relations that enables determination of the coefficients of the cut polynomials which, in turn, allows for the calculation of their roots.As demonstrated in a subsequent section, these roots provide complete information for the construction of closed-form expressions for the BP distribution functions.

Continued Fraction
where M is the full tridiagonal birth-death process generator matrix for the BP problem as given in (6), and let M 1 be the matrix obtained from M 0 by omitting the first row and column.We then define the general matrix M k recursively as the one obtained by omitting the first row and column from M k−1 .This can be continued until one reaches the scalar Cramer's rule yields By applying the cofactor expansion in the first row, we are led to the recurrence relation where we have introduced the notation on applying the recurrence once, and then successively applying the recurrence, we arrive at the continued fraction representation One may note that if we set η N (s) ≡ rψ − (s), and for k = N − 1, N − 2, . . ., 1, then we have ϕ 1 (s) = η 1 (s)/r.Setting b(s) ≡ s + r + 1, we can write it follows that rψ from which we also obtain the identity With the aid of these relationships, (76) leads directly to the following explicit forms linear in ψ − (s): For N = 3, with For N = 1, we trivially have ϕ 1 (s) = ψ − (s).While the continued-fraction representation can be used to derived explicit expressions for the MGF when the number of servers is small, there does not seem to be any obvious way to extract an explicit closed-form result for the general case.A finite continued fraction representation was previously considered by Daley and Servi [9], who make the same observation that 'there is no simple non-recursive closed-form solution'.Prior to this, Conolly [7] had also derived a continued fraction representation.

Cut Polynomials
Further progress for the general case can be made by observing the general structure where N 0,1 and D 0,1 are polynomials in s.This representation follows from applying the cofactor expansion for det M 0 and det M 1 to the last rather than the first row.At a pole, which occurs when σ N −1 − rψ − = −D 1 /D 0 , the numerator becomes rational: It is convenient to introduce the functions so that ϕ 1 (s) = µ 1 N − (s)/D − (s).One can show that which is a polynomial in s.We also have for some polynomial A(s), and one can show that Thus, we see that the representation has polynomial denominator and is linear in ψ − (s).
On the cut in the s-plane, that is due to square-root term in ψ − (s) = [b(s) − |∆(s)|]/(2r), the function ϕ 1 (s) contributes only the component where |∆(s)| = 4r − b 2 (s) on the cut.This allows us to identify the cut function introduced previously as where we find it convenient to introduce the cut polynomials C(σ), defined as functions of the variable σ ≡ σ 1 = s + r + µ 1 via C N (σ(s)) ≡ −µ 1 D + (s)D − (s).The first few of these are given as follows: For N = 2, C 2 (σ) = σ − r/µ 1 .For N = 3, For N = 4, In general, For consistency with (92), one may set C 1 (σ) ≡ −1/µ 2 1 for the single server case.The residue of ϕ 1 (s) at a given pole s, defined by where the prime in (D + D − ) ′ denotes that the pole factor has been divided out of the polynomial, i.e.
We can derive the general relationship and note that D − = 0 at a pole, to obtain that, at a pole, Consequently, at a pole s, the residue term is It may be observed that the polynomials N 0 , N 1 have completely dropped out of the calculations.
The full BP distribution is obtained from the sum of the pole and cut contributions to the MGF.
For the PDF, N denote the the tridiagonal matrix with main diagonal ) and constant upper sub-diagonal with common value −r.Then, Let M (1) N denote the the tridiagonal matrix with main diagonal (σ ) and constant upper sub-diagonal with common value −r.Then, In terms of the matrix M in (6), M N = sI − M and M (1) N with the first row and column removed.We have the recurrences det M (j) for each j = 0, 1.Only the X N ≡ det M (0) N are significant for general computation of the MGF.Setting σ n = σ 1 + μn with μn ≡ µ n − µ 1 , we see that X N is obtained from the recurrence for n = 2, 3, . . ., N , subject to X 0 = 0, X 1 = 1.To obtain an explicit expansion in powers of σ 1 , we write for n = 0, 1, 2, . ... This yields the simple two-dimensional recursion for 0 ≤ ℓ ≤ n − 3. It is straightforward to show that for n = 1, 2, . .., which is required to seed the recursion.Alternatively, one may achieve this by introducing x (n) n ≡ 0, n ≥ 0. This scheme enables us, in principle, to calculate all the coefficients Table 1 provides expressions for N 0,1 , D 0,1 for server numbers N = 2, 3, 4.
The critical value of the total traffic intensity is the value r = r c at which the pole contribution disappears.This occurs when D 0 (σ and for N = 4, it solves the cubic In Table 2, we present numerical values for r c corresponding to the first few values of N .

Asymptotic Limits
There are two natural ways to scale the problem in the limit of a large number of servers N → ∞.The first is to keep the total traffic intensity r constant as N → ∞.The second is to keep the mean arrival and service rates (namely λ, µ, respectively) constant.As either one of the two time-dimensional parameters simply serves to calibrate the clock in the model, we can choose to set µ = 1.This is equivalent to measuring model time in units of the mean service time, and physical time can be recovered by appropriately rescaling the time axis of the resulting distributions.Thus, the model essentially depend on only two parameters, which can be taken to be the number of servers, N , and either r or the partial traffic intensity ρ ≡ λ/µ [12].An advantage of the former parametrization is that the ergodic region of interest here is characterized by a common bounded interval (viz.0 ≤ r < 1) independent of N .
The limit of an infinite number of servers with given constant ρ = r/N is referred to as the M/M/∞ queue, and corresponds to our second scaling option.It has been previously extensively studied by Guillemin and Simonian [13], and also earlier by Morrison et al. [20].However, while this model is useful in some applications, it is also rather trivial as a queueing model given that it does not actually contain a queue.We shall proceed to study the BP distributions of both scaling options.
Having explicit results for the N → ∞ limits, and comparing these with full calculations for finite given values of N , enables one to decide, in creating and optimizing numerical algorithms, how far one needs to go before simply being able to invoke asymptotic results.

Constant-ρ Model
To study this limit, it is convenient to consider the transpose problem.Let A ≡ sI − M , where M is given by (6) 2 .In the large-N limit, M becomes countably infinite dimensional and the boundary condition for (3), viz.ϕ N = ψ − (s)ϕ N −1 , disappears as it gets pushed to infinity.In Section 3, we considered a linear system which, in the large-N limit, now reads A|ϕ⟩ = |e 1 ⟩ since we can set |w⟩ = |e 1 ⟩, noting that µ 1 = 1 here.This can be inverted to yield the BP MGF ϕ 1 (s) = ⟨e 1 |A −1 |e 1 ⟩.However, we have Thus, we can also consider the transpose problem and compute φ 1 (s) Clearly A T = sI − M T and, in the N → ∞ limit, the infinite tridiagonal matrix M T is given by with µ n ≡ n being adopted here.On setting φ 0 ≡ 1/ρ, the implied recurrence relation for φ n reads for n = 1, 2, . ... Following [13], we next introduce the generating function so that g(0) = φ 1 = ϕ 1 .Summation over the recurrence (113) followed by some standard manipulations gives rise to the differential equation Solution via the integrating factor method yields The vanishing of the residue at z = 1 implies that Guillemin and Simonian [13] have defined the integrals In terms of these, we can write the result for the BP MGF as This result coincides with the MGF of [13] for their congestion duration (denoted θ) in the case C = 0.These authors, however, did not proceed to recover the distribution from the MGF.We shall also introduce the functions 1 (s).We note from the defining matrix equation A|ϕ⟩ = |e 1 ⟩, that ϕ 1 (s) has poles at the eigenvalues of the matrix M T .Similarly, ϕ (1) 1 (s) has poles at the eigenvalues of the truncated problem obtained by removing the first row and column from the matrix M T , and these correspond to the zeros of I 1 (s − 1, ρ).
At this point, it is useful to recall that the Kummer hypergeometric function [23], also known as the hypergeometric function of the first kind, M (a, b, z) = F 1 1 (a; b; z) has integral representation In terms of this function, we can write Thus, the requirement that ϕ 1 (0) = 1 is recovered and, noting that we see that the mean BP for N → ∞ is (e ρ − 1)/ρ as expected from the explicit exact result (56).Furthermore, a function that is entire in all arguments is obtained via With this definition and the relationship we arrive at the result ϕ In this expression, both numerator and denominator are entire functions, and we may consequently observe that the zeros of M(s, s + 1, ρ) are identified with the eigenvalues of the full problem, while the zeros of M(s, s + 2, ρ) are identified with the eigenvalues of the truncated problem.In fact, it follows directly from Cramer's rule that where A 11 denotes A with the first row and column removed.
In view of the foregoing discussion, we may write the BP MGF as where an exact expression is produced when L = ∞.This product form of linear factors is directly implied by ( 126), but we also know that the χ ℓ , χ (1) ℓ correspond to the eigenvalues of the matrices A, A 11 , respectively.The two sequences of eigenvalues (ordered according to increasing modulus) are interleaved, and both tend to negative integral values as ℓ increases.Moreover, they can be paired up in such a way that χ ℓ+1 − χ (1) ℓ → 0 as ℓ → ∞.The following inequalities hold: for ℓ = 1, 2, . ... Therefore, the individual ratios in the product above eventually become indistinguishable from unity, implying that only a finite number, L, needs to be retained.The residue theorem may be applied to recover the SF as an exponential mixture with the weights being given by for k = 2, 3, . ... The weights W k are positive, sum to unity and tend to zero for large k.It follows that, if we set U k ≡ −1/χ k , then the node/weight pairs (U k , W k ) define a texture distribution of the sort discussed in Section 3. The χ ℓ are bracketed according to ℓ − 1 < −χ ℓ < ℓ for all ℓ = 1, 2, . . .while, for the χ ℓ , we may write ℓ < ℓ + 1 for ℓ = ceil(ρ), . . ., −χ (1)  ρ = ρ for integral ρ ≥ 1 . ( This behaviour can be observed in Figure 4 for the case ρ = 1.9.The function values on the vertical axis have been logarithmically compressed such the values indicate the exponent of the order of magnitude.The bijective compression function that was employed is given by where sgn(x) denotes the sign function.Thus, y = f (x) can be inverted according to x = sgn(y)(10 |y| − 1).A bracketed root-finding method, such as bisection or a secant method, can be used to compute the desired number of zeros of the relevant Kummer functions.Alternatively, an eigenvalue problem can be solved for the matrices A, A 11 truncated to a sufficiently large finite dimension (which will be greater than L).An initial dimension can be estimated and subsequently iterated until convergence is achieved.The eigenvalue of least modulus will dominate the large-t behaviour of the distribution and provide a lower bound to the hazard function H(t) ≡ P (t)/ F (t) > |χ 1 |.This asymptotic level is indicated in Figure 5 as the dotted line.The hazard function for the M/M/∞ model corresponds to the dashed black curve.It is clear that this curve is approached rapidly as the number of servers increases.

Constant-r Model
We shall first consider the boundary case r = 1 when N ≫ 1.Let m ≡ (Λ/µ 1 ) 2 with Λ given by (17), noting that m → ∞ as N → ∞ since, by Stirling's formula, m ∼ N ≫1 e 2N /(2πN ).We already know that F (mt) ∼ t→∞ 1/ √ πt for any N , and we can write Thus, F (mt)    But also which leads to the result where erfcx(t) ≡ e t 2 • erfc(t) is the scaled complementary error function, and erfc(t) denotes the standard complementary error function.In completing this derivation, we have used the identity for the incomplete gamma function Γ( The former integral expression is amenable to efficient numeral evaluation via Gauss-Laguerre quadrature.The asymptotic BP PDF is given by The log-mean is given by ⟨ln This result follows easily from direct integration of the representation In the ergodic region 0 < r < 1, the BP distribution tends to a mixture of two exponentials as N become large: one that dominates in the tail for large times and one that dominates for very small times.For a fixed number of servers N , there is always a critical total traffic intensity r = r c beyond which there is no contribution to the distribution from discrete poles.On the other hand, for any fixed traffic intensity r, there exists at least one discrete pole when the number of servers is sufficiently large, and the position of smallest-magnitude pole ζ 1 tends to zero as the number of servers increases.The exponential term generated by smallest-magnitude pole also dominates the large-t behaviour of the distribution.Given that ζ 1 approaches zero for N ≫ 1, it follows that this exponential asymptotic behaviour will set in ever earlier in t as N continues to increase.Hence, we must have that for all t bounded away from zero, where m bp ≡ ⟨T bp ⟩ T bp denotes the mean BP.A corollary of this is that ζ 1 ∼ N ≫1 −1/m bp , which is easily verified numerically by computing ζ 1 from (66) and comparing with the exact result for m bp .The constant ν is derived from the residue of the pole at where the function J(s) has been given in (52).It also holds that ν → 1 as N → ∞.On the timescale that corresponds to measuring t in units of the mean BP, we have a PDF of the mixture form where δ(t) denotes the Dirac delta function.The SF for the two-exponential mixture that leads to this PDF is given by with m ′ = m bp − (1 − ν)m ′′ chosen to reproduce the correct mean, and m ′′ is a constant of order unity estimated as , which strengthens the relationship given above.It is easily confirmed numerically that this result remains accurate down to quite small values of N .In Figure 7, we have plotted the SF exponents − log 10 F (t) for an increasing sequence of server numbers N , where values on the time axis correspond directly to the model definition in the introductory section with the choice µ = 1/N adopted there.Thus, time is being measured in units of 1/N times the mean treatment time.If we were to consider a situation where both r and µ were to remain constant with respect to N , in which case λ would be linear in N , then any given value on the time axis would be a different point in physical time for each curve.It can be seen from the figure that, in the short-time regime being plotted, agreement between the numerical exact curves and the analytical asymptotic curves improves rapidly with increasing N .
The long-time behaviour is presented in Figure 8, where values on the time axis are measured in units of the mean BP.The black dashed curve is the analytical asymptotic curve evaluated here for N = 80, and is seen to coincide with its numerical exact counterpart to within the line-width of the graph.The black dotted line is the limiting exponential that represents the strict N → ∞ limit F (m bp •t) ∼ N →∞ e −t .It constitutes a lower bound to the tail of the SFs for all N , and is approached from above as N → ∞.

Complex-Pole Method
The explicit expression for the cut function developed in the previous section allows us to bypass the spectral decomposition and derive explicit closed-form results for the cut contribution to the BP PDF for an arbitrary number of servers.These comprise sums of Bessel and Marcum Q-functions where the individual terms are generated by poles in a complex plane that are induced by the zeros of the cut polynomial.In fact, it will become clear that knowledge of the cut polynomial is all that is required to completely determine the BP distribution for any given number of servers.While this method is less robust numerically than the spectral approach, it provides more insight into the analytic structure of the problem.Busy period (number of means) The general form of the cut contribution to the BP PDF van be expressed as where the contour C denotes the anti-clockwise traversed unit circle about the origin, and R cut (s) is the cut function, with This follows from setting v = cos θ in (38), and then writing cos θ = −(z + 1/z)/2 where z lies on the unit circle in the complex z-plane3 .For N = 1 (µ 1 = 1), the cut function is the trivial constant while, for N = 3 (µ 1 = 1/3), By writing the cut polynomial in terms of its roots, ), we have the general formula The Schlaefli formula for the modified Bessel function of the first kind is given by where the closed contour is taken around the origin.This may be used to derive the generating function from which it follows that Thus, suppose that g(z) is a meromorphic function such that g(z) = g(1/z) and that C denotes the anti-clockwise unit circle.Then we have the integral Suppose now that G(z) ≡ g(z)/z has only simple poles within the unit circle at z = β 1 , β 2 , . . ., β M .Then we have where the coefficients, for n = 0, 1, . .., are the residue sums For N = 1, (144) and the Schaefli formula (149) reproduce the well-known result [4] P cut (t) = e −(1+r In the general case, for N > 1, the relationships lead to the form where Hence, (144) becomes where (152) applies with which allows us to the identify the c n in the relevant instance of (153).We note that for N ≥ 3 there is no pole at the origin.Next, we write where β ℓ is chosen to be the unique root of the pair of roots β ℓ = −α ℓ ± α 2 ℓ − 1 that lies within the unit circle |β ℓ | < 1.In each pair, one always lies inside and one outside the unit circle, and they are complex conjugates.With this in mind, we can express (154) as with In terms of these quantities, the cut PDF for N ≥ 3 reads It is useful to recall that c 0 = − 2N −3 k=1 γ k .The generalized Marcum Q-function for order ν is defined by the integral which allows it to be extended to complex-valued arguments a, b.The foregoing series is useful for numerical computation when |b/a| < 1.Otherwise, one may appeal to the alternative form The Marcum Q-function was originally introduced in radar detection theory for non-fluctuating targets [19], and has subsequently found applications in communications and signal processing.Abate and Whitt [1] were first to use it in the context of queueing theory.Here we adopt it to write (163) in the more compact form where and these are in general complex valued.A robust numerical algorithm for the computation of the Marcum Q-function 4 that allows for complex arguments has been given in [11].
For N = 2, there is an additional pole at the origin.In this case, we have from which it is evident that a pole at the origin contributes when n = 0, 1.We find that c and it turns out that β 1 = min( √ 4r, 1/ √ 4r), so that Busy period (number of treatment times) Note.BP distribution via the complex-pole method for r = 0.9, N = 15.The NW graph is the cut function Rcut(•) on a logarithmic scale as a function of τ .The NE graph is the unit circle in the complex z-plane displaying the locations of α-poles (blue) and β-poles (red) within the unit circle.The single real α-pole lies just outside the unit circle.The SW graph compares the survival function from the spectral method with that from the complexmethod.The SE graph is the same comparison for the negative base-10 logarithm of the survival function and emphasizes agreement in the tail.This is similar to the expression for the BP distribution as an infinite series of Bessel functions that was previously derived by Arora [3] for the two-server case.In terms of the complementary Marcum Q-function, the two-server result becomes and we may note that In this instance, the arguments of the Marcum Q-function are real.

Summary Statistics
7.1.Analytical In order to test how the numerical algorithms for computing the BP distribution, developed in this work, perform across a wide range of input parameters (r, N ) in the ergodic region, we plot various summary statistics against a grid of these parameters.We can then ascertain (i) whether the observed variations are sufficiently regular and continuous, and (ii) how well the values approach the analytically calculated asymptotic limits as the number of servers grows.The results presented here are also useful in determining the values of N beyond which the asymptotic limit can serve as a proxy for the exact distribution, given a desired level of accuracy.
Since the mean and higher moment diverge rapidly as the traffic intensity r approaches unity, where heavy-tailed behaviour sets in, we seek summary statistics that remain finite at r = 1 and can be computed conveniently.For this purpose, we have chosen the differential entropy H ≡ ⟨− ln(P (T ))⟩ T and the log-mean L ≡ ⟨ln(T )⟩ T .The log-mean is widely used in radar detection theory to estimate the shape parameters of candidate heavy-tailed distributions for high-resolution radar clutter from collected experimental data [27].
Figure 10 presents the BP differential entropy plotted against number of servers N for various values of the traffic intensity r.It includes a comparison with asymptotic limits as detailed in the caption.Figure 11 presents the BP log-mean plotted against number of servers N for various values of the traffic intensity r.It also includes a comparison with asymptotic limits as detailed in the caption.In Figure 12, the BP mean is plotted against traffic intensity r as its base-10 logarithm, for various numbers of servers N , and is compared with the known exact results.All these graphs demonstrate consistent behaviour for the numerical computations, agreement with exact results and the expected approach to asymptotic limits.

Simulation
We have implemented a discrete event simulation (DES) for the M/M/c system.Adhering to the ergodic region (0 < r < 1), we run the simulation for a long time period as a steady-state simulation that is initialized to the empty state, and collect the times (epochs) at which an empty system becomes non-empty and at which a non-empty system becomes empty.These data are used to generate the empirical partial BP distribution from which summary statistics can be estimated.The goal is to compare the analytical results for various summary statistics with their empirical estimates across a wide range of the input-parameter pairs (r, N ) in order to verify the theoretical results derived in this work.A simulation is run for a given parameter pair whenever the expected number of BPs generated in the maximum tolerated simulation time T stop is at least 10.The maximum simulation time was taken to be 10 6 multiples of the mean treatment time.According to the theory of regenerative processes [8] the expected number of regeneration cycles (empty periods followed by a BP) in time T stop is simply given by n reg = T stop /[(m bp + 1/r)•T scl ] where T scl is a timescale needed to ensure that time in both the simulation and the analysis are measured in the same units.In the simulation, we set the mean treatment time to unity so that the simulation clock runs in units of the mean treatment time.In this case we must set T scl = 1/N .
The nearest-neighbour estimator [17,30] is used for the empirical differential entropy.This choice avoids calculation of a kernel density estimator (KDE) for the empirical PDF that is required by other empirical entropy estimators [5].The KDE is problematic for distributions with semi-infinite support because undesirable artefacts are inevitably produced at the boundary.For a sample of size M {T i : 1 ≤ i ≤ M }, the nearest-neighbour estimator is given by where ρ i ≡ min j̸ =i ∥T i − T j ∥ is the nearest-neighbour Euclidean distance of T i from all other members of the sample.Due to the small sample numbers in many cases (viz.r close to unity, or large N ) and lack of alternatives for the entropy, the bootstrap method [14] has been used to generate the confidence intervals that provide the error bars in the graphs.The significance level is taken to be α = 0.01.Thus, the displayed error bars in Figures 13 and 14, for the differential entropy and log-mean, respectively, indicate the central 99% confidence intervals.In cases where only a dot appears, the length of the confidence interval is less than the dot size.
The nearest-neighbour entropy estimator is somewhat problematic for the bootstrap method as it generates duplicates, for which the nearest-neighbour distance is zero, leading to logarithmic divergence.We have adopted two strategies to address this difficulty: (i) removal of duplicates after bootstrap resampling but prior to evaluating the estimator, and (ii) constructing an ensemble of nearest-neighbour distances from the entire sample and bootstrapping on the nearest-neighbour ensemble.While these mitigations have different draw-backs, they are observed to perform equally well when tested on a large variety of simple distributions and compared with accurate results for the confidence intervals obtained from repeated Monte Carlo simulation trials.
Figure 13 presents the BP differential entropy versus traffic intensity r for various server numbers N , with overlaid simulation data as indicated by the dots with error bars.The green dashed curve is the extreme large-N limit based on a single exponential given by F (t) ∼ N ≫1 e −t/m bp with the mean m bp computed for N = 60.The dashed black curve is the finer large-N asymptotic form based on the two-exponential mixture given in (143) and computed for N = 60.It aligns very closely with the exactly calculated result.Figure 14 presents the BP log-mean versus traffic intensity r for various server numbers N , with overlaid simulation data as indicated by the dots with error bars.The green dashed curve is the extreme large-N limit based on a single exponential given by F (t) ∼ N ≫1 e −t/m bp with the mean m bp computed for N = 50.The dashed black curve is the finer large-N asymptotic form based on the two-exponential mixture given in (143) and computed for N = 50.It aligns very closely with the exactly calculated result.

Conclusions
This paper has developed two distinct methods for generating explicit exact results for the distribution of the partial BP pertaining to the M/M/c queue and a variety of models that generalize it with priority levels and distinct arrival classes.The spectral method allows for a robust and efficient numerical implementation.The algebraic method furnishes closed-form results for the PDF and SF that elucidate the analytical structure of the problem.In particular, for any given number of servers, it has identified a unique polynomial that completely characterizes the distribution.The present discussion has also served to connect previous diverse approaches to the problem.

Figure 1 .
Figure 1.Note.BP pole locations for case r = 0.5, N = 30 at s = ζ ℓ , ℓ = 1, 2, 3 are indicated by the crosses on the horizontal axis.The smallest pole lies very close to zero.

9 Figure 2 .
Figure 2. Note.BP survival function for various traffic intensities (r) and Nsrv = 30 servers, showing long-time behaviour.The negative base-10 logarithm of the survival function is plotted.

9 Figure 3 .
Figure 3. Note.BP survival function for various traffic intensities (r) and Nsrv = 30 servers, showing short-time behaviour.The negative base-10 logarithm of the survival function is plotted.

Figure 4 .
Figure 4. Note.Zeros of the Kummer M-functions in the numerator and denominator of the M/M/∞ BP MGF for partial traffic intensity ρ = 1.9.Denominator zeros are indicated by crosses.Numerator zeros are indicated by dots.

Figure 5 .
Figure 5. Note.BP hazard function for various numbers of servers (Nsrv) and partial traffic intensity ρ = 0.75.The dashed black curve is the hazard function for the M/M/∞ model.The dotted black line is large-t limit and lower bound of the asymptotic hazard function.

Figure 6 .
Figure 6.Note.Deviation of the BP survival function from its large-N asymptotic form for various numbers of servers (Nsrv) and traffic intensity r = 1.The difference of the survival function logarithms is plotted on the vertical axis.Time on the horizontal axis is measured in units of the scale m = (Λ/µ1) 2 .

Figure 7 .
Figure 7.Note.BP survival functions for various numbers of servers (N ) and traffic intensity r = 0.9, compared with their respective large-N asymptotic forms.These curves elucidate the very short-time behaviour.The asymptotic forms are evaluated at the corresponding values of N .

Figure 8 .
Figure 8. Note.BP survival functions for various numbers of servers (Nsrv) and traffic intensity r = 0.75, showing long-time behaviour.The dashed black line is the asymptotic form evaluated for N = 80.The black dotted line is the strict N → ∞ limiting exponential.
where I ν (x) is the modified Bessel function of the first kind of order ν.It constitutes a cumulative distribution function (CDF) in the variable b.The original Marcum Q-function Q(a, b) is the special case of the generalized Q-function Q ν (a, b) for ν = 1.Its complementary function Q(a, b) ≡ 1 − Q(a, b) can be represented by the infinite Neumann expansion Q(a, b) = e −(a 2 +b 2 ) Figure 9.Note.BP distribution via the complex-pole method for r = 0.9, N = 15.The NW graph is the cut function Rcut(•) on a logarithmic scale as a function of τ .The NE graph is the unit circle in the complex z-plane displaying the locations of α-poles (blue) and β-poles (red) within the unit circle.The single real α-pole lies just outside the unit circle.The SW graph compares the survival function from the spectral method with that from the complexmethod.The SE graph is the same comparison for the negative base-10 logarithm of the survival function and emphasizes agreement in the tail.

Figure 12 .
Figure 12.Note.BP mean versus traffic intensity r, plotted as its base-10 logarithm.For each number of servers N , the computed coloured dashed curve overlays a solid black curve representing the exact analytical result.Discrepancies are less than the linewidth.