Stationary Distributions and Convergence for M/M/1 Queues in Interactive Random Environment

TL;DR

This paper analyzes a Markovian single-server queue in an interactive random environment, deriving explicit stationary distributions and convergence rates for two types of environments, enhancing understanding of their long-term behavior.

Contribution

It provides explicit stationary distributions and convergence rate estimates for M/M/1 queues in interactive random environments, including jump and jump-diffusion types.

Findings

Stationary distribution is explicitly derived as a weighted geometric distribution.

Exponential convergence rates to stationarity are explicitly estimated.

Results apply to both pure jump and reflected jump-diffusion environments.

Abstract

A Markovian single-server queue is studied in an interactive random environment. The arrival and service rates of the queue depend on the environment, while the transition dynamics of the random environment depends on the queue length. We consider in detail two types of Markov random environments: a pure jump process and a reflected jump-diffusion. In both cases, the joint dynamics is constructed so that the stationary distribution can be explicitly found in a simple form (weighted geometric). We also derive an explicit estimate for exponential rate of convergence to the stationary distribution via coupling.

Equations

$$\mathcal{L}f(n,z)={\mathcal{M}}_{z}f(n,z)+{\mathcal{A}}_{n}f(n,z),$$

$${\mathcal{M}}_{z}f(n,z)=\lambda (z)(f(n+1,z)-f(n,z))+1_{\{n\ne 0\}}\mu (z)(f(n-1,z)-f(n,z)).$$

$$\mathcal{A}_{n}f(n,z)=b_{n}(z)\cdot\nabla_{z}f(n,z)+\frac{1}{2}\operatorname{tr}\big{(}\Sigma_{n}(z)\nabla^{2}_{z}f(n,z)\big{)}$$

$${\mathcal{A}}_{n}f(n,z)={\beta }_n{\rho }^{-n}(z)\mathcal{A}f(n,z)$$

$$\begin{array}{rl}R[(n,z),(n+1,z)]& =\lambda (z),R[(n,z),(n-1,z)]=\mu (z),\\ R[(n,z),(n,z^{\prime })]& ={\rho }^{-n}(z){\tau }_n(z,z^{\prime }),R[(n,z),(n^{\prime },z^{\prime })]=0,n\ne n^{\prime },z\ne z^{\prime }.\end{array}$$

$$v(z)\sum _{z^{\prime }\in D}{\tau }_n(z,z^{\prime })=\sum _{z^{\prime }\in D}v(z^{\prime }){\tau }_n(z^{\prime },z).$$

$$\sum _{z^{\prime }\in D}{\tau }_n(z,z^{\prime })=\sum _{z^{\prime }\in D}{\tau }_n(z^{\prime },z),$$

$${\tau }_n(z,z^{\prime })={\beta }_n\tau (z,z^{\prime })$$

$$\rho (z)<1\text{mbox}forz\in D,$$

$$\Xi :=\sum _{z\in D}\frac{v(z)}{1-\rho (z)}=\sum _{n=0}^{\infty }\sum _{z\in D}{\rho }^n(z)v(z)<\infty .$$

$$\pi (n,z):=\eta (n,z)/\Xi ,\forall (n,z)\in {\mathbb{Z}}_{+}\times D,$$

$$\eta (n,z):={\rho }^n(z)v(z),\forall (n,z)\in {\mathbb{Z}}_{+}\times D,$$

$$\Vert P_t(x,\cdot )-\pi (\cdot ){\Vert }_{\mathrm{TV}}\to 0\text{mbox}ast\to \infty ,\text{mbox}forallx\in {\mathbb{Z}}_{+}\times D.$$

$$\sum _{(n,z)}\eta (n,z)=\sum _{(n,z)}{\rho }^n(z)v(z)=\sum _z\frac{v(z)}{1-\rho (z)}<\infty ,$$

$$\begin{array}{rl}-\eta (n,z)R[(n,z),(n,z)]& =\eta (n-1,z)R[(n-1,z),(n,z)]+\eta (n+1,z)R[(n+1,z),(n,z)]\\ & +{\sum }_{z^{\prime }\ne z}\eta (n,z^{\prime })R[(n,z^{\prime }),(n,z)],\\ -\eta (0,z)R[(0,z),(0,z)]& =\eta (1,z)R[(1,z),(0,z)]+{\sum }_{z^{\prime }\ne z}\eta (0,z^{\prime })R[(0,z^{\prime }),(0,z)].\end{array}$$

$$\eta (n,z)$$

$$\displaystyle=\rho^{n}(z)v(z)\bigg{(}R[(n,z),(n+1,z)]+R[(n,z),(n-1,z)]+\sum_{z^{\prime}\neq z}R[(n,z),(n,z^{\prime})]\bigg{)}$$

$$\displaystyle=\rho^{n}(z)v(z)\bigg{(}\lambda(z)+\mu(z)+\sum_{z^{\prime}\neq z}\rho^{-n}(z)\tau_{n}(z,z^{\prime})\bigg{)}$$

$$={\rho }^n(z)v(z)(\lambda (z)+\mu (z))+v(z)\sum _{z^{\prime }\ne z}{\tau }_n(z,z^{\prime });$$

$${\rho }^{n-1}(z)$$

$$=\lambda (z)v(z){\rho }^n(z)(\lambda (z)+\mu (z))+\sum _{z^{\prime }\ne z}v(z^{\prime }){\tau }_n(z^{\prime },z).$$

$${\tau }_n(z,z^{\prime })=\frac{{\beta }_n}{m_n-1},\forall z,z^{\prime }\in D_{n}withz\ne z^{\prime }.$$

$${\tau }_n(z(i),z(j))=\frac{1}{2}\delta (j,i\pm 1),\forall i,j\in \{1,\dots ,m_n\}.$$

$${\tau }_n(z_0,z_j)={\tau }_n(z_i,z_{i+j})={\beta }_n\delta (j,n),\forall i,j\in \mathbb{Z}.$$

$$\mathrm{d}Z(t)=g(Z(t))\mathrm{d}t+{\Sigma }^{1/2}(Z(t))\mathrm{d}B(t)+\mathcal{N}(t)+\sum _{i=1}^m{r}_i(Z(t))\mathrm{d}{\ell }_i(t),t\ge 0.$$

$$\mathcal{A}f(z)=g(z)\cdot \nabla f(z)+\frac{1}{2}\mathrm{tr}(\Sigma (z){\nabla }^{2}f(z))+{\int }_D(f(z^{\prime })-f(z))\varpi (z,\mathrm{d}{z}^{\prime }),$$

$$r_i(z)\cdot \nabla f(z)=0,z\in F_i,i=1,\dots ,m.$$

$$\rho (z):=\frac{\lambda (z)}{\mu (z)}\le 1.$$

$${\mathcal{M}}_{z}f(n)=\lambda (z)(f(n+1)-f(n))+1_{\{n\ne 0\}}\mu (z)(f(n-1)-f(n)).$$

$$\mathcal{L}f(n,z)={\mathcal{M}}_{z}f(n,z)+{\beta }_n{\rho }^{-n}(z)\mathcal{A}f(n,z),f\in \mathcal{D}.$$

$$\mathcal{D}$$

$${\mathcal{D}}_D$$

$$\Vert f{\Vert }_{D,2}:=\underset{z\in D}{\sup }\left(|f(z)|+\Vert \nabla f(z)\Vert +\Vert {\nabla }^{2}f(z)\Vert \right).$$

Full text

1. Introduction

In this paper we propose a tractable modeling approach to studying queues in an interactive random environment, where the arrival and/or service rates are modulated by a Markov process and the dynamics of the environment also depends on the state of the queue. Such models may be used in the following setting: In a service system (for example, on-demand service platforms), the demand may be affected by the service quality as indicated by dynamic ‘ratings’ which may be modeled by a Markov chain, while the ratings dynamics may depend on the congestion level in the system.

For an $M/M/1$ queue in an interactive random environment, let $N(t)$ be the queue-length process (the number of customers in the system) and $Z(t)$ be the random process for the environment. The joint process $(N(t),Z(t))$ can be modeled as a continuous-time Markov process on ${\mathbb{N}}\times\mathcal{Z}$ ($\mathcal{Z}$ representing the range of $Z(t)$), with a generator

[TABLE]

where $\mathcal{M}_{z}$ describes the queueing dynamics depending on the environment state $z$ and $\mathcal{A}_{n}$ describes the environment dynamics depending on the queueing state $n$. Specifically, given the arrival and service rates $\lambda(z)$ and $\mu(z)$, we can write

[TABLE]

On the other hand, the generator $\mathcal{A}_{n}$ can be for a general Markov process, depending on the queue length $n$. For example, for a given $n$, $\mathcal{A}_{n}$ may represent the generator of a diffusion process,

[TABLE]

or a continuous-time jump Markov chain with a transition rate matrix depending on the queueing state $n$. In the utmost generality, one can impose mild conditions on the generators $\mathcal{M}_{z}$ and $\mathcal{A}_{n}$ to guarantee the existence of an invariant measure for the joint process $(N(t),Z(t))$. However, it seems difficult to go beyond that without any structural assumptions on the joint generator, especially $\mathcal{A}_{n}$. In many applications, it is convenient to have an explicit invariant measure to work with. In general, it is hard to find an explicit form for stationary distributions of multidimensional Markov processes. (For example, in [42] it is shown that an obliquely reflected Brownian motion (RBM) in a polyhedral domain in ${\mathbb{R}}^{d}$ has a product-of-exponentials stationary distribution under the skew symmetry condition, the only case with an explicit stationary measure.)

Therefore, in order to provide an explicit expression for the invariant measure of the joint process, we study a particular multiplicative (scaled) form in the generator component $\mathcal{A}_{n}$, that is,

[TABLE]

Here $\beta_{n}$ is a positive constant, $\rho(z)=\lambda(z)/\mu(z)$ is the traffic intensity in the queue, and $\mathcal{A}f(n,z)$ is a generator corresponding to a Markov process whose transition dynamics does not depend on $n$. (In the case of reflected processes, the boundary conditions should be treated carefully; see Section 3 for details.) The scaling factors not only depend on the queue length $n$, but also include the traffic intensity $\rho(z)$. For an environment state $z$, $\rho^{-n}(z)>1$ for all queue state $n$, but the factor $\beta_{n}$ gives more flexibility (slowing down or speeding up) to the scaling of the generator $\mathcal{A}$. Our approach is motivated by applications where the environment dynamics may be sped up or slowed down by the congestion. For example, in on-demand service systems, the transitions among the different service quality ‘ratings’ may simultaneously change faster when many customers experience more congestion due to higher response rates.

We discuss two types of random environment: a pure jump Markov chain taking values in a discrete state space $D$ (finite or countable), and a reflected (jump) diffusion in a piecewise smooth domain, also denoted by $D$. Each type of environment is of its own interest. Under certain assumptions, we prove the existence of the joint invariant measure, derive its explicit expression and establish the exponential rate of convergence to the steady state (in the total variation norm). The explicit expression of the invariant measure can be regarded as a weighted geometric form (or some “product form”, although not exactly in the same sense as in the literature on stochastic networks [22, 10, 23]). Specifically, we have the joint invariant measure for $(N(t),Z(t))$ of the form $\pi(\{n\},\mathrm{d}z)=\Xi^{-1}\rho^{n}(z)\nu(\mathrm{d}z),$ where $\Xi$ is some normalization constant, and $\nu(\cdot)$ is the invariant measure associated with the generator $\mathcal{A}$. Recall that the steady-state distribution of the $M/M/1$ queue itself given an environment state $z$ is geometric ($P(N(\infty)=n)=(1-\rho(z))\rho^{n}(z)$). The product of the terms “$\rho^{n}(z)$” and “$\nu(\mathrm{d}z)$” mixes the invariant measures for the queue and the environment, despite $\rho(z)$ depending on $z$. Here the scaling factor $\rho^{-n}(z)$ in $\mathcal{A}_{n}$ is critical. For the two types of environment processes we are able to establish the exponential rate of convergence.

With a diffusive environment, our work introduces new stochastic models. The simple models include: (a) an $M/M/1$ queue with an interactive diffusive arrival rate: the arrival rate is a one-dimensional reflected (jump) diffusion process in $[0,1]$ under a fixed service rate $1$; (b) an $M/M/1$ queue with an interactive diffusive service rate: the service rate is a one-dimensional reflected (jump) diffusion process in $[1,\infty)$ under a fixed arrival rate $1$; (c) the arrival and service rates form a two-dimensional RBM in an open convex cone (with arrival rate strictly lower than service rate). RBMs have been extensively studied in the queueing (network) literature as scaling limits. However, RBMs as arrival and/or service rates have not been carefully studied. When there is no interactive behavior, the $M/M/1$ queue with a RBM arrival rate can be regarded as a special case of queues of the so-called doubly stochastic Poisson arrival processes with the arrival rate being an independent stochastic process (see, e.g., [3, 4, 5]). Our first model extends such existing interesting studies to include feedback loop from queue to environment. The second and third models with RBM being the service rate or the RBM in the wedge for both arrival and service rates are new, even in the setting of no interactive behavior. Such models are worth further careful investigation. Of course, our models go beyond RBMs, to general reflected (jump) diffusion models.

We aim to find the explicit rate of convergence to the stationary distribution in these models. For standard $M/M/1$ queues, it is well known that the rate of convergence is exponential, see, e.g., [36, Proposition 5.8]. However, for diffusion processes (solutions of SDEs), reflected diffusions, and their versions with jumps, the characterization of an explicit rate of convergence to steady state (as opposed to simply proving that there exists an exponential rate of convergence) is quite a challenging problem. See, e.g., [11, 38, 37, 20]. Thus, it is a much more difficult problem to study the rate of convergence for the joint Markov process with a generator in the general form in (1.1) due to the complicated interactive behavior of the two processes (one being discrete and the other being continuous). We attempt to solve this problem via a coupling technique for the joint process $(N,Z)$. We provide a novel way to construct the coupling time for the joint process in order to prove that the convergence rate is exponential, and more importantly, provide good estimates of the rate of convergence via careful studies of the exponential bounds for the coupling time. This appears to be the first work in the literature to carefully find the estimates of the coupling times of joint processes for queueing processes in random environments.

Although our main focus is on the multiplicative (scaling) form in the generator $\mathcal{A}_{n}$, we have also considered a setup where the the environment jump diffusion described above depend on the queue length $n$ via its domain $D_{n}\subseteq D$. In particular, the drift vector field, covariance matrix field, and the jump measure remain the same for all $n$, but reflection vector fields may depend on the queueing state $n$. The entire domain $D$ is the union of these $D_{n},\,n=0,1,2,\ldots$. We assume that this reflected jump-diffusion in $D_{n}$ has a unique invariant probability measure $\nu^{(n)}_{D_{n}}$ inside the domain $D_{n}$, which is the projection of a certain finite measure on $D$ to $D_{n}$. (The corresponding boundary measures may depend on $n$.) See Assumptions 3.3–3.6. We prove similar results as above in this setting. We construct two special examples: an $M/M/1$ queue with a fixed service rate and a reflected diffusion arrival rate, controlled based on a threshold of queue length (Example 3.1) and an $M/M/1$ queue with a fixed arrival rate and a diffusive service rate, controlled similarly (Example 3.2).

When the random environment is a Markov chain taking discrete values, our results also extend to the generator $\mathcal{A}_{n}$ of the form $\rho^{-n}(z)\tau_{n}(z,z^{\prime})$, where the generator rate $\tau_{n}$ may depend on the queueing state $n$ unlike the multiplicative case. However, it is assumed that an invariant measure associated with the transition rate $\tau_{n}(z,z^{\prime})$ exists such that it is independent of the queueing state $n$ (Assumption 2.1). This is slightly more general than the multiplicative case, so we state the model and results in section 2 in this setup. We also give an example to illustrate how this slightly more general setup is used (see Examples 2.1 and 2.2).

1.1. Literature review on queues in interactive environments

Queues in random environments (e.g., Markov modulated models) have been extensively studied in the literature. Most of the literature assumes that the queueing dynamics is affected by the environment, but not interactive. For example, the paper [35] studies Markov-modulated arrival and service rates with finite environment space, and finds expressions of waiting times. The paper [41] deals with similar questions by comparing this queue with an appropriate $M/M/1$ queue. Optimization of service rate for the case when arrival rate is a Markov process is studied in [26]. See also, a birth–death process in random environment [13] and a Markov chain in Markov environment, studied in [12, 16, 33]. A particular case of a Markov–modulated setting is when the service dynamics is subject to interruptions. In this case, the random environment only affects service rate $\mu$. The survey [25] summarizes the existing literature on this topic.

In the Markov-modulated queueing literature, the arrival or service rates under modulation take finite or countable number of values. However, in practice, the rates under modulation can possibly take continuous values. Our work thus goes beyond the existing frameworks and develops new queueing models.

In [19], the authors study a random particle (a distinguished customer) walking randomly over the sites of a symmetric Jackson network (open or closed), where the arrival rate of a station/node or the transition of customers from it to other stations/nodes is affected if the particle occupies it, while the jump rate of the particle depends on the state of the station/node it currently occupies. An explicit steady state distribution for the joint process is derived. In [24], Jackson networks in interactive random environments are studied, where the service capacities are affected by the environment, while customer departing may enforce the environment to jump immediately. An explicit expression of the product form is derived for the joint queueing and environment processes. Inspired by [19], a different construction of Markov processes in random environments resulting in product-form invariant measure is provided. In [6], various Markov processes with interactive random environment are constructed. This paper is in the same flavor as that in [6]. The paper [17] deals with feedback loop created by blocking some channels in a multi-server queue, and finds a product-form stationary distribution for the joint process. None of these papers investigate the rate of convergence to stationarity. Our model of a single-server queue is also constructed in a more general manner.

The papers [14, 44] study birth-death processes in random environment with feedback. This is a more general setup than in our paper, because an $M/M/1$ queue is a particular case of a birth-death process. However, [14] is concerned with explosion questions, rather than stationary distributions and convergence rates, and [44] focuses on generating function approach and achieves only partial results for the steady-state distribution.

1.2. Notation

The integral with respect to the measure $\nu$ applied to the function $f$ is written as $\langle\nu,f\rangle$. Exponential distribution with rate $\alpha$ is denoted by $\operatorname{Exp}(\alpha)$. The arrow $\Rightarrow$ indicates weak convergence. The dot product of two vectors $a$ and $b$ is denoted by $a\cdot b$. We say two finite measures $\mu,\nu$ on $\mathbb{R}$ satisfy $\mu\preceq\nu$ if for all $u\in\mathbb{R}$ we have $\mu(-\infty,u]\leq\nu(-\infty,u]$, but $\mu(\mathbb{R})=\nu(\mathbb{R})$. We say that $\mu$ is stochastically dominated by $\nu$. We transfer this concept to random variables: $X$ is stochastically dominated by $Y$ if the distribution of $X$ is stochastically dominated by the distribution of $Y$. Let ${\mathbb{Z}}_{+}=\{0,1,2,\ldots\}$ and ${\mathbb{R}}_{+}:=[0,\infty)$. Define the total variation norm: For a signed measure $\mathbf{\nu}$, let $\lVert\mathbf{\nu}\rVert_{\operatorname{TV}}:=\sup_{A}|\mathbf{\nu}(A)|$. Throughout this article, we consider continuous-time random processes (unless otherwise noted) on a filtered probability space $(\Omega,\mathcal{F},(\mathcal{F}_{t})_{t\geq 0},\mathbb{P})$ with the filtration satisfying the usual conditions.

1.3. Organization of the paper

In Section 2, we study the model in an interactive jump environment. In Section 3, we study the single-server queue with a reflected jump diffusion environment. In Section 4, we estimate the explicit rate of exponential convergence for the case of compact environment state space, for both models in Sections 2 and 3. In Section 5 we state and prove some auxiliary lemmata. We make some concluding remarks in Section 6.

2. $M/M/1$ queue in an interactive jump environment

Consider an $M/M/1$ queue with an infinite waiting space operating in an interactive jump environment described as follows. Let $D$ be a finite or countable state space. For every $n\in\mathbb{Z}_{+}$, let $\mathbf{T}_{n}=(\tau_{n}(z,z^{\prime}))_{z,z\in D}$ be the generator of an irreducible continuous-time Markov chain on $D$; this (finite or countable-sized) matrix is called nominal jump intensity matrix for the jump process $Z$ in the queueing state $n$. We define a two-component Markov process $(N,Z)$ taking values in the countable state space ${\mathbb{Z}}_{+}\times D$ with the following generator matrix \mathbf{R}=\big{(}R[(n,z),(n^{\prime},z^{\prime})]\big{)}:

[TABLE]

where $\rho(z):=\lambda(z)/\mu(z)$ for each $z\in D$. Here $N=\{N(t):t\geq 0\}$ represents the number of jobs in the system (including those in queue and in service), taking values in ${\mathbb{Z}}_{+}$, and $Z=\{Z(t):t\geq 0\}$ represents a jump process taking values in $D$. When the environment is in state $z$, the arrival and service rates for the queueing process are $\lambda(z)$ and $\mu(z)$, respectively, both depending on state $z$.

When the queue size is in state $n$, the transition of the environment $Z$ from state $z$ to state $z^{\prime}$ occurs at the rate $\rho^{-n}(z)\tau_{n}(z,z^{\prime})$. Note that the fourth equation in (2.1) does not allow simultaneous jumps for $N$ and $Z$. It is evident that the pair $(N,Z)$ is a well-defined Markov process in ${\mathbb{Z}}_{+}\times D$ with the generator $\mathbf{R}$.

Remark 2.1**.**

We do not multiply this transition rate $\tau_{n}$ by a factor $\beta_{n}$: Dependence on $n$ is already enshrined in the rate $\tau_{n}$. We impose a condition (2.2) to guarantee the product form of the steady state.

We first make the following assumption on the nominal jump intensity matrix $\mathbf{T}_{n}$.

Assumption 2.1**.**

For each $n\in{\mathbb{Z}}_{+},\,z\in D$, and for some function $v:D\rightarrow{\mathbb{R}}_{+}$,

[TABLE]

For fixed $n\in{\mathbb{Z}}_{+}$, if we define a Markov process $\tilde{Z}_{n}:=\{\tilde{Z}_{n}(t):t\geq 0\}$ on $D$ with the nominal jump intensity matrix $\mathbf{T}_{n}$ as the generator, then (2.2) implies that $v(\cdot)$ defines an invariant measure for $\tilde{Z}_{n}$. If $\sum_{z\in D}v(z)<\infty$, then this measure can be normalized to a probability distribution. If

[TABLE]

then the counting measure is invariant for $\tilde{Z}_{n}$; if $D$ is a finite set, then it is normalized to a uniform distribution on $D$. It is important to note that the invariant measure $v(\cdot)$ does not depend on $n$, although the jump intensity matrix $\mathbf{T}_{n}$ depends on $n$.

Remark 2.2**.**

A simplest example is when $\tau_{n}(z,z^{\prime})$ has a multiplicative form:

[TABLE]

for some transition rate matrix $\tau(z,z^{\prime})$ satisfying $v(z)\sum_{z^{\prime}\in D}\tau(z,z^{\prime})=\sum_{z^{\prime}\in D}v(z^{\prime})\tau(z^{\prime},z)$. However, we provide examples below in which $\tau_{n}(z,z^{\prime})$ depends on $n$ in a nontrivial manner while the existence of $v$ independent of $n$ is guaranteed. See Examples 2.1 and 2.2.

Assumption 2.2**.**

The functions $\rho,v$ satisfy

[TABLE]

[TABLE]

Note that the constant $\Xi$ is the normalization constant in the joint invariant measure $\pi$ in (2.5).

Theorem 2.1**.**

Under Assumptions 2.1 and 2.2, the Markov process $(N,Z)$ is irreducible, aperiodic, and positive recurrent. It has an invariant probability measure

[TABLE]

where $\Xi$ is given in (2.4), and

[TABLE]

This process has transition kernel $P_{t}(x,\cdot)$ which converges to this invariant measure:

[TABLE]

Proof. We first show that the process $(N,Z)$ is irreducible and aperiodic. It follows from the observation that for every $t>0$, $(n,z),(n^{\prime},z^{\prime})\in\mathbb{Z}_{+}\times D$, one can with positive probability get from $(n,z)$ to $(n^{\prime},z^{\prime})$ in time $t$. For the measure $\eta$ from (2.6) to be finite, we need

[TABLE]

which is implied by (2.3)–(2.4) in Assumption 2.2. If we prove that $\eta$ from (2.6) is indeed an invariant measure, the positive recurrent property follows from [34, Theorem 3.5.3], [40, Theorem 2.7.18], and then the ergodicity, as in (2.7), follows from [32]. To verify that $\eta(n,z)$ in (2.6) is an invariant measure, we show that $\eta^{\prime}\mathbf{R}=0$. Let us show that for all $n=1,2,\ldots$ and $z\in D$,

[TABLE]

By (2.1), the left- and right-hand sides of the first equation in (2.8) are equal to, respectively,

[TABLE]

From (2.2) in Assumption 2.1, the last terms in the right-hand side of these two last equations are equal. This proves the first equation in (2.8); the second one is similar. This completes the proof.

Example 2.1**.**

($D$ as a union of finite sets)* In Examples 2.1 and 2.2, $\delta(i,j)$ stands for the Kronecker delta. Given $n\in{\mathbb{Z}}_{+}$, let $D_{n}$ be a finite set in $(0,1)$ with cardinality $m_{n}$. For definiteness, assume that $1<m_{n}<M$ where $M\in{\mathbb{Z}}_{+}$ is a fixed value. Introduce an enumeration of points in each $D_{n}$: $D_{n}=\{z(1),\ldots,z(m_{n})\}$ (say, in a increasing order) and make a convention that $z(0)=z(m_{n})$, $z(m_{n}+1)=z(1)$. Sets $D_{n}$ can have common points for different $n$ or be pair-wise disjoint. Set $D=\operatornamewithlimits{\cup}\limits_{n}D_{n}$ and $\upsilon(z)=1$ for $z\in D$. Set $D$ can be finite or countable.*

Next, take a subset ${\mathbb{L}}\subseteq{\mathbb{Z}}_{+}$ (${\mathbb{L}}$ or ${\mathbb{Z}}_{+}\setminus{\mathbb{L}}$ can be empty). For $n\in\mathbb{L}$, set

[TABLE]

For $n\in{\mathbb{Z}}_{+}\setminus{\mathbb{L}}$, set

[TABLE]

Here $\beta_{n}\in(0,\infty)$ are scaling constants depending on $n$ (which is irrelevant for the invariant measure of the process $(N,Z)$). Pictorially, $\tau_{n}$ for $n\in\mathbb{L}$ describes uniform jumps on $D_{n}$ while for $n\in{\mathbb{Z}}_{+}\setminus{\mathbb{L}}$, $\tau_{n}$ yields a ‘nearest-neighbor’ walk with cyclic (periodic) boundary condition. Either way, the counting measure $\upsilon$ is invariant; cf. Assumption 2.1. Thus, (2.2) holds true.

Then \mathbf{T}_{n}=\big{(}\tau_{n}(z,z^{\prime})\big{)} generates a Markov chain $\tilde{Z}_{n}$ with an invariant probability measure ${\mathbf{1}}_{D_{n}}(z)/m_{n}$, $z\in D$. The invariant measure $\eta$ is then given in (2.6) with $\eta(n,z)=z^{n}$.

Example 2.2**.**

($D$ as a countable set, $\tau_{n}$ as a null-recurrent jump chain.)* Assume that $D\subset(0,1)$ is countable, and can be enumerated by $i=0,\pm 1,\pm 2$, so that $\rho_{i}:=\rho_{z_{i}}$ for $i\geq 0$ satisfies $\rho_{0}<\rho_{1}<\cdots<1$ and $\lim_{i\rightarrow\infty}\rho_{i}=1$. (Enumeration with labels $i=-1,-2,\ldots$ does not matter.) Set $\upsilon(z)=1$ and*

[TABLE]

Here, as earlier, $\beta_{n}$ is a scaling constant depending on $n$ (again irrelevant for the invariant measure of the process $(N,Z)$). Then \mathbf{T}_{n}=\big{(}\tau_{n}(z,z^{\prime})\big{)} generates a null-recurrent Markov chain $\tilde{Z}_{n}$ with the invariant measure $\upsilon(z)=1$, $z\in D$. Thus, the random traffic intensity $\rho_{\tilde{Z}_{n}}$, depending on both the state of the queue and the environment, will approach the critical value $1$ infinitely often. However, under the condition (2.4) the resulting Markov process $(N,Z)$ is positive recurrent, with an invariant measure $\eta(n,z)=\rho^{n}(z)$ for $(n,z)\in{\mathbb{Z}}_{+}\times D$.

3. $M/M/1$ queue in an interactive diffusive environment

3.1. Reflected jump-diffusions

In this section we consider the queue with $\lambda$ and $\mu$ dependent on a diffusive environment process $Z(t)$. First, let us define the dynamics of this environment process as a reflected (jump) diffusion in a certain domain in $\mathbb{R}^{d}$.

It is instrumental to recapitulate some basic notion. A domain in ${\mathbb{R}}^{d}$ is the closure of an open connected subset. A domain $D$ is called smooth if its boundary $\partial D$ is a $(d-1)$-dimensional $C^{2}$ manifold. Take $m$ smooth domains $D_{1},\ldots,D_{m}$ in ${\mathbb{R}}^{d}$. Assume $D=\cap_{i=1}^{m}D_{i}$ has boundary $\partial D$ with $m$ faces: $F_{i}:=\partial D\cap\partial D_{i}$ which are $(d-1)$-dimensional manifolds with an edge, and such that all $m$ domains are essential: Removal from the intersection of any domain will change the result. Then $D$ is called a piecewise smooth domain in ${\mathbb{R}}^{d}$. Define by $\mathbf{n}_{i}(z)$ the inward unit normal vector to $\partial D_{i}$ at $z\in F_{i}$. Inward in this case is defined as pointing inside $D_{i}$, even if this is not inside $D$. An important example is a convex polyhedron with $D_{i}$ being half-spaces. Of particular interest is the positive orthant $D={\mathbb{R}}^{d}_{+}$. Of course, smooth domains also belong to this class of domains, with $m=1$.

Take continuous functions $g:D\to\mathbb{R}^{d}$ and $\Sigma:D\to\mathbb{R}^{d\times d}$ such that the matrix $\Sigma(z)=(a_{ij}(z))$ is symmetric and positive definite for all $z\in D$, and there exists a $\delta>0$ such that $\Sigma(z)v\cdot v\geq\delta\lVert v\rVert^{2}$ for all $v\in\mathbb{R}^{d}$ and $z\in D$. For every $z\in D$, define a finite measure $\varpi(z,\cdot)$ on $D$ such that $\varpi(z,\cdot)\Rightarrow\varpi(z^{0},\cdot)$ as $z\to z^{0}$ in $D$. Recall that $\Rightarrow$ denotes weak convergence. Take $r_{i}:F_{i}\to\mathbb{R}^{d}$: continuous functions, pointing inside $D$; that is, $r_{i}(z)\cdot\mathbf{n}_{i}(z)>0$ for $i=1,\ldots,m,\,z\in F_{i}$. Let us define a reflected jump-diffusion: a process $Z=\{Z(t):\,t\geq 0\}$ in $D$ with drift vector field $g$, diffusion matrix field $\Sigma$, jump measures $\varpi(z,\cdot)$ and reflection vector fields $r_{1},\ldots,r_{m}$.

This process will be adapted and right-continuous with left limits. Take a $d$-dimensional Brownian motion $B=\{B(t):\,t\geq 0\}$, adapted to the filtration. Take continuous nondecreasing processes $\ell_{i}=\{\ell_{i}(t):\,t\geq 0\}$ for $i=1,\ldots,m$ such that $\ell_{i}$ can grow only when $Z(t)\in F_{i}$, another right-continuous process with left limits $Z=(Z(t),\,t\geq 0)$ with values in $D$, and yet another process $\mathcal{N}=\{\mathcal{N}(t):t\geq 0\}$ which is right-continuous piecewise constant, with jump measure $\varpi(Z(t-),\cdot)$, and such that

[TABLE]

We assume the equation (3.1) has a well-defined unique weak solution, and forms a Feller continuous strong Markov semi-group, with generator

[TABLE]

which consists of a nondegenerate uniformly elliptic diffusion and a state-dependent finite jump measure. This existence and uniqueness were proved under Lpischitz conditions on vector field $g(\cdot)$ and the matrix $(a_{ij}(\cdot))$, as well as continuity of $r_{i}(\cdot)$ for each $i=1,\ldots,m$, and some additional technical conditions. The case without jumps was proved in [28]; the general case follows from the standard construction by piecing out, [39]. The reflection at the boundary translates into boundary conditions for (3.2):

[TABLE]

The dynamics of this process can be described as follows:

• •

As long as it is strictly inside $D$, this process behaves as a jump-diffusion in $d$ dimensions with drift vector field $g$ diffusion matrix field $\Sigma$, and family $\varpi$ of jump measures. These jump measures are such that the process does not jump out of $D$.

• •

At a point $z\in F_{i}$, $i=1,\ldots,m$, it is reflected back inside the domain $D$, according to the vector $r_{i}(z)$.

• •

If it hits the lower-dimensional edges: intersections of two or more faces $F_{1},\ldots,F_{m}$, it is reflected back inside $D$ according to a positive linear combination of reflection vectors corresponding to these intersecting faces.

Normal reflection corresponds to the case when $r_{i}(z)=\mathbf{n}_{i}(z)$, where $z\in F_{i}$ and $i=1,\ldots,m$.

Remark 3.1**.**

In the case of a diffusion without reflection, the state space may be $\mathbb{R}^{d}$, or still some subset $D$. The latter happens if the drift coefficient is sufficiently large to compel the process to stay in a certain domain. An example of this is the drift for a Bessel process on the half-line, see [21, Chapter 3, Problem 3.23].

In the case $d=1$, for a reflection on $[a,b]$, we have a normal reflection, and the boundaries consisting of two pieces $\{a\}$ and $\{b\}$. For a reflection on $[a,\infty)$, we have a normal reflection again, with the boundary $\{a\}$.

3.2. Construction of the joint Markov process

Let us now use symbol $z$ for a point in $D$ (instead of $x$). Take continuous functions $\lambda,\mu:D\to(0,\infty)$ with $\lambda(z)\leq\mu(z)$ for $z\in D$. Define the traffic intensity:

[TABLE]

For every $z\in D$, consider an $M/M/1$ queue with arrival intensity $\lambda(z)$ and service intensity $\mu(z)$, where $n$ is the state of this queue. The process $\tilde{N}$ counting the number of jobs in the system, called the queueing process in the sequel, is a continuous-time Markov process on $\mathbb{Z}_{+}$ with generator

[TABLE]

We now consider a $(1+d)$-dimensional Markov process $(N,Z)=\{(N(t),Z(t)):t\geq 0\}$ with values in $\mathbb{Z}_{+}\times D$ which evolves as follows:

(a) If $N(t)=n\in\mathbb{Z}_{+}$, then $Z$ behaves as a reflected jump-diffusion in $D$ with generator $\rho^{-n}(z)\beta_{n}\mathcal{A}$ and reflection fields $r_{1},\ldots,r_{m}$.

(b) if $Z(t)=z$, then $N(t)$ jumps from $n$ to $n+1$ with intensity $\lambda(z)$, and (if $n\neq 0$) to $n-1$ with intensity $\mu(z)$.

Here $\beta_{n}$ is the variability coefficient for the diffusive environment, depending on the queueing state $n$, and $\rho^{-n}(z)$ is the queueing impact factor, capturing the impact from the traffic intensity (congestion) from the queueing process. The component $N$ can be informally described as the queueing process of an $M/M/1$ queue with arrival and service rates, $\lambda(z)$ and $\mu(z)$, respectively. These rates depend on an auxiliary process $Z$. The dynamics of $Z$, however, depends on the current position of this queueing process. Therefore, we call such a system as an $M/M/1$ queue in an interactive diffusive environment.

The joint dynamics is described via a combined Markov process $(N,Z)$ with the following generator:

[TABLE]

Here, $\mathcal{D}$ stands for the following subspace of the domain of $\mathcal{L}$:

[TABLE]

(Note that we were intentionally loose on the domains of $f$ in (3.2) and (3.5), but they are clear from this definition.) From the general theory of piecing out it follows that this is a Feller process, see [39]. We denote by $C^{2}_{b}(D)$ the set of twice continuously differentiable functions $D\to\mathbb{R}$ which are bounded with their first and second derivatives (the last condition, automatically fulfilled for bounded $D$). This is a separable Banach space with the norm

[TABLE]

Denote by $P^{t}(y,\cdot)$ the transition kernel of $(N,Z)$ where $y=(n,z)\in\mathbb{Z}_{+}\times D$.

We give three special cases to illustrate the construction above.

• (a)

$M/M/1$ queue with an interactive diffusive arrival rate.

Assume that $\lambda(z)=z$ and $\mu\equiv 1$. Let $D=[0,1]$, and the generator $\mathcal{A}$ in (3.2) be that of a reflected diffusion in $(0,1)$ without jumps. The reflections at [math] and $1$ correspond to the Neumann boundary conditions:

[TABLE]

• (b)

$M/M/1$ queue with an interactive diffusive service rate.

Assume that $\lambda\equiv 1$ and $\mu(z)=z$. Let $D=[\mu_{0},\infty)$ for some $\mu_{0}\geq 1$, and the generator $\mathcal{A}$ in (3.2) be that of a simple RBM on $[\mu_{0},\infty)$ without jumps. The reflection at $\mu_{0}$ satisfies the Neumann boundary condition.

• (c)

$M/M/1$ queue with both diffusive arrival and service rates.

Take $D=\{(z_{1},z_{2})\in{\mathbb{R}}^{2}_{+}:z_{2}\leq z_{1}\}$ be a cone in the positive orthant, and the generator $\mathcal{A}$ in (3.2) be that of a two-dimensional Brownian motion in $D$ with normal reflections at the boundary. Let $(\lambda(z),\mu(z))=z$. Then the arrival and service rates of the $M/M/1$ queue follow the dynamics of a reflected RBM in $D$ in the interactive manner described above.

3.3. Invariant measures

We need the following assumptions on some properties of the reflected jump-diffusion process. The first assumption states that for each level, there exists a steady-state distribution. The second assumption ensures that the whole process has a steady-state distribution.

Assumption 3.1**.**

Assume the (reflected) jump-diffusion with generator $\mathcal{A}$ is positive recurrent, and has a unique stationary/invariant measure $\nu_{D}$, together with boundary measures $\nu_{F_{i}},\,i=1,\ldots,m$. This means that the stationary copy of this process $\tilde{Z}^{*}=\{\tilde{Z}^{*}(t):\,t\geq 0$}* with $\tilde{Z}^{*}(t)\sim\nu_{D}$ for $t\geq 0$, satisfies the following condition: For every $t\geq 0$, each $i=1,\ldots,m$, and every bounded function $f:F_{i}\to\mathbb{R}$,*

[TABLE]

where $\ell_{i}(s)$ is the nondecreasing process in (3.1).

Assumption 3.2**.**

The measure $\nu_{D}(\cdot)$ satisfies

[TABLE]

Note that $\Xi$ is the normalization constant in the joint invariant measure of $(N,Z)$ in (3.10). This invariant measure on each boundary $F_{i}$ has value zero.

Remark 3.2**.**

Similarly to (3.8), we can define the concept of boundary measures for the joint process $(N,Z)$. First, construct the boundary process $\ell_{i}=(\ell_{i}(t),\,t\geq 0)$ for the component $Z$ and face $F_{i}$ of the boundary $\partial D$. Assume $0=\rho_{0}<\rho_{1}<\ldots$ are jump times for $N$. Then $Z(\rho_{k}+t)$ for $t\in[0,\rho_{k+1}-\rho_{k}]$ behaves as a reflected jump-diffusion on $D$ with generator $\rho^{-n_{k}}(z)\beta_{n_{k}}\mathcal{A}$ and reflection fields $r_{1},\ldots,r_{m}$, with $N(t)=n_{k}$ for $t\in[\rho_{k},\rho_{k+1})$. Thus there exist a continuous nondecreasing process $\ell_{i}^{(k)}(t),\,t\in[0,\rho_{k+1}-\rho_{k}]$ such that (3.1) holds with adjusted drift vector field, diffusion matrix field, and jump measures family. Define

[TABLE]

using induction over $k$. This defines $\ell_{i}=(\ell_{i}(t),\,t\geq 0)$ for $i=1,\ldots,m$. Next, define a boundary measure $\nu_{F_{i}}$ on the face $F_{i}$ corresponding to a stationary distribution $\pi$ for this joint process $(N,Z)$: Take the corresponding stationary copy $(N^{*},Z^{*})$ with $(N^{*}(t),Z^{*}(t))\sim\pi$ for $t\geq 0$. For a bounded function $f:\mathbb{Z}_{+}\times F_{i}\to\mathbb{R}$ and a $t\geq 0$,

[TABLE]

Now we are ready to state and prove the main result of this section.

Theorem 3.1**.**

Under Assumptions 3.1 and 3.2, there is a unique invariant measure for $(N,Z)$:

[TABLE]

The corresponding boundary measures $\pi_{i}$ for $F_{i}$ (if there is reflection) are given by

[TABLE]

Finally, this Markov process is ergodic: for every $y\in\mathbb{Z}_{+}\times D$,

[TABLE]

Proof.

From stationarity we immediately get: for all $f\in C^{2}_{b}(D)$,

[TABLE]

This is called the basic adjoint relationship in the literature. We refer to [43] for its deduction in the case of a convex polyhedron; the same is true for a general piecewise smooth domain $D$, as in our case. Apply [27, Theorem 1.7, Theorem 2.2, Lemma 2.4, Remark 2.5] using their notation, with the state space $E=\mathbb{Z}_{+}\times D$; $U=\{0,1,\ldots,m\}$, where the point [math] corresponds to the domain $D$ itself, and $i=1,\ldots,m$, correspond to faces $F_{1},\ldots,F_{m}$ of the boundary; for all $z\in D$, $n\in\mathbb{Z}_{+}$, and $u\in U$,

[TABLE]

We need to check [27, Condition 1.2] on the absolutely continuous generator $A$ and the singular generator $B$. Let

[TABLE]

Part (i) requires that $A,B:\mathcal{D}\subset C_{b}(E)\to C(E\times U)$, and the unity function $\mathbf{1}(n,z)=1$ for $(n,z)\in E$ satisfies $\mathbf{1}\in\mathcal{D}$, $A\mathbf{1}=0$, and $B\mathbf{1}=0$. This is trivially satisfied.

Part (ii) requires that there exist $\psi_{A}(n,z)$ and $\psi_{B}(n,z)$ in $C(E\times U)$, $\psi_{A},\psi_{B}\geq 1$ and constants $a_{f},b_{f}$, $f\in\mathcal{D}$ such that

[TABLE]

where $\mathcal{U}$ is any closed set of $E\times U$. We can take $a_{f}=b_{f}:=\lVert f\rVert_{D,2}$, and

[TABLE]

Part (iii) requires the following: Defining $(A_{0},B_{0})=\{(f,\psi_{A}^{-1}Af,\psi_{B}^{-1}Bf):f\in\mathcal{D}\}$, $(A_{0},B_{0})$ is separable in the sense that there exists a countable collection $\{g_{k}\}\subset\mathcal{D}$ such that $(A_{0},B_{0})$ is contained in the bounded, pointwise closure of the linear span of $\{(g_{k},A_{0}g_{k},B_{0}g_{k})=(g_{k},\psi_{A}^{-1}Ag_{k},\psi_{B}^{-1}g_{k})\}$. This is proved by taking a dense countable subset $\Upsilon$ of $C^{2}_{b}(D)$ in the norm $\lVert\cdot\rVert_{2}$, and then taking a countable subset

[TABLE]

This subset is dense in the sense of pointwise convergence.

Part (iv) requires that for each $u\in U$, the operators $A_{u}$ and $B_{u}$ defined by $A_{u}f(x)=Af(x,u)$ and $B_{u}f(x)=Bf(x,u)$ are pre-generators. This follows from [27, Remark 1.1], because all these operators satisfy the positive maximum principle.

Part (v) requires that $\mathcal{D}$ is closed under multiplication and separates points. This follows directly from the definition.

Finally, we need to prove the main condition as in [27, Theorem 1.7, (1.17)]:

[TABLE]

From (3.14) and (3.6), canceling $\beta_{n}$ and $\rho^{n}(z)$ when appropriate, we rewrite the left-hand side of (3.15) as follows:

[TABLE]

The first line in (3.16) is equal to zero; this follows from (3.13). Let us show that the second line in (3.16) is equal to zero, too. For every $z\in D$, $\mathcal{M}_{z}$ is the generator of the $M/M/1$ queue with arrival and service rates $\lambda(z)$ and $\mu(z)$. This queue has geometric stationary distribution $(1-\rho(z))\rho^{n}(z),\,n\in\mathbb{Z}_{+}$. Thus

[TABLE]

Integrating (3.17) with respect to $\mu_{D}(\mathrm{d}z)$, we get: The second line in (3.16) is equal to zero. We interchanged integration and series, which we can do by uniform boundedness of $f$ combined with Assumption 3.2. This completes the proof of (3.15), and with it [27, (1.17)]. Next, $K_{1}:=\partial D$ is the closed support for $\mu_{1}^{E}$. By [27, Remark 2.5], the results of [27, Lemma 2.4] hold, and we can apply [27, Theorem 2.2 (f)], and obtain the stationary copy of our process $(N,Z)$.

We have written the proof for reflected diffusions. For non-reflected ones, it is simpler: we can simply verify (3.13), which in our case then becomes

[TABLE]

This is done similarly to the computation above, but without all boundary terms. The lack of reflection obviates the need to apply results cited above from [27].

Finally, ergodicity follows from [32, Theorem 6.1] in the following way (for terminology, we refer the reader to this cited article [32]). Our process is positive Harris recurrent, since the invariant measure is finite. Meanwhile, every skeleton chain is irreducible, because of the following irreducibility property. Define a Lebesgue measure on $\mathbb{Z}_{+}\times D$ as a sum of Lebesgue measures on each layer of this set.

Lemma 3.2**.**

For every $n\in\mathbb{Z}_{+},\,z\in D$, and a subset $G\subseteq\mathbb{Z}_{+}\times D$ of positive Lebesgue measure,

[TABLE]

Proof.

Without loss of generality, assume $G=\{m\}\times E$ for a subset $E\subseteq D$ of positive Lebesgue measure, and $m\geq n$. We prove the statement (3.19) by induction over $m$.

Induction Base: $m=n$. Consider the probability

[TABLE]

that, starting from $y=(n,z)$, the joint process $(N,Z)$ at time $t$ will be in $\{n\}\times E$. This probability is bounded from below by

[TABLE]

This probability $Q^{t}_{n}(z,E)$, in turn, is estimated from below by (with $z_{*}>0$ fixed later):

[TABLE]

Here, $q_{*}$ is the probability that, starting from $Z_{n}(0)=z$, the reflected jump-diffusion $Z_{n}$ in $D$ with generator $\rho^{-n}(z)\mathcal{L}$ and reflection vector fields $r_{1},\ldots,r_{m}$ ends at $Z_{n}(t)\in E$ and $\lVert Z_{n}(s)\rVert\leq z_{*}$ for $s\in[0,t]$. It follows from known properties of reflected jump-diffusions with nonsingular covariance matrix $\Sigma(\cdot)$ that $q_{*}>0$ for large enough $z_{*}>0$. This, together with (3.20) and (3.21), proves that

[TABLE]

Thus we have proved the statement (3.19) for $m=n$.

Induction Step: First, consider the case $m=n+1$. This probability $P^{t}(y,G)$ is estimated from below by the probability that for some time $\tau\in[0,t]$, the process $N$ will stay at level $n$, then jump at time $\tau$ at level $n+1$ and stay there until time $t$, and $Z(t)\in E$. If $\hat{\mu}$ is the distribution of $\tau$ (which is a positive measure on $[0,t]$), then

[TABLE]

It suffices to show that the double integral in the right-hand side of (3.23) is positive. Indeed, from (3.22) we get: $Q^{s}_{n}(y,E^{\prime})>0$ for $E^{\prime}\subseteq D$ of positive Lebesgue measure, and $\quad Q^{t-s}_{n+1}(w,E)>0$. In addition, $\hat{\mu}$ is a positive measure on $[0,t]$. Use twice the observation that the integral of a positive function over a positive measure is positive, and complete the proof that the right-hand side (and therefore the left-hand side) in (3.23) is positive.

Assuming we proved (3.19) for $m=n+k$, $k\geq 0$, let us prove this for $m=n+k+1$:

[TABLE]

This follows from the same logic: The function $P^{t/2}((n+k,w),\{n+k+1\}\times E)$ is positive by the previous part of the induction step, applied to $n+k$ instead of $n$, and to $n+k+1$ instead of $m=n+1$. The measure $P^{t/2}(y,(n+k,\mathrm{d}w))$ is positive by the induction hypothesis. This completes the proof of this lemma. ∎

Using Lemma 3.19, we have shown ergodicity as in (3.12). Earlier, we have proved (3.10) and (3.11). Thus we have completed the proof of Theorem 3.12. ∎

Remark 2.3. The crucial property is that for each $n\in{\mathbb{Z}}_{+},\,z\in D,\,t>0,V,V^{\prime}\subseteq D$,

[TABLE]

In fact, further generalizations depend on whether an analog of this equality can be established. Here ${\tt p}^{t}$ stands for the transition density for the diffusion with generator $\mathcal{A}$ in (3.2).

3.4. A more general setup

We offer a similar result under a more general feedback scheme. For $n=1,2,\ldots$, fix a piecewsie smooth domain $D_{n}\subseteq D$ with $m_{n}$ faces of the boundary $\partial D_{n}$:

[TABLE]

and corresponding reflection vector fields

[TABLE]

For each $n\in\mathbb{Z}_{+}$, this domain $D_{n}$, its boundary $\partial D_{n}$ with faces (3.26), and reflection vector fields (3.27) satisfy the same assumptions enunciated at the very beginning of Section 3, as the original domain $D$ and reflection vector fields $r_{1},\ldots,r_{m}$. In addition, we impose the following assumptions on domains $D$ and $D_{n}$.

Assumption 3.3**.**

For all $n\in\mathbb{Z}_{+}$, $D_{n}\cap D_{n+1}$ contains an open subset of $D$; and $D=\cup_{n\in\mathbb{Z}_{+}}D_{n}$.

For every level $N(t)=n$ of the queue-size component, the environment variable $z\in D$ is kept fixed when $z\in D\setminus D_{n}$ and follows a reflected jump-diffusion process in $D_{n}$ as in (3.1) where parameters vary with $n$. In other words, the process ${\widetilde{Z}}_{n}$ lives in $D$ but its mechanism depends on $n$. The generator ${\mathcal{A}}_{n}$ of ${\widetilde{Z}}_{n}$ has the form

[TABLE]

for $z\in D_{n}$, and $\mathcal{A}_{n}f(z)=0$ for other $z$. The generator $\mathcal{L}$ of the joint process, instead of (3.6), has the following form:

[TABLE]

Here $\widetilde{\mathcal{D}}$ is the following domain, defined similarly to (3.7):

[TABLE]

Note that in this setup, the dependence of the generator $\mathcal{A}_{n}$ on the queueing state $n$ is only through the domain $D_{n}$ while the drift vector field $g(\cdot)$, covariance matrix field $\Sigma(\cdot)$, and jump measure family $\varpi(\cdot,\cdot)$ are all independent of $n$; see also Examples 3.1 and 3.2.

Let us impose assumptions on $\mathcal{A}_{n}$, similar to Assumptions 3.1 and 3.2.

Assumption 3.4**.**

For every $n\in\mathbb{Z}_{+}$, the above (reflected) jump-diffusion in $D_{n}$ has a unique invariant distribution $\nu^{(n)}_{D_{n}}$, with corresponding boundary measures $\nu^{(n)}_{F^{(n)}_{i}},\,i=1,\ldots,m_{n}$.

Assumption 3.5**.**

There exists a finite measure $\upsilon$ on $D$ whose restriction $\nu^{(n)}_{D_{n}}$ on $D_{n}$ is a stationary measure for $\widetilde{Z}_{n}$, for every $n$.

This independence of the invariant measure $\upsilon$ of $n$ is similar to Assumption 2.1 in Section 3.

Assumption 3.6**.**

We have:

[TABLE]

Under these assumptions, we obtain the following theorem, analogous to Theorem 3.12.

Theorem 3.3**.**

Under Assumptions 3.3–3.6, the combined proces $(N,Z)$ with the generator $\mathcal{L}$ from (3.28) has a unique invariant probability distribution $\pi$ given by

[TABLE]

The corresponding boundary measures $\nu_{F_{i}}$ for $F_{i}$ (if there is reflection) are given by

[TABLE]

Finally, this process is ergodic in the sense of (3.12).

Proof.

For the proof of the stationary measure, we proceed very similarly to the proof of Theorem 3.12, except that we change (3.14)

[TABLE]

To prove ergodicity as in (3.12), similarly to Theorem 3.12, we show an analogue of Lemma 3.19:

Lemma 3.4**.**

For all $n,m\in\mathbb{Z}_{+}$, $z\in D$, and a subset $G\subseteq\mathbb{Z}_{+}\times D$ of positive Lebesgue measure:

[TABLE]

Proof.

Similarly to Lemma 3.36, without loss of generality, assume $G=\{m\}\times E$ for a subset $E\subseteq D$ of positive Lebesgue measure, and $m\geq n$.

Case (a). $z\in D_{n},\,E\subseteq D_{m}$. We prove this statement similarly to Lemma 3.19, using induction over $m$. Induction base ($m=n$): can be shown as in (3.20). Induction step: for $m=n+1$ we prove this as in (3.23) (using the same notation), but we integrate over $D_{n}\cap D_{n+1}$ instead of $D$:

[TABLE]

Assuming we proved this for $m=n+k$, let us prove this for $m=n+k+1$. Similarly to (3.24), but integrating over $D_{n+k}\cap D_{n+k+1}$, we get:

[TABLE]

This completes the proof of the induction step, and with it the proof of (3.36) in case (a).

Case (b). $z\in D_{n},\,E\cap D_{m}=\varnothing$. (Clearly, we can reduce the case of a general $E$ to these two cases (a) and (b).) Since $E\subseteq D=\cup_{k}D_{k}$, there exists a $k$ such that $E\cap D_{k}$ has positive Lebesgue measure. Take the $k$ with such a property which is closest to $m$. The process can get from $(n,z)$ to $\{k\}\times(E\cap D_{k})$ with positive probability in time $t/2$, using the path described in case (a) above. Afterwards, for every $z\in E\cap D_{k}$, the process $(N,Z)$ can jump from $(k,z)$ to $(m,z)$ in time $t/2$ with positive probability. Indeed, for $l$ between $k$ and $m$ we have $z\notin D_{l}$; thus the component $N$ will jump from $k$ to $m$, and the environment component $Z$ will stay constant at $z$.

Case (c). $z\notin D_{n}$. There exists a $k$ such that $z\in D_{k}$, since $z\in D=\cup_{k}D_{k}$. Find such $k$ wich is closest to $n$. The process $(N,Z)$ can get from $(n,z)$ to $(k,z)$ in time $t/2$ with positive probability: The queue component $N$ will jump from $n$ to $k$, and the environment component $Z$ will stay constant at $z$, since $z\notin D_{l}$ for $l$ between $n$ and $k$. Starting the process from $(k,z)$ instead of $(n,z)$ now, we are back to cases (a) and (b). Applying results from these cases for $t/2$ instead of $t$, we prove (3.36) for $z\notin D_{n}$. ∎

We proved Lemma 3.36, and with it we proved ergodicity (3.12), and thus Theorem 3.3. ∎

Now we provide examples in which the generator of the diffusive component in the joint process depends on the queueing state in a nontrivial manner.

Example 3.1**.**

Assume $D=[0,1]$, $D_{n}=[0,\alpha_{n}]$, $\lambda(z)=z$, $\mu(z)=1$. Assume ${\mathcal{A}}_{n}$ is a reflected diffusion (without jumps):

[TABLE]

The functions $a,b\in C^{2}_{b}([0,1])$ are given, describing the local diffusion coefficient and the local drift of the processes ${\widetilde{Z}}_{n}$ in $D_{n}$, with $a(z)>0$ for $z\in(0,1)$. Standard formulas from [8] guarantee that the measure $\upsilon$ on $D$ has Lebesgue density

[TABLE]

assuming that

[TABLE]

Assumption 3.6 becomes

[TABLE]

In particular, for $a(z)\equiv 1$ and $b(z)=\theta/(z-1)$ with $\theta>0$, we get $q(z)=2(1-z)^{-2\theta}$. If we choose

[TABLE]

for some $\alpha_{*}\in(0,1)$ and $n_{0}\in{\mathbb{N}}$, then (3.40) holds for $\theta\in(0,1/2)$. If $a\equiv 1$ and $b\equiv 0$, then the driving process for the environment is a reflected Brownian motion, with $\upsilon$ being the Lebesgue measure, and $q(z)\equiv 1$.

This example can be interpreted as follows. We keep the service rate fixed: $\mu=1$, while the arrival rate $\lambda$ varies as a reflected diffusion on $[0,1]$ if the queue size $n$ is less than an agreed threshold $n_{0}$. However, if $n$ reaches level $n_{0}$ while $\lambda<\alpha^{*}$, we allow $\lambda$ to vary only in a “safety range” $[0,\alpha^{*}]$. If $n$ attains level $n_{0}$ while $\lambda\geq\alpha^{*}$, we simply “freeze” $\lambda$ until the queue size becomes $n_{0}-1$, at which time $\lambda$ is again allowed to follow the diffusion on $[0,1]$.

Example 3.2**.**

Fix the arrival rate $\lambda=1$ while the service rate $\mu_{n}$ is subject to a reflected diffusion on the interval $D_{n}:=[\alpha_{n},\alpha^{*}]\subset[1,\alpha^{*}]=:D$ and kept unchanged in $D\setminus D_{n}$. Here $\alpha_{*}>1$ is a fixed constant. Here again, the generator ${\mathcal{A}}_{n}$ is given by (3.37), with $a,b\in{C}_{b}^{2}([1,\alpha^{*}])$; this operator from (3.37) acts on $f\in{C}^{2}([\alpha_{n},\alpha^{*}])$ with boundary conditions $f^{\prime}(\alpha_{n})=f^{\prime}(\alpha^{*})=0$. Instead of (3.38), we have

[TABLE]

and instead of assumption (3.39), we have

[TABLE]

Assumption 3.6 becomes

[TABLE]

As in Example 3.1, if $a\equiv 1$, $b\equiv 0$, then the driving process for the environment is a reflected Brownian motion, with $\upsilon$ being the Lebesgue measure, and $q(z)\equiv 1$.

4. Explicit rates of exponential convergence

4.1. A brief summary of results and methods

In this section, we prove (for both discrete-space and reflected diffusion environments) that for some constants $C,\varkappa>0$, we have

[TABLE]

and estimate the constant $\varkappa$. We do this by coupling: Take two copies $(N_{1},Z_{1})$ and $(N_{2},Z_{2})$ of this process starting from $x_{1}=(n_{1},z_{1})$ and $x_{2}=(n_{2},z_{2})$. Couple them (that is, construct them on the same probability space) such that the coupling time

[TABLE]

satisfies $\mathbb{E}\left[e^{\varkappa\tau}\right]<\infty$ for some constant $\varkappa>0$. By the standard Lindvall inequality we get

[TABLE]

We need only to integrate (4.2) with respect to $x_{2}\sim\pi$ to get (4.1). To obtain such a coupling, we apply the following method. We wait until the queue component hits [math] for both copies. Thus these queue components become coupled, that is, they are at the same point. Then we wait until: (a) either one of these queue components jumps back to $1$, or (b) the environment components become coupled. In case of (b), we have coupled both copies. In case of (a), we have failed, and need to repeat this procedure. Each time, we succeed with positive probability (bounded from below). Thus the number of tries is dominated by a geometric distribution.

To couple the environment components, we use the results of [37]; however, it is well-known how to find hitting time of zero by the $M/M/1$ queue [36]. Note that assuming exponential rates of convergence of $\mathcal{A}_{n}$ given each queue state $n$ does not immediately imply the exponential rate of convergence of the joint process $(N,Z)$. The particular multiplicative structure we consider in $\mathcal{A}_{n}$ enables us to obtain exponential estimates for the coupling time constructed for the joint processes $(N,Z)$ under the mild conditions imposed on $\mathcal{A}$ as well as the arrival and service rates.

4.2. Main statements

We impose two assumptions. The first assumes exponential bounds on the coupling time (uniform in state variables) associated with the generator $\mathcal{A}$.

Assumption 4.1**.**

The domain $D\subseteq\mathbb{R}^{d}$ is bounded. There exist constants $\alpha>1$ and $\gamma>0$ such that for all $z_{1},z_{2}\in D$ we can couple two processes $Z_{1},Z_{2}$ with generator $\mathcal{A}$, starting from $Z_{1}(0)=z_{1}$ and $Z_{2}(0)=z_{2}$, in time $\tau_{z_{1},z_{2}}:=\inf\{t\geq 0\mid Z_{1}(t)=Z_{2}(t)\}$, with

[TABLE]

The other assumption is a stronger condition on the traffic intensity: In previous sections, we assumed it is less than $1$, but now it has to be uniformly bounded away from $1$.

Assumption 4.2**.**

There exist constants $\overline{\lambda},\overline{\mu}>0$ which satisfy

[TABLE]

From this Assumption 4.2,

[TABLE]

Next, define the function

[TABLE]

This function is concave, increasing on $[1,c^{*}]$ and decreasing on $[c^{*},\infty)$, with $c^{*}:=\overline{\rho}^{-1/2}$, and

[TABLE]

Finally, define the function

[TABLE]

for any $\alpha>1$, $\beta,\gamma>0$ and $a\geq 0$.

Theorem 4.1**.**

Fix an initial condition $x_{0}=(n_{0},z_{0})\in\mathbb{Z}_{+}\times\mathbb{R}_{+}$. Under Assumptions 4.3 and 4.2, for some constants $C>0$ and $c\in(1,c_{*})$,

[TABLE]

where we can take any $\varkappa=(1-\varepsilon)m(c)$ for $\varepsilon\in(0,1)$ and $c\in(1,c^{*})$ such that

[TABLE]

The proof of the theorem is given at the end of this section. The only condition on the environment process is the Assumption 4.3 on coupling time with (uniformly) exponential tail for the environment process corresponding to $N(t)=0$. It is natural to assume this condition also holds for finite environment space.

Note that there exists a $c\in(1,c^{*})$ such that (4.8) is satisfied. Indeed, the left-hand side of (4.8) is continuous with respect to $c$, and is equal to $1$ for $c=1$. Whereas the right-hand side of (4.8) is larger than one for any $\epsilon\in(0,1)$. However, to find a maximal rate of convergence, one needs to maximize $\varkappa$ over the space of two parameters $(\varepsilon,c)$ which satisfy (4.8). Possible values of $\varkappa$ form an interval $[0,\varkappa_{*})$, which does not contain its upper endpoint; therefore, we cannot claim that $\varkappa_{*}$ is itself a rate of convergence.

Compare this with the simple $M/M/1$ queue with constant rates: arrival rate $\overline{\lambda}$ and service rate $\overline{\mu}$, which has an exact rate of convergence $e^{-m(c)t}$ for $c>0$ such that $m(c)>0$ from (4.5). [36, Proposition 5.8] states that the upper bound, restricting to only the queueing process, is

[TABLE]

The constant in the exponent does not depend on $n$. Our result matches this rate.

After some modifications, this theorem is applicable not only for reflected diffusions from Section 2, but for discrete environment space from Section 3. Here is its version:

Assumption 4.3**.**

There exist constants $\alpha>1$ and $\gamma>0$ such that for all $z_{1},z_{2}\in D$ we can couple two continuous-time Markov chains $Z_{1},Z_{2}$ with common generator $\sigma(\cdot)\mathbf{T}_{0}$, starting from $Z_{1}(0)=z_{1}$ and $Z_{2}(0)=z_{2}$, in time $\tau_{z_{1},z_{2}}$, such that (4.3) holds.

Theorem 4.2**.**

Under Assumptions 4.2 and 4.3, the result (4.7) for $(c,\varepsilon)$ satisfying (4.8) holds.

4.3. On the Assumptions 4.3 or 4.3

Below we give examples of discrete and continuous environment processes which satisfy Assumptions 4.3 or 4.3.

4.3.1. Coupling of jump processes

First, let us start with discrete-space Markov chains. The relation between coupling times and mixing times (for $P^{t}(x,\cdot)$ to converge within a fixed $\operatorname{TV}$ distance from the stationary distribution) is partially explored in [18]. There is a lot of existing literature on mixing times. For example, an extensive treatment of mixing times is given by [30]. The literature on coupling times is sparse. Much of the existing research is focused on $\gamma$ from Assumption 4.3, see for example [9], but we need to know both $\alpha$ and $\gamma$. We could not find articles which estimate both of them. Thus we present an elementary result, which we hope will be useful. The proof is in the Appendix.

Lemma 4.3**.**

Take a pure jump Markov process on the state space $D$ (finite, countable, or a domain in $\mathbb{R}^{d}$) such that the family of jump measures $(\nu(x,\cdot))_{x\in D}$ obeys

[TABLE]

and the family of probability measures

[TABLE]

satisfies the following condition:

[TABLE]

Then the coupling times $\tau_{x,y}$ satisfy the following uniform estimate:

[TABLE]

Remark 4.1**.**

The same result is true if the process is a reflected jump-diffusion with jump measures satisfying conditions of Lemma 4.3.

Example 4.1**.**

The condition (4.9) is not true if at least two measures $\overline{\nu}(x,\cdot)$ and $\overline{\nu}(y,\cdot)$ are mutually singular; that is, there exists a set $D_{0}\subseteq D$ such that $\overline{\nu}(x,D_{0})=0$ but $\overline{\nu}(y,D_{0})=1$. Indeed, we then have

[TABLE]

Example 4.2**.**

Assume that for all $x\in D$, $\nu(x,\cdot)\ll\mu(\cdot)$ for some $\sigma$-finite Borel measure $\mu$ on $D$. It can be the Lebesgue measure if $D$ is a domain in $\mathbb{R}^{d}$, or the counting measure for discrete $D$. Define the Radon-Nikodym derivative

[TABLE]

Then condition (4.9) is equivalent to

[TABLE]

For example, take a finite $D$ (with $m$ elements). Let $\mu$ be the counting measure, then $\nu(\cdot,\cdot)$ can be given by an $m\times m$ matrix $(\nu_{ij})$ (with zero diagonal elements). Each $i^{\rm th}$ row gives Radon-Nikodym derivative of $\nu(i,\cdot)$ with respect to $\mu$. Thus we obtain

[TABLE]

4.3.2. Coupling of reflected diffusions

Now consider a reflected diffusion on $[0,a]$. It is stochastically ordered, so every $\tau_{x,y}$ is stochastically dominated by $\mathcal{T}$: hitting time of $a$ starting from [math]. Thus

[TABLE]

Let us estimate the tail of $\mathcal{T}$. Take a non-reflected diffusion $Z^{*}=\{Z^{*}(t):\,t\geq 0\}$ on the real line, with drift and diffusion coefficients

[TABLE]

Let $\mathcal{T}^{*}:=\inf\{t\geq 0:|Z^{*}(t)|=a\}$. Then the laws of $Z(\cdot\wedge\mathcal{T})$ and $Z^{*}(\cdot\wedge\mathcal{T}^{*})$ are the same, and the laws of $\mathcal{T}$ and $\mathcal{T}^{*}$ are the same. Thus we have reduced this to tail estimation for an exit time of a diffusion process from a strip $[-a,a]$.

Denote by $u^{*}(t,x)$ the probability that $Z^{*}$ stays in $(-a,a)$ until at least time $t$, if $Z^{*}(0)=x$. Denote by $G(t,x,y)$ the transition density of this diffusion killed at $\pm a$, otherwise known as Green’s function (or heat kernel) of the infinitesimal generator $\mathcal{A}^{*}$ of $Z^{*}$. Then the function $u^{*}$ satisfies the initial-boundary value problem

[TABLE]

with initial and boundary conditions $\left.u^{*}\right|_{t=0}=1$ and $\left.u\right|_{x=\pm a}=0$. Thus we can express

[TABLE]

Knowing spectral decomposition of $G$ gives us the exponent in (4.3). To find the constant $A$ is a little harder, since it requires some information on the function $G$ itself, or its eigenvalues. In some simple cases, however, it can be found explicitly. For example, for a RBM $Z$ on $[0,a]$, the process $Z^{*}$ is also a Brownian motion, and [21, Chapter 2, Problem 8.2] gives us an exact estimate.

4.4. Proof of Theorem 4.2

We proceed in seven steps.

Step 1. It suffices to prove the following version of (4.2): For $(n_{1},z_{1}),(n_{2},z_{2})\in\mathbb{Z}_{+}\times D$,

[TABLE]

for some constant $C_{*}$ (which will be determined below). Indeed, then we can rewrite (4.10) as follows: For every Borel subset $A\subseteq\mathbb{Z}_{+}\times D$,

[TABLE]

Integrate (4.11) with respect to $(n_{2},z_{2})\sim\pi$. Note that the function $(n,z)\mapsto c^{n}$ is integrable with respect to $\pi$. Indeed, this integral is equal to

[TABLE]

From (4.4), $\nu_{D}(D)=1$, and $c<c_{*}=\overline{\rho}^{-1/2}$,

[TABLE]

Combining (4.11) and (4.12), we get (4.7).

Step 2. To get (4.10), we use coupling: As explained in the beginning of this section, we take on the same filtered probability space two copies $X_{1}=(N_{1},Z_{1})$ and $X_{2}=(N_{2},Z_{2})$ of this queue, starting from $x_{1}=(n_{1},z_{1})$ and $x_{2}=(n_{2},z_{2})$. Assume $\tau\equiv\tau(x_{1},x_{2})$ is a stopping time such that $X_{1}(t)=X_{2}(t)$ for $t\geq\tau$ a.s. Then $\tau$ is called a coupling time. For every $t\geq 0$ and a function $f:\mathbb{Z}_{+}\times D\to\mathbb{R}$ with $|f|\leq 1$, we can write

[TABLE]

In other words, we get the classic Lindvall inequality

[TABLE]

Next, assuming that we prove that $\mathbb{E}e^{\varkappa\tau}<\infty$, then

[TABLE]

Combining (4.14) with (4.15), we get (4.11). In the proof below, we shall see that the constant before $e^{-\varkappa t}$ turns out to be of the same form as required in (4.11).

Step 3. Let us now describe the coupling in detail.

(a) First, we couple the queue components. Both $N_{1}$ and $N_{2}$ are stochastically dominated by $\overline{N}$, which is described as the $M/M/1$ queue with arrival rate $\overline{\lambda}$ and service rate $\overline{\mu}$, starting from $\overline{N}(0)=n_{1}\vee n_{2}$. Therefore, we can take copies of $N_{1},N_{2},\overline{N}$ such that

[TABLE]

From (4.16) it follows that for $\tau_{0}:=\inf\{t\geq 0\mid\overline{N}(t)=0\}$, we have $N_{1}(\tau_{0})=N_{2}(\tau_{0})=0$.

(b) At $\tau_{0}$, we start two competing clocks. The first one is an exponential clock $\eta_{0}\sim\operatorname{Exp}(\overline{\lambda})$, which measures the time until arrival of the process $\overline{N}$ to $1$ from [math]. The second one is $\zeta_{0}$, a coupling time of $Z_{1}(\tau_{0}+\cdot)$ and $Z_{2}(\tau_{0}+\cdot)$. This time $\zeta_{0}$ exists by Assumption 4.3, since these two processes are copies of the environment process with generator $\mathcal{A}$ (recall $\beta_{0}=1$) starting from $Z_{1}(\tau_{0})$ and $Z_{2}(\tau_{0})$, respectively. At least (importantly for us here), this is true until $\eta_{0}$, when those drift and diffusion coefficients change.

(c) If $\zeta_{0}<\eta_{0}$, then $Z_{1}$ and $Z_{2}$ have time to couple while $\overline{N}(t)=0$. By stochastic domination, $N_{1}(t)=N_{2}(t)=0$. Thus $S_{0}:=\tau_{0}+\zeta_{0}$ is a coupling time for $X_{1}$ and $X_{2}$.

(d) If, however, $\zeta_{0}\geq\eta_{0}$, then the coupling did not work. The process $\overline{N}$ has jumped at time $\tau_{0}+\eta_{0}$ back to $1$, and we need to repeat this procedure. Let

[TABLE]

Let $\zeta_{1}$ be a coupling time of $Z_{1}(\tau_{1}+\tau_{0}+\eta_{0}+\cdot)$ and $Z_{2}(\tau_{1}+\tau_{0}+\eta_{0}+\cdot)$. If $\zeta_{1}<\eta_{1}$, then for $S_{1}:=\tau_{0}+\eta_{0}+\tau_{1}+\zeta_{1}$ we have $\overline{N}(S_{1})=0$, and thus $N_{1}(S_{1})=N_{2}(S_{1})=0$. But since $\zeta_{1}$ is also a coupling time for environment components, $Z_{1}(S_{1})=Z_{2}(S_{1})$. Thus $S_{1}$ is a coupling time for $(N_{1},Z_{1})$ and $(N_{2},Z_{2})$.

(e) If $\zeta_{1}\geq\eta_{1}$, then this coupling did not work, and we need to repeat this procedure, with $\zeta_{2},\eta_{2},S_{2}$, and so on. Let $\mathcal{J}:=\min\{j\geq 0\mid\zeta_{j}<\eta_{j}\}$. Then the ultimate coupling time is

[TABLE]

where we define the following random times:

[TABLE]

Next, we estimate the MGF of $\tau$ from (4.17).

Step 4. First, we estimate the MGF for each $\tau_{k}$. The generator of $\overline{N}$ is

[TABLE]

Therefore, letting $f(n)=c^{n}$ for a constant $c>1$, we get

[TABLE]

with the constant $m(c)$ defined in (4.5). The following process

[TABLE]

is a local supermartingale, because the function $W_{N}:n\mapsto c^{n}$ satisfies

[TABLE]

In the terminology of [37, Section 4], this is a modified Lyapunov function for $\overline{N}$. Then the derivation is similar to [37, Section 5]. By Fatou’s lemma, $L$ is a true supermartingale. Let

[TABLE]

Because $e^{ms}\geq 0$, this process is also a supermartingale. Consider the process

[TABLE]

By an elementary calculation, $\mathrm{d}L^{*}(t)=\mathrm{d}L_{*}(t)$. Therefore $L^{*}(t)=L_{*}(t)+\mathrm{const}$, and $L^{*}$ is itself a supermartingale. Thus, for every $t\geq 0$,

[TABLE]

Let $t\to\infty$ in (4.19). By Fatou’s lemma with the observation that $\overline{N}(\tau_{0})=0$, we get

[TABLE]

Similarly to (4.20), we get estimates for the MGFs of $\tau_{1},\tau_{2},\ldots$, with the difference that the initial state becomes $1$ instead of $n_{1}\vee n_{2}$. Therefore,

[TABLE]

Step 5. By Assumption 4.3, we have $P(\zeta_{k}>t)\leq\alpha e^{-\gamma t}$ for $t>0$, and recall that $\eta_{k}\sim\operatorname{Exp}(\bar{\lambda})$. Also, $\zeta_{k}$ and $\eta_{k}$ are independent. Thus, by Lemma 5.1, we have for all $k\in\mathbb{Z}_{+}$,

[TABLE]

Thus the number of ‘tries’, $\mathcal{J}$, is stochastically dominated by a geometric random variable $\widetilde{\mathcal{J}}$, which is the number of trials that one needs to get to the first success if the probability of success of each trial is $p$. It has the distribution and generating function (with $q:=1-p$)

[TABLE]

Step 6. Let us estimate the MGF of $\xi_{k}$, defined in (4.18). By Assumption 4.3 and Lemma 5.1 applied to $a:=m(c)$ for $c\in[1,c_{*}]$,

[TABLE]

The expression for $\theta(\alpha,\beta,\gamma,a)$ is given in (4.6). Combining (4.21) and (4.23), we get

[TABLE]

The same holds if we do conditional expectation

[TABLE]

Combining (4.20) and (4.23), we get

[TABLE]

Step 7. Finally, recall (4.17). We estimate from above the MGF for appropriate $\varkappa>0$: $\mathbb{E}\left[e^{\varkappa S_{\mathcal{J}}}\right]=\mathbb{E}\left[e^{\varkappa\tau}\right]$. By (4.24), the process $(M_{k})_{k\in\mathbb{Z}_{+}}$ defined by

[TABLE]

is an $(\mathcal{F}_{S_{k}})_{k\in\mathbb{Z}_{+}}$-supermartingale. It is positive, and $\mathcal{J}$ is an $(\mathcal{F}_{S_{k}})_{k\in\mathbb{Z}_{+}}$-stopping time. Applying the optional stopping theorem and using (4.25), we obtain

[TABLE]

By Hölder’s inequality,

[TABLE]

Since $\kappa(c)>0$ for $c\in[1,c_{*}]$, and $\mathcal{J}$ is stochastically dominated by a geometric random variable $\widetilde{\mathcal{J}}$ as in (4.22), we have

[TABLE]

Here we require that $\kappa(c)^{(1-\varepsilon)/\varepsilon}<q^{-1}$, which is exactly the condition for $c$ in (4.8). Combining (4.26),  (4.27) and  (4.28), we get

[TABLE]

From (4.17), this completes the proof of (4.10) for $\varkappa:=(1-\varepsilon)m(c)$, and Theorem 4.1.

5. Appendix

5.1. Proof of Lemma 4.3

Alternatively we can describe such pure jump process $X=(X(t),\,t\geq 0)$ as follows: Run an exponential clock $\eta_{1}\sim\operatorname{Exp}(\Lambda)$, and then let $X(t)=X(0)$ for $t<\eta_{1}$, and $X(\eta_{1})\sim\overline{\nu}(X(0),\cdot)$ (independently of $\eta_{1}$). Run another exponential clock $\eta_{2}\sim\operatorname{Exp}(\Lambda)$ independent of those random variables, then $X(S_{2}),Y(S_{2})$ with $S_{2}:=\eta_{1}+\eta_{2}$, and repeat the process. Thus we couple these processes $X=\{X(t):\,t\geq 0\}$ and $Y=\{Y(t):\,t\geq 0\}$ starting from $X(0)=x$ and $Y(0)=y$ as follows: We use the same exponential clocks $\eta_{1},\eta_{2},\ldots$, and couple $X(S_{k})$ and $Y(S_{k})$ with $S_{k}:=\eta_{1}+\ldots+\eta_{k}$, using the maximal coupling from [30, Proposition 4.7]:

[TABLE]

The coupling time then becomes

[TABLE]

Combining (4.9) and (5.1), we get

[TABLE]

Therefore, $\mathcal{J}$ is stochastically dominated by a geometric random variable $\tilde{\mathcal{J}}$ (the number of tries until the first success in a sequence of independent Bernoulli trials with individual success probability $p$), independent of $\eta_{1},\eta_{2},\ldots$. From (5.2) we get

[TABLE]

The MGF of each of these exponential random variables is

[TABLE]

and the generating function for this geometric random variable is

[TABLE]

Therefore, the MGF for $\tilde{S}$ from the right-hand side of (5.4) is the composition:

[TABLE]

Thus $\tilde{S}\sim\operatorname{Exp}(p\Lambda)$, and it satisfies $\mathbb{P}(\tilde{S}\geq t)\leq e^{-p\Lambda t}$. The rest is trivial.

5.2. A technical comparison lemma

Lemma 5.1**.**

Fix constants $\alpha>1$, $\beta,\gamma>0$. Take two independent random variables $\xi\sim\operatorname{Exp}(\beta)$ and $\eta>0$ which satisfies $\mathbb{P}(\eta>u)\leq\alpha e^{-\gamma u}$ for $u\geq 0$. Then

[TABLE]

For $a\in[0,\beta+\gamma)$, the moment generating function for $\xi\wedge\eta$ satisfies

[TABLE]

where the function $\theta$ is defined in (4.6).

Proof.

Let us first show (5.5). We have $\alpha e^{-\gamma u}<1$ for $u>u_{0}:=\gamma^{-1}\ln(\alpha)$. Then we can rewrite our tail estimate for $\eta$ as follows:

[TABLE]

Therefore, we have

[TABLE]

From here (5.5) immediately follows. Next, let us show (5.6). For every $u\geq 0$,

[TABLE]

Calculate the integral in the right-hand side by splitting it into two integrals: from [math] to $u_{0}$ (where $u_{0}$ is defined above), and from $u_{0}$ to $u$. If $u\in[0,u_{0}]$, this integral is equal to

[TABLE]

If $u>u_{0}$, then this integral is equal to

[TABLE]

Combining all these computations, we get

[TABLE]

Now integrate (5.7) with respect to the exponential distribution of $\xi$: $\beta e^{-\beta u}\,\mathrm{d}u$:

[TABLE]

This completes the proof. ∎

6. Concluding Remarks

We have found the explicit invariant measure for the joint interactive queueing and environment process, and estimated the exponential rate of convergence for the compact environment case. One interesting question would be to consider unbounded environment domains, but with environment process being exponentially ergodic. This will require much finer estimates, because Assumption 4.3 will hold only with $\alpha$ dependent on $z_{1}$ and $z_{2}$. One way to find such coupling was developed in [7, 20, 31, 37] via Lyapunov functions. Subgeometric rates of convergence seem interesting. Some work was done in [15, 29] for general Markov processes and in [1, 2] for some SDEs arising from many-server queues; but to the best of our knowledge none for our setup.

Acknowledgments

G. Pang was supported in part by NSF grants CMMI-1635410, and DMS/CMMI-1715875 and in part by an Army Research Office grant W911NF-17-1-0019. Y. Belopolskaya was supported in part by RSF 17-11-01136. Y. Suhov thanks Department of Mathematics at Pennsylvania State university for hospitality and support.