{'fraction_non_alphanumeric': 0.061808684965178204, 'fraction_numerical': 0.01590024580090127, 'mean_word_length': 3.6875900144023044, 'pattern_counts': {'":': 0, '<': 19, '<?xml version=': 0, '>': 26, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Martingale approach for tail asymptotic problems in the generalized Jackson network

We study the tail asymptotic of the stationary joint queue length distribution for a generalized Jackson network (GJN for short), assuming its stability. For the two station case, this problem has been recently solved in the logarithmic sense for the marginal stationary distributions under the setting that arrival processes and service times are of phase-type. In this paper, we study similar tail asymptotic problems on the stationary distribution, but problems and assumptions are different. First, the asymptotics are studied not only for the marginal distribution but also the stationary probabilities of state sets of small volumes. Second, the interarrival and service times are generally distributed and light tailed, but of phase type in some cases. Third, we also study the case that there are more than two stations, although the asymptotic results are less complete. For them, we develop a martingale method, which has been recently applied to a single queue with many servers by the author.


INTRODUCTION
Asymptotic analyses have been actively studied in the recent queueing theory. This is because queueing models, particularly, queueing networks, become very complicated and their exact analyses are getting harder. We are interested in asymptotic analyses for large queues in a generalized Jackson network and aim to understand their asymptotic behaviors through its modeling primitives.
There are two different types of asymptotic analyses for large queues. One is for a given fixed model. Large deviations are typically studied for this. Another is to study them through an approximating model. For example, such a model is obtained as the limit of a sequence of models under heavy traffic by scaling time, space and/or modeling primitives. It is called heavy traffic approximation (e.g., see [20], [22]). Here, large queues are caused by heavy traffic. In this paper, we focus on the large deviations for a fixed model. Among them, we are particularly interested in the logarithmic tail asymptotics of the stationary distribution for a generalized Jackson network, GJN for short . This problem has been studied by the standard approach of large deviations, but the decay rates are hard to obtain analytically by using modeling primitives (e.g., see [14]). The author [17] recently studied it by a matrix analytic method, and derived the decay rates for the marginal stationary distributions in an arbitrary direction for a two-station GJN, assuming phase-type distributions for service times and arrival processes, called a phase-type setting. We aim to generalize this result under a more general setting by a different approach.
Let d 2 be the number of stations in the GJN. For d = 2, we relax the phase-type assumption, and consider the decay rates of the stationary probabilities for state sets of small volumes in addition to those of the marginal stationary distribution. For d 3, we derive upper and lower bounds for those decay rates.
Our basic idea is to simplify the derivation of those asymptotic results in such a way that they are obtained in a similar manner to a reflecting random walk on a multidimensional orthant. This simplification greatly benefits for analysis although the decay rate problems for the reflecting random walk have not been fully solved even for d = 2. To this end, we take an approach studied for a single queue with heterogeneous servers in [18], and modify it for a queueing network. In this approach, we first describe the GJN by a piecewise deterministic Markov process, PDMP for short. We then derive martingales for change of measures, and formulate the asymptotic problems under a new measure.
PDMP is a continuous time Markov process whose sample path is composed of two parts, a continuous part, which is deterministic, and a discontinuous part, called jumps, by which randomness is created. Thus, PDMP is particularly suitable for queueing models. However, jump instants are random, and state changes at them are complicated. Because of this, PDMP is hard for analysis. So, other methods have been employed in queueing theory. For example, the state space is discretized by using phase-type distributions, and a Markov chain is obtained. Then, matrix analysis is applicable. This phase-type approach is numerically powerful but analytically less explicit because of state space description. Furthermore, it is getting harder to apply as a queueing model becomes complicated like a queueing network. We will not use such a matrix analysis. Nevertheless, it turns out that the phase-type assumption is helpful in our asymptotic analysis in some cases.
Contrary to the analytical difficulty, the PDMP has a simple sample path. Its time evolution is easily presented by a stochastic integral equation using a test function, which maps the states of the PDMP to real values (see (2.7) below). In this stochastic equation, state changes at the jump instants cause difficulty for analysis as we mentioned above. Davis [9], who introduced PDMP, replaces those state changes at jump instants by a martingale and the so-called boundary condition on the test function.
However, it is not easy to find a good class of the test functions which characterize a distribution on the state space of the PDMP. The idea of [18] is to choose a smaller class of test functions to overcome those difficulties. We then have a semi-martingale which cannot characterize a distribution on the state space, but still retains full information to study large queues. Once the semi-martingale is obtained, we use the standard technique for change of measure through constructing an exponential martingale, called a multiplicative functional. In applying this martingale approach to the GJN, we need to know how the network model is changed under the new measure. Intuitively, some of its stations must be unstable for the tail asymptotic analysis to work. To study this instability problem, we will use the fact that the network structure is unchanged under the change of measure, and therefore the stability of each station is characterized by the traffic intensity at that station. These traffic intensities are obtained from the traffic equations, but they are non-linear because of unstable stations. Thus, the instability problem is not obvious. We challenge it, and find some sufficient conditions for the GJN to be partly unstable under the new measure, which depends on the choice of a martingale for change of measure.
This paper is made up of four sections. In Section 2, the GJN (generalized Jackson network) is described by a PDMP, and a martingale for change of measure is derived. This section also considers geometric interpretations of the stability condition of the GJN, and present main results for the asymptotic problems. Section 3 discusses the method of change of measure, and considers how the GJN is changed under a new measure. In Section 4, the main results are proved. For this, we first list major steps for deriving upper and lower bounds, then prepare several lemmas to complete the proofs.
In this paper, we will use real vectors in the following way. Column and row vectors and their dimensions are not specified as long as they can be identified in the context where they are used. Their inequality holds in entrywise; e k is the unit vector whose k-th entry is one while all the other entries vanish; 1 is the vector all of whose entries are one. The inner product ∑ i x i y i of vectors x, y of the same dimension is denoted by ⟨x, y⟩, and ∥x∥ = √ ⟨x, x⟩; x is said to be a unit direction vector if x 0 and ∥x∥ = 1. We denote the set of all unit direction vectors in R d + by − → U d . For x in a finite-dimensional vector space S and its subset B, we will use the convention that x + B = {x + y ∈ S; y ∈ B}. For a finite set A, its cardinality is denoted by |A|.
year. I have benefited from those conferences about not only academic collaborations but also personal interactions. This paper is one of such outcomes, and I am grateful to Tomasz's great academic contribution as well as his warm hospitality.

GENERALIZED JACKSON NETWORK
We are concerned with a queueing network which has a finite number of stations with single servers and single class of customers. At each station, there is an infinite buffer, exogenous customers arrive subject to a renewal process if any, and customers are served in first-come-first-served manner by independent and identically distributed service times. Furthermore, the renewal process and service times are independent of everything else. Customers who complete service at a station are independently routed to the next stations or leave the network according to a given probability. We refer to this queueing network as a GJN (generalized Jackson network).

2.1.
Notation and assumptions. Let us introduce the notation for a GJN. Let d be the total number of stations. We index stations by elements in J ≡ {1, 2, . . . , d}, and let E be the set of the stations which have exogenous arrivals. For each station i, let F e,i for i ∈ E be the interarrival time distribution of exogenous customers, and let F s,i for i ∈ J be the service time distribution. Let p ij be the probability that a customer completing service at station i is routed to station j for i, j ∈ J , where those customers leave the outside of the network with probability To exclude trivial cases, we assume that d × d matrix P ≡ {p ij ; i, j ∈ J } is strictly substochastic, and (d + 1) where p 00 = 0, and p 0i > 0 only if i ∈ E, where the value of p 0i is specified later. We call P a routing matrix, while P is called an over all routing matrix.
At time t, let L i (t) be the number of customers in station i ∈ J , and let R s,i (t) be the residual service time of a customer being served there if any, where we set up a new service time just after service completion, and this service time is unchanged as long as station i is empty. Thus, all R s,i (t) are always positive because of the right continuity, and R s,i (t−) vanishes only at service completion instants. For i ∈ E, let R e,i (t) be the residual time to the next exogenous arrival at station i.
Denote the vectors whose i-th entries are (t), respectively, and define X (t) and R(t) as where Z + and R + are the sets of all nonnegative integers and all nonnegative real numbers, respectively. As usual, we assume that X(t) is right-continuous and has left-hand limits. Let {F t ; t 0} be a filtration generated by histories of all the sample paths of X(·); then F t is right-continuous, and {X(t); t 0} is an F t -Markov process.
Let F e,i and F s,i be the moment generating functions, MGF for short, of the distributions F e,i and F s,i , respectively. We define β w,i and θ w,i for w = e, s as We will assume that β w,i > 0 and θ w,i = ∞ for all w = e, s and i. That is, all the distributions, F w,i , have light tail and their moment generating functions diverge at their singular points. These conditions are assumed for technical simplicity.
For some important cases, we have to restrict these distributions in the following class. A positive random variable T or its distribution is said to have a conditional MGF with a uniform bound if there is a function h of θ > 0 such that E(e θT ) < ∞ implies that as long as P(T > t) > 0. Obviously, if T is bounded, it satisfies the condition (2.1).
Another obvious example is an NBU distribution, which is characterized by P(T > s + t|T > s) P(T > t) for s, t > 0. An important class for our application is of phase-type, which is defined as where a is a finite-dimensional probability vector, and U is a defective transition rate matrix with the same dimension as a such that (−U ) −1 is finite.
LEMMA 2. 1. A phase-type distribution has a conditional MGF with a uniform bound.
P r o o f. Assume that F is given by (2.2). Let T be a random variable subject to F , and let b(s) = ae sU /(ae sU 1). Then b(s) is a probability vector, and and therefore which is finite as long as E(e θT ) = a(−U )(θI + U ) −1 1 is finite. Hence, we have (2.1) by letting h(θ) be the maximum of all the entries of the vector (−U )(θI + U ) −1 1.
Thus, we consider the tail asymptotic problem for the GJN assuming the distributions of T e,i , T s,j to have light tails, and, in some cases, we assume the following: (A1) All the T e,i for i ∈ E and T s,i for j ∈ J have conditional MGF with uniform bounds, that is, satisfy (2.1).
For convenience, we put λ i = 0 for i ∈ J \ E. We now put p 0i = λ i / ∑ j∈J λ j for i ∈ J in P . Let α i for i ∈ J be the solutions of the following traffic equation: It is easy to see that the solutions uniquely exist because the routing matrix P is strictly substochastic and P is irreducible. Let ρ i = α i E(T s,i ), and assume the stability condition, that is, In Section 2.5, we will consider the case where some of stations are unstable. This case occurs under change of measure, and ρ i is no longer a right traffic intensity.

Piecewise deterministic Markov process (PDMP).
In this paper, we consider {X(t); t 0} as a piecewise deterministic Markov process, PDMP for short, introduced by Davis [9]. PDMP is a Markov process with piecewise deterministic and continuously differentiable sample path and finitely many discontinuous epochs in each finite time interval. Its randomness arises at discontinuous epochs, which are uniquely determined by hitting times when the deterministic sample path gets into a specified state set. The set of those discontinuous epochs constitutes a counting process, and the piecewise deterministic sample path is randomly changed at those times. We here assume that there is no other discontinuous state change. This slightly changes the standard description of PDMP due to [9], but it is a matter of formulation since Davis' PDMP can be described by the present formulation as well.
We now introduce the notation to describe X(t) as a PDMP. Let N be a counting process for the expiring times of all the remaining times. That is, where ∆ is the difference operator such that ∆f (t) = f (t) − f (t−) for a function f which is right-continuous and has left-hand limits. Clearly, N counts all the discontinuous points of X(t). However, it may multiply counts at the same instant, and therefore ∆N (t) may be greater than one. To avoid this, we define a simplification of N as We then let t 0 = 0, and inductively define t n = inf{u > t n−1 ; ∆N * (u) > 0} for n = 1, 2, . . . Thus, t n is the n-th discontinuous epoch of X (t), and a stopping time with respect to F t . Between times t n−1 and t n , X(t) is linearly changing, so continuously differentiable in such a way that This differentiation can be described by an operator A on C 1 (S), which is the set of all continuously differentiable functions from S to R. Namely, A is defined as Since the PDMP is a strong Markov process, the conditional distribution of X(t n ) given F tn− is a function of X(t n −) for each n 1, which is characterized by the transition kernel K given by ∃i ∈ E, y e,i = 0 or ∃i ∈ J , z i 1, y s,i = 0.
This Γ is referred to as a terminal set, while K is referred to as a jump kernel.

2.3.
Martingale decomposition of the PDMP. From (2.5) and (2.6) and the counting process N * , we have a time evolution equation for f ∈ C 1 (S), that is, .
We refer to f as a test function, as is typically called. We apply the same martingale method as discussed in [18]. We repeat it here briefly for this paper to be self-contained. We first note that We define M (·) and D(·) as we have the following fact.

LEMMA 2.2. For the PDMP X(·), if the condition
We refer to (2.12) as a terminal condition following the terminology of [18].

Terminal condition for the GJN.
A key of our arguments is to find a set of test functions satisfying the terminal condition (2.12). For this, we mainly apply the following test function, parameterized by θ ∈ R d , (2.13) f θ (x) = e ⟨θ,z⟩−⟨γ e (θ),y e ⟩−⟨γ s (θ),y s ⟩ , x ≡ (z, y e , y s ) ∈ S, using some vector-valued functions γ e (θ) ∈ R E and γ s (θ) ∈ R J , where we recall that ⟨a, b⟩ is the inner product of vectors a, b of the same dimensions. In some cases, it needs to truncate some of y e,i and y s,j as y e,i ∧ v and y s, where y = (y e , y s ) and Our first task is to determine functions γ e,i (v, ·) and γ s, For this, we first consider a prototype for them as we have done in Section 2.3 of [18]. Let T be a positive-valued random variable, and denote its distribution by F . We truncate T by a positive number v as T (v) ≡ T ∧ v, and denote the distribution of T (v) by F (v) . We denote the moment generating functions of F and F (v) by F and F (v) , respectively.
Note that F (v) (θ) exists and is finite for all θ ∈ R, but this may not be true for Obviously, ξ uniquely exists for each θ ∈ R and v > 0. We denote it by ξ F (v, θ). It has the following properties, which are proved in Lemma 2.4 of [18].
We define ξ which exist and are finite. These functions have some nice properties. For them, we cite Lemma 2.5 of [18] in which θ = ∞ in the present case.
F (θ) has the following properties: Throughout the paper, we assume that which means that F e,i and F s,j have light tails, and their moment generating functions diverge at the upper boundaries of their convergence domains. This assumption can be removed by using ξ (△) F as shown in [18] for a single queue. However, it will be complicated for a queueing network. This is the reason why we assume (2.20). Let Then q i (θ) > 0 and is convex in θ, and it is easy to see that is a convex function of θ ∈ R d because ∑ j∈J p ij e θ j + p i0 is a sum of convex functions (see the lemma of [12]) . We now define γ e,i (v, θ i ) and γ s, As informally mentioned, we let γ e,i ( Due to the assumption (2.20), these functions are well defined for all θ ∈ R d . Clearly, their definitions are equivalent to These equations mean that ∆R e,i (t) and ∆R s,j (t) at the jump instants are compensated by the change of the queue lengths so that the terminal condition (2.12) is satisfied. This is an intuitive background for the definitions of γ e,i , γ s,j .
REMARK 2. 1. The reader may wonder why the minus signs are needed in (2.21) because γ e,i , γ s,j in the test functions f θ , and f θ also have the minus signs, and they can be cancelled. The reason for this is that they have nice interpretations for large deviations. For example, let N e,i (t) be the number of arrivals at station i ∈ E by time t. Then N e,i (·) is a renewal process, and Glynn and Whitt [10] show that for any initial distribution for N e,i (·). This suggests that γ e,i (θ i ) must be one of the key information for the tail asymptotic of our problem. However, we will not use this property of γ e,i because the definition (2.21) is sufficiently informative for our analysis.
Furthermore, γ J(v) (θ) converges to γ(θ) for each θ ∈ R d as v → ∞, which is uniform on a compact set of θ. The next lemma is a key for our arguments, and easily follows from Lemma 3.2 in [4]. We also remarked its intuitive meaning below (2.22). So far, its proof is omitted.
, the test function f J(v),θ of (2.14) satisfies the terminal condition (2.12) with equality for all θ ∈ R d , respectively.
We next consider a martingale for the test functions f J(v),θ . Denote the probability measure for X(·) with the initial state x ∈ S by P x , and let E x stand for the expectation under P x . We first note that always holds for each x ∈ S and θ ∈ R d because the total number of exogenous arrivals and service completions in each finite time interval has a super-light tail (lighter than any exponential decay); see, e.g., Lemma 4.1 of [18] for the single queue case. Hence, Lemmas 2.2 and 2.5 immediately imply the following fact.

Stability condition and geometric interpretation.
As we mentioned in Section 2.1, the GJN (generalized Jackson network) is stable if the stability condition (2.4) holds. Except for trivial cases, it is also necessary. We will consider this network under change of measure, which is generally unstable, and it is important to see under what condition which station is unstable. To make these arguments clear, we formally define stability and instability for each station. Station i is said to be weakly stable (stable) if L i (t) is recurrent (positive recurrent, respectively), and to be weakly unstable (unstable) if L i (t) is null recurrent or transient (transient, respectively).
In this subsection, we so far do not assume the stability condition (2.4), and consider conditions for a station to be unstable (or stable). For this, we first need to compute an arrival rate at each station, which is obtained as the maximal solution { α i ; i ∈ F} of the following traffic equation (e.g., see [5], [6]): It is easy to see that α i α i for all i ∈ J , where recall that α i is the solution of the standard traffic equation (2.3). The α i can be numerically obtained from (2.27) in finite steps, but it is hard to get its analytic expression. For us, it is particularly important to give sufficient conditions in terms of γ e,i , γ s,i for a station to be unstable or weakly unstable because these functions are well handled under change of measure. We first give sufficient conditions for instability in terms of λ i , µ i , p ij and α i .
Unfortunately, this is generally not true. A counterexample is easily obtained, for example, for a three-station tandem queue (see Section 4 of [5] for some related discussions). We need to update α j using the information on the unstable station to be available to get such a generalization, but it would be less analytically tractable. Thus, we will not pursue it in this paper.
P r o o f. (a) Since α j α j , it follows from α j µ j that ρ j 1. If (2.28) holds, (2.27) implies that α j µ j , which is equivalent to ρ j 1.
(b) Suppose that ρ j 1, contrary to the claim; then ρ k 1 for all k ∈ J by the second assumption. Hence, α k µ k for all k ∈ J , and therefore the nonlinear traffic equation (2.27) is identical with the linear traffic equation (2.3). Thus, α k = α k µ k for all k ∈ J . This contradicts the assumption that µ i < α i , and 27) can be written as Hence, (2.29) implies that We then sum up both sides of this inequality for all j ∈ A, which yields which contradicts the irreducibility of the over all routing matrix P , and therefore A = ∅. This proves the first half of (c). It also implies that α j = µ j for all j ∈ J . Hence, if (2.29) holds with strict inequality, then (2.27) implies that This proves the remaining part of (c).
We next characterize the conditions in Lemma 2.7 by the gradient vector of γ(θ) and γ s,i (θ) at θ = 0. Define the gradient operator ∇ as Note that this t i is uniquely determined except for its length ∥t i ∥.
Thus, µ k − α k and ⟨∇γ(0), t k ⟩ have the same sign, which proves the claim.  The left panel shows that n < 0 and ⟨n, t i ⟩ > 0 for i = 1, 2, and both stations are stable by Lemma 2. 9. The right panel shows that n 1 > 0, n 2 < 0, ⟨n, t 1 ⟩ < 0 and ⟨n, t 2 ⟩ > 0, and station 1 is unstable while station 2 is stable by Lemmas 2.7 and 2. 9. Lemma 2.9 together with Lemmas 2.7 and 2.8 provides us geometric interpretations of the stability and instability conditions for stations through curves γ(θ) = 0 and γ s,i (θ) = 0 for i ∈ J . We will use them for the GJN before and after change of measure.

Tail asymptotics.
We now return to the assumption that the GJN is stable, that is, (2.4) holds. Under this assumption, we will use the following sets for considering the tail asymptotics of the stationary distribution. Let For A ⊂ J , let Those sets are open and connected. We denote their boundaries by putting the is a non-empty bounded and convex set because γ(θ) is convex and diverges as ∥θ∥ goes to infinity in any direction, and therefore Γ OUT is also not empty. We check below that Γ IN A is not empty for A ̸ = ∅. LEMMA 2. 10. Assume that the GJN is stable, and let A ⊂ J . (b) Define We now present main results, which are proved in Section 4. 3. For this, we use the following notation. For x ∈ R d and A ⊂ J , let x A be the |A|-dimensional vector which is obtained from x by dropping its i-entry of x for all i ∈ J \ A. Let Note that Θ k = R 2 for d = 2, and therefore r * (e k ) = r {k} (e k ). For A ⊂ J and c ∈ − → U d , that is, the unit direction vector c, let Note that r A (c) m A (c) because ∥c∥ = 1.

define a convex corn as
For d = 2, we can get bounds explicitly. For this, let Then we have the following corollary.

COROLLARY 2.1. Assume the stable GJN has two stations (d = 2).
(a) For k = 1, 2, It is notable that D 2 have been obtained as the convergence domain of φ(θ) ≡ E(e ⟨θ,L⟩ ), and used to derive (2.47) for the two-station JGN with phase-type F e,i , F s,j in Theorem 4.2 of [17]. The asymptotic (2.45) in the coordinate directions is not derived in [17], but can be obtained from Theorem 3.2 there. We here have the asymptotic (2.45) without the phase-type assumption. We conjecture that the assumption (A1) can be removed from all the results, but it seems a quite hard problem.
For d 3, there is very little known about the tail asymptotics of the stationary distribution not only for the GJN but also a reflecting random walk and SRBM. There are some studies in the framework of sample path large deviations, but those results need to solve certain optimization problems, which are hard to solve even numerically (e.g., see [14]). Contrary to them, (2.38) and (2.39) may be applied to get explicit bounds, by using ideas for a reflecting random walk (see Theorem 6.1 of [16]).

CHANGE OF MEASURE FOR GJN
In this section, we present some preliminary results for proving Theorems 2.1 and 2.2 and Corollary 2.1. A change of measure is typically used in the theory of large deviations. We also use it, and construct a new measure using a multiplicative functional, which is obtained from the martingale M J(v),θ (·) in Section 2. 4. However, we assume J (v) = ∅ in this section for making arguments simpler. It also suffices for major applications in the later sections.
Thus, the new measure is constructed from M θ (·) ≡ M ∅,θ (·). For this, we first drive a multiplicative functional. Its derivation is rather standard, but will be detailed because it is crucial for our arguments. Our major interest in this section is to see how the GJN is modified under the new measure. It is important for us to specifically identify its modeling parameters, which has not been studied in the literature except for the single queue case (see [18]), and may have an independent interest.
We define the integral of Y (·) with respect to martingale M (·) by where integration is a natural extension of a Riemann-Stieltjes integral (see Section 4d of Chapter I of [11]). For a positive-valued test function f , choose Y (t) as which is obviously positive and continuous in t and adapted to F t . Hence, Y · M (·) is a martingale. We denote it by E f (·). Thus, it follows from (2.9) that which is an F t -martingale under P x . On the other hand, E f (·) is a multiplicative functional because it is right- Thus, we can define a probability measure P f x for an initial state x ∈ S by is a martingale (see [13] for details). We refer to (3.2) as exponential change of measure. Let P ν and P f ν be probability measures such that P ν (C) = ∫ S P x (C)ν(dx) and P f ν (C) = ∫ S P f x (C)ν(dx) for X(0) to have a probability distribution ν on S, (3.2) implies that, for a nonnegative F t -measurable random variable U with finite expectation, we have where E ν and E f ν represent the expectations concerning P ν and P f ν , respectively. Similarly, for conditional expectations, we have, for 0 s < t, One can easily check this equation from the definition of a conditional expectation (see, e.g., Section III.3 of [11]).
When f = f θ of (2.13) and M = M θ of (2.26) with

GJN under the new measure.
Let us consider how the GJN is modified under the new measure P (θ) x . A general principle for change of measure is considered for a PDMP in [19], but we need to compute specific modeling parameters. For this, we follow the method of [18] studied for a single queue with many heterogeneous servers. We modify it here for the GJN. Since the differential operator A is unchanged because it works on a deterministic part of the sample path of X(·), we only need to consider the jump kernel K. Denote it under P (θ) x by K (θ) . We first write E f θ (t) explicitly as where w(θ, y) = w ∅ (θ, y) (see (2.15)), that is, Our first task is to compute the distributions of T e,i , T s,j under P (θ) x . These distributions (moment generating functions) are denoted, respectively, by F Since F where E x just before time when external arrivals and service completion, respectively, at station i occur. Then, by Lemma 3.1, we have The jump kernel K is changed to K (θ) as where θ 0 = 0. Hence, the routing probability from station i to j under P (θ) ij ≡ e −θ i +θ j p ij /q i (θ). (3.11) Thus, the GJN keeps the same network structure under the new probability measure P (θ) x , but its modeling primitives, F e,i , F s,j and p ij , are changed to F s,j and p (θ) ij , respectively, which do not depend on the initial state x. Let which corresponds to q i (η) in the original model, where η ∈ R d is a variable. From this definition and (3.11), we have As in the original network model, we define γ s,j (η) as the unique solutions of the following equations: for v ∈ (0, ∞]. These definitions yield and define γ (θ) (η) as Let α (θ) = λ (θ) (I − P (θ) ) −1 , and α (θ) is defined similarly to α.
We immediately see from these formulas that Similarly to (2.33) and (2.24), we have , j ∈ J . (3.14) Hence, we can update Lemmas 2.7 and 2.9 in exactly the same way as for the network model under P j , respectively. The following lemma is almost immediate from (3.11) and (3.14), but will be useful to check the conditions in Lemma 2.9 in which µ k , α k are updated to µ (θ) k , α (θ) k , respectively. As t i of (2.34), we define t Hence, the following lemma can be obtained similarly to Lemma 2.8. LEMMA 3.2. Let T (θ) be the matrix whose i-th column is t (θ) i ; then T (θ) is non-singular and not positive, that is, t 4

. PROOFS
The goal of this section is to prove the theorems and their corollary. A main idea is to use the new measure introduced in Section 3.2 by appropriate choosing the parameter θ. Some of its arguments are parallel to those in Section 4 of [18], but we require more lemmas because the state space for the queue lengths is multidimensional. We start to represent the stationary tail probability under the new measure.

4.1.
A procedure for deriving tail asymptotics. Recall the notation w(θ, y), and, for Then, it follows from (3.3) and (3.5) that, for a given initial distribution ν, where we recall that f θ ( X(0) ) = e ⟨θ,L(0)⟩−w(θ,R(0)) . We take the initial distribution ν in the following way. Let A be the first exit from and return times to Let ν − A be the distribution of X(τ EX A −) given that X(0) is subject to the normalized stationary distribution limited on ∂ L A S ≡ ∂ A S L × S R . This ν − A is taken for ν in (4.1). Denote a random vector subject to the stationary distribution of X(t) by X ≡ (L, R e , R s ). Then, the cycle formula yields, for x > 0 and B(x) ⊂ S L \ ∂ A S L , We are interested here in the asymptotic of P ( L ∈ B(x) ) as x → ∞. We apply change of measure to (4.2). For this, let τ x be a stopping time which satisfies We do not make any extra assumption for this τ x in this section, but it will be either τ IN B(x) or τ IN C(x) for some C(x) ⊂ S L \ ∂ A S L in Sections 4.2 and 4. 3. In what follows, except for Steps 5 and 5 ′ , we fix A ⊂ J . Let Then it follows from (4.1) with ν = ν − A and (4.2) that We are now ready to consider the asymptotic of P ( L ∈ B(x) ) as x → ∞. For its upper bound, we take the following steps.
S t e p 1. Choose θ ∈ Γ IN A ; then γ(θ) < 0 and γ s,i (θ) > 0 for all i ∈ J \ A, and therefore In Steps 2-4 below, we keep the same θ ∈ Γ IN A . S t e p 2. Verify that there is a constant S t e p 4. Find finite real-valued functions a 0 (θ) and a 1 (θ) > 0 such that then e −⟨θ,L(τx)⟩ is bounded above by e a 0 (θ)−a 1 (θ)x .
S t e p 5. Derive an inequality from (4.4) using Steps 1-4. Take the logarithms of both sides of this inequality, and divide by x; then letting x → ∞ yields the upper bound −a 1 (θ) for each θ and A ⊂ J . We then take the infimum of −a 1 (θ) for all θ and A ⊂ J for which Steps 1-4 work well.
To derive the lower bounds, we will use the martingale M J(v),θ of (2.26) in Lemma 2.6 for change of measure, choosing the index set for truncation: ) for each fixed θ ∈ R d as and we choose a sufficiently large v such that γ e,i (v, θ i ) < 0 for all i ∈ J e (v) and γ s,j (v, θ) < 0 for all j ∈ J s (v), which is possible by Lemma 2.4 and the assumption (2.20). Then, w J(v) (θ, y) of (2.15) is bounded below for all y = (y e , y s ) ∈ R E × R d . Namely, Then, using the martingale (2.26), we see that (4.4) is changed as . Note that the first integration term with minus sign in the exponent of (4.9) is bounded below by zero by the choice of J (v). We now take the following steps for the lower bounds.
S t e p 1 ′ . Choose θ ∈ Γ OUT such that γ s,i (θ) < 0 for all i ∈ J \ A, which implies i ∈ J s (v) and, for sufficiently large v > 0, The lower bounds are only used for Theorem 2.2. Thus, A = J for general d and A = {k} for d = 2. We also keep the same θ in Steps 2 ′ -4 ′ below.
S t e p 2 ′ . Verify that there is a constant C 2 such that, if τ x < ∞, then This inequality can be considered opposite to (4.5) in Step 2 because of (4.8).
S t e p 3 ′ . Find finite-valued functions a 0 (θ), a 1 (θ) such that In this procedure, we first need to find appropriate B(x) and τ x so that (4.11) and (4.12) hold, then go through steps. Among them, (4.12) is technically most demanding.

Lemmas for tail asymptotics.
We start with the following lemma.
We evaluate the integrand in the last integral using change of measure by H(·) similar to L(·). Let and we choose a sufficiently large v such that γ e,i (v, θ i ) > 0 for all i ∈ J e (v) and γ s,j (v, θ) > 0 for all j ∈ J s (v), which is possible by the same reason as that used for (4.7).
For change of measure, we use the test function f J(v),θ of (2.14) and the martingale M J(v),θ of (2.26), where L(t) is replaced by H(t). Then, the exponential martingale E f (·) is obtained as and define the new measure P We choose this v; then it follows from its conditional expectation version (3.4) that We note here that the condition (4.14) implies (4.18) from which the last term in (4.17) is bounded by Hence, the last term in (4.15) is bounded by This proves the lemma.
In the proof of Lemma 4.1, the condition (4.14) may be weakened as long as (4.18) holds. However, we also require the conditions (4.3) and (4.6) for τ x = τ IN B(x) to get an upper bound. In the view of these conditions, (4.14) is close to be necessary.
To see this, let N e,i (·) be the counting process for the exogenous arrivals at station i, and let N d,j,i (·) be the counting process for the customers who completed service at station j and routed to station i. Then the Palm formulas for stationary point processes yield where N s,j,i (·) is the stationary counting process for the number of service completions routed to station i when the server at station j is always busy. Since This proves the claim that for i ∈ E, j ∈ J are independent and have super-exponential distributions, that is, their tails are asymptotically faster than any exponential function (see, e.g., Lemma 4.1 of [18]). We now prove (4.19). Note that its terms multiplied by γ e,i (θ i ) 0, which is equivalent to θ i 0, and γ s,i (θ) for i ∈ J \ {k} can be dropped to bound the second expectation term in (4.19) Thus, it follows from the equation in (4.19) that (4.21)  (t). Let E e,i and E d,j represent the expectations concerning the stationary embedded distributions at exogenous arrivals at station i and at departure instants at station j, respectively. Then, they are known as Palm distributions (e.g., see [3]), and obtained as . From a similar bound in (4.20), the inequality (4.22) is obtained by using Both can be proved similarly to Lemma 4.8 of [18] (see Section S2 of its supplemented version for its proof), which is originally in Lemma 4.2 of [21].  Then lim inf x→∞ P ∥ is sufficiently small and for j ̸ = k and j ∈ {1, 2}.
P r o o f. For notational symmetry, we consider only the case for k = 1. Clearly, the lemma is obtained if lim inf x→∞ P hard to see that this is proved if station 1 is unstable and station 2 is stable under P . By Lemma 2.9 and ( 3.12), this holds if Otherwise, assume that θ (θ (CP 1 ) ) ̸ = θ (1,max) , and let f be a function from R to R such that θ 2 = f (θ 1 ) is determined by γ s,2 (θ) = 0. We then observe that f (θ 1 ) is increasing convex in θ 1 , and its derivative is smaller than that of the curve γ(θ) = 0 at θ = θ (CP k ) because θ (CP k ) is only one cross point of those two curves for θ 1 > 0 and Γ IN is not empty. Again the tangent vector t (θ (CP 1 ) ) 1 < 0 by Lemma 3.2, and therefore the second inequality of (4.26) must hold.   We only need to verify Step 3, that is, E ν − i ( f θ ( X(0) )) < ∞ for all i ∈ A because Lemma 4.2 cannot be used. We use here the assumption (A1); then it is not hard to see that, for i ∈ A, φ i (θ) < ∞ implies E ν i ( f θ ( X(0) )) < ∞. The latter finiteness implies that E ν − i ( f θ ( X(0) )) < ∞, as shown in the proof of Lemma 4.2. Thus, Step 3 is verified, and (4.27) is obtained. Taking the minimum of the right-hand side of (4.27) for A ⊂ J and A ̸ = ∅, we obtain (2.38).
We next prove (2.39 Thus, we complete the proof by taking the minimum of the right-hand side of the above inequality over A ⊂ J \ ∅.  . We choose θ such that γ(θ) > 0 and γ s,2 (θ) < 0; then Step 1 ′ works.
Step 2 ′ is obviously verified because Y (τ x ) does not decrease as x gets large.
Step 3 ′ is also obvious because B 0 is compact. For Step 4 ′ , we can take any bounded set for U ⊂ ∂ A S L . Then, if we take θ which is sufficiently close to θ (CP k ) , then (4.12) holds by Lemma 4.3, while (4.13) obviously holds. This completes the procedure, and (2.40) is obtained.
(b) We restrict the initial state in a bounded set C so that C ⊂ ∂ J S L × S R and Let c ∈ Corn( ← − Γ IN ∩∂Γ IN ), which implies ∇g(uc) 0 for uc ∈ Corn( ← − Γ IN ∩∂Γ IN ) by the convexity of Γ IN . Let is choose z 0 such that x ≡ (z, y) ∈ C implies that max i∈J z i < z 0 . We let B(x) = {z ∈ Z d + ; ⟨c, z⟩ > x}, x > z 0 , and let τ x = τ IN B(x) . Then, Step 2 ′ is obviously valid. Because the initial state is in C, we have ux < ⟨uc, L(τ x )⟩ ux + 1.
Hence, the condition (4.11) in Step 3 ′ is satisfied for x z 0 . Furthermore, if we take θ = uc for the change of measure, then all the stations are weakly unstable by Lemma 4.4, which implies that (4.28) and therefore (4.12) is satisfied for A = J . Thus, all the steps work well, and the proof is completed.