Scaling limits for infinite-server systems in a random environment

This paper studies the effect of an overdispersed arrival process on the performance of an infinite-server system. In our setup, a random environment is modeled by drawing an arrival rate $\Lambda$ from a given distribution every $\Delta$ time units, yielding an i.i.d. sequence of arrival rates $\Lambda_1,\Lambda_2, \ldots$. Applying a martingale central limit theorem, we obtain a functional central limit theorem for the scaled queue length process. We proceed to large deviations and derive the logarithmic asymptotics of the queue length's tail probabilities. As it turns out, in a rapidly changing environment (i.e., $\Delta$ is small relative to $\Lambda$) the overdispersion of the arrival process hardly affects system behavior, whereas in a slowly changing random environment it is fundamentally different; this general finding applies to both the central limit and the large deviations regime. We extend our results to the setting where each arrival creates a job in multiple infinite-server queues.


Introduction
Empirical studies show that the number of arrivals in customer contact centers, hospital emergency departments and cloud computing systems typically varies strongly over time [8,16].This motivates modeling such arrival processes by a non-homogeneous Poisson process (NHPP) with time-dependent arrival rate λ(t), see e.g.[9].At the same time, various studies show that in a broad variety of real-life systems the intensity of the fluctuations in the arrival rate is so severe that the Poisson assumption ceases to hold [2,8].The observed level of overdispersion urges the need to develop stochastic models that can capture such persistent fluctuations.
Starting from the classical Poisson process, it is common practice to increase dispersion by using a mixed Poisson process [2,13], to that end replacing a deterministic parameter λ by a random parameter Λ.This leads to the idea of modeling overdispersed arrival processes by a mixed version of NHPPs, so-called Cox processes [5], where the time-dependent rate λ(t) of the classical NHPP is replaced by a stochastic process Λ(t).For instance, one could use Markov-modulated Poisson processes in which the arrival rate Λ(t) = λ J(t) is a function of a continuous-time Markov chain J(•) on a finite state space S and non-negative rates λ i for i ∈ S (see e.g.[1,3]).To also include, say, diurnal patterns, one could work with the arrival rate Λ(t) + λ(t) for some function λ(t).Although the Markov-modulated Poisson process is versatile and has various attractive properties, it has considerable drawbacks as well.First, while a substantial body of results for single-server queues with Markov modulation has been established, considerably less is known about their many-server and infinite-server counterparts; see e.g. an account of this issue for the infinite-server system in [4].Second, due to the fact that the process J(•) is not observed, estimating the parameters of a Markovmodulated Poisson process from data is a non-trivial task [14].
The main objective of this paper is to develop an arrival process simpler than a Markovmodulated Poisson process -arguably the simplest in terms of analysis -that fits the overdispersed and time-dependent setting, and to assess the impact of these characteristics on a corresponding system's performance.The model we propose is a mixed Poisson arrival process in a random environment.It is defined as follows.Let Λ a non-negative random variable with finite first two moments and density f Λ (•).Introduce a sampling frequency 1  ∆ ; then the arrival rate at time t is given by Λ j when t ∈ [j∆, (j + 1)∆), where the Λ j are independent and distributed as a non-negative random variable Λ, for j ∈ Z.In other words, this arrival process is a special case of a stationary Cox process where the arrival rate at time t is given by Λ(t) = j Λ j 1 [j∆,(j+1)∆) (t). (1.1) To add nonstationarity in the arrivals, one could include a deterministic component λ(t) without intrinsically complicating the analysis; for ease of presentation we omit the extra component here.The resulting process can be viewed as an extension of the classical mixed Poisson setting, which is enriched by (independently) resampling the arrival rate after every time slot of length ∆ > 0. The intuition is that the arrival rate changes every ∆ time units, so that the number over a large time slot fluctuates more severely than standard Poisson data would, as can be made explicit via an elementary computation.Let the number of arrivals up to time t be given by N t ∼ Pois( t 0 Λ s ds) and let t be some multiple of ∆ (for simplicity).Then EN t = tEΛ, whereas Conclude that, as desired, the variance-to-mean ratio is strictly larger than 1 for non-deterministic Λ, i.e., Var(N t ) EN t = 1 + ∆ Var(Λ) EΛ .
Observe that the level of overdispersion is determined by the interval length ∆ and the level of overdispersion in Λ (through its variance-to-mean ratio).Given this model for the arrival process, various queueing models can be studied; in this paper we focus on single-class infinite-server systems with exponential service times.The proposed arrival process being overdispersed, the main objective of this paper is to reveal, in a compact manner, the impact of overdispersion on system performance.Infinite-server systems are a natural choice when the system at hand is designed to (almost) immediately serve all customers [15], but it may also serve as a tractable proxy for the more complicated multi-server systems, which is for instance exploited in the modified offered-load (MOL) and pointwise stationary approximation (PSA) methods for staffing large-scale service systems in a time-varying setting [10,16].
Contributions.Infinite-server systems with overdispersed arrivals are, as described above, very tractable.As shown in Section 2, it is fairly straightforward to compute the probability generating function (pgf) of the stationary and time-dependent queue length processes in terms of transforms.This is due to the fact that customers are served immediately upon arrival, independently of each other; as a result, when analyzing the queue length at a given point in time, we can separately consider the individual (independent!)contributions that correspond to each of the preceding intervals of length ∆.
The queue length distribution can be characterized in terms of its pgf, which effectively means that evaluation of the accompanying performance measures requires numerical inversion.However, by imposing a scaling on both the time and scale parameters, ∆ and Λ, we succeed in identifying an asymptotic regime in which the distribution can be explicitly given.We inflate the arrival rate and sampling frequency in the following way: where we let N → ∞.Importantly, Λ and ∆ −1 do not necessarily grow at the same rate under scaling (1.2).The value of α determines the asymptotic behavior of the resulting scaled system, giving rise to a trichotomy.For α > 1, in which case the arrival rate is resampled relatively frequently, we find that the system behaves as a standard infinite-server queue (no overdispersion), whereas for α < 1 the overdispersion remains present in the asymptotic regime.
The case α = 1 essentially reflects a superposition of the two distinct types of behavior.
For preparatory purposes, we show in Section 2 that the centered and normalized stationary queue length is asymptotically normal under the scaling in (1.2).Next, in Section 3 we consider a multidimensional setting with correlated arrivals: an arrival triggers jobs in multiple queues.Hence, we work with a coupled system in which d parallel queues are fed by a single arrival process; cf.[11,12].With U (N ) (•) denoting the vector of centered and normalized queue length processes, the asymptotic normality now translates to the corresponding limiting process U (•) being Gaussian: U (•) is a d-dimensional process of the Ornstein-Uhlenbeck type with parameters that depend on the scaling regime.Following the approach in [1], we show this by applying a lemma due to Kurtz and a martingale central limit theorem (mclt) to a suitable stochastic integral equation.
Subsequently, in Section 4 we carry out a large deviations analysis to obtain the logarithmic tail asymptotics corresponding to the queue length distribution.The crucial observation in this analysis is that rare events can essentially be realized in two ways: (i) the random arrival rate attains an exceptionally high value, (ii) the Poisson process generates an unusually large number of arrivals given the (not so rare) value of the random parameter.Again, the value of α determines what type of tail behavior dominates: for α < 1 this is effect (i), for α > 1 effect (ii), and for α = 1 a combination of effects (i) and (ii).These findings complement similar results that have been established for an infinite-server system with Markov-modulated input, where it is noted that the slow regime (α ∈ (0, 1)) was not covered in that setting [3,6].We conclude Section 4 by pointing out how the large deviations results can be extended to the multidimensional setting.

Overdispersion in an infinite-server context
In this section we present a stationary and transient analysis of the single-class Markovian infinite-server system in a random environment just introduced.A crucial role is played by Λ(t), the arrival rate at time t given in (1.1).Remember that we assumed that the arrival rates are i.i.d. and distributed as a random variable Λ 0 with finite first two moments and density f Λ (•).The corresponding service times are assumed i.i.d.(and in addition independent of the arrival process) exponentially distributed random variables with mean 1/µ.
First, in Section 2.1, we analyze the stationary system behavior, in terms of its pgf and the corresponding moments, which we then extend to the associated transient behavior.We then study the stationary behavior in a central limit regime under parameter scaling (1.2) in Section 2.2.This exposition serves as an illustration for the reader, and is intended to create intuition as for why the scaled stationary queue length is asymptotically normal and why the three different limiting regimes appear; in addition, in Section 4 we need a result that is proven along the same lines.We remark that in Section 3 the normality is generalized in several directions: we establish a functional central limit theorem (fclt) for the (scaled) transient process M (N ) (•) corresponding to the d-dimensional parallel system as defined in the introduction.

Pre-limit results
This subsection presents 'pre-limit results'; later we study their counterparts in the limiting regime after imposing a parameter scaling.
Transform of stationary queue length.Let M be the random variable associated with the stationary number of jobs (also sometimes referred to as 'customers') in the system.Exploiting 'thinning' properties, we can identify the pgf φ(z) := Ez M of M.
In the sequel we write p t := e −µt for the probability that a job present at kt is still present at (k + 1)t and q t := (1 − e −µt )/(µt) for the probability that a job arriving at a uniform epoch in [kt, (k + 1)t) is still present at (k + 1)t.Denote pt := 1 − p t .
Note that M can be written as the sum of M 0 , M 1 , M 2 , . .., where M k represents the number of jobs that arrived in [−(k + 1)∆, −k∆) and are still present at time 0. Furthermore, observe that these 'thinned' random variables M k are independent.A job that arbitrarily arrived in [−(k + 1)∆, −k∆) (i.e., having arrived at a uniform epoch in this interval) is still in the system at time 0 with probability As a consequence, with r t := tq t , Observe that φ k (z) is a pgf of 'mixed Poisson' type: conditional on Λ k = λ the pgf corresponds with that of a Poisson random variable with mean λr ∆ p k ∆ .We conclude that M k is distributed as a mixed Poisson random variable with random parameter with Λ k the value of the arrival rate in the interval [−(k + 1)∆, −k∆) (note that, in fact, we should have written Λ −(k+1) rather than Λ k , but due to the i.i.d.assumption the processes {Λ(s)} s 0 and {Λ(−s)} s 0 have the same finite-dimensional distributions).Therefore, M is mixed Poisson as well and its random parameter is given by (Note that κ(•) is defined as a functional; κ(Λ) should be interpreted as κ(Λ(•)).) There is an alternative way to obtain this result.Indeed, since we observe the system in stationarity, where g Λ,t (z) is defined by Applying an iteration argument to (2.3) yields (2.4) In the factors g Λ,∆ 1 − (1 − z)p k ∆ we recognize the expression for φ k (z) as in (2.1).
First two moments.We now evaluate the first two moments of M. This is an interesting computation in its own right, but it also provides useful results that can be exploited when considering this system under the central limit scaling (as is done in the next subsection).Differentiating (2.3) and letting z ↑ 1, we obtain a fixed-point equation, Hence EM = φ ′ (1) = r ∆ EΛ/(1 − e −µ∆ ) = EΛ/µ.This quantity could have been computed more directly as well, using a standard identity for conditional means: Then observe that (M k | Λ k ) is Poisson, and hence its mean equals its parameter.As a result, (2.5) equals For the variance we use that , and hence It thus follows that, after some algebra, where Alternatively, one could use the 'law of total variance' to identify VarM: Observe that, because of the 'mixed Poisson property', and as a result the first term at the right-hand side of (2.7) equals EM.The second term, which is inherently non-negative, gives rise to 'overdispersion', i.e., the effect that the variance of the stationary queue length exceeds the corresponding mean.This is a distinguishing feature compared to the analogous system in which the Poissonian arrival rate is deterministic: the stationary queue length in an M/M/∞ system is Poisson, and cannot accommodate any overdispersion.In order to evaluate the second term in the right-hand side of (2.7), we note that Substituting (2.8) in the second term in the right-hand side of (2.7), we find that Var M equals (2.6), as desired.Formula (2.6) lends itself to a nice interpretation: the term EΛ/µ is the contribution to the variance that one would have if the arrival rate would have had the deterministic value EΛ, whereas the term C VarΛ/µ 2 needs to be added in order to deal with the non-Poisson variability due to the stochasticity of the arrival rate.
Transient behavior.As the analysis of the transient system behavior strongly resembles its stationary counterpart, we restrict ourselves to a short account of this.We let the system start empty (for ease of presentation; a non-empty initial condition can be analyzed without any additional difficulty).Denote by M(t) the number of jobs present at time t.Then, for n the smallest integer such that t − n∆ < ∆, where Mj ( M[n∆,t) ) represents the number of jobs that have arrived between in [j∆, (j + 1)∆) ([n∆, t)) and are still around at n∆ (t).As before, these have pgfs As the individual random variables M1 , M2 , . . .and M[n∆,t) are independent, M(t) is mixed Poisson with random parameter κ t (Λ) := t 0 Λ(s)e −µs ds. (2.9)

Limit results
This section focuses on the central limit regime that results from simultaneously scaling, in a controlled way, both the arrival rate Λ and the sampling frequency 1 ∆ as in (1.2).Let the scaled counterpart of Λ(t) be NΛ (N ) (t), with (2.10) That is, the sampling frequency and the arrival rates are both inflated as we let N tend to ∞, but, importantly, at rates that are not necessarily identical.As mentioned in the introduction, depending on the value of α, we obtain fundamentally different behavior.
We consider a sequence of systems indexed by N, where the N-scaled system uses a mixed Poisson arrival process with time-dependent random rate NΛ (N ) (t).Let M (N ) denote the stationary queue length in the N-scaled system, with parameter Nκ(Λ (N ) ) (cf. (2.2)).We start our exposition by a preliminary calculation, in which we compute the mean and variance of M (N ) and study their behavior for large N, which indeed reveals the announced trichotomy.Then, after centering and normalizing M (N ) , we derive a central limit theorem.
Qualitative behavior of first two moments: trichotomy in variance.First, we identify the steady-state mean and variance in our scaling regime, using (2.5), (2.7) and (2.8).We find that (2.11) where it is noted that for large N, (2.12) behaves approximately as (the ratio of the two converges to 1).We thus observe the trichotomy (2.13) For α > 1, the sampling frequency dominates the variability of Λ.Consequently, the model behaves essentially as an M/M/∞ system, with the variance of M (N ) being linear in N and equal to EM (N ) , for large N.For α < 1, we find a superlinear relation between N and Var Λ, and both the sampling frequency (i.e., the reciprocal of the interval length ∆) and the variance of Λ play a role.Hence, the asymptotic variance indeed grows faster than the asymptotic mean for α < 1; in this regime the system is overdispersed.For α = 1, the variance is 'slightly larger' than for α > 1, but it is still linear in N. In this case the sampling frequency and the variance of Λ grow at the same rate, so that the variance for M (N ) combines the effects observed in the two former cases.
As observed from the above computation, the variance of M (N ) is essentially proportional to N γ with γ := max{1, 2 − α}.As a consequence, one may expect that, under (1.2), converges to a (zero-mean) normally distributed random variable.It is this property that we verify now.
Asymptotic normality.We show how to establish asymptotic normality for the centered and normalized version of M (N ) in (2.14) via evaluation of the corresponding Laplace transform.
Appealing to Lévy's convergence theorem, we establish the desired convergence in distribution.
For simplicity, the proof of Thm.2.1 assumes that all moments of Λ are finite; however, as will appear from the proof of Thm.3.2 only finiteness of the first two moments is necessary.
Theorem 2.1 (clt).As N → ∞, M(N) converges to a zero-mean normally distributed random variable with variance Proof.Let φ (N ) (z) be the counterpart of (2.4) under scaling as in (1.2); likewise g We are interested in the behavior of M (N ) in the central limit regime, hence we need to analyze the limiting distribution of M(N) .To this end, we evaluate the logarithm of the corresponding Laplace transform: (2.15) We now use that log Observe that (2.16) is the sum of cumulant generating functions (which is again a cumulant generating function), each of them related to the random variable Λ but evaluated at different arguments.Let m ℓ denote the ℓ-th cumulant of Λ (for ℓ ∈ N); in particular m 1 = EΛ and m 2 = Var Λ.In addition, we define Then it follows that (2.17) Let us first consider the contribution of the term corresponding to ℓ = 1.Observe that, as Note that the first term in the right-hand side of (2.18) is canceled by the first term in the right-hand side of (2.15), so that we are left with the second term, i.e., The second term in (2.17) corresponding to ℓ = 2 gives  1  2 s 2 ∆ VarΛ/(2µ).Finally, if α = 1, we find that both terms converge to the expected finite positive limit.
We have therefore established that, as as claimed.

Functional central limit theorem
In this section we generalize the central limit result of Thm.2.1 in two ways.First, we establish the functional version: the centered and normalized transient queue length process converges to a limiting process of Ornstein-Uhlenbeck type with parameters that depend on the value of α.Second, we extend this to a multidimensional setting with correlated arrivals: every arrival triggers jobs in multiple queues.The correlation structure of the resulting multidimensional Gaussian limiting process is explicitly identified.
Let us start by describing the mechanics of the generalized setting.We consider a parallel system in which d queues are fed by a single arrival process that was constructed in the same way as the one in the previous section: a Markovian process with arrival rate Λ(t) as in (1.1).The service times in queue i are i.i.d.exponential random variables with mean µ −1 i ; the service processes of the individual queues are independent, and also independent of the arrival process.We perform the same scaling as before: the sampling frequency is sped up by a factor N α , while the (random) arrival rate is blown up by a factor N. This results in a mixed Poisson arrival process with time-dependent rate Λ (N ) (t) as in (2.10).Let where M (N ) i (t) is the queue length at time t in the i-th queue of the N-scaled system, for i ∈ {1, . . ., d}.Note that the M (N ) i (t) are mixed Poisson with time-dependent random parameter Nκ t,i (Λ (N ) ), with κ t,i (Λ (N ) ) as defined in (2.9) but with µ replaced by µ i .
We now present an alternative way of writing M (N ) i (t), which facilitates the use of a martingale central limit theorem.We introduce the functional ](t) is to be interpreted as the 'cumulative service capacity' in queue i over the interval [0, t].In addition, for the 'cumulative arrival rate' we have Ψ[Λ](t), with scaled counterpart NΨ[Λ (N ) ](t).By the law of large numbers, Ψ[Λ](t)/t converges a.s. to EΛ as t → ∞ and for fixed t, Ψ[Λ (N ) ](t) converges a.s. to t EΛ as N → ∞.
With Y 0 (•), . . ., Y d (•) denoting independent unit-rate Poisson processes, Our objective is to derive a d-dimensional fclt for M (N ) (•).This result characterizes the timedependent queue length in the scaled system and makes explicit how the correlated arrivals lead to correlation between the individual queue length processes.It will be stated and proved in subsection 3.2; first we study the stationary behavior by presenting the corresponding first two moments of M (N ) (including covariances between the individual queues).

Qualitative behavior of first two moments in stationarity
Note that the individual queue lengths are only coupled through the arrival process, so under (1.2), the mean and variance of M (N ) are, as in (2.11) and (2.12), given by for i = 1, . . ., d.Hence, we find the same behavior as in (2.13).Interestingly, the same trichotomy is observed for the covariances, as stated in the next lemma.
Lemma 3.1 (Covariance in M (N ) ).For i, k ∈ {1, . . ., d} with i = k, and for large N, Proof.Without loss of generality, we take i = 1 and k = 2.We first consider the non-scaled model, by studying the joint probability generating function, where ξ jn (z 1 , z 2 ) is the contribution due to the slot between j∆ and (j + 1)∆; as z 1 and z 2 are held fixed for the moment, we suppress them.Now we introduce functions (for ℓ = 1, 2) where it is noted that g j (µ, n) behaves as e −µ(n−j)∆ for small ∆.In addition, we define the quantities Using arguments similar to those we have used before, Because the contributions to M 1 (n∆) and M 2 (n∆) resulting from different time intervals are independent, we obtain that . Now imposing scaling (1.2) and considering the stationary behavior by letting n → ∞, it is readily derived that (for large N) observe that, for reasons of symmetry, it is allowed to replace n − j by j in the definition of the g j (µ, n).We thus arrive at which behaves in accordance with (3.2) for N large.
Recall that γ = max{1, 2 − α}; the above computation shows that the covariance matrix of M (N ) is essentially proportional to N γ .Therefore, we expect that the centered and normalized version of the joint stationary queue length process converges to a (zero-mean) d-dimensional Gaussian random vector with covariance matrix C such that for i, k ∈ {1, . . ., d}.This is verified in the next subsection.

Proof of functional central limit theorem based on mclt
The main objective of this subsection is to derive a functional limit theorem for M (N ) (t), the vector describing the queue lengths of the scaled system at time t.To this end, we consider the process M(N) We will need the following lemma, which uses the law of large numbers for Poisson processes; see [1].
Lemma 3.2.Let Y be a unit-rate Poisson process.Then for any U > 0, almost surely The uniform convergence in Lemma 3.2 entails that (3.3) converges almost surely to the solution of the functional equation as N → ∞, under the proviso that M(N) i (0) converges a.s. to some value ̺ i,0 for i = 1, . . ., d.The solution is given by a convex mixture of the initial position ̺ i (0) = ̺ i,0 and the limiting value EΛ/µ i : Having identified this fluid limit, the next objective is to establish an fclt for the centered and normalized process U (N ) (•) given by with β := 2 − γ = min{1, α}.Here we closely follow the approach in [1], where the idea is to use an mclt, so as to obtain weak convergence to a (generalized) Ornstein-Uhlenbeck process.
The version of the mclt that we need in our setting is stated below.
Theorem 3.1 (mclt, [1]).Let {R (N ) } N ∈N be a sequence of R d -valued martingales.Assume that the following condition on the jump sizes is met: in addition, assume that, as N → ∞, for a deterministic function C ik (t), continuous in t for all t > 0 and for i, k = 1, . . ., d.
Stated below is the main theorem of this section: an fclt for U (N ) (•), the process defined via (3.6).In line with earlier findings, three regimes need to be distinguished: α > 1 (the fast regime), α < 1 (the slow regime) and α = 1 (the intermediate regime).Theorem 3.2 (fclt).As N → ∞, U (N ) (•) converges weakly to a zero-mean d-dimensional Gaussian process with covariance matrix given by ) Proof.Using (3.3), we write for i = 1, . . ., d. Adding and subtracting ̺ i,0 , (3.12) is equivalent to which, by filling out the implicit form of ̺ i (t) as in (3.4), simplifies to We consider the three individual components separately.
(i) Component U (N ) 0,i (t) consists of the starting value of the process, which is assumed to converge to some value U i (0), minus a reverting term.It is now straightforward that, as (ii) Then consider U (N ) 1 (•).For α ≤ 1 (and hence β 2 = α 2 ), the standard functional central limit theorem for partial sums of i.i.d.random variables entails that, as N → ∞, with V (•) a standard Brownian motion.On the other hand, for α > 1 the limiting process is identical to 0, as a consequence of (iii) Finally, from Lemma 3.3, we conclude that U (N ) 2 (•) converges weakly to a d-dimensional zero-mean Brownian motion with covariance matrix K 0 (t) + K(t) for α 1, and to 0 else.
Using the above observations, we can now complete the proof.Each of the three regimes will be considered separately.
1. Fast regime (α > 1).We have obtained above that U (N ) (•) converges weakly to the solution U (•) of the d-dimensional stochastic integral equation given by with W0 (•), W1 (•), . . ., Wd (•) independent standard Brownian motions.It takes a routine calculation to derive that All linear combinations of the U i (•) are Gaussian processes, so we conclude that this ddimensional process is Gaussian.It is readily seen that E U i (t) = U i (0)e −µ i t .For the variance, an elementary computation gives Likewise, for the covariance, with This shows (3.10) for α > 1.
2. Slow regime (α < 1).In the slow regime, U (N ) (•) converges to the solution of which can be written as Therefore the U i (•) are Ornstein-Uhlenbeck processes given by: As before, we can conclude that this d-dimensional process is Gaussian with expectation vector given by U i (0)e −µ i t .Computations as above reveal that for α < 1, as claimed in (3.10), 3. Intermediate regime (α = 1).In this regime, a combination of the processes from the other cases appears: The marginal solutions U i (t) are, for i = 1, . . ., d, equal to Again, we conclude that this d-dimensional process is Gaussian with expectation vector given by U i (0)e −µ i t ; routine computations yield the desired covariance matrix, as given in (3.10) and (3.11).This completes the proof.
It is interesting to study the impact of the scaling parameter α on the correlation between the individual queue lengths.Remarkably, it turns out that for α = 1 this correlation depends on the service rates only, whereas for α = 1 also the first and second moment of Λ play a role; see the following corollary for a result on the stationary regime.
Copying the approach -using cumulant generating functions -underlying the proof of Theorem 2.1, it is readily derived that Hence, lim inf Now it remains to show that k t (Λ (N ) ) can again be replaced by κ t (Λ (N ) ) (which we abbreviate for compactness to k t and κ t ).Note that, as N → ∞, The rate function in (4.6) is continuous in η, so now letting η ↓ 0 yields (4.2).
As in the central limit regime, we distinguish between the cases α > 1, α = 1, and α < 1.For all three cases we derive the logarithmic asymptotics.
1. Fast regime (α > 1).We can bound (4.1) from below by for some ε ∈ (0, a − ρ(t)).As N tends to infinity, it is directly shown that the second factor in (4.7) converges to 1, and hence has exponential decay rate 0. Now an application of Cramér's theorem [7] yields lim inf On the other hand, (4.1) is majorized by As this holds for all ε > 0, we conclude that lim Recognizing the decay rate of a Poisson distribution with mean ρ(t), we observe that the essential behavior in the fast regime is again of M/M/∞ type.
2. Slow regime (α < 1).In this regime we need to distinguish between the situation in which the random variable Λ almost surely results in a κ t (Λ (N ) ) below a, and the situation in which this is not the case.The proof of the following lemma is straightforward hence omitted.
Lemma 4.2.Given Λ, let y = inf{x > 0 : The cases u(t) a and u(t) < a should be treated differently, as follows from the following intuitive explanation that is based on the decomposition (4.1).If u(t) a, then the random variable Λ can be 'large' with respect to a, which enables M (N ) (t) to reach Na without the Poisson random variable attaining an unlikely value.If on the contrary u(t) < a, then Λ is 'small' with respect to a; M (N ) (t) can only exceed level Na by the Poisson random variable attaining an extraordinarily large value.
We first consider the case u(t) < a.For ease we assume that Λ attains values in a discrete set of positive values, of which y is the largest (occurs with probability p ∈ (0, 1)) and y ′ < y the one-but-largest.It is directly seen that, for θ > 0, As α < 1, this leads to lim In addition, Ee θM (N) (t) is majorized by which converges to the right-hand side of (4.9) on an exponential scale (use y > y ′ ).Applying 'Cramér', we thus find that the probability of interest decays exponentially: whereas for α = 1 it turns out to equal e −µ i s (e θ i − 1) + 1 − 1 ds.
For α < 1, as before, the decay is either exponential in N (if the the multi-dimensional random Poisson parameter cannot attain values that are contained in A), or exponential in N α .The latter regime being the more complicated one, we here include the corresponding decay rate.
The probability of our interest can be rewritten as θ i e −µ i s ds .

Discussion and future research
In this paper we propose to model an overdispersed arrival process by a mixed Poisson process in a random environment.We assess the impact of overdispersion on system performance when feeding such an arrival process into an infinite-server system.Under a specific scaling, we derive (functional) central limit results and large deviations asymptotics.Various extensions can be explored, a few of which are mentioned here.To start with, many results seem to carry over to the setting with generally distributed service times.In addition, systematically studying the effect of adding a deterministic trend λ(•) to the random environment Λ(•), the results could be generalized to a setting with nonstationary Cox arrival processes.Another challenge lies in refining the logarithmic asymptotics, as obtained in Section 4, to exact asymptotics.
In all of the results obtained, we revealed a trichotomy in system behavior depending on the imposed scaling on system size and sampling frequency.Here the scaling primarily serves to change the level of overdispersion in the system.The combination of tunable sampling and tunable overdispersion provides a rich framework for modeling real-world arrival processes.One could imagine that in a rapidly changing environment, the inherent overdispersion of the arrival process hardly plays a role, whereas in a slowly changing random environment, overdispersion is expected to be more dominant.This interplay between sampling frequency and overdispersion is a convenient feature of our model, which could be used to calibrate the model to real-world data.The latter could be a promising direction for future research, involving challenging statistical issues.
Another application of our model would be in the area of dimensioning service systems or staffing, and in particular square-root staffing in many-server systems.The general idea behind square-root staffing is as follows: a finite-server system is modeled as a system in heavy traffic, where the number of servers s is large and at the same time, the system is critically loaded.Under Markovian assumptions this can be achieved by setting s = λ + β √ λ (denoting the load on the system by λ) and letting λ → ∞ while keeping β > 0 fixed.The system then reaches the desirable Quality-and-Efficiency-Driven (QED) regime, in which the system load approaches 100% while the delays experienced by customers remain limited.In such large-scale service operations, it is natural to use an infinite-server system as a proxy to the many-server system.Infinite-server models are extremely useful because of their tractability; this can be exploited by translating detailed knowledge of the infinite-server system state to the finite-server setting.This returns rather good estimates of future arrivals, even in situations of time-varying arrival processes [15,16].The model developed in this paper provides a new way of modeling such large-scale service systems, with the additional feature of a tunable level of overdispersion, essentially replacing a deterministic λ by a stochastically fluctuating Λ.The possibility to design asymptotic dimensioning schemes compatible with our new model -for both static and time-varying overdispersed arrival processes -is currently investigated by the authors.

Proof.
We start by checking the conditions of Thm.3.1.First, observe that for each N, Y (N ) 0 (•) and Y (N ) (•) are d-dimensional real-valued martingales.Also, condition (3.7) is met, as both for R (N ) = Y (N ) 0 and R