# The $X(6900)$ peak could be a molecular state

###### Abstract

The analyses of the LHCb data on $X(6900)$ found in di-$J/\psi$ and
$J/\psi\psi(3686)$ systems are performed using a momentum-dependent Flatté-
like parameterization. The use of the pole counting rule and spectral density
function sum rule give consistent conclusions that $X(6900)$ may not be a
molecule of $J/\psi\psi(3686)$. Nevertheless it is still possible that
$X(6900)$ be a molecule of higher states, such as $J/\psi\psi(3770)$,
$\chi_{c0}\chi_{c2}$, etc.

Ye Lu${}^{\dagger\,}$, Chang Chen${}^{\dagger\dagger\,}$, Kai-Ge
Kang${}^{\dagger\dagger\,}$, Guang-you Qin${}^{\dagger\,}$, Han-Qing
Zheng${}^{\heartsuit\,}$
† Institute of Particle Physics and Key Laboratory of Quark and Lepton Physics
(MOE), Central China Normal University, Wuhan, Hubei 430079, China
†† Department of Physics and State Key Laboratory of Nuclear Physics and
Technology, Peking University, Beijing 100871, China
♡ College of Physics, Sichuan University, Chengdu, Sichuan 610065, China

The LHCb Collaboration has observed a structure named as $X(6900)$ in the
di-$J/\psi$ invariant mass spectrum [1], with the signal statistical
significance above 5$\sigma$. It is probably composed of four (anti)charm
quarks ($c\bar{c}c\bar{c}$) and its width [1] are determined to be $80\pm
19~{}({\rm stat.})\pm 33~{}({\rm sys.})$ and $168\pm 33~{}({\rm stat.})\pm
69~{}({\rm sys.})$ MeV in two fitting scenarios of Breit-Wigner
parameterizations with constant widths. Additionally, a broad bump and a
narrow bump exist in the low and high sides of the di-$J/\psi$ mass [1],
respectively, where the former might be a result from a lower broad resonant
state (or several lower states) or interference effect, and the latter is
found to be a hint of a state located at $\sim$7200 MeV, named as $X(7200)$.
The $X(6900)$ peak in di-$J/\psi$ channels is also found by the CMS
Collaboration [2]. More recently the $X(6900)$ peak is reported to be observed
in the $J/\psi\psi(3686)$ invariant mass spectrum [3].

The experimental observation has triggered tremendous studies, see for example
Refs. [4] –[30]. 111For an incomplete list of earlier studies, see for example
references listed in Ref. [4].  Generally speaking, a molecular state may
locate near the threshold of two color singlet hadrons, like deuteron,
$Z_{b}(10610)$ [31], $Z_{c}(3900)$ [32, 33, 34], $P_{c}(4470)$ [35, 36],
$etc.$ The $X(6900)$ state is close to the threshold of $J/\psi\psi(3770)$,
$J/\psi\psi_{2}(3823)$, ${J/\psi\psi_{3}(3842)}$, and $\chi_{c0}\chi_{c1}$;
and the $X(7200)$ is close to the threshold of $J/\psi\psi(4160)$ and
$\chi_{c0}$$\chi_{c1}(3872)$. Inspired by this, it is studied in this paper,
as an extension to the work of Ref. [4], the properties of $X(6900)$ and
$X(7200)$, by assuming for example the $X(6900)$ coupling to $J/\psi J/\psi$,
$J/\psi\psi(2s)$, $J/\psi\psi(3770)$, $J/\psi\psi_{2}(3823)$,
${J/\psi\psi_{3}(3842)}$ and $\chi_{c0}\chi_{c2}$ channels, etc.; and
$X(7200)$ to $J/\psi J/\psi$, $J/\psi\psi(2s)$ and $J/\psi\psi(4160)$
channels, etc. Notice that in this paper we limit ourselves discussing only
the $S$-wave ($l=0$) couplings, since only in this situation we are able to
distinguish a molecular state from an ‘elementary’ or a confining state. Such
assignment corresponds to $J^{P}=0^{++}$, $2^{++}$ spin quantum numbers of
$X(6900)$. For the $S$-wave $J/\psi J/\psi$ coupling, the pole counting rule
(PCR) [37], which has been applied to the studies of “$XYZ$” physics in Refs.
[38, 39, 41, 40], and spectral density function sum rule (SDFSR) [41, 42, 43,
44, 45] are employed to analyze the nature of the two structures in Ref. [4].
It is found that the di-$J/\psi$ data alone is not enough to judge the
intrinsic properties of the two states. It is also pointed out that the
$X(6900)$ is unlikely a molecule of $J/\psi\psi(3686)$ [4] – a conclusion
drawn before the new data from Ref. [3].

In this paper we make an upgraded analysis on $X(6900)$ by using the new
$J/\psi\psi(3686)$ data. What we conclude from this reanalysis is that, even
though the $X(6900)$ is very unlikely a molecule of $J/\psi\psi(3686)$, it is
not necessarily an “elementary state”, i.e., a compact $\bar{c}\bar{c}cc$
tetraquark state. Considering that the tetra-quark idea is rather attractive
in the literature (see for example, refs. [11], [12], [18], [19], [30] ), it
is important to make a more careful re-analysis on the nature of $X(6900)$.
The present analysis points out that, there still exists the possibility that
the $X(6900)$ state be a molecular state composed of particles which
thresholds are closer to $X(6900)$ comparing with $J/\psi\psi(3686)$, such as
$J/\psi\psi(3770)$, $J/\psi\psi_{2}(3823)$, ${J/\psi\psi_{3}(3842)}$ and
$\chi_{c0}\chi_{c2}$, etc.

We start with a two channel re-analysis (i.e., di-$J/\psi$,
$J/\psi\psi(3686)$), using the neural network program which has been recently
developed in the study of ‘XYZ’ states in Ref. [25]. Then we make a further
numerical three channel studies by including another nearby threshold and
point out that a strong coupling to the third channel is not excluded by the
current data, and hence the $X(6900)$ may still be of molecule nature.

Examination of the couple channel situation by using neural network

For the purpose, a supervised learning scheme is adopted as in Ref. [25]. It
means that we should get some ‘molecule’ type and ‘elementary’ type samples
beforehand. Then these samples are taken into the machine learning program.
The samples named ‘molecule’ are with label ‘0’ and named ‘elementary’ are
with label ‘1’. Our purpose is to let program find out the different
characteristics of these input line shapes between ‘0’ and ‘1’. We use PCR
[37] as our criteria. That is, if there exists one pole near one threshold in
$s$ plane, the S-wave resonance can be interpreted as the molecule in that
channel. On the other side, if there exists one more nearby pole (on different
sheets), the $S$-wave resonance is labeled as ‘elementary’. This method has
been used in many works to study the nature of exotic hadron states, see for
example, Refs. [41, 40, 38, 16].

To simulate a couple channel resonance amplitude, we chose the Flatté-like
parametrization to generate training data. For a two channel parametrization
of invariant mass spectrum, it is written as

$\displaystyle\mathcal{R}_{i}(s)$
$\displaystyle=\rho_{i}(s)\left|\frac{e^{i\phi_{i}}}{s-M^{2}+iM\left(g_{1}\rho_{1}(s)+g_{2}\rho_{2}(s)\right)}+bg_{1i}\right|^{2}$
(1) $\displaystyle+\rho_{i}(s)bg_{2i},$

where $\rho_{i}$ ($i=1,2$) represents two-body phase space factor for final
state 1 (FS1) and final state 2 (FS2), respectively. Parameter $M$ is the bare
mass of the resonance. Coherent and incoherent background contributions are
also considered as the noise which should be weak in order not to influence
the judgement on lineshapes. In the region we care about, the background
terms, $bg_{1i}$ and $bg_{2i}$ are all first order polynomials of $\sqrt{s}$.
All the parameters are tuned for generating a resonance near FS2 threshold and
there is one or two poles near FS2 representing ‘0’ and ‘1’ samples,
respectively. In practice, the coupling constants $g_{1},g_{2}$ are important
to lineshapes as well as the pole positions in the complex s plane [25].

Considering that the experimental data is affected by the energy resolution,
to match the real signal better, we add a gaussian convolution to the model.
According to the experiment by LHCb [1], we fix $\Delta=5~{}\mathrm{MeV}$ and

$\frac{d\sigma_{i}}{d\sqrt{s}}\sim\frac{1}{\sqrt{2\pi}\Delta}\int_{\sqrt{s}-3\Delta}^{\sqrt{s}+3\Delta}\mathcal{R}_{i}\left(s^{\prime}\right)e^{-\frac{\left(\sqrt{s}-\sqrt{s^{\prime}}\right)^{2}}{2\Delta^{2}}}d\sqrt{s^{\prime}}\
.$ (2)

According to above equation, when the parameter space is determined, all the
training data can be generated. Some pre-treatment is made to fit our machine
learning program. The first is to set the window size of the signal. For FS1
di-$J/\psi$ signal, we set the window size from $E_{th2}-100$ MeV to
$E_{th2}+300$ MeV to cover the relevant data and keep away from the noise from
the peak X(6200), where $E_{th2}$ is the threshold energy for FS2. For
$J/\psi\psi(2S)$ FS2 signal, we set the window size from $E_{th2}$ to
$E_{th2}+200$ MeV. See Figs. 1 and 2 for two typical examples of the generated
training data. The energy interval is fixed at 1 MeV uniformly so it can meet
the demand that all the inputs sent to a special neural network should have
the same size. The values at every energy point are calculated using a linear
extrapolation of the experimental data, and their errors are taken randomly as
5, 10, 15 percent of their values. Finally, the normalization is made before
sending them as learning input. In Fig. 3 examples of pole locations of two
‘molecules’ and an ‘elementary’ state are drawn and are labeled appropriately
using PCR.

Figure 1: An example of ‘0’ signal from training data for 2 channel case.
Figure 2: An example of ‘1’ signal from training data for 2 channel case.
Figure 3: The pole locations of ‘molecule 1’ (blue dots), ‘molecule 2’ (green
dots), and ‘elementary state’ (yellow dots) of the training data. Here
‘molecule 1’ corresponds to a bound state of channel 2, whereas ‘molecule 2’
corresponds to a virtual pole of channel 2.

Before sending the experimental data into the neural network, it is noted that
the energy interval is wider than $1$ MeV. We have to manage these data to
make the input the same size as the neural network. To be specific, if the
energy interval of experimental data is larger than 1 MeV then the method of
linear interpolation will be employed to supply extra points.

Fig. 4 shows the loss function of our two channel machine for both training
data and test data. One can see that the data are trained well and can be used
to judge whether a peak stands for a molecule of a given threshold. As shown
in Fig. 5, the 100 data points got from the experimental data are all judged
around ‘1’ by the trained machine. It tells us that $X(6900)$ is not a
molecule composed of $J/\psi\psi(3686)$.

Figure 4: left: loss function of the training set (blue), and test set
(orange); right: the output result of test set.
Figure 5: Test result of true data, it indicates that $X(6900)$ is not a
molecule of $J/\psi\psi(2S)$.

A triple channel study:

From above discussion, it is concluded that $X(6900)$ is not a molecule of
$J/\psi\psi(2S)$, but it does not necessarily mean that $X(6900)$ is
‘elementary’ since it may still be a molecule in some other higher channels,
e.g., $J/\psi\psi(3770)$, $\chi_{c0}\chi_{c2}$, etc. To examine this
possibility, we extend the Flatté-like parametrization to three channels:

$\displaystyle\mathrm{Event}_{i}(s)$
$\displaystyle=\mathrm{C}_{i}\rho_{i}(s)\left|\frac{e^{i\phi_{i}}}{s-M^{2}+iM\left(g_{1}\rho_{1}(s)+g_{2}\rho_{2}(s)+g_{3}\rho_{3}(s)\right)}+{\left(a+b(\sqrt{s}-E_{thi})\right)}\right|^{2}.$
(3)

In this part, we try to fit $J/\psi J/\psi$ and $J/\psi\psi(2S)$ signal
directly by using above Flatté-like parametrization equation. $E_{thi}$,
$i=1,2$ are the threshold energies of FS1 and FS2. All the parameters $g_{1}$,
$g_{2}$, $g_{3}$, $\phi_{i}$, $a$, $b$, $\mathrm{C}_{i}$ are considered as fit
parameters where $\mathrm{C}_{i}$ is a normalization constant. Here the ‘third
channel’ is taken as $\chi_{c0}\chi_{c2}$ for example (corresponding to
$2^{++}$ assignment of $X(6900)$).222The $J/\psi\psi(3770)$ or
$\chi_{c1}\chi_{c1}$ molecule ($0^{++}$ or $2^{++}$) are also possible, though
slightly less preferred. In the former situation, $X(6900$ is a virtual state
composed of $J/\psi\psi(3770)$. As expected, the current data leave room for
totally different scenarios. See table 1, apparently there exist multi-
solutions. One is an ‘elementary’ solution, which is similar to the couple
channel solution discussed previously. However, when the coupling to the
‘third’ channel, $g_{3}$, tuned large, a molecule (of the third channel)
solution appears. See figs. 6, 7 and figs 8, 9 for illustration.

Parameter | $g_{1}$ | $g_{2}$ | $g_{3}$ | M | a | b | $\phi_{1}$ | $\phi_{2}$ | $\mathrm{C}_{1}$ | $\mathrm{C}_{2}$ |  $\chi^{2}$/d.o.f
---|---|---|---|---|---|---|---|---|---|---|---
Elementary | 0.64 | 0.07 | 0.21 | 6.93 | 0.88 | -0.39 | 1.49 | 3.82 | 276 | 275 | 0.69
Molecule | 0.69 | 0.48 | 1 | 6.91 | 0.82 | -0.46 | 1.31 | 3.78 | 477 | 263 | 0.73
Table 1: Fit parameters for two different solutions.
Figure 6: Fit to the invariant mass spectrum, for the ‘elementary solution’ in
table 1. The data are from Ref. [2].
Figure 7: Pole positions of the ‘elementary solution’ in table 1.
Figure 8: Fit to the invariant mass spectrum, for the molecule solution in
table 1. The data are from Ref. [2].
Figure 9: Pole positions of the molecule solution in table 1.

Based on the standard PCR analysis we suggest that it is possible the
$X(6900)$ peak may be a molecule composed of, for example,
$\chi_{c0}\chi_{c_{2}}$. In principle, one may also make use of neural network
to test the three channel situation. However in the absence of data of the
third channel, it is very difficult to make a reliable test. Finally we would
like to comment on the work of Ref. [46], in which the authors also adopted a
triple channel study (di-$J/\psi$, $J/\psi\psi(2S)$, $J/\psi\psi(3770)$). They
conclude from their fit results (without including the new $J/\psi\psi(2S)$
data of Ref. [3]) that the $X(6900)$ is likely a true tetra quark state.
However it is clear indicated from our analysis that the fit program contains
redundancy of fit parameters and it certainly leaves the room for a ‘molecule’
solution. So additional experimental information in the
$\chi_{c0}\chi_{c_{2}}$, $J/\psi\psi(3770)$ channels, etc., are needed to
further clarify the issue on whether $X(6900)$ is a molecular state or a
compact tetra-quark state.

Acknowledgements: This work is supported in part by National Nature Science
Foundations of China under Contract Number 11975028 and 10925522.

## References

  * [1] Roel Aaij et al. (LHCb Collaboration), Sci.Bull. 65 (2020) 23, 1983.
  * [2] CMS collaboration, CMS-PAS-BPH-21-003 (2022).
  * [3] ATLAS collaboration, ATLAS-CONF-2022-040 (2022).
  * [4] Q. F. Cao $et$ $al.$, Chin. Phys. C 45 (2021) 10, 103102.
  * [5] C. Gong, M. C. Du, Qiang Zhao, X. H. Zhong, and B. Zhou, Phys. Lett. B 824 (2022) 136794.
  * [6] X. Dong, F. K. Guo and B. S. Zou, Phys. Rev. Lett. 126 (2021) 152001.
  * [7] M. Z. Liu and L. S. Geng, Eur. Phys. J. C 81 (2021) 2, 179.
  * [8] V. P. Gonçalves and B. Moreira, Phys. Lett. B 816 (2021) 136249.
  * [9] H. W. Ke, X. Han, X. H. Liu and Y. L. Shi, Eur. Phys. J. C 81 (2021) 5, 427.
  * [10] Z. R. Liang, X. Y. Wu and D. L. Yao, Phys. Rev. D 104 (2021) 3, 034034.
  * [11] H. Mutuk, Eur. Phys. J. C 81 (2021) 4, 367.
  * [12] G. J. Wang, L. Meng, M. Oka, and S. L. Zhu, Phys. Rev. D 104 (2021) 3, 036016.
  * [13] Q. N. Wang, Z. Y. Yang, W. Chen, and H. X. Chen, Phys. Rev. D 104 (2021) 1, 014020.
  * [14] A. V. Nefediev, Eur. Phys. J. C 81 (2021) 8, 692.
  * [15] Q. F. Lü, D. Y. Chen, Y. B. Dong and E. Santopinto, Phys. Rev. D 104 (2021) 5, 054026.
  * [16] H. Chen, H. R. Qi and H. Q. Zheng, Eur. Phys. J. C 81 (2021) 9, 103102.
  * [17] Q. N. Wang, Z. Y. Yang and W. Chen, Phys. Rev. D 104 (2021) 11, 114037.
  * [18] A. Esposito, C. Andrea Manzari, A. Pilloni and A. D. Polosa, Phys. Rev. D 104 (2021) 11, 114029.
  * [19] F. X. Liu, M. S. Liu, X. H. Zhong and Q. Zhao, Phys. Rev. D 104 (2021) 11, 116029.
  * [20] Z. J. Zhuang, Y. Zhang, Y. Z. Ma and Q. Wang, Phys. Rev. D 105 (2022) 5, 054026.
  * [21] S. Z. Chen et al., arXiv: 2111.14360 [hep-ph].
  * [22] R. F. Lebed, Moscow Univ.Phys. Bull. 77 (2022) 2, 458.
  * [23] Z. G. Wang, Int. J. Mod. Phys. A 37 (2022) 31n32, 2250189.
  * [24] R. H. Wu et al., JHEP 11 (2022) 023.
  * [25] C. Chen, H. Chen, W. Q. Niu, H. Q. Zheng, Eur. Phys. J. C 83 (2023)1, 52.
  * [26] Z. R. Liang and De-Liang Yao, Rev. Mex. Fis. Suppl. 3 (2022) 3, 0308042.
  * [27] C. Gong, M. C. Du and Q. Zhao, Phys. Rev. D 106 (2022) 5, 054011.
  * [28] H. X. Chen, Y. X. Yan and W. Chen, Phys. Rev. D 106 (2022) 9, 094019.
  * [29] G. J. Wang, Q. Meng and M. Oka, Phys. Rev. D 106 (2022) 9, 096005.
  * [30] H. Mutuk, Phys. Lett. B 834 (2022) 137404.
  * [31] A. Bondar et al. (Belle Collaboration), Phys. Rev. Lett. 108 (2012) 122001.
  * [32] M. Ablikim et al. (BESIII Collaboration), Phys. Rev. Lett. 110 (2013) 252001.
  * [33] Z. Q. Liu et al. (Belle Collaboration), Phys. Rev. Lett. 110 (2013) 252002.
  * [34] T. Xiao, S. Dobbs, A. Tomaradze and K. K. Seth, Phys. Lett. B 727 (2013) 366.
  * [35] R. Aaij et al. (LHCb Collaboration), Phys. Rev. Lett. 115 (2015) 072001.
  * [36] R. Aaij et al. (LHCb Collaboration), Phys. Rev. Lett. 122 (2019) 222001.
  * [37] D. Morgan, Nucl. Phys. A 543 (1992) 632.
  * [38] O. Zhang, C. Meng and H. Q. Zheng, Phys. Lett. B 680 (2009) 453.
  * [39] L. Y. Dai, M. Shi, G. Y. Tang and H. Q. Zheng, Phys. Rev. D 92 (2015) 1, 014020.
  * [40] Q. F. Cao, H. R. Qi, Y. F. Wang and H. Q. Zheng, Phys. Rev. D 100 (2019) 5, 054040.
  * [41] Q. R. Gong, Z. H. Guo, C. Meng, G. Y. Tang, Y. F. Wang and H. Q. Zheng, Phys. Rev. D 94 (2016) 11, 114019.
  * [42] V. Baru, J. Haidenbauer, C. Hanhart, Y. Kalashnikova and A. E. Kudryavtsev, Phys. Lett. B 586 (2004) 53.
  * [43] S. Weinberg, Phys. Rev. 130 (1963), 776.
  * [44] S. Weinberg, Phys. Rev. 137 (1965), B672.
  * [45] Y. S. Kalashnikova and A. V. Nefediev, Phys. Rev. D 80 (2009) 074004.
  * [46] Q. Zhou, D. Guo, S. Q. Kuang, Q. H. Yang and L. Y. Dai, Phys. Rev. D 106 (2022) 11, L111502.


# On the Persistent-Idle Load Distribution Policy Under Batch Arrivals and
Random Service Capacity

Rami Atar(2), Isaac Keslassy(2), Gal Mendelson(1), Ariel Orda(2), Shay
Vargaftik(3)
1\. Graduate School of Business, Stanford University
2\. Viterbi Faculty of Electrical Engineering, Technion
3\. VMware Research

###### Abstract

The Persistent-Idle (PI) load distribution policy was recently introduced as
an appealing alternative to current low-communication load balancing
techniques. In PI, servers only update the dispatcher when they become idle,
and the dispatcher always sends jobs to the last server that reported being
idle. PI is unique in that it does not seek to push the server queue lengths
towards equalization greedily. Rather, it aggressively pulls the servers away
from starvation. As a result, PI’s analysis requires different tools than
other load balancing approaches. So far, PI was proven to achieve the
stability region for Bernoulli arrivals and deterministic and constant service
capacities.

Our main contribution is proving that PI achieves the stability region in a
model with batch arrivals and random service capacities. Proving this result
requires developing tighter bounds on quantities of interest and proving the
fulfillment of a generalized version of the state-dependent drift criteria
previously considered. We also present PI-Split, which allows batches to be
split among idle servers and prove its stability. Finally, we conduct
simulations comparing PI and PI-Split to other load balancing policies.
Simulation results indicate that both achieve excellent performance for the
model under consideration for a wide range of parameters.

Keywords— Parallel Server Model; Load Balancing; Persistent-Idle; Stability;
Lyapunov Drift

##  1 Introduction

The Persistent Idle (PI) load distribution policy was recently introduced in
[1] as an appealing reduced-state load balancing scheme. It is designed for
the parallel server model, in which jobs arrive at a dispatcher and are routed
to one of $n$, possibly heterogeneous, servers. In PI, a job is routed to an
idle server if there is one or to the last server a job was sent to,
otherwise. Thus the dispatcher only uses idleness information to make its
routing decisions making the communication cost at most a single message per
job.

Despite this low communication overhead, PI was shown to achieve the stability
region when servers are heterogeneous; namely, with possibly different
processing rates. In their considered model the following assumption were
made: at most one job may arrive at a time slot; each job’s workload depends
on the server it is assigned to and equals the number of time slots required
to process it on that server; each server completes one unit of work at each
time slot. Thus at most one departure can occur in each server at each time
slot.

In this paper, we study PI in a considerably different, well-established model
with batch arrivals and random service capacities (e.g. [2]). An additional
variant of the policy is considered, called PI-Split, in which the dispatcher
may split jobs belonging to the same batch to several servers. Our main
results are that PI and PI-Split achieve the stability region in this model.

Proving that PI achieves the stability region in this model presents a
significant challenge. In particular, the proof in [1] relies on sampling the
Markov chain describing the state of the system at the random times at which
new servers become idle, and proving a state-dependent drift criteria is
fulfilled by a carefully chosen Lyapunov function. This proof heavily relies
on the fact that, given the current state, the next sampling time is
deterministic. This is no longer the case in the model we consider.

Our main contribution is proving the fulfillment of a generalized version of a
state-dependent drift criteria (we discuss this at length in Section 3.2.1)
and proving it is sufficient for stability. Our proof includes developing
tighter bounds on quantities of interest and a method to relate the current
queue length state to the time interval’s duration until the next sampling
time.

Finally, we present simulation results comparing PI and PI-Split to other
load-distribution policies that were considered in [1]. Simulation results
indicate that PI has a comparable performance to Join-the-Shortest-Queue (JSQ)
up to high loads, and favorable performance compared to other well-known
reduced-state load balancing schemes, i.e., Join-the-Idle-Queue (JIQ) [3],
Power-of-choice (JSQ(2)) [4] and Power-of-Memory (JSQ(1,1)) [5].

The remainder of this paper is organized as follows. Section 2 contains the
system model and the definition of the PI and PI-Split policies. Section 3
contains the stability results for PI and PI-Split. Section 4 contains the
simulation results, comparing PI and PI-Split to other load-distribution
policies.

##  2 Model

###  2.1 Persistent Idle

We consider a parallel-server system that evolves in discrete time
$t\in\\{0,1,2\ldots\\}$ with a single dispatcher and $n$ heterogeneous work-
conserving servers. Each of the $n$ servers has an unbounded FIFO queue
holding pending jobs.

Order of events. The order of events during each time slot is as follows. (1)
New jobs may arrive at the dispatcher; (2) the dispatcher immediately decides
on a server to which the jobs will be assigned to and routes them accordingly;
(3) servers provide service. In particular, servers can process jobs that have
just arrived; (4) communication between the dispatcher and the servers may
occur. As a result, the dispatcher may update its state (e.g. which servers
have no pending jobs), on which it may base the assignment decision at the
next time slot.

Batch arrivals. The job arrival process is given by a sequence of _i.i.d._
random variables (RVs) $\\{a(t)\\}_{t=1}^{\infty}$. Namely, at each time slot
$t>0$, $a(t)$ jobs arrive as a batch at the dispatcher. We assume
$\mathbb{P}\left(a(1)=0\right)>0$ to avoid dealing with a periodic Markov
chain. We denote the mean and standard deviation of $a(1)$ by $\lambda$ and
$\sigma_{a}$, respectively.

Random service capacities. The service process of server
$i\in\\{1,\ldots,n\\}:=[n]$ is given by a sequence of _i.i.d._ RVs
$\\{s_{i}(t)\\}_{t=1}^{\infty}$, representing the potential number of jobs
that a server can complete in a time slot. We assume that these RVs are
positive and bounded, _i.e.,_ there exists $s_{max}\in\mathbbm{N}$ such that

$1\leq s_{i}(t)\leq s_{max},\quad\mbox{ for all }i\in[n]\mbox{ and }t>0.$ (1)

We denote the mean and standard deviation of $s_{i}(1)$ by $\mu_{i}$ and
$\sigma_{s_{i}}$, respectively. Assumption (1) implies

$1\leq\mu_{i}\leq s_{max},\quad\mbox{ for all }i\in[n].$ (2)

Queue lengths and idleness. Let $Q_{i}(t)$ denote the number of pending jobs
in queue $i$ at the end of time slot $t$, after arrival and departure. Let
$Q(t)=(Q_{1}(t),\ldots,Q_{n}(t))$, and $Q=(Q(t))_{t=0}^{\infty}$. We refer to
$Q(0)$ as the initial condition. Server $i$ is idle at time slot $t$ if
$Q_{i}(t)=0$. A server becomes idle at time slot $t>0$ if it had jobs to
process during time slot $t$ and processed all of them. For example, if
$Q_{i}(t-1)=0$, jobs arrived to server $i$ at the beginning of time slot $t$
and it served all of them, i.e. $Q_{i}(t)=0$, then it became idle at time slot
$t$. But, if $Q_{i}(t-1)=0$ and $Q_{i}(t)=0$ and no jobs arrived at server
$i$, it did not become idle at time slot $t$.

Tokens. A single token is associated with each server. Thus, the total number
of tokens in the system is $n$. Each token can be held by either the
corresponding server or the dispatcher. When a server becomes idle, it sends
its token to the dispatcher. Similarly, when the dispatcher chooses an idle
server as a destination for a job, it implicitly sends the token back to this
server along with the job. At time slot $t=0$ the dispatcher holds the tokens
corresponding to the servers that are initially idle.

Example 1. If $Q_{i}(t-1)>0$, then server $i$ must hold its token at the
beginning of time slot $t$. If, in addition, $Q_{i}(t)=0$, server $i$ becomes
idle at time slot $t$ and sends the token at the end of time slot $t$.

Example 2. If $Q_{i}(t-1)=0$, the dispatcher must hold server $i$’s token at
the beginning of time slot $t$. If, in addition, jobs arrive to server $i$ at
the beginning of time slot $t$ then the dispatcher sends the token to server.
If server $i$ now completes processing all jobs, it sends the token at the end
of time slot $t$. Thus, in this case, the token travels from the dispatcher
and back at the same time slot.

PI routing policy. The routing is encoded in a process that takes values in
$[n]$, denoted by $LI=(LI(t))_{t=0}^{\infty}$. If a job arrives at time $t>0$,
it joins queue $LI(t)$. We refer to the corresponding server as the Last-Idle
server. Roughly speaking, the Last-Idle server is the one that was idle most
recently. If there are multiple idle servers, the Last-Idle server is chosen
by some tie-breaking rule. For concreteness we assume that the tie-breaking
rule is that the Last-Idle server is chosen uniformly at random among the tied
servers. The exact mathematical definition of the process $LI$ appears below.
Thus, the policy persistently sends the arriving jobs to the Last-Idle server
until other servers become idle.

Markov chain formulation. The precise mathematical formulation is given in
terms of a Markov chain, that we now define. The Markov chain, denoted
$\\{X(t)\\}=\\{(Q(t),LI(t))\\}$, has state space

$\hat{\mathcal{A}}=\mathbb{N}^{n}\times[n].$

A member $\alpha$ of the state space $\hat{\mathcal{A}}$ is of the form

$\alpha=(\alpha_{1},\ldots,\alpha_{n},l),$

where $\alpha_{1},\ldots,\alpha_{n}$ are queue lengths, and $l$ is the
identity of the Last-Idle server.

Define the set of idle servers at the end of time slot $t$ by

$I(t)=\\{i:Q_{i}(t)=0\\}.$ (3)

The queue length and Last-Idle processes satisfy the following recursion for
$t>0$,

$\displaystyle\begin{cases}Q_{i}(t)=[Q_{i}(t-1)+a(t)\mathbbm{1}_{\\{LI(t)=i\\}}-s_{i}(t)]^{+}\\\
LI(t)=\begin{cases}LI(t-1),&\text{if }|I(t-1)|=0\\\
\xi_{t},&\text{otherwise,}\end{cases}\end{cases}$ (4)

where for each $t>0$, $\xi_{t}$ is a RV whose conditional distribution given

$(X(0),X(1),\ldots,X(t-1),a(t),s_{1}(t),\ldots,s_{n}(t))$

is uniform on the set $I(t-1)$ on the event $\\{|I(t-1)|>0\\}$, and
$\xi_{t}=1$ on the event $\\{|I(t-1)|=0\\}$ (note that the value of $\xi_{t}$
on the latter event is immaterial; it is specified only for concreteness). It
is clear from the definition that $X$ forms a Markov chain.

We define the set ${\cal{A}}\subset\hat{\cal{A}}$ as the collection of all
states $\alpha\in\hat{\cal{A}}$ such that $\alpha$ can be reached with a
positive probability from any of the states of the form $(0,\ldots,0,i)$,
$i\in[n]$, corresponding to an empty system. The initial condition
$X(0)=(Q(0),LI(0))$ is assumed to lie in $\cal{A}$ with probability 1. Under
this assumption, we may and will assume without loss of generality that the
state space of $X$ is ${\cal A}$.

By (4), the definition of $\cal{A}$ and the fact that the process can reach
$(0,\ldots,0,i)$ from $(0,\ldots,0,j)$ for any $i,j$ (thanks to the
assumptions that $\mathbb{P}\left(a(1)=0\right)>0$ and the uniformly-at-random
tie-breaking rule), it follows that $X$ is an irreducible, aperiodic,
homogeneous Markov chain on ${\cal A}$.

Departures. Let $d_{i}(t)$ denote the number of jobs completed by server $i$
during time slot $t$. Clearly, $d_{i}(t)\leq s_{i}(t)$. Moreover, $d_{i}(t)$
is given by

$d_{i}(t)=\min\\{s_{i}(t),Q_{i}(t-1)+a(t)\mathbbm{1}_{\\{LI(t)=i\\}}\\},$

and the first equation in recursion (4) can also be written as

$Q_{i}(t)=Q_{i}(t-1)+a(t)\mathbbm{1}_{\\{LI(t)=i\\}}-d_{i}(t).$ (5)

###  2.2 Persistent Idle with batch splitting

The PI policy sends all of the jobs belonging to a batch to a single server.
The batch-splitting PI policy (referred to as PI-Split) is a variation in
which whenever there are at least two idle servers, the dispatcher is allowed
to split an arriving batch into sub-batches and send them to multiple idle
servers. If there are none, the dispatcher sends the batch to the server
corresponding to the value of $LI(t)$.

Under PI-Split, the tokens are handled exactly as under PI: when a server
becomes idle, it sends its token to the dispatcher and when a dispatcher sends
jobs to a server (either a batch or a sub-batch) it implicitly returns the
token. Thus, under PI-Split, whenever a batch is split into sub-batches,
multiple tokens are returned to the corresponding servers.

The updating rule for $LI$ remains the same as under PI: If there are idle
servers, $LI$ is chosen uniformly at random. If there are none, its value
remains unchanged. We add a superscript $S$ for the batch splitting case such
that the queue lengths, idleness and Last-Idle processes are given by $Q^{S}$,
$I^{S}$ and $LI^{S}$ respectively. The state of the system is now given by the
process

$X^{S}=:(Q^{S},LI^{S}).$

The only difference between PI and PI-split is how the queue lengths are
updated whenever $|I^{S}(t-1)|>0$. To keep the Markov chain formulation simple
for this case, we assume that the dispatcher determines how to split the batch
between the idle servers based only on the current queue length state and the
current batch size.

Formally, let $\\{B(t)\\}_{t=1}^{\infty}$ be an i.i.d process taking values in
$\mathbbm{R}$. Let $v[i]$ denote the $i$th member of a vector $v$ in
$\mathbbm{N}^{n}$. Define a family of deterministic functions

$f:\mathbbm{N}^{n}\times\mathbbm{N}\times\mathbbm{R}\rightarrow\mathbbm{N}^{n}$
(6)

satisfying the following conditions:

$\displaystyle\sum_{i=1}^{n}f(\alpha_{1},\ldots,\alpha_{n},a,b)[i]=a,$ (7)
$\displaystyle f(\alpha_{1},\ldots,\alpha_{n},a,b)[i]=0\quad\text{ if
}\min_{1\leq j\leq n}{\alpha_{j}}=0\text{ and }\alpha_{i}>0.$ (8)

The inputs to $f$ are the queue lengths $Q^{S}$, the batch size $a(t)$ and the
value of the i.i.d process $B(t)$. The $i$th member of the output vector
equals the number of jobs the dispatcher should send to server $i$. The first
condition in (7) states that the total number of assigned jobs is exactly the
batch size. The second condition states that if there are idle servers
($\min_{1\leq j\leq n}\alpha_{j}=0$), then any server that is not idle should
not receive any jobs. The behaviour of $f$ in the case where there are no idle
servers is immaterial.

Following the definitions above, we can now state the updating rule for
$X^{S}$:

$\displaystyle\begin{cases}Q^{S}_{i}(t)=[Q^{S}_{i}(t-1)+a_{i}(t)-s_{i}(t)]^{+}\\\
LI^{S}(t)=\begin{cases}LI^{S}(t-1),&\text{if }|I^{S}(t-1)|=0\\\
\xi^{S}_{t},&\text{otherwise,}\end{cases}\\\
a_{i}(t)=\begin{cases}a(t),&\text{if }|I^{S}(t-1)|=0\text{ and }i=LI^{S}(t)\\\
f(X^{S}_{i}(t-1),a(t),B(t))[i],&\text{if }|I^{S}(t-1)|>0\text{ and }i\in
I^{S}(t-1)\\\ 0,&\text{otherwise}\end{cases}\end{cases}$ (9)

where $\xi^{S}_{t}$ is a RV distributed uniformly on the idle servers and is
defined analogously to $\xi_{t}$ in (4) and the function $f$ encodes the batch
splitting rule and belongs to the family of functions defined by (6) and (7).
Thus, $X^{S}$ is an irreducible, aperiodic, time homogeneous Markov chain.

##  3 Main results

###  3.1 Communication overhead

We begin with a simple result regarding the communication overhead of PI and
PI-Split. Both only rely on tokens passed from the servers to the dispatcher.

######  Proposition 3.1.

The communication overhead of PI (resp. PI-Split) consists of at most a single
message per batch (resp. sub-batch).

###### Proof.

A necessary condition for a transmission of a token is the completion of a
job. Since PI does not split batches, a completion of a job which results in a
server becoming idle must be a completion of a batch. Similarly, for PI-Split,
a message can be sent only when sub-batches are completed and not all
completions result in a message. ∎

###  3.2 Stability results

We introduce the following notation. For $\alpha\in\mathcal{A}$, we write
$\mathbb{E}_{\alpha}[\ {\cdot}\ ]$ and $\mathbb{P}_{\alpha}\left(\ {\cdot}\
\right)$ for the conditional expectation $\mathbb{E}[\ {\cdot}\
|X(0){=}\alpha]$ and probability $\mathbb{P}\left(\ \cdot\
|X(0){=}\alpha\right)$, respectively.

####  3.2.1 Overview

A useful approach to proving stability for irreducible Markov processes relies
on the construction of a Lyapunov function $\cal{L}$, mapping the state space
to $[0,\infty)$, that satisfies a drift condition within a single time step,
of the form

$\mathbb{E}_{\alpha}[{\cal{L}}(X(1))]-{\cal{L}}(\alpha)<-\varepsilon<0,$ (10)

such that the drift is negative and bounded away from zero at all states
$\alpha$ outside of some finite set $A\subset\cal{A}$, and a finite
expectation at all states $\alpha\in A$, namely
$\mathbb{E}_{\alpha}[{\cal{L}}(X(1))]<\infty$. The existence of such a
function ensures that provided the chain starts in the finite set $A$, the
expected return time of the chain to $A$ is finite. This, in turn, implies
several stability properties, such as positive recurrence and convergence to a
unique invariant measure (Lemma I.3.10 of [6]).

However, this approach does not appear to be directly applicable for the model
under consideration. As discussed at length in [1], the main difficulty arises
when considering states where one queue length is large while the others are
close to zero. It is possible that the large queue length belongs to the Last-
Idle server, which receives all of the work.

The proof in [1] relies on sampling the underlying Markov chain at (random)
times at which the identity of the Last-Idle server may change, namely, the
times at which a token is available at the dispatcher and it is not the last
used token. If the service capacities are deterministic (as in [1]), the next
such time is a function of the current state alone. Then, the calculation of
the drift is performed with respect to time intervals that depend
(deterministically) on the states the process traverses (as in [7], [8]).

At these sampling times the queue length vector is always at the boundary of
the positive orthant. The drift condition then only needs to be verified at
these special states, rather than at the whole space. Denote the set of these
states by ${\cal{A}}_{b}$ (where $b$ is mnemonic for boundary). Because the
Lyapunov function is observed at sampling times, the negative drift condition
must be stated in terms of the duration of the intervals. That is, in place of
(10), it is shown in [1] that

$\mathbb{E}_{\alpha}[{\cal{L}}(X(\tau))]-{\cal{L}}(\alpha)<-\varepsilon\tau,$
(11)

for all but finitely many $\alpha$ in the set ${\cal{A}}_{b}$, where
$\varepsilon>0$ and $\tau$ denotes the next sampling time.

The random service capacities we consider in this paper further complicate the
analysis. In this case, given the current state $X(0)=\alpha$, the next
sampling time $\tau$ is random. Instead of (11), our method is based on
showing that

$\mathbb{E}_{\alpha}[{\cal{L}}(X(\tau))]-{\cal{L}}(\alpha)<-\varepsilon\mathbb{E}_{\alpha}[\tau],$
(12)

for all but finitely many $\alpha$ in the set ${\cal{A}}_{b}$.

We analyze the stability of PI and PI-Split in Sections 3.3 and 3.4
respectively.

###  3.3 Stability of PI

######  Theorem 3.2.

Assume $\lambda<\sum_{i=1}^{n}\mu_{i}$. Then:
(i) $X$ is positive recurrent. Consequently,
(ii) $X$ has a unique stationary distribution, denoted by $\pi_{X}$, and
(iii) For any initial state $\alpha\in\mathcal{A}$ and any
$B\subset\mathcal{A}$, $\mathbb{P}_{\alpha}(X(t)\in B)\rightarrow\pi_{X}(B)$
as $t\rightarrow\infty$.

Proof. For $A\subset\cal{A}$ define the hitting time

$\tau^{A}=\inf{\\{t\geq 1;\,X(t)\in A\\}}.$ (13)

We prove that there exists a nonempty finite set $A\subset{\cal{A}}$ such that
for all $\alpha_{0}\in A$,

$\mathbbm{E}_{\alpha_{0}}[\tau^{A}]<\infty.$ (14)

By I.3.10 of [6], the irreducibility of $X$ and (14) imply positive recurrence
and hence part (i) of Theorem 3.2. Part (ii) then follows by I.3.6 of [6].
Finally, since the chain is aperiodic, part (iii) follows by the ergodic
theorem for Markov chains, Theorem I.4.2 of [6]. It thus remains to show (14).
To this end, we define a sampled chain $Y$ of $X$, sampled at the times when
the value of $LI$ may change.

To define the random sampling times, we first define the following
filtrations. The first is the natural filtration of $X$, namely

${\cal{F}}_{t}=\sigma(X(0),\ldots,X(t))$ (15)

and the second also includes the RVs $\xi_{1},\ldots,\xi_{t+1}$ which are used
to break ties between the idle servers, namely

${\cal{F}}^{\xi}_{t}=\sigma(X(0),\ldots,X(t),\xi_{1},\ldots,\xi_{t+1}).$ (16)

The filtration ${\cal{F}}^{\xi}_{t}$ in (16) is useful when knowing the
identity of the newly used token simplifies the analysis.

By (4), the value of $LI$ can only change at time $t$ if there is an idle
server which is not the Last-Idle server, corresponding to the event

$\Omega_{t}:=\\{\exists i\mbox{ s.t. }Q_{i}(t)=0,i\neq LI(t)\\}.$ (17)

Define the sequence of stopping times $\tau_{0},\tau_{1},\ldots$ by

$\displaystyle\tau_{0}=\inf\\{t\geq 1:\mathbbm{1}_{\Omega_{t}}=1\\}$ (18)
$\displaystyle\tau_{k+1}=\inf\\{t>\tau_{k}:\mathbbm{1}_{\Omega_{t}}=1\\},\quad
k=0,1,2\ldots$ (19)

Since the values of $\tau^{A},\tau_{0},\tau_{1},\ldots$ only depend on the
values of $X$, the former are stopping times with respect to both filtrations,
${\cal{F}}_{t}$ and ${\cal{F}}^{\xi}_{t}$. By the definition of the policy,
the finiteness of the initial condition and the finiteness of the arrival
process, the stopping times $\tau_{k}$ are finite a.s. for all $k\geq 0$ (the
proof follows the same argument as in Lemma 2 of [1] and hence is omitted).
Define the sampled chain

$Y_{k}=X({\tau_{k}}),\quad k\geq 0.$ (20)

Denote by ${\cal{A}}^{*}\subset{\cal{A}}$ the collection of states that can be
reached by $Y$, namely

$\displaystyle{\cal{A}}^{*}=\\{\alpha=(\alpha_{1},\ldots,\alpha_{n},l)\in{\cal{A}}:\alpha_{i}=0\text{
for some }i\neq l\\}.$

The following lemma asserts that a certain drift criterion is fulfilled by
$Y$.

######  Lemma 3.3.

There exist a constant $\gamma$ such that $0<\gamma\leq 1$, a function
${\cal{L}}:{\cal{A}}\rightarrow\mathbbm{R}_{+}$ of the form

${\cal{L}}(\alpha)=\sum_{i=1}^{n}\alpha_{i}^{1+\gamma},$ (21)

a positive constant $\epsilon$ and a finite non-empty set
$A\subset{\mathcal{A}}$, such that for all $k\geq 0$ and
$\alpha\in{\cal{A}}^{*}\setminus A$, the following inequality holds:

$\displaystyle\mathbb{E}_{\alpha}[({\cal{L}}(Y_{k+1})-{\cal{L}}(Y_{k}))\mathbbm{1}_{\\{k<\sigma\\}}\mid{\cal{F}}^{\xi}_{\tau_{k}}]\leq-\epsilon\mathbb{E}_{\alpha}[(\tau_{k+1}-\tau_{k})\mathbbm{1}_{\\{k<\sigma\\}}\mid{\cal{F}}^{\xi}_{\tau_{k}}],$
(22)

where

$\sigma=\inf\\{k\geq 0:Y_{k}\in A\\}.$ (23)

The stopping time $\sigma$ counts the number of samples needed until $Y$ hits
the set $A$ for the first time. Inequality (22) states that, given the initial
condition is outside of $A$, as long as $k<\sigma$, the drift of the function
$\cal{L}$ applied to $Y$ is strictly negative and is proportional to the
expected duration of the intervals between samples. The proof is provided
below the proof of the theorem.

Now, fix an initial state $\alpha_{0}$ such that $X(0)=\alpha_{0}\in A$. By
the definition of the sampling times $\sigma,\tau_{0},\tau_{1},\ldots$ in
(18), we have

$\tau_{\sigma}=\inf\\{t\geq 1:X(t)\in{\cal{A}}^{*}\cap A\\},$ (24)

such that $X(\tau_{\sigma})=Y(\sigma)\in{\cal{A}^{*}}\cap A$. The value of
$\tau_{\sigma}$ only depends on the values of $X$, hence $\tau_{\sigma}$ is a
stopping time with respect to ${\cal{F}}_{t}$ and ${\cal{F}}^{\xi}_{t}$. While
$X$ reaches $A$ at time slot $\tau_{\sigma}$, it is not necessarily the first
time slot in $[1,\tau_{\sigma}]$ in which this occurs. Since $\tau^{A}$,
defined in (13), equals the number of time slots it takes $X$ to reach the set
$A$ for the first time, we must have $\tau^{A}\leq\tau_{\sigma}$ and therefore

$\mathbbm{E}_{\alpha_{0}}[\tau^{A}]\leq\mathbbm{E}_{\alpha_{0}}[\tau_{\sigma}].$
(25)

Our goal is now to prove that $\mathbbm{E}_{\alpha_{0}}[\tau_{\sigma}]$ is
finite. We have:

$\displaystyle\mathbbm{E}_{\alpha_{0}}[\tau_{\sigma}]$
$\displaystyle=\sum_{\alpha\in{\cal{A}^{*}}\cap
A}\mathbbm{E}_{\alpha_{0}}[\tau_{\sigma}\mid
X(\tau_{0})=\alpha]\mathbbm{P}_{\alpha_{0}}(X(\tau_{0})=\alpha)$ (26)
$\displaystyle\quad+\sum_{\alpha\in{\cal{A}^{*}}\setminus
A}\mathbbm{E}_{\alpha_{0}}[\tau_{\sigma}\mid
X(\tau_{0})=\alpha]\mathbbm{P}_{\alpha_{0}}(X(\tau_{0})=\alpha).$ (27)

Conditioned on the event $\\{X(\tau_{0})=\alpha\in{\cal{A}^{*}}\cap A\\}$, by
the definitions of $\tau_{0}$ in (18) and $\tau_{\sigma}$ in (24), we have
$\tau_{\sigma}=\tau_{0}$ and therefore

$\mathbbm{E}_{\alpha_{0}}[\tau_{\sigma}\mid
X(\tau_{0})=\alpha]=\mathbbm{E}_{\alpha_{0}}[\tau_{0}\mid
X(\tau_{0})=\alpha],\quad\alpha\in{\cal{A}^{*}}\cap A.$ (28)

We proceed with analyzing the term $\mathbbm{E}_{\alpha_{0}}[\tau_{\sigma}\mid
X(\tau_{0})=\alpha]$ for the case $\alpha\in{\cal{A}^{*}}\setminus A$. First,
we prove that

$\mathbbm{E}_{\alpha_{0}}[\tau_{\sigma}\mid{\cal{F}}_{\tau_{0}}]=\tau_{0}+g(X(\tau_{0})),$
(29)

where $g(x):=\mathbbm{E}_{x}[\tau_{\sigma}]$.

Since $\tau_{0}$ and $\tau_{\sigma}$ are stopping times with respect to
${\cal{F}}_{t}$ and $\tau_{0}\leq\tau_{\sigma}$ a.s., by the Strong Markov
Property (Theorem 5.2.5. of [9]), we have

$\mathbbm{E}_{\alpha_{0}}[(\tau_{\sigma}-\tau_{0})\wedge
r\mid{\cal{F}}_{\tau_{0}}]=\phi^{r}(X(\tau_{0})),$ (30)

where $r$ is a positive constant and
$\phi^{r}(x):=\mathbbm{E}_{x}[\tau_{\sigma}\wedge r]$. By the monotone
convergence theorem we have

$\mathbbm{E}_{x}[\tau_{\sigma}\wedge
r]\xrightarrow{r\rightarrow\infty}\mathbbm{E}_{x}[\tau_{\sigma}]$ (31)

and

$\mathbbm{E}_{\alpha_{0}}[(\tau_{\sigma}-\tau_{0})\wedge
r\mid{\cal{F}}_{\tau_{0}}]\xrightarrow{r\rightarrow\infty}\mathbbm{E}_{\alpha_{0}}[(\tau_{\sigma}-\tau_{0})\mid{\cal{F}}_{\tau_{0}}].$
(32)

Taking the limit w.r.t. $r$ in both sides of (30) and using (31) and (32)
yields (29).

Since $X(\tau_{0})$ is measurable with respect to ${\cal{F}}_{\tau_{0}}$,
taking the conditional expectation of both sides of (29), conditioned on
$X(\tau_{0})$, yields

$\mathbbm{E}_{\alpha_{0}}[\mathbbm{E}_{\alpha_{0}}[\tau_{\sigma}\mid{\cal{F}}_{\tau_{0}}]\mid
X(\tau_{0})]=\mathbbm{E}_{\alpha_{0}}[\tau_{\sigma}\mid
X(\tau_{0})]=\mathbbm{E}_{\alpha_{0}}[\tau_{0}\mid
X(\tau_{0})]+g(X(\tau_{0})),$

which implies

$\mathbbm{E}_{\alpha_{0}}[\tau_{\sigma}\mid
X(\tau_{0})=\alpha]=\mathbbm{E}_{\alpha_{0}}[\tau_{0}\mid
X(\tau_{0})=\alpha]+\mathbbm{E}_{\alpha}[\tau_{\sigma}].$ (33)

Substituting (28) and (33) in (26) yields

$\displaystyle\mathbbm{E}_{\alpha_{0}}[\tau_{\sigma}]$
$\displaystyle=\sum_{\alpha\in{\cal{A}^{*}}}\mathbbm{E}_{\alpha_{0}}[\tau_{0}\mid
X(\tau_{0})=\alpha]\mathbbm{P}_{\alpha_{0}}(X(\tau_{0})=\alpha)$ (34)
$\displaystyle\quad+\sum_{\alpha\in{\cal{A}^{*}}\setminus
A}\mathbbm{E}_{\alpha}[\tau_{\sigma}]\mathbbm{P}_{\alpha_{0}}(X(\tau_{0})=\alpha)$
(35)
$\displaystyle=\mathbbm{E}_{\alpha_{0}}[\tau_{0}]+\sum_{\alpha\in{\cal{A}^{*}}\setminus
A}\mathbbm{E}_{\alpha}[\tau_{\sigma}]\mathbbm{P}_{\alpha_{0}}(X(\tau_{0})=\alpha),$
(36)

where the last equality is due to fact that the value of $X(\tau_{0})$ is
almost surely in ${\cal{A}}^{*}$. We proceed by showing that the two terms in
the right hand side of (34) are finite, which, combined with (25) proves (14)
and concludes the proof.

First, we consider the term $\mathbbm{E}_{\alpha_{0}}[\tau_{0}]$. By (18), the
stopping time $\tau_{0}$ is the first time after time 0 such that one of the
queue lengths, corresponding to a server different from $LI(0)$, equals zero.
Under PI, the servers different from $LI(0)$ do not receive any work during
$[1,\tau_{0}]$. By (1), a non-idle server must serve at least one job during a
time slot. Thus,

$\tau_{0}\leq\max_{i}\\{Q_{i}(0):i\neq LI(0)\\}\leq\max_{i}\\{Q_{i}(0)\\}$
(37)

and therefore

$\mathbbm{E}_{\alpha_{0}}[\tau_{0}]\leq\max_{i}\\{Q_{i}(0)\\}<\infty,$ (38)

where the last inequality is due to the finitness of the initial condition.

Next, we consider the second term in the right hand side of (34). The proof
proceeds by deriving an upper bound on $\mathbbm{E}_{\alpha}[\tau_{\sigma}]$
for $\alpha\in{\cal{A}^{*}}\setminus A$ by using the sampled chain $Y$,
defined in (20), and Lemma 3.3. Taking expectation on both sides of (22)
yields

$\mathbb{E}_{\alpha}[({\cal{L}}(Y_{k+1})-{\cal{L}}(Y_{k}))\mathbbm{1}_{\\{k<\sigma\\}}]\leq-\epsilon\mathbb{E}_{\alpha}[(\tau_{k+1}-\tau_{k})\mathbbm{1}_{\\{k<\sigma\\}}].$
(39)

We use (39) twice: to prove that $\tau_{\sigma}$ is well defined, namely, that
$\sigma<\infty$ a.s., and then to derive a bound for
$\mathbbm{E}_{\alpha}[\tau_{\sigma}]$.

By the definition of the sampling times in (18), we have
$\tau_{k+1}-\tau_{k}\geq 1$. Thus (39) implies that

$\mathbb{E}_{\alpha}[({\cal{L}}(Y_{k+1})-{\cal{L}}(Y_{k}))\mathbbm{1}_{\\{k<\sigma\\}}]\leq-\epsilon\mathbb{E}_{\alpha}[\mathbbm{1}_{\\{k<\sigma\\}}].$

Summing over $k\in[0,m-1]$ for $m\geq 1$ yields

$\mathbb{E}_{\alpha}[{\cal{L}}(Y_{m\wedge\sigma})]-\mathbb{E}_{\alpha}[{\cal{L}}(Y_{0})]\leq-\epsilon\sum_{k=0}^{m-1}\mathbb{E}_{\alpha}[\mathbbm{1}_{\\{k<\sigma\\}}]=-\epsilon\sum_{k=0}^{m-1}\mathbbm{P}_{\alpha}(\sigma>k),$

After rearranging we obtain

$\sum_{k=0}^{m-1}\mathbbm{P}_{\alpha}(\sigma>k)\leq\mathbb{E}_{\alpha}[{\cal{L}}(Y_{0})]/\epsilon-\mathbb{E}_{\alpha}[{\cal{L}}(Y_{m\wedge\sigma})]/\epsilon\leq\mathbb{E}_{\alpha}[{\cal{L}}(Y_{0})]/\epsilon=\mathbb{E}_{\alpha}[{\cal{L}}(X(\tau_{0}))]/\epsilon,$

where in the second inequality we used the fact that
$\mathbb{E}_{\alpha}[{\cal{L}}(Y_{m\wedge\sigma})]\geq 0$. Hence

$\mathbb{E}_{\alpha}[\sigma]=\lim_{m\rightarrow\infty}\sum_{k=0}^{m-1}\mathbbm{P}_{\alpha}(\sigma>k)\leq\mathbb{E}_{\alpha}[{\cal{L}}(X(\tau_{0}))]/\epsilon.$
(40)

We now argue that $\mathbb{E}_{\alpha}[{\cal{L}}(X(\tau_{0}))]$ is finite for
all initial states $\alpha\in{\cal{A}}$ (and therefore, specifically, for
$\alpha\in{\cal{A}}^{*}\setminus A$). Using the definition of $\cal{L}$ in
(21), we have

$\mathbbm{E}_{\alpha}[{\cal{L}}(X(\tau_{0}))]=\mathbbm{E}_{\alpha}\Big{[}\sum_{i=1}^{n}(Q_{i}(\tau_{0}))^{1+\gamma}\Big{]}\leq\mathbbm{E}_{\alpha}\Big{[}\sum_{i=1}^{n}(Q_{i}(\tau_{0}))^{2}\Big{]},$
(41)

where the last inequality is due to the fact that $0<\gamma\leq 1$. Denote the
maximum queue length at time $t$ by $Q_{M}(t)$, namely,
$Q_{M}(t)=\max_{i}\\{Q_{i}(t)\\}$. The value of $Q_{M}(0)$ is finite for all
initial states $\alpha\in{\cal{A}}$. The quantities $Q_{i}(\tau_{0})$ for
servers different from the Last Idle server are each bounded from above by
$Q_{M}(0)$. The queue length of the Last Idle server at $\tau_{0}$ is bounded
from above by $Q_{M}(0)$ plus the number of jobs that arrived to the system
during $[1,\tau_{0}]$. Therefore, we have

$\displaystyle\mathbbm{E}_{\alpha}\Big{[}\sum_{i=1}^{n}(Q_{i}(\tau_{0}))^{2}\Big{]}$
$\displaystyle\leq\mathbbm{E}_{\alpha}\Bigg{[}(n-1)Q_{M}^{2}(0)+\Big{(}Q_{M}(0)+\sum_{t=1}^{\tau_{0}}a(t)\Big{)}^{2}\Bigg{]}$
(42)
$\displaystyle\leq\mathbbm{E}_{\alpha}\Bigg{[}(n-1)Q_{M}^{2}(0)+\Big{(}Q_{M}(0)+\sum_{t=1}^{Q_{M}(0)}a(t)\Big{)}^{2}\Bigg{]}$
(43) $\displaystyle\leq
B\mathbbm{E}_{\alpha}\Big{[}Q_{M}^{2}(0)\Big{]}=B\Big{(}\max_{i\in[n]}\alpha_{i}\Big{)}^{2},$
(44)

for some constant $B>0$, where the second inequality is due to (37) and the
third inequality is due to the finite first and second moments of the arrival
process and the fact that $Q_{M}(0)$ is an integer. Using (42) in (41) and
recalling that $Q_{max}$ is finite yields

$\mathbbm{E}_{\alpha}[{\cal{L}}(X(\tau_{0}))]\leq
B\Big{(}\max_{i\in[n]}\alpha_{i}\Big{)}^{2}<\infty,\quad\alpha\in{\cal{A}}.$
(45)

Equations (40) and (45) establish that $\mathbbm{E}_{\alpha}[\sigma]<\infty$
and consequently

$\sigma<\infty\quad\mbox{a.s.}$ (46)

Now, revisiting (39), summing over $k\in[0,m-1]$ yields

$\mathbb{E}_{\alpha}[{\cal{L}}(Y_{m\wedge\sigma})]-\mathbb{E}_{\alpha}[{\cal{L}}(Y_{0})]\leq-\epsilon\mathbb{E}_{\alpha}[\tau_{m\wedge\sigma}].$

Using the fact that $\mathbb{E}_{\alpha}[{\cal{L}}(Y_{m\wedge\sigma})]\geq 0$
and rearranging yields

$\mathbb{E}_{\alpha}[\tau_{m\wedge\sigma}]\leq\mathbb{E}_{\alpha}[{\cal{L}}(X(\tau_{0}))]/\epsilon.$
(47)

Using (46), the monotone convergence theorem and (47), we obtain

$\mathbb{E}_{\alpha}[\tau_{\sigma}]=\mathbb{E}_{\alpha}[\lim_{m\rightarrow\infty}\tau_{m\wedge\sigma}]=\lim_{m\rightarrow\infty}\mathbb{E}_{\alpha}[\tau_{m\wedge\sigma}]\leq\mathbb{E}_{\alpha}[{\cal{L}}(X(\tau_{0}))]/\epsilon\leq\bar{B}\Big{(}\max_{i\in[n]}\alpha_{i}\Big{)}^{2},$
(48)

where $\bar{B}>0$ is some constant and in the last inequality we used (45).
Using (25), (34) and (48) yields

$\displaystyle\mathbbm{E}_{\alpha_{0}}[\tau^{A}]\leq\mathbbm{E}_{\alpha_{0}}[\tau_{\sigma}]$
$\displaystyle\leq\mathbbm{E}_{\alpha_{0}}[\tau_{0}]+\bar{B}\sum_{\alpha\in{\cal{A}^{*}}\setminus
A}\Big{(}\max_{i\in[n]}\alpha_{i}\Big{)}^{2}\mathbbm{P}_{\alpha_{0}}(X(\tau_{0})=\alpha)$
(49)

Recalling the definition of $Q_{M}$, the second term on the right hand side of
(49) satisfies

$\sum_{\alpha\in{\cal{A}^{*}}\setminus
A}\Big{(}\max_{i\in[n]}\alpha_{i}\Big{)}^{2}\mathbbm{P}_{\alpha_{0}}(X(\tau_{0})=\alpha)\leq\sum_{\alpha\in{\cal{A}^{*}}}\Big{(}\max_{i\in[n]}\alpha_{i}\Big{)}^{2}\mathbbm{P}_{\alpha_{0}}(X(\tau_{0})=\alpha)=\mathbbm{E}_{\alpha_{0}}\Big{[}Q_{M}^{2}(\tau_{0})\Big{]}.$

Therefore

$\mathbbm{E}_{\alpha_{0}}[\tau^{A}]\leq\mathbbm{E}_{\alpha_{0}}[\tau_{0}]+\bar{B}\mathbbm{E}_{\alpha_{0}}\Big{[}Q_{M}^{2}(\tau_{0})\Big{]}.$
(50)

Combining (38), (45) and (50) proves (14), which concludes the proof. ∎

Proof of Lemma 3.3.

First, we discuss the term $\mathbbm{1}_{\\{k<\sigma\\}}$ in (22). By (20), we
have $(Y_{0},\ldots,Y_{k})=(X({\tau_{0}}),\ldots,X({\tau_{k}}))$. By (16),
${\cal{F}}^{\xi}_{\tau_{k}}$ contains the entire evolution of the process $X$
during $[0,\tau_{k}]$ and therefore the values of the sampled chain $Y$,
$Y_{0},\ldots,Y_{k}$, are measurable with respect to
${\cal{F}}^{\xi}_{\tau_{k}}$. Since $\sigma$ (defined in (23)) corresponds to
the first sample which lies in the finite set $A$, the values
$Y_{0},\ldots,Y_{k}$ are enough to determine whether or not the event
$\\{k<\sigma\\}$ occurred. Thus, the random variable
$\mathbbm{1}_{\\{k<\sigma\\}}$ is also measurable with respect to
${\cal{F}}^{\xi}_{\tau_{k}}$, and we can write (22) as

$\displaystyle\mathbbm{1}_{\\{k<\sigma\\}}\mathbb{E}_{\alpha}[({\cal{L}}(Y_{k+1})-{\cal{L}}(Y_{k}))\mid{\cal{F}}^{\xi}_{\tau_{k}}]\leq-\epsilon\mathbbm{1}_{\\{k<\sigma\\}}\mathbb{E}_{\alpha}[(\tau_{k+1}-\tau_{k})\mid{\cal{F}}^{\xi}_{\tau_{k}}],$
(51)

Clearly, on the event $\\{k\geq\sigma\\}$, the two sides of inequality (51)
equal zero and (22) holds. The rest of proof is dedicated to proving (22) on
the event $\\{k<\sigma\\}$, namely, given that the values of
$Y_{0},\ldots,Y_{k}$ are all outside of $A$. For simplicity, we drop the
subscript $\alpha$ from the expected value notation and write $\mathbbm{E}$
for $\mathbbm{E}_{\alpha}$ for the rest of the proof.

We begin by fixing $\epsilon>0$ and a function $\cal{L}$. We then consider a
finite set $A$ of the form

$A=\\{\alpha=(\alpha_{1},\ldots,\alpha_{n},l)\in{\cal{A}}\mid\sum_{i}\alpha_{i}\leq\cal{C}\\},$
(52)

and show that if $\cal{C}$ is large enough, (22) holds.

Since $\lambda<\sum_{i=1}^{n}\mu_{i}$, there exists $\epsilon_{0}>0$ such that

$\lambda-\sum_{i=1}^{n}\mu_{i}=-\epsilon_{0}.$

Denote $\mu_{min}=\min_{i}\mu_{i}$ and let

$\epsilon=(\epsilon_{0}\wedge\mu_{min})/4.$ (53)

Fix $0<\gamma\leq 1$ such that

$\max_{j\in\\{1,\ldots,n\\}}\Big{\\{}\big{(}[\lambda-\mu_{j}]^{+}\big{)}^{1+\gamma}-\sum_{i\neq
j}\mu_{i}\Big{\\}}<-2\epsilon.$ (54)

To prove the existence of such a $\gamma$, define
$g_{j}(\gamma)=([\lambda-\mu_{j}]^{+})^{1+\gamma}-\sum_{i\neq j}\mu_{i}$. Now,

$\displaystyle g_{j}(0)=[\lambda-\mu_{j}]^{+}-\sum_{i\neq
j}\mu_{i}=-\big{(}\epsilon_{0}\wedge\sum_{i\neq
j}\mu_{i}\big{)}\leq-\big{(}\epsilon_{0}\wedge\mu_{min}\big{)}=-4\epsilon.$
(55)

Define $h(\gamma)=\max_{j\in\\{1,\ldots,n\\}}g_{j}(\gamma)$. By (55),
$h(0)\leq-4\epsilon$. Since $g_{j}(\gamma)$ are continuous functions of
$\gamma$, so is $h(\gamma)$. Thus, there exists $\gamma\in(0,1]$ such that
$h(\gamma)<-2\epsilon$, and (54) holds.

Recall (by (21)) the function ${\cal{L}}:{\cal{A}}\rightarrow\mathbbm{R}_{+}$,
for $\alpha=(\alpha_{1},\ldots,\alpha_{n},l)\in{\cal{A}}$, is defined by

${\cal{L}}(\alpha)=\sum_{i=1}^{n}\alpha_{i}^{1+\gamma}.$ (56)

Starting with the left hand side of (22), we use (56) to obtain

$\displaystyle\mathbb{E}[{\cal{L}}(Y_{k+1})-{\cal{L}}(Y_{k})\mid{\cal{F}}^{\xi}_{\tau_{k}}]=\mathbb{E}\Big{[}\sum_{i=1}^{n}\Big{(}Q_{i}^{1+\gamma}(\tau_{k+1})-{Q_{i}^{1+\gamma}(\tau_{k})}\Big{)}\mid{\cal{F}}^{\xi}_{\tau_{k}}\Big{]}.$
(57)

The proof proceeds by decomposing the sum in the right hand side of (57) into
two parts and analyzing them separately: the first part corresponds to the
Last Idle server at time $\tau_{k}+1$, which we denote by

$l_{k}=LI(\tau_{k}+1).$ (58)

The second part corresponds to the rest of the servers.

By (4), the value of $l_{k}$ is determined using $X(\tau_{k})$ and
$\xi_{\tau_{k}+1}$, which, by (16), are included in
${\cal{F}}^{\xi}_{\tau_{k}}$. Thus, the random variable $l_{k}$ is measurable
with respect to ${\cal{F}}^{\xi}_{\tau_{k}}$. By the definition of the
sampling times in (17) and (18), the identity of the Last Idle server remains
unchanged during $[\tau_{k}+1,\tau_{k+1}]$ and equals to $l_{k}$. Also,

$Q_{l_{k}}(\tau_{k})=0.$ (59)

Thus, during $[\tau_{k}+1,\tau_{k+1}]$, the server corresponding to $l_{k}$
receives all incoming jobs and its queue length is given by a reflected random
walk starting at zero. The rest of the queue lengths decrease until the first
reaches zero at time $\tau_{k+1}$. Therefore, using the balance equation (4),
we obtain

$Q_{i}(\tau_{k+1})=\Bigg{[}Q_{i}(\tau_{k})-\sum_{t=\tau_{k}+1}^{\tau_{k+1}}s_{i}(t)\Bigg{]}^{+},\quad
i\in\\{1,\ldots,n\\}\setminus\\{l_{k}\\}.$ (60)

For simplicity, denote

$\displaystyle Q_{i}^{k}=Q_{i}(\tau_{k}),\quad i\in\\{1,\ldots,n\\},$ (61)
$\displaystyle R_{i}^{k}=\sum_{t=\tau_{k}+1}^{\tau_{k+1}}s_{i}(t),\quad
i\in\\{1,\ldots,n\\}\setminus\\{l_{k}\\}.$ (62)

Using (58), (59), (60) and (61), we decompose the sum in the right hand side
of (57) as follows:

$\displaystyle\mathbb{E}\Big{[}\sum_{i=1}^{n}\Big{(}Q_{i}^{1+\gamma}(\tau_{k+1})-{Q_{i}^{1+\gamma}(\tau_{k})}\Big{)}\mid{\cal{F}}^{\xi}_{\tau_{k}}\Big{]}$
$\displaystyle=\mathbb{E}\Big{[}\sum_{i\neq
l_{k}}^{n}\Big{(}([Q_{i}^{k}-R_{i}^{k}]^{+})^{1+\gamma}-{(Q_{i}^{k})^{1+\gamma}}\Big{)}\mid{\cal{F}}^{\xi}_{\tau_{k}}\Big{]}$
(63)
$\displaystyle\quad+\mathbb{E}\Big{[}\Big{(}{Q_{l_{k}}(\tau_{k+1})}\Big{)}^{1+\gamma}\mid{\cal{F}}^{\xi}_{\tau_{k}}\Big{]}.$
(64)

Denote

$\Delta_{k}=\tau_{k+1}-\tau_{k}.$ (65)

We proceed by bounding the two terms in the right hand side of (63)
separately. We begin with
$\mathbb{E}\Big{[}\Big{(}{Q_{l_{k}}(\tau_{k+1})}\Big{)}^{1+\gamma}\mid{\cal{F}}^{\xi}_{\tau_{k}}\Big{]}$.
Conditioned on ${\cal{F}}^{\xi}_{\tau_{k}}$, the duration of the interval
$[\tau_{k}+1,\tau_{k+1}]$, given by $\Delta_{k}$, depends only on the queue
lengths at the servers other than $l_{k}$ and is therefore independent of the
reflected random walk in server $l_{k}$. Denote

$\displaystyle\beta_{i}$ $\displaystyle:=\lambda-\mu_{i}$ (66)
$\displaystyle\sigma_{i}^{2}$
$\displaystyle:=\sigma_{a}^{2}+\sigma_{s_{i}}^{2}.$ (67)

Using a modification of Lemma 6 of [1], we obtain

$\mathbb{E}\Big{[}\Big{(}{Q_{l_{k}}(\tau_{k+1})}\Big{)}^{1+\gamma}\mid{\cal{F}}^{\xi}_{\tau_{k}}\Big{]}\leq([\lambda-\mu_{l_{k}}]^{+})^{1+\gamma}\mathbb{E}[\Delta_{k}^{1+\gamma}\mid{\cal{F}}^{\xi}_{\tau_{k}}]+C\mathbb{E}[\Delta_{k}^{0.75(1+\gamma)}\mid{\cal{F}}^{\xi}_{\tau_{k}}],$
(68)

where
$C=\max_{i}\\{\big{(}16\sigma_{i}+8(\sigma_{i})^{1/2}(\beta_{i})^{+}\big{)}^{(1+\gamma)/2}\\}$.
The proof is deferred to the Appendix.

Next, we consider the first term in the right hand side of (63). For ease of
exposition, we suppress the index $k$ and write $Q_{i}$ and $R_{i}$ for
$Q_{i}^{k}$ and $R_{i}^{k}$ respectively. We have

$\displaystyle([Q_{i}-R_{i}]^{+})^{1+\gamma}-{Q_{i}^{1+\gamma}}$
$\displaystyle=Q_{i}^{1+\gamma}(([1-R_{i}/(Q_{i}\vee 1)]^{+})^{1+\gamma}-1)$
(69) $\displaystyle\quad\leq Q_{i}^{1+\gamma}([1-R_{i}/(Q_{i}\vee 1)]^{+}-1)$
(70) $\displaystyle=-Q_{i}^{\gamma}(Q_{i}\wedge R_{i}),$ (71)

where the inequality is due to the fact that $R_{i}/(Q_{i}\vee 1)\geq 0$ and
thus $0\leq[1-R_{i}/(Q_{i}\vee 1)]^{+}\leq 1$. Combining (57), (63), (68) and
(69) we obtain

$\displaystyle\mathbb{E}[{\cal{L}}(Y_{k+1})-{\cal{L}}(Y_{k})\mid{\cal{F}}^{\xi}_{\tau_{k}}]\leq-\epsilon\mathbb{E}[\Delta_{k}\mid{\cal{F}}^{\xi}_{\tau_{k}}]+\Phi_{k},$
(72)

where

$\Phi_{k}:=\mathbb{E}\Big{[}\epsilon\Delta_{k}+([\lambda-\mu_{l_{k}}]^{+})^{1+\gamma}\Delta_{k}^{1+\gamma}+C\Delta_{k}^{0.75(1+\gamma)}-\sum_{i\neq
l_{k}}^{n}Q_{i}^{\gamma}(Q_{i}\wedge
R_{i})\mid{\cal{F}}^{\xi}_{\tau_{k}}\Big{]}.$ (73)

Thus, using (72), (22) will follow once we prove that $\Phi_{k}\leq 0$ a.s. To
this end, denote

$\displaystyle Q_{min}=\min_{1\leq i\leq n,i\neq l_{k}}Q_{i},$ $\displaystyle
Q_{max}=\max_{1\leq i\leq n}Q_{i}.$

Denote by $i_{max}$ the server with the maximum queue length at time
$\tau_{k}$, namely

$i_{max}=\mbox{argmax}_{1\leq i\leq n}\\{Q_{i}\\}$

Since the queue length at time $\tau_{k}$ of the Last Idle server $l_{k}$
equals zero, we have $i_{max}\neq l_{k}$.

The key to analyzing (73) is to relate $\Delta_{k}$, the duration of the
interval $[\tau_{k}+1,\tau_{k+1}]$, to the value of $Q_{min}$, which, by (16),
is measurable with respect to ${\cal{F}}^{\xi}_{\tau_{k}}$. By the definition
of the sampling times in (18) we have $\Delta_{k}\geq 1$. By (1), any strictly
positive queue lengths of the servers different from $l_{k}$ must decrease by
at least one and by at most $s_{max}$ at each time slot. Therefore, given
${\cal{F}}^{\xi}_{\tau_{k}}$, we have

$\displaystyle(Q_{min}/s_{max})\vee 1\leq\Delta_{k}\leq Q_{min}\vee 1$ (74)

The argument proceeds by analyzing (73) in two cases, corresponding to the
value of $Q_{min}$. The idea is that if $Q_{min}$ is small, then by (74), so
is $\Delta_{k}$. Thus, if the finite set in (52) is chosen large enough so
that $Q_{max}$ must be very large, one of the members of $\sum_{i\neq
l_{k}}^{n}Q_{i}^{\gamma}(Q_{i}\wedge R_{i})$ is very large compared to the
other terms in the right hand side of (73), resulting in $\Phi_{k}\leq 0$.
Otherwise, if $Q_{min}$ is large, then by (74) so is $\Delta_{k}$. In this
case, we prove that

$\mathbb{E}\Big{[}\sum_{i\neq l_{k}}^{n}Q_{i}^{\gamma}(Q_{i}\wedge
S_{i})\mid{\cal{F}}^{\xi}_{\tau_{k}}\Big{]}\geq\mathbb{E}\Big{[}\Big{(}\sum_{i\neq
l_{k}}^{n}\mu_{i}\Big{)}\Delta_{k}^{1+\gamma}-(n-1)s_{max}\Delta_{k}^{\gamma}\mid{\cal{F}}^{\xi}_{\tau_{k}}\Big{]}.$

Thus, if we use this to further bound $\Phi_{k}$ from above, by (54), the
coefficient of the leading term $\Delta_{k}^{1+\gamma}$ is negative. If
$\Delta_{k}$ is large enough we again obtain $\Phi_{k}\leq 0$.

Define

$f(x)=\epsilon x-2\epsilon
x^{1+\gamma}+Cx^{0.75(1+\gamma)}+(n-1)s_{max}x^{\gamma},$ (75)

where $C$ is from (68). Since the coefficient of the leading term
$x^{1+\gamma}$ in $f(x)$ is negative, there exists $u>0$ such that

$u>s^{2}_{max}\quad\mbox{ and }\quad f(x)<0\quad\mbox{ if }x>u/s_{max}.$ (76)

Fix such $u$.
Choose the constant $\cal{C}$ in the definition of the finite set $A$ in (52)
such that

${\cal{C}}>ns_{max}u+n\Big{(}\epsilon
u+(\max_{i}\\{[\lambda-\mu_{i}]^{+}\\})^{1+\gamma}u^{1+\gamma}+Cu^{0.75(1+\gamma)}\Big{)}^{1/\gamma}+n\Big{(}\epsilon+n(\sigma_{a}^{2}+\lambda^{2})\Big{)}^{1/\gamma}.$
(77)

By (52), for any state outside of $A$ we have ${\cal{C}}<\sum_{i}Q_{i}\leq
nQ_{max}$. Thus, by (77), we obtain

$Q_{max}>s_{max}u+\Big{(}\epsilon
u+(\max_{i}\\{[\lambda-\mu_{i}]^{+}\\})^{1+\gamma}u^{1+\gamma}+Cu^{0.75(1+\gamma)}\Big{)}^{1/\gamma}.$
(78)

We now consider two cases that correspond to the value of $Q_{min}$, which, by
(74), is related to the duration of the interval. The first two members of the
right hand side of (77) correspond to our requirements of the finite set $A$
for dealing with the aforementioned two cases. The third member is used for
the proof of Theorem (3.4) on PI-Split.

Case 1: (short interval) $Q_{min}\leq u$.

In this case, by using (74), we obtain

$\Delta_{k}\leq u\quad\mbox{if}\quad Q_{min}\leq u,$ (79)

namely, the duration of the interval is at most $u$. Since by (78)
$Q_{max}>s_{max}u$ and $s_{max}$ is the maximum number of jobs that can be
served during a single time slot, the largest queue length at $\tau_{k}$ (in
server $i_{max}$) does not reach zero at $\tau_{k+1}$. By (1) the queue length
in server $i_{max}$ must decrease by at least one. Thus, using (61), we
conclude that $R_{i_{max}}\geq 1$ and $Q_{max}\wedge R_{i_{max}}\geq 1$.
Therefore, we have

$\sum_{i\neq l_{k}}^{n}Q_{i}^{\gamma}(Q_{i}\wedge R_{i})\geq
Q_{max}^{\gamma}(Q_{max}\wedge R_{i_{max}})\geq Q_{max}^{\gamma}.$ (80)

Using (79) and (80) to upper bound the right hand side of (73) yields

$\Phi_{k}\leq\mathbb{E}\Big{[}\epsilon
u+([\lambda-\mu_{l_{k}}]^{+})^{1+\gamma}u^{1+\gamma}+Cu^{0.75(1+\gamma)}-Q_{max}^{\gamma}\mid{\cal{F}}^{\xi}_{\tau_{k}}\Big{]}<0,$

where the last inequality is due to (78).

Case 2: (long interval) $Q_{min}>u$.

By the definition of the sampling times, $R_{i}$ can be larger than $Q_{i}$ by
at most $s_{i}(\tau_{k+1})$ (the number of jobs potentially served by server
$i$ at the last time slot of the interval $[\tau_{k}+1,\tau_{k+1}]$), which we
assumed in (1) is bounded from above by $s_{max}$. Thus

$Q_{i}\wedge R_{i}\geq(R_{i}-s_{max})\wedge R_{i}=R_{i}-s_{max}.$

Therefore

$\mathbb{E}\Big{[}\sum_{i\neq l_{k}}^{n}Q_{i}^{\gamma}(Q_{i}\wedge
R_{i})\mid{\cal{F}}^{\xi}_{\tau_{k}}\Big{]}\geq\mathbb{E}\Big{[}\sum_{i\neq
l_{k}}^{n}Q_{i}^{\gamma}(R_{i}-s_{max})\mid{\cal{F}}^{\xi}_{\tau_{k}}\Big{]}=\sum_{i\neq
l_{k}}^{n}Q_{i}^{\gamma}(\mathbb{E}[R_{i}\mid{\cal{F}}^{\xi}_{\tau_{k}}]-s_{max}),$
(81)

where the last equality is due to the fact that $Q_{i}$ is measurable with
respect to ${\cal{F}}^{\xi}_{\tau_{k}}$ for all $i$.

Next, we calculate $\mathbb{E}[R_{i}\mid{\cal{F}}^{\xi}_{\tau_{k}}]$. By the
definitions of $R_{i}$ in (61) and $\Delta_{k}$ in (65), this is the expected
value of a random walk with strictly positive step sizes, distributed as
$s_{i}(0)$, starting at time $\tau_{k}+1$ and stopped after $\Delta_{k}$ time
slots. If $\Delta_{k}$ only depended on the queue length at server $i$, then
it would have been a stopping time, and by Wald’s Equation and the finitness
of the sampling times the expected value would have been equal to
$\mu_{i}\mathbb{E}[\Delta_{k}\mid{\cal{F}}^{\xi}_{\tau_{k}}]$. However, this
is not the case here, since $\Delta_{k}$ depends on all servers different from
$l_{k}$. Even so, we are able to prove the same result holds. Using the
definition of $R_{i}$ in (61) and the fact that
$\tau_{k+1}=\tau_{k}+\Delta_{k}$ we obtain

$\displaystyle\mathbb{E}[R_{i}\mid{\cal{F}}^{\xi}_{\tau_{k}}]$
$\displaystyle=\mathbb{E}\Bigg{[}\sum_{t=\tau_{k}+1}^{\tau_{k+1}}s_{i}(t)\mid{\cal{F}}^{\xi}_{\tau_{k}}\Bigg{]}=\mathbb{E}\Bigg{[}\sum_{t=\tau_{k}+1}^{\tau_{k}+\Delta_{k}}s_{i}(t)\mid{\cal{F}}^{\xi}_{\tau_{k}}\Bigg{]}$
(82)
$\displaystyle=\mathbb{E}\Bigg{[}\sum_{t=\tau_{k}+1}^{\tau_{k}+Q_{min}}s_{i}(t)\mathbbm{1}_{\\{t\leq\tau_{k}+\Delta_{k}\\}}\mid{\cal{F}}^{\xi}_{\tau_{k}}\Bigg{]}$
(83)
$\displaystyle=\mathbb{E}\Bigg{[}\sum_{m=1}^{Q_{min}}s_{i}(\tau_{k}+m)\mathbbm{1}_{\\{m\leq\Delta_{k}\\}}\mid{\cal{F}}^{\xi}_{\tau_{k}}\Bigg{]},$
(84)

where the second to last equality is due to the fact that by (74), we have
$\Delta_{k}\leq Q_{min}$, and the last equality is a simple variable change.

Since the random variables $Q_{min}$ and $\tau_{k}$ are both measurable with
respect to ${\cal{F}}^{\xi}_{\tau_{k}}$ and are finite a.s., we can change the
order of the summation and the expectation, yielding

$\displaystyle\mathbb{E}\Bigg{[}\sum_{m=1}^{Q_{min}}s_{i}(\tau_{k}+m)\mathbbm{1}_{\\{m\leq\Delta_{k}\\}}\mid{\cal{F}}^{\xi}_{\tau_{k}}\Bigg{]}=\sum_{m=1}^{Q_{min}}\mathbb{E}[s_{i}(\tau_{k}+m)\mathbbm{1}_{\\{m\leq\Delta_{k}\\}}\mid{\cal{F}}^{\xi}_{\tau_{k}}].$
(85)

The event $\\{m\leq\Delta_{k}\\}$ is equivalent to the event
$\\{\Delta_{k}\leq m-1\\}^{c}$. Thus, given ${\cal{F}}^{\xi}_{\tau_{k}}$,
whether or not the event $\\{m\leq\Delta_{k}\\}$ occurred is determined by the
values of all the service processes, in all servers, during the interval
$[\tau_{k}+1,\ldots,\tau_{k}+m-1]$, which are independent of
$s_{i}(\tau_{k}+m)$. Therefore

$\displaystyle\sum_{m=1}^{Q_{min}}\mathbb{E}[s_{i}(\tau_{k}+m)\mathbbm{1}_{\\{m\leq\Delta_{k}\\}}\mid{\cal{F}}^{\xi}_{\tau_{k}}]$
$\displaystyle=\sum_{m=1}^{Q_{min}}\mathbb{E}[s_{i}(\tau_{k}+m)\mid{\cal{F}}^{\xi}_{\tau_{k}}]\mathbb{E}[\mathbbm{1}_{\\{m\leq\Delta_{k}\\}}\mid{\cal{F}}^{\xi}_{\tau_{k}}]$
(86)
$\displaystyle=\sum_{m=1}^{Q_{min}}\mu_{i}\mathbb{E}[\mathbbm{1}_{\\{m\leq\Delta_{k}\\}}\mid{\cal{F}}^{\xi}_{\tau_{k}}]=\mu_{i}\mathbb{E}[\Delta_{k}\mid{\cal{F}}^{\xi}_{\tau_{k}}],$
(87)

where the second equality is due to $s_{i}(\tau_{k}+m)$ being independent of
${\cal{F}}^{\xi}_{\tau_{k}}$ and the last equality is due to the fact that
$\sum_{m=1}^{Q_{min}}\mathbbm{1}_{\\{m\leq\Delta_{k}\\}}=\Delta_{k}$.

Combining (81), (82), (85) and (86) yields

$\mathbb{E}\Big{[}\sum_{i\neq l_{k}}^{n}Q_{i}^{\gamma}(Q_{i}\wedge
R_{i})\mid{\cal{F}}^{\xi}_{\tau_{k}}\Big{]}\geq\sum_{i\neq
l_{k}}^{n}Q_{i}^{\gamma}(\mu_{i}\mathbb{E}[\Delta_{k}\mid{\cal{F}}^{\xi}_{\tau_{k}}]-s_{max})=\sum_{i\neq
l_{k}}^{n}\mathbb{E}[Q_{i}^{\gamma}(\mu_{i}\Delta_{k}-s_{max})\mid{\cal{F}}^{\xi}_{\tau_{k}}].$
(88)

Since we consider the case where $Q_{min}>u$, and by (76) $u>s^{2}_{max}$, by
(74) we obtain

$\Delta_{k}>Q_{min}/s_{max}>u/s_{max}>s_{max}.$ (89)

Thus, using the fact that by (2), $\mu_{i}\geq 1$ for all $i$, we obtain

$\mu_{i}\Delta_{k}-s_{max}>0.$ (90)

Using (88), (90) and the fact that $\Delta_{k}\leq Q_{min}\leq Q_{i}$ for
$i\neq l_{k}$, we obtain

$\displaystyle\mathbb{E}\Big{[}\sum_{i\neq
l_{k}}^{n}Q_{i}^{\gamma}(Q_{i}\wedge
R_{i})\mid{\cal{F}}^{\xi}_{\tau_{k}}\Big{]}$ $\displaystyle\geq\sum_{i\neq
l_{k}}^{n}\mathbb{E}[\Delta_{k}^{\gamma}(\mu_{i}\Delta_{k}-s_{max})\mid{\cal{F}}^{\xi}_{\tau_{k}}]=\sum_{i\neq
l_{k}}^{n}\mathbb{E}[\mu_{i}\Delta_{k}^{1+\gamma}-s_{max}\Delta_{k}^{\gamma}\mid{\cal{F}}^{\xi}_{\tau_{k}}]$
(91) $\displaystyle=\mathbb{E}\Big{[}\Big{(}\sum_{i\neq
l_{k}}^{n}\mu_{i}\Big{)}\Delta_{k}^{1+\gamma}\mid{\cal{F}}^{\xi}_{\tau_{k}}\Big{]}-(n-1)s_{max}\mathbb{E}[\Delta_{k}^{\gamma}\mid{\cal{F}}^{\xi}_{\tau_{k}}].$
(92)

Using (91) to bound the right hand side of (73) from above, we obtain

$\displaystyle\Phi_{k}$
$\displaystyle\leq\mathbb{E}\Big{[}\epsilon\Delta_{k}+\Big{(}([\lambda-\mu_{l_{k}}]^{+})^{1+\gamma}-\sum_{i\neq
l_{k}}^{n}\mu_{i}\Big{)}\Delta_{k}^{1+\gamma}+C\Delta_{k}^{0.75(1+\gamma)}+(n-1)s_{max}\Delta_{k}^{\gamma}\mid{\cal{F}}^{\xi}_{\tau_{k}}\Big{]}$
$\displaystyle\leq\mathbb{E}\Big{[}\epsilon\Delta_{k}-2\epsilon\Delta_{k}^{1+\gamma}+C\Delta_{k}^{0.75(1+\gamma)}+(n-1)s_{max}\Delta_{k}^{\gamma}\mid{\cal{F}}^{\xi}_{\tau_{k}}\Big{]}=\mathbb{E}[f(\Delta_{k})]\mid{\cal{F}}^{\xi}_{\tau_{k}}],$

where we have used the definitions of $\gamma$ in (54) and $f(\cdot)$ in (75).
Finally, by (89) $\Delta_{k}>u/s_{max}$, so by (76) we have $f(\Delta_{k})<0$
yielding $\Phi_{k}\leq 0$. This concludes the proof.

∎

###  3.4 Stability of PI-Split

######  Theorem 3.4.

Assume $\lambda<\sum_{i=1}^{n}\mu_{i}$. Then:
(i) $X^{S}$ is positive recurrent. Consequently,
(ii) $X^{S}$ has a unique stationary distribution, denoted by $\pi^{S}_{X}$,
and
(iii) For any initial state $\alpha\in\mathcal{A}$ and any
$B\subset\mathcal{A}$, $\mathbb{P}_{\alpha}(X^{S}(t)\in
B)\rightarrow\pi^{S}_{X}(B)$ as $t\rightarrow\infty$.

Proof. The proof is similar to that of Theorem (3.2). Specifically, we prove
that the state dependent drift in Lemma (3.3) is also fulfilled by the sampled
chain $Y^{S}$, with the same finite set $A$ (52), constants $\epsilon$ (53),
$\gamma$ (54) and the Lyapunov function ${\cal{L}}$ (56), but with stopping
times and a filtration which correspond to the chain $X^{S}$. Theorem (3.4)
then follows from the same arguments made in the proof of Theorem (3.2).

Recall that ${\cal{A}}^{*}$ denotes the collection of states in which there is
at least one queue length equal to zero that does not correspond to the Last
Idle server. For initial states $\alpha\in{\cal{A}}^{*}\setminus A$ in which
there is only one such queue length, the behavior of PI and PI-Split is
identical since batches are not split until (at least) the next sampling time.
Thus, the analysis in the proof of Lemma (3.3) proves that the desired state
dependent drift holds for these initial states.

It remains to consider initial states $\alpha\in{\cal{A}}^{*}\setminus A$ such
that there are at least two queue lengths equal to zero which do not
correspond to the Last-Idle server. By the definition of the sampling times,
the next sampling time must be after exactly 1 time slot. Thus, it remains to
prove that

$\mathbbm{E}_{\alpha}\Big{[}{\cal{L}}(X^{S}(1))-{\cal{L}}(X^{S}(0))\Big{]}=\sum_{i=1}^{n}\mathbbm{E}_{\alpha}\Big{[}(Q_{i}^{S}(1))^{1+\gamma}-(Q_{i}^{S}(0))^{1+\gamma}\Big{]}\leq-\epsilon.$
(93)

Decompose the sum after the first equality into three parts:

Part 1: the single member of the sum corresponding to the server with the
maximum queue length at time slot 0 (if there are several, choose one
arbitrarily). Denote the server with the maximal queue length by
$j^{*}=\mbox{argmax}_{1\leq i\leq n}Q^{S}_{i}(0)$. By the definition of PI-
Split, server $j^{*}$ does not receive jobs at the beginning of time slot $1$
because its queue length is positive and there are other serves which are
idle. By the definition of the finite set in (52) and (77),
$Q^{S}_{j^{*}}(0)\geq s_{max}$. Hence, the queue length at server $j*$ cannot
reach zero after one time slot and we have
$Q^{S}_{j^{*}}(1)=Q^{S}_{j^{*}}(0)-s_{j^{*}}(1)$.

Part 2: the members of the sum corresponding to servers with queue length
equal to zero at time slot 0. The queue lengths of these servers at time slot
1 is upper bounded by the number of jobs that arrive at the beginning of time
slot 1, namely, $a(1)$. Thus, since $\gamma\leq 1$, $a^{2}(1)$ is an upper
bound for the contribution of each member of part 2 to the sum in (93).

Part three: the members of the sum not included in parts 1 and 2. By the
definition of PI-Split, members of part 3 do not receive any work. Thus, their
contribution to the sum in (93) is non positive. Specifically, it is upper
bounded by $a^{2}(1)$.

We conclude that

$\displaystyle\mathbbm{E}_{\alpha}\Big{[}{\cal{L}}(X^{S}(1))-{\cal{L}}(X^{S}(0))\Big{]}$
$\displaystyle\leq\mathbbm{E}_{\alpha}\Big{[}(Q_{j^{*}}^{S}(1))^{1+\gamma}-(Q_{j^{*}}^{S}(0))^{1+\gamma}\Big{]}+(n-1)\mathbbm{E}_{\alpha}[a^{2}(1)]$
(94)
$\displaystyle=\mathbbm{E}_{\alpha}\Big{[}(Q_{j^{*}}^{S}(0)-s_{j^{*}}(1))^{1+\gamma}-(Q_{j^{*}}^{S}(0))^{1+\gamma}\Big{]}+n(\sigma_{a}^{2}+\lambda^{2})$
(95)
$\displaystyle\leq(Q_{j^{*}}^{S}(0)-1)^{1+\gamma}-(Q_{j^{*}}^{S}(0))^{1+\gamma}+n(\sigma_{a}^{2}+\lambda^{2})$
(96)
$\displaystyle\leq-(Q_{j^{*}}^{S}(0))^{\gamma}+n(\sigma_{a}^{2}+\lambda^{2}),$
(97)

where in the second inequality we used the fact that by (1), $s_{j^{*}}(1)\geq
1$ and that given the initial state $\alpha$, $j^{*}$ and $Q_{j^{*}}^{S}(0)$
are known. The last inequality follows from the same calculation as is done in
(69). Finally, by the choice of the finite set in (77) (specifically, the
third member in the right hand side of (77)), we obtain

$\mathbbm{E}_{\alpha}\Big{[}{\cal{L}}(X^{S}(1))-{\cal{L}}(X^{S}(0))\Big{]}\leq-\epsilon,$

which concludes the proof. ∎

##  4 Simulations

10 slow : 90 fast 50 slow : 50 fast 90 slow : 10 fast

Figure 1: Performance of different policies with 100 heterogeneous servers and
different slow:fast server ratios. The policies we consider here cannot split
arriving batches of jobs. We emphasize this by adding ’u’ before each policy
name (which stands for ’unsplittable’).

10 slow : 90 fast 50 slow : 50 fast 90 slow : 10 fast

Figure 2: Performance of different policies with 100 heterogeneous servers and
different slow:fast server ratios. The policies we consider here can split
arriving batches of jobs. We emphasize this by adding ’s’ before each policy
name (which stands for ’splittable’).

The main purpose of this section is to complement the simulation results
presented in [1] for a model with batch arrivals and random service
capacities. To this end, we compare the performance of PI and PI-Split with
that of the token-based JIQ as well as several other load-distribution
policies, namely JSQ, JSQ($2$) and JSQ($1$,$1$).

We consider systems with 100 heterogeneous servers and study two cases
corresponding to whether arriving batches can be split or not. When batches
cannot be split, the aforementioned policies route the entire batch to a
single queue. When batches can be split, it is not straight-forward how the
policies should split them. There are many sensible ways in which this can be
done. In our simulations we chose the following splitting schemes, based on
the “water-filling” approach (see, for example, [2]):

  * •

JSQ: Uses the water-filling approach where the jobs belonging to the batch are
sent to the servers with the smallest queue lengths in an attempt to equalize
them as much as possible. For example, if there are 3 servers with queue
lengths $(1,2,4)$ and a batch arrives with 3 jobs, server 1 will get 2 jobs
and server 2 will get 1 job, resulting in the queue lengths $(3,3,4)$. Ties
are broken randomly.

  * •

JSQ(2): Uses the same water-filling approach as JSQ, but only for the two
servers that were chosen uniformly at random.

  * •

JSQ(1,1): Same as JSQ(2), but for the “remembered” server and the one chosen
uniformly at random.

  * •

JIQ: Uses water-filling only for the idle servers. Specifically, a batch is
split as evenly as possible between the idle servers. If there are no idle
servers, each job in the batch is sent to a randomly chosen server.

  * •

PI: If there are idle servers, the splitting is the same as in JIQ. Otherwise,
the entire batch is sent to the last idle server.

As shall be seen throughout this section, simulation results indicate that PI
is indeed stable and offers appealing properties in comparison to JIQ and
other reduced-state policies.

###  4.1 Simulation settings

We consider a system with $100$ heterogeneous servers, some being slow
servers, working at a rate $\mu^{-1}=1$, and the others being fast servers,
working at a rate $(10\mu)^{-1}=1/10$. The service time distribution of all
servers is uniform and lower bounded by 1. Note that these service time
distributions concur with Assumption (1)). We implement three scenarios, which
correspond to different partitions of the servers into the two classes. More
specifically, we consider the following slow:fast ratios: $10{:}90$, $50{:}50$
and $90{:}10$. The number of arrivals at each time slot is drawn from a
Poisson process of rate $\lambda$. We are interested in values of $\lambda$
such that $0<\lambda<\sum_{i=1}^{n}\mu_{i}$, _i.e.,_ any strictly positive and
sub-critical arrival rates.

For each scenario, we run simulations for $10^{6}$ time slots for different
loads, where all policies receive the same input process, but possibly make
different routing decisions. The number of time slots was chosen such that the
difference in the outputs of different runs at the maximal load was
negligible. We present the following graphs:

Steady state average queue length vs. load. For each simulation run with a
given load, we calculate the time averaged sum of queue lengths in the system.
Under all policies, whenever the chain is positive recurrent, it is also
ergodic. Thus the average queue length converges, and one simulation run is
enough to calculate the steady-state average queue length. Also, there are
points that do not appear on the graphs since they are vertically truncated.
However, since the truncation is at very high average queue lengths, these
missing points indicate either poor performance or instability.

CDF of JCT in steady state for 50$\%$ and 99$\%$ load. We keep track of the
job completion times of all jobs that entered the system during the
simulations, with $\lambda{=}0.5\sum_{i=1}^{n}\mu_{i}$ and
$\lambda{=}0.99\sum_{i=1}^{n}\mu_{i}$. We calculate the cumulative
distribution function (CDF). Policies that are unstable under a given load are
not shown in the graph.

Messages per time slot vs. load. For JIQ and PI, for each simulation run with
a given load, we calculate the number of times a token is passed from a server
to the dispatcher (the reverse path is not counted because the token is sent
along with a job). We then divide this by the number of time slots to obtain
the average number of messages per time slot. Counting messages for the other
policies is more straightforward. At each time slot with non-zero arrivals,
for JSQ it is the number of servers needed to be probed, _i.e.,_ $100$.
Similarly, for JSQ(2) and JSQ(1,1) it is $2$.

######  Remark 4.1.

Instead of calculating messages per time slot, one can choose to look at
messages per batch, per sub batches (in the splittable case) or per job.
However, measuring the number of messages per time slot is more informative
for the amount of needed communication, whereas the other options are affected
by the arrival and departure rates rather than by their ratio (_i.e.,_
normalized load).

###  4.2 The unsplittable case: simulation results

Stability and average queue lengths. The results are depicted in Figure 1. As
can be seen, for all the different slow:fast ratios, PI appears to achieve the
stability region. Also, the average queue length under PI is favorable to that
of other reduced-state policies and even comparable to that of JSQ for low and
moderate loads. JSQ(1,1) also appears to achieve the stability region, but
with a much larger average queue length. JSQ(2) appears to be unstable for
moderately-high loads for the $50{:}50$ and $90{:}10$ ratios. Also, whenever
JSQ(2) is stable, its average queue length is much higher than under PI. JIQ
appears unstable for loads close to $100\%$ for all ratios.

Job completion time. We measure the completion time of all jobs during the
simulation at a moderate load of 50$\%$ and a high load of 99$\%$. It is
evident how both PI and JIQ are competitive to JSQ at 50% for all slow:fast
ratios. This is because at such load there are many idle servers and both of
these token-based policies resemble JSQ. JSQ(2) and JSQ(1,1) show lesser
performance, and at a 50:50 ratio JSQ(2) is missing as its completion times
diverge due to apparent instability. At a load of 99$\%$, PI shows the best
overall performance after JSQ. Both JIQ and JSQ(2) do not appear in 2 out of 3
scenarios due to what appears to be instability at this load. JSQ(1,1) again
shows lesser performance.

Messages per time slot. Since the probability of zero arrivals in a time slot
is negligible even at low loads (_i.e.,_ $e^{-\lambda}$), the message overhead
of JSQ, JSQ(2) and JSQ(1,1) is approximately constant. It is 100 for JSQ and 2
for JSQ(2) and JSQ(1,1). On the other hand, the number of messages per time
slot for PI and JIQ is 1 up to a high load. As the load increases, there is an
increasing number of time slots without any servers becoming idle and the
message overhead for both policies decreases.

###  4.3 Splittable case: simulation results

Stability and average queue length. The results are depicted in Figure 2. As
can be seen, the results show a similar trend to those of the unsplittable
case. It is also evident how the advantage of JSQ in comparison to the other
policies increases.

Job completion time and Messages per time slot. In the splittable case, the
communication overhead of token based policies has different characteristics
than in the unsplittable case. This is because in the splittable case, the
overhead may reach 1 message per job rather than 1 message per batch of jobs.
Indeed, this would be the case for moderate loads at which we have both a
sufficient number of jobs in each batch and they are split among several/many
idle servers.

As can be seen from the graphs, at such moderate loads, both PI ans JIQ are
competitive with JSQ but with often lower message overhead. JSQ(2) and
JSQ(1,1) perform significantly worse but also have a lower message overhead in
this case. At the high load of 99$\%$ the message overhead of the token based
policies rapidly decreases due to the decreased number of idled servers. In
this case the results are similar to those of the unspittable case with
apparent instability of JSQ(2) in 2 out of 3 cases and JIQ in 1 out of 3.
JSQ(1,1) appears to be stable but again performs poorly.

######  Remark 4.2.

We have found in our simulations that the stability region of JIQ increases
when batches are allowed to be split. To the best of our knowledge, this
phenomenon was not pointed out previously. We remark it as a future research
challenge to provide formal analysis for this phenomenon and explore how it
can be leveraged to improve the performance of token based policies.

##  5 Appendix

Proof of Equation (68). To reiterate (68), we prove that

$\mathbb{E}\Big{[}\Big{(}{Q_{l_{k}}(\tau_{k+1})}\Big{)}^{1+\gamma}\mid{\cal{F}}^{\xi}_{\tau_{k}}\Big{]}\leq([\lambda-\mu_{l_{k}}]^{+})^{1+\gamma}\mathbb{E}[\Delta_{k}^{1+\gamma}\mid{\cal{F}}^{\xi}_{\tau_{k}}]+C\mathbb{E}[\Delta_{k}^{0.75(1+\gamma)}\mid{\cal{F}}^{\xi}_{\tau_{k}}],$
(98)

where $l_{k}\in[n]$ is the Last-Idle server during the interval
$[\tau_{k}+1,\tau_{k+1}]$, $\Delta_{k}=\tau_{k+1}-\tau_{k}$ is the duration of
the interval and $Q_{l_{k}}(\tau_{k+1})$ is the queue length at server $l_{k}$
at time $\tau_{k+1}$. Recall that the RVs $l_{k}$, $\tau_{k}$ and
$\\{Q_{i}(\tau_{k})\\}_{i\in[n]}$ are measurable w.r.t
${\cal{F}}^{\xi}_{\tau_{k}}$, and that by the definition of the sampling
times, we must have $Q(\tau_{k})=0$ and that server $l_{k}$ receives all
incoming jobs during $[\tau_{k}+1,\tau_{k+1}]$.

The left hand side of (98) is therefore the $1+\gamma$ moment of the position,
after $\Delta_{k}$ time slots, of a reflected random walk starting at zero.
Recall the definition in (66) of $\beta_{i}:=\lambda-\mu_{i}$ and
$\sigma_{i}^{2}:=\sigma_{a}^{2}+\sigma_{s_{i}}^{2}$. Given the value of
$l_{k}$, the mean and variance of the step size of this random walk are given
by $\beta_{l_{k}}$ and $\sigma_{l_{k}}^{2}$ respectively.

Lemma 6 of [1] provides an upper bound on the $1+\gamma$ moment of a very
similar reflected random walk. There are two differences between that result
and what we want to prove here. First, the random walk in [1] is of the
workload process, not the queue length process as we consider here.
Specifically, the distribution of the step size has a different dependence on
the model parameters. However, the proof in [1] only uses the mean and
variance of the step size. Therefore, the same derivation is applicable here
and we just plug in our mean and variance.

Second, the bound in [1] is derived for a fixed number of steps. In our case,
given ${\cal{F}}^{\xi}_{\tau_{k}}$, the number of steps is random. This is due
to the random service capacities. But, the number of steps until the next
sampling time is independent of the random walk at server $l_{k}$. It only
depends on the other servers and how many time slots it takes one of their
queue lengths to reach zero.

Thus, a direct application of Lemma 6 in [1] (after plugging in
$\beta_{l_{k}}$ and $\sigma_{l_{k}}$) yields

$\mathbb{E}\Big{[}\Big{(}{Q_{l_{k}}(\tau_{k+1})}\Big{)}^{1+\gamma}\mid{\cal{F}}^{\xi}_{\tau_{k}},\Delta_{k}=t\Big{]}\leq([\beta_{l_{k}}]^{+})^{1+\gamma}t^{1+\gamma}+Ct^{0.75(1+\gamma)},$

where

$C=\max_{i}\\{\big{(}16\sigma_{i}+8(\sigma_{i})^{1/2}(\beta_{i})^{+}\big{)}^{(1+\gamma)/2}\\}.$

Thus, we have

$\displaystyle\mathbb{E}\Big{[}\Big{(}{Q_{l_{k}}(\tau_{k+1})}\Big{)}^{1+\gamma}\mid{\cal{F}}^{\xi}_{\tau_{k}}\Big{]}$
$\displaystyle=\sum_{t=1}^{\infty}\mathbb{E}\Big{[}\Big{(}{Q_{l_{k}}(\tau_{k+1})}\Big{)}^{1+\gamma}\mid{\cal{F}}^{\xi}_{\tau_{k}},\Delta_{k}=t\Big{]}\mathbbm{P}(\Delta_{k}=t\mid{\cal{F}}^{\xi}_{\tau_{k}})$
$\displaystyle\leq\sum_{t=1}^{\infty}\Big{(}([\beta_{l_{k}}]^{+})^{1+\gamma}t^{1+\gamma}+Ct^{0.75(1+\gamma)}\Big{)}\mathbbm{P}(\Delta_{k}=t\mid{\cal{F}}^{\xi}_{\tau_{k}})$
$\displaystyle=([\lambda-\mu_{l_{k}}]^{+})^{1+\gamma}\mathbb{E}[\Delta_{k}^{1+\gamma}\mid{\cal{F}}^{\xi}_{\tau_{k}}]+C\mathbb{E}[\Delta_{k}^{0.75(1+\gamma)}\mid{\cal{F}}^{\xi}_{\tau_{k}}].$

∎

## References

  * [1] R. Atar, I. Keslassy, G. Mendelson, A. Orda, S. Vargaftik, Persistent-idle load-distribution, Stochastic Systems 10 (2) (2020) 152–169.
  * [2] L. Ying, R. Srikant, X. Kang, The power of slightly more than one sample in randomized load balancing, Mathematics of Operations Research 42 (3) (2017) 692–722.
  * [3] Y. Lu, Q. Xie, G. Kliot, A. Geller, J. R. Larus, A. Greenberg, Join-idle-queue: A novel load balancing algorithm for dynamically scalable web services, Performance Evaluation 68 (11) (2011) 1056–1071.
  * [4] M. Mitzenmacher, The power of two choices in randomized load balancing, IEEE Transactions on Parallel and Distributed Systems 12 (10) (2001) 1094–1104.
  * [5] D. Shah, B. Prabhakar, The use of memory in randomized load balancing, in: Proceedings IEEE International Symposium on Information Theory, 2002, p. 125.
  * [6] S. Asmussen, Applied probability and queues, Vol. 51, Springer Science & Business Media, 2008.
  * [7] V. A. Malyshev, M. V. Men’shikov, Ergodicity, continuity and analyticity of countable markov chains, Trudy Moskovskogo Matematicheskogo Obshchestva 39 (1979) 3–48.
  * [8] S. P. Meyn, R. Tweedie, State-dependent criteria for convergence of markov chains, The Annals of Applied Probability (1994) 149–168.
  * [9] R. Durrett, Probability: theory and examples, Vol. 49, Cambridge university press, 2019.