Spectral gap of the Erlang A model in the Halfin-Whitt regime

We consider a hybrid diffusion process that is a combination of two Ornstein-Uhlenbeck processes with different restraining forces. This process serves as the heavy-traffic approximation to the Markovian many-server queue with abandonments in the critical Halfin-Whitt regime. We obtain an expression for the Laplace transform of the time-dependent probability distribution, from which the spectral gap is explicitly characterized. The spectral gap gives the exponential rate of convergence to equilibrium. We further give various asymptotic results for the spectral gap, in the limits of small and large abandonment effects. It turns out that convergence to equilibrium becomes extremely slow for overloaded systems with small abandonment effects.


Introduction
Within the fields of stochastic processes and queueing theory, the Halfin-Whitt regime refers to a mathematical way of establishing economies-of-scale in many-server queueing systems like call centers (see [13]). The Halfin-Whitt regime in fact prescribes a scaling under which the manyserver systems converge to limiting processes, which are for most systems diffusion processes. This paper deals with many-server systems in the Halfin-Whitt regime with the additional feature that customers are impatient, and may abandon the system without being served. For such systems with abandonments, we are interested in the spectral gap, which is inversely related to the relaxation time or the speed at which a system reaches stationarity. A large relaxation time in general indicates that replacing time-dependent characteristics by their stationary counterparts might lead to poor approximations. As it turns out, the rate at which customers renege (abandon the system) greatly influences the relaxation time.
In recent years, a large number of papers have dealt with the influence of reneging or abandonments on the system behavior (see, e.g., [8,15,35,36,37,41,42]), and it is widely accepted that reneging is indeed one of the main factors driving the system performance. One of the key insights is that the system behavior strongly depends on whether it is stable or overloaded. By stable we mean that it can serve all customers, even if none of the customers would abandon the system. The spectral gap (or relaxation time) is also very different for stable or overloaded systems. In fact, we find that stable systems have a relatively short relaxation time, whereas the relaxation time of overloaded systems can become extremely large, particularly when the reneging rate is small.
The model we shall consider is the M/M/s + M system, better known as the Erlang A model. This model is a standard Markovian many-server queueing system with Poisson arrivals, exponential service times, s servers, and with the additional feature that customers that are waiting in the queue abandon the system after exponentially distributed reneging times. The queue length process in the Erlang A model, denoted by (Q(t)) t≥0 , is a birth-death process. Whitt [36] (see also [38] and [19]) derived a fluid approximation for the the steady-state behavior of the overloaded Erlang A model, and he further showed that a diffusion limit might provide refined approximations. Garnett, Mandelbaum and Reiman [15] proved a diffusion limit for the Erlang A model in the critical regime. In particular, they showed that under certain conditions a sequence of normalized queue length processes converges to a certain diffusion process (X(t)) t≥0 . These conditions are in fact the ones that correspond to the Halfin-Whitt regime, in which the arrival rate λ and the numbers of servers s are scaled such that, while both λ and s increase toward infinity, the traffic intensity ρ 0 = λ/s approaches one, with ( The diffusion process (X(t)) t≥0 is a combination of two Ornstein-Uhlenbeck (OU) processes with different restraining forces, depending on whether the process is below or above zero. The number of customers in the Erlang A model can be roughly expressed as s + √ sX(t) for s sufficiently large. The diffusion process is generally easier to study than the birth-death process, and can thus be employed to obtain simple approximations for the system behavior. The steady-state distribution of the diffusion can be easily obtained (see (6) below), but less is known about the time-dependent behavior. In this paper we shall present an explicit and asymptotic characterization of the spectral gap of (X(t)) t≥0 . The spectral gap of the diffusion process provides an understanding of the relaxation times for the Erlang A model in the Halfin-Whitt regime. Most importantly, we shall study in detail the impact on the spectral gap of the capacity parameter β and the reneging rate η, which shall enhance our understanding of how the Erlang A model behaves for positive/negative β and small/large values of η.
The diffusion process (X(t)) t≥0 also applies to the G/M/s + M system, in the same asymptotic limit, which was proven by Whitt [37]. Stochastic processes for more general systems with abandonments were obtained recently by Dai, He and Tezcan [8] for the G/P h/s + M system. In this case, the limiting process is still a diffusion process, but it becomes multi-dimensional. Zeltyn and Mandelbaum [42] derived approximations for the M/M/n + G in the Halfin-Whitt regime. In case of general service times, the limiting process is not even a diffusion process (see e.g. [27,37] for cases without reneging). Therefore, the one-dimensional diffusion process (X(t)) t≥0 strikes the proper balance between simplicity and tractability, while retaining the essential features of abandoning customers in many-server systems.
The Erlang A model is particularly interesting, as it incorporates three classical queueing systems as special cases. In the case of no reneging (with η = 0) the Erlang A model reduces to the Erlang C model, or M/M/s system. Halfin and Whitt [17] established that the limiting process behaves as a Brownian motion above zero and an OU process below zero. In [23] we have referred to this process as the Halfin-Whitt diffusion. For η = 1 the Erlang A model becomes an infinite server queue or M/M/∞ system, for which the stochastic-process limit is known to be an OU process [18]. For η → ∞ the Erlang A model becomes the Erlang B model or M/M/s/s system, in which case the stochastic-process limit is a reflected OU process (see [12,24,35]). For the diffusion approximations, we show the η → ∞ reduction in Section 3. 3. Our analysis of the spectral gap of (X(t)) t≥0 provides results for each of these three cases.
Mathematically, determining the transient distribution for the present diffusion process involves analyzing a Schrödinger type equation (see Section 4) with a piecewise parabolic potential function, or, equivalently, a Fokker-Planck equation with a piecewise linear drift (see (2) and (3) below). Such problems arise in a variety of other applications, such as linear systems driven by white noise [3,2], the Kramers' problem [26] and escape over potential barriers [22]. Invariably, the solution involves the parabolic cylinder functions, and these we discuss in detail in Section 5. The key in determining the spectral gap is in fact determining the Laplace transform of the transient probability distribution over time. The spectral gap then follows from the dominant singularity of the Laplace transform. The main results are presented in Section 2, the three special cases (η = 0, 1, ∞) are discussed in Section 3 and the proofs are given in Section 6. Before the proofs we give some basic background on Schrödinger equations (Section 4) and parabolic cylinder functions (Section 5), whose properties are heavily used later. In Section 7 we establish monotonicity properties of the spectral gap.

Main results
The diffusion process (X(t)) t≥0 is a Markov process on the real line with continuous paths and density p = p(x, t) = p(x, t; x 0 ; β, η) that satisfies the forward Kolmogorov equation where and (with δ(·) the Dirac function and p x = ∂p/∂x) p(0 + , t) = p(0 − , t), p x (0 + , t) = p x (0 − , t), and p(x, t) must decay as x → ±∞. The limiting distribution of the diffusion process is (see [15]) p(x, ∞; x 0 ; β, η) = C e − 1 2 ηx 2 e −βx , x > 0, e − 1 2 x 2 e −βx , x < 0, where C −1 = ∞ 0 e − 1 2 ηx 2 e −βx dx + 0 −∞ e − 1 2 x 2 e −βx dx. As shall be discussed in Section 4, this problem has a purely discrete spectrum for all η > 0, and it is confined to the real axis. The spectral gap can thus be defined as the absolute value of the least negative eigenvalue of the operator in the right-hand side of (2). It governs the asymptotic rate of convergence to the stationary distribution. An alternative description of the spectral gap is the absolute value of the singularity closest to the imaginary axis in the range Re(θ) < 0 of the Laplace transformp. Denote this dominant singularity byθ and the spectral gap by r(β, η).
The relaxation time, which measures the time it takes for the system to approach its steady-state behavior, is defined as (see [4,7]) and hence τ −1 = −Re(θ) = r(β, η). For this problemθ is real, so that −θ = r(β, η). Our definition of the relaxation time in (7) assumes the initial condition p(x, 0) = δ(x − x 0 ) in (4), and then the approach to equilibrium is governed by λ 1 = r. But we could certainly have initial conditions that would lead to a faster approach. For example, if p(x, 0) = p(x, ∞) then p(x, t) = p(x, ∞) for all t and equilibrium is attained instantaneously. We could also have initial distributions p(x, 0) that have zero projections on, say, the first L eigenfunctions, and then the sums in (31) and (32) below would be replaced by ∞ n=L+1 e −λnt c n φ n (x), where the c n may be computed in terms of p(x, 0). Then the approach to equilibrium would be governed by eigenvalue λ L+1 .
Here is the main result: The spectral gap of the diffusion process (X(t)) t≥0 is given by r(β, η) = −θ whereθ is the least negative solution to V(θ; η, β) = 0 with D ν (z) the parabolic cylinder function with index ν and argument z, and D ′ We also note that the result in (8) corresponds to taking the limit of the discrete queueing model with the scaling in (1) and then (7) examines what happens for large times. We show in Appendix D that (8) may also be obtained from the exact solution of the M/M/s + M queue. We have Proposition 2. The spectral gap in the discrete M/M/s + M model is the least negative solution to where F n and H n are the contour integrals in (199) and (201). For s → ∞, with ρ 0 = 1 − β/ √ s and θ = O(1), the roots of ∆(θ) may be approximated by those of V in (8). This shows that the exchange of the limit in (1) and of large time is permissible in this particular case.
Theorem 1 is an implicit description of the spectral gap, and it can be used to calculate r(β, η) numerically or asymptotically. For some values of β and η the spectral gap is shown in Figure  1. We observe that the spectral gap decreases with β if η > 1 and increases with β if η < 1. If η = 1 the spectral gap is r(β, 1) = 1 for all β, since then the problem reduces to a standard Ornstein-Uhlenbeck process (see Section 3.2). This suggests that for systems with a large reneging rate, increasing the load (increasing ρ 0 and decreasing β) leads to shorter time scales for achieving equilibrium, while the opposite it true for small reneging rates. We also see that r increases as a function of η. Later, we establish the monotonicity of r with β, which is suggested by the numerical results in Figure 1 (see Section 7). To further substantiate our findings, we accompany the observations from Figure 1 by results for the spectral gap in various asymptotic regimes. In order to do so we assume that η → 0 (small abandonment rate). The asymptotics for η → ∞ can be obtained from the following important symmetry result, which we establish in Section 6.

Proposition 4.
For β < 0 the spectral gap behaves asymptotically as r(β, η) ∼ η with the correction term and hence r − η is exponentially small as η → 0. Proposition 4 describes the part at the far right end of Figure 1, where r increases linearly with η. Since β < 0, the diffusion process is mostly in the positive part of the state space, since the process has a positive drift for 0 < x < −β/η and an equilibrium point at x = −β/η = |β|/η ≫ 1. Hence, particularly when there is little reneging, one has to be far up in the state space before the process starts stabilizing. For the underlying queueing model, this scenario corresponds to large queues building up until enough customers renege so that the situation stabilizes. For this scenario, the spectral gap r = O(η) suggests large relaxation times. Note also from (6) that the steady-state distribution concentrates about x = −β/η. Table 2 compares exact and asymptotic results for β = −1.  We next consider β positive and sufficiently large, where we obtain a very different result for r.
The equation (16) corresponds to the discrete part of the spectrum of the Halfin-Whitt diffusion with no reneging (i.e., with η = 0 in (3)); see the discussion after Proposition 13. The various solution branches of (16) are demonstrated in Figure 2, where we plot the implicit functionṼ(p, β) = 0 for β, p > 0.
Proposition 5 applies to the flat part in Figure 1. Indeed, when β is large enough, the spectral gap is hardly influenced by η. The diffusion process will spend most time below zero, near x = −β. A likely queueing scenario would be that queues hardly ever build up, which makes the impact of reneging customers negligible. As the spectral gap is O(1), we expect relaxation times that are O(1). While asymptotically r(β, η) ranges from r 0 (β) to 1, numerically this corresponds to the interval (β 2 * /4, 1) = (.86231, 1), which is quite small, leading to the flatness of the surface in Figure  1 for β > β * . Table 3 compares exact and asymptotic results for β = 2.
We next consider β > 0 but with β < β * , in which case (16) has no positive solutions.
From Proposition 6 we see that the asymptotic series now involves powers of η 1/3 , which illustrates the lack of analyticity of r(β, η) at η = 0. Now r(β, η) ∼ β 2 /4 and for η = 0 the spectral gap is in fact exactly 1 4 β 2 (see Section 3.1). Table 4 compares exact and asymptotic results for β = 1. When β becomes small, both (14) and (18) become invalid, as the correction terms become larger than the leading term. Then a separate analysis leads to the following result.  Table 3: Results for β = 2.  Table 4: Results for β = 1.
The results in (21) and (22) are also consistent with Figure 1, which suggests that r(0, η) increases as a concave function of η. The queueing counterpart is such that the load is one, and hence the reneging is necessary to alleviate the system. As η becomes larger, more customers will leave the system, which reduces the queue lengths and, as seen from Proposition 7, shortens the relaxation times. Table 5 displays numerical results for γ = 1, 0, −1.  Table 5: Results for γ = 1, 0, −1.
As discussed in (57)-(59) in Section 4, the Sturm-Liouville ODE theory guarantees that there are infinitely many real solutions to (24). The solution branches of (24) are illustrated in Figure 4.
We also note that if we order the roots of Ai(z) = 0 as 0 > a 0 > a 1 > · · · and the roots of Ai ′ (z) = 0 as 0 > b 0 > b 1 > · · · , these roots interlace as 0 > b 0 > a 0 > b 1 > a 1 > · · · , and this fact can be used to establish more directly that (24) has infinitely many solution branches, for any fixed W .
We have thus obtained the asymptotic connection between Propositions 5 and 6. Numerical results for the case β = β * are given in Table 6.
This concludes our asymptotic analysis of the spectral gap in Theorem 1. The proof of Theorem 1 follows immediately from an explicit expression for the Laplace transformp of the transient Defining the auxiliary function M by we have the following result: Theorem 10. Consider x 0 < 0, with V in (8) and M in (27), and assume that Re(θ) > 0.
The proof of Theorem 10 is presented in Section 6.2. Note that the results for x 0 > 0 follow immediately from the symmetry relation (11).
Hence, the large-time behavior of the diffusion process is dominated by the least negative zero of V, which gives the result on the spectral gap in Theorem 1.
From (28)- (30), by evaluating the contour integral for the inversion of the Laplace transform p(x; θ), we can obtain a spectral expansion of the form and where and The eigenfunctions ψ ± n then satisfy the orthogonality relation When η = 1 we have λ n = n, k n = D 2 n (β)/(n! √ 2π), and D n (β) = e −β 2 /4 He n (β) so that (36) reduces to Note that the pole at θ = 0 of (28)- (30) corresponds to the steady state behavior p(x, ∞), while the poles at θ = −λ N and their residues lead to the decaying terms in (31) and (32). However, the spectral expansion does not yield any more insight than (28)-(30).

Three special cases
We shall now consider the three special cases of the diffusion process that arise by setting η equal to zero, one and infinity.

The Halfin-Whitt diffusion
As η → 0 we end up with a process that behaves like a Brownian motion with drift above zero and like an Ornstein-Uhlenbeck process below zero. In [23] we have called this diffusion process the Halfin-Whitt diffusion, after Halfin and Whitt [17] who identified this process as a heavy-traffic limiting process for the GI/M/s system. The mean hitting time of the Halfin-Whitt diffusion was obtained in Maglaras and Zeevi [25]. Gamarnik and Goldberg [14] were the first to identify the spectral gap of the M/M/s system, asymptotically in the Halfin-Whitt regime.
To establish Theorem 12, Gamarnik and Goldberg used the framework of Karlin and McGregor [20] for birth-death processes, and the result of Van Doorn [9] on the spectral gap of the M/M/s system. In [14] the starting point is the discrete M/M/s model, and its spectral gap is then analyzed in the Halfin-Whitt regime (1). An alternative proof of Theorem 12 was given by the authors in [23] by deriving the expression for the Laplace transformp of the transient density in the diffusion limit, which shows that the limits of large time and (1) may be, in this case, interchanged. Below we summarize the main result in [23].

Free-space OU process
When η = 1 it immediately follows from the process description that the diffusion process (X(t)) t≥0 reduces to a free-space OU process, for which it is known that (with Indeed, this result also follows from Theorem 10 using the Wronskian identity in (69), which shows that Expression (40) is obtained for example in [34], in the context of the harmonic oscillator (see (46)-(50) below). Also, M → 0 (cf. (27)) as η → 1 and then (40) follows from (28)- (30). It is easy to invert the Laplace transform (40), as its poles are at zero and at all negative integers. Hence, (see, e.g., [21]) Here D n (z) = e −z 2 /4 He n (z) where He n (·) is the nth Hermite polynomial. Alternatively, there is the closed-form expression

Reflected OU process
As η → ∞, the process will spend all its time below zero, and hence (X(t)) t≥0 reduces to a reflected OU process (see Ward and Glynn [35], Linetsky [24] and Fricker et al. [12]). In this limit we have , and then, using (67), which can be used to simplify (29) and (

Schrödinger equations and spectral properties
Here we give some basic background on spectral properties that are relevant to PDE's such as (2).
In particular we show that the discreteness of the spectrum for any η > 0 follows from classic results on the Schrödinger equation. We set p = e −λt φ(x) where λ is a spectral or eigenvalue parameter. Then (2) and (3) lead to and the interface conditions are φ(0 − ) = φ(0 + ) and φ ′ (0 − ) = φ ′ (0 + ). Furthermore, we can transform (46) and (47) into the self-adjoint form by setting which leads to the Schrödinger equation where E and λ are related by and the "potential" function V(x) is We also require the eigenfunctions ψ(x) to satisfy ψ(0 + ) = ψ(0 − ), and ψ ′ (0 + ) = ψ ′ (0 − ). Since the problem is defined over the entire real line, additional conditions must be imposed at x = ±∞, and most often it is required that However, for parabolic and piecewise parabolic potentials, such as the one in (51), this condition is equivalent to simply rejecting solutions of (49) that have Gaussian growth as x → ±∞.
We note that if the potential V(x) were exactly quadratic, say V(x) = x 2 /4, then the problem is just the quantum harmonic oscillator (or, for our application, the standard Ornstein-Uhlenbeck process), and then the eigenvalues are E N = N + 1/2 and the corresponding eigenfunctions are where c N is a normalizing constant and He N (x) is the N th Hermite polynomial. Thus the spectrum is purely discrete for quadratic potentials.
We can also view the differential equation in (49) as constituting a singular Sturm-Liouville boundary value problem. The study of such problems dates back to the work of Sturm in the nineteenth century, and they are discussed in detail in the books of Titchmarsh [34], Stakgold [29], Reid [28] and Coddington and Levinson [6,. The problem in (49) is singular since it is defined over the infinite interval x ∈ (−∞, ∞).
Singular Sturm-Liouville problems are classified as either of limit circle or limit point type. For limit point problems the condition that the solution be square integrable is sufficient to determine it, while limit circle problems require a more explicit boundary condition at the singular point(s) (which are at x = ±∞ for (49)).
Since Sturm-Liouville problems and Schrödinger equations are self-adjoint, their spectra are confined to the real axis. Singular problems may have both discrete and continuous spectra. However, there is a general result, originally due to Weyl, with simplified proofs by Titchmarsh appearing in [32,33] (see also the book [34]), that guarantees that (49) will have a purely discrete spectrum. This needs only the conditions that Our potential in (51) clearly satisfies these conditions and thus has a purely discrete spectrum, for any η > 0. The fact that V(x) has a jump discontinuity at x = 0 does not affect the spectrum; it only means that some jump conditions must be specified at x = 0. However, if η = 0, then the potential does not grow at x = +∞, and then in fact, as we discussed in [23], the problem has a continuous spectrum in the range λ > β 2 /4 (E > β 2 /4 + 1/2), and may also have any number of discrete eigenvalues, depending on the value of β. Much of the asymptotic work here assumes that η → 0 + , so we are looking at a very singular limit where the discrete spectrum begins to resemble a continuous one, in certain ranges of λ.
For the problem in (51) the smallest eigenvalue is E 0 = 1/2 (thus λ 0 = 0) with the corresponding eigenfunction being the piecewise Gaussian and this corresponds to the steady state distribution in our model. Given the discrete spectrum we order the eigenvalues E N as with λ N = E N − 1/2. By general results for Sturm-Liouville problems the sequence {E N } satisfies E N → ∞ as N → ∞. Also, for every eigenvalue there is only one linearly independent eigenfunction, so all eigenvalues are simple. This can be shown directly from (49), for if ψ(x) and ψ(x) corresponded to the same eigenvalue E, then ψ( and thus ψ must be a multiple of ψ.
There are two other singular Sturm-Liouville problems that are relevant to the analysis here. First consider with the boundary condition This problem has a regular point at X = 0, where a standard boundary condition is applied, and a singular point at X = ∞, where we require that Ψ(X) ∈ L 2 (0, ∞). This is a singular problem of limit point type at X = ∞ which may be explicitly solved in terms of parabolic cylinder functions, with Ψ(X) = D E (X + γ). Then (55) leads to the eigenvalue condition The results of Weyl and Titchmarsh again guarantee that the problem has a purely discrete spectrum and thus an infinite sequence of eigenvalues { E N }. Also, E 0 = 0 is the lowest eigenvalue with Ψ 0 (X) = e −(X+γ) 2 /4 . Note that (56) is essentially the same as equation (20) Thus the existence of infinitely many real solutions to (20) follows from Sturm-Liouville ODE theory, though in the next section we shall also show that it follows from the oscillatory nature of the parabolic cylinder functions, as functions of their index E.
Another singular Sturm-Liouville problem is with the boundary conditions Here ω is a real parameter. Again, since V (x) = x grows linearly as x → ∞ and x = 0 is a regular point, we have a purely discrete spectrum. But (57) is related to the Airy equation, with solutions proportional to Ai(x − E), and the eigenvalues are determined by (58), hence But (59) is equivalent to equation (24) in Proposition 9, so again ODE theory can be used to establish the existence of infinitely many solutions. Note also that if ω = 0 the eigenvalues are the roots of Ai ′ (·), while if ω = ∞ the eigenvalues are the roots of Ai(·).
To summarize we have given some basic background on Sturm-Liouville theory, and on Schrödinger equations and their spectral properties, that are useful in the present investigations. In particular this theory guarantees infinitely many discrete solutions to the equations that arise in Propositions 6 and 7, and in Theorem 1.

Parabolic cylinder functions and their properties
The parabolic cylinder equation is the second order ordinary differential equation where z is a complex variable and p is a parameter. Since (60) has no singular points (except at z = ∞) its solutions are entire functions of z (see [6]). One solution of (60) is denoted by D p (z), which is called a parabolic cylinder function of order p, and it is defined by the integral representation Here Br is a vertical Bromwich contour on which Re(u) > 0, and the branch of u p will be defined by u p = |u| p e i arg(u) where −π < arg(u) ≤ π. Then the integrand in the right-hand side of (61) is analytic exterior to the branch cut where Im(u) = 0 and Re(u) ≤ 0. The function D p (−z) provides a second linearly independent solution to (60), so that is the general solution, with c 1 and c 2 being complex constants. When p = 0, 1, 2, . . . is a nonnegative integer we can obtain D p (z) in a closed form as where He N (·) is the N th Hermite polynomial. Note that if p = N the integrand in (61) becomes an entire function of u. Here we use the notation He(·) for the Hermite polynomials, so that He 0 (z) = 1, another linearly independent solution must be used in (62), but we shall not need it in the present analysis. The function D p (z) is real valued when z and p are real. As discussed in [31], D p (z) is an entire function of both p and z, and indeed one can easily compute derivatives of all orders from the integral representation in (61). For example, we have and In (65) log u is real for u real and positive, and analytic exterior to the cut Im(u) = 0, Re(u) ≤ 0. If z = 0 the integrals in (61) and (64) may be expressed in terms of the Gamma function, with Since Γ(z) has simple poles at z = 0, −1, −2, . . . it follows that D p (0) has simple zeros at p = 1, 3, 5, . . . , while D ′ p (0) has simple zeros at p = 0, 2, 4, . . . . This also shows that the functions in (66), as functions of p, tend to oscillate for p > 0, but have one sign for p < 0 (actually, for all p < 1 for D p (0)). Also, in view of the growth of Γ(z) as z → +∞, the functions in (66) decay very rapidly for p → −∞.
Using (61) and (64) we can easily derive the recurrence relations which we shall use in the present analysis. The Wronskian of D p (z) and and it has a very simple form, with which vanishes if p = 0, 1, 2, . . . . We also note that D p (z) and D ′ p (z) cannot both vanish simultaneously. For if D p (z * ) = D ′ p (z * ) = 0 for some p and z * , then (60) shows that D ′′ p (z * ) = 0. Then repeated differentiation of (60) would show that all derivatives of D p (z) vanish at z = z * . Then we could expand D p (z) in Taylor series about z = z * to conclude that D p (z) = 0 in some neighborhood of z = z * . But since D p (z) is an entire function of z this would imply that D p (z) is identically zero, which is clearly not the case.
To better understand the behavior of these special functions, many asymptotic formulas have been derived for z and/or p large. We summarize some of these below, since they are used to establish our main results. First, for z large and positive, we have (see for example [16, p. 1093 A more general result, which allows z to be negative, is the following (see for example [16, p. 1094]): Here we let z = |z|e i arg(z) where | · | denotes the complex modulus.
The limit of z large and negative corresponds to setting arg(z) = π in (71), and then the leading term becomes In some of the analysis that follows we will need to consider cases where p is not exactly a non-negative integer, but is very close to one. For p = N the first series in (71) has Gaussian decay as z → −∞ (O(e −z 2 /4 )), while the second series in (71) (cf. also (72)) has Gaussian growth as z → −∞. But if p is very close to an integer these two terms may be of comparable magnitude. For example if p = ε is small Γ(−p) ∼ −1/ε and if ε → 0 and z → −∞ in such a way that e z 2 /2 ε is O(1), then the two parts of (71) are roughly comparable.
The asymptotic results in (70)-(72) follow easily by expanding (61), using techniques for the asymptotic evaluation of integrals, such as the saddle point method and singularity analysis. General references for such methods are the books of Bleistein and Handelsman [5], Wong [39], Szpankowski [30], and Flajolet and Sedgewick [11]. The integrand in (61) has a saddle point at u = z and a branch point at u = 0, and one of these (or both) determines the asymptotic behavior of D p (z) as z → ∞, for any direction arg(z) in the complex plane. The saddle leads to (70) and the first part of (71), while the branch point leads to (72) and the second part of (71).
We next consider a fixed (real) z and expand D p (z) in the limits of p → ±∞. Then From (74) we see faster than exponential growth with p, coupled with oscillations, in view of the trigonometric factor. Near a zero of the cosine the O(p −1 ) error term may become important, and it may be also explicitly obtained. Thus D p (z) has an infinite number of zeros as p increases toward +∞, for any fixed real z. This is in sharp contrast to viewing D p (z) for a fixed p as a function of z, in which case it has at most finitely many zeros. It is known [1, p. 696] that D p (z) has no zeros in the range p + 1/2 < z 2 /4. Also, (74) shows that the large zeros can be estimated by The results in (73) and (74) may be obtained, for example, by expanding the integral in (61). When p → −∞ the asymptotics are governed by a single saddle point at u = √ −p, while as p → +∞ two saddle points, at u = ±i √ p, contribute.
We next consider asymptotic limits where p and z are simultaneously large, restricting ourselves to real p and z. The results quoted below are taken out of Abramowitz and Stegun [1], where we note that in [1] the results are given for the function U(a, z), which is related to D p (z) by U(a, z) = D −a−1/2 (z), or D p (z) = U(−p − 1/2, z). A complete summary of the asymptotics of U is also given in Temme [31].
When p + 1/2 > 0 and z 2 − 4p is large and positive, the so-called Darwin's expansions apply, where [1, pp. 689-690] , and for p + 1/2 < 0 and z 2 − 4p large and positive The results in (76) and (77) Here sin −1 (·) ∈ (−π/2, π/2). The result in (78) is uniform for p → ∞ and z → ±∞ for intervals where z 2 /(4p) ∈ [0, 1 − ε] for any ε > 0, except if we are at or near a zero of the cosine function, in which case the correction term(s) to (78) must be considered. For a fixed large z, p has to increase past approximately z 2 /4 in order for the zeros of D p (z) to become evident. The expansions in (76) and (78) develop non-uniformities when z 2 /(4p) ≈ 1 and there yet other expansions apply. The following result [1, p. 689] is more uniform and applies for all z 2 /(4p) ∈ [0, ∞], as long as p and |z| are large: and Here Ai(·) is the Airy function, which has the following asymptotic behaviors as z → ±∞ (see [1, p. 448 For p → ∞ with a fixed τ > 0, we can simplify (79) by using (83) to approximate the Airy function, and then we obtain (76) as a special case, when ξ > 1.
We have thus summarized the various "uniform" asymptotic approximations to D p (z), where both z and p become large. Despite the seeming complexity of these results, they are easily obtained from (61) via the saddle point method. Indeed, setting u = zv (with z > 0) in (61) leads to The saddle point equation is Φ ′ (v) = 0 which is the quadratic equation To expand the integral in (87) for z → ∞, p → ∞ with p/z 2 fixed we find that for 4p/z 2 < 1 the real saddle at v = v + determines the asymptotic behavior, and we ultimately obtain (76), up to a Stirling approximation of Γ(p + 1) and the equivalence z 2 − 4p − 2 ∼ z 2 − 4p. In contrast, when 4p/z 2 > 1 two complex saddles, at 1 2 1 ± i 4pz −2 − 1 , contribute to the asymptotics and then we obtain (78). The transition range in (86) corresponds to 4pz −2 ≈ 1 and then the two saddles coalesce to form a higher order saddle, and such transitions invariably involve Airy functions (see Chapter 9 in [5]) . We have discussed here only approximations to D p (z), but some of our main results involve also the derivative D ′ p (z) (see, for example, Theorem 1). Its asymptotics follow from the integral in (64), but the same results can be obtained by formally differentiating the results for D p (z), as in this case term by term differentiation of the asymptotic series is permissible. For example, the logarithmic derivatives of D p (z) and Ai(z) satisfy and later we shall make use of these results.

Proof of Proposition 3
Here we establish the symmetry relations (11)-(13). These may be obtained without solving explicitly for p(x, t). Consider the problem in (2) with p = p(x, t; x 0 ; β, η) and set Then (2) becomes where R x ′ = ∂R/∂x ′ and the initial condition becomes where we have used the scaling law of the delta function.

Proof of Theorem 10
We letp(x; θ) = ∞ 0 e −θt p(x, t)dt and note thatp will be analytic in the right half-plane Re(θ) > 0. If p satisfies (2) its Laplace transform satisfies Assume that x 0 < 0 and x > 0, so that δ(x − x 0 ) = 0. By writingp = e −ηx 2 /4 e −βx/2 v(x; θ), (96) reduces to the differential equation where v ′ = dv/dx. This is the parabolic cylinder equation (Erdelyi [10], p. 116) and as we discussed in Section 5 two linearly independent solutions (at least for Re(θ) > 0) are given by as the second solution in (99) must be rejected due to its Gaussian growth (see (72)). For x < 0 we must use the second expression in (97), and then solve (98) with η = 1. Now we reject solutions with Gaussian growth as x → −∞, so that for x < x 0 the appropriate solution to the second equation in (97) (using (98) with η = 1) iŝ But in the range x 0 < x < 0 the solution will involve both of the parabolic cylinder functions The functions γ j (θ) are determined from continuity conditions at x = 0 and x = x 0 (cf. (4) and (5)). Continuity ofp and d dxp at x = 0 leads to Continuity ofp at x = x 0 yields and the jump condition of d dxp at x = x 0 , i.e., Here we used (101) to computep(x − 0 ), (102) to computep(x + 0 ), multiplied (106) by e  (107) give a 4 × 4 linear system for the γ j , whose solution leads to Theorem 10. We note that the calculations assumed that Re(θ) > 0. If Re(θ) ≤ 0 the 4 × 4 system may become singular, and in fact this occurs when θ = 0 and at the eigenvalues −θ = λ N , N ≥ 1. Theorem 10 thus gives the Laplace transformp for Re(θ) > 0, and then the expression can be analytically continued to the left half-plane, since we know how to continue the parabolic cylinder functions, which are entire functions of θ. After the continuation, locating the singularities in (28)- (30) in the range Re(θ) ≤ 0 can be used, for example, to obtain the spectral representation in (31).

General considerations for establishing Propositions 4-9
Here we discuss some general principles about solving V = 0 in Theorem 1, for the minimal root r(β, η), in the limit of η → 0 + . As discussed in Section 4, general results for Schrödinger equations and Sturm-Liouville problems show that the roots of V = 0 are all on the real axis, and that the sequence of roots (or eigenvalues) −θ N = λ N satisfies λ N → ∞ as N → ∞, for any fixed β and η > 0. Also, we know from Section 4 that the roots are all simple, and thus ∂V/∂θ = 0 when θ = −λ N .
Consider r as a function of β and η. Then V(−r(β, η); η, β) = 0 and by implicit differentiation we obtain ∂V ∂θ θ=−r · ∂r ∂β + ∂V ∂β = 0 (108) and since ∂V/∂θ| θ=−r = 0 we can use this relation to compute ∂r/∂β (for any β and any η > 0). By taking higher order derivatives of V = 0 with respect to β, a similar argument shows that r has derivatives of all orders with respect to β. Also, by differentiating V = 0 implicitly with respect to η, we conclude that r has derivatives of all orders with respect to η, for any η > 0. Thus r(β, η) is infinitely smooth for all real β and for η > 0, and this is true for the higher roots also. Note that since D p (z) is an entire function of both p and z (see [31]), V is an entire function of θ and β, and real analytic for η > 0. However, the limit η → 0 + is quite singular, as we shall show.
The above discussion shows that the roots of V = 0 vary smoothly with β and η, and a root cannot simply appear/disappear, say at some critical value η c . Thus for η → 0 + the roots have to lie in some range(s) of θ. In Section 5 we gave detailed results of the different asymptotic expansions of the parabolic cylinder functions D p (z), for different ranges of p, z. Applying these results to the equation V = 0, the function V can be approximated by simpler functions in the limit of η → 0 + , but these approximations are different in different ranges.
In what follows we shall use the following principle: suppose the equation F (u, ǫ) = 0 has roots u j = u j (ǫ) which depend on the small parameter ǫ, and these roots are smooth functions of ǫ. Also, suppose that F (u, ǫ) is an analytic function of both u and ǫ, with an expansion of the form Then if F 0 (u) has a simple root at u * then F (u, ǫ) has a root close to u * for ǫ → 0. The same conclusion holds if F is not analytic in ǫ, but has an asymptotic expansion of the form (109), where the expansion holds uniformly on some (finite) u interval that contains u * . In our case −θ plays the role of u and ǫ will correspond to η or a fractional power of η, such as √ η or η 1/3 .

Proof of Proposition 7
For β = O( √ η) and η → 0 + , we shall show that V = 0 has solutions in the range −θ = O(η). We use the facts that D 0 (0) = 1 and D ′ 0 (x) ∼ −x/2 as x → 0. We let η → ∞ and from (8) we obtain The error term in (110) is uniformly O(η −1 ) on finite θ intervals, and we note that the right-hand side of (110) is an entire function of θ, as will be the error terms. Then the symmetry relation for V in Proposition 3 implies that if we scale θ = ηS and β = γ √ η we obtain with an error term that is uniformly O(η) on finite intervals of S and γ. Proposition 7 follows upon setting S = −R and using the identity D 1−S (γ) + D ′ −S (γ) = 1 2 γD −S (γ) (see [16], p. 1066).
Then both θ * η M and θ * η M −1 are small and expanding the parabolic cylinder functions in (112) in Taylor series in both index and argument yields where the error term is o(1) for η → 0 + and this is uniform in finite θ * intervals. But then we conclude that θ * = 0, contradicting our assumption that there is a root in the range θ = Θ(η M ) for M > 1. To obtain the last expression in (113) we also used the recurrence (67). The calculation that led to (113) only used the fact that θ and θ/η are both small. Indeed, for any θ = o(η) we obtain (113) with θ * η M replaced by θ. Thus if there is a root in any range where θ = o(η), we again conclude that θ = 0. Hence, there can be no roots in ranges where θ = Θ( η log(1/η) ), θ = Θ( η log log(1/η) ), etc.
The error terms in (136) and (137) are uniform on finite δ intervals.

Monotonicity of the spectral gap
The surface sketched in Figure 1 suggested certain monotonicity properties of r(β, η), and these were partially confirmed by the various asymptotic formulas in Section 2. We now establish these analytically, for all values of (β, η). We shall obtain: Hence, for a fixed η < 1 the spectral gap r is an increasing function of β, while it decreases with β for fixed η > 1. If η = 1, r(β, 1) = 1 is constant.
To establish this result it is useful to set and then in view of Theorem 1 V (r(β, η); η, β) = 0.
By implicit differentiation of (151) we have By definition, r is the minimal positive solution of V (P ; η, β) = 0 and we also note that V (0; η, β) = 0, as θ = 0 is a simple pole of p(x; θ) in Theorem 10, which corresponds to the steady state limit in (6). Thus P = 0 and P = r are consecutive zeros of V (P ; η, β) = 0. To determine the sign of ∂r/∂β in (152) requires that we know the signs of ∂V /∂P and ∂V /∂β at P = r. For the former we can compute ∂V /∂P using the expression in Theorem 1 and the integral respresentations in (61), (64) and (65). However, an indirect argument leads immediately to the value of sgn( ∂V /∂P | P =r ). The solutions of V = 0 for P ≥ 0 correspond to poles of the Laplace transform p(x; θ) in Theorem 10 and these are the eigenvalues λ N for N ≥ 0, with λ 0 = 0 and λ 1 = r. From the general theory of the one-dimensional Schrödinger equation, the eigenvalues are all simple (see the discussion below (52)) and the equation V = 0 has simple zeros. Hence, ∂V /∂P | P =r = 0. We also note that if V had, say, a double zero at P = r, then (28) would imply that the spectral expansion of p(x, t) would involve the terms e −λ 1 t = e −rt and also te −rt , and this would contradict the selfadjointness of the Schrödinger equation in (49). Now, since P = 0 and P = r are consecutive simple zeros on the real axis of the entire function V (as a function of P ) we must have Computing the right-hand side of (153) is much easier than computing the left-hand side, as we show below. We define the functions I(z) and J(z) by Br u(log u)e −zu e u 2 /2 du.
Then we expand V in (150) in Taylor series about P = 0 to obtain .
But an integration by parts shows that Here we also used the fact that a parabolic cylinder function of order P = −1 can be expressed in terms of the standard error function, or probability integral. Using (157) in (156) we have so that ∂V /∂P | P =0 < 0 and hence, in view of (152) and (153), Now, Using the parabolic cylinder equation D ′′ P (z) = 1 4 z 2 − P − 1 2 D P (z) we can simplify (160) to When P = r we can further use the fact that V (r; η, β) = 0 to simplify the right side of (161). We will need to separately consider the two cases D r (−β) = 0 (a degenerate case that occurs rarely) and D r (−β) = 0 (which is typical).
In the degenerate case we have and thus We note that if D r (−β) = 0 then certainly D ′ r (−β) = 0, as discussed in Section 5 below (69). But if both V = 0 and D r (−β) = 0 then certainly D r/η (β/ √ η) = 0. In the non-degenerate case we can Using (164) in (161) and (159) we conclude that We proceed to determine the signs of the various terms in (163) and (165). It proves useful to understand the behaviors of r(β, η) as β → ±∞ for a fixed η. By a calculation completely analogous to that used to establish Proposition 4, we find that Whereas Proposition 4 applies for η → 0 with fixed β < 0, (166) applies for fixed η as β → −∞. By expanding (14) for β → −∞ and (166) as η → 0 we see that the two agree in this intermediate limit. Thus for β large and negative, r is exponentially close to η, as could be expected since then almost all of the probability mass in the model migrates to the range x > 0. A similar analysis as β → +∞ shows that which can also be obtained simply by using (166) and the symmetry relation in (13). For β large and positive the probability mass migrates to the region x < 0. Note that (166) and (167) suggest that ∂r/∂β has the oppositive sign as η − 1, at least for |β| sufficiently large, and this we proceed to establish for any β.
Note also that for β → −∞, (166) and (70) lead to and this verifies the conclusion about opposite signs. We have thus simplified (165) to in the non-degenerate case. In the degenerate case we conclude from (169) that D ′ rc (−β c ) and D ′ rc/η (β c / √ η) have the same sign, and then (163) shows that sgn(∂r/∂β) = − sgn(η − 1), which establishes Proposition 15. It remains to show that the last factor in (172) has always positive sign. Let us define and we consider H as a function of both P and z. We clearly have H(P, 0) = P D 2 P (0)+[D ′ P (0)] 2 > 0 for P > 0, with H(0, 0) = 0. Also, H(0, z) = 0 for all z, in view of (63). We consider P > 0 and z > 0.
We shall show that H(P, z) > 0 for all z ≥ 0 when P > 0. For z → ∞ the estimate in (89) leads to and then from (173) and (70) so that H is positive for z sufficiently large. By differentiating (173) with respect to z we obtain Here we also used (60) and the recurrence (67). From (176) we conclude that H has maximum or minimum values at roots of D P (z) and D P −1 (z), as functions of z. As discussed in Section 5, there are at most finitely many of these. But if D P (z * ) = 0 for some z * then D ′ P (z * ) = 0 and H(P, z * ) = [D ′ P (z * )] 2 > 0. If D P −1 ( z) = 0 for some z then (173) and (67) show that which is again positive for P > 0. Note that we cannot have simultaneously D P −1 ( z) = 0 = D P ( z), for then (67) would imply that D ′ P ( z) = 0 also. We have thus shown that H(P, z) is (for P > 0) positive at z = 0 and as z → +∞, and also H > 0 at any maximum/minimum value of H. We then conclude that H(P, z) > 0 for P > 0 and z ≥ 0.
Note that if H becomes negative at some z = z ′ then H would need to reach a minimum value at a point z ′′ where H < 0, since for sufficiently large z we again have H > 0. Now we let P = r(β, η) > 0 and z = −β and use (173) and (178) in (172) to conclude that sgn ∂r ∂β = − sgn(η − 1), β ≤ 0 and we have thus established Proposition 15 for β ≤ 0 and all η > 0. To show the result holds also for β > 0, we need only use the symmetry relation in (13), which shows that if r increases with β for β < 0 and 0 < η < 1, (resp. η > 1) then r will decrease with β for β > 0 and η > 1 (resp. 0 < η < 1). Alternately, we can use the relation (168) (in the non-degenerate case) and (172) to conclude that and apply (178) with P = r/η > 0 and z = β/ √ η > 0. This concludes the proof of Proposition 15, which was suggested by our numerical and asymptotic results.

B Appendix B
We discuss the singularities of (28)- (30) in the complex θ-plane and thus establish Proposition 11. As discussed in Section 5, D −θ (β) is an entire function of θ, so that the only singularities of (28) are the zeros of V(θ; η, β). The existence of an infinite sequence of zeros and the fact that they lie on the real axis (Im(θ) = 0) follows from standard ODE theory, which was discussed in Section 4. Now consider (29) and (30). The factor Γ(θ) has simple poles at θ = 0, −1, −2, −3, . . . . If θ = −M , M ≥ 0 we can simplify V by using (63), so that Similary, (27) and (63) lead to But the first equality in (185), along with the Wronkskian identity in (68), leads to .
we see that as η → 0 the expression in (28) becomes that in (39). We also have which can be used to obtain the limits of (29) and (30), and this agrees with the results we obtained in [23]. Throughout this calculation we divided several times by D −θ/η β/ √ η , which was permissible since, for θ > 0, we are outside of the oscillatory range of the special function, as we discussed in Section 5.

D Appendix D
Here we discuss the discrete M/M/m + M model. We shall obtain an explicit, albeit complicated, expression for the Laplace transform (over time) of p n (t) = Prob[N (t) = n | N (0) = n 0 ], where N (t) is the number of customers in the system. Then we will give an alternate derivation of Theorem 1, by evaluating the discrete model in the limit m → ∞ with ρ = λ/µ = m + O( √ m). The analysis here closely parallels the proof of Theorem 10, so we just give the main points.
We solve the following infinite system of ODEs (we assume time is scaled to make the service rate µ = 1, so that ρ = λ): − [m + (n − m)η] p n (t), n ≥ m + 1 with the initial condition p n (0) = δ(n, n 0 ). We need to consider the cases n 0 < m and n 0 > m separately, as the discrete model has no analog of the symmetry relation in Proposition 3.
Here we used the integral representation in (61) for the parabolic cylinder function D. In (212) Br is a vertical contour with Re(ξ) > 0, which can be used to approximate C 0 in (199) with this scaling of z. A completely analogous expansion of (200) leads to and, after some calculation, we obtain from (201) Using (201) and Stirling's formula we also have Next we consider the limiting form of (207). Noting that F m − F m−1 can be computed from (199) by multiplying the integrand by 1 − z and setting n = m, we have Then using (212)-(215) and (217), we see that the expressions in (203)-(205) reduce to those in (28)-(30), up to a factor of 1/ √ m in the former, which arises due to the fact that p n (t) is normalized by a sum over n while p(x, t) is normalized by an integral over x. We have thus given an alternate derivation of Theorem 10. Note that all of the asymptotic calculations do not involve scaling time t or the transform variable θ.
Finally, we discuss the uniformity of the approximation in (216) for small values of θ. This will show that ∆(θ) = ∆(θ; m, β, η) ≡ F m H m−1 − F m−1 H m can have no roots for m → ∞ in ranges where θ = o(1). We set ρ = m − β √ m and consider finite intervals of β and η, with η > 0. First we note that, using (199)  The pole at θ = 0 corresponds to the steady state limit, which exists for all m and β, for η > 0. We expand ∆ for fixed finite β and fixed η > 0, as m → ∞, which will refine the leading order result in (216) and show that the higher order terms remain smaller than the leading term, for θ = o(1) as m → ∞. Setting z = 1 − ξ/ √ m in the integral in (199) leads to m −θ/2 F m (θ) = e ρ √ m C ′ e ξ 2 /2 e βξ ξ −θ F(ξ; m)dξ, where C ′ is the image of the contour C 0 and where P j (ξ) is a polynomial in ξ of degree 3j. We have P 1 (ξ) = ξ + 1 3 ξ 3 , P 2 (ξ) = ξ 2 + 7 12 ξ 4 + 1 18 ξ 6 , etc. Using (220) and some contour deformation (as |z| < 1 in (199) implies asymptotically that Re(ξ) > 0), we obtain the asymptotic series where f j (θ, β) = 1 2πi Br + P j (ξ)e ξ 2 /2 e βξ ξ −θ dξ and Re(ξ) > 0 on the vertical contour Br + . Since the P j are polynomials, the integral in (222) is a finite sum of parabolic cylinder functions of different orders (for example f 1 involves D 1−θ (·) and D 3−θ (·)), and thus each f j is an entire function of θ. A completely analogous calculation shows that Next consider H m (θ) (setting n = m in (201)) in the same asymptotic limit. Now we scale z = 1 + ξ η/m and obtain whereP j are again polynomials in ξ. Then we obtain the asymptotic expansion of (224) as where h j (θ, β, η) = 1 2πi Br +P j (ξ, η)e ξ 2 /2 e −βξ/ √ η ξ −θ/η dξ where Here V 0 = V and f 0 and h 0 can be identified from (221) and (226). For (230) we see that each V j is a finite sum of products of two parabolic cylinder functions, possibly of different indices and arguments. Thus each V j is an entire function of θ, and the expansion in (229) is uniform in finite β, θ, η intervals, with η > 0. Setting θ = 0 and using the fact that ∆(θ), which is itself and entire function of θ, vanishes for all values of m, β, η, we conclude that V j (0; η, β) = 0 for each j. From (158), V has a simple zero at θ = 0 and then V j for j ≥ 1 have zeros of orders ≥ 1. Thus for any θ = o(1) as m → ∞ the correction terms in (229) remain smaller than the leading term, and thus ∆(θ) cannot have a zero for θ = o(1), except the one at θ = 0. We have thus shown that for fixed β, fixed η > 0, and m → ∞ the zeros of ∆ in (216) in the range θ = Θ(1) must approach the zeros of V, and that ∆ has no zeros where θ = o(1), except for θ = 0.