{'fraction_non_alphanumeric': 0.05782357574166448, 'fraction_numerical': 0.01785245471252297, 'mean_word_length': 3.7603624244948968, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 35, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Heavy-traffic approximations for a layered network with limited resources

Motivated by a web-server model, we present a queueing network consisting of two layers. The first layer incorporates the arrival of customers at a network of two single-server nodes. We assume that the inter-arrival and the service times have general distributions. Customers are served according to their arrival order at each node and after finishing their service they can re-enter at nodes several times (as new customers) for new services. At the second layer, active servers act as jobs which are served by a single server working at speed one in a Processor-Sharing fashion. We further assume that the degree of resource sharing is limited by choice, leading to a Limited Processor-Sharing discipline. Our main result is a diffusion approximation for the process describing the number of customers in the system. Assuming a single bottleneck node and studying the system as it approaches heavy traffic, we prove a state-space collapse property. The key to derive this property is to study the model at the second layer and to prove a diffusion limit theorem, which yields an explicit approximation for the customers in the system.


Introduction
We consider a network with a two-layered architecture. The first layer models the processing of customers by a network of two nodes. Each node can have multiple (but finitely many) servers. Customers are served according to their order of arrival and after finishing their service, they can re-enter at nodes several times for new services. The servers of the first layer act as jobs in the second layer, where they are simultaneously served by a common server working at speed one according to Processor-Sharing (PS) with rates depending on the number of customers in the system. Our goal is to derive an explicit approximation of the process describing the number of customers in the system. We analyze the system as it approaches heavy traffic. Under the assumption that there is a single bottleneck, we derive explicit results for the joint distribution of the number of customers in the system by proving a diffusion limit theorem. To achieve this, we look at the system in the second layer. In this way, we can aggregate the whole system since the total workload of the system (including the future workload due to customers re-entering the queues) acts as if were that of a single server queue with two independent renewal inputs.
To derive our diffusion limit theorem, we carry out a program inspired by the work of Barmson [1] and Williams [14], which consists of two main steps. First, we consider a critical fluid model, which can be thought of as a formal law of large numbers approximation under appropriate scaling. We identify the invariant states for the critical fluid model and we study the convergence to equilibrium of critical fluid model solutions as time goes to infinity. Our analysis has some similarities with the head-of-the-line processor sharing discipline as studied in [1], but there are differences. In particular, as the degree of resource sharing at each node is finite in our case, we need to define appropriate spatial regions in which the fluid model solutions have qualitatively different behavior. Our main result is to show that a solution of the fluid model converges to equilibrium uniformly (in terms of the initial condition) on compact sets. To achieve this, we perform a time change that facilitates our analysis.
The second main step is to show a state-space collapse property for the joint queue length vector process in heavy traffic. For an appropriately defined sequence of stochastic processes, we show that the difference between this vector and an appropriate deterministic mapping of the one-dimensional total workload process vanishes. The latter process is shown to converge to a one-dimensional reflected Brownian motion.
Our work can be seen as a partial network extension of the limited processor sharing (LPS) queue, of which fluid, diffusion, and steady-state heavy traffic limit theorems have been derived [18]- [20]. In our model, we assume that the inter-arrival and the service times have general distributions, but we consider that only one customer at each node can receive service at any time. In case that the inter-arrival and the service times are exponential, the service discipline at each station becomes irrelevant. An extension in the direction of general service times, using processor sharing at each node would require measure valued processes and is beyond the scope of the present paper. A mostly heuristic description of the results in this paper has appeared in [11]. In the classical applied probability literature, a version of our model has been investigated in a steady-state setting using boundary value techniques [6]; the solution in that paper may be used for numerical purposes, and is complementary to our heavy traffic limit, which yields explicit formulae, both for time-dependent as well as steady-state results.
In addition, our work is a contribution to the performance analysis of layered queueing networks. These are queueing networks where some entities in the system have a dual role (e.g., servers become customers to a higher-layer). In such systems, the dynamics in layers are correlated and the service speeds vary over time. Layered queueing networks can be characterised by separate layers (see [8] and [15]) or simultaneous layers. In the first case, customers receive service with some delay. An application where layered networks with separate layers appear is the manufacturing systems e.g., [3] and [4]. On the other hand, in layered networks with simultaneous layers, customers receive service from the different layers simultaneously. Layered networks with simultaneous layers have applications in communications networks. An application example where layered networks with simultaneous layers (such as our model) appear naturally are web-based multi-tiered system architectures. In such environments, different applications compete for access to shared infrastructure resources, both at the software level (e.g., mutex and database locks, thread-pools) and at the hardware level (e.g., bandwidth, processing power, disk access). For background, see [9] and [10].
The paper is organized as follows. We provide a detailed model description in Section 2 and we introduce the systems dynamics. In Section 3, we derive the fluid model and analyze it under the assumption of a single bottleneck and heavy traffic in the network. As we see, the assumption of the single bottleneck allows us to prove a State Space Collapse (SSC) property. Then, we show that a fluid model solution converges to equilibrium uniformly on compact sets. The main result of this paper is contained in Section 4. Namely, we provide a diffusion limit theorem for the joint customer population process for this two-layered queueing network. First, we prove that the diffusion scaled total workload process converges in distribution to a reflected Brownian motion. This result together with results in [1], lead to the main theorem..

Model
We assume a network with two layers. In layer 1, there are 2 single-server nodes indexed by i. Customers arrive at node i ∈ {1, 2} randomly one by one and have a random service requirement. A customer completing service at node i may be routed at node l, l ∈ {1, 2} for another service. It is assumed that customers are served according to their arrival order at each node; i.e., First In First Out. Only the first customer at each node can receive service at any time; i.e., the network is a Head of the Line network (HL).
In layer 2, there is a single server working at speed one. The servers of layer 1 are served by this single server simultaneously and at a rate which depends on the number of customers in the system. The model is illustrated in the following figure. In Section 2.1, we give a formal description of the model and in Section 2.2, we introduce the dynamics describing the model. In the sequel, we use the subscript i to refer to processes or quantities pertaining to each node and by convention, we omit the subscript to denote the 2-dimensional vector of these processes or quantities.

Preliminaries and model description
In this section, we give a formal model description. Let (Ω, F, P) be a probability space. For T > 0, let D[0, T ] 2 be the Skorokhod space; i.e., the space of 2-dimensional real-valued functions on [0, T ] that are right continuous with left limits endowed with the J 1 topology (as all candidate limit objects we consider are continuous, we actually only need to work with the uniform topology); cf. [2]. We denote by B(D[0, T ] 2 ) the Borel σ−algebra of D[0, T ] 2 . All the processes are defined from (Ω, F, P) to . For a process X(t), we denote the uniform norm by X(t) T = sup 0≤t≤T |X(t)|, where |X(·)| = max i |X i (·)|. We adopt the convention that all mentioned vectors are 2-dimensional columns and use a T to denote the transpose of a vector or a matrix a. We use A −1 to denote the inverse of a square matrix A, A k its k th power, and A the maximum element of A. Furthermore, I represents the identity matrix and e and e 0 are the vectors consisting of 1's and 0's, respectively, the dimensions of which are clear from the context. Also, e i is the vector whose i th element is 1 and the rest are all 0. Last, for a real number x, its integer part is represented by [x].
We start by describing the first layer. Let u i (j) for j = 2, 3, . . ., be the time between the (j − 1) th and j th external arrival at node i and u i (1) > 0 be the residual arrival time of the first customer entering at node i after time 0. We assume that the sequence {u i (j)} for i = 1, 2 and j = 2, 3, . . . is a sequence of positive i.i.d. random variables with mean 1/λ i , λ i > 0, and that u i (1) is independent of this sequence but sampled from an arbitrary distribution with the same mean. For i = 1, 2, define the cumulative arrival time process U i (·), as follows: U i (0) := 0 and U i (m) := m j=1 u i (j) for m ∈ N. The number of the external arrivals at node i until time t > 0 is given by the external arrival process In order to be able to define the total workload process in the system, including future service requirements due to routing, we need to introduce a sequence of random variables for any customer j. For any fixed time t ≥ 0 and j > 1, let v (1) ii (j), be the immediate service requirement of the j th customer (external or routed) at node i. Also, we define v (k+1) il (j) to be the service requirement of the j th customer (external or routed) at node i at the k th future time he visits node l, for i, l ∈ {1, 2} and j, k ∈ N. The sequence {v (k) il (j)}, indexed by j, is a sequence of i.i.d. random variables for any fixed i, l, k and for j > 1 and has mean β l := 1/µ l , µ l > 0. The random variable v (1) ii (1) denotes the residual service time for the first customer being served at node i at time 0; it is independent of the sequence {v (k) il (j)} but sampled from an arbitrary distribution with the same mean. In addition, we assume that all the above-mentioned random variables have finite second moments (more precisely, we need a Lindeberg-type condition to hold to make sure that the exogenous input processes satisfy a functional central limit theorem; see Section 4 for more details). We define the cumulative service time process as V i (0) := 0 and for m ∈ N, and the counting process We shall use the random variables {v (k+1) il (j)} to count the future workload in the system at time t. For a fixed t > 0, v (k+1) il (j) represents the k th future service requirement of the j th customer waiting to being served at node i and routed at node l. This event will occur after time t and after the completion of his service at node i.
Customers can move between queues according to Markovian routing. To describe the routing process, we define the following quantities. Let P be the (square) routing matrix of dimension 2. It is assumed that it is substochastic with a spectral radius less than one; i.e., its largest eigenvalue is less than one. In other words, the network is open, so the following relations hold: For any customer j at node i (external or routed), we define the random variables ϕ (k) il (j) = 1 if the j th departing customer from node i is routed to node l in k steps. The probability of this event is given by il denotes the (i, l) th element of the matrix P k . For i = 1, 2, we define the 2-dimensional random vector ϕ Note that ϕ . Now, we can define the routing process, which counts the number of customers who are routed from node l to node i as The total arrival rate at node i, γ i , is given by the solution of the following traffic equations In vector form, this can be written as It is shown in [2,Theorem 7.3] that under the assumptions described above, (2.2) has a unique solution γ = (γ 1 , γ 2 ) T . The traffic intensity of node i is ρ i := γ i /µ i . Now, we describe the service discipline at the second layer. Here, there is a single server. The servers of layer 1 become jobs at layer 2 in the sense they are served by the server of layer 2 simultaneously and at a rate that depends on the number of customers in layer 1 (at any time). The rate that each node receives is given by the service allocation function R(·) : R 2 + → R 2 + , with R(·) := (R 1 (·), R 2 (·)) T and for i = 1, 2, 3) The quantity q i represents the number of customers at node i. The 2-dimensional vector K = (K 1 , K 2 ) T is constant and we call it the degree of resource sharing. We assume that it is always finite and the user can choose it as a parameter of the system. Observe that K in this service allocation function guarantees a minimum service rate for each customer in the system. Also, note that the above function is Lipschitz continuous for q = (0, 0). We make the additional assumption that there exists a unique bottleneck in our system, which w.l.o.g. we let it be node 1. The definition of bottleneck is given below. By the previous definition, a straightforward inequality follows Observe that, if K 1 = K 2 , an intuitive explanation of the above definition is that the average occupancy of the server at node 1 is strictly greater than the server at node 2. In case of multi-server nodes where K i represents the number of servers an node i, the fraction ρi Ki is the average occupancy of a server at node i.

System dynamics
In this section, we introduce the dynamics that describe our model. We denote by Q i (t) the number of customers at node i at time t. This is given by where Q i (0) denotes the number of customers initially at node i. We define the cumulative service time of the server at node i as This quantity can be viewed as the effort that the server of node i has put in processing customers during [0, t]. Note that as the allocation function might be less than one, the above process is not necessarily equal to the amount of time that the server at node i is busy during [0, t]. In case the other node is empty during [0, t], (2.6) coincides with the busy time at node i. Recall that E i (t) is the number of external arrivals at node i up to time t. Observe that S i T i (t) , which is a composition of the renewal process (2.1) and the process T i (t), represents the number of departures at node i until time t. Furthermore, the total arrival process is given by The amount of time that both servers at the nodes are idle during [0, t] is given by the 1-dimensional process Alternatively, we can see this quantity as the idle time of the server in layer 2 during [0, t]. Further, the immediate workload at node i at time t is defined as Observe that W i (t) is nonnegative for any t ≥ 0. Last, due to the work-conserving property in layer 2, the following relation holds: Recall that when we omit the subscript i, we refer to the 2-dimensional column vector of the corresponding process/quantity; for example A(·) = (A 1 (·), A 2 (·)) T and W (·) = (W 1 (·), W 2 (·)) T . All the essential information of the evolution of the system is contained is the following 6-tuple In addition, the total (immediate and future) workload of the system plays a key role in our analysis. First, we define the remaining service requirement of the j th customer waiting to be served at node i = 1, 2 as is the future service requirement of the above-mentioned customer. Observe that, for an external arrival, s i (j) is the total service requirement. The first and the second moments of (2.11) are given (in vector form) by and Now, we can define the (1-dimensional) total workload of the system as (2.14) In case that , there are no customers at node i, we understand the second sum of the last equation as zero. Obviously, the total workload is not a Markov process as it is dependent on future service requirements. In Section 4.1, we shall see that under an appropriate scaling (i.e., the diffusion scaling) the dependence of the total workload on the future vanishes. Last, as our network is HL, only one customer can be in service at node i at any time. This property gives an upper and a lower bound for the cumulative service time (2.6) at node i, which is given in [1, We have so far defined the system dynamics for the above-mentioned two layered network and stated all the assumptions we need for our analysis. We are now ready to study the fluid model of this network, which is the first essential step to show a SSC property.

Fluid analysis
In this section, we study a critical fluid model, which is a deterministic model and can be thought of as a formal law of large numbers approximation under appropriate scaling. We shall give a rigorous proof of the last statement in the next section.
The main goal is to prove uniform convergence (w.r.t. the initial condition) on compact sets for the fluid model under the critical loading assumption; i.e., the traffic intensity of the network is one. First, we find the invariant points (or equilibrium states) and define an appropriate lifting map which describes these points. Then, we define a time-changed version of the original fluid model and we show that it is enough to prove the convergence for the time-changed function. As the time-changed function is given by a piece-wise linear ODE, we are able to find the solution and to show the convergence. Because the degree of resource sharing (the vector K) is finite we need to separate the state space in suitable regions and to distinguish cases depending on initial conditions.

Definition and invariant points
The traffic intensity of the network is given by ρ : We make the critical loading assumption, i.e., To derive the fluid model equations we replace any random quantity in (2.5)-(2.10) with its mean. The fluid model equations are given bȳ We can show that the immediate workload in the fluid model can be written asW i (t) = β iQi (t).

Definition 3.1 (Fluid model). We say that a 2-dimensional vectorQ(·) with non-negative components is a solution of the fluid model if it is continuous and satisfies
We define an auxiliary quantity, which can be interpreted as the total workload in the fluid model. It is defined by functionQ(·) as follows: A useful result in our analysis is that the fluid total workload in the system remains constant under the critical loading assumption.
Proposition 3.1. For any fluid model solutionQ(·), we have that Proof. IfQ(0) = 0, thenW T ot (t) ≡ 0, for t ≥ 0. LetQ(0) > 0. Assume nowQ(0) = 0. By definition,Q(t) is continuous, so let t be such thatQ(t) is positive in a neighborhood of t. Calculating the derivative of the total workload at time t, we derivē The last equation holds due to (3.1) and the property of the service allocation function; i.e.
Combining this with the continuity ofQ(t), we see that, necessarily, δ = ∞. Thus, the result extends to all positive t.
In the following lemma, we show that there exists a solution to the fluid model equations for all non-zero initial states and it is unique. Proof. LetQ(0) > 0. A by-product of the previous lemma is that δ = ∞. Thus, we can discard the origin and define the function Ψ(·) : where µ • R(·) indicates the Hadamard product; i.e., µ • R(·) = (µ 1 R 1 (·), µ 2 R 2 (·)) T . This function is Lipschitz continuous because R(·) is. Now, note that (3.2) can be written as where the prime denotes the derivative with respect to time. The existence and uniqueness follows directly by [12, Section 10, Theorem IV]. Now, we characterize the invariant points x ∈ R 2 + of the fluid model. Before we state our result, we first proceed in an informal manner. Equate the total rate into node i with the total rate out of node i. That is, Thus, we have that for the points on the invariant manifold (i.e. the set of the invariant points), ρ i = R i (x). Using the definition of a bottleneck and keeping in mind that we assume node 1 to be the bottleneck, we now describe the invariant points. We know by (2.4) that ρ1 K1 > ρ2 K2 , which yields Thus, by combining the last inequality and the definition of the service allocation function (2.3), we have that the following inequality holds The last inequality implies that for all invariant points x = (x 1 , x 2 ) of the fluid model, we have that The invariant manifold is thus given by We make the previous arguments rigorous by showing that a sufficient and necessary condition of the fluid queue length to remain constant in time is the initial state lies on the invariant manifold. Now, supposing thatQ(t) =Q(0) ∈ R 2 + for all t ≥ 0, then it follows from the previous discussion thatQ(0) ∈ I.
Having found the invariant (or equilibrium) points of the fluid model, we now turn to its stability property, namely the convergence of the solutions of fluid model equations to the invariant manifold as time goes to infinity.

Convergence to the invariant manifold for the fluid model
For this point, we define the critical workload as (cf. (3.1)) In order to prove a SSC property based on the critical workload level w * , we define a lifting map, ∆ : R + → R 2 + , as follows: Note that the lifting map is Lipschitz continuous with constant In the sequel, we show that fluid model solution converges to the invariant manifold as t goes to infinity. Theorem 3.4 (Convergence to the invariant manifold for the fluid model).
for t > t 0 and ∆W T ot (0) is an invariant state.
Sketch of proof. Here, we give a sketch of the proof. The complete proof is extended in the rest of this section. The first step is to define a function y(·) and a function G(·), and to show that y(·) can be interpreted as a time-change ofQ(t), namelyQ(t) = y(G(t)). Then, we show that the convergence of the time changed version implies the convergence of the original functionQ(·). To this end, let Ξ = {ξ ij }, i, j ∈ {1, 2} be the matrix and (3.14) We shall show that the above-defined function can be interpreted as a time-change ofQ(t). Let G(·) : Note that G(·) is continuous and that This means that the function G(·) is strictly increasing and unbounded in time, which implies that G(·) is also invertible. The original functionQ(t) can be interpreted asQ i (t) = y i (G(t)), for t > 0. To see this,Q where the function Ψ(·) is defined in (3.7). The idea now is to prove (3.13) by showing that y(t) converges as t → ∞. We will formally do so in the remainder of Section 3. Assuming that y(t) converges, we now show thatQ(t) converges. To see this, if for any > 0 there exists t 0 (independent of M ) such that where u = G(t). For u > G −1 (t 0 ), the last term becomes smaller that .
The remainder of the current section is devoted to the proof of (3.15). To do it, we first define appropriate spatial regions in which the fluid model solutions have qualitatively different behavior and we solve (3.14) in these regions. This is done in Section 3.3. In Section 3.4, these solutions are used to show that both y 1 (t) and y 2 (t) are monotone in t. A crucial observation is that we only need to look at y 2 (t), as τ 1 y 1 (t) + τ 2 y 2 (t) =W T ot (0) by Proposition 3.1. This paves the way for a global convergence analysis of y(t), also establishing the desired uniformity. This is done in Section 3.5.

Explicit local solutions of time-changed ODE
Because we assume that the degree of resource sharing K is finite, the form of (3.14) depends on the value of y(t). For this reason, we need to define the following regions. For x = (x 1 , x 2 ) ∈ R 2 + , we define The following picture shows these regions and the invariant manifold as defined in (3.10). Let y(0) = (y 1 (0), y 2 (0)) T ∈ R 2 + . We solve the time-changed ODE, which is given by (3.14), in regions Π 3 and Π 4 (considering these two regions only is sufficient for our purposes). It is useful to observe the relations between the coefficients of matrix Ξ in (3.14), which we use later. We know by the definition of the total arrival rate that The constant ξ 22 can be expressed as In a similar way, we can obtain By definition (2.12), τ 2 can be written as τ 2 (1 − p 22 ) = β 2 + p 21 τ 1 . Also, by (2.2) and (2.12), we have that τ T λ = β T (I − P T ) −1 λ = β T γ = 1. Combining the previous two equations, we get Last, by (3.17) and (3.18) we have that ξ 21 = − τ 1 τ 2 ξ 11 .

Solution in region Π 3
Assuming that y(s) is in region Π 3 for s ∈ [0, t], we can directly solve the second equation of system (3.14) since it is independent of y 1 (t): Then, by using (3.17), the solution is given by Now, we can easily obtain the solution of the first equation of system (3.14). By the relations between the coefficients of (3.14), this solution is given by
Using the relations between the coefficients of matrix Ξ, the solution of (3.20) is given by Having found the solution of (3.14) in each region Π 3 and Π 4 , we observe that the 2-dimensional equation can be reduced to a 1-dimensional equation sinceW T ot (0) = τ 1 y 1 (t) + τ 2 y 2 (t), for t ≥ 0. Now, it is enough to show the convergence of this equation. To see this, define x(·) as and its derivative as Using the observation that we can reduce the dimension by one and the solutions to system (3.14), we show that the fluid model solutions converge to an equilibrium state uniformly for all initial states within a compact set. First, we find the sign of the derivative of the above reduced equation. Then, as the limiting point depends on the sign of the quantityW T ot (0) − w * , we have to distinguish between the following three cases: the total workload in the fluid model is i) greater than, ii) less than or iii) equal to the critical workload.

Local analysis: establishing monotonicity
By (3.22), it is clear that we need to study only the behaviour of y 1 (t). In the sequel, we find the sign of (3.23) in each region Π i , for i = 1, . . . , 4.
In region Π 1 , we know that y 1 (t) ≤ K 1 . By (3.14) and (3.18), we have that , which is strictly positive by (2.4) and the fact that ξ 12 > 0. Therefor we conclude that the derivative of x(t) is strictly negative in region Π 1 .

Global analysis: convergence to Invariant Manifold
We are now ready to connect all pieces. From the previous section, we know that x(t) = y 2 (t) must be smaller than K 2 after a finite time, thus exiting regions Π 1 and Π 2 . Therefore, we can focus on the remaining two regions. In order to do so, we need to consider whether y 1 (t) will eventually be larger than, smaller than, or equal to K 1 . This leads to three cases, treated separately in the remainder of this section.
This means that the function x(t) lies in region Π 3 after that time. It is enough to prove that the equation y 1 (t) = K 1 , has a positive solution. Note that by (3.21), we have that Now, setting y 1 (t) = K 1 , the previous equation becomes We argue that (3.26) is satisfied by a t ≥ 0 as follows. Observe that by (3.11) the quantity ρ1 τ1ρ1+τ2ρ2 is equal to K1 w * , and by the assumption that y 1 (0) ∈ Π 4 , we have that y 1 (0) ≤ K 1 . Now, we can obtain the following inequality which proves the statement. Note that if the total workload is equal to the critical workload, then (3.26) does not have a positive solution since (3.26) would imply that exp{α 2 t} = 0. It would only be satisfied by t = 0. This means, that if y 1 (0) ∈ Π 4 andW T ot (0) − w * = 0, then the function x(t) remains in region Π 4 for ever. Since in this first caseW T ot (0) − w * > 0, we have thatW T ot (0) w * K 1 − K 1 > 0. Thus, by combining the fact that y 1 (0) ≤ K 1 and the previous display, we have shown thatW T ot (0) − w * K1 y 1 (0) > W T ot (0) − w * > 0. Combining these arguments leads to the conclusion that the equation y 1 (t * ) = K 1 has a (unique) positive solution, say t * . Therefore, y(t) ∈ Π 3 for t > t * . In other words, it is enough to prove that forW T ot (0) − w * > 0, the function x(t) converges to a point in region Π 3 .
We now show that it converges to the invariant manifold in region Π 3 . By (3.12), we have that ∆ 2WT ot (0) = ρ2 ρ1 K 1 . By (3.19) and (3.22), Case 2:W T ot (0) − w * < 0. Adapting the previous case, we first show that ifW T ot (0) − w * < 0 and y 1 (0) ∈ Π 4 , the the function x(t) remains for ever in region Π 4 . To see this, by (3.21) we have that Observing that y 1 (0) ≤ K 1 in Π 4 , we derive the following inequality is positive by the assumptionW T ot (0) − w * < 0 and that exp{α 2 t} ≤ 1. Combining these three facts, we have that y 1 (t) ≤ K 1 . Now, we show that if the process starts in region Π 3 andW T ot (0) − w * < 0, then there exists a finite time t * * such that x(t) ∈ Π 4 after that time. Again, here we prove that the equation y 1 (t * ) = K 1 has a positive solution. Then, the result follows by observing that y 1 (0) ≥ K 1 > ∆ 1WT ot (0), and in that case x(t) is an increasing function. By (3.19), we have that Setting y 1 (t) = K 1 , we obtain To show that the previous equation has a positive solution, it is enough to show that Recall that y 1 (0) ≥ K 1 and thatW T ot (0) − w * < 0. We can now derive the previous inequality by observing that Analogously with the previous case, we note that if the total workload is equal to the critical workload, and recall that α 2 < 0. The uniform convergence in a compact set for any initial state follows for the same reason as in previous case.
Case 3:W T ot (0) − w * = 0. In this case, the convergence follows from the comments we made in the previous two cases. IfW T ot (0) − w * = 0, then the function x(·) always stays in the region where y 1 (0) lies (see comments after (3.26) and (3.27)). As we see, the function converges in regions Π 3 and Π 4 .
This concludes the proof of Theorem 3.4, which will be applied to prove a diffusion theorem for the queue length process in the next section.

Diffusion approximations
The main objective in this section is to show a state-space collapse property (SSC) for the diffusion queue length process. This yields a diffusion limit theorem for the diffusion-scaled process. To do it, we follow the strategy set up in [1]. Let us consider a family of single-server systems indexed by n ∈ N, where n tends to infinity, with the same basic structure as that of the network described in Section 2. To indicate the position in the sequence of networks, a superscript n will be appended to the network parameters and processes. Diffusion (or central limit theorem) scaling is indicated by placing a hat over a process. Thus, the well-known diffusion scaling is given by X n (·) = 1 n X n (n 2 ·). Let ρ n = 2 i=1 ρ n i = (γ n ) T β. We set γ n = (I − P T ) −1 λ n , λ n = λ(1 − θ n ), µ n ≡ µ, P n ≡ P , and K n i = nK i , where θ is a positive real number. Thus, we have that ρ n = 1 − θ n . It is clear that under the critical loading assumption, λ n → λ and n(1 − ρ n ) → θ, as n → ∞. These are our heavy traffic assumptions. Furthermore, we assume that Q n (0) = 1 n Q n (0) →Q(0), whereQ(0) is a positive constant. The service allocation function for the n th model is given by where we observe that R 1 i (·) = R i (·). Recall that R(·) is a Lipschitz-continuous function on R 2 + \{0} and observe that the following scaling property holds: R n (n·) = R(·). In the sequel, we state the technical assumptions that allow us to apply the functional central limit theorem and Bramson's weak law estimates. We assume that for i = 1, 2, in probability as n → ∞. In addition, we assume that there exists a function η(·), with zero limit at infinity such that for j > 1, i, l ∈ {1, 2} and k ∈ N, More details about these assumption can be found in [1] and [17]. In the rest of this section, we assume that the previous assumptions are satisfied without refer to them again. The main result in this section is the following Theorem 4.1. Assume that the diffusion-scaled initial state converges in distribution as n → ∞, i.e.
denotes convergence in distribution. Then, the diffusion-scaled stochastic process converges in distribution as n → ∞, i.e.
where W T ot (t) is a 1-dimensional Brownian motion with drift −θ and variance si denotes the variance of s i (j) and c ui denotes the coefficient of variation of u i (j).
The proof of this theorem is given in the end of this section. This theorem can be applied to develop heavy traffic approximations for the joint queue length process, as they are both a piece-wise linear function of a one-dimensional RBM, of which the time-dependent distribution can be expressed in closed form (in terms of the Gaussian cdf and pdf); cf. [2]. To see this, note that the functions ∆ 1 and ∆ 2 are invertible with inverses We know that W T ot (t) is a RBM(−θ, σ 2 ). Let Q(·) be the diffusion limit. We have that for x, y ≥ 0, where z = max{∆ −1 1 (x), ∆ −1 2 (y)}. The last expression can be written in terms of the Gaussian cdf and pdf; cf. [2]. Also, using a similar coupling argument as in [20], it can be shown that one can interchange the steady-state and heavy traffic limits in this case. For space considerations we will leave this as detail to the reader.
The rest of this section is devoted to a proof of Theorem 4.1. It is organized as follows.
1. We first prove a heavy traffic limit theorem for the total workload process.
2. After that, we define a family of shifted fluid-scaled processes in Section 4.2 and we show that they are stochastically bounded in Section 4.3.
3. In Section 4.4, we establish some technical auxiliary estimates and tightness of these families. Moreover, we establish that limit points of these fluid scaled processes, which are called fluid limits, are in fact fluid model solutions as defined in Section 3. The development in this section is very similar to those in Bramson [1] and is therefore kept concise.
4. In Section 4.5, we establish a similar tightness property for a family of shifted fluid-scaled workload processes.
5. The proof is then completed by showing a state-space collapse result in Section 4.6.

Convergence of the total workload
Lemma 4.2. Under the critical loading assumption, the diffusion-scaled total workload, W n T ot (t) = 1 n W n T ot (n 2 t), converges in distribution to a RBM(−θ, σ 2 ).
Proof. By (2.11) for i, l = 1, 2 and i = l, we have that the total service requirement of the j th external customer (including customers who already are in queue i at time zero) who enter at queue i is given by We define the following process We recall that E i (·) denotes the external arrivals at queue i and by construction of the model, is a sequence of positive i.i.d. random variables for i = 1, 2. The process given in (4.1) represents the workload of a single queue with input given by two independent renewal process that have independent service requirements from each other. The busy time of this system is t 0 1 {W G (s)>0} ds = T 1 (t) + T 2 (t), and it represents the busy time of server in layer 2. Note that the busy time is zero if and only if both queues are empty. The diffusion-scaled process (after subtracting and adding the means of the random quantities in (4.1)), is given by where Y n (t) = t 0 1 { W n G (s)=0} ds. By the time change theorem in [13] and the functional central limit theorem (also see [2, Theorem 6.8]), we have that W n G (·) d → W G (·), as n → ∞. Furthermore, the limit can be described as where W i (t), 1 ≤ i ≤ 4 are independent 1-dimensional standard Brownian motions and Y (t) can be increased only if W G (t) = 0. Thus, the process W G (t) satisfies a 1-dimensional Skorokhod problem. That is, W G (t) is a reflected Brownian motion starting at point 2 i=1 τ iQi (0) with drift −θ and variance σ 2 which is given by where σ 2 si denotes the variance of s i (j) and c ui denotes the coefficient of variation of u i (j). The second moment of the random variables s i (j) is given by (2.13). In case of Poisson external arrivals, this result is reduced to the well-known heavy-traffic limit (see e.g. [7,Theorem 2.3]). Now, we shall prove that We do it by showing that we can change the label of how we count the service requirements of the customers in the system. Counting the total service requirements of the external arrivals until time t is the same as counting the immediate and remaining service requirements of the total arrivals in the system until time t. Recall that s i (j) = v (1) ii (j) + s i (j) for i = 1, 2, where s i (j) are the future service requirements. If t = 0, then we have nothing to prove as E i (0) = A i (0) = 0. If t > 0 and Φ li S l T l (t) = 0 for l, i ∈ {1, 2}, then (4.2) holds as all the departures until time t leave the system and so s i (j) = 0 for j = 1, . . . , S i T i (t) and E i (t) = A i (t).
For the general case, we first assume that all customers are routed only one time until the time t. The right-hand side in (4.1) can be written as In order to separate the customers who depart form node i (if they are routed or leave the system), we define the following sets for i, l ∈ {1, 2} and i = l, A i = {j : 1 ≤ j ≤ S i T i (t) }, which includes all the customers who depart from node i. Customers can be routed at the same node and we denote this set by A i 1 = {j ∈ A i : ϕ ii (j) = 1} or to the other node and we denote this set by A i 2 = {j ∈ A i : ϕ il (j) = 1}. Last, we define the following set , which represents all customers who leave the system. If j ∈ A i 1 , then there exist natural numbers k i j such that Last, j ∈ A i 3 means that the j th customer leaves the system after his first service and so s i (j) = 0. The quantities k i j and h i j denote the number of customers in node i = 1, 2 (including that one in service) who the j th customer meets after his departure from node l = 1, 2. Therefore, (4.3) can be written as Also, we define the number k i 0 = max j k i 0 (j), and the number of maximum routes for any external arrival in the system until time t, k 0 = max i=1,2 {k i 0 }. Observe that for any t ≥ 0, k 0 ∈ N and it is finite because we assume Markov routing. Take the following partition [0, t] = k0 m=0 [t m , t m+1 ], where t 0 = 0 and t k0+1 = t. We take the previous partition in such way, so that in each interval [t m , t m+1 ], any customer in the system can be routed only one time and so that in the interval [t k0 , t] there is no routing. Now, we can write the right-hand side of (4.1) as Applying the previous idea where customers are routed only one time per interval, we have that the above quantity can be written as Split again the last term of the previous quantity until time t 2 , and apply the previous idea when customers are routed only one time to obtain Adapting the previous steps until time t k0 , and recalling that the total arrival process is given by (2.7), we derive (4.2).
In the sequel, this result plays a key role. The next step is to define the so-called shifted fluid-scaled processes and to show that they are stochastically bounded.

Shifted fluid-scaled processes
We introduce the shifted fluid scaling, which is an extension of the classical fluid scaling. Let T > 0 and m ≤ nT . We defineĀ n,m i (t) = 1 n A n i (nm + nt) − A n i (nm) , and the analogous scaling for the processes T (·), Y L2 (·), and E(·). For the departure process, the cumulative service time process, and the routing process we have that Last, the queue length process is scaled as followsQ n,m i (t) = 1 n Q n i (nm + nt) and analogously for the scaling of the immediate and total workload. The system dynamics (2. In the sequel, we shall be referring to shifted fluid scaling, and shifted fluid process as shifted scaling, and shifted process for simplicity. The main step of SSC is to show the shifted process can be approximated by a solution of the fluid model. This is done in Section 4.4. We first need to prove that the shifted workload and shifted queue length are bounded at zero, which we do in the following section. Using these bounds and some properties of the cumulative service time (2.6), we can apply the results in [1, Sections 4 and 5].

Bounding the shifted processes
First, we find the relation between the diffusion scaling and the shifted scaling. Although this relation is easily obtained and is already known in the literature (e.g. [19]), we provide it here for completeness. We denote the event {ω ∈ Ω : max m≤nTW n,m T ot (0) ≤ B 1 } by G n 1 (B 1 ). Using now (4.10), it can be shown that the shifted queue length process is stochastically bounded at zero. (4.11) We denote the event {ω ∈ Ω : max m≤nTQ n,m i Proof. We prove the result by deriving a contradiction, so suppose that (4.11) does not hold. Thus, assume there exists at least one i such that Q n,m i (0) is not stochastically bounded. In other words, there exists a δ > 0 such that for any B > 0, Suppose that m n is such that it optimises the quantity max m≤nT Q n i (nm). We can choose in (4.10), = δ 3 and a (large enough) constant B 1 . Also, we choose a constant B such that B > 2 B1 βi . By the definition of the total workload (2.14) and (2.15), we have that We know that v (1) ii (j) are i.i.d. with mean β i . Also, in the previous summation j > S n i (T n i (nm n )), which means that v (1) ii (j) are independent of the process S n i (T n i (nm n )). Define the the following event G n = {ω ∈ Ω : 1 Bn S n i (T n i (nmn))+Bn By the weak law of large numbers (which we can apply due to the independence of the v (1) ii (j) and S n i (T n i (nm n ))), we have that for large n, P(G n ) ≥ 1 − δ 3 . In the sequel, we assume that ω ∈ G n ∩ G n 1 (B 1 ) ∩ (G n 2 (B)) c , and note that P(G n ∩ G n 1 (B 1 ) ∩ (G n 2 (B)) c ) ≥ δ 3 . Applying (4.12) and dividing by n we derive ii (j).
By (4.13) and the last inequality, we obtain for sufficiently large n, This yields a contradiction.
Having proved that the shifted processes are bounded, we can show that the shifted processes can be approximated by a solution of the fluid model. This is the topic of next section, in which we use a very similar approach as is Bramson [1, Sections 4 and 5].

Uniform fluid approximation
By [1, Proposition 5.1], we have that for any > 0, (4.14) Also, by [1,Proposition 5.2] it is known that the shifted arrival process is almost Lipschitz continuous, which means that for some N 1 > 0, Furthermore, using the definition of the cumulative service time (2.6), the property 2 i=1 R n i (q) = 1, and the observation that Y n L2 (·) and T n (·) are increasing functions in time, we conclude that the shifted process, T n,m (·), and the shifted idle time are Lipschitz continuous with constant equal to 1. Proof. It is shown in [1] that for a renewal process S(·), which is equivalent to (the process can start anywhere in the interval [0, Applying (4.15) to the last inequality we obtain (4.16).
In the following proposition, we show that all the shifted processes are almost Lipschitz continuous.

The scaled shifted total workload process
In this section, we see that we can approximate the scaled shifted total workload process by a solution to the fluid model. We begin with a preliminary result. Proof. The random variables s i (j) depend on n, but to keep the notation simple we omit the index n.
Note that (4.21) can be written as Adapting Remark 4.1, we can replace in the proposition above by (n) such that (n) → 0 as n → ∞. This will be done in the next result, where we combine all technical estimates so far to construct a "good" event.
Proposition 4.8. Let > 0 and G n 0 ⊆ Ω be as in Remark 4.1. Let G n 6 ⊆ Ω be the event such that the complement of (4.21) holds if we replace by (n). Define the event G n (B) = G n 0 (B) ∩ G n 6 . Then lim n→∞ P G n (B) ≥ 1 − .
Proof. Note that G n (B) = 6 i=1 G n i ∩ G n 1 (B) ∩ G n 2 (B). We denote (G i ) c the complement of the event G i . For 1 ≤ i ≤ 6, we have by construction that P((G n i ) c ) ≤ (n). Also, by (4.10) and (4.11) we can choose a constant B such that for i = 1, 2 lim n→∞ P G n i (B) ≥ 1 − 2 .
Combining the above inequalities, we obtain as n → ∞, In the sequel, we assume that ω ∈ G n (B). In other words, Proposition 4.8 allows us to use a samplepath approach. As a final step towards proving state-space collapse, we use this approach to show that there exists a fluid approximation for the total workload of the system. Proposition 4.9. There exists a solution of the fluid model equations, W T ot (·), such that for m ≤ nT W n,m T ot (·, ω) − W T ot (·) L ≤ (n) ∀ ω ∈ G n (B). Proof. Take > 0 and let ω ∈ G n (B) as defined in Proposition 4.8. Let the functions Q(t) and W (t) be such so that they satisfy (4.20). Define the function which is a solution to the fluid model equations because Q(t) is. Omitting again the index n in the quantity s i (j), by the definition of the total workload, we have that

State-space collapse
Now, we can state and prove SSC for the diffusion queue length process. as n → ∞.