# Stabilizing Queuing Networks with Model Data-Independent Control

Qian Xie and Li Jin Q. Xie is with the Tandon School of Engineering, New York
University, USA. L. Jin is with the UM Joint Institute and the School of
Electronic Information and Electrical Engineering, Shanghai Jiao Tong
University, China and with the Tandon School of Engineering, New York
University, USA (emails<EMAIL_ADDRESS>li.jin@sjtu.edu.cn).

###### Abstract

Classical queuing network control strategies typically rely on accurate
knowledge of model data, i.e. arrival and service rates. However, such data
are not always available and may be time-variant. To address this challenge,
we consider a class of model data-independent (MDI) control policies that only
rely on traffic state observation and network topology. Specifically, we focus
on the MDI control policies that can stabilize multi-class Markovian queuing
networks under centralized and decentralized policies. Control actions include
routing, sequencing, and holding. By expanding the routes and constructing
piecewise-linear test functions, we derive an easy-to-use criterion to check
the stability of a multi-class network under a given MDI control policy. For
stabilizable multi-class networks, we show that a centralized, stabilizing MDI
control policy exists. For stabilizable single-class networks, we further show
that a decentralized, stabilizing MDI control policy exists. In addition, for
both settings, we construct explicit policies that attain maximal throughput
and present numerical examples to illustrate the results.

Keywords: Multi-class queuing networks, Dynamic routing, Lyapunov function,
Stability

##  1 Introduction

Control on multi-class queuing networks has been studied in numerous contexts
of transportation, logistics, and communication systems [1, 2, 3, 4]. Most
existing analysis and design approaches rely on full knowledge of model data,
i.e., arrival and service rates, to ensure stability and/or optimality [5].
However, in many practical settings, such data may be unavailable or hard to
estimate, and may be varying over time. Such challenges motivate the idea of
_model data-independent (MDI)_ control policies. MDI control policies select
control actions, including routing, sequencing, and/or holding, according to
state observation and network topology but independent of arrival/service
rates. Such policies are easy to implement and, if appropriately designed, can
resist modeling error or non-stationary environment. However, the stability of
general open multi-class queuing networks with centralized or decentralized
MDI control policies has not been well studied.

In this paper, we consider the stability of multi-class queuing networks with
throughput-maximizing MDI control policies. Particularly, we focus on acyclic
open queuing networks with Poisson arrivals and exponential service times.
Jobs (customers) are classified according to their origin-destination (OD)
information. Service rates are independent of job classes. A network is
stabilizable if there exists a control policy that ensures positive Harris
recurrence of the queuing process, whether the network is open-loop or closed-
loop, centralized or decentralized [6]. By standard results on Jackson
networks, stabilizability is equivalent to the existence of a (typically model
data-dependent) stabilizing Bernoulli routing policy [7]. We assume that the
class-specific arrival rates and the server-specific service rates are unknown
to the controller. The main results are as follows:

  1. 1.

An easy-to-use criterion to check the stability of a multi-class network under
a given MDI control policy (Proposition 1).

  2. 2.

For a multi-class network, a stabilizing centralized MDI control policy exists
if and only if the network is stabilizable (Theorem 1).

  3. 3.

For a single-class network, a stabilizing decentralized MDI control policy
exists if and only if the network is stabilizable (Theorem 2).

Previous works on stability of queuing networks are typically based on full
knowledge of model data [8, 1, 9, 10, 11, 12, 13]. So far, the best-studied
MDI control policy is the join-the-shortest-queue (JSQ) routing policy for
parallel queues [14, 15, 16, 17, 18, 19, 20, 21] or simple networks [22],
which requires only the queue lengths and does not rely on model data [23].
When and only when the network is stabilizable, i.e., the demand is less than
capacity (service rate), the JSQ policy guarantees the stability of parallel
queues/simple networks [24, 22] and the optimality of homogeneous servers
[15]. However, JSQ routing does not guarantee stability of more complex
networks [23]. MDI routing for general networks has been numerically evaluated
[25], but no structural results are available. Besides, most studies on MDI
routing for general networks are not aimed for stability [26, 27, 28, 29]. In
addition, decentralized dynamic routing has been considered for single origin-
destination networks [30, 31] but not in MDI settings.

To design stabilizing MDI control policies, we first develop a stability
criterion (Proposition 1) based on route expansion for queuing networks and
explicit construction of a piecewise-linear test function. The expanded
network is essentially a parallel connection of all routes from the set of
origins to the set of destinations. With this expansion, we use insights on
the behavior of parallel queues and of tandem queues to construct the test
function and derive the stability criterion. The test function can be used to
obtain a smooth Lyapunov function verifying a negative drift condition. The
piecewise-linear test function technique was proposed by Down and Meyn [6];
however, their implementation relies on linear programming formulations to
determine parameters of the test function, which depends on model data. We
will extend this technique to the MDI setting using explicitly constructed
test functions.

Based on the stability criterion, we design control policies in centralized
and decentralized settings. First, for multi-class networks, we present a
stabilizing centralized MDI control policy requiring dynamic routing and
preemptive sequencing named JSR policy (Theorem 1). The control policy is
obtained by minimizing the mean drift of the piecewise-linear test function,
and the mean drift is guaranteed to be negative if and only if the network is
stabilizable. The JSR policy, which is centralized and MDI, maximizes
throughput among all control policies. Compared with other centralized
policies, it does not require knowledge of model data, and compared with other
MDI policies (e.g., JSQ), it guarantees stability for any stabilizable
networks. Second, for single-class networks, we present a decentralized
routing and holding policy that guarantees stability (Theorem 2). Such
policies can also maximize the throughput since the stabilizability of the
network implies that the throughput can be as large as close to the capacity.
The results are closely related to the theory on the classical JSQ routing
policy [23] and the decentralized max-pressure control policy [32].

The rest of this paper is organized as follows. Section 2 defines the multi-
class queuing network model. Section 3 presents the stability criterion based
on route expansion and piecewise-linear test function. Section 4 and Section 5
consider the control design problem in centralized and decentralized settings
respectively. Section 6 gives concluding remarks.

##  2 Multi-class queuing network

Consider an acyclic network of queuing servers with infinite buffer spaces.
Let $\mathcal{N}$ be the set of _servers_. Each server $n$ has an exponential
_service rate_ $\bar{\mu}_{n}$. The network has a set $\mathcal{S}$ of
_origins_ and a set $\mathcal{T}$ of _destinations_. jobs are classified
according to their origins and destinations. That is, we can use an origin-
destination (OD) pair $(S,T)\in\mathcal{C}$ to denote a _job class_ , or
simply _class_. For notational convenience, classes (OD pairs) are indexed by
$c=(S_{c},T_{c})$. jobs of class $c$ arrive at $S_{c}$ according to a Poisson
process of rate $\lambda_{c}\geq 0$. We assume that service rates are
independent of job class.

The _topology_ of the network is characterized by _routes_ between origins and
destinations. We use $|r|$ to denote the number of servers on route $r$. Let
$\mathcal{R}_{c}$ be the set of routes between $S_{c}$ and $T_{c}$, and define
$\mathcal{R}=\bigcup_{c\in\mathcal{C}}\mathcal{R}_{c}$. Below is an example
network to illustrate the notations.

###### Example 1

Consider the Wheatstone bridge network in Fig. 1.

Figure 1: A two-class queuing network.

Two classes of jobs arrive at $S_{1}$ (resp. $S_{2}$) with $\lambda_{1}>0$
(resp. $\lambda_{2}>0$). The set of servers is
$\mathcal{N}=\\{1,2,\ldots,5\\}$ and the set of OD-specific routes are

$\displaystyle\mathcal{R}_{1}=\\{(1,3),(4)\\},\
\mathcal{R}_{2}=\\{(2),(3,5)\\}.$

The _state_ of the network is defined as follows. Let
$\bar{x}=[\bar{x}_{n}^{c}]_{n\in\mathcal{N},c\in\mathcal{C}}$ be the vector of
class-specific _job numbers_ , where $\bar{x}_{n}^{c}$ is the number of jobs
of class $c$ in server $n$, either waiting or being served. Let
$\bar{\mathcal{X}}$ be the space of $\bar{x}$. We use $\bar{X}(t)$ to denote
the state of the queuing process at time $t$.

We consider three types of _control actions_ , viz. routing, sequencing, and
holding. All control actions are essentially Markovian (in terms of $\bar{x}$
plus additional auxiliary states) and are applied at the instant of
_transitions_ , which include the arrival of a job at an origin or the
completion of service at a server. _Routing_ refers to allocating an incoming
job to a server downstream to the origin or allocating a job discharged by a
server to another downstream server. _Sequencing_ refers to selecting a job
from the waiting queue to serve. The default sequencing policy is the first-
come-first-serve (FCFS) policy. For the multi-class setting, we consider the
preemptive-priority that can terminate an ongoing service and start serving
jobs from another class, while the job with incomplete service is sent back to
the queue. _Holding_ refers to holding a job that has completed its service in
the server while blocking the other jobs in the queue from accessing the
server.

Following [33], we say that a queuing network is _stable_ if the queuing
process is positive Harris recurrent. For details about the notion of positive
Harris recurrence for queuing networks, see [8, 33, 6]. Finally, we say that
the network is _stabilizable_ if a stabilizing control exists. One can check
the stabilizability using the following result:

###### Lemma 1

An open acyclic queuing network is stabilizable if and only if there exists a
vector $[\xi_{r}]_{r\in\mathcal{R}}$ such that

$\displaystyle\xi_{r}\geq 0,\quad\forall r\in\mathcal{R},$
$\displaystyle\lambda_{c}=\sum_{r\in\mathcal{R}_{c}}\xi_{r},\quad\forall
c\in\mathcal{C},$ $\displaystyle\sum_{r\in\mathcal{R}:n\in
r}\xi_{r}<\bar{\mu}_{n},\quad\forall n\in\mathcal{N}.$

The proof and implementation are straightforward.

##  3 Stability criterion

In this section, we derive a stability criterion for multi-class networks
under given control policies. The techniques that we use include the route
expansion of the original network and the explicit construction of a
piecewise-linear test function based on the network topology. In Section 3.1,
we construct an expanded network based on the original network. In Section
3.2, we apply a piecewise-linear test function to the expanded network to
obtain a stability criterion (Proposition 1) for both the expanded and the
original networks.

###  3.1 Route expansion

For the convenience of constructing test function, we first introduce the
route expansion. _Route expansion_ refers to the construction of an _expanded
network_ based on the topology of _original network_ (defined in Section 2).
The high-level idea is to decompose the network into routes, and the specific
procedures are:

  1. 1.

Place all routes $\mathcal{R}$ in the original network in parallel.

  2. 2.

Add two-way connections between duplicates of servers in the original network.

For example, Fig. 2 shows the expanded network constructed from the original
network in Fig. 1.

Figure 2: Route expansion of the network in Fig. 1.

We call ”servers” in the expanded network as _subservers_ , since they are
obtained by duplicating actual servers in the original network. Subservers are
indexed by $k$, $c_{k}\in\mathcal{C}$ is the class index,
$r_{k}\in\mathcal{R}$ is the route index, and
$i_{k}\in\\{1,2,\ldots,|r_{k}|\\}$ is the numbering of subserver $k$ on route
$r_{k}$. We use $k\in r$ to refer to that subserver $k$ is on route $r$. Let
$\mathcal{K}$ be the set of all subservers and $\mathcal{K}_{c}$ be the set of
subservers with $c_{k}=c$. We use $n_{k}\in\mathcal{N}$ to denote the actual
server that corresponds to subserver $k$. In addition, let $k_{p}$ (resp.
$k_{s}$) denote the subserver immediately upstream (resp. downstream) to
subserver $k$.

The _state_ of the expanded network is
$x=\\{x^{c}_{k};k\in\mathcal{K},c\in\mathcal{C}\\}$, denoting the vector of
number of class-$c$ jobs in subserver $k$. Let
$x_{k}:=\sum_{c\in\mathcal{C}}x_{k}^{c},k\in\mathcal{K}.$ The expanded state
space is $\mathcal{X}=\mathbb{Z}_{\geq 0}^{|\mathcal{C}|\times|\mathcal{K}|}$.
Note that the states of the expanded network and the states of the original
network are related by

$\displaystyle\bar{x}^{c}_{n}=\sum_{k\in\mathcal{K}:n_{k}=n}x^{c}_{k},\quad
k\in\mathcal{K},$ (1)

for each $n\in\mathcal{N}$.

The routing policy is characterized by
$\pi:\mathcal{X}\to[0,1]^{|\mathcal{C}|\times|\mathcal{K}|^{2}}$, where
$\pi^{c}_{k,k^{\prime}}$ is the probability that a class-$c$ job is routed
from subserver $k$ to subserver $k^{\prime}$.

The holding policy is characterized by
$\zeta:\mathcal{X}\to\\{0,1\\}^{|\mathcal{K}|}$, where $\zeta_{k}$ specifies
whether subserver $k$ is holding ($\zeta_{k}(x)=0$) or not holding
($\zeta_{k}(x)=1$) when the current state is $x$.

Two subservers $k$ and $k^{\prime}$ are _duplicating_ if
$n_{k}=n_{k^{\prime}}$. Note that the service rates of duplicating subservers
are coupled in the sense that for each server $n\in\mathcal{N}$, at a given
time, at most one subserver $k$ such that $n_{k}=n$ can be actively serving
jobs, or _active_. This can be modeled as an _imaginary service rate control
policy_ $\mu:\mathcal{X}\to\mathbb{R}^{|\mathcal{K}|}$ such that the service
rate $\mu_{k}(x)$ of subserver $k$ satisfies

$\displaystyle\sum_{k:n_{k}=n}\mu_{k}(x)\leq\bar{\mu}_{n},\quad\forall
x\in\mathcal{X}.$

Such control policy is essentially equivalent to the class-based preemptive
sequencing in the original network.

Note that $\\{X(t):t\geq 0\\}$ is a Markov process, and the positive Harris
recurrence refers to that there exists a unique invariant measure $\nu$ on
$\mathcal{X}$ such that for every measurable set $D\subseteq\mathcal{X}$ with
$\nu(D)>0$ and for every initial condition $x\in\mathcal{X}$,

$\Pr\\{\tau_{D}<\infty|X(0)=x\\}=1,$

where $\tau_{D}=\inf\\{t\geq 0:X(t)\in D\\}$. Also, though
$\\{\bar{X}(t):t\geq 0\\}$ is not a Markov process, it will eventually
converge to a steady state distribution.

The route expansion technique not only expands the network but also decomposes
the state variables. Jobs can move along the expanded network using two
transition mechanisms. One is _actual transition_ , referring to moving a job
from subserver $k$ (or an origin) to its downstream subserver $k_{s}$ (or a
destination). The other is _imaginary transition_ that moves a job from one
subserver $k$ to a duplicating subserver $k^{\prime}$ thereof, see imaginary
switch in Section 5. Imaginary transitions always occur instantaneously. Note
that an actual transition corresponds to a transition in the original network,
while an imaginary transition does not; this is also revealed in (1).

One can always map a control action in the expanded network to the original
network. However, an MDI control policy may not exist on the state space of
the original network; we do need an expanded state space for MDI control. In
addition, we allow imaginary control actions in the expanded network,
including _imaginary service rate control_ and _imaginary switch_ ; see
Section 4 and Section 5. Such imaginary actions only make sense in the
expanded network and do not correspond to actual service rate control or
switch in the original network.

###  3.2 Stability of the expanded network

After introducing the expanded network and the mathematical definition of the
control policy, the main question of this paper can be expressed in a formal
way as follows.

Given an expanded queuing network and a control policy $\phi=(\pi,\mu,\zeta)$,
how can we tell if the control policy $\phi$ is stabilizing, i.e., the
expanded network is stable under $\phi$?

The answer of this question will be given in Proposition 1. Before that, we
need to introduce the test function technique first. As opposed to linear
programming-based construction in [6], we provide an explicit construction,
where parameters of the test function do not rely on solving any optimization
problems. The high level idea is to identify the bottlenecks and their
upstream subservers. Our construction is based on the route expansion
described in the previous subsection.

  1. 1.

For each class $c\in\mathcal{C}$ and expanded state $x\in\mathcal{X}$, define

$\displaystyle
g_{c}(x):=\max_{\begin{subarray}{c}K_{c}\subseteq\mathcal{K}_{c}:\\\ \kappa\in
K_{c}\Rightarrow\kappa_{p}\in\mathcal{K}_{c}\end{subarray}}\sum_{k\in
K_{c}}a_{k}x_{k},$

where $a_{k}\in(0,1)$ is a parameter.

  2. 2.

Define a piecewise-linear _test function_

$\displaystyle V(x):=\max_{C\subseteq\mathcal{C}}\sum\limits_{c\in
C}b_{c}g_{c}(x),$

where $b_{c}\in(0,1)$ is a parameter.

We call $V(x)$ the test function rather than the Lyapunov function, since
strictly speaking, a smooth Lyapunov function should be developed based on the
piecewise-linear test function to verify the Foster-Lyapunov stability
criterion. Down and Meyn [6] showed that as long as a piecewise-linear test
function can be determined, one can always smooth it to obtain a qualified
$C^{2}$ Lyapunov function.

###### Remark 1

The test functions we proposed in this work are MDI. But generally speaking,
they do not need to be MDI since it does not affect the control policies to be
MDI.

######  Definition 1 (Dominance)

Consider state $x\in\mathcal{X}$.

  1. 1.

We call $C^{*}$ a set of _dominant_ classes if

$C^{*}\in\operatorname*{argmax}_{C\subseteq\mathcal{C}}\sum_{c\in
C}b_{c}g_{c}(x).$

Each class $c\in C^{*}$ is a _dominant_ class.

  2. 2.

We call $K^{*}_{c}$ a set of _dominant_ class-c subservers if

$K^{*}_{c}\in\operatorname*{argmax}\limits_{\begin{subarray}{c}K_{c}\subseteq\mathcal{K}_{c}:\\\
k\in K_{c}\Rightarrow k_{p}\in\mathcal{K}_{c}\end{subarray}}\sum_{k\in
K_{c}}a_{k}x_{k}.$

Each subserver $k\in K^{*}_{c}$ is a _dominant_ class-c subserver.

  3. 3.

A route $r\in\mathcal{R}_{c}$ is _dominant_ if it includes dominant class-c
subservers, i.e. there exists dominant class-c subserver $k\in K^{*}_{c}$ such
that $k\in r$.

Let $R_{c}$ be the set of dominant class-c routes.

  4. 4.

A subserver $b\in K^{*}_{c}$ is called a _bottleneck_ if it is a dominant
class-c subserver while its immediate downstream subserver $b_{s}\notin
K^{*}_{c}$ is not.

###### Remark 2

A route or server is dominant if changes in its traffic state immediately
affect the test function $V$.

A _regime_ $X$ of the piecewise-linear test function is a subset of
$\mathcal{X}$ such that there exist $C^{X}\subseteq\mathcal{C}$,
$K^{X}=\bigcup_{c\in C^{X}}K_{c}^{X}\subseteq\mathcal{K}$, and
$R^{X}=\bigcup_{c\in C^{X}}R_{c}^{X}\subseteq\mathcal{R}$ where
$C^{X},K_{c}^{X},R_{c}^{X}$ are dominant for each $x\in X$, i.e., the test
function is linear over $X$. Let $\mathscr{X}$ be the set of regimes; note
that $\bigcup_{X\in\mathscr{X}}X=\mathcal{X}$.

######  Definition 2 (Mean velocity and drift)

Consider a multi-class network with state $x\in\mathcal{X}$ under an expanded
control policy $\phi=(\pi,\mu,\zeta)$.

  1. 1.

The _mean velocity_ at state $x$ is a function
$v:\mathcal{X}\to\mathbb{R}^{|\mathcal{K}|}$ such that for each
$k\in\mathcal{K}$,

$\displaystyle
v_{k}(x):=\sum_{c\in\mathcal{C}}\lambda_{c}\pi_{S_{c},k}^{c}(x)+\mu_{k_{p}}(x)\zeta_{k_{p}}(x)-\mu_{k}(x)\zeta_{k}(x).$

where $\pi^{c}_{S_{c},k}$ is the probability that a class-$c$ job is routed
from origin $S_{c}$ to subserver $k$, while $\mu_{k}$ and $\zeta_{k}$ are the
controlled service rate and the holding status of the subserver $k$
respectively.

  2. 2.

Given $X\in\mathscr{X}$ such that $x\in X$, the _mean drift over $X$_ is given
by

$\displaystyle D^{X}(x):=\sum_{c\in C^{X}}b_{c}\sum_{k\in
K^{X}_{c}}a_{k}v_{k}(x).$

###### Remark 3

In our subsequent analysis, the mean drift $D^{X}(x)$ of the test function
will play the role of infinitesimal generator applied to a Lyapunov function;
see [6] for the connection between the test function and the Lyapunov
function.

The main result of this section is as follows:

###### Proposition 1

Consider a multi-class network under the expanded control policy $\phi$.
Suppose there exist constants $M<\infty$, $\epsilon>0$, and
$a_{k},b_{c}\in(0,1)$ ($\forall$ $c\in\mathcal{C}$, $k\in K_{c}$), such that
for each $X\in\mathscr{X}$ and each $x\in X$ where
$|x|=\sum\limits_{k\in\mathcal{K}}x_{k}>M$,

$\displaystyle\sum_{c\in C^{X}}b_{c}\sum_{k\in
K^{X}_{c}}a_{k}v_{k}(x)\leq-\epsilon.$ (2)

Then, the network is stable.

_Proof._ Consider the test function $V(x)$. By its definition, if $x\in X$, we
have

$V(x)=\sum_{c\in C^{X}}b_{c}\sum_{k\in K^{X}_{c}}a_{k}x_{k}.$

The mean drift is given by

$\displaystyle D^{X}(x)$ $\displaystyle=\sum_{c\in C^{X}}b_{c}\sum_{k\in
K^{X}_{c}}a_{k}v_{k}(x)$
$\displaystyle\stackrel{{\scriptstyle\footnotesize\eqref{eq_betaalpha}}}{{\leq}}-\gamma^{|C^{X}|-1}\epsilon\leq-\gamma^{|\mathcal{C}|-1}\epsilon,\quad
x:|x|>M.$

One can then apply [6, Theorem 1] and [6, Lemma 5] to obtain the stability of
the network. $\square$

As a benchmark, the approach in [6, Theorem 1] requires solving linear
programs to obtain parameters of the test functions in Proposition 1, while
our approach explicitly constructs the parameters (see Section 4 and Section
5). Moreover, the proposed control, which is independent of model data,
guarantees stability if and only if the network is stabilizable (see Theorem 1
and Theorem 2), while the approach in [6] relies on knowledge of model data.

##  4 Centralized control for multiple classes

In this section, we consider the “join-the-shortest-route (JSR)” policy (a
joint routing and sequencing policy) for centralized control. The JSR policy
is MDI and constructed based on the expanded network. We will show that it is
stabilizing if and only if the network is stabilizable.

The test functions are constructed as follows.

  1. 1.

For each class $c\in\mathcal{C}$, each route $r\in\mathcal{R}_{c}$, and each
expanded state $x\in\mathcal{X}$, let

$\displaystyle f_{r}(x):=\max_{k\in r}\alpha^{i_{k}-1}\sum_{j:i_{j}\leq
i_{k}}x_{j},$ $\displaystyle
g_{c}(x):=\max_{R_{c}\subseteq\mathcal{R}_{c}}\beta^{|R_{c}|-1}\sum_{r\in{R_{c}}}f_{r}(x),$

where $\alpha\in(0,1),\beta\in(0,1)$ are constant parameters.

  2. 2.

The piecewise-linear _test function_ is given by

$\displaystyle V(x):=\max_{C\subseteq\mathcal{C}}\gamma^{|C|-1}\sum_{c\in
C}g_{c}(x),$

where $\gamma\in(0,1)$ is a constant parameter.

Let the parameters be such that

$\displaystyle\alpha=\beta\geq\frac{|\mathcal{R}|-1}{|\mathcal{R}|},\quad\gamma\geq\frac{|\mathcal{C}|-1}{|\mathcal{C}|},$
(3)

and follow the notions of dominance accordingly (see Definition 1). Note that
such MDI parameters $\alpha,\beta,\gamma$ always exist. The control that we
consider in this subsection only depends on $\alpha,\beta,\gamma$ and is thus
MDI. Specifically, we define the JSR policy as follows:

######  Definition 3 (Join-the-shortest-route (JSR) policy)

  1. 1.

(Routing) At an origin $S$, an incoming job of class $c$ is allocated to the
route $r^{*}\in\mathcal{R}_{c}$ such that

$\displaystyle r^{*}\in\operatorname*{argmin}_{r\in\mathcal{R}_{c}}f_{r}(x).$

If there is only one minimum, then $r^{*}$ must be a non-dominant route.
Otherwise, let $b^{*}$ be the bottleneck on route $r^{*}$. Then, an incoming
job of class $c$ is allocated to the route $r^{*}\in\mathcal{R}_{c}$ with the
largest $i_{b^{*}}$, which is denoted by $i_{c}$. Further ties are randomly
broken.

  2. 2.

(Imaginary service rate control) Let $\mathcal{K}_{n}$ be the set of
subservers corresponding to server $n$ and let $\mathcal{B}$ be the set of
bottlenecks for a given $x$. Then, a subserver $k\in\mathcal{K}_{n}$ is
activated if $k\in\mathcal{B}$. If multiple subservers are in
$\mathcal{K}_{n}\cap\mathcal{B}$, then activate the subserver $k^{*}$ such
that

$\displaystyle
k^{*}=\operatorname*{argmin}_{k\in\mathcal{K}_{n}\cap\mathcal{B}}\\{i_{c_{k}}+|\mathcal{R}_{c_{k}}|\\};$

ties are randomly broken. This is to ensure that the bottlenecks are active to
discharge jobs and only one of the duplicating subservers can be active.

The main result of this section is the following:

######  Theorem 1 (Stability of JSR policy)

The JSR policy stabilizes a multi-class network if and only if the network is
stabilizable.

This theorem implies that the JSR policy is also throughput-maximizing, as
long as the network is stabilizable, i.e., the demand is less than the total
capacity. Note that the stabilizability can be easily checked using Lemma 1.

In the rest of this section, we apply Theorem 1 to study the stability of the
Wheatstone bridge network under the JSR policy (Subsection 4.1) and then prove
this theorem (Subsection 4.2).

###  4.1 Numerical Example

Consider the network in Fig. 1 and suppose that
$\lambda_{1}=\lambda_{2}=\lambda=1$ and $\bar{\mu}_{n}=\mu=1$ for $n=1,2,4,5$
and $\bar{\mu}_{3}=\frac{1}{4}$. This example is for illustrating the route
expansion and the test function construction.

Note that under the above model parameters, the decentralized JSQ policy is
destabilizing. To see this, $\bar{\mu}_{1}=\bar{\mu}_{4}$ implies that on
average, class-1 jobs are evenly distributed between server 1 and server 4.
Thus, the average departure rate of class-1 jobs from server 1 is
$\frac{1}{2}$, which exceeds the service rate of server 3. Therefore, the
queue at server 3 is unstable. The main reason that the JSQ policy is
destabilizing is the ignorance of downstream congestion. As $\bar{X}_{3}(t)$
gets large, a reasonable action is to allocate fewer class-1 jobs to server 1.
However, the JSQ policy disallows such far-sighted decisions.

An alternative centralized stabilizing routing policy can be the following JSR
policy:

  1. 1.

A class-1 job arriving at $S_{1}$ is routed to server 1 if
$\bar{X}_{1}^{1}(t)+\bar{X}_{3}^{1}(t)<\bar{X}_{4}^{1}(t)$, to server 4 if
$\bar{X}_{1}^{1}(t)+\bar{X}_{3}^{1}(t)>\bar{X}_{4}^{1}(t)$, and randomly
otherwise.

  2. 2.

A class-2 job arriving at $S_{2}$ is routed to server 3 if
$\bar{X}_{3}^{2}(t)+\bar{X}_{5}^{2}(t)<\bar{X}_{2}^{2}(t)$, to server 2 if
$\bar{X}_{3}^{2}(t)+\bar{X}_{5}^{2}(t)>\bar{X}_{2}^{2}(t)$, and randomly
otherwise.

  3. 3.

The dominant class has a higher priority.

That is, when jobs are routed at $S_{1}$, the decision is based on not only
the local state ($\bar{X}_{1}(t)$ and $\bar{X}_{4}(t)$), but also the state
further downstream ($\bar{X}_{3}(t)$).

The expanded network is shown in Fig. 2. Each block in the figure represents a
subserver. In particular, subservers 3a and 3b are decomposed from server 3;
the other servers are remained. Solid arrows correspond to actual transitions
in an original network, while dashed arrows correspond to imaginary
transitions between duplicating subservers.

In the expanded network, a job can move along both solid and dashed arrows.
The color of an arrow shows which class can move along it: blue means class
$(S_{1},T_{1})$, red means $(S_{2},T_{2})$, and purple means both. For ease of
presentation, we label $(S_{1},T_{1})$ as class 1 and $(S_{2},T_{2})$ as class
2. For example, a job of class $(S_{1},T_{1})$ can visit subservers 4, 1, 3a,
3b and the destination $T_{1}$.

In the expanded network, the JSR policy works as follows.

  1. 1.

A class-1 job arriving at $S_{1}$ is routed to subserver 4 if
$X_{4}(t)<X_{1}(t)+X_{3a}(t)$, to subserver 1 if
$X_{4}(t)>X_{1}(t)+X_{3a}(t)$, and randomly otherwise.

  2. 2.

A class-2 job arriving at $S_{2}$ is routed to subserver 3b if
$X_{3b}(t)+X_{5}(t)<X_{2}(t)$, to subserver 2 if
$X_{3b}(t)+X_{5}(t)>X_{2}(t)$, and randomly otherwise.

  3. 3.

If subserver 3a is dominant while subserver 3b is non-dominant, and server 3
is serving a class-2 job, then server 3 preempts the class-2 job being served
in 3b to the class-1 job in 3a, and vice versa. If both subserver 3a and
subserver 3b are dominant, then server 3 gives priority to class-2 job since
the index of 3b is smaller.

By Theorem 1, the network can be stabilized by the JSR policy if and only if

$\displaystyle\lambda_{1}<2,\ \lambda_{2}<2,\
\lambda_{1}+\lambda_{2}<\frac{9}{4}.$

We use the following parameters for the test function:

$\alpha=\beta=\gamma=\frac{3}{4},\ \epsilon=\Big{(}\frac{3}{4}\Big{)}^{5}.$

One can verify that the above parameters satisfy (3) and Proposition 1 by
considering the following cases:

  1. 1.

Only one route is dominant. In this case, an incoming job is always allocated
to a non-dominant route, leading to non-positive contribution to the mean
drift:

$\displaystyle
D^{X}(x)\leq-\gamma\beta\alpha\mu=-\Big{(}\frac{3}{4}\Big{)}^{3}\leq-\epsilon.$

  2. 2.

Two routes with different OD pairs are dominant. This case is analogous to the
previous case:

$\displaystyle
D^{X}(x)\leq-\gamma\beta\alpha\mu=-\Big{(}\frac{3}{4}\Big{)}^{3}\leq-\epsilon.$

  3. 3.

Two routes with the same OD pair or more than two routes are dominant. In such
cases, the mean drift satisfies

$\displaystyle
D^{X}(x)\leq\gamma\beta^{3}(\lambda-\mu-\alpha\mu)=-\Big{(}\frac{3}{4}\Big{)}^{5}\leq-\epsilon.$

Consequently, the network is stable under the MDI JSR policy.

###  4.2 Proof of Theorem 1

In this subsection, we will the sufficiency and the necessity respectively,
based on the connection between the sign of the mean drift and the
stabilizability of the network. When analyzing the mean drift, we consider two
parts: external arrivals and internal transmission. We first show that any
internal transmission does not positively contribute to the mean drift and
then show that any positive contribution from external arrivals can always be
compensated by internal transmissions.

####  4.2.1 Internal transmissions

Note that under the JSR policy, every job remains on the route assigned to the
job when it enters the network. Hence, internal transmissions only occur
between subservers on the same route.

Given $x$, consider an internal transmission from subserver $k$ to subserver
$j$; this implicitly requires $x_{k}\geq 1$. The definition of dominance
ensures that if $j$ is dominant, then so is $k$. Hence, we need to consider
the following cases:

  1. 1.

If $k$ and $j$ are both dominant, the transmission leads to zero contribution
to the mean drift $D^{X}(x)$ for all $X$ such that $x\in X$.

  2. 2.

If $k$ is dominant and $j$ is non-dominant, the transmission leads to the
following contribution to the mean drift:

$\displaystyle-\alpha^{i_{k}-1}\mu_{k}(x)\leq 0.$

Hence, internal transmissions never lead to positive contribution to the mean
drift.

####  4.2.2 External arrivals

Given $x\neq 0$, consider a regime $X\in\mathscr{X}$ such that $x\in X$. For
each $c\in\mathcal{C}$, the JSR policy ensures that if there exists a non-
dominant route in $\mathcal{R}_{c}$, then an incoming job must be allocated to
a non-dominant route in $\mathcal{R}_{c}$, leading to non-positive
contribution to the mean drift. Hence, we only need to consider dominant
classes $c$ such that every route in $\mathcal{R}_{c}$ is dominant, i.e.
$R^{X}_{c}=\mathcal{R}_{c}$. Recall that $C^{*}\subseteq\mathcal{C}$ is the
set of dominant classes. The part of the mean drift associated with $c\in
C^{*}$ satisfies

$\displaystyle D^{X}_{c}(x)$
$\displaystyle\leq\gamma^{|C^{*}|-1}\beta^{|\mathcal{R}_{c}|-1}\Big{(}\alpha^{i_{c}-1}\lambda_{c}-\sum_{b\in\mathcal{B}^{X}:c_{b}=c}\alpha^{i_{b}-1}\mu_{b}(x)\Big{)}$
$\displaystyle:=\gamma^{|C^{*}|-1}\Delta_{c}^{X}(x)$

over any regimes of the piecewise-linear test function, where $i_{c}$ is given
in Definition 3.

###### Lemma 2

When $x\neq 0$, there is no empty bottleneck, i.e.

$\displaystyle x_{b}\geq 1.$ (4)

_Proof._

Since $x\neq 0$ and $r_{b}$ is dominant, we have

$\sum_{k:i_{k}\leq i_{b}}x_{k}>0.$

If $i_{b}=1$, then the above inequality directly implies (4).

Now consider the case that $i_{b}\geq 2$. Since $b$ is a bottleneck, we have

$\displaystyle\alpha^{i_{b}-1}\sum_{k\in r_{b}:i_{k}\leq
i_{b}}x_{k}\geq\alpha^{i_{b}-2}\sum_{k\in r_{b}:i_{k}\leq i_{b}-1}x_{k},$

which implies

$\displaystyle x_{b}\geq(1-\alpha)\sum_{k:i_{k}\leq i_{b}}x_{k}>0$

and thus we have (4). $\square$

Lemma 2 is to ensure that the bottlenecks are none-empty to discharge jobs and
thus contribute negative terms to the drift.

Next, we show the sufficiency of Theorem 1. Based on the definition of the
routing policy (see Definition 3), $\forall b\in\mathcal{B}^{X}$, we have
$i_{b}\leq i_{c}$ when the incoming job is allocated to a dominant route. Then

$\displaystyle\sum_{c\in{C^{*}}}\Delta_{c}^{X}(x)$
$\displaystyle\leq\sum_{c\in{C^{*}}}\beta^{|\mathcal{R}_{c}|-1}\alpha^{i_{c}-1}\Big{(}\lambda_{c}-\sum_{b\in\mathcal{B}^{X}:c_{b}=c}\mu_{b}(x)\Big{)}$
$\displaystyle=\sum_{c\in{C^{*}}}\alpha^{i_{c}+|\mathcal{R}_{c}|-2}\Big{(}\lambda_{c}-\sum_{b\in\mathcal{B}^{X}:c_{b}=c}\mu_{b}(x)\Big{)}$
$\displaystyle:=\sum_{c\in{C^{*}}}\alpha_{c}\beta_{c}$

Without loss of generality assume $C^{*}=\\{1,2,\cdots,m\\}$ and

$i_{1}+|\mathcal{R}_{1}|\leq i_{2}+|\mathcal{R}_{2}|\leq\cdots\leq
i_{m}+|\mathcal{R}_{m}|.$

Then by using Abel transformation (summation by parts), the right hand side of
the above inequality (abbr. RHS):

$\displaystyle
RHS=\sum_{i=1}^{m-1}(\alpha_{i}-\alpha_{i+1})\sum_{j=1}^{i}\beta_{j}+\alpha_{m}\sum_{j=1}^{m}\beta_{m}.$

Based on the assumption, we have $\alpha_{i}\geq\alpha_{i+1}$ and

$\displaystyle\sum_{j=1}^{i}\beta_{j}$
$\displaystyle=\sum_{j=1}^{i}\Big{(}\lambda_{j}-\sum_{b\in\mathcal{B}^{X}:c_{b}=j}\mu_{b}(x)\Big{)}$
$\displaystyle=\sum_{j=1}^{i}\lambda_{j}-\sum_{n_{b}:b\in\mathcal{B}_{i}^{X}}\bar{\mu}_{n_{b}}$
$\displaystyle=\sum_{j=1}^{i}\lambda_{j}-\sum_{n\in\mathcal{N}_{i}}\bar{\mu}_{n}$
$\displaystyle<0,$

where $\mathcal{B}_{i}$ is the set of bottlenecks in the first $i$ classes and
$\mathcal{N}_{i}$ is the min-cut of the original network with the first $i$
classes. Here we use the definition of the imaginary service rate control (see
Definition 3) and Lemma 1.

Since $RHS<0$, we have $\sum_{c\in{C^{*}}}\Delta_{c}^{X}(x)<0$ and thus
$D_{c}^{X}(x)<0$. Then by noting that internal transmissions lead to non-
positive contributions to the mean drift, we have

$\displaystyle D^{X}(x)\leq\sum_{c\in{C^{*}}}D^{X}_{c}(x)<0,$

which implies stability. $\square$

Finally, the necessity is apparent: if a network is not stabilizable, then
there exists no MDI control that can stabilize the network.

##  5 Decentralized control for a single class

For a single-class network, we can drop the class index and use $x_{k}$ to
denote the number of jobs in subserver $k$. Note that such network has a
single origin and a single destination. Again we can do route expansion on
such network.

We consider a decentralized MDI control policy as follows.

######  Definition 4 (JSQ with artificial spillback)

The JSQ with artificial spillback (JSQ-AS) policy is as follows:

  1. 1.

(Routing) A discharged job is routed to the shortest downstream queue, with
ties randomly broken.

  2. 2.

(Holding) For each subserver $k$, any job which has finished the service will
be held if and only if $X_{s_{k}}(t)\geq X_{k}(t)$.

  3. 3.

(Imaginary switch) When a dominant subserver $k$ is inactive while its non-
dominant duplicate $k^{\prime}$ is active, and both are not in the holding
status, then the job in $k^{\prime}$ is moved to (and discharged from) $k$
after service, and then routed to the downstream of $k$ (i.e., $s_{k}$).

Note that under the holding policy, the process $\\{X(t);t\geq 0\\}$ admits an
invariant set $\mathcal{Q}\subseteq\mathcal{X}$ given by

$\displaystyle\mathcal{Q}:=\\{x\in\mathcal{X}:x_{s_{k}}\leq
x_{k},k\in\mathcal{K}\\}.$ (5)

Since we consider the long-time stability of the network, it suffices to
consider the states in an invariant set. The above result indicates that in
the invariant set $\mathcal{Q}$, the queue size of any subserver is upper-
bounded by the queue size of its immediate upstream subserver.

The JSQ-AS policy is decentralized in the sense that control actions on
subserver $k$ only depend on local traffic information: the number of jobs in
duplicate subservers $\\{x_{k^{\prime}}:n_{k}=n_{k^{\prime}}\\}$ and that in
immediate downstream subservers
$\\{x_{s_{k^{\prime}}}:n_{k}=n_{k^{\prime}}\\}$. A key characteristic of such
policies is that congestion information can propagate through the network via
the forced holding: if a subserver becomes congested (i.e. $x_{k}$ gets
large), the congestion will propagate to the upstream subservers in a
cascading manner (“artificial spillback”). Importantly, such artificial
spillback does not undermine throughput like the natural spillback caused by
the limited buffer size. The reason is that though congestion can propagate,
the queue size in any downstream subserver is not upper-bounded. Artificial
spillback is the main difference between the JSQ-AS policy and the classic JSQ
policies.

Note that though the JSQ-AS policy is constructed based on the expanded
network, its actions can always be converted to the ones in the original
network. Importantly, the decentralized control in the expanded network must
also be decentralized in the original network. Also note that the imaginary
switch has no impact on the original network or the test function.

The main result of this section is as follows:

######  Theorem 2 (Stability of JSQ-AS policy)

For the route expansion of a single-class network, the JSQ-AS policy is
stabilizing if and only if

$\displaystyle\lambda<\bar{\mu}^{mc},$ (6)

where $\bar{\mu}^{mc}$ is the min-cut service rate of the original network.

This theorem implies that JSQ-AS policy is also a throughput-maximizing policy
since we allow any throughput that satisfies (6).

In the rest of this section, we apply Theorem 2 to study the stability of the
Wheatstone bridge network under the JSQ-AS policy (Subsection 5.1) and then
prove this theorem (Subsection 5.2).

###  5.1 Numerical Example

Figure 3: A single-class queuing network and its expanded network.

Consider the original network with route expansion in Fig. 3. Again suppose
that $\lambda=1$, $\bar{\mu}_{n}=\frac{3}{4}$ for $n=1,2,\cdots,5$. Similarly,
with the above parameters, the JSQ policy is destabilizing since the queue at
server 5 is unstable. However, in the decentralized setting, the control
actions can only depend on the local state, say the routing decision at the
origin can be based on $\bar{X}_{1}(t)$ and $\bar{X}_{4}(t)$, but not
$\bar{X}_{5}(t)$. A remedy is to introduce the holding policy (artificial
spillback) to the JSQ policy so that the downstream congestion can be relieved
and the local state can somehow reflect the states further downstream.

In the expanded network, server 1 is decomposed into subserver 1a and 1b,
server 5 is decomposed into subserver 5a and 5b. The states in the original
network and those in the expanded network satisfy
$\bar{X}_{1}(t)=X_{1a}(t)+X_{1b}(t)$ and $\bar{X}_{5}(t)=X_{5a}(t)+X_{5b}(t)$.
The initial states of the expanded network can be not unique. Say the initial
queue size of server 1 is 2, then the initial queue sizes of subserver 1a and
1b can be 2, 0 or 1, 1 or 0, 2 respectively. Then the states are updated based
on the model and our JSQ-AS policy. For example, the routing decision at the
origin is based on $X_{1a}(t)$, $X_{1b}(t)$ and $X_{4}(t)$ rather than
$\bar{X}_{1}(t)$ and $\bar{X}_{4}(t)$; a job which has just finished the
service at server 3 will be held if $X_{5a}(t)\geq X_{3}(t)$, once released,
it will be routed to the shorter downstream queue by comparing $X_{5a}(t)$ and
$X_{5b}(t)$.

###  5.2 Proof of Theorem 2

This proof uses the connection between the stabilizability condition (6) and
the sign of the mean drift. But before showing the mean drift is negative, we
first present the explicit MDI piecewise-linear test function and several key
lemmas that can help us analyze the mean drift.

The piecewise-linear test function is constructed as follows:

$\displaystyle V(x):=\max\limits_{\begin{subarray}{c}K\subseteq\mathcal{K}:\\\
\kappa\in K\Rightarrow
p_{\kappa}\in\mathcal{K}\end{subarray}}\bigg{\\{}\frac{1+(|K|-1)\delta}{|K|}\sum_{k\in
K}x_{k}\bigg{\\}},$

where $\delta$ can be any small value such that $0<\delta<1$.

The following lemmas are useful in proving Theorem 2 where we consider the
regime $X\subseteq\mathcal{Q}$ containing $x$.

###### Lemma 3

A bottleneck can not be in the holding status.

_Proof._ Otherwise, the bottleneck must have at least one downstream
subserver. By (5), $x_{s_{k}}\geq x_{b}$. Since $b$ is a bottleneck, we have

$\displaystyle\frac{1+(|K^{X}|-2)\delta}{|K^{X}|-1}\sum_{k\neq
b}x_{k}\leq\frac{1+(|K^{X}|-1)\delta}{|K^{X}|}\sum_{k\in K^{X}}x_{k},$ (7)

which implies

$\displaystyle(1-\delta)\sum_{k\in
K^{X}}x_{k}\leq|K^{X}|[1+(|K^{X}|-2)\delta]x_{b}.$ (8)

Since $x_{b}\leq x_{s_{k}}$, we have

$\displaystyle(1-\delta)\sum_{k\in
K^{X}}x_{k}<|K^{X}|(1+|K^{X}|\delta)x_{s_{k}},$

which is equivalent to

$\displaystyle\frac{1+(|K^{X}|-1)\delta}{|K^{X}|}\sum_{k\in
K^{X}}x_{k}<\frac{1+|K^{X}|\delta}{|K^{X}|+1}\Big{(}\sum_{k\in
K^{X}}x_{k}+x_{s_{k}}\Big{)},$

contradicting with the fact that subserver $b$ is dominant and subserver
$s_{k}$ is non-dominant. $\square$

###### Corollary 1

Based on (8), we have $x_{b}>0$, i.e., any bottleneck $b$ must be non-empty.

This corollary and Lemma 3 ensure that all bottlenecks can discharge customers
and contribute negative terms to the drift.

###### Lemma 4

Let $k_{r}^{1}$ be the first subserver on route $r$, then either the route
with the smallest $x_{k_{r}^{1}}$ is non-dominant or every route is dominant.

_Proof._ If there is only one route, then that route must be dominant. Now
assume there are at least two routes and route $\hat{r}$ has the smallest
$x_{k_{r}^{1}}$, i.e. $\forall\ r\in\mathcal{R}$, $x_{k_{\hat{r}}^{1}}\leq
x_{k_{r}^{1}}$. Suppose $k_{\hat{r}}^{1}\in K^{X}$ and $\exists\
r\in\mathcal{R}$ s.t. $k_{r}^{1}\notin K^{X}$. Note that by (5), $x_{b}\leq
x_{k_{\hat{r}}^{1}}\leq x_{k_{r}^{1}}$, then from (7) we have

$\displaystyle\frac{1+(|K^{X}|-1)\delta}{|K^{X}|}\sum_{k\in
K^{X}}x_{k}<\frac{1+|K^{X}|\delta}{|K^{X}|+1}\Big{(}\sum_{k\in
K^{X}}x_{k}+x_{k_{r}^{1}}\Big{)},$

contradicting with our supposition. Therefore, either $\hat{r}$ is non-
dominant or every route in $\mathcal{R}$ is dominant. $\square$

###### Lemma 5

If $x\in\mathcal{X}$ makes every route $r\in\mathcal{R}$ dominant, then we
have

$\displaystyle\sum_{k\in\mathcal{B}^{X}}\mu_{k}(x)=\sum_{n:n=n_{k},k\in\mathcal{B}^{X}}\bar{\mu}_{n}$

_Proof._ Once there is an inactive bottleneck $k$ and a non-dominant but
active duplicate subserver $k^{\prime}$, the imaginary switch mechanism will
move the job being served in $k^{\prime}$ to $k$ and move one job in $k$ to
$k^{\prime}$. This is allowed since both $k$ and $k^{\prime}$ contain at least
one job due to the fact that a job being served in $k^{\prime}$ and the
bottleneck $k$ must be non-empty. $\square$

Similar to the proof of Theorem 1, we first analyze the internal transmissions
and then the external arrivals.

####  5.2.1 Internal transmissions

In the proof of Theorem 1, we have already discussed the case where internal
transmissions between subservers are on the same route. However, unlike the
JSR policy, the JSQ-AS policy allows internal transmissions between subservers
on different routes. Hence, we also need to consider the internal transmission
from subserver $k$ to subserver $j$ where $r_{k}\neq r_{j}$.

The definition of dominance ensures that if $k$ is non-dominant, so is
$s_{k}$. According to the routing policy, $x_{s_{k}}\geq x_{j}$. Let $\ell$ be
the first non-dominant subserver on route $r_{k}$ and $b$ be the bottleneck on
route $r_{j}$. If $j$ is non-dominant, then by (5), we have $x_{\ell}\geq
x_{s_{k}}\geq x_{j}\geq x_{b}$. Now from (7) we can obtain

$\displaystyle\frac{1+(|K^{X}|-1)\delta}{|K^{X}|}\sum_{k\in
K^{X}}x_{k}<\frac{1+|K^{X}|\delta}{|K^{X}|+1}\Big{(}\sum_{k\in
K^{X}}x_{k}+x_{\ell}\Big{)},$

contradicting with the definition of dominant subservers.

Thus, it cannot be the case that $k$ is non-dominant and $j$ is dominant,
which implies that any internal transmission does not positively contribute to
the mean drift.

####  5.2.2 External arrivals

According to Lemma 4, if a non-dominant route exists, then the routing policy
guarantees that an arriving job must be routed to the first subserver on a
non-dominant route $r$; this leads to non-positive contribution to the mean
drift. Otherwise, every route is dominant. Then for any $x\in\mathcal{Q}$
($x\neq 0$), the drift satisfies

$\displaystyle D^{X}(x)$ $\displaystyle\overset{\text{Corollary
}\ref{cor_bottleneck}}{\leq}\frac{1+(|K^{X}|-1)\delta}{|K^{X}|}\Big{(}\lambda-\sum_{b\in\mathcal{B}^{X}}\mu_{b}(x)\zeta_{b}(x)\Big{)}$
$\displaystyle\overset{\eqref{invariant}}{=}\frac{1+(|K^{X}|-1)\delta}{|K^{X}|}\Big{(}\lambda-\sum_{b\in\mathcal{B}^{X}}\mu_{b}(x)\Big{)}$
$\displaystyle\overset{\text{Lemma
}\ref{lmm_cut}}{=}\frac{1+(|K^{X}|-1)\delta}{|K^{X}|}\Big{(}\lambda-\sum_{n:n=n_{k},k\in\mathcal{B}^{X}}\bar{\mu}_{n}\Big{)}$
$\displaystyle\overset{\text{Lemma }\ref{lmm_stabilizability}}{<}0,$

which completes the proof. $\square$

The JSQ-AS policy cannot be directly applied to multi-class network, because
the imaginary switch mechanism may move a job to the subserver of a different
class with a different destination. Although the imaginary service rate
control in the JSR policy can be used for multiple classes, it needs global
information such as the information of dominance and bottlenecks for the
preemption, so it is not suitable for the decentralized setting. The design of
a decentralized MDI control policy for multi-class network can be a future
work.

##  6 Concluding remarks

We study the stability of open queuing networks under a class of model data-
independent control policies. In addition, we derive an easy-to-use stability
criterion based on route expansion of the network and explicit piecewise-
linear test functions. With the stability criterion, we generalize the
classical join-the-shortest-queue policy to ensure stability and attain
maximum throughput under centralized/decentralized settings. Our analysis and
design can also be applied to specific network control problems with stability
issues.

## References

  * [1] P. Kumar and S. P. Meyn, “Stability of queueing networks and scheduling policies,” _IEEE Transactions on Automatic Control_ , vol. 40, no. 2, pp. 251–260, 1995.
  * [2] S. P. Meyn, “Sequencing and routing in multiclass queueing networks part i: Feedback regulation,” _SIAM Journal on Control and Optimization_ , vol. 40, no. 3, pp. 741–776, 2001.
  * [3] S. L. Smith, M. Pavone, F. Bullo, and E. Frazzoli, “Dynamic vehicle routing with priority classes of stochastic demands,” _SIAM Journal on Control and Optimization_ , vol. 48, no. 5, pp. 3224–3245, 2010.
  * [4] R. Zhang, F. Rossi, and M. Pavone, “Analysis, control, and evaluation of mobility-on-demand systems: a queueing-theoretical approach,” _IEEE Transactions on Control of Network Systems_ , vol. 6, no. 1, pp. 115–126, 2018\.
  * [5] D. Bertsimas, I. C. Paschalidis, and J. N. Tsitsiklis, “Optimization of multiclass queueing networks: Polyhedral and nonlinear characterizations of achievable performance,” _The Annals of Applied Probability_ , pp. 43–75, 1994.
  * [6] D. Down and S. P. Meyn, “Piecewise linear test functions for stability and instability of queueing networks,” _Queueing Systems_ , vol. 27, no. 3-4, pp. 205–226, 1997.
  * [7] R. G. Gallager, _Stochastic processes: theory for applications_. Cambridge University Press, 2013.
  * [8] J. G. Dai, “On positive Harris recurrence of multiclass queueing networks: A unified approach via fluid limit models,” _The Annals of Applied Probability_ , pp. 49–77, 1995.
  * [9] S. Foss and N. Chernova, “On the stability of a partially accessible multi-station queue with state-dependent routing,” _Queueing Systems_ , vol. 29, no. 1, pp. 55–73, 1998.
  * [10] Y. Tang and L. Jin, “Analysis and control of dynamic flow networks subject to stochastic cyber-physical disruptions,” _arXiv preprint arXiv:2004.00159_ , 2020.
  * [11] Y. Sarikaya, T. Alpcan, and O. Ercetin, “Dynamic pricing and queue stability in wireless random access games,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 6, no. 2, pp. 140–150, 2011.
  * [12] P. Dube and R. Jain, “Bertrand games between multi-class queues,” in _Proceedings of the 48h IEEE Conference on Decision and Control (CDC) held jointly with 2009 28th Chinese Control Conference_. IEEE, 2009, pp. 8588–8593.
  * [13] K. Savla and E. Frazzoli, “A dynamical queue approach to intelligent task management for human operators,” _Proceedings of the IEEE_ , vol. 100, no. 3, pp. 672–686, 2011.
  * [14] G. Foschini and J. Salz, “A basic dynamic routing problem and diffusion,” _IEEE Transactions on Communications_ , vol. 26, no. 3, pp. 320–327, 1978\.
  * [15] A. Ephremides, P. Varaiya, and J. Walrand, “A simple dynamic routing problem,” _IEEE transactions on Automatic Control_ , vol. 25, no. 4, pp. 690–693, 1980.
  * [16] N. D. Vvedenskaya, R. L. Dobrushin, and F. I. Karpelevich, “Queueing system with selection of the shortest of two queues: An asymptotic approach,” _Problemy Peredachi Informatsii_ , vol. 32, no. 1, pp. 20–34, 1996.
  * [17] P. Eschenfeldt and D. Gamarnik, “Join the shortest queue with many servers. the heavy-traffic asymptotics,” _Mathematics of Operations Research_ , vol. 43, no. 3, pp. 867–886, 2018.
  * [18] V. Gupta, M. H. Balter, K. Sigman, and W. Whitt, “Analysis of join-the-shortest-queue routing for web server farms,” _Performance Evaluation_ , vol. 64, no. 9-12, pp. 1062–1081, 2007.
  * [19] A. Mukhopadhyay and R. R. Mazumdar, “Analysis of randomized join-the-shortest-queue (jsq) schemes in large heterogeneous processor-sharing systems,” _IEEE Transactions on Control of Network Systems_ , vol. 3, no. 2, pp. 116–126, 2015.
  * [20] S. Mehdian, Z. Zhou, and N. Bambos, “Join-the-shortest-queue scheduling with delay,” in _2017 American Control Conference (ACC)_. IEEE, 2017, pp. 1747–1752.
  * [21] Y. Tang, Y. Wen, and L. Jin, “Security risk analysis of the shorter-queue routing policy for two symmetric servers,” in _2020 American Control Conference (ACC)_. IEEE, 2020, pp. 5090–5095.
  * [22] M. Bramson _et al._ , “Stability of join the shortest queue networks,” _The Annals of Applied Probability_ , vol. 21, no. 4, pp. 1568–1625, 2011\.
  * [23] J. Dai, J. J. Hasenbein, and B. Kim, “Stability of join-the-shortest-queue networks,” _Queueing Systems_ , vol. 57, no. 4, pp. 129–145, 2007.
  * [24] R. D. Foley and D. R. McDonald, “Join the shortest queue: stability and exact asymptotics,” _The Annals of Applied Probability_ , vol. 11, no. 3, pp. 569–607, 2001.
  * [25] X. Ling, M.-B. Hu, R. Jiang, and Q.-S. Wu, “Global dynamic routing for scale-free networks,” _Physical Review E_ , vol. 81, no. 1, p. 016113, 2010\.
  * [26] F. Ren, T. He, S. K. Das, and C. Lin, “Traffic-aware dynamic routing to alleviate congestion in wireless sensor networks,” _IEEE Transactions on Parallel and Distributed Systems_ , vol. 22, no. 9, pp. 1585–1599, 2011.
  * [27] C. H. Papadimitriou and J. N. Tsitsiklis, “The complexity of optimal queueing network control,” in _Proceedings of IEEE 9th Annual Conference on Structure in Complexity Theory_. IEEE, 1994, pp. 318–322.
  * [28] D. Towsley, “Queuing network models with state-dependent routing,” _Journal of the ACM (JACM)_ , vol. 27, no. 2, pp. 323–337, 1980.
  * [29] F. Kelly and C. Laws, “Dynamic routing in open queueing networks: Brownian models, cut constraints and resource pooling,” _Queueing systems_ , vol. 13, no. 1-3, pp. 47–86, 1993.
  * [30] P. Sarachik and U. Ozguner, “On decentralized dynamic routing for congested traffic networks,” _IEEE Transactions on Automatic Control_ , vol. 27, no. 6, pp. 1233–1238, 1982.
  * [31] J. Gregoire, X. Qian, E. Frazzoli, A. De La Fortelle, and T. Wongpiromsarn, “Capacity-aware backpressure traffic signal control,” _IEEE Transactions on Control of Network Systems_ , vol. 2, no. 2, pp. 164–173, 2014\.
  * [32] P. Varaiya, “Max pressure control of a network of signalized intersections,” _Transportation Research Part C: Emerging Technologies_ , vol. 36, pp. 177–195, 2013.
  * [33] J. G. Dai and S. P. Meyn, “Stability and convergence of moments for multiclass queueing networks via fluid limit models,” _IEEE Transactions on Automatic Control_ , vol. 40, no. 11, pp. 1889–1904, 1995.


#  A Mass Ordering Sum Rule for the Neutrino Disappearance Channels
in T2K, NOvA and JUNO

Stephen J. Parke<EMAIL_ADDRESS># orcid: 0000-0003-2028-6782 Theoretical
Physics Department, Fermi National Accelerator Laboratory, Batavia, IL 60510,
USA  Renata Zukanovich Funchal<EMAIL_ADDRESS># orcid:
0000-0001-6749-0022 Instituto de Física, Universidade de São Paulo, 05315-970
São Paulo, Brazil

(April 12, 2024)

###### Abstract

We revisit a method for determining the neutrino mass ordering by using
precision measurements of the atmospheric $\Delta m^{2}$’s in both electron
neutrino and muon neutrino disappearance channels, proposed by the authors in
2005 Nunokawa  _et al._ (2005). The mass ordering is a very important
outstanding question for our understanding of the elusive neutrino and
determination of the mass ordering has consequences for other neutrino
experiments. The JUNO reactor experiment will start data taking this year, and
the precision of the atmospheric $\Delta m^{2}$’s from electron anti-neutrino
measurements will improve by a factor of three from Daya Bay’s 2.4% to 0.8%
within a year. This measurement, when combined with the atmospheric $\Delta
m^{2}$’s measurements from T2K and NOvA for muon neutrino disappearance, will
contribute substantially to the $\Delta\chi^{2}$ between the two remaining
neutrino mass orderings. In this paper we derive a mass ordering sum rule that
can be used to address the possibility that JUNO’s atmospheric $\Delta
m^{2}$’s measurement, when combined with other experiments in particular T2K
and NOvA, can determine the neutrino mass ordering at the 3 $\sigma$
confidence level within one year of operation. For a confidence level of 5
$\sigma$ in a single experiment we will have to wait until the middle of the
next decade when the DUNE experiment is operating.

††preprint: FERMILAB-PUB-24-0117-T

Introduction — We have known for more than a quarter of century that neutrinos
are massive Fukuda  _et al._ (1998) but we still do not know whether the
neutrino with the least amount of $\nu_{e}$, usually labelled $\nu_{3}$, is at
the top or bottom of the neutrino mass spectrum. This is the neutrino mass
ordering question. The SNO experiment Aharmim  _et al._ (2013) determined that
the mass ordering of the other two neutrino mass eigenstates was such that the
neutrino with the most $\nu_{e}$, usually labelled $\nu_{1}$, was lighter than
the other mass state, $\nu_{2}$ which has a smaller $\nu_{e}$ fraction than
$\nu_{1}$ but a larger $\nu_{e}$ fraction than $\nu_{3}$. Thus, the remaining
possible mass ordering for the neutrino mass states is, either
$m_{1}<m_{2}<m_{3}$ which is known as the normal ordering (NO) or
$m_{3}<m_{1}<m_{2}$ which is known as the inverted ordering (IO), see Fig. 1.
The mass squared splitting between $\nu_{2}$ and $\nu_{1}$ was measured with
good precision by the KamLAND experiment to be Gando  _et al._ (2013)

$\displaystyle\Delta m^{2}_{21}\equiv m^{2}_{2}-m^{2}_{1}\approx+7.5\times
10^{-5}~{}{\rm eV}^{2}\,.$ (1)

Whereas the magnitude of the mass-squared-splitting between $\nu_{3}$ and
$\nu_{1}$ has been determined by a number of experiments to be 30 times
larger, i.e.

$\displaystyle\Delta m^{2}_{31}\equiv m^{2}_{3}-m^{2}_{1}\approx\pm\,2.5\times
10^{-3}~{}{\rm eV}^{2}\,,$ (2)

where the ambiguity in the sign comes from the undetermined mass
ordering111For NO $\Delta m^{2}_{atm}>0$ and for IO $\Delta m^{2}_{atm}<0$ for
any atmospheric $\Delta m^{2}$, atm =(ee, $\mu\mu$, 31, 32).. There exist
numerous ways to determine the mass ordering and hence the above sign in the
literature Dighe and Smirnov (2000); Minakata and Nunokawa (2001a, b);
Lunardini and Smirnov (2001); Barger  _et al._ (2002a, b); Petcov and Piai
(2002); Huber  _et al._ (2003); Dighe  _et al._ (2003a, b); Lunardini and
Smirnov (2003); Mena and Parke (2004); Barger  _et al._ (2005); Nunokawa  _et
al._ (2005).

Figure 1: The two remaining mass orderings for the 3 neutrino mass states,
from Parke (2020). Dark blue is the $\nu_{e}$ fraction, cyan the $\nu_{\mu}$
fraction and red the $\nu_{\tau}$ fraction. If the mass state with the least
fraction of $\nu_{e}$, labelled $\nu_{3}$, is at the top of the spectrum this
is called the normal ordering (labelled NO) whereas if it is at the bottom of
the spectrum it is called the inverted ordering (IO). SNO Aharmim  _et al._
(2013) determined the mass ordering of the other two mass states using solar
neutrinos. The set $(\nu_{1},\nu_{2})$, usually called the solar pair, has the
state with most $\nu_{e}$, labelled $\nu_{1}$, below the other member of the
pair, labelled $\nu_{2}$.

Nevertheless, the use of matter effects in neutrino oscillations have been
guiding most of the experimental efforts. The long-baseline accelerator
neutrino experiments NOvA Acero  _et al._ (2023) and T2K Abe  _et al._ (2023a)
as well as the atmospheric neutrinos experiments Super-Kamiokande Wester  _et
al._ (2023) and Ice-Cube Aartsen  _et al._ (2015) operate in a regime where
neutrino oscillations are mostly driven by the higher mass-squared-splitting
and matter effects are significant in the $\nu_{\mu}\to\nu_{e}$ (and
antineutrino) appearance probabilities. They are responsible for the current
status of our understanding of the mass ordering question. At present T2K and
NOvA, individually, have a slight preference for NO, while the combined fit,
flips this preference to IO Kelly  _et al._ (2021); Esteban  _et al._ (2020),
although the preference is weak. Ice-Cube has no preference for either NO or
IO. On the other hand, Super-Kamiokande data seems to favor NO by 92.3% CL
even though their Monte Carlo simulation indicates they should not be able to
discriminate the ordering better than $\sim 80$% CL, see Wester  _et al._
(2023). It is the statistical weight of their data, corresponding to an
exposure of 364.8 kiloton-years, that makes the final global data fit prefer
NO Esteban  _et al._ (2020). So, at the current time, we do not have 3
$\sigma$ or more preference for one ordering over the other.

In this article we revisit one way proposed by the authors of this paper
Nunokawa  _et al._ (2005) which requires precision measurements of $|\Delta
m^{2}_{31}|$ or $|\Delta m^{2}_{32}|$ by both $\nu_{e}$/$\bar{\nu}_{e}$ and
$\nu_{\mu}$/$\bar{\nu}_{\mu}$ disappearance experiments. The reason for
writing this paper now is that within the next year the JUNO experiment An
_et al._ (2016); Abusleme  _et al._ (2022) is expected to improve the
precision of these measurements for $\bar{\nu}_{e}$ disappearance by a factor
of approximately 3, i.e. from Daya Bay’s 2.4 % to better than 0.8 %, a very
significant improvement.

There are a number of papers, such as Blennow and Schwetz (2012); Li  _et al._
(2013); Cabrera  _et al._ (2022), that use the method first presented in
Nunokawa  _et al._ (2005) to address the mass ordering question. Here, we
revisit the combination of long baseline experiments and reactor experiments
by simplifying and extend all of these analyses. First we derive a mass
ordering sum rule, eq. 7, for the electron and muon neutrino disappearance
channels. Then we use NuFIT’s combined T2K and NOvA analysis to generate a
contour plot of the $\Delta\chi^{2}$ between the mass orderings in the plane
of JUNO’s measurement of $\Delta m^{2}_{\rm 31}[\rm NO]$ and its precision,
Fig. 2. This plot allows the reader to immediately estimate the
$\Delta\chi^{2}$ between the two mass orderings as soon as JUNO presents their
$\Delta m^{2}_{\rm 31}|^{\rm NO}$ measurement with its corresponding
precision. Neither the sum rule or the contour plot has appeared elsewhere.

In Nunokawa  _et al._ (2005), with more details available in Parke (2016), we
have shown that the effective atmospheric $\Delta m^{2}$ ($\Delta m^{2}_{\rm
atm}$) for $\nu_{e}$ and $\bar{\nu}_{e}$ disappearance at a baseline divided
by neutrino energy of 0.5 km/GeV, in vacuum, is given by222We use the PDG
conventions for the neutrino mixing matrix.

$\Delta m^{2}_{\rm ee}=\Delta m^{2}_{31}\cos^{2}\theta_{12}+\Delta
m^{2}_{32}\sin^{2}\theta_{12}\,.$ (3)

This is the only $\Delta m^{2}$ measured by Daya Bay An  _et al._ (2017) or
RENO Bak  _et al._ (2018) without additional information from other
experiments. The final precision obtained by Daya Bay (DB) was 2.4%, see An
_et al._ (2023):

$\displaystyle|\Delta m^{2}_{\rm ee}|_{\rm DB}=2.519\pm 0.060~{}~{}{\rm
eV}^{2}\,.$ (4)

Since the Daya Bay experiment is insensitive to the neutrino mass ordering,
the magnitude of $\Delta m^{2}_{\rm ee}$ is the same for both orderings and
the sign is undetermined. Matter effects on the magnitude of $\Delta
m^{2}_{\rm ee}$ are smaller than one tenth of one per cent level, see Khan
_et al._ (2020), and are therefore much smaller than the measurement
uncertainties.

Similarly, for $\nu_{\mu}$ (and $\bar{\nu}_{\mu}$) disappearance, in vacuum,
$\Delta m^{2}_{\mu\mu}$ is given by

$\displaystyle\Delta m^{2}_{\rm\mu\mu}\approx$ $\displaystyle\;\Delta
m^{2}_{31}\sin^{2}\theta_{12}+\Delta m^{2}_{32}\cos^{2}\theta_{12}$
$\displaystyle+\sin\theta_{13}\cos\delta~{}\Delta m^{2}_{21}\,,$ (5)

where $\delta$ is the CP phase and we have set $\sin
2\theta_{12}\tan\theta_{23}\approx 1$ for simplicity in the coefficient of the
$\cos\delta$ term. For the current long baseline experiments such as T2K and
NOvA, even though matter effects are significant for the appearance channels,
it was recently shown by Denton and Parke, in Denton and Parke (2024), that
for the $\nu_{\mu}$ disappearance channels matter effects are at the tenths of
one per cent level. Therefore the $\nu_{\mu}$ disappearance channels are
effectively in vacuum as there is a cancellation between the matter effects in
$\nu_{\mu}\rightarrow\nu_{e}$ and in $\nu_{\mu}\rightarrow\nu_{\tau}$ for the
long baseline experiments T2K and NOvA. As a result these disappearance
channels are also insensitive to the neutrino mass ordering. Similar to Daya
Bay and RENO, $|\Delta m^{2}_{\rm\mu\mu}|$ is the only $\Delta m^{2}$ that is
measurable in T2K and NOvA without information from other experiments.

For NO $\Delta m^{2}_{31}>\Delta m^{2}_{32}$ whereas for IO $|\Delta
m^{2}_{31}|<|\Delta m^{2}_{32}|$, therefore, since $\sin^{2}\theta_{12}\approx
0.3$ we can determine the mass ordering by comparing the magnitude of $|\Delta
m^{2}_{\rm ee}|$ to $|\Delta m^{2}_{\rm\mu\mu}|$:

$\displaystyle|\Delta m^{2}_{\rm ee}|>|\Delta m^{2}_{\rm\mu\mu}|\text{ for
NO}$ $\displaystyle\hskip 42.67912pt{\rm whereas}~{}|\Delta m^{2}_{\rm
ee}|<|\Delta m^{2}_{\rm\mu\mu}|~{}\text{for IO},$ (6)

although the size of the difference between $|\Delta m^{2}_{\rm ee}|$ and
$|\Delta m^{2}_{\rm\mu\mu}|$ is at the couple of percent level and therefore
precise measurements are required 333If T2K and NOvA reported the $|\Delta
m^{2}_{\rm\mu\mu}|$ fits to their disappearance only data, this comparison
could be made more directly without any information needed on $\cos\delta$..
Fortunately we are about to enter an era where very precise measurements will
be made for $\bar{\nu}_{e}$ disappearance by the JUNO experiment.

The above observations can be converted into a “mass ordering sum rule for the
neutrino disappearance channels” :

$\displaystyle(\Delta m^{2}_{31}|^{\rm NO}_{\mu~{}\rm disp}-\Delta
m^{2}_{31}|^{\rm NO}_{e~{}\rm disp})~{}~{}+~{}~{}(|\Delta m^{2}_{32}|^{\rm
IO}_{e~{}\rm disp}-|\Delta m^{2}_{32}|^{\rm IO}_{\mu~{}\rm disp})~{}~{}=$
$\displaystyle~{}~{}(2\cos
2\theta_{12}-2\sin\theta_{13}\,\widehat{\cos\delta})\Delta m^{2}_{21}\,,$ (7)

where the RHS can also be written as
$(2.4-0.9\,\widehat{\cos\delta})\%~{}|\Delta m^{2}_{\rm atm}|$. The subscript
“$\mu$ disp” means the results from $\nu_{\mu}$ disappearance measurements in
T2K and NOvA, whereas “e disp” means the result from $\bar{\nu}_{e}$
disappearance experiments such as Daya Bay and RENO. The symbol
$\widehat{\cos\delta}$ is the average $\cos\delta$ for the NO and IO fits. If
one changes which $\Delta m^{2}$ one uses for both experiments for a given
mass ordering, the RHS of this sum rule is unchanged. For detailed derivation
of this sum rule see the appendix/supplemental material.

From this sum rule it is clear that if NO is Nature’s choice then the first
term will be zero, within measurement uncertainties, and if IO is Nature’s
choice the second term will be zero but in both scenarios the sum of the two
must add up to the RHS, independent of the mass ordering. Consequently, if NO
is Nature’s ordering, the measurements of the $\Delta m^{2}_{32}$’s assuming
IO will not align between $\nu_{e}$ and $\nu_{\mu}$ disappearance within
$\approx$ 2.4% and similar for the IO ordering. The measurement uncertainties
of the experiments Daya Bay, T2K and NOvA are now small enough that this
method is already contributing to the global fits on the neutrino mass
ordering. This can be perceived in the latest NuFIT figure on the synergies
for the $\Delta m^{2}_{3\ell}$’s, see Esteban  _et al._ (2022), one can
observe a preference for NO although the precision of the current measurements
is not sufficient for a 3 $\sigma$ determination of the neutrino mass
ordering.

The question of immediate current interest is how will the precision
measurements of the $\Delta m^{2}_{3i}$’s by JUNO affect the determination of
the mass ordering as JUNO measurements are expected to have an uncertainty
smaller than one third of Daya Bay’s. This measurement is expected to come
very quickly after JUNO turns on, most likely in the first year of operation.

Reactor Measurement ($\bar{\nu}_{e}$ Disappearance) — JUNO is a medium
baseline ($\sim 50$ km) high precision reactor antineutrino oscillation
experiment aiming to determining the neutrino mass ordering by a careful
measurement of the $\bar{\nu}_{e}$ energy spectrum using an idea first
proposed in Petcov and Piai (2002) and further investigated in Choubey  _et
al._ (2003); Bilenky  _et al._ (2017). It was shown in Minakata  _et al._
(2007) that medium baseline reactor experiments can, in principle, determine
the ordering by precisely measuring the effective combination $\Delta
m^{2}_{\rm ee}$ and the sign of a phase ($\pm\Phi_{\odot}$; $+$ for NO, $-$
for IO). This is a very challenging measurement due to various systematic
effects (energy resolution, non-linear detector response etc.) that have to be
tamed and understood to a very high level. The JUNO collaboration claimed, in
their 2015 paper An  _et al._ (2016), that it will take 6 years to determine
the mass ordering at 3 $\sigma$, although more recent papers, see e.g. Forero
_et al._ (2021), have questioned that claim444Two 4.6 GWth were not built at
Taishan and the best fit values of the relevant parameters
$\sin^{2}\theta_{13}$ and $\Delta m^{2}_{atm}$ by Daya Bay, both moved in the
unfavorable directions, see Fig 8. of Forero  _et al._ (2021), compared to
what was used in An  _et al._ (2016).  . JUNO’s recent update Abusleme  _et
al._ (2022) does not contain an update on their expected mass ordering
sensitivity.

On the other hand, JUNO is expected to reach, after a few months of operation,
and much sooner than they can start to be sensitive to the mass ordering,
unprecedented sub-percent precision on the determination of $|\Delta
m^{2}_{\rm atm}|$ where $|\Delta m^{2}_{\rm atm}|$ could be any one of
$|\Delta m^{2}_{\rm ee}|$, or $|\Delta m^{2}_{31}|$, or $|\Delta m^{2}_{32}|$
depending on the experiment’s analysis choice for both mass orderings but all
are related to one another. JUNO claims that after 100 days of data taking
they will be able to determine $|\Delta m^{2}_{\rm atm}|$ at $0.8\%$ precision
and will continue to improve ultimately reaching 0.2% precision, see Abusleme
_et al._ (2022). In contrast, the expected $\Delta\chi^{2}$ between the NO and
IO in JUNO’s fits is expected to grow quite slowly, at no more than 1.5 units
per year.

It was shown in Forero  _et al._ (2021), that due to the phase advance (NO) or
retardation (IO) the best fits to the spectrum at the far JUNO (JU) detector
will give a $|\Delta m^{2}_{\rm ee}|$ for IO which is 0.7% larger than
$|\Delta m^{2}_{\rm ee}|$ fit for NO, i.e.

$\displaystyle|\Delta m^{2}_{\rm ee}|^{\rm IO}_{\rm JU}$ $\displaystyle=\Delta
m^{2}_{\rm ee}|^{\rm NO}_{\rm JU}+1.8\times 10^{-5}~{}{\rm eV}^{2}\,.$ (8)

Note this shift is a fraction of $\Delta m^{2}_{21}$ and comes from the fact
that in the IO the atmospheric oscillations are retarded with respect to the
NO and therefore $\Delta m^{2}_{ee}$ for IO is slightly larger than that for
NO to compensate for this phase shift. We have used that $0.007\times|\Delta
m^{2}_{\rm ee}|\approx 0.018\times 10^{-3}{\rm eV}^{2}$ and the value of
$\Delta m^{2}_{21}$ is given by eq. 1. Using eq. 3 for each mass ordering we
can relate this result for any of the other possible $|\Delta m^{2}_{\rm
atm}|$ as given in Table 1.

$\Delta(\Delta m^{2})$ | NO ee | NO 31 | NO 32
---|---|---|---
IO ee | 1.8 | -0.5 | 6.9
IO 31 | -0.5 | -2.7 | 4.7
IO 32 | 6.9 | 4.7 | 12.1
Table 1: This Table relates the different $\Delta m^{2}_{\rm atm}$’s for the
best fits of the JUNO experiment, such that each entry gives the difference
between $\Delta(\Delta m^{2})\equiv|\Delta m^{2}_{ij}|^{\rm IO}-\Delta
m^{2}_{kl}|^{\rm NO}$ in units of $10^{-5}$ eV2 where (ij) and (kl) are (ee,
31, 32). A value of 2.5 in this table represents a 1.0 % difference. The
difference between $|\Delta m^{2}_{32}|^{\rm IO}$ and $\Delta m^{2}_{32}|^{\rm
NO}$ is $\sim 5$%, much larger than the measurement error expected in JUNO,
whereas the difference between $|\Delta m^{2}_{31}|^{\rm IO}$ and $\Delta
m^{2}_{31}|^{\rm NO}$ is $\sim 1$%.

This modifies the “mass ordering sum rule for the neutrino disappearance
channels” , given in eq. 7, by increasing the RHS by 0.7%, so that for T2K,
NOvA (LBL) and JUNO (JU) we have the following sum rule

$\displaystyle(\Delta m^{2}_{31}|^{\rm NO}_{\rm LBL}-\Delta m^{2}_{31}|^{\rm
NO}_{\rm JU})+(|\Delta m^{2}_{32}|^{\rm IO}_{\rm JU}-|\Delta m^{2}_{32}|^{\rm
IO}_{\rm LBL})$ $\displaystyle\hskip
56.9055pt\approx(3.1-0.9\,\widehat{\cos\delta})\%~{}|\Delta m^{2}_{\rm
atm}|\,.$ (9)

With enough precision on the measurements of the $\Delta m^{2}_{\rm atm}$’s
between $\nu_{e}$ and $\nu_{\mu}$ disappearance, the $\Delta\chi^{2}$ between
the two mass ordering fits, can contribute significantly to the determination
of the mass ordering.

We can describe for JUNO the $\chi^{2}$ fit to data to determine $|\Delta
m^{2}_{\rm atm}|$, as the parabola which will depend on the assumed best fit
value $\Delta m^{2}_{\rm atm}|^{\rm NO}$ or $|\Delta m^{2}_{\rm atm}|^{\rm
IO}$ and $\sigma_{\rm JU}$, the precision of the measurement, which does not
depend on the ordering

$\displaystyle\chi^{2}_{\rm NO}(\Delta m^{2},\sigma_{\rm JU})_{\rm JU}$
$\displaystyle=\left(\frac{\Delta m^{2}-\Delta m^{2}_{\rm atm}|^{\rm
NO}}{\sigma_{\rm JU}}\right)^{2}$ (10) $\displaystyle\chi^{2}_{\rm IO}(|\Delta
m^{2}|,\sigma_{\rm JU})_{\rm JU}$ $\displaystyle=\left(\frac{|\Delta
m^{2}|-|\Delta m^{2}_{\rm atm}|^{\rm IO}}{\sigma_{\rm JU}}\right)^{2}$ (11)

Note that the subscript ‘atm’ can be any one of (ee, 31, 32) and they do not
need to be the same for NO and IO. However, the best fit points must be
related by the numbers given in Table 1. Different experiments and different
global fit groups may make different choices. Here, we will use the choice
used by the NuFIT collaboration, that is $\Delta m^{2}_{31}$ for NO and
$\Delta m^{2}_{32}$ for IO and the relationship that for the JUNO experiment

$\displaystyle|\Delta m^{2}_{32}|^{\rm IO}_{\rm JU}$ $\displaystyle=\Delta
m^{2}_{31}|^{\rm NO}_{\rm JU}+4.7\times 10^{-5}~{}{\rm eV}^{2}\,.$ (12)

The physics conclusions will be independent of the (31, 32) arbitrary choices.

Long-baseline Accelerator Measurement ($\nu_{\mu}/\bar{\nu}_{\mu}$
Disappearance) — T2K is a long-baseline ($\sim 295$ km) accelerator neutrino
oscillation experiment in Japan that has collected a total of $1.97\times
10^{21}$ and $1.63\times 10^{21}$ protons on target in neutrino and
antineutrino modes, respectively. Similarly, NOvA is a long-baseline ($\sim
810$ km) accelerator neutrino oscillation experiment in the US that also has
collected $1.36\times 10^{21}$ and $1.25\times 10^{21}$ protons on target of
data in neutrino and antineutrino modes, respectively. Both experiments
operate as a $\nu_{\mu}\to\nu_{\mu}/\bar{\nu}_{\mu}\to\bar{\nu}_{\mu}$
disappearance experiment as well as a
$\nu_{\mu}\to\nu_{e}/\bar{\nu}_{\mu}\to\bar{\nu}_{e}$ appearance experiment.
Their disappearance measurements have no sensitivity to the mass ordering but
are the responsible for the precise determination of $|\Delta m^{2}_{32}|$ (or
equivalently $|\Delta m^{2}_{31}|$). T2K results given in Abe  _et al._
(2023b) are $\Delta m^{2}_{32}|^{\rm NO}=(2.49\pm 0.05)\times 10^{-3}\,\rm
eV^{2}$ and $|\Delta m^{2}_{31}|^{\rm IO}=(2.46\pm 0.05)\times 10^{-3}\,\rm
eV^{2}$, note that uncertainty is $\sim$ 2%. NOvA’s results are given by
$\Delta m^{2}_{32}|^{\rm NO}=(2.39\pm 0.06)\times 10^{-3}\,\rm eV^{2}$ and
$|\Delta m^{2}_{32}|^{\rm IO}=(2.44\pm 0.06)\times 10^{-3}\,\rm eV^{2}$, see
Acero  _et al._ (2023). Given the consistency of T2K and NOvA disappearance
measurements, they can be combined , as in Esteban  _et al._ (2020) , for both
a NO and a IO fit. Using the update of these fits, given in Esteban  _et al._
(2022), we have for NO

$\displaystyle\chi^{2}_{\rm NO}(\Delta m^{2}_{31})_{\rm LBL}$
$\displaystyle=\left(\frac{\Delta m^{2}_{31}-\Delta m^{2}_{\rm 31}|^{\rm
NO}}{\sigma^{\rm NO}_{\rm LBL}}\right)^{2}$ (13)

with $(\Delta m^{2}_{31}|^{\rm NO},\sigma^{\rm NO}_{\rm
LBL})=(2.516,0.031)\times 10^{-3}$ eV2. Similarly for the IO fit,

$\chi^{2}_{\rm IO}(|\Delta m^{2}_{32}|)_{\rm LBL}=\left(\frac{|\Delta
m^{2}_{32}|-|\Delta m^{2}_{\rm 32}|^{\rm IO}}{\sigma^{\rm IO}_{\rm
LBL}}\right)^{2}$ (14)

with $(|\Delta m^{2}_{32}|^{\rm IO},\sigma^{\rm IO}_{\rm
LBL})=(2.485,0.031)\times 10^{-3}$ eV2. Note since the uncertainty for both NO
and IO are the same, $\sigma^{\rm IO}_{\rm LBL}=\sigma^{\rm NO}_{\rm LBL}$, we
can drop the mass ordering on this symbol. Also it is important that this
combined uncertainty is $\sim$ 1.2%, much less than the RHS of eq. 7 and 9.
T2K’s and NOvA’s results are not expected to change significantly in the next
few years due to the large statistics already collected by these experiments.

Figure 2:  Iso-contours of $\Delta\chi^{2}_{\rm min}\equiv\chi^{2}_{\rm
IO}|_{\rm min}-\chi^{2}_{\rm NO}|_{\rm min}$ in the plane of JUNO’s $(\Delta
m^{2}_{31}|^{\rm NO}\otimes\sigma)$ fit to the far detector spectrum assuming
normal ordering (NO). The precision is expressed as a % of $\Delta
m^{2}_{31}|^{\rm NO}$. The blue region favors NO whereas the red region favors
IO. The current values from the Daya Bay (DB) experiment and from Daya Bay and
RENO combined according to NuFIT (React. Comb.) are also shown. Table 1 can be
used to translate the horizontal axis to any other $\Delta m^{2}_{\rm atm}$.
The white dashed lines mark the precision achievable by JUNO after 100 days
(0.8%) and 6 years (0.2%).

Combining JUNO with T2K and NOvA — We combine the results of T2K, NOvA and
JUNO by just adding the $\chi^{2}$ for the combine LBL results to that of
JUNO, using the best fit for NO and the measurement precision of JUNO as
variables as follows:

$\displaystyle\chi^{2}_{\rm NO}(\Delta m^{2}_{31},\sigma_{\rm JU})$
$\displaystyle=\chi^{2}_{\rm NO}(\Delta m^{2}_{31},\sigma_{\rm JU})_{\rm
JU}+\chi^{2}_{\rm NO}(\Delta m^{2}_{31})_{\rm LBL}$
$\displaystyle\chi^{2}_{\rm IO}(|\Delta m^{2}_{32}|,\sigma_{\rm JU})$
$\displaystyle=\chi^{2}_{\rm IO}(|\Delta m^{2}_{32}|,\sigma_{\rm JU})_{\rm
JU}+\chi^{2}_{\rm IO}(|\Delta m^{2}_{32}|)_{\rm LBL}$

where for JUNO we use the relationship given in eq. 12 for the best fit
values. Then the difference in the $\Delta\chi^{2}$ is given by

$\displaystyle\Delta\chi^{2}(\Delta m^{2}_{31}|^{\rm NO}_{\rm JU},\sigma_{\rm
JU})_{\rm min}$ $\displaystyle=\chi^{2}_{\rm IO}|_{\rm min}-\chi^{2}_{\rm
NO}|_{\rm min}\,.$ (15)

Now everything is determined except for the best fit value $\Delta
m^{2}_{31}|^{\rm NO}$ and precision of this measurement $\sigma_{\rm JU}$ by
the JUNO experiment. In Fig. 2 we have plotted the $\Delta\chi^{2}_{\rm min}$
defined as the minimum for IO minus the minimum for NO as a function of JUNO’s
$\Delta m^{2}_{31}|^{\rm NO}$ and the fractional precision of the measurement
by JUNO555 We fix the solar parameters at their best fit values determined by
NuFIT Esteban  _et al._ (2020) $\Delta m^{2}_{21}=7.41\times 10^{-5}$ eV2 and
$\sin^{2}\theta_{12}=0.303$. Small changes in these parameters when measured
by JUNO will not effect our results as the experimental uncertainty on $\Delta
m^{2}_{\rm atm}$ from the LBL experiments will dominate. . On this fan shaped
figure we have shown for reference Daya Bay’s result as well as the one from
the combined fit of Daya Bay and RENO data given by the NuFIT group. Clearly
as JUNO’s precision on the measurement gets better and better one moves down
this plot which increases the $\Delta\chi^{2}$, however the central value
could also move right (left) thereby increasing (decreasing) the
$\Delta\chi^{2}$. This plot assumes that the disappearance $\Delta m^{2}_{\rm
atm}$ results from T2K and NOvA will not significantly change from what they
are now. So, as soon as JUNO presents results on the $|\Delta m^{2}_{\rm
atm}|$ fit to the spectrum at the far detector, one can read off from this
plot the contribution that $\nu_{\mu}/\bar{\nu}_{\mu}$ and $\bar{\nu}_{e}$
disappearance measurements from T2K, NOvA and JUNO make to the
$\Delta\chi^{2}$ for the mass ordering determination. If JUNO has the same
central value as the NuFIT (React. Comb.) value but with a precision below 1%,
then the contribution from these disappearance measurements will be greater
than 9 units of $\Delta\chi^{2}$ , i.e. greater than 3$\sigma$.

Conclusions — We re-examine an idea we have had on how to determine the
neutrino mass ordering, almost 20 years ago, using only neutrino disappearance
data in vacuum. This is auspicious today in the light of the current precision
on $|\Delta m^{2}_{\rm atm}|$ (31 for NO and 32 for IO) achieved by long-
baseline experiments ($\sim$2% individually and $\sim$1.2% combined) and of
the imminent few per mil determination of $|\Delta m^{2}_{\rm atm}|$ by JUNO.
This kind of unprecedented accuracy allows one to discuss a mass ordering sum
rule for the neutrino disappearance channels (see eq. 9) which may be used to
determine the mass ordering in the near future solely using data from these
disappearance experiments. To show this with more clarity, we combined in a
$\chi^{2}$ function the present results from T2K and NOvA $\nu_{\mu}$ and
$\bar{\nu}_{\mu}$ disappearance measurements with that expected from
$\bar{\nu}_{e}$ disappearance at JUNO, as a function of the assumed best fit
value ($\Delta m^{2}_{31}|^{\rm NO}_{\rm JU}$) and fractional accuracy of the
JUNO measurement ($\sigma_{\rm JU}$). In Fig. 2 one can see the main result,
the values of $\Delta\chi^{2}_{\rm min}=\chi^{2}_{\rm IO}|_{\rm
min}-\chi^{2}_{\rm NO}|_{\rm min}$ in the plane of JUNO’s $(\Delta
m^{2}_{31}|^{\rm NO}\otimes\sigma)$. This is a very useful figure because as
soon as JUNO presents its first result on their $|\Delta m^{2}_{\rm atm}|$
measurement, one can read from it if NO is preferred, and if so, how much it
is preferred in terms of how many units of $\Delta\chi^{2}$. In this manner it
is conceivable, if JUNO measures $|\Delta m^{2}_{\rm atm}|$ close to the one
given by combining Daya Bay and RENO data, that NO could be soon (in a year or
so) determined, by the combined (T2K, NOvA and JUNO) disappearance
measurements alone, to better than 3 $\sigma$ i.e. a confidence level of
99.73%.

###### Acknowledgements.

SJP would like to thank Peter Denton for many useful discussions on matter
effects for $\nu_{\mu}$ disappearance. SJP acknowledges support by the United
States Department of Energy under Grant Contract No. DE-AC02-07CH11359. R. Z.
F. is partially supported by Fundação de Amparo à Pesquisa do Estado de São
Paulo (FAPESP) under Contract No. 2019/04837-9, and Conselho Nacional de
Desenvolvimento Científico e Tecnológico (CNPq). We would both like to thank
our long time collaborator, Hiroshi Nunokawa, for many enlightening
discussions on this topic over many years. Its a pity he chose not to join us
for this follow on paper of Nunokawa  _et al._ (2005). This project has also
received support from the European Union’s Horizon 2020 research and
innovation programme under the Marie Sklodowska-Curie grant agreement No
860881-HIDDeN as well as under the Marie Skłodowska-Curie Staff Exchange grant
agreement No 101086085 - ASYMMETRY.

## References

  * Nunokawa  _et al._ (2005) H. Nunokawa, S. J. Parke, and R. Zukanovich Funchal, Phys. Rev. D 72, 013009 (2005), arXiv:hep-ph/0503283 .
  * Fukuda  _et al._ (1998) Y. Fukuda _et al._ (Super-Kamiokande), Phys. Rev. Lett. 81, 1562 (1998), arXiv:hep-ex/9807003 .
  * Aharmim  _et al._ (2013) B. Aharmim _et al._ (SNO), Phys. Rev. C 88, 025501 (2013), arXiv:1109.0763 [nucl-ex] .
  * Gando  _et al._ (2013) A. Gando _et al._ (KamLAND), Phys. Rev. D 88, 033001 (2013), arXiv:1303.4667 [hep-ex] .
  * Note (1) For NO $\Delta m^{2}_{atm}>0$ and for IO $\Delta m^{2}_{atm}<0$ for any atmospheric $\Delta m^{2}$, atm =(ee, $\mu\mu$, 31, 32).
  * Dighe and Smirnov (2000) A. S. Dighe and A. Y. Smirnov, Phys. Rev. D 62, 033007 (2000), arXiv:hep-ph/9907423 .
  * Minakata and Nunokawa (2001a) H. Minakata and H. Nunokawa, Phys. Lett. B 504, 301 (2001a), arXiv:hep-ph/0010240 .
  * Minakata and Nunokawa (2001b) H. Minakata and H. Nunokawa, JHEP 10, 001 (2001b), arXiv:hep-ph/0108085 .
  * Lunardini and Smirnov (2001) C. Lunardini and A. Y. Smirnov, Nucl. Phys. B 616, 307 (2001), arXiv:hep-ph/0106149 .
  * Barger  _et al._ (2002a) V. Barger, D. Marfatia, and K. Whisnant, Phys. Rev. D 65, 073023 (2002a), arXiv:hep-ph/0112119 .
  * Barger  _et al._ (2002b) V. Barger, D. Marfatia, and K. Whisnant, Phys. Rev. D 66, 053007 (2002b), arXiv:hep-ph/0206038 .
  * Petcov and Piai (2002) S. T. Petcov and M. Piai, Phys. Lett. B 533, 94 (2002), arXiv:hep-ph/0112074 .
  * Huber  _et al._ (2003) P. Huber, M. Lindner, and W. Winter, Nucl. Phys. B 654, 3 (2003), arXiv:hep-ph/0211300 .
  * Dighe  _et al._ (2003a) A. S. Dighe, M. T. Keil, and G. G. Raffelt, JCAP 06, 006 (2003a), arXiv:hep-ph/0304150 .
  * Dighe  _et al._ (2003b) A. S. Dighe, M. T. Keil, and G. G. Raffelt, JCAP 06, 005 (2003b), arXiv:hep-ph/0303210 .
  * Lunardini and Smirnov (2003) C. Lunardini and A. Y. Smirnov, JCAP 06, 009 (2003), arXiv:hep-ph/0302033 .
  * Mena and Parke (2004) O. Mena and S. J. Parke, Phys. Rev. D 70, 093011 (2004), arXiv:hep-ph/0408070 .
  * Barger  _et al._ (2005) V. Barger, P. Huber, and D. Marfatia, Phys. Lett. B 617, 167 (2005), arXiv:hep-ph/0501184 .
  * Parke (2020) S. Parke, in  _28th International Symposium on Lepton Photon Interactions at High Energies_ (WSP, Singapur, 2020) pp. 116–127, arXiv:1801.09643 [hep-ph] .
  * Acero  _et al._ (2023) M. A. Acero _et al._ (NOvA), (2023), arXiv:2311.07835 [hep-ex] .
  * Abe  _et al._ (2023a) K. Abe _et al._ (T2K), Phys. Rev. D 108, 072011 (2023a), arXiv:2305.09916 [hep-ex] .
  * Wester  _et al._ (2023) T. Wester _et al._ (Super-Kamiokande), (2023), arXiv:2311.05105 [hep-ex] .
  * Aartsen  _et al._ (2015) M. G. Aartsen _et al._ (IceCube), Phys. Rev. D 91, 072004 (2015), arXiv:1410.7227 [hep-ex] .
  * Kelly  _et al._ (2021) K. J. Kelly, P. A. N. Machado, S. J. Parke, Y. F. Perez-Gonzalez, and R. Z. Funchal, Phys. Rev. D 103, 013004 (2021), arXiv:2007.08526 [hep-ph] .
  * Esteban  _et al._ (2020) I. Esteban, M. C. Gonzalez-Garcia, M. Maltoni, T. Schwetz, and A. Zhou, JHEP 09, 178 (2020), arXiv:2007.14792 [hep-ph] .
  * An  _et al._ (2016) F. An _et al._ (JUNO), J. Phys. G 43, 030401 (2016), arXiv:1507.05613 [physics.ins-det] .
  * Abusleme  _et al._ (2022) A. Abusleme _et al._ (JUNO), Chin. Phys. C 46, 123001 (2022), arXiv:2204.13249 [hep-ex] .
  * Blennow and Schwetz (2012) M. Blennow and T. Schwetz, JHEP 08, 058 (2012), [Erratum: JHEP 11, 098 (2012)], arXiv:1203.3388 [hep-ph] .
  * Li  _et al._ (2013) Y.-F. Li, J. Cao, Y. Wang, and L. Zhan, Phys. Rev. D 88, 013008 (2013), arXiv:1303.6733 [hep-ex] .
  * Cabrera  _et al._ (2022) A. Cabrera _et al._ , Sci. Rep. 12, 5393 (2022), arXiv:2008.11280 [hep-ph] .
  * Parke (2016) S. Parke, Phys. Rev. D 93, 053008 (2016), arXiv:1601.07464 [hep-ph] .
  * Note (2) We use the PDG conventions for the neutrino mixing matrix.
  * An  _et al._ (2017) F. P. An _et al._ (Daya Bay), Phys. Rev. D 95, 072006 (2017), arXiv:1610.04802 [hep-ex] .
  * Bak  _et al._ (2018) G. Bak _et al._ (RENO), Phys. Rev. Lett. 121, 201801 (2018), arXiv:1806.00248 [hep-ex] .
  * An  _et al._ (2023) F. P. An _et al._ (Daya Bay), Phys. Rev. Lett. 130, 161802 (2023), arXiv:2211.14988 [hep-ex] .
  * Khan  _et al._ (2020) A. N. Khan, H. Nunokawa, and S. J. Parke, Phys. Lett. B 803, 135354 (2020), arXiv:1910.12900 [hep-ph] .
  * Denton and Parke (2024) P. B. Denton and S. J. Parke, Phys. Rev. D 109, 053002 (2024), arXiv:2401.10326 [hep-ph] .
  * Note (3) If T2K and NOvA reported the $|\Delta m^{2}_{\rm\mu\mu}|$ fits to their disappearance only data, this comparison could be made more directly without any information needed on $\mathop{cos}\nolimits\delta$.
  * Esteban  _et al._ (2022) I. Esteban, M. C. Gonzalez-Garcia, M. Maltoni, T. Schwetz, and A. Zhou, “Nufit figure on the synergies for the $\delta m^{2}_{3\ell}$’s,”  (2022), http://www.nu-fit.org/sites/default/files/v52.fig-chisq-dma.pdf.
  * Choubey  _et al._ (2003) S. Choubey, S. T. Petcov, and M. Piai, Phys. Rev. D 68, 113006 (2003), arXiv:hep-ph/0306017 .
  * Bilenky  _et al._ (2017) S. M. Bilenky, F. Capozzi, and S. T. Petcov, Phys. Lett. B 772, 179 (2017), [Erratum: Phys.Lett.B 809, 135765 (2020)], arXiv:1701.06328 [hep-ph] .
  * Minakata  _et al._ (2007) H. Minakata, H. Nunokawa, S. J. Parke, and R. Zukanovich Funchal, Phys. Rev. D 76, 053004 (2007), [Erratum: Phys.Rev.D 76, 079901 (2007)], arXiv:hep-ph/0701151 .
  * Forero  _et al._ (2021) D. V. Forero, S. J. Parke, C. A. Ternes, and R. Z. Funchal, Phys. Rev. D 104, 113004 (2021), arXiv:2107.12410 [hep-ph] .
  * Note (4) Two 4.6 GWth were not built at Taishan and the best fit values of the relevant parameters $\mathop{sin}\nolimits^{2}\theta_{13}$ and $\Delta m^{2}_{atm}$ by Daya Bay, both moved in the unfavorable directions, see Fig 8. of Forero  _et al._ (2021), compared to what was used in An  _et al._ (2016).
  * Abe  _et al._ (2023b) K. Abe _et al._ (T2K), Eur. Phys. J. C 83, 782 (2023b), arXiv:2303.03222 [hep-ex] .
  * Note (5) We fix the solar parameters at their best fit values determined by NuFIT Esteban  _et al._ (2020) $\Delta m^{2}_{21}=7.41\times 10^{-5}$ eV2 and $\mathop{sin}\nolimits^{2}\theta_{12}=0.303$. Small changes in these parameters when measured by JUNO will not effect our results as the experimental uncertainty on $\Delta m^{2}_{\rm atm}$ from the LBL experiments will dominate.

##  Appendix A Supplemental Material

Here we derive the “mass ordering sum rule for neutrino disappearance
channels” from the observations of Nunokawa  _et al._ (2005). Daya Bay and
Reno measure the same $|\Delta m^{2}_{ee}|$ for both mass orderings. Therefore
from eq. 3 we have the following:

$\displaystyle\Delta m^{2}_{31}|^{\rm NO}_{e~{}\rm disp}$
$\displaystyle=|\Delta m^{2}_{ee}|+\sin^{2}\theta_{12}\Delta m^{2}_{21}$
$\displaystyle|\Delta m^{2}_{32}|^{\rm IO}_{e~{}\rm disp}$
$\displaystyle=|\Delta m^{2}_{ee}|+\cos^{2}\theta_{12}\Delta m^{2}_{21}$
$\displaystyle|\Delta m^{2}_{32}|^{\rm IO}_{e~{}\rm disp}-\Delta
m^{2}_{31}|^{\rm NO}_{e~{}\rm disp}$ $\displaystyle=\cos 2\theta_{12}\Delta
m^{2}_{21}.$ (16)

Daya Bay’s reported measurements given in An  _et al._ (2023) of $\Delta
m^{2}_{32}|^{\rm NO}_{\rm DB}=2.466\pm 0.060\times 10^{-3}$ eV2 and $|\Delta
m^{2}_{32}|^{\rm IO}_{\rm DB}=2.571\pm 0.060\times 10^{-3}$ eV2 satisfy eq.
16, after adding $\Delta m^{2}_{21}$ to $\Delta m^{2}_{32}|^{\rm NO}$ to
obtain $\Delta m^{2}_{31}|^{\rm NO}$. The agreement is much, much smaller than
the quoted uncertainty, demonstrating that these two results are connected via
$|\Delta m^{2}_{ee}|$.

Similarly for the disappearance channels in T2K and NOvA, they both measure
the same $|\Delta m^{2}_{\mu\mu}|$ for both mass orderings and matter effects
are negligible in this channel for both T2K and NOvA, see Denton and Parke
(2024). Therefore from eq. 5 we have the following:

$\displaystyle\Delta m^{2}_{31}|^{\rm NO}_{\mu~{}\rm disp}$
$\displaystyle=|\Delta
m^{2}_{\mu\mu}|+(\cos^{2}\theta_{12}-\sin\theta_{13}\cos\delta^{\rm NO})\Delta
m^{2}_{21}$ $\displaystyle|\Delta m^{2}_{32}|^{\rm IO}_{\mu~{}\rm disp}$
$\displaystyle=|\Delta
m^{2}_{\mu\mu}|+(\sin^{2}\theta_{12}+\sin\theta_{13}\cos\delta^{\rm IO})\Delta
m^{2}_{21}$ $\displaystyle\Delta m^{2}_{31}|^{\rm NO}_{\mu~{}\rm disp}-|\Delta
m^{2}_{32}|^{\rm IO}_{\mu~{}\rm disp}$ $\displaystyle=(\cos
2\theta_{12}-2\sin\theta_{13}\,\widehat{\cos\delta})\Delta m^{2}_{21}.$ (17)

where $\widehat{\cos\delta}\equiv\frac{1}{2}(\cos\delta^{\rm
NO}+\cos\delta^{\rm IO})$. T2K provides enough significant figures such that
again the identity eq. 17 can be checked. T2K results given in Table 13 of Abe
_et al._ (2023b) are reported as $\Delta m^{2}_{32}|^{\rm NO}_{\rm
T2K}=2.494\pm 0.054\times 10^{-3}$ eV2 and $|\Delta m^{2}_{31}|^{\rm IO}_{\rm
T2K}=2.463\pm 0.049\times 10^{-3}$ eV2 again after correcting for the
$31\leftrightarrow 32$ for both mass ordering, we get excellent agreement when
using $\widehat{\cos\delta}\approx 0$. The agreement is much, much smaller
than the quoted uncertainty, demonstrating that these two results are
connected via $|\Delta m^{2}_{\mu\mu}|$. NOvA, unfortunately, does not provide
enough significant figures for their measurements to demonstrate the $|\Delta
m^{2}_{\mu\mu}|$ connection as convincingly as T2K, but their measurements are
in agreement with 17.

Now, adding and rearranging eq. 16 and 17 we have the “mass ordering sum rule
for neutrino disappearance channels” given in eq. 7:

$\displaystyle(\Delta m^{2}_{31}|^{\rm NO}_{\mu~{}\rm disp}-\Delta
m^{2}_{31}|^{\rm NO}_{e~{}\rm disp})~{}~{}+~{}~{}(|\Delta m^{2}_{32}|^{\rm
IO}_{e~{}\rm disp}-|\Delta m^{2}_{32}|^{\rm IO}_{\mu~{}\rm disp})~{}~{}$
$\displaystyle=~{}~{}(2\cos
2\theta_{12}-2\sin\theta_{13}\,\widehat{\cos\delta})\Delta m^{2}_{21}\,,$ (18)
$\displaystyle=~{}~{}(2.4-0.9\,\widehat{\cos\delta})\%~{}|\Delta m^{2}_{\rm
atm}|\,,$

where for the last line we have used $\Delta m^{2}_{21}/|\Delta m^{2}_{\rm
atm}|=0.03$, $\sin^{2}\theta_{12}=0.3$ and $\sin\theta_{13}=0.15$. Note the
interchange between “$\mu$ disp” and “e disp” when going from NO to IO. Also
this sum is invariant if we replace $\Delta m^{2}_{31}$ with $\Delta
m^{2}_{32}$ for NO or/and $\Delta m^{2}_{32}$ with $\Delta m^{2}_{31}$ for IO.
We use the above choice because that is the choice made by NuFIT, but our
physics conclusions are independent of this choice. In passing it is worth
noting that the NuFIT results given in Esteban  _et al._ (2022) suggest that
$\widehat{\cos\delta}\leq 0$.

Due to the phase advance (NO) or retardation (IO) of the atmospheric
oscillations, JUNO does not measure exactly the same $|\Delta m^{2}_{ee}|$ for
both mass orderings, see Forero  _et al._ (2021), infact

$\displaystyle|\Delta m^{2}_{ee}|^{\rm IO}\approx 1.007\times\Delta
m^{2}_{ee}|^{\rm NO}=\Delta m^{2}_{ee}|^{\rm NO}+1.8\times 10^{-5}{\rm
eV}^{2}.$ (19)

This changes eq. 16 and also the “mass ordering sum rule for neutrino
disappearance channels” to what is given in eq. 9 by adding 0.7 % times
$|\Delta m^{2}_{\rm atm}|$ to the RHS:

$\displaystyle(\Delta m^{2}_{31}|^{\rm NO}_{\rm LBL}-\Delta m^{2}_{31}|^{\rm
NO}_{\rm JU})+(|\Delta m^{2}_{32}|^{\rm IO}_{\rm JU}-|\Delta m^{2}_{32}|^{\rm
IO}_{\rm LBL})\approx(3.1-0.9\,\widehat{\cos\delta})\%~{}|\Delta m^{2}_{\rm
atm}|\,.$ (20)

Here the label “e disp” (“ $\mu$ disap”) has been replaced with the label “JU”
(“LBL”) as this sum rule is specific for T2K, NOvA (LBL) and JUNO (JU). The
additional 0.7 % has an impact on Fig. 2 because of the precision of the
JUNO’s measurements on $\Delta m^{2}_{\rm atm}$’s.

If Nature’s choice for the mass ordering is NO, then

$\displaystyle(\Delta m^{2}_{31}|^{\rm NO}_{\rm LBL}-\Delta m^{2}_{31}|^{\rm
NO}_{\rm JU})\approx 0\quad{\rm and}\quad(|\Delta m^{2}_{32}|^{\rm IO}_{\rm
JU}-|\Delta m^{2}_{32}|^{\rm IO}_{\rm
LBL})\approx(3.1-0.9\,\widehat{\cos\delta})\%~{}|\Delta m^{2}_{\rm atm}|$ (21)

whereas for IO

$\displaystyle(\Delta m^{2}_{31}|^{\rm NO}_{\rm LBL}-\Delta m^{2}_{31}|^{\rm
NO}_{\rm JU})\approx(3.1-0.9\,\widehat{\cos\delta})\%~{}|\Delta m^{2}_{\rm
atm}|\quad{\rm and}\quad(|\Delta m^{2}_{32}|^{\rm IO}_{\rm JU}-|\Delta
m^{2}_{32}|^{\rm IO}_{\rm LBL})\approx 0\,,$ (22)

where “$\approx$” should be interpreted to mean “within measurement
uncertainties”. The current measurement uncertainty for the combined T2K and
NOvA measurement is $\sim$ 1.2%. JUNO’s measurement uncertainty will reach 0.8
% within one year of data taking and reach 0.2 % within 6 years. If T2K and
NOvA can improve their measurement uncertainties on the $\Delta m^{2}_{\rm
atm}$’s even modestly, in the future, this can impact the $\Delta\chi^{2}$
between the two mass orderings substantially.