2. The model and preliminaries

Let $X=\{(X_1(t),\ldots,X_n(t))\mid t\ge 0\}$ be a multivariate Lévy process with $X(0)=0$ , where $X_1$ is a Lévy process with no negative jumps which is not a subordinator and $X_i$ is a subordinator which is not identically zero for each $2\le i\le n$ . We have $\mathbb{E}[ {\mathrm{e}}^{-\alpha^T X(t)}] = {\mathrm{e}}^{\varphi(\alpha)t}$ , where the Laplace exponent $\varphi(\alpha)$ has the form

\begin{equation*}\varphi(\alpha)=-a^T\alpha+\frac{\sigma^2}{2}\alpha_1^2+\int_{\mathbb{R}_+^n}({\mathrm{e}}^{-\alpha^T x}-1+\alpha^T x1_{(0,1]^n}(x))\nu({\mathrm{d}} x) ,\end{equation*}

where, if we denote $\nu_i(A)=\nu(\mathbb{R}_+^{i-1}\times A\times\mathbb{R}_+^{n-i})$ for $1 \leq i \leq n$ , then

\begin{equation*}\int_{\mathbb{R}_+}(x^2_1\wedge 1)\nu_1({\mathrm{d}} x_1)<\infty\quad\text{and}\quad\int_{\mathbb{R}_+}(x_i\wedge 1)\nu_i({\mathrm{d}} x_i)<\infty\quad\text{for $2\le i\le n$.}\end{equation*}

For this set-up we can rewrite the Laplace exponent in the following way:

\begin{equation*}\varphi(\alpha)=-c^T\alpha+\frac{\sigma^2}{2}\alpha_1^2+\int_{\mathbb{R}_+^n}({\mathrm{e}}^{-\alpha^T x}-1+\alpha_1 x_1 1_{(0,1]}(x_1))\nu({\mathrm{d}} x) ,\end{equation*}

where

\begin{equation*}c_1=a_1+\int_{(0,1]\times(\mathbb{R}_+^{n-1}\setminus(0,1]^{n-1})}x_1\nu({\mathrm{d}} x_1) \quad\text{and}\quad c_i=a_i-\int_{(0,1]^n}x_i\nu({\mathrm{d}} x)\quad\text{for $2\le i\le n$.}\end{equation*}

Letting

\begin{align*}\varphi_1(\alpha_1)&=\varphi(\alpha_1,0,\ldots,0) \notag \\[3pt]&=-c_1\alpha_1+\frac{\sigma^2}{2}\alpha_1^2+\int_{\mathbb{R}_+}({\mathrm{e}}^{-\alpha_1x_1}-1+\alpha_1x_11_{(0,1]}(x_1))\nu_1 ({\mathrm{d}} x_1)\end{align*}

gives

(1) \begin{equation}\varphi(\alpha)=\varphi_1(\alpha_1)-\sum_{i=2}^nc_i\alpha_i-\int_{\mathbb{R}_+^n}\,{\mathrm{e}}^{-\alpha_1x_1}(1-{\mathrm{e}}^{-\sum_{i=2}^n\alpha_ix_i})\nu({\mathrm{d}} x) .\end{equation}

It should be noted that the integral on the right is finite for every choice of $\alpha\in \mathbb{R}_+^n$ . This follows from the fact that

\begin{equation*}{\mathrm{e}}^{-\alpha_1x_1}(1-{\mathrm{e}}^{-\sum_{i=2}^n\alpha_ix_i})\le \biggl(\sum_{i=2}^n\alpha_ix_i\biggr)\wedge 1 ,\end{equation*}

so that

\begin{align*}\int_{\mathbb{R}_+^n} \,{\mathrm{e}}^{-\alpha_1x_1}(1-{\mathrm{e}}^{-\sum_{i=2}^n\alpha_ix_i})\nu({\mathrm{d}} x) &\le\sum_{i=2}^n\alpha_i\int_{(0,1]^n}x_i\nu({\mathrm{d}} x)+\nu(\mathbb{R}_+^n\setminus(0,1]^n) \notag \\[3pt] &\le \sum_{i=2}^n \alpha_i \int_{(0,1]}x_i\nu_i({\mathrm{d}} x)+\nu(\mathbb{R}_+^n\setminus(0,1]^n) \notag \\[3pt] & <\infty .\end{align*}

We note that we may also write

\begin{equation*}\varphi(\alpha)=\varphi_1(\alpha_1)+\varphi(0,\alpha_2,\ldots,\alpha_n)+\int_{\mathbb{R}_+^n}(1-{\mathrm{e}}^{-\alpha_1x_1})(1-{\mathrm{e}}^{-\sum_{i=2}^n\alpha_ix_i})\nu({\mathrm{d}} x) ,\end{equation*}

observing that

(2) \begin{equation}\varphi(0,\alpha_2,\ldots,\alpha_n)=-\sum_{i=2}^nc_i\alpha_i-\int_{\mathbb{R}_+^n}(1-{\mathrm{e}}^{-\sum_{i=2}^n\alpha_i x_i})\nu({\mathrm{d}} x)\equiv-\eta(\alpha_2,\ldots,\alpha_n)\end{equation}

is the Laplace exponent of $(X_2,\ldots,X_n)$ , which implies (since it is an ( $n-1$ )-dimensional subordinator) that necessarily $c_i\ge 0$ for $2 \le i \le n$ .

In a similar manner, letting

\begin{equation*}\varphi_k(\beta)=\varphi(\underbrace{\beta,\ldots,\beta}_k,\underbrace{0,\ldots,0}_{n-k})\end{equation*}

for $\beta\in\mathbb{R}_+$ and $2\le k\le n$ , we have that $\varphi_k$ is the Laplace exponent of $\sum_{i=1}^kX_i$ , which (because of $X_1$ ) is not a subordinator, and, from (1),

(3) \begin{align}&\varphi(\underbrace{\beta,\ldots,\beta}_k,\alpha_{k+1},\ldots,\alpha_n)\notag \\[3pt]&\quad =\varphi_k(\beta)-\sum_{i=k+1}^nc_i\alpha_i-\int_{\mathbb{R}_+^n}\,{\mathrm{e}}^{-\beta\sum_{i=1}^kx_i}(1-{\mathrm{e}}^{-\sum_{i=k+1}^n\alpha_ix_i})\nu({\mathrm{d}} x) .\end{align}

3. The main results

In this section we determine the Laplace–Stieltjes transform of the steady-state buffer content vector for the system of n parallel queues (Theorem 1 and Corollary 1). The special structure of this buffer content transform allows us to rewrite it as a product of n joint Laplace–Stieltjes transforms. Each term of the product is given a probabilistic interpretation. We are thus able to interpret the buffer content vector as a sum of n independent random vectors. We end the section with a few remarks concerning the relations between the model under consideration and some other multivariate stochastic models.

Lemma 1. For $k=1,\ldots,n-1$ , there is a unique positive

\begin{equation*} \beta=\psi_k(\alpha_{k+1},\ldots,\alpha_n) \end{equation*}

with $\alpha_{k+1},\dots,\alpha_n \geq 0$ such that

(4) \begin{equation}\varphi(\beta,\ldots,\beta,\alpha_{k+1},\ldots,\alpha_n)=0 .\end{equation}

Proof. Since $X_i$ is not identically zero for each $2\le i\le n$ , then, when $\alpha_{k+1},\ldots,\alpha_n$ are not all zero, it follows from (3) and the fact that $c_2,\ldots,c_n \ge 0$ that

(5) \begin{equation} \varphi(0,\ldots,0,\alpha_{k+1},\ldots,\alpha_n)=-\sum_{i=k+1}^nc_i\alpha_i-\int_{\mathbb{R}_+^n} (1-{\mathrm{e}}^{-\sum_{i=k+1}^n\alpha_ix_i})\nu({\mathrm{d}} x)<0 . \end{equation}

Since $\sum_{i=1}^kX_i$ is not a subordinator, as $\beta \rightarrow \infty$ we have $\varphi_k(\beta)\to \infty$ and

\begin{equation*} \int_{\mathbb{R}_+^n}\,{\mathrm{e}}^{-\beta\sum_{i=1}^kx_i}(1-{\mathrm{e}}^{-\sum_{i=k+1}^n\alpha_ix_i})\nu({\mathrm{d}} x)\to 0 \end{equation*}

(dominated convergence), and thus (see (3)) $\varphi(\beta,\ldots,\beta,\alpha_{k+1},\ldots,\alpha_n)\to\infty$ . Therefore, since $\varphi$ is convex (hence $\varphi(\beta,\ldots,\beta,\alpha_{k+1},\ldots,\alpha_n)$ is convex in $\beta$ ), and $\varphi(0,\ldots,\break 0,\alpha_{k+1},\ldots,\alpha_n) <0$ according to (5), the statement of the lemma follows.□

Now assume that $Y_i=\sum_{j=1}^i X_j$ for $1\le i\le n$ . Then, for each $t\ge 0$ , we have for each $s<t$ that

\begin{equation*}Y_1(t)-Y_1(s)\le Y_2(t)-Y_2(s)\le\cdots\le Y_n(t)-Y_n(s) .\end{equation*}

Clearly $Y =\{(Y_1(t),\ldots,Y_n(t))\mid t \ge 0\}$ is also a Lévy process with Laplace exponent

\begin{equation*}\tilde\varphi(\alpha)=\varphi(\alpha_1+\cdots+\alpha_n,\alpha_2+\cdots+\alpha_n,\ldots,\alpha_n).\end{equation*}

Lemma 1 implies that for

(6) \begin{equation} \beta=\psi_k(\alpha_{k+1}+\cdots+\alpha_n,\ldots,\alpha_{n-1}+\alpha_n,\alpha_n)-\sum_{i=k+1}^n\alpha_i ,\end{equation}

$\tilde\varphi(0,\ldots,0,\beta,\alpha_{k+1},\ldots,\alpha_n)$ is zero.

We are now ready to study the buffer content of a system of n parallel fluid queues, with net input Y. Let

\begin{equation*}L_i(t)=-\inf_{0\le s\le t}Y_i(s), \quad Z_i(t)=Y_i(t)+L_i(t) .\end{equation*}

Note that $Z_i(t)$ can be viewed as the buffer content of a queue with net Lévy input $Y_i(t) = X_1(t)+\cdots + X_i(t)$ , $i=1,\ldots,n$ , and $L_i(t)$ is the local time at level 0 of that queue.

We necessarily have (see Theorem 6 of [Reference Kella and Whitt10]) that $Z_i(t)\le Z_{i+1}(t)$ for $1\le i\le n-1$ . Thus, if we assume that $\mathbb{E}Y_n(1)<0$ (hence $\mathbb{E}Y_i(1)<0$ for all $1\le i\le n$ ), then our system of n parallel fluid queues has a stationary distribution. We shall now determine the Laplace–Stieltjes transform (LST) of the steady-state workload vector, to be denoted by $Z^*$ .

Theorem 1. The LST of the steady-state workload vector $Z^*$ is given by

(7) \begin{equation}\tilde\varphi(\alpha)\mathbb{E}\,{\mathrm{e}}^{-\alpha^TZ^*}=\sum_{k=1}^{n-1}\alpha_k\,f_k(\alpha_{k+1},\ldots,\alpha_n)+\alpha_n\,f_n ,\end{equation}

where $f_n=-\mathbb{E}Y_n(1)$ is constant and $f_k(\alpha_{k+1},\ldots,\alpha_n)$ are recursively given by

(8) \begin{equation}f_k(\alpha_{k+1},\ldots,\alpha_n)=\frac{\sum_{i=k+1}^{n-1}\alpha_i\,f_i(\alpha_{i+1},\ldots,\alpha_n)+\alpha_n\,f_n}{\sum_{i=k+1}^n\alpha_i-\psi_k(\alpha_{k+1}+\cdots+\alpha_n,\ldots,\alpha_{n-1}+\alpha_n,\alpha_n)} ,\end{equation}

where an empty sum is defined to be zero.

Proof. If $Z^*$ has this distribution, then, since $Z_1^*\le\cdots\le Z_n^*$ (hence $Z_i^*=0$ implies that $Z_j^*=0$ for $1\le j\le i$ ), as in (2.12) of [Reference Kella6], (7) is satisfied for all $\alpha_1,\ldots,\alpha_n$ satisfying $\sum_{i=k}^n\alpha_i\ge 0$ for all $k=1,\ldots,n$ . Below we first determine $f_n$ , and thereafter we show how $f_{n-1},\ldots,f_1$ can successively be determined: $f_{n-1}$ is expressed in $f_n$ , then $f_{n-2}$ is expressed in $f_{n-1}$ and $f_n$ , etc. Letting $\alpha_1=\cdots=\alpha_{n-1}=0$ , we have, upon dividing by $\alpha_n$ and letting $\alpha_n\downarrow 0$ , that

\begin{equation*}f_n=\frac{\partial \tilde\varphi}{\partial \alpha_n}(0)=-\mathbb{E}Y_n(1)=-\sum_{i=1}^n \mathbb{E}X_i(1) .\end{equation*}

Now note that if we let $\alpha_1=\cdots=\alpha_{k-1}=0$ , then we have the formula

(9) \begin{equation}\tilde\varphi(0,\ldots,0,\alpha_k,\ldots,\alpha_n)\mathbb{E}\,{\mathrm{e}}^{-\sum_{i=k}^n\alpha_iZ^*_i}=\sum_{i=k}^{n-1}\alpha_i\kern0.5pt f_i(\alpha_{i+1},\ldots,\alpha_n)+\alpha_n\,f_n .\end{equation}

Assume that $f_i(\alpha_{i+1},\ldots,\alpha_n)$ are known for $i=k+1,\ldots,n-1$ . Set $\alpha_k$ to be the right-hand side of (6). Recall that $Z^*_{i-1}\le Z^*_i$ for $2\le i\le n$ . Therefore

\begin{equation*} \sum_{i=k}^n\alpha_i Z_i^*= {\psi_k(\alpha_{k+1}+\cdots+\alpha_n,\ldots,\alpha_{n-1}+\alpha_n,\alpha_n)Z^*_k} +\sum_{i=k+1}^n\alpha_i(Z_i^*-Z_k^*) \geq 0,\end{equation*}

so that $\mathbb{E}\,{\mathrm{e}}^{-\sum_{i=k}^n\alpha_iZ_i^*}\le 1$ for our choice of $\alpha_k$ (remember that we do not demand that all $\alpha_i \geq 0$ ). This implies (8).□

In principle we can thus compute the right side of (7) and hence $\mathbb{E}{\mathrm{e}}^{-\alpha^TZ^*}$ for all $\alpha_1,\ldots,\alpha_n$ satisfying $\sum_{i=k}^n\alpha_i\ge 0$ for all $k=1,\ldots,n$ . It follows from (2.12) of [Reference Kella6] that $f_k(\alpha_{k+1},\ldots,\alpha_n)$ is a constant times some (joint) Laplace–Stieltjes transform for every $1\le k\le n-1$ . We shall also see this in Corollary 1.

It should be mentioned that the above theorem generalizes Section 3 of [Reference Kella6] in two ways. The first is that it is not necessary to assume that $X_1$ is a subordinator minus a drift and the second is that it is not necessary to assume that $X_1,\ldots,X_n$ are independent. It also generalizes results from [Reference Badila, Boxma, Resing and Winands2] from the compound Poisson setting to the more general Lévy subordinator setting. The latter paper considers a system of n queues which simultaneously receive input from an n-dimensional compound Poisson process, and the jump sizes of the simultaneous jumps are stochastically ordered. The steady-state joint workload LST of that system is determined in [Reference Badila, Boxma, Resing and Winands2].

In [Reference Badila, Boxma, Resing and Winands2], furthermore, a decomposition is presented for the n-dimensional workload LST, and the terms of this decomposition are given an interpretation. In the remainder of this section, our aims are also to decompose the n-dimensional LST of $Z^*$ into a product of terms and to give interpretations of these terms. In particular, it would be nice to understand the meaning of $\varphi_k(\alpha_{k+1},\ldots,\alpha_n)$ . For this it would suffice to consider $k=1$ . We follow ideas from [Reference Kella6], where some unnecessary independence assumptions were made, but with a new observation which was missed at the time. The first thing to observe is the following. Let

\begin{equation*}T_1(x)=\inf\{t\mid X_1(t)=-x\}.\end{equation*}

Then it is well known that $T_1(x)$ is a.s. finite for each $x\ge 0$ and in fact $\{T_1(x)\mid x\ge 0\}$ is a subordinator with Laplace exponent $-\varphi^{-1}_1(\cdot)$ , where $\varphi^{-1}_1$ is the inverse of the strictly increasing and continuous function $\varphi_1$ (from $\mathbb{R}_+$ to $\mathbb{R}_+$ ).

Now, since ${\mathrm{e}}^{-\sum_{i=1}^n \alpha_i X_i(t)-\varphi(\alpha)t}$ is a mean-one martingale (Wald martingale) for each $\alpha\in \mathbb{R}^n_+$ , this holds in particular for $\alpha_1=\psi_1(\alpha_2,\ldots,\alpha_n)$ , in which case $\varphi(\alpha)=0$ . Since $T_1(x)$ is a stopping time, the optional stopping theorem implies that

\begin{equation*} \mathbb{E}\exp\biggl(-\psi_1(\alpha_2,\ldots,\alpha_n)X_1(T_1(x)\wedge t)-\sum_{i=2}^n\alpha_iX_i(T_1(x)\wedge t)\biggr)=1 .\end{equation*}

Noting that $\psi_1(\alpha_2,\ldots,\alpha_n)\ge 0$ and $X_1(T_1(x)\wedge t)\ge -x$ , and that $X_i(T_1(x)\wedge t)\ge 0$ for each $2\le i\le n$ , this implies, by bounded convergence and the fact that $X_1(T_1(x))=-x$ , that

(10) \begin{equation}\mathbb{E}{\mathrm{e}}^{-\sum_{i=2}^n\alpha_iX_i(T_1(x))}={\mathrm{e}}^{-\psi_1(\alpha_2,\ldots,\alpha_n)x} .\end{equation}

In [Reference Kella6] it was explained why $\{(X_2(T_1(x)),\ldots,X_n(T_1(x)))\mid x\ge 0\}$ is a (non-decreasing) Lévy process, and it was stated that it is important that $X_2,\ldots,X_n$ are independent. This statement about independence was an oversight, as with a similar argument it holds that it is a Lévy process without any independence assumptions. Formula (10) implies that $-\psi_1(\alpha_2,\ldots,\alpha_n)$ is in fact the Laplace exponent of this Lévy process. If $X_1$ is independent of $(X_2,\ldots,X_n)$ then, recalling (2),

\begin{equation*}\mathbb{E}\exp\biggl(-\sum_{i=2}^n\alpha_iX_i(T_1(x))\biggr)= \mathbb{E} \,{\mathrm{e}}^{-\eta(\alpha_2,\ldots,\alpha_n)T_1(x)}={\mathrm{e}}^{-\varphi_1^{-1}(\eta(\alpha_2,\ldots,\alpha_n))x} ,\end{equation*}

and thus for this case

\begin{equation*}\psi_1(\alpha_2,\ldots,\alpha_n)=\varphi^{-1}_1(\eta(\alpha_2,\ldots,\alpha_n)) ,\end{equation*}

which is what appears in [Reference Kella6] (with different notation). In fact, since $X_1(T_1(x))=-x$ , we clearly have that $\{X(T_1(x))\mid x\ge 0\}$ is an n-dimensional Lévy process with Laplace exponent $\alpha_1-\psi_1(\alpha_2,\ldots,\alpha_n)$ . In particular, this implies that the Laplace exponent of $\{Y(T_1(x))\mid x\ge0\}$ is

\begin{equation*}\sum_{i=1}^n\alpha_i-\psi_1(\alpha_2 + \cdots + \alpha_n,\ldots,\alpha_{n-1}+\alpha_n,\alpha_n) .\end{equation*}

At this point it would be good to refer to (6).

Now we note that, for $2\le i\le n$ ,

\begin{equation*}\inf_{0\le s\le T_1(x)}Y_i(s)=\inf_{0\le y\le x}Y_i(T_1(y))\end{equation*}

and thus

\begin{equation*}Z_i(T_1(x))=Y_i(T_1(x))-\inf_{0\le y\le x}Y_i(T_1(y)) .\end{equation*}

In addition to the formal proof given in [Reference Kella6], note that if $0\le s\le T_1(x)$ and

\begin{equation*}s\not\in \{T_1(y)\mid 0\le y\le x\}\end{equation*}

then, since $T_1(\cdot)$ is right continuous and $T_1(0)=0$ , there must be some $y\in[0,x]$ for which $T_1(y)<s$ and $X_1(T_1(y))=-y<X_1(s)$ . Since $X_2,\ldots,X_n$ are non-decreasing, $Y_i(T_1(y))<Y_i(s)$ . This means that the only contenders to minimize $Y_i$ on the interval $[0,T_1(x)]$ are $\{T_1(y)\mid 0\le y\le x\}$ .

Similar ideas imply that, with $T_k(x)=\inf\{t\mid Y_k(t)=-x\}$ ,

\begin{equation*}\{(X_{k+1}(T_k(x)),\ldots,X_n(T_k(x)))\mid x\ge0\}\end{equation*}

is a subordinator with Laplace exponent $-\psi_k(\alpha_{k+1},\ldots,\alpha_n)$ , and for every $k+1\le i\le n$ it follows that

\begin{equation*}Z_i(T_k(x))=Y_i(T_k(x))-\inf_{0\le y\le x}Y_i(T_k(y)),\end{equation*}

and we also observe that

\begin{equation*}Y_i(T_k(x))=\sum_{j=k+1}^iX_j(T_k(x))-x\quad\text{for $k+1 \le i \le n$.}\end{equation*}

That is, for each $k+1 \le i \le n$ , $Y_i(T_k(x))$ is a subordinator minus a unit drift.

Letting $Z^{k*}$ denote a random vector having the limiting distribution of $Z(T_k(x))$ , noting that necessarily $Z_1^{k*}=\cdots=Z_k^{k*}=0$ , then from Corollary 2.3 of [Reference Kella6] we have that

(11) \begin{equation}f_k(\alpha_{k+1},\ldots,\alpha_n)=-\mathbb{E}Y_k(1)\, \mathbb{E}\exp\biggl(-\sum_{i=k+1}^n\alpha_iZ^{k*}_i\biggr) .\end{equation}

The following corollary is now implied by Theorem 1.

Corollary 1.

\begin{equation*}\tilde\varphi(\alpha)\mathbb{E}\,{\mathrm{e}}^{-\alpha^TZ^*}=\sum_{k=1}^{n}\alpha_k (-\mathbb{E}Y_k(1))\, \mathbb{E}\exp\biggl(-\sum_{i=k+1}^n\alpha_iZ^{k*}_i\biggr) ,\end{equation*}

with an empty sum being zero.

This confirms the statement, made in the proof of Theorem 1, that all $f_k$ are LSTs of some (joint) distribution, up to a multiplicative constant.

We emphasize that $Z^{k*}$ is a vector of workloads having the steady-state distribution associated with the vector of workloads process embedded at time instants where station k (hence also stations $1,\ldots,k-1$ ) is empty. For what follows, we recall that if $-\xi(\beta)$ is the Laplace exponent of some one-dimensional subordinator, then for some $b\in\mathbb{R}_+$ and Lévy measure $\mu$ satisfying $\int_{\mathbb{R}_+}(u\wedge 1)\mu({\mathrm{d}} u)<\infty$ we have

(12) \begin{equation}\xi(\beta)=b\beta+\int_{\mathbb{R}_+}(1-{\mathrm{e}}^{-\beta u})\mu({\mathrm{d}} u) ,\end{equation}

and when the mean $\xi'(0)$ is finite, then ${{\xi(\beta)}/{(\beta\xi'(0))}}$ is the LST of the mixture, with weight factors

(13) \begin{equation}\biggl(\frac{b}{b+\int_0^\infty\mu(u,\infty)\,{\mathrm{d}} u} ,\,\frac{\int_0^\infty\mu(u,\infty)\,{\mathrm{d}} u}{b+\int_0^\infty\mu(u,\infty)\,{\mathrm{d}} u}\biggr)\end{equation}

of zero and a (residual lifetime) distribution having the density $\mu(x,\infty)/\int_0^\infty\mu(u,\infty)\,{\mathrm{d}} u$ . See e.g. (4.6) of [Reference Kella7].

We first discuss the workload decomposition for the case $n=2$ , exposing the key ideas; thereafter we briefly treat the case of a general n via a repetition of the argument.

3.1. Workload decomposition for the two-dimensional case

For the case $n=2$ , we first note that from $\varphi(\psi_1(\alpha_2),\alpha_2)=0$ it follows that with

\begin{equation*}\eta_2(\alpha_2)=-\varphi(0,\alpha_2)=c_2\alpha_2+\int_{\mathbb{R}_+}(1-{\mathrm{e}}^{-\alpha_2 x})\nu_2({\mathrm{d}} x) ,\end{equation*}

we have

\begin{equation*}\psi^{\prime}_1(0)=\frac{\eta^{\prime}_2(0)}{\varphi^{\prime}_1(0)} .\end{equation*}

It may be seen after some very simple manipulations that, for $\alpha_2\ge 0$ and $\alpha_1\ge -\alpha_2$ , with $\alpha_1\not=\psi(\alpha_2)$ ,

\begin{equation*} \mathbb{E}\,{\mathrm{e}}^{-\alpha_1Z_1^*-\alpha_2Z_2^*}=\frac{\varphi^{\prime}_1(0)(\alpha_1+\alpha_2-\psi_1(\alpha_2))}{\varphi(\alpha_1+\alpha_2,\alpha_2)}\cdot\frac{(1-\psi^{\prime}_1(0))\alpha_2}{\alpha_2-\psi_1(\alpha_2)} .\end{equation*}

An identical formula is given, employing different methods, in Proposition 1 of [Reference Badila, Boxma, Resing and Winands2]. That paper restricts itself to the special case where, with J a two-dimensional compound Poisson process with non-negative (two-dimensional) jumps, we have $X(t) = J(t)-(t,0)$ (so that $Y(t) = (J_1(t)-t,J_1(t)+J_2(t)-t)$ ).

If we denote $W_1^*=Z_1^*$ and $W_2^*=Z_2^*-Z_1^*$ (non-negative), then

(14) \begin{equation}\mathbb{E}\,{\mathrm{e}}^{-\alpha_1W_1^*-\alpha_2W_2^*}=\frac{\varphi^{\prime}_1(0)(\alpha_1-\psi_1(\alpha_2))}{\varphi(\alpha_1,\alpha_2)}\cdot\frac{(1-\psi^{\prime}_1(0))\alpha_2}{\alpha_2-\psi_1(\alpha_2)} .\end{equation}

The second expression in the product on the right-hand side is the LST of $W_2^{1*} = Z_2^{1*}$ (which has the steady-state distribution of the workload in station 2 observed only at times when station 1 is empty) and from its form (generalized Pollaczek–Khinchin formula, e.g. (2.5) of [Reference Kella7], among many others), it is indeed the steady-state LST of a reflected process of a Lévy fluid queue with subordinator input having the Laplace exponent $-\psi_1(\alpha_2)$ and a processing rate of one. That is, it is the steady-state LST of the reflected process associated with the driving process $\{\kern2pt\hat{\kern-2pt J}(x)-x\mid x\ge 0\}$ , where $\kern2pt\hat{\kern-2pt J}(x)=X_2(T_1(x))$ . In the compound Poisson setting of [Reference Badila, Boxma, Resing and Winands2], it was observed to be the steady-state LST of the workload in an M/G/1 queue with service times distributed as the extra (compared with queue 1) workload accumulated in station 2 during a busy period of station 1.

Setting $\alpha_1=0$ , we can rewrite the resulting equation as follows:

(15) \begin{equation}\frac{\eta_2(\alpha_2)}{\alpha_2\eta^{\prime}_2(0)}\mathbb{E}\,{\mathrm{e}}^{-\alpha_2W^*_2}=\frac{\psi_1(\alpha_2)}{\alpha_2\psi^{\prime}_1(0)}\mathbb{E}\,{\mathrm{e}}^{-\alpha_2 W_2^{1*}} .\end{equation}

That is, if $\xi_2^*$ and $\xi_2^{1*}$ have LSTs

\begin{equation*}\frac{\eta_2(\alpha)}{\alpha_2\eta^{\prime}_2(0)}\quad\text{and}\quad\frac{\psi_1(\alpha_2)}{\alpha_2\psi^{\prime}_1(0)}\end{equation*}

(see the paragraph that includes (12) and (13)) and are respectively independent of $W_2^*$ and $W_2^{1*}$ , then we have the following decomposition:

(16) \begin{equation}\xi^*_2+W^*_2\sim \xi_2^{1*}+W_2^{1*} .\end{equation}

We note that this is different from the decomposition described in Theorem 4.2 of [Reference Kella6], which holds when $X_1$ and $X_2$ are independent processes but not otherwise. It is easy to check that with this independence, this decomposition holds here as well, where, unlike there, $X_1$ need not be a subordinator minus some drift.

Our next goal is to identify a joint distribution with LST given by the first expression of the product on the right-hand side of (14). First we observe from Corollary (2.1) of [Reference Kella6] and the facts that $Z_1(T_1(x))=0$ , $Z_1(t)=0$ for points of (right) increase of $L_1$ and $Z_1(t)=Z_2(t)=0$ for points of (right) increase of $L_2$ (e.g. [Reference Kella8]), that

(17) \begin{align} &\tilde \varphi(\alpha)\mathbb{E}\int_0^{T_1(x)}\,{\mathrm{e}}^{-\alpha^TZ(s)}\,{\mathrm{d}} s \notag \\ &\quad =\mathbb{E}\,{\mathrm{e}}^{-\alpha_2 Z_2(T_1(x))}-1+\alpha_2\mathbb{E}L_2(T_1(x)) +\alpha_1 \mathbb{E}\int_0^{T_1(x)}\,{\mathrm{e}}^{-\alpha_2 Z_2(s)}\,{\mathrm{d}} L_1(s) .\end{align}

This holds in particular when we substitute $\alpha_1 = \psi_1(\alpha_2) - \alpha_2$ (see (6)), in which case $\tilde\varphi(\alpha)=0$ . Therefore, subtracting (17) with $\alpha_1=\psi_1(\alpha_2)-\alpha_2$ from (17) and noting (see also Corollary 2.1 of [Reference Kella6]) that

\begin{equation*}\int_0^{T_1(x)}\,{\mathrm{e}}^{-\alpha_2 Z_2(s)}\,{\mathrm{d}} L_1(s)=\int_0^x\,{\mathrm{e}}^{-\alpha_2 Z_2(T_1(y))}\,{\mathrm{d}} y,\end{equation*}

we have

\begin{equation*}\tilde \varphi(\alpha)\mathbb{E}\int_0^{T_1(x)}\,{\mathrm{e}}^{-\alpha^TZ(s)}\,{\mathrm{d}} s=(\alpha_1+\alpha_2-\psi_1(\alpha_2))\int_0^x\mathbb{E}\,{\mathrm{e}}^{-\alpha_2 Z_2(T_1(y))}\,{\mathrm{d}} y .\end{equation*}

Dividing by $\mathbb{E}T_1(x)=x/\varphi^{\prime}_1(0)$ and recalling that $\tilde\varphi(\alpha)=\varphi(\alpha_1+\alpha_2,\alpha_2)$ now gives

(18) \begin{equation}\frac{1}{\mathbb{E}T_1(x)}\mathbb{E}\int_0^{T_1(x)}\,{\mathrm{e}}^{-\alpha^TZ(s)}\,{\mathrm{d}} s=\frac{\varphi^{\prime}_1(0)(\alpha_1+\alpha_2-\psi_1(\alpha_2))}{\varphi(\alpha_1+\alpha_2,\alpha_2)}\frac{1}{x}\int_0^x\mathbb{E}\,{\mathrm{e}}^{-\alpha_2 Z_2(T_1(y))}\,{\mathrm{d}} y .\end{equation}

We now observe two facts. The first is that by regenerative theory the left-hand side is the LST of the steady-state distribution of a corresponding regenerative process where at time $T_1(x)$ the process is restarted. At this time $Z_1(T_1(x))=0$ and the remaining quantity at station 2 is lost. The second fact is that as $x\downarrow 0$ , by bounded convergence and the fact that $Z(T_1(\cdot))$ is right continuous with $Z(T_1(0))=0$ (since $T_1(0)=0$ ), we have

\begin{equation*}\lim_{x\downarrow 0}\frac{1}{x}\int_0^x \mathbb{E}\,{\mathrm{e}}^{-\alpha_2Z_2(T_1(y))}\,{\mathrm{d}} y=\lim_{x\downarrow 0}\mathbb{E}\,{\mathrm{e}}^{-\alpha_2 Z_2(T_1(x))}=1 .\end{equation*}

Using (18), this implies that

\begin{equation*}\lim_{x\downarrow 0}\frac{1}{\mathbb{E}T_1(x)}\mathbb{E}\int_0^{T_1(x)}\,{\mathrm{e}}^{-\alpha^TZ(s)}\,{\mathrm{d}} s=\frac{\varphi^{\prime}_1(0)(\alpha_1+\alpha_2-\psi_1(\alpha_2))}{\varphi(\alpha_1+\alpha_2,\alpha_2)} ,\end{equation*}

where the right-hand side converges to 1 as $(\alpha_1,\alpha_2)\to (0,0)$ . By the continuity theorem for LSTs we have that necessarily the right-hand side is an LST of a non-negative random vector (for any $\alpha_2\ge 0$ and $\alpha_1\ge -\alpha_2$ ). This also implies that for $\alpha_1,\alpha_2\ge 0$ , with $W(s) =(W_1(s),W_2(s)) := (Z_1(s),Z_2(s)-Z_1(s))$ ,

(19) \begin{equation}\lim_{x\downarrow 0}\frac{1}{\mathbb{E}T_1(x)}\mathbb{E}\int_0^{T_1(x)}\,{\mathrm{e}}^{-\alpha^TW(s)}\,{\mathrm{d}} s=\frac{\varphi^{\prime}_1(0)(\alpha_1-\psi_1(\alpha_2))}{\varphi(\alpha_1,\alpha_2)} ,\end{equation}

so this is also an LST. For the compound Poisson case, this implies that this is the LST of a version of this process such that, whenever $Z_1(t)=0$ , any quantity available in queue 2 is lost. For this case, this interpretation was discovered in [Reference Badila, Boxma, Resing and Winands2] and it continues to be valid if $X_1$ is a compound Poisson process with a negative drift and $X_2$ is a general subordinator ( $X_1,X_2$ can be dependent).

Remark 1. To obtain moments of $W_1^*,W_2^*$ or $Z_1^*=W_1^*$ , $Z_2^*=W_1^* + W_2^*$ requires slightly tedious but straightforward calculations. First (see (14)),

\begin{equation*}\mathbb{E}({\mathrm{e}}^{-\alpha_1 W_1^*}) = \frac{\varphi^{\prime}_1(0)\alpha_1}{\varphi(\alpha_1,0)}\end{equation*}

yields

\begin{equation*}\mathbb{E}W_1^* = \frac{\varphi^{\prime\prime}_1(0)}{2 \varphi^{\prime}_1(0)} = \frac{\text{Var} X_1(1)}{-2 \mathbb{E}X_1(1)} .\end{equation*}

Second (again from (14)),

\begin{equation*}\mathbb{E}( {\mathrm{e}}^{-\alpha_2 W_2^*}) =- \frac{\varphi^{\prime}_1(0)\psi_1(\alpha_2)}{\varphi(0,\alpha_2)}\cdot\frac{(1-\psi^{\prime}_1(0))\alpha_2}{\alpha_2-\psi_1(\alpha_2)}\end{equation*}

yields

\begin{equation*}\mathbb{E}W_2^* =\frac{\bigl(\tfrac{{\mathrm{d}} ^2}{{\mathrm{d}} y^2}\bigr) \varphi(0,y)}{2 \bigl(\tfrac{{\mathrm{d}} }{{\mathrm{d}} y}\bigr) \varphi(0,y)}\bigg|_{y=0}-\frac{\psi^{\prime}_1(0)}{2 \psi^{\prime}_1(0)}- \frac{\psi^{\prime}_1(0)}{2(1-\psi^{\prime}_1(0))} ,\end{equation*}

which term by term corresponds to (see (15) and (16))

\begin{equation*}\mathbb{E}W_2^* = - \mathbb{E} \xi_2^{*} + \mathbb{E}\xi_2^{1*} + \mathbb{E}W_2^{1*} .\end{equation*}

Insertion of

\begin{gather*}\frac{{\mathrm{d}} }{{\mathrm{d}} y} \varphi(0,y)|_{y=0} = -\mathbb{E}X_2(1), \quad \frac{{\mathrm{d}} ^2}{{\mathrm{d}} y^2} \varphi(0,y) = \text{Var}(X_2(1)),\\\psi^{\prime}_1(0) = - \frac{\mathbb{E}X_2(1)}{\mathbb{E}X_1(1)},\\\begin{aligned} \psi^{\prime}_1(0) &=-\frac{\varphi^{\prime}_1(0)}{\varphi^{\prime}_1(0)} (\psi^{\prime}_1(0))^2 -\frac{\bigl(\tfrac{{\mathrm{d}} ^2}{{\mathrm{d}} y^2}\bigr) \varphi(0,y)|_{y=0}}{\varphi^{\prime}_1(0)} - 2\biggl(\frac{{\mathrm{d}} ^2}{{\mathrm{d}} x\, {\mathrm{d}} y}\biggr) \varphi(x,y)|_{x=y=0} \frac{\psi^{\prime}_1(0)}{\varphi^{\prime}_1(0)}\notag\\&=\frac{\text{Var}(X_1(1))}{\mathbb{E}X_1(1)} \biggl(\frac{\mathbb{E}X_2(1)}{\mathbb{E}X_1(1)}\biggr)^2 + \frac{\text{Var}( X_2(1))}{\mathbb{E}X_1(1)}- 2 \text{Cov}(X_1(1),X_2(1)) \frac{\mathbb{E}X_2(1)}{(\mathbb{E}X_1(1))^2}\end{aligned}\end{gather*}

allows us to express the moments of $W_1^*$ and $W_2^*$ in terms of the first two moments and covariance of $X_1(1)$ and $X_2(1)$ . Compared to formula (4.6) of [Reference Kella6], the expression for $\mathbb{E}W_2^*$ contains an extra term that includes $\text{Cov}(X_1(1),X_2(1))$ .

3.2. Workload decomposition for the n-dimensional case

The above argument may be repeated for the n-dimensional case. Indeed, it follows from (9), (8), and (11) that for $\alpha$ such that $\sum_{i=k}^n\alpha_i \ge 0$ for all $k=1,\ldots,n$ ,

\begin{align*} &\mathbb{E}\,{\mathrm{e}}^{-\alpha^TZ^*}\\ &\quad =\frac{\varphi^{\prime}_1(0)(\sum_{i=1}^n\alpha_i-\psi_1(\alpha_2+\cdots+\alpha_n,\ldots,\alpha_{n-1}+\alpha_n,\alpha_n))}{\tilde\varphi(\alpha)}\ \mathbb{E}\exp\biggl(-\sum_{i=2}^n\alpha_iZ_i^{1*}\biggr) ,\end{align*}

or equivalently, introducing $\alpha_i := \alpha_i - \alpha_{i+1}$ for $1 \leq i \leq n-1$ , we have for $\alpha_i\ge 0$ , for $i=1,\ldots,n$ ,

\begin{align*}\mathbb{E}{\mathrm{e}}^{-\alpha^TW^*}&=\frac{\varphi^{\prime}_1(0)(\alpha_1-\psi_1(\alpha_2,\ldots,\alpha_{n-1},\alpha_n))}{\varphi(\alpha)}\ \mathbb{E}\exp \biggl(-\sum_{i=2}^n\alpha_i W_i^{1*}\biggr) .\end{align*}

By induction, the right-hand side may be written as a product of n (joint) LSTs where, for $1\le k\le n-1$ , the kth multiple is

\begin{equation*}\frac{\bigl(1-\bigl(\tfrac{\partial\psi_{k-1}}{\partial\alpha_k}\bigr)(0)\bigr)(\alpha_k-\psi_k(\alpha_{k+1},\ldots,\alpha_n))}{\alpha_k-\psi_{k-1}(\alpha_k,\ldots,\alpha_n)}\end{equation*}

and the last one is

\begin{equation*}\frac{(1-\psi^{\prime}_{n-1}(0))\alpha_n}{\alpha_n-\psi_{n-1}(\alpha_n)} .\end{equation*}

From (4) it follows that, for $2\le k\le n$ ,

\begin{equation*}1-\frac{\partial\psi_{k-1}}{\partial\alpha_k}(0)=\frac{\varphi^{\prime}_1(0)-\sum_{i=2}^k\eta^{\prime}_i(0)}{\varphi^{\prime}_1(0)-\sum_{i=2}^{k-1}\eta^{\prime}_i(0)}=\frac{-\mathbb{E}Y_k(1)}{-\mathbb{E}Y_{k-1}(1)} ,\end{equation*}

where for $k=2$ the empty sum in the denominator is defined to be zero.

For the case where $X_1$ is a subordinator minus a drift and $X_1,\ldots,X_n$ are independent, it seems that this decomposition can be inferred from more general results reported in Theorem 6.1 of [Reference Debiçki, Dieker and Rolski5].

Finally, the ideas that led to (19) may be repeated to conclude that each multiplicative factor participating in this decomposition is indeed a (joint) LST.

Remark 2. The results of this section are also of immediate relevance for a tandem fluid model, i.e. a model of n stations in series, in which material or work leaves each station as a fluid. Such a connection between parallel stations and stations in series was already pointed out in [Reference Kella6]. Tandem fluid models were introduced in [Reference Kella and Whitt9]. The following n-station tandem fluid model is introduced and studied in that paper. The input process of station 1 is a non-decreasing Lévy process and the jth station receives the output of the $(j-1)$ th station at a constant rate $r_{j-1}$ , as long as that station is not empty. It was assumed that $r_1 \ge \cdots \ge r_n$ , to avoid the trivial case that a station is always empty. That tandem fluid model was generalized in [Reference Kella6] by allowing additional external inputs to stations $2,\dots,n$ . Those inputs were assumed to be subordinators and, for most of the results, they were assumed to be independent from each other and from the input to station 1. It was observed in [Reference Kella6] that there is a direct connection between the workloads in this tandem fluid model and the workloads in a model of n parallel stations. The same observation, but for the case of compound Poisson inputs (and allowing dependence between the input processes), was also made in [Reference Badila, Boxma, Resing and Winands2]. That connection also extends to the case of dependent external Lévy inputs. More precisely, let $X_1,\dots,X_n$ be the external inputs to stations $1,\dots,n$ of the tandem fluid model, with $X_2,\dots,X_n$ being subordinators, and let $W_1,\dots,W_n$ denote their buffer content level processes. Then we can identify $W_j(t)$ with $ Z_j(t)-Z_{j-1}(t)$ , $1 \le j \le n$ (hence $W_1(t)+\cdots+W_j(t) = Z_j(t)$ ), where $Z_j(t)$ as before denotes the buffer content level of station j in the system of parallel stations studied in the present paper, and $Z_0(t) \equiv 0$ .

Remark 3. In [Reference Kella6] an equivalence between the tandem fluid model and a particular single server priority queue is also pointed out. Assume a compound Poisson input vector X of customer classes $1,\ldots,n$ . Class i has preemptive resume priority over class j if $i<j$ ; the total workload in the first j classes can now be identified with $Z_j^*$ .

In [Reference Badila, Boxma, Resing and Winands2] a multivariate duality relation is established between (i) the model of n parallel M/G/1 queues with simultaneous arrivals, with stochastic ordering of the n service times of each arriving batch, and (ii) a Cramér–Lundberg insurance risk model featuring n insurance companies with simultaneous claim arrivals and stochastic ordering of those claims. In particular, when the arrival processes in both systems have the same distribution, the joint steady-state workload distribution $\mathbb{P}(V_1 \le x_1,\dots,V_n \le x_n)$ in the queueing model equals the survival probability for all companies, with initial capital vector $(x_1,\dots,x_n)$ . This is a generalization of a well-known duality that is discussed, for example, on page 46 of [Reference Asmussen and Albrecher1]. Using the sample path argument presented there, one should also be able to generalize this duality to the case of the Lévy input process of the present paper, thus also obtaining the survival probability for n insurance companies with the same n-dimensional input process as the n parallel stations of the present paper.