Two coupled Levy queues with independent input

We consider a pair of coupled queues driven by independent spectrally-positive Levy processes. With respect to the bi-variate workload process this framework includes both the coupled processor model and the two-server fluid network with independent Levy inputs. We identify the joint transform of the stationary workload distribution in terms of Wiener-Hopf factors corresponding to two auxiliary Levy processes with explicit Laplace exponents. We reinterpret and extend the ideas of Cohen and Boxma (1983) to provide a general and uniform result with a neat transform expression.

We consider a pair of coupled queues driven by independent spectrally-positive Lévy processes. With respect to the bi-variate workload process this framework includes both the coupled processor model and the two-server fluid network with independent Lévy inputs. We identify the joint transform of the stationary workload distribution in terms of Wiener-Hopf factors corresponding to two auxiliary Lévy processes with explicit Laplace exponents. We reinterpret and extend the ideas of Cohen and Boxma (1983) to provide a general and uniform result with a neat transform expression.
1. Introduction. In the queueing literature, several studies have been devoted to a queueing model of two servers, each with their own customer arrival process, with the special feature that the speed of one server changes when the other server becomes idle. This has become known as the coupled processor model. A possibly even more popular model of two servers is a fluid network with independent arrival processes, where fixed fractions of fluid exiting one queue are routed into the same and the other queue, as well as out of the system. These models are intimately related and in the case of Lévy input both can be put in our framework below; see Section 1.1 and Section 1.3 for possible interpretations of our model and for related literature respectively.
More specifically, we assume that our queues are driven by two independent Lévy processes X 1 (t) and X 2 (t) without negative jumps. We model a pair of workload processes (W 1 (t), W 2 (t)) as a 2-dimensional reflected process, see e.g. Harrison and Reiman (1981); Kella (2006), where W i (t) are nonnegative, L i (t) are nonnegative and nondecreasing with L i (0) = 0, and, in addition, it is required that if t is a point of increase of L i (t) then W i (t) = 0. Sometimes the latter condition is replaced by an equivalent integral condition or minimality requirement. We assume that r 1 , r 2 ≥ 0 and r 1 r 2 < 1, in which case workload processes (W 1 (t), W 2 (t)) (with given initial values) and unused capacity processes (L 1 (t), L 2 (t)) satisfying the above conditions exist and are unique, see (Kella, 2006, Sec. 5).
1.1. Interpretations. It is the easiest to understand the model given by (1) in the case of compound Poisson inputs and constant service rates c i > 0, i.e. when each X i (t) is a compound Poisson process (CPP) minus c i t. Note that when W i (t) hits zero it stays at zero until the arrival of the next customer, which leads to the following four cases. While W 1 (t), W 2 (t) > 0 these workload processes evolve according to X 1 (t) and X 2 (t). Whereas while W 1 (t) = 0 and W 2 (t) > 0 the process L 1 (t) evolves as c 1 t resulting in an additional service rate r 2 c 1 in the second queue, i.e. the fraction r 2 of the first server capacity is used to help the second. Similarly, while W 1 (t) > 0, W 2 (t) = 0 the service rate in the first queue is c 1 + r 1 c 2 . Finally, when both queues are empty, the processes L i (t) evolve as certain linear drifts canceling the negative drifts of X i (t) and each other's influence, which is possible if and only if r 1 r 2 < 1. It is noted that compound Poisson input allows for a formulation of the coupled processor model, which goes beyond our assumption of r 1 r 2 < 1. One simply replaces (1) by an explicit description of the workload processes in the above four cases, see Cohen and Boxma (1983).
As mentioned above, our model includes two-dimensional fluid networks with independent Lévy input, where the column vector of workloads is a reflected process of the form: see e.g. Kella (1996). HereX(t) is a column vector of external non-decreasing input processes into each queue, P is a routing matrix (a substochastic matrix) with P n → 0 for n → ∞, P ′ is its transpose, I is the identity matrix, and c is a column vector of (maximal) service rates. One usually interprets ct −L(t) as a vector of cumulative outflows from the queues, which are routed according to P . We can write 1−p jj we obtain (1) and guarantee the above conditions. We remark that commonlyX i (t) is a non-decreasing Lévy process and hence X i (t) is a Lévy process without negative jumps having bounded variation (on finite intervals). We allow X i (t) to be general Lévy processes without negative jumps, which may lead to a certain debate about an appropriate model for the fluid network, because cumulative outflows (if defined at all) are not necessarily non-decreasing in this general setup. Nevertheless, such models have appeared in the literature, see e.g. Kella and Ramasubramanian (2012). Finally, one can go the other way around and produce a network from the model (1), which is immediate if r 1 , r 2 ≤ 1. If r 1 > 1 (and similarly for r 2 > 1) then consider (W 1 (t), r 1 W 2 (t)) and note that it corresponds to a pair of workload processes in a network with routing matrix given by p 11 = p 22 = 0, p 12 = r 1 r 2 , p 21 = 1 and driving processes X 1 (t), r 1 X 2 (t), see also (Kella, 1996, Lem. 4.1).
1.2. Stationary distribution. Let us note that (I − P ′ ) −1 EX(1) < 0, where X(t) is a multidimensional driving process, is a sufficient stability condition for a general Lévy network, i.e. under this condition the joint workload process has a limiting distribution (independent of initial conditions), which is also a unique stationary distribution, see (Kella and Ramasubramanian, 2012, Thm. 2.4). Furthermore, if none of X i (t) is a zero process then this condition is also necessary. Stronger limiting results are available in Kella (1996) for the case when both X i (t) have bounded variation. Stability of (1) can be easily related to the stability of the corresponding network yielding the following condition where d i = −EX i (1). Assuming that (2) holds we let a pair of random variables (W 1 , W 2 ) refer to the joint stationary distribution of (W 1 (t), W 2 (t)). Our main result is an expression for the transform Ee −α 1 W 1 −α 2 W 2 in terms of Wiener-Hopf factors corresponding to two auxiliary processes with explicit Laplace exponents, see Theorem 1. We reinterpret and extend the ideas from Chapter III.3 of Cohen and Boxma (1983) to provide a general and uniform result. Its derivation is rather compact, and is based on a number of identities and observations from the fluctuation theory for Lévy processes.
Let us shortly discuss a special case, when X 1 (t) is a subordinator, i.e. a non-decreasing Lévy process. In this case L 1 (t) can increase only when L 2 (t) increases, hence both queues should be empty. This feature allows for a rather simple analysis of the joint transform similarly to Kella and Whitt (1992a). So we can assume in the following that each X i (t) is a spectrallypositive Lévy process, i.e. it is a Lévy process which is not a subordinator, and which can have only positive jumps.
1.3. Related literature and motivation. The main application/motivation of the coupled processor model is provided by the fact that, in a network of work stations, a user may use other machines than its own when those machines are idle; this is often referred to as cycle stealing. Another application occurs in integrated-service communication networks. Differentiated quality-of-service among different traffic flows is achieved in such networks via scheduling algorithms such as Weighted Fair Queueing. Mathematically, such scheduling algorithms may often be represented by a form of Generalized Processor Sharing, where traffic flow i gets a weight factor w i ∈ (0, 1), with w i = 1. If all traffic flows are backlogged, then flow i is served at rate w i . If some of the flows are not backlogged, then the excess capacity is redistributed among the backlogged flows proportionally to their weights. Again, this may be viewed as a form of cycle stealing. A pioneering paper on the mathematical analysis of coupled processors is Fayolle and Iasnogorodski (1979). They consider two M/M/1 queues with service speeds c 1 and c 2 , respectively, unless the other queue is empty; then the speeds are c * 1 and c * 2 , respectively. They study the two-dimensional queue length process, and show how the generating function of the joint steady-state queue length distribution can be obtained via the solution of a Riemann-Hilbert boundary value problem. Konheim, Meilijson and Melkman (1981) provide an elegant solution of the special, slightly easier, case of two symmetric M/M/1 queues in which a server doubles its speed when the other server is idle (one might say that the idle server helps the other one). Cohen and Boxma (1983) have generalized the model of Fayolle and Iasnogorodski (1979) to the case of general service time distributions. They consider the two-dimensional workload process. See Cohen (1992) for a further extension to the case that arrivals may also simultaneously occur at both queues.
The analytic approach of Fayolle and Iasnogorodski (1979); Cohen and Boxma (1983), exploiting a relation to boundary value problems in complex function theory, seems to be limited to two dimensions. This has led to work in the following directions. (i) Application of a numerical-analytic method, the Power Series Algorithm, gives numerical results for more than two coupled processors Blanc (1987); Hooghiemstra, Keane and van de Ree (1988). (ii) Osogami, Harchol-Balter and Scheller-Wolf (2003) have developed an approximation method which yields mean response times; the approximation can be made as accurate as desired. (iii) Several studies (see, e.g., Jelenkovic (2000, 2003)) consider tail asymptotics of workloads for coupled processors and multi-queue systems with some form of Generalized Processor Sharing.
The body of literature concerning fluid networks with Lévy input is huge. The joint transform of the stationary workload in such networks is not known apart from a few special cases. The transform can be obtained for tandem and feed-forward networks with decreasing service rates (in the direction of flow), see e.g. Kella and Whitt (1992a) and Debicki, Dieker and Rolski (2007). In such networks, if one queue is empty then all the queues preceding it are empty as well. This is the main feature facilitating the computation, which can be also guaranteed in some other models, see Badila et al. (2012). Otherwise, the only tractable examples concern networks of two queues, which are closely related to the coupled processor model discussed above. The main result of the present paper yields an exact expression for the transform in a two-node network with independent Lévy input. This model generalizes, e.g., a tandem queue of Miyazawa and Rolski (2009), for which the authors study tail asymptotics.
1.4. Organization of the paper. Section 2 summarizes basic facts about spectrally-positive Lévy processes and about Wiener-Hopf factorization of Lévy processes. In Section 3 we relate a spectrally-positive Lévy process to a certain pure-jump subordinator which plays a fundamental role in our main result, Theorem 1, formulated in Section 4. Section 5 contains the proof of Theorem 1. It basically consists of three steps. In Subsection 5.1 we derive a functional equation for the joint workload transform; in Subsection 5.2 the kernel of that functional equation is studied, and the functional equation is solved via Wiener-Hopf factorization assuming certain bounds; these bounds are established in Subsection 5.3. Some special cases are considered in Section 6, where we also discuss the result of Cohen and Boxma (1983) in the case of CPP inputs.
2. Basic facts. For ease of reference let us recall the Lévy-Khintchine formula for a spectrally-positive Lévy process X(t) (cf. Kyprianou (2006)): The process X(t) has bounded variation (on finite intervals) if and only if σ = 0 and 1 0 xν(dx) < ∞, in which case we have an alternative representation and µ can be interpreted as a linear drift. The case of µ = 0 corresponds to a pure-jump subordinator. This subordinator is either a CPP or an infinite activity subordinator according to ν(0, ∞) being finite or infinite.
Differentiating under the integral sign in (3), which can be justified, we get for α > 0. This shows that X(t) has bounded variation if and only if lim α→∞ φ ′ (α) is finite.
Finally, let us recall the celebrated Wiener-Hopf factorization for a general (two-sided) Lévy process X(t) and some positive constant p > 0, see e.g. (Kyprianou, 2006, Thm. 6.16 and comments on p. 167 about the CPP case). Letting ℜ(α) be the real part of α ∈ C, consider the Laplace transforms where e p is an independent exponentially distributed r.v. with rate p and X(t), X(t) denote supremum and infimum processes respectively. Note that Ψ ± (α) are analytic in the corresponding half-planes and continuous on the imaginary axis. They satisfy the following factorization for w ∈ iR: Let us finally note that identification of the Wiener-Hopf factors is a difficult but well-studied problem with some numerical evaluation techniques available, see e.g. Den Iseger, Gruntjes and Mandjes (2013).
3. Fundamental subordinators. Consider a spectrally-positive Lévy process X(t), which will serve as a driving process in our model. The goal of this section is to associate to X(t) a certain pure-jump subordinator Y (t), which will play a fundamental role in our main result. Recall that φ(α) denotes the Laplace exponent of X(t), and d = −EX(1) = φ ′ (0) ∈ (−∞, ∞). Note that we have excluded only d = −∞, which is allowed because of the stability condition (2).
It is known that −α/Φ(α) is the Laplace exponent of a certain exponentially killed subordinator (ascending ladder time process, see e.g. (Kyprianou, 2006, p. 170)). Note that if d ≥ 0 then lim α↓0 If, however, d < 0 then this limit is 0. Hence d + = d ∨ 0 is the rate of killing, which we remove to obtain a non-killed subordinator Y (t) with the Laplace exponent This is a pure-jump subordinator, which follows from lim α→∞ φ Y (α)/α = 0 and representation (4). Note also that (8) holds for all α = 0 with ℜ(α) ≥ 0 by analyticity and continuity of Laplace exponents, see Section 5.2. Let us consider the dichotomy of bounded and unbounded variation for the process X(t): This can be seen by considering , which is finite in the first case and is −∞ in the second as was discussed in Section 2. So in the first case we have P(Y (1) = 0) > 0 and in the second P(Y (1) = 0) = 0, which correspond to a CPP and an infinite activity subordinator respectively.
4. Transform of the stationary workload. Consider the model specified by (1), where r i ≥ 0 and r 1 r 2 < 1. Recall that X 1 (t) and X 2 (t) are two independent spectrally-positive Lévy processes with Laplace exponents φ i (α), d i = −EX i (1) ∈ (−∞, ∞), and assume that stability condition (2) holds. Let Y i (t) be a pure-jump subordinator associated to X i (t), whose Laplace exponent φ Y i (α) is given in (8). Define two Lévy processes and two positive constants: Their corresponding Laplace exponents for w ∈ iR, w = 0 are given by Finally, we let Ψ ± L (α) be the Wiener-Hopf factors corresponding to X L (t) and rate parameter p L . Similarly, Ψ ± R (α) are the Wiener-Hopf factors corresponding to X R (t) and p R .
Theorem 1. The joint transform of the stationary workloads is given by It is noted that if d 1 , d 2 ≥ 0 then p 0 L = p L and p 0 R = p R . Moreover, if r i = 0 then d i > 0 according to (2) implying d − i = 0. In this case 0/0 in the definition of p 0 R , p 0 L is interpreted as 0. Consider the above systems of queues for r 1 = r 2 = 0. Then X R (t) = Y 1 (t) and X R (t) = 0 for all t, and p R = d 1 . From the definition of the W-H factors we have Ψ − R (α) = 1 and Ψ + . Plugging in α = φ 1 (α 1 ) we obtain Ψ + R (φ 1 (α 1 )) = d 1 α 1 φ 1 (α 1 ) , and similarly we obtain expressions for the other terms. Putting things together we see that (11) becomes which is indeed the transform of the workload in two independent queues. Another verification of Theorem 1 is given in Section 6.1, where we assume that X 2 (t) = −d 2 t for d 2 > 0. For such a (degenerate) system we first provide a quick alternative derivation of the joint transform and then check it against Theorem 1.

Proof.
In this section we prove Theorem 1. The proof consists of three steps. In Subsection 5.1 we derive a functional equation for the joint workload transform; in Subsection 5.2 the kernel of that functional equation is studied, and the functional equation is solved via Wiener-Hopf factorization assuming certain bounds; these bounds are established in Subsection 5.3.
Remark 1. Our proof of Theorem 1 is mainly analytic. The result is however formulated in terms of Wiener-Hopf factors corresponding to processes constructed from Y i (t), the ascending ladder time processes. This hints that there may be a direct probabilistic proof based on fluctuation theory of Lévy processes, which may provide better insight into the problem and help in solving other important problems concerning coupled queues and risk processes. This type of proof is given in Debicki, Dieker and Rolski (2007) for a feed-forward network, and remains an open challenge for a general network as considered in the present paper.
5.1. The functional equation. In this section we derive an equation for the two-dimensional joint workload transform, which involves two unknown functions. Identification of these functions is the main problem, which will be addressed in Subsections 5.2 and 5.3. The following result is based on a by now standard argument using the Kella-Whitt martingale, see Kella and Whitt (1992b).
Remark 2. The functional equation (13) can be used to derive the following identity for the means: One way is to put α 1 = r 2 α, α 2 = α, express F 2 (r 2 α), differentiate it at α = 0, and then to do the same for α 1 = α, α 2 = 0. Then the above identity follows by expressing F ′ 2 (0) from these equations. Note also that if d 1 , d 2 > 0 then the right side of (15) is r 2 d 1 EV 1 + r 1 d 2 EV 2 , where V i refers to the stationary workload in queue i considered alone. It may be an interesting exercise to prove this relation probabilistically from first principles, at least for Poisson inputs.
We use a similar uniformization approach as in Cohen and Boxma (1983) for the special case of compound Poisson input: we consider (13) for α 1 = Φ 1 (w) and α 2 = Φ 2 (−w), where w ∈ iR lies on the imaginary axis, and obtain Assuming w = 0 we multiply this equation by which immediately translates into according to (10). Finally, from the Wiener-Hopf factorization (6) we have which also holds for w = 0 by continuity. The left side of (17) is analytic in the left-half plane and the right side is analytic in the right-half plane, and both are continuous and coincide on the boundary. So one is an analytic continuation of the other, see (Lang, 1999, Thm IX.1.1). Assume for a moment that the so-obtained entire function is bounded; this will be proven in Subsection 5.3. Then by Liouville's theorem (Lang, 1999, Thm III.7.5) it is a constant, call it C.
Let us determine the constant C, by plugging w = 0 in (17). According to the stability condition (2) at least one of d i is positive. If d 1 > 0 then Φ 1 (0) = 0 and hence C = p R E * L 2 (1), whereas C = p L E * L 1 (1) if d 2 > 0. Note also that for a stationary system which yields and provides the expression for C. Furthermore, where p 0 L , p 0 R are given in (12). This can be checked by considering three scenarios d i ≥ 0, d 2 < 0, d 1 < 0 separately.
This together with (13) completes the proof of Theorem 1 under the assumption that the entire function defined by (17) is a constant.
5.3. Bounds on the entire function. In this section we show that the entire function defined by (17) is a constant. Consider (14) and observe that F 1 (α) and F 2 (α) are bounded for ℜ(α) ≥ 0. LetX L (t) andX R (t) be the processes X L (t) and X R (t), see (9), in the model defined by (1) but with interchanged indices. That is, we consider the same system with reversed indexing. Observe that Therefore, it is sufficient to analyze Ψ + L (w)/Ψ + R (w), ℜ(w) ≥ 0. Recall the following Spitzer-type identity for a Lévy process X L (t): see Thm. 6.16 and comments on p. 168 in Kyprianou (2006). Observe also that for ℜ(w) ≥ 0 we have (1 − e −wx )P(X L (t) ∈ dx) dt. 5.3.1. Bounded variation case. If both X 1 (t) and X 2 (t) have bounded variation then X L (t) and X R (t) are CPPs and hence This immediately shows that A(w) and hence Ψ + L (w)/Ψ + R (w) are bounded.

5.3.2.
The general case. The proof in the general case is based on the following proposition.
Proposition 2 together with (18) implies that the entire function defined by (17) is bounded by a polynomial and hence it is a polynomial itself, see (Lang, 1999, Cor. III.7.4). Taking limit along the reals shows that this polynomial is just a constant. The proof of Proposition 2 relies on the following technical lemma.
For positive z the bound can be replaced by 1 + ln(|z| ∨ 1).
Proof. For |z| > 1 we write The result is immediate for |z| < 1 and for positive z.
Proof of Proposition 2. For positive w the function A(w) defined in (20) is bounded from above by Recall that φ Y 1 (w) ≤ 0, and use Lemma 1, to bound A(w) by 1 + ln(|φ Y 1 (w)| ∨ 1). Finally, exponentiate and use (8) to establish that Ψ + L (w)/Ψ + R (w) = o(w). For w with |w| > 1, ℜ(w) ≥ 0 we mimic the steps in (21) to show that A similar bound (with the constant 2r 2 EY 1 (1)) can be obtained for the term involving X L (t). Assuming that EY 1 (1) < ∞ we have |A(w)| ≤ C + 4 ln |w| for |w| > 1 and some constant C, which yields the result. If EY 1 (1) = ∞ then we can easily reduce our problem to the one, where the large jumps of Y 1 (t) are removed and hence EY 1 (1) < ∞.
6.2. Queues fed by Brownian motion. Suppose φ i (α) = d i α + α 2 /2, i.e. the driving processes are standard Brownian motions with drifts. Then Φ i (α) = −d i + d 2 i + 2α and so which corresponds to an inverse Gaussian subordinator. Hence the processes X L and X R can be seen as differences of two inverse Gaussian processes. Their respective Laplace exponents are given by The final step according to Theorem 1 is to identify the Wiener-Hopf factors corresponding to these Laplace exponents.
6.3. Compound Poisson input. This subsection briefly examines the relation between the general result given in Theorem 1 and the result of Cohen and Boxma (1983) for CPP inputs. Assume that customers arrive into queue i with intensity λ i and bring iid amount of work distributed as B i , and the server speed is s i . Suppressing the index i, the Laplace exponent of the driving process X(t) is where ρ = λEB and R has the stationary residual life distribution associated to B. This further leads to α where R is assumed to be independent of the driving process X(t) and hence of its first passage time τ − x . Note that τ − R has the interpretation of the length of the busy period in a queue driven by X(t) and started with workload R. For simplicity we assume that s − ρ > 0 and hence τ − R is a proper positive random variable, which we denote by U .

Consider (23)
− w Φ 2 (−w) + r 2 w Φ 1 (w) = s 2 − ρ 2 Ee wU 2 + r 2 s 1 − r 2 ρ 1 E −wU 1 appearing in (16). Note that this expression can be rewritten as (s 2 + r 2 s 1 )(1 − zEe −wŨ ), where z = (ρ 2 + r 2 ρ 1 )/(s 2 + r 2 s 1 ) ∈ (0, 1) andŨ is a mixture of U 1 and −U 2 . This allows to apply Wiener-Hopf factorization for the random walk induced byŨ to decompose (23) into a product of functions, which are analytic in different half-planes, see Cohen and Boxma (1983) for details. Finally, we mention that a similar technique can be used, when both X i (t) are Lévy processes of bounded variation (a subordinator minus linear drift). In this case R can be interpreted as an asymptotic overshoot of the corresponding subordinator. This technique fails to generalize further. In general one can use Wiener-Hopf factorization for Lévy processes as is demonstrated in the present paper, which furthermore provides a uniform and neat solution, see Theorem 1.