text:
\begin{document}

\title{On the Benefits of Network-level Cooperation in IoT Networks with Aggregators}

\author{Nikolaos Pappas~\IEEEmembership{Member,~IEEE}, Ioannis Dimitriou,
         Zheng Chen~\IEEEmembership{Member,~IEEE}
\thanks{ N. Pappas is with the Department of Science and Technology, Link\"{o}ping University, Norrk\"{o}ping SE-60174, Sweden. (e-mail: nikolaos.pappas@liu.se). I. Dimitriou is with the Department of Mathematics, University of Patras, Patra, Peloponnese, Greece.
(e-mail: idimit@math.upatras.gr). Z. Chen is with the Department of Electrical Engineering, Link\"{o}ping University, 58183 Link\"{o}ping, Sweden.  (e-mail: zheng.chen@liu.se). This work was supported in part by ELLIIT, Center for Industrial Information Technology (CENIIT), and the Swedish Foundation for Strategic Research. This work was presented in part in the 30th International Teletraffic Congress (ITC) 2018 \cite{ITC30}.}
}

\maketitle

\begin{abstract}
In this work, we consider a random access Internet of Things IoT wireless network assisted by two aggregators collecting information from two disjoint groups of sensors. The nodes and the aggregators are transmitting in a random access manner under slotted time, the aggregators perform network-level cooperation for the data collection. The aggregators are equipped with queues to store data packets that are transmitted by the network nodes and relaying them to the destination node. We characterize the throughput performance of the IoT network and we obtain the stability conditions for the queues at the aggregators. We apply the theory of boundary value problems to analyze the delay performance. Our results show that the presence of the aggregators provides significant gains in the IoT network performance, in addition, we provide useful insights regarding the scalability of the IoT network.
\end{abstract}

\begin{IEEEkeywords}
Random Access, IoT, Throughput, Delay, Queueing, Stability, Boundary Value Problems.
\end{IEEEkeywords}

\IEEEpeerreviewmaketitle

\section{Introduction}
The Internet of Things (IoT) is one of the most attractive concepts in the area of information and communication technology. IoT is expected to play an important role in our daily life by supporting massive connectivity with seamless service. It involves the interconnection of many, and possibly heterogeneous objects through the Internet using different communication technologies. The objects are equipped with communications capabilities and can vary from sensors, smart objects, etc. \cite{StankovicIoT2014, KhanWirComMag2017, LenkaACCESS2018}.

The total number of IoT connections is expected to grow tremendously the next years. The application of random access protocols in the IoT networks can potentially mitigate the co-channel interference caused by massive amount of IoT devices with low signaling overhead \cite{HasanComMag2013}.
To support the massive connectivity in future IoT networks, practical techniques are required to collect data from a large set of devices and the traditional orthogonal multi-access schemes are not efficient. The works in \cite{MalakTCOM2016, HsuICC2017, GuoTCOM2017, LopezTWC2018, SalamAccess2018, ZhangTVT2018, HattabGC2018, Rahman2018QoS} have considered data aggregation under different scenarios and have evaluated the benefits of such technique. The work in \cite{Globecom2017Latency} proposes a cloud-assisted priority-based channel access and data aggregation scheme for irregularly deployed sensor nodes to minimize the network latency and to enhance the system reliability of IoT networks. The work in \cite{ZhangTWC2019} considers an uplink IoT network where a cellular user operating as aggregator can assist IoT devices. IoT can leverage the unlicensed spectrum for machine-to-machine communications. The problem of spectrum allocation and device association in uplink two-hop IoT networks is studied in \cite{IbrahimIoT2019}.

Network-level cooperation introduced in \cite{SadekTIT2007} and \cite{RongTIT12} can be an effective alternative method for data aggregation. It is plain relaying without any physical layer considerations and it has been shown to provide large gain in terms of throughput and delay performance.
Recently, several works have investigated relaying at the network level \cite{PappasTWC2015,PappasGlobalsip2013, pap, PappasJCN2016, DimitriouASMTA2016, NomikosSurvey2016,dim3, PappasJCN2018, DIMITRIOU2018, CristianGC2018, Dimitriou2019AdHoc}. The deployment of aggregators under network-level cooperation can improve the throughput and delay in IoT networks \cite{KimTVT2018}. However, due to the queueing delay, the stability conditions at the aggregator queues need to be considered.

In this work, we consider a random access IoT network with two aggregators which can receive and forward data in form of packets from the network nodes.
The nodes transmit in a random access manner, which is a common assumption in IoT technologies such as LoRa \cite{LoRa2016}. The aggregators use network-level cooperation and their receiving and transmitting modules are operating in different frequency bands, thus, they have out-of-band full duplex capabilities. Furthermore, we assume multi-packet reception (MPR) capabilities at the receivers. MPR is suitable to capture the SINR (Signal-to-Interference-plus-Noise-Ratio) based model and is more realistic to model the wireless transmission. Similar assumptions to our work can be found in \cite{MalakTCOM2016, SalamAccess2018, KimTVT2018}.

\subsection{Contributions}

The contributions of our work can be summarized as follows.
Our primary goal is to provide support for data collection in IoT networks by applying network-level cooperation. We provide the throughput analysis of an IoT network consisting of sensors that are assisted by two aggregators, from which we can gain insights on the scalability of the considered network. In addition, we study the stability conditions for the queues at the aggregators, which guarantee finite queueing delay. The main difficulty for characterizing the stability conditions lies on the interaction of the queues, which arises when the service rate of a queue depends on the state of the other. The technique of stochastic dominance is utilized to bypass this difficulty and characterize necessary and sufficient stability conditions for the queues at the aggregators.

Furthermore, we study the average delay of the packets possibly received and forwarded by the aggregators. Our system is modeled as a two-dimensional discrete time Markov chain, and we show that the generating function of the stationary joint queue length distribution can be obtained by solving a fundamental functional equation with the theory of boundary value problems. For the general MPR case, we obtain a lower and an upper bound for the average delay. In addition, we characterize in closed form expression the gap between the lower and the upper bound. These bounds as it is also seen in the numerical results appear to be tight. For the model with capture effect, i.e., a subclass of MPR model, we provide the explicit expression for the average delay. The analysis in this work can act as a framework for other research directions that involve multiple aggregators with interacting queues.

The remainder of the paper is organized as follows: Section \ref{sec:model} describes the system model and in Section \ref{sec:anal1} we present the analysis related to the network-wide throughput. In Sections \ref{anal} and \ref{sec:delaysym}, we provide the analysis for the average delay per packet. The numerical evaluation of the theoretical results is presented in Section \ref{sec:Res}, and Section \ref{sec:conclusions} concludes the paper.

\section{System Model} \label{sec:model}

\subsection{Network Layer Model}
We consider a wireless network consisting of IoT nodes/sensors/objects, which intend to communicate to a common destination/sink $D$, and two aggregators, denoted by $R_{1}$ and $R_{2}$, which can help aggregating and relaying messages from the IoT nodes to $D$. The network model is depicted in Fig. \ref{fig:model}. The nodes are located in two non-overlapping regions. There are in total $M_1$ and $M_2$ nodes within the service range of the aggregators $R_1$ and $R_2$, respectively.
Note that the transmissions in the first region cannot be overheard by the nodes in the second region and vice versa. This can be done by carefully planning for the placement of the aggregators in order to increase the coverage area without interfering with each other. \textit{In the following, we will use the terms nodes, sensors, and objects interchangeably}.

The sensors intend to transmit packets to the destination node $D$ and they are assumed to be saturated, i.e., they always have packets to transmit. In case the transmission from a sensor to $D$ fails, the aggregator can help relaying the message to $D$ and the aggregators do not generate their own traffic. We consider using \textit{network-level cooperation} at the aggregators \cite{SadekTIT2007, RongTIT12}, which means that the aggregators are cooperating as relays in a decode-and-forward manner. The packets are assumed to have equal length and the time is divided into slots, which corresponds to the transmission time of a packet. We assume that the sensors access the wireless channel randomly without any coordination among them. We consider a full multi-packet reception (MPR) channel model, which allows the receivers to successfully receive more than one packets when there are multiple transmissions in the same slot \cite{verdu}. As the result, the sink node $D$ can receive information simultaneously from the sensors and the aggregators. Note that when all nodes are transmitting, we can have in total up to $M_1+M_2$ interfering devices at the sink. This assumption is common in the literature \cite{MalakTCOM2016, KimTVT2018, ZhangTVT2018}.

We assume that different frequency bands are allocated for the transmissions from the sensors and the aggregators, thus, there is no interference between them. On the other hand, the transmission of a node creates interference to the other nodes of the same kind, i.e., there is interference between the sensors, and between the aggregators. The transmitting and receiving units of the aggregators are operating in different channels/frequency bands to avoid self-interference, which can be considered as \textit{out-of-band full duplex mode} \cite{ramakrishnan2016design}.
The aggregators are equipped with queues that store possible packets from the sensors that failed to reach the destination. The queues are assumed to have infinite capacity, thus there is no packet dropping.\footnote{In practice, the buffers have limited size, which is usually quite large. Our analysis based on the infinite buffer size assumption can capture this scenario with minor modifications and it is more general. Furthermore, our analysis can provide design guidelines for selecting the appropriate queue size in order to achieve specific performance requirements regarding throughput and delay. The arrival and service rates for the queues are defined in Section \ref{sec:anal1}.} This is a common assumption in the literature, in the IoT context \cite{KimTVT2018}.

\begin{figure}[!ht]
\centering
\includegraphics[scale=0.25]{system_model}\vspace{-0.1in}
\caption{The IoT network considered in this work. The nodes are assisted by two aggregators that are equipped with queues, and they are operating in an out-of-band full duplex mode.}
\label{fig:model}
\end{figure}

\subsection{Physical Layer Model} \label{sec:PHY}
At the beginning of a timeslot, sensor nodes that belong to the coverage area of $R_i$, attempt to transmit with probability $t_i$, $i=1,2$. The aggregator $R_{i}$ will attempt to transmit a packet with probability $\alpha_{i}$, if it has a non-empty queue. Note that we assume that all the sensors in the same area transmit with the same probability. Our analysis can be easily extended to handle the general case where each node has different transmit probabilities however, the expressions will be cumbersome, and it will be difficult to extract meaningful insights.

We assume that a packet transmitted by sensor $s$ in the first coverage area is successfully received by its aggregator $R_1$ if and only if ${\rm SINR}(s,R_1)\geq \gamma$, where $\gamma$ is the SINR threshold. The wireless channel is subject to fading; let $P_{tx}(s)$ be the transmit power at sensor $s$ and $r(s,R_1)$ be the distance between the sensor and the aggregator. The received power at $R_1$, when $s$ transmits is $P_{rx}(s,R_1)=A(s,R_1)h(s,R_1)$ where $A(s,R_1)$ is a random variable representing small-scale fading. Under Rayleigh fading assumption, $A(s,R_1)$ is exponentially distributed. The received power factor $h(s,R_1)$ is given by $h(s,R_1) = P_{tx}(s)(r(s,R_1))^{-\theta}$ where $\theta$ is the path loss exponent. Then, the success probability between a sensor in the first coverage area and its aggregator is denoted by $P^{R_1}_{{i}}$, when there are $i$ sensors from the same area transmitting in a timeslot, and the expression is given by \cite{PappasTWC2015}
\begin{equation}
\begin{aligned}
\label{eq:succprob}
P^{R_1}_{{i}}=\exp\left(-\frac{\gamma \eta_{R_1}}{v(s,R_1)h(s,R_1)}\right) \prod_{k\in \mathcal{T}\backslash \left\{s\right\}}{\left(1+\gamma\frac{v(k,R_1)h(k,R_1)}{v(s,R_1)h(s,R_1)}\right)}^{-1},
\end{aligned}
\end{equation}
where $\mathcal{T}$ is the set of transmitting sensors in the first area and $|\mathcal{T}|=i$; $v(s,R_1)$ is the parameter of the Rayleigh fading random variable, $\eta_{R_1}$ is the receiver noise power at $R_1$.

Note that transmitting sensors from the other coverage area do not create interference at the aggregator.
However, the concurrently transmitting sensors from both areas interfere with each other at the destination $D$. Denote by $P^{1D}_{{i},{j}}$ the success probability to the sink from a sensor in coverage area $1$ when there are $i$ active transmitters from area 1 and $j$ active transmitters from area 2. Similarly we can define $P^{2D}_{{i},{j}}$.

A packet transmission from a sensor in the first area fails to reach the destination with probability $\bar{P}^{1D}_{{i},{j}}=1-P^{1D}_{{i},{j}}$, when there are $i$ active sensors in the first area and $j$ active sensors in the second area. In this case, that packet will be stored in the queue of aggregator $R_{1}$ with probability $P^{R_1}_{{i}}$. Otherwise, with probability $\bar{P}^{R_1}_{{i}}$, the aggregator fails to decode that packet and it has to be re-transmitted by the sensor in a future time slot.

Recall that if there are stored packets in the queues of the aggregator $R_{i}$, $i=1,2$, then $R_i$ transmits a packet with probability $\alpha_{i}$. If only one aggregator $R_{i}$ is active, then the packet will be successfully transmitted to $D$ with probability $p_{R_{i},\{R_{i}\}}^{D}$. If both aggregators transmit simultaneously, then with probability $p_{R_i/R_1,R_2}^{D}$ the packet from $R_i$ is successfully received by node $D$. If a transmitted packet from an aggregator fails to reach the $D$, that packet remains in the queue and will be retransmitted in a later time slot.

\section{Throughput and Stability Analysis} \label{sec:anal1}

In this section, we characterize the network throughput performance and provide the stability conditions for the queues at the aggregators.

\subsection{Throughput Analysis}
The throughput per node consists of the direct throughput from each sensor to the destination and the throughput contributed by the aggregator.
Recall that the devices that are in coverage from the first aggregator cannot cause interference at the receiver of the second aggregator. Moreover, the devices that are covered by the $i$-th aggregator, $R_i$, are transmitting with probability $t_{i}$, for $i=1,2$.

The direct throughput from a sensor in the first coverage area to the sink is given by
\begin{equation}
\begin{aligned}
T_{1,D}=\sum_{i=0}^{M_1-1}\sum_{j=0}^{M_2}{M_1-1 \choose i}{M_2 \choose j} {t_1^{i+1}\overline{t}_1^{M_1-i-1}} {t_2^{j}\overline{t}_2^{M_2-j}} P^{1D}_{{i+1},{j}}.
\end{aligned}
\end{equation}
Note that the direct throughput in this setup is equivalent to the throughput in the network without aggregators.

The contributed throughput from a sensor to the aggregator in the first coverage area is given by
\begin{equation}
\begin{aligned}
T_{1,R_1}=\sum_{i=0}^{M_1-1}\sum_{j=0}^{M_2}{M_1-1 \choose i}{M_2 \choose j} {t_1^{i+1}\overline{t}_1^{M_1-i-1}} {t_2^{j}\overline{t}_2^{M_2-j}} \overline{P}^{1D}_{{i+1},{j}}P^{R_1}_{i+1}.
\end{aligned}
\end{equation}
The total throughput seen by a sensor in the first coverage area is
$T_1=T_{1,D}+T_{1,R_1}$. Similarly, we can obtain the throughput seen by a sensor located in the second coverage area which is assisted by the second aggregator $R_2$.
The percentage of throughput

\begin{remark}
The term $T_{1,D}$ can also be interpreted as the probability that a transmitted packet from a sensor in the first coverage area reaches the destination directly. Furthermore, $T_{1,R_1}$ is the probability of unsuccessful transmission from a sensor to the destination while the packet is received by the aggregator. The percentage of a sensor's traffic that is being relayed is $\frac{T_{1,R_1}}{T_1}$.
\end{remark}

In addition, we need to characterize the average arrival rates at the aggregators, denoted by $\lambda_1$ and $\lambda_2$. Since we assume full MPR capability at the receivers, the $i$-th aggregator can receive up to $N_i$ packets in a timeslot.
We define $l_{k,i}$ as the probability that $k$ packets will arrive in a timeslot at the $i$-th aggregator.
The average arrival rate at the $i$-th aggregator is given by
\begin{equation} \label{eq:lambdaju}
\lambda_{i} = \sum_{k=1}^{M_i} k l_{k,i}, \text{ }i=1,2.
\end{equation}
The probability $l_{k,1}$ where $1 \leq k \leq N_1$ is given by
\begin{equation}
l_{k,1}=\sum_{s=k}^{M_1}\sum_{m=0}^{M_2}{M_1 \choose s}{s \choose k}{M_2 \choose m} {t_1^{s}\overline{t}_1^{M_1-s}} {t_2^{m}\overline{t}_2^{M_2-m}}  \left( \overline{P}^{1D}_{{s},{m}}P^{R_1}_{s} \right)^{k} \left(P^{1D}_{{s},{m}} + \overline{P}^{1D}_{{s},{m}}\overline{P}^{R_1}_{s} \right)^{s-k}.
\end{equation}

Similarly we can obtain $l_{k,2}$.
The network-wide throughput is $T=M_1 T_1+M_2 T_2$. The previous expressions for the throughput are valid when the queues at the aggregators are stable.

\begin{remark}
Previously we considered the case where the queues are stable, but in order to obtain the throughput when the queues are not stable we needs to replace the sum of expressions $T_{i,R_i}$ with the service rate of the aggregator. In this case the network-wide throughput will be given by $T=M_1 T_{1,D}+M_2 T_{2,D}+\mu_1+\mu_2$.
\end{remark}

\subsection{Stability Analysis at the Aggregators} \label{sec:stability_conditions}
The average service rate for the aggregator $i$ is given by
\begin{equation} \label{eq:mu}
\begin{array}{l}
\mu_i=Pr(N_j \neq 0) \left[\alpha_i\bar{\alpha}_{j}p_{R_{i},\{R_{i}\}}^{D}+\alpha_i \alpha_j p_{R_i/R_i,R_j}^{D}\right] + Pr(N_j = 0) \alpha_{i} p_{R_{i},\{R_{i}\}}^{D},j=i\text{ }\mathrm{mod}2+1,
\end{array}
\end{equation}
where $N_{j}$ is the queue size at queue $j$.
The notation for the success probabilities used in (\ref{eq:mu}) is the one introduced in Section \ref{sec:model}.
The theorem below, provides the stability conditions for the considered setup.

\begin{thm}\label{thm1}
The stability conditions for the queues at the aggregators for fixed transmission probabilities $\alpha_{i}, i=1,2$, are described by the region $\mathcal{R}=\mathcal{R}_1 \cup \mathcal{R}_2$ where $\mathcal{R}_i$ is given by
\begin{align} \label{eq:stabCond}
\mathcal{R}_i = \left\lbrace(\lambda_1,\lambda_2):\lambda_i \leq \alpha_i p_{R_{i},\{R_{i}\}}^{D} - \frac{\alpha_i \alpha_j(p_{R_{i},\{R_{i}\}}^{D}-p_{R_i/R_i,R_j}^{D})}{\alpha_j\bar{\alpha}_ip_{R_{j},\{R_{j}\}}^{D}+\alpha_i \alpha_j p_{R_j/R_i,R_j}^{D}}\lambda_j, \right. \\ \nonumber \left.  \lambda_j \leq \alpha_j\bar{\alpha}_ip_{R_{j},\{R_{j}\}}^{D}+\alpha_i \alpha_jp_{R_j/R_i,R_j}^{D} \right\rbrace, j=i\text{ }\mathrm{mod}2+1.
\end{align}
\end{thm}

\begin{proof}
The proof is given in Appendix \ref{sec:stability_proof}.
\end{proof}

\section{Preliminary Analysis}\label{anal}
Let $N_{k,n}$ be the number of packets in the buffer of aggregator $R_{k}$, $k=1,2$, at the beginning of the $n$th slot. Then, $Y_{n}=(N_{1,n},N_{2,n})$ is a discrete time Markov chain with state space $E=\{(i,j):i,j=0,1,2,...\}$. The queues of both aggregators evolve as follows:
\begin{equation}
\begin{array}{c}
N_{k,n+1}=[N_{k,n}+F_{k,n}]^{+},\,k=1,2,
\end{array}
\label{x}
\end{equation}
where $F_{k,n}$ is the the difference of the number of packets that enter the buffer of the $k$th aggregator at the beginning of slot $n$ ($F_{k,n}$ equals $0$ or $\pm 1$), and $[x]^{+}:=max(0,x)$.

Before proceeding with the analysis, we will slightly modify the notation for the success probabilities presented in Section \ref{sec:model} in order to be more convenient for the delay analysis. If two sensors transmit a packet simultaneously, $p_{i,\{1,2\}}^{D}$ represents the probability that the packet from sensor $i$ is successfully received by node $D$ and the packet from sensor $j=i mod 2+1$ failed to be received by $D$. $p_{1,2,\{1,2\}}^{D}$ represents the probability that the packets from both nodes are successfully received by $D$. Then we have $\bar{p}_{\{1,2\}}^{D}=1-p_{1,2,\{1,2\}}^{D}-p_{1,\{1,2\}}^{D}-p_{2,\{1,2\}}^{D}$,which denotes the probability that both packets fail to be received by the node $D$. If both aggregators transmit simultaneously, then with probability $p_{R_{k},\{R_{1},R_{2}\}}^{D}$ the packet from $R_{k}$ is successfully received by node $D$, with probability $p_{R_{1},R_{2},\{R_{1},R_{2}\}}^{D}$ the packets from both aggregators are successfully received by $D$, while with probability $\bar{p}_{\{R_1,R_2\}}^{D}=1-p_{R_{1},R_{2},\{R_{1},R_{2}\}}^{D}-p_{R_{1},\{R_{1},R_{2}\}}^{D}-p_{R_{2},\{R_{1},R_{2}\}}^{D}$, both of them failed to be received by $D$.

Let $H(x,y)$ be the generating function of the joint stationary queue process,
\begin{displaymath}
\begin{array}{c}
H(x,y)=\lim_{n\to\infty}E(x^{N_{1,n}}y^{N_{2,n}}),\,|x|\leq1,|y|\leq1.
\end{array}
\end{displaymath}
Then, by exploiting (\ref{x}) we obtain after lengthy calculations
\begin{equation}
\begin{array}{r}
R(x,y)H(x,y)=A(x,y)H(x,0)+B(x,y)H(0,y)+C(x,y)H(0,0),
\end{array}
\label{we}
\end{equation}
where,
\begin{displaymath}
\begin{array}{rl}
R(x,y):=1-L(x,y)[1-\alpha_{1}\widehat{\alpha}_{2}(1-\frac{1}{x})-\widehat{\alpha}_{1}\alpha_{2}(1-\frac{1}{y})-\alpha_{1}\alpha_{2}p_{R_{1},R_{2},\{R_{1},R_{2}\}}^{D}(1-\frac{1}{xy})],
\end{array}
\end{displaymath}
is the \textit{kernel} of functional equation (\ref{we}) and
\begin{displaymath}
\begin{array}{rl}
A(x,y)=&L(x,y)[d_{1}(1-\frac{1}{x})+\alpha_{2}\widehat{\alpha}_{1}(1-\frac{1}{y})+\alpha_{1}\alpha_{2}p_{R_{1},R_{2},\{R_{1},R_{2}\}}^{D}(1-\frac{1}{xy})],\\
B(x,y)=&L(x,y)[d_{2}(1-\frac{1}{y})+\alpha_{1}\widehat{\alpha}_{2}(1-\frac{1}{x})+a_{1}a_{2}p_{R_{1},R_{2},\{R_{1},R_{2}\}}^{D}(1-\frac{1}{xy})],\\
C(x,y)=&L(x,y)[d_{1}(\frac{1}{x}-1)+d_{2}(\frac{1}{y}-1)+\alpha_{1}\alpha_{2}p_{R_{1},R_{2},\{R_{1},R_{2}\}}^{D}(\frac{1}{xy}-1)],\\
L(x,y)=&1-(1-x)t_{1}[\bar{t}_{2}\bar{p}_{1,\{1\}}^{D}p_{1,\{1\}}^{R_{1}}+t_{2}(\bar{p}_{\{1,2\}}^{D}+p_{2,\{1,2\}}^{D})p_{1,\{1,2\}}^{R_{1}}]+(1-y)t_{2}[\bar{t}_{1}\bar{p}_{2,\{2\}}^{D}p_{2,\{2\}}^{R_{2}}\\&+t_{1}(\bar{p}_{\{1,2\}}^{D}+p_{1,\{1,2\}}^{D})p_{1,\{1,2\}}^{R_{2}}]+(1-x)(1-y)t_{1}t_{2}\bar{p}_{\{1,2\}}^{D}p_{1,\{1,2\}}^{R_{1}}p_{2,\{1,2\}}^{R_{2}},
\end{array}
\end{displaymath}

\begin{displaymath}
\begin{array}{rl}
\widehat{\alpha}_{k}=&\bar{\alpha}_{k}p_{R_k,\{R_k\}}^{D}+\alpha_{k}p_{R_k,\{R_1,R_2\}}^{D},\,k=1,2,\\
d_{1}=&\alpha_{1}(\widehat{\alpha}_{2}-p_{R_{1},\{R_{1}\}}^{D}),
d_{2}=\alpha_{2}(\widehat{\alpha}_{1}-p_{R_{2},\{R_{2}\}}^{D}).
\end{array}
\end{displaymath}
Clearly, for every fixed $y$ with $|y|\leq 1$, $H(x,y)$ it is regular in $x$ for $|x|<1$, and continuous in $x$ for $|x|\leq1$; similar observation hold for the variable $y$.

Some interesting relations are directly obtained using (\ref{we}). In particular, by setting in (\ref{we}) $y=1$, dividing with $x-1$, and taking the limit $x\to 1$, by using the L'Hospital rule, and vice-versa we obtain the following conservation of flow relations:
\begin{equation}
\begin{array}{rl}
\lambda_{1}=&\alpha_{1}\tilde{\alpha}_{2}[1-H(1,0)-H(0,1)+H(0,0)]+\alpha_{1}p_{R_{1},\{R_{1}\}}^{D}[H(1,0)-H(0,0)],\\
\lambda_{2}=&\alpha_{2}\tilde{\alpha}_{1}[1-H(1,0)-H(0,1)+H(0,0)]+\alpha_{2}p_{R_{2},\{R_{2}\}}^{D}[H(0,1)-H(0,0)],
\end{array}
\label{con}
\end{equation}
where, $\tilde{\alpha}_{k}:=\widehat{\alpha}_{k}+\alpha_{k}p_{R_{1},R_{2},\{R_{1},R_{2}\}}^{D}$, $k=1,2,$
\begin{displaymath}
\begin{array}{l}
\lambda_{1}:=t_{1}\bar{t}_{2}\bar{p}_{1,\{1\}}^{D}p_{1,\{1\}}^{R_{1}}+t_{1}t_{2}(\bar{p}_{1,2,\{1,2\}}^{D}+p_{2,\{1,2\}}^{D})p_{1,\{1,2\}}^{R_{1}},\\
\lambda_{2}:=t_{2}\bar{t}_{1}\bar{p}_{2,\{2\}}^{D}p_{2,\{2\}}^{R_{2}}+t_{1}t_{2}(\bar{p}_{1,2,\{1,2\}}^{D}+p_{1,\{1,2\}}^{D})p_{2,\{1,2\}}^{R_{2}}.
\end{array}
\end{displaymath}
Note that the previous equations are the same with the ones obtained in the previous section, but here we use the more convenient notation for the delay analysis regarding the success probabilities. In order to facilitate presentation we present the case of two nodes, but clearly the analysis holds for the general case of $N$ nodes, just by replacing the right parts of $\lambda_{1}$ and $\lambda_{2}$.

The expressions in (\ref{con}), equate the flow of jobs into an aggregator, with the flow of jobs out of the aggregator. Looking carefully at (\ref{con}) it is readily seen that the following analysis is distinguished in two cases:
\begin{enumerate}
\item $\frac{\alpha_{1}\tilde{\alpha}_{2}}{\alpha_{1}p_{R_{1},\{R_{1}\}}^{D}}+\frac{\alpha_{2}\tilde{\alpha}_{1}}{\alpha_{2}p_{R_{2},\{R_{2}\}}^{D}}=1$. Then, $H(0,0)=1-\frac{\lambda_{1}}{\alpha_{1}p_{R_{1},\{R_{1}\}}^{D}}-\frac{\lambda_{2}}{\alpha_{2}p_{R_{2},\{R_{2}\}}^{D}}=1-\rho$.

\item $\frac{\alpha_{1}\tilde{\alpha}_{2}}{\alpha_{1}p_{R_{1},\{R_{1}\}}^{D}}+\frac{\alpha_{2}\tilde{\alpha}_{1}}{\alpha_{2}p_{R_{2},\{R_{2}\}}^{D}}\neq1$. Then, (\ref{con}) yields $H(1,0)=\frac{\alpha_{1}\tilde{\alpha}_{2}(\lambda_{2}-\alpha_{2}p_{R_{2},\{R_{2}\}}^{D})-\tilde{d}_{2}(\lambda_{1}+H(0,0)\alpha_{1}p_{R_{1},\{R_{1}\}}^{D})}{\tilde{d}_{1}\tilde{d}_{2}-\alpha_{1}\alpha_{2}\tilde{\alpha}_{1}\tilde{\alpha}_{2}},$ and $H(0,1)=\frac{\alpha_{2}\tilde{\alpha}_{1}(\lambda_{1}-\alpha_{1}p_{R_{1},\{R_{1}\}}^{D})-\tilde{d}_{1}(\lambda_{2}+H(0,0)\alpha_{2}p_{R_{2},\{R_{2}\}}^{D})}{\tilde{d}_{1}\tilde{d}_{2}-\alpha_{1}\alpha_{2}\tilde{\alpha}_{1}\tilde{\alpha}_{2}},$

where $\tilde{d}_{1}=\alpha_{1}(\tilde{\alpha}_{2}-p_{R_{1},\{R_{1}\}}^{D})$, $\tilde{d}_{2}=a\alpha_{2}(\tilde{\alpha}_{1}-p_{R_{2},\{R_{2}\}}^{D})$.
\end{enumerate}
\subsection{Kernel analysis}
The kernel $R(x,y)$ plays a crucial role in the following analysis and here we provide some important properties. For convenience, assume in the following that $p_{R_{1},R_{2},\{R_{1},R_{2}\}}^{D}=0$. It is readily seen that
\begin{displaymath}
R(x,y)=a(x)y^{2}+b(x)y+c(x)=\widehat{a}(y)x^{2}+\widehat{b}(y)x+\widehat{c}(y),
\end{displaymath}
where, for $L=t_{1}t_{2}\bar{p}_{\{1,2\}}^{D}p_{1,\{1,2\}}^{R_{1}}p_{2,\{1,2\}}^{R_{2}}$,
\begin{displaymath}
\begin{array}{rl}
\widehat{a}(y)=&L(1-\alpha_{1}\widehat{\alpha}_{2}-\alpha_{2}\widehat{\alpha}_{1})y^{2}+y[\lambda_{1}(a_{1}\widehat{\alpha}_{2}+\alpha_{2}\widehat{\alpha}_{1}-1)\\&+L(\alpha_{1}\widehat{\alpha}_{2}+2\alpha_{2}\widehat{\alpha}_{1}-1)]-\alpha_{2}\widehat{\alpha}_{1}(\lambda_{1}+L),\\
\widehat{b}(y)=&y^{2}[\lambda_{2}(\alpha_{1}\widehat{\alpha}_{2}+\alpha_{2}\widehat{\alpha}_{1}-1)+L(2\alpha_{1}\widehat{\alpha}_{2}+\alpha_{2}\widehat{\alpha}_{1}-1)]\\
&+y[\lambda_{1}(1-2\alpha_{1}\widehat{\alpha}_{2}-\alpha_{2}\widehat{\alpha}_{1})+\lambda_{2}(1-2\widehat{\alpha}_{1}\alpha_{2}-\alpha_{1}\widehat{\alpha}_{2})\\
&+L(1-2(\widehat{\alpha}_{1}\alpha_{2}+\alpha_{1}\widehat{\alpha}_{2}))+\alpha_{1}\widehat{\alpha}_{2}+\alpha_{2}\widehat{\alpha}_{1}]\\&+\widehat{\alpha}_{1}\alpha_{2}(\lambda_{1}+\lambda_{2}+L-1),\\
\widehat{c}(y)=&\alpha_{1}\widehat{\alpha}_{2}y[\lambda_{1}-1+(\lambda_{2}+L)(1-y)],\\
a(x)=&L(1-\alpha_{1}\widehat{\alpha}_{2}-\alpha_{2}\widehat{\alpha}_{1})x^{2}+x[\lambda_{2}(\alpha_{1}\widehat{\alpha}_{2}+\alpha_{2}\widehat{\alpha}_{1}-1)\\&+L(2\alpha_{1}\widehat{\alpha}_{2}+\alpha_{2}\widehat{\alpha}_{1}-1)]-\alpha_{1}\widehat{\alpha}_{2}(\lambda_{2}+L),\\
b(x)=&x^{2}[\lambda_{1}(\alpha_{1}\widehat{\alpha}_{2}+\alpha_{2}\widehat{\alpha}_{1}-1)+L(2\alpha_{2}\widehat{\alpha}_{1}+\alpha_{1}\widehat{\alpha}_{2}-1)]\\
&+x[\lambda_{1}(1-2\alpha_{1}\widehat{\alpha}_{2}-\alpha_{2}\widehat{\alpha}_{1})+\lambda_{2}(1-2\widehat{\alpha}_{1}\alpha_{2}-\alpha_{1}\widehat{\alpha}_{2})\\
&+L(1-2(\widehat{\alpha}_{1}\alpha_{2}+\alpha_{1}\widehat{\alpha}_{2}))+\alpha_{1}\widehat{\alpha}_{2}+\alpha_{2}\widehat{\alpha}_{1}]\\&+\widehat{\alpha}_{2}\alpha_{1}(\lambda_{1}+\lambda_{2}+L-1),\\
c(x)=&\alpha_{2}\widehat{\alpha}_{1}x[\lambda_{2}-1+(\lambda_{1}+L)(1-x)].
\end{array}
\end{displaymath}
The roots of $R(x,y)=0$ are $X_{\pm}(y)=\frac{-\widehat{b}(y)\pm\sqrt{D_{y}(y)}}{2\widehat{a}(y)}$, $Y_{\pm}(x)=\frac{-b(x)\pm\sqrt{D_{x}(x)}}{2a(x)}$, where $D_{y}(y)=\widehat{b}(y)^{2}-4\widehat{a}(y)\widehat{c}(y)$, $D_{x}(x)=b(x)^{2}-4a(x)c(x)$.

\begin{lemma}\label{LEM}
For $|y|=1$, $y\neq1$, $R(x,y)=0$ has exactly one root $x=X_{0}(y)$ such that $|X_{0}(y)|<1$. When $\lambda_{1}<\alpha_{1}\widehat{\alpha}_{2}$, $X_{0}(1)=1$. Similarly, $R(x,y)=0$ has exactly one root $y=Y_{0}(x)$, such that $|Y_{0}(x)|\leq1$, for $|x|=1$.
\end{lemma}
\begin{proof}
See Appendix \ref{a0}.
\end{proof}

The lemma below provides information about the location of the branch points of the two-valued functions $Y(x)$, $X(y)$, its proof is based on algebraic arguments, thus it is omitted.

\begin{lemma}\label{lem1}
The algebraic function $Y(x)$, defined by $R(x,Y(x)) = 0$, has four real branch points $0< x_{1}<x_{2}\leq1<x_{3}<x_{4}$. Moreover, $D_{x}(x)<0$, $x\in(x_{1},x_{2})\cup(x_{3},x_{4})$ and $D_{x}(x)>0$, $x\in(-\infty,x_{1})\cup(x_{2},x_{3})\cup(x_{4},\infty)$. Similarly, $X(y)$, is defined by $R(X(y),y) = 0$, it has four real branch points $0\leq y_{1}<y_{2}\leq1<y_{3}<y_{4}$, and $D_{x}(y)<0$, $y\in(y_{1},y_{2})\cup(y_{3},y_{4})$ and $D_{x}(y)>0$, $y\in(-\infty,y_{1})\cup(y_{2},y_{3})\cup(y_{4},\infty)$.
\end{lemma}
Let $\tilde{C}_{x}=\mathbb{C}_{x}-([x_{1},x_{2}]\cup[x_{3},x_{4}]$, $\tilde{C}_{y}=\mathbb{C}_{y}-([y_{1},y_{2}]\cup[y_{3},y_{4}]$, where $\mathbb{C}_{x}$, $\mathbb{C}_{y}$ the complex planes of $x$, $y$, respectively. In $\tilde{C}_{x}$ (resp. $\tilde{C}_{y}$), denote by $Y_{0}(x)$ (resp. $X_{0}(y)$) the root of $R(x,Y(x))=0$ (resp. $R(X(y),y)=0$) with the smallest modulus, and $Y_{1}(x)$ (resp. $X_{1}(y)$) the other one.
Define the image contours, $\mathcal{L}=Y_{0}[\overrightarrow{\underleftarrow{x_{1},x_{2}}}]$, $\mathcal{M}=X_{0}[\overrightarrow{\underleftarrow{y_{1},y_{2}}}]$, where $[\overrightarrow{\underleftarrow{u,v}}]$ stands for the contour traversed from $u$ to $v$ along the upper edge of the slit $[u,v]$ and then back to $u$ along the lower edge of the slit. In the following lemma we provide exact characterization for the smooth and closed contours $\mathcal{L}$, $\mathcal{M}$ respectively.
\begin{lemma}\label{SQ} The algebraic function $X(y)$, $y\in[y_{1},y_{2}]$ lies on a closed contour $\mathcal{M}$, which is symmetric with respect to the real line and defined by
\begin{displaymath}
\begin{array}{l}
|x|^{2}=m(Re(x)),\,m(\delta)=\frac{\widehat{c}(\zeta(\delta))}{\widehat{a}(\zeta(\delta))},\,|x|^{2}\leq\frac{\widehat{c}(y_{2})}{\widehat{a}(y_{2})},
\end{array}
\end{displaymath}
where, $\zeta(\delta)=\frac{k_{2}(\delta)-\sqrt{k_{2}^{2}(\delta)-4k_{3}(\delta)k_{1}(\delta)}}{2k_{1}(\delta)}$, and
\begin{displaymath}
\begin{array}{rl}
k_{1}(\delta):=&a_{1}\widehat{a}_{2}(\lambda_{2}+2L(1-\delta))-(\lambda_{2}+L(1-2\delta))(1-\widehat{a}_{1}a_{2}),\\
k_{2}(\delta):=&2\delta[a_{1}\widehat{a}_{2}(\lambda_{1}+L)+\widehat{a}_{1}a_{2}(\lambda_{1}+2L)-\lambda_{1}-L]+\lambda_{1}+\lambda_{2}+L+a_{1}\widehat{a}_{2}(1-2\lambda_{1})\\
&+\widehat{a}_{1}a_{2}(1-\lambda_{1}-2(\lambda_{2}+L)),\\
k_{3}(\delta):=&-[\lambda_{2}(\widehat{a}_{1}a_{2}+a_{1}\widehat{a}_{2})+\widehat{a}_{1}a_{2}(1+(\lambda_{1}+L)(1+2\delta))].
\end{array}
\end{displaymath}
Moreover, $\beta_{0}:=\sqrt{\frac{\widehat{c}(y_{2})}{\widehat{a}(y_{2})}}$, $\beta_{1}:=-\sqrt{\frac{\widehat{c}(y_{1})}{\widehat{a}(y_{1})}}$ are the extreme right and left points of $\mathcal{M}$, respectively. Similarly, $Y(x)$, $x\in[x_{1},x_{2}]$ lies on a closed contour $\mathcal{L}$. Its exact representation is derived as for $\mathcal{M}$, and further details are omitted.
\end{lemma}
\begin{proof} See Appendix \ref{a1}.
\end{proof}
\subsection{The boundary value problems}\label{bound}
Here, we proceed with the derivation of the probability generating function of the joint stationary queue length distribution at relays. The analysis is distinguished in two cases according to the values of the parameters.
\subsubsection{A Dirichlet boundary value problem}\label{dir}
Let $\frac{a_{1}\tilde{a}_{2}}{a_{1}p_{R_{1},\{R_{1}\}}^{D}}+\frac{a_{2}\tilde{a}_{1}}{a_{2}p_{R_{2},\{R_{2}\}}^{D}}=1$. It can be easily seen that
$A(x,y)=\frac{d_{1}}{\alpha_{1}\widehat{\alpha}_{2}}B(x,y)\Leftrightarrow A(x,y)=\frac{\alpha_{2}\widehat{\alpha}_{1}}{d_{2}}B(x,y)$.

Therefore, for $y\in \mathcal{D}_{y}=\{y\in\mathcal{C}:|y|\leq1,|X_{0}(y)|\leq1\}$,
\begin{equation}
\begin{array}{l}
\alpha_{2}\widehat{\alpha}_{1}H(X_{0}(y),0)+d_{2}H(0,y)+\frac{\alpha_{2}\widehat{\alpha}_{1}(1-\rho)C(X_{0}(y),y)}{A(X_{0}(y),y)}=0.
\end{array}
\label{cons}
\end{equation}
Both $H(X_{0}(y),0)$, and $H(0,y)$, where $y\in \mathcal{D}_{y}-[y_{1},y_{2}]$, are analytic functions. Using analytic continuation arguments we consider (\ref{cons}) for $x\in\mathcal{M}$
\begin{equation}
\begin{array}{c}
\alpha_{2}\widehat{\alpha}_{1}H(x,0)+d_{2}H(0,Y_{0}(x))+\frac{\alpha_{2}\widehat{\alpha}_{1}(1-\rho)C(x,Y_{0}(x))}{A(x,Y_{0}(x))}=0.
\end{array}
\label{con2}
\end{equation}
By noticing that $H(0,Y_{0}(x))$ is real for $x\in\mathcal{M}$, i.e., $Y_{0}(x)\in[y_{1},y_{2}]$, we have
\begin{equation}
\begin{array}{c}
Re(iH(x,0))=w(x),\,x\in\mathcal{M},
\end{array}
\label{p1}
\end{equation}
where $w(x):=Re\left(-i\frac{C(x,Y_{0}(x))}{A(x,Y_{0}(x))}\right)(1-\rho)$. Clearly, some technical requirements are needed to be everything well defined. In particular, we have to investigate the possible poles of $H(x,0)$, $x\in S:=G_{\mathcal{M}}\cap\bar{D}_{x}^{c}$, where $G_{\mathcal{U}}$ be the interior domain bounded by $\mathcal{U}$, and $D_{x}=\{x:|x|<1\}$, $\bar{D}_{x}=\{x:|x|\leq1\}$, $\bar{D}_{x}^{c}=\{x:|x|>1\}$. This is equivalent with the investigation of the zeros of $A(x,Y_{0}(x))$, $x\in S_{x}$.
Moreover, in order to solve (\ref{p1}) we first transform the problem from $\mathcal{M}$ to the unit circle; see Appendix \ref{conf} for details. Let the conformal mapping, $z=\gamma(x):G_{\mathcal{M}}\to G_{\mathcal{C}}$, and its inverse $x=\gamma_{0}(z):G_{\mathcal{C}}\to G_{\mathcal{M}}$.

By applying the transformation, the problem is reduced to the determination of function $\tilde{T}(z)=H(\gamma_{0}(z),0)$ regular for $z\in G_\mathcal{C}$, continuous for $z\in\mathcal{C}\cup G_\mathcal{C}$ such that, $Re(i\tilde{T}(z))=w(\gamma_{0}(z))$, $z\in\mathcal{C}$. The solution of the Dirichlet problem with boundary condition (\ref{p1}) is:
\begin{equation}
\begin{array}{c}
H(x,0)=-\frac{1-\rho}{2\pi}\int_{|t|=1}f(t)\frac{t+\gamma(x)}{t-\gamma(x)}\frac{dt}{t}+C,\,x\in\mathcal{M},
\end{array}
\label{sol1}
\end{equation}
where $f(t)=Re\left(-i\frac{C(\gamma_{0}(t),Y_{0}(\gamma_{0}(t)))}{A(\gamma_{0}(t),Y_{0}(\gamma_{0}(t)))}\right)$, $C$ a constant to be defined by setting $x=0\in G_{\mathcal{M}}$ in (\ref{sol1}) and using the fact that $H(0,0)=1-\rho$, $\gamma(0)=0$. In case $H(x,0)$ has a pole, say $\bar{x}$, we still have a Dirichlet problem for the function $(x-\bar{x})H(x,0)$.

Setting $t=e^{i\phi}$, $\gamma_{0}(e^{i\phi})=\rho(\psi(\phi))e^{i\psi(\phi)}$, we obtain after some algebra, $f(e^{i\phi})=\frac{d_{1}\alpha_{2}\sin(\psi(\phi))(1-Y_{0}(\gamma_{0}(e^{i\phi}))^{-1})}{\rho(\psi(\phi))k(\phi)}$,

which is an odd function of $\phi$, and
\begin{displaymath}
\begin{array}{c}
k(\phi)=\left[\alpha_{2}\widehat{\alpha}_{1}(1-Y_{0}^{-1}(\gamma_{0}(e^{i\phi})))+d_{1}(1-\frac{\cos(\psi(\phi))}{\rho(\psi(\phi))})\right]^{2}+\left(d_{1}\frac{\sin(\psi(\phi))}{\rho(\psi(\phi))}\right)^{2}.
\end{array}
\end{displaymath}
Thus, $C=1-\rho$. Substituting in (\ref{sol1}) we obtain after simple calculations an integral representation of $H(x,0)$ on a real interval for $x\in G_{\mathcal{M}}$, i.e.,
\begin{equation}
\begin{array}{c}
H(x,0)=(1-\rho)\left\lbrace 1+\frac{2\gamma(x)i}{\pi}\int_{0}^{\pi}\frac{f(e^{i\phi})\sin(\phi)d\phi}{1-2\gamma(x)\cos(\phi)-\gamma(x)^{2}}\right\rbrace.
\end{array}
\label{soll}
\end{equation}

Similarly, we can determine $H(0,y)$ by solving another Dirichlet boundary value problem on the closed contour $\mathcal{L}$. Then, using the fundamental functional equation (\ref{we}), we uniquely obtain $H(x,y)$.
\subsubsection{A homogeneous Riemann-Hilbert boundary value problem}\label{rh}
In case $\frac{a_{1}\tilde{a}_{2}}{a_{1}^{*}p_{R_{1},\{R_{1}\}}^{D*}}+\frac{a_{2}\tilde{a}_{1}}{a_{2}^{*}p_{R_{2},\{R_{2}\}}^{D*}}\neq1$,
consider the following transformation:
\begin{displaymath}

G(x):=H(x,0)+\frac{\alpha_{1}p_{R_{1},\{R_{1}\}}^{D} d_{2}H(0,0)}{d_{1}d_{2}-\alpha_{1}\widehat{\alpha}_{2}\alpha_{2}\widehat{\alpha}_{1}},
L(y):=H(0,y)+\frac{\alpha_{2}p_{R_{2},\{R_{2}\}}^{D} d_{1}H(0,0)}{d_{1}d_{2}-\alpha_{1}\widehat{\alpha}_{2}\alpha_{2}\widehat{\alpha}_{1}}.

\end{displaymath}
Then, for $y\in \mathcal{D}_{y}$,
\begin{equation}
\begin{array}{c}
A(X_{0}(y),y)G(X_{0}(y))=-B(X_{0}(y),y)L(y).
\end{array}
\label{po31}
\end{equation}
Using similar arguments as in previous subsection, we have for $x\in\mathcal{M}$
\begin{equation}
\begin{array}{c}
A(x,Y_{0}(x))G(x)=-B(x,Y_{0}(x))L(Y_{0}(x)).
\end{array}
\label{za1}
\end{equation}
Clearly, $G(x)$ is holomorphic in $D_{x}$, continuous in $\bar{D}_{x}$, but it might have poles in $S_{x}$, based on the values of the system parameters. These poles, if exist, coincide with the zeros of $A(x,Y_{0}(x))$ in $S_{x}$. For $y\in[y_{1},y_{2}]$, let $X_{0}(y)=x\in\mathcal{M}$ and realize that $Y_{0}(X_{0}(y))=y$ so that $y=Y_{0}(x)$. Taking into account the possible poles of $G(x)$, and noticing that $L(Y_{0}(x))$ is real for $x\in\mathcal{M}$, since $Y_{0}(x)\in[y_{1},y_{2}]$, we have
\begin{equation}
\begin{array}{c}
Re[iU(x)\tilde{G}(x)]=0,\,x\in\mathcal{M},\\
U(x)=\frac{A(x,Y_{0}(x))}{(x-\bar{x})^{r}B(x,Y_{0}(x))},\,\tilde{G}(x)=(x-\bar{x})^{r}G(x),
\end{array}
\label{df3}
\end{equation}
where $r=0,1$, whether $\bar{x}$ is zero or not of $A(x,Y_{0}(x))$ in $S_{x}$. Thus, $\tilde{G}(x)$ is regular for $x\in G_{\mathcal{M}}$, continuous for $x\in\mathcal{M}\cup G_{\mathcal{M}}$, and $U(x)$ is a non-vanishing function on $\mathcal{M}$. By conformally transform the problem (\ref{df3}) from $\mathcal{M}$ to the unit circle, using $z=\gamma(x):G_{\mathcal{M}}\to G_{\mathcal{C}}$, and its inverse given by $x=\gamma_{0}(z):G_{\mathcal{C}}\to G_{\mathcal{M}}$, the problem in (\ref{df3}) is reduced to the following: find a function $F(z):=\tilde{G}(\gamma_{0}(z))$, regular in $G_{\mathcal{C}}$, continuous in $G_{\mathcal{C}}\cup\mathcal{C}$ such that, $Re[iU(\gamma_{0}(z))F(z)]=0,\,z\in\mathcal{C}$.

We now need additional information in order to derive a solution for the above problem, i.e., we must determine the index $\chi=\frac{-1}{\pi}[arg\{U(x)\}]_{x\in \mathcal{M}}$, where $[arg\{U(x)\}]_{x\in \mathcal{M}}$, denotes the variation of the argument of the function $U(x)$ as $x$ moves along the closed contour $\mathcal{M}$ in the positive direction, provided that $U(x)\neq0$, $x\in\mathcal{M}$. Following the methodology in \cite{fay}, the value of the index is closely related to the stability conditions. The following lemma provides the necessary information.
\begin{lemma}\begin{enumerate}
\item If $\lambda_{2}<\alpha_{2}\widehat{\alpha}_{1}$, then $\chi=0$ is equivalent to
\begin{displaymath}
\begin{array}{l}
\frac{d A(x,Y_{0}(x))}{dx}|_{x=1}<0\Leftrightarrow\lambda_{1}<\alpha_{1}p_{R_{1},\{R_{1}\}}^{D}+\frac{d_{1}\lambda_{2}}{\alpha_{2}\widehat{\alpha}_{1}},\\ \frac{d B(X_{0}(y),y)}{dy}|_{y=1}<0\Leftrightarrow\lambda_{2}<\alpha_{2}p_{R_{2},\{R_{2}\}}^{D}+\frac{d_{2}\lambda_{1}}{\alpha_{1}\widehat{\alpha}_{2}}.
\end{array}
\end{displaymath}
\item If $\lambda_{2}\geq \alpha_{2}\widehat{\alpha}_{1}$, $\chi=0$ is equivalent to
\begin{displaymath}
\begin{array}{c}
\frac{d B(X_{0}(y),y)}{dy}|_{y=1}<0\Leftrightarrow \lambda_{2}<\alpha_{2}p_{R_{2},\{R_{2}\}}^{D}+\frac{d_{2}\lambda_{1}}{\alpha_{1}\widehat{\alpha}_{2}}.
\end{array}
\end{displaymath}
\end{enumerate}
\end{lemma}
Thus, under stability conditions the problem defined in (\ref{df3}) has a unique solution for $x\in G_{\mathcal{M}}$ given by,
\begin{equation}

H(x,0)=K(x-\bar{x})^{-r}\exp\left\lbrace\frac{1}{2i\pi}\int_{|t|=1}\frac{\log J(t)}{t-\gamma(x)}dt \right\rbrace-\frac{\alpha_{1}p_{R_{1},\{R_{1}\}}^{D} d_{2}H(0,0)}{d_{1}d_{2}-\alpha_{1}\widehat{\alpha}_{2}\alpha_{2}\widehat{\alpha}_{1}},

\label{sool1}
\end{equation}
where $K$ is a constant and $J(t)=\frac{\overline{U_{1}(t)}}{U_{1}(t)}$, $U_{1}(t)=U(\gamma_{0}(t))$, $|t|=1$. Setting $x=0$ in (\ref{sool1}) we derive a relation between $K$ and $H(0,0)$. Then, for $x=1\in G_{\mathcal{M}}$, and using the first in (\ref{con}) we can obtain $K$ and $H(0,0)$. Substituting back in (\ref{sool1}) we obtain for $x\in G_{\mathcal{M}}$,
\begin{equation}
\begin{array}{l}
H(x,0)=\frac{\lambda_{1}d_{2}+\alpha_{1}\widehat{\alpha}_{2}(\bar{t}_{1}\bar{t}_{2}\alpha_{2}p_{R_{2},\{R_{2}\}}^{D}-\lambda_{2})}{(\bar{x}-1)^{r}\bar{t}_{1}\bar{t}_{2}(\alpha_{1}\widehat{\alpha}_{2}\alpha_{2}\widehat{\alpha}_{1}-d_{1}d_{2})} \times\\\left((\bar{x}-x)^{r}\exp\left[\frac{\gamma(x)-\gamma(1)}{2i\pi}\int_{|t|=1}\frac{\log J(t)}{(t-\gamma(x))(t-\gamma(1))}dt\right]+\frac{\alpha_{1}p_{R_{1},\{R_{1}\}}^{D} d_{2}\bar{x}^{r}}{\alpha_{1}\widehat{\alpha}_{2}\alpha_{2}p_{R_{2},\{R_{2}\}}^{D}}\exp\left[\frac{-\gamma(1)}{2i\pi}\int_{|t|=1}\frac{\log J(t)}{t(t-\gamma(1))}dt \right]\right).
\end{array}
\label{fin}
\end{equation}
The other unknown function $H(0,y)$ is determined similarly, by solving another Riemann-Hilbert boundary value problem on the closed contour $\mathcal{L}$. Then, using the fundamental functional equation (\ref{we}) we uniquely obtain $H(x,y)$. After that, we can derive several performance metrics. In particular, using Little's law, the expected delay at each aggregator is
\begin{displaymath}
\begin{array}{c}
D_{1}=\frac{\lambda_{1}+d_{1}H_{1}(1,0)}{\lambda_{1}a_{1}\widehat{a}_{2}},\,D_{2}=\frac{\lambda_{2}+d_{2}H_{2}(0,1)}{\lambda_{2}a_{2}\widehat{a}_{1}},
\end{array}
\end{displaymath}
where $H_{1}(1,0):=\frac{d}{dx}H(x,0)|_{x=1}$, and $H_{2}(0,1):=\frac{d}{dy}H(0,y)|_{y=1}$.

\section{Explicit bounds for the mean delay in the symmetrical system} \label{sec:delaysym}
In the following we consider the symmetrical model and obtain explicit bounds for the average delay. In particular, for $k=1,2$ let $\alpha_{k}=\alpha$, $t_{k}=t$, $p_{k,\{k\}}^{R_{k}}=s_{1}^{R}$, $p_{k,\{k\}}^{D}=s_{1}^{D}$, $p_{k,\{1,2\}}^{R_{k}}=s_{2}^{R}$, $p_{k,\{1,2\}}^{D}=s_{2}^{D}$, $p_{1,2,\{1,2\}}^{D}=s_{0}^{D}$, $p_{R_k,\{R_k\}}^{D}=r_{1}^{D}$, $p_{R_k,\{R_1,R_{2}\}}^{D}=r_{2}^{D}$, $p_{R_1,R_2,\{R_1,R_2\}}^{D}=r_{0}^{D}$, and as a result $d_{k}=d:=\alpha^{2}(r_{2}^{D}-r_{1}^{D})$. Note that due to the symmetry, $H(0,1)=H(1,0)$. Using the fact that $H(1,1)=1$, and setting in (\ref{we}) $y=1$, dividing with $x-1$, and taking the limit $x\to 1$, by using the L'Hospital rule,
\begin{equation}
\begin{array}{rl}
\alpha[\widehat{\alpha}+\alpha r_{0}^{D}]-\lambda=&(d+\alpha\widehat{\alpha}+2\alpha^{2}r_{0}^{D})H(1,0)-(d+\alpha^{2}r_{0}^{D})H(0,0),
\end{array}
\label{xb}
\end{equation}
where now due to the symmetry $\lambda_{k}=\lambda$, $k=1,2,$ $\widehat{\alpha}=\bar{\alpha}r_{1}^{D}+\alpha r_{2}^{D}$, and $\lambda=t^{2}s_{2}^{R}(\bar{s}_{0}^{D}+s_{2}^{D})+t\bar{t}\bar{s}_{1}^{D}s_{1}^{R}$.

Note that $\alpha[\widehat{\alpha}+\alpha r_{0}^{D}]>\lambda$ due to the stability condition.

Denote by $H_{1}(x,y)$, $H_{2}(x,y)$ the derivatives of $H(x,y)$ with respect to $x$, $y$. Due to the symmetry, $H_{1}(1,1)=H_{2}(1,1)$, $H_{1}(1,0)=H_{2}(0,1)$. Denote by $M_{k}=H_{1}(1,1)$, $k=1,2$. Differentiating (\ref{we}) with respect to $x$, setting $(x,y)=(1,1)$, and using (\ref{xb}) we obtain
\begin{equation}
\begin{array}{l}
M_{1}=\frac{\lambda\bar{\lambda}+(d+\alpha^{2}r_{0}^{D})H_{1}(1,0)}{\alpha[\widehat{\alpha}+\alpha r_{0}^{D}]-\lambda},
\end{array}
\label{d1}
\end{equation}
where $\bar{\lambda}=1-\lambda$. Now set $x=y$ in (\ref{we}) we obtain
\begin{equation}
\begin{array}{l}
2(\alpha[\widehat{\alpha}+\alpha r_{0}^{D}]-\lambda)\frac{d}{dx}H(x,x)|_{x=1}=2(1+2\lambda)(\alpha[\widehat{\alpha}+\alpha r_{0}^{D}]-\lambda)+\alpha^{2}r_{0}^{D}P(N_{1}>0,N_{2}>0)\\
+2H_{1}(1,0)(d+\alpha\widehat{\alpha}+2\alpha^{2}r_{0}^{D})-2\alpha(\widehat{\alpha}+\alpha r_{0}^{D})
+t^{2}\bar{s}_{0}^{D}(s_{2}^{R})^{2}+4\lambda(1-\alpha\widehat{\alpha}-\alpha^{2}r_{0}^{D}).
\end{array}
\label{d2}
\end{equation}
Using (\ref{d1}), (\ref{d2}), and realizing that due to the symmetry $\frac{d}{dx}H(x,x)|_{x=1}=2M_{1}$, we obtain
\begin{equation}
\begin{array}{l}
M_{1}=M_{2}=\frac{\lambda\bar{\lambda}(d+\alpha\widehat{a}+2\alpha^{2}r_{0}^{D})}{2\alpha(\alpha[\widehat{\alpha}+\alpha r_{0}^{D}]-\lambda)}
-\frac{(d+\alpha^{2}r_{0}^{D})[2\lambda(1+2\bar{\lambda})+t^{2}\bar{s}_{0}^{D}(s_{2}^{R})^{2}+\alpha^{2}r_{0}^{D}P(N_{1}>0,N_{2}>0)]}{2 \alpha(\alpha[\widehat{\alpha}+\alpha r_{0}^{D}]-\lambda)}.
\end{array}
\end{equation}
Using Little's law, the average delay at each aggregator is
\begin{equation}
\begin{array}{c}
D_{1}=D_{2}:=S-\phi,
\end{array}
\label{del}
\end{equation}
where $S=\frac{\lambda\bar{\lambda}(d+\alpha\widehat{a}+2\alpha^{2}r_{0}^{D})-(d+\alpha^{2}r_{0}^{D})[2\lambda(1+2\bar{\lambda})+t^{2}\bar{s}_{0}^{D}(s_{2}^{R})^{2}]}{2\lambda \alpha(\alpha[\widehat{\alpha}+\alpha r_{0}^{D}]-\lambda)}$ and $\phi=\frac{(d+\alpha^{2}r_{0}^{D})\alpha^{2}r_{0}^{D}P(N_{1}>0,N_{2}>0)}{2\lambda \alpha(\alpha[\widehat{\alpha}+\alpha r_{0}^{D}]-\lambda)}$.

Note that in case the destination can successfully receive at most one packet from the aggregators even if both of them transmit, i.e., $r_{0}^{D}=0$, the exact average queueing delay in a node is given by (\ref{del}) for $\phi=0$. Moreover, if $r_{0}^{D}=r_{2}^{D}=0$, then (\ref{del}) provides also the exact average delay for the case of a \textit{collision} channel. For the case of the general MPR model where $r_{0}^{D}\neq 0$, we can easily obtain bounds for the average delay at aggregators based on the sign of $\phi$. Since $Pr(N_1 > 0, N_2 > 0) > 0$, the sign of $\phi$ coincides with the
sign of $d+\alpha^{2}r_{0}^{D}$. To proceed, we distinguish the analysis in the following two cases:
\begin{enumerate}
\item If $d+\alpha^{2}r_{0}^{D}<0$, then $0\leq\phi\leq -\frac{\alpha^{2}r_{0}^{D}(d+\alpha^{2}r_{0}^{D})}{2\lambda \alpha(\alpha[\widehat{\alpha}+\alpha r_{0}^{D}]-\lambda)}$.
Thus, the upper and lower delay bound, say $D_{1}^{low}$, $D_{1}^{up}$ respectively are,
\begin{displaymath}
\begin{array}{lr}
D_{1}^{low}=S,&D_{1}^{up}=D_{1}^{low}-\frac{\alpha^{2}r_{0}^{D}(d+\alpha^{2}r_{0}^{D})}{2\lambda \alpha(\alpha[\widehat{\alpha}+\alpha r_{0}^{D}]-\lambda)}.
\end{array}
\end{displaymath}
\item If $d+\alpha^{2}r_{0}^{D}>0$, then $-\frac{\alpha^{2}r_{0}^{D}(d+\alpha^{2}r_{0}^{D})}{2\lambda \alpha(\alpha[\widehat{\alpha}+\alpha r_{0}^{D}]-\lambda)}\leq\phi\leq 0$.
In such a case,
\begin{displaymath}
\begin{array}{rl}
D_{1}^{up}=S,&
D_{1}^{low}=D_{1}^{up}-\frac{\alpha^{2}r_{0}^{D}(d+\alpha^{2}r_{0}^{D})}{2\lambda \alpha(\alpha[\widehat{\alpha}+\alpha r_{0}^{D}]-\lambda)}.
\end{array}
\end{displaymath}
\end{enumerate}

The following remark shows that the gap between the upper and the lower bound of the average delay has a closed form expression.

\begin{remark}
We observe that the gap between the lower and the upper average delay is given by $D_{1}^{up} - D_{1}^{low} =\left|\frac{\alpha^{2}r_{0}^{D}(d+\alpha^{2}r_{0}^{D})}{2\lambda \alpha(\alpha[\widehat{\alpha}+\alpha r_{0}^{D}]-\lambda)} \right|$. Note that when the arrival rate $\lambda$ increases, then the gap tends to zero. Furthermore, when $\lambda\to 0$, then the average delay is very close to zero thus, the lower bound on the delay is tight.
\end{remark}
\begin{remark}
In the case of a collision channel model assumption, which is equivalent to $r_{0}^{D}=r_{2}^{D}=0$, we have that $D_{1}^{up} = D_{1}^{low}$. This is also true for the case where $r_{0}^{D}=0$, which means that the destination can successfully receive at most one packet from the aggregators when both of them transmit. This is known as the capture effect which is common in LoRaWAN \cite{LoRaCapture2017}.
\end{remark}

\section{Numerical Results} \label{sec:Res}
In this section, we provide numerical results to evaluate the presented theoretical performance analysis. For exposition convenience, we consider the case where all sensors have the same link characteristics and transmission probabilities.

We consider the following network setup, the distance between the sensors and the sink is $130m$, the distance between the sensors and the aggregator is $60m$, and the distance between the aggregators and the destination is $80m$. The path loss exponent is assumed $4$, and the transmission power for each sensor is $1mW$ and for each aggregator is $10mW$. We assume Rayleigh block fading channel model.
The values for the SINR threshold are considered in this section are $0.2$, $0.5$, $1.2$, and $2$. In addition, regarding the transmission probabilities, for the sensors we consider the values $t=0.1$, $t=0.2$, and for the aggregators $\alpha=0.8$.

\subsection{Stability Region}

Here we present the closure of the stability region for the queues at the aggregators for all the possible random access probabilities $\alpha$. We consider two cases for the SINR threshold, the case where SINR$=0.5$, which is depicted in Fig. \ref{fig:stability_low}. \textit{We observe that this region is a convex set, thus, the system performs better than time sharing schemes and this is an indication of strong MPR capabilities at the destination}. The stability region for SINR$=1.2$ is depicted in Fig. \ref{fig:stability_high}. \textit{In this case, a time sharing scheme has higher performance since the destination has weak MPR capabilities}.

\begin{figure}[!ht]
	\centering
	\subfigure[SINR$=0.5$]{
		\includegraphics[scale=0.38]{stability_low_SINR}
		\label{fig:stability_low}
	}
	\subfigure[SINR$=1.2$]{
		\includegraphics[scale=0.38]{stability_high_SINR}
		\label{fig:stability_high}
	}\vspace{-0.1in}
	\caption{The stability region for the queues at the aggregators.}
\end{figure}

\subsection{Network Throughput}

In this subsection, we present the performance in terms of throughput per sensor and the aggregate throughput of the considered system. We consider four cases for the SNR/SINR threshold,
two cases for the transmission probabilities of the sensors, $t=0.1$ and $t=0.2$, and for the aggregators $\alpha=0.8$. We also plot the throughput of the IoT network without the presence of the aggregators. The case where SINR$ < 1$, is depicted in Fig. \ref{fig:thlow}, in this case, we have strong MPR capabilities at the receivers, thus it is more likely to have more concurrent successful transmissions. Furthermore, we observe that the performance of the IoT network without the assistance of the aggregators can be sufficient for the case with low number of sensors.

\begin{figure}[!ht]
	\centering
	\subfigure[Throughput per sensor]{
		\includegraphics[scale=0.45]{thr-low}
		\label{fig:thrlow}
	}
	\subfigure[Aggregate Throughput]{
		\includegraphics[scale=0.45]{athr-low}
		\label{fig:athrlow}
	}\vspace{-0.1in}
	\caption{Throughput performance of the IoT network for SINR$ < 1$.} \label{fig:thlow}
\end{figure}

The case where SINR$ > 1$, is depicted in Fig. \ref{fig:thhigh}. In this case, the benefits of the aggregators deployment are prominent. The throughput of the IoT network without the aggregators is almost zero. \textit{The presence of the aggregators that deploy network-level cooperation provides significant gains in the throughput performance of the network.}
\begin{figure}[!ht]
	\centering
	\subfigure[Throughput per sensor]{
		\includegraphics[scale=0.45]{thr-high}
		\label{fig:thrhigh}
	}
	\subfigure[Aggregate Throughput]{
		\includegraphics[scale=0.45]{athr-high}
		\label{fig:athrhigh}
	}\vspace{-0.1in}
	\caption{Throughput performance of the IoT network for SINR$ > 1$.} \label{fig:thhigh}
\end{figure}

In Fig. \ref{fig:rth}, we present the percentage of the traffic that is being relayed. We observe that in almost all cases the majority of the traffic is relayed.
Furthermore, in Fig. \ref{fig:rthrhigh} it becomes apparent that almost all the traffic that ends to the common destination is relayed traffic.

\begin{figure}[!ht]
	\centering
	\subfigure[SINR$ < 1$]{
		\includegraphics[scale=0.45]{rthr-low}
		\label{fig:rthrlow}
	}
	\subfigure[SINR$ > 1$]{
		\includegraphics[scale=0.45]{rthr-high}
		\label{fig:rthrhigh}
	}\vspace{-0.1in}
	\caption{Percentage of relayed traffic per sensor.} \label{fig:rth}
\end{figure}

Table \ref{tab:stable_sensors} provides what is the number of sensors per aggregator that can cause its queue stable or unstable.
As expected, when SINR$<1$, we have a clear transition from a stable queue to unstable queue as the number of sensors increase. However, the more interesting case is when SINR$>1$. We observe that for low number of sensors, we start with a stable queue, then the queue becomes unstable, and finally becomes stable again. The reason behind this is that initially, the incoming traffic increases up to a point the queue becomes unstable. Then due to the increased interference, the incoming traffic is reduced significantly thus, the queues are stable.
The shape of the curves presented in Fig. \ref{fig:rth} can be explained by Table \ref{tab:stable_sensors}. These results provide useful guidelines on the number of sensors that can be supported by an aggregator in order to achieve a required network performance regarding throughput.

\begin{table}[h]
\begin{center}
	\renewcommand{\arraystretch}{0.7}
	\caption{State of the queue per aggregator and the number of sensors.}
\begin{tabular}{|c|c|c|}
\hline
 & Stable & Unstable \\\hline
$SINR=0.2, t=0.2$ & $\{1,...,6\}$ & $\{7,...,30\}$ \\\hline
$SINR=0.5, t=0.2$ & $\{1,...,4\}$ & $\{5,...,30\}$ \\\hline
$SINR=0.2, t=0.1$ & $\{1,...,13\}$ & $\{14,...,30\}$ \\\hline
$SINR=0.5, t=0.1$ & $\{1,...,9\}$ & $\{10,...,30\}$ \\\hline \hline
$SINR=1.2, t=0.2$ & $\{1,...,3\},\{19,...,30\} $ & $\{4,...,18\}$ \\\hline
$SINR=2, t=0.2$ & $\{1,2\}, \{14,...,30\}$ & $\{3,...,13\}$ \\\hline
$SINR=1.2, t=0.1$ & $\{1,...,7\}$ & $\{8,...,30\}$ \\\hline
$SINR=2, t=0.1$ & $\{1,...,6\}, \{28,...,30\}$ & $\{7,...,27\}$ \\\hline
\end{tabular}
\label{tab:stable_sensors}
\end{center}
\end{table}

\subsection{Delay Performance}
In Fig. \ref{fig:delay}, the lower and the upper bounds for the average delay per packet versus the arrival rate at the aggregator are depicted for the four cases of SINR threshold, $0.2$, $0,5$, $1.2$ and $2$ respectively. The values of arrival rates are obtained by the stability conditions and we can connect the value of arrival rate with the number of users and their transmission probability through expression \eqref{eq:lambdaju}. Furthermore, we observe that the average delay increases with the increase of the SINR threshold. It is also important to notice that the upper and the lower bounds are very close, as an indication that the obtained bounds are tight.

\begin{figure}[!ht]
\centering
\includegraphics[scale=0.5]{delay}\vspace{-0.1in}
\caption{Average delay performance versus average arrival rate for various values of the SINR threshold for the case $\alpha=0.8$.}
\label{fig:delay}
\end{figure}

\section{Conclusions and Future Directions} \label{sec:conclusions}

In this work, we considered a random access IoT wireless network assisted by two aggregators providing support for data collection by applying network level cooperation. We characterized the throughput performance of the IoT network and obtained the stability conditions for the queues at the aggregators, which guarantee finite queueing delay. Furthermore, by applying the theory of boundary value problems we provided a detailed analysis for the delay. We showed that the benefits for the deployment of aggregators are more profound for higher SINR threshold and when the MPR capability of the network is weak. Our analytical results provide useful design guidelines for deploying aggregators in random access wireless IoT networks.

The suggested framework in this work can be extended to capture the case of cognitive aggregators that are adapting according to the incoming traffic from the sensors and also the channel conditions.
Considering adaptation of the operation of the aggregator according to the incoming traffic, can provide also performance improvements related to power consumption and energy efficiency, which is also another important factor for the design of IoT networks. In addition, it will be of further interest to consider delay-aware operation protocols for the IoT network and using successive interference cancellation to mitigate interference among the aggregators at the destination. The case of multiple aggregators in the same coverage area is an interesting future direction, our approach can also be applied but with some modifications.

\appendices

\section{Proof of Theorem \ref{thm1}} \label{sec:stability_proof}
To determine the stability region of the aggregators in our network, we apply the stochastic dominance technique \cite{Rao_TIT1988}. More specifically, we construct hypothetical dominant systems, in which an aggregator transmits dummy packets for the packet queue that is empty, while for the non-empty queue it transmits according to its traffic. Under this approach, we consider the $\mathcal{D}_{1}$, and $\mathcal{D}_{2}$-dominant systems. In the $\mathcal{D}_{k}$ dominant system, whenever the queue of aggregator $k$, $k=1,2$ empties, it transmits a dummy packet instead.

Thus, in $\mathcal{D}_{1}$, the first aggregator never empties, and hence, the second aggregator sees a constant service rate, while the service rate of aggregator $1$ depends on the state of aggregator $2$, i.e., empty or not. We proceed with the queue at aggregator $1$, its service rate is given by
\begin{equation} \label{eq:mu1}
\begin{array}{l}
\mu_1=Pr(N_2 \neq 0) \left[\alpha_1\bar{\alpha}_{2}p_{R_{1},\{R_{1}\}}^{D}+\alpha_1 \alpha_2 p_{R_1/R_1,R_2}^{D}\right] + Pr(N_2 = 0) \alpha_{1} p_{R_{1},\{R_{1}\}}^{D}.
\end{array}
\end{equation}
The service rate of the second aggregator is given by
\begin{equation}
\begin{array}{c}
\mu_2=\alpha_2\bar{\alpha}_{1}p_{R_{2},\{R_{2}\}}^{D}+\alpha_1 \alpha_2 p_{R_2/R_1,R_2}^{D}.
\end{array}
\end{equation}
By applying Loyne's criterion, the second node is stable if and only if the average arrival rate is less that the average service rate, $\lambda_{2} < \alpha_2\bar{\alpha}_{1}p_{R_{2},\{R_{2}\}}^{D}+\alpha_1 \alpha_2 p_{R_2/R_1,R_2}^{D}$. We can obtain the probability that the second node is empty and is given by $Pr(N_2 = 0)=1-\frac{\lambda_{2}}{\mu_2}$. After replacing $Pr(N_2 = 0)$ into \eqref{eq:mu1}, and applying Loynes criterion we can obtain the stability condition for the first node. Then, we have the stability region $\mathcal{R}_1$ given by \eqref{eq:stabCond} for $i=1$.
Similarly, we can obtain the stability region for the second dominant system $\mathcal{R}_2$. For a detailed treatment of dominant systems please refer to \cite{Rao_TIT1988}.

An important observation made in \cite{Rao_TIT1988} is that the stability conditions obtained by the stochastic dominance technique are not only sufficient but also necessary for the stability of the original system. The \emph{indistinguishability} argument \cite{Rao_TIT1988} applies here as well. Based on the construction of the dominant system, we can see that the queue sizes in the dominant system are always greater than those in the original system, provided they are both initialized to the same value and the arrivals are identical in both systems. Therefore, given $\lambda_{2}<\mu_{2}$, if for some $\lambda_{1}$, the queue at the first user is stable in the dominant system, then the corresponding queue in the original system must be stable. Conversely, if for some $\lambda_{1}$ in the dominant system, the queue at the first node saturates, then it will not transmit dummy packets, and as long as the first user has a packet to transmit, the behavior of the dominant system is identical to that of the original system since dummy packet transmissions are eliminated as we approach the stability boundary. Therefore, the original and the dominant system are indistinguishable at the boundary points.

\section{Proof of Lemma \ref{LEM}}\label{a0}
It is readily seen that $R(x,y)=\frac{xy-\Psi(x,y)}{xy}$, where $\Psi(x,y)=L(x,y)[xy+y(1-x)\alpha_{1}\widehat{\alpha}_{2}+x(1-y)\alpha_{2}\widehat{\alpha}_{1}]$, where for $|x|\leq1$, $|y|\leq1$, $\Psi(x,y)$ is a generating function of a proper probability distribution. Now, for $|y|=1$, $y\neq1$ and $|x|=1$ it is clear that $|\Psi(x,y)|<1=|xy|$. Thus, a direct application of Rouch\'e's theorem states that, $xy-\Psi(x,y)$ has exactly one zero inside the unit circle. Therefore, $R(x,y)=0$ has exactly one root $x=X_{0}(y)$, such that $|x|<1$. For $y=1$, $R(x,1)=0$ implies

$(1-x)\left(\lambda_{1}+\lambda_{1}\frac{\alpha_{1}\widehat{\alpha}_{2}(1-x)}{x}-\frac{\alpha_{1}\widehat{\alpha}_{2}}{x} \right)=0.$

Therefore, for $y=1$, and since $\lambda_{1}<\alpha_{1}\widehat{\alpha}_{2}$, the only root of $R(x,1)=0$ for $|x|\leq1$, is $x=1$.\hfill$\square$
\section{Proof of Lemma \ref{SQ}}\label{a1}
For $y\in[y_{1},y_{2}]$, $D_{y}(y)$ is negative, so $X_{0}(y)$, $X_{1}(y)$ are complex conjugates. Therefore, $|X(y)|^{2}=\frac{\widehat{c}(y)}{\widehat{a}(y)}=g(y)$. Clearly, $g(y)$ is an increasing function for $y\in[0,1]$ and thus, $|X(y)|^{2}\leq g(y_{2})=\beta_{0}$. Using simple algebraic considerations we can prove that, $X_{0}(y_{1}):=\beta_{1}=-g(y_{1})$ is the extreme left point of $\mathcal{M}$. Finally, $\zeta(\delta)$ is derived by solving $Re(X(y))=-\widehat{b}(y)/2\widehat{a}(y)$ for $y$ with $\delta = Re(X(y))$, and taking the solution such that $y\in[0,1]$.\hfill$\square$
\section{Construction of the conformal mappings}\label{conf}
Our goal in order to to obtain expressions for the basic performance metrics is to construct the mapping $\gamma_{0}(z)$. To proceed, we need a representation of $\mathcal{M}$ in polar coordinates, i.e., $\mathcal{M}=\{x:x=\rho(\phi)\exp(i\phi),\phi\in[0,2\pi]\}.$

In the following, we summarize the basic steps: Since $0\in G_{\mathcal{M}}$, for each $x\in\mathcal{M}$, a relation between its absolute value and its real part is given by $|x|^{2}=m(Re(x))$ (see Lemma \ref{SQ}). Given the angle $\phi$ of some point on $\mathcal{M}$, the real part of this point, say $\delta(\phi)$, is the zero of $\delta-\cos(\phi)\sqrt{m(\delta)}$, $\phi\in[0,2\pi].$ Since $\mathcal{M}$ is a smooth, egg-shaped contour, the solution is unique. Clearly, $\rho(\phi)=\frac{\delta(\phi)}{\cos(\phi)}$, and the parametrization of $\mathcal{M}$ in polar coordinates is fully specified.

 Then, the mapping from $z\in G_{\mathcal{C}}$ to $x\in G_{\mathcal{M}}$, where $z = e^{i\phi}$ and $x= \rho(\psi(\phi))e^{i\psi(\phi)}$, satisfying $\gamma_{0}(0)=0$, $\gamma_{0}(z)=\overline{\gamma_{0}(\bar{z})}$ is uniquely determined by (see \cite{coh}, Section I.4.4),
\begin{equation}
\begin{array}{rl}
\gamma_{0}(z)=&z\exp[\frac{1}{2\pi}\int_{0}^{2\pi}\log\{\rho(\psi(\omega))\}\frac{e^{i\omega}+z}{e^{i\omega}-z}d\omega],\,|z|<1,\\
\psi(\phi)=&\phi-\int_{0}^{2\pi}\log\{\rho(\psi(\omega))\}\cot(\frac{\omega-\phi}{2})d\omega,\,0\leq\phi\leq 2\pi,
\end{array}
\label{zx}
\end{equation}
i.e., $\psi(.)$ is uniquely determined as the solution of a Theodorsen integral equation with $\psi(\phi)=2\pi-\psi(2\pi-\phi)$. This integral equation has to be solved numerically by an iterative procedure. For the numerical evaluation of the integrals we split the interval $[0,2\pi]$ into $M$ parts of length $2\pi/M$, by taking $M$ points $\phi_{k}=\frac{2k\pi}{M}$, $k=0,1,...,M-1$. For the $M$ points given by their angles $\left\{\phi_{0},...,\phi_{M-1}\right\}$ we should solve the second in (\ref{zx}) to obtain the corresponding points $\left\{\psi(\phi_{0}),...,\psi(\phi_{M-1})\right\}$, iteratively from,
\begin{equation}

\psi_{0}(\phi_{k})=\phi_{k},
\psi_{n+1}(\phi_{k})=\phi_{k}-\frac{1}{2\pi}\int_{0}^{2\pi}\log\left\{\frac{\delta(\psi_{n}(\omega))}{\cos(\psi_{n}(\omega))}\right\}\cot\left(\frac{\omega-\phi_{k}}{2}\right)d\omega,

\label{uip}
\end{equation}
where $\lim_{n\to\infty}\psi_{n+1}(\phi)=\psi(\phi)$, and $\delta(\psi_{n}(\omega))$ is determined by, $\delta(\psi_{n}(\omega))=cos(\psi_{n}(\omega))\sqrt{m(\delta(\psi_{n}(\omega)))}$, using the Newton-Raphson root finding method. For each step, the integral in (\ref{uip}) is
numerically determined by again using the trapezium rule with $M$ parts of equal length
$2\pi/M$. For the iteration, we have used the following stopping criterion $\max_{k\in\left\{0,1,...,M-1\right\}}\left|\psi_{n+1}(\phi_{k})-\psi_{n}(\phi_{k})\right|<10^{-6}$

After obtaining $\psi(\phi)$ numerically, the values of the conformal mapping $\gamma_{0}(z)$, $\left|z\right|\leq 1$, can be calculated by applying the Plemelj-Sokhotski formula to the first in (\ref{zx}) for $0\leq\phi\leq 2\pi$,
\begin{displaymath}
\begin{array}{rl}
\gamma_{0}(e^{i\phi})=\frac{e^{i\psi(\phi)}\delta(\psi(\phi))}{\cos(\psi(\phi))}=\delta(\psi(\phi))[1+i \tan(\psi(\phi))].
\end{array}
\end{displaymath}
We can find $\gamma(1)$, $\gamma^{\prime}(1)$ by applying the Newton's method and solving $\gamma_{0}(z_{0})=1$, in $[0,1]$, i.e., $z_{0}$ is the zero in $[0,1]$ of $\gamma_{0}(z)=1$. Then, $\gamma(1)=z_{0}$. Moreover, using the first in (\ref{zx})
\begin{equation}
\begin{array}{rl}
\gamma^{\prime}(1)=&(\gamma_{0}^{\prime}(z_{0}))^{-1}=\left\lbrace\frac{1}{\gamma(1)}+\frac{1}{2\pi i}\int_{0}^{2\pi}\frac{\log\{\rho(\psi(\omega))\}2e^{i\omega}}{(e^{i\omega}-\gamma(1))^{2}}d\omega\right\rbrace^{-1},
\end{array}
\label{cv}
\end{equation}
which can be obtained numerically by using the trapezoidal rule for the integral on the right-hand side of (\ref{cv}).

\bibliographystyle{IEEEtran}
\bibliography{bibliography}
\end{document}

config: arxiv


TITLE: Diffgeo question on curves and properties of dot product
QUESTION [2 upvotes]: Let $\gamma$ be a curve such that $\gamma(a)=p$ and $\gamma(b)=q$ and let $u$ be a unit vector.
A) Prove $\gamma'(t) \cdot u$ $\leq$ $\|\gamma’(t)\|$
For this one I use the properties of dot product, noting that $u=\frac{\gamma’}{\|\gamma\|}$ but I always end up just showing that the quantities are equal.
B) Use path integrals to show ($\gamma(b)-\gamma(a)) \cdot u$ $\leq$ $\int_{a}^{b} \|\gamma'(t)\| dt$
For this one I observe that $u=\frac{(q-p)}{\|(q-p)\|}$ and sloppily reach a conclusion.
C) Prove the length of gamma from $p$ to $q$ is larger than or equal to the length of the straight line between $p$ and $q$.
This one is intuitively obvious, but for the proof I observe the straight line dist is just $\|p-q\|$ and the length of gamma is the arc length. I had myself convinced I had this part the other day, but now I’m confused.
Hopefully someone can make this all more clear for me. I apologize for the terrible (lack of) formatting
Thx
Nick

REPLY [2 votes]: Here are three hints.
A) Use the Cauchy-Schwarz inequality $|\langle v, w \rangle|\leq \|v\| \|w\|$. What do you choose for $v$ and $w$?
B) Use the fact that for two real functions $f$ and $g$ with $f\leq g$, one has
$$
   \int_a^b f(t)\,dt \leq  \int_a^b g(t)\,dt.
$$
What $f$ and $g$ do you choose?
C) You have now proven the inequality
$$(\gamma(b)-\gamma(a))\cdot u \leq \int_a^b \|\gamma'(t)\|\,dt$$
for every unit vector $u$. Note that the right hand side $\int_a^b \|\gamma'(t)\|\,dt$ is the length of $\gamma$. Choose an appropriate $u$ to show $\|p-q\|$ is smaller than the length of $\gamma$.


TITLE: Write statements using quantifiers
QUESTION [0 upvotes]: I want to know how to properly write below statements using quantifiers and simple relations ($\in, =, <, \leq, |$):

Between any 2 rational numbers exists number different than these 2.
Every natural number greater than 1 has at least one prime factor.

1st statement seems to be easy, I think it's right:
$$\forall x, y\in \mathbb{R} \exists c\in \mathbb{R}: x < c < y \lor y < c < x$$
but I have problems writing properly 2nd statement. I tried something like that, but I think it's wrong:
$$\forall n\in\mathbb{N}, n > 1 \exists p \in \mathbb{N} : p|n \land \forall x \in \mathbb{N} x < p : x\not| p \land p \not= 1$$
Also, I'd be interested how to write "definition" of NWD(a,b) using quantifiers and basics relations. Can you point me to some books/tutorials where I can excercise "rewriting" statemants in natural language to mathematics language?
Thanks for help

REPLY [1 votes]: Between any 2 rational numbers exists number different than these 2.

$\forall x, y\in \mathbb{Q}, \exists c\in \mathbb{Q}, x < c < y \lor y < c < x$
By the way this statement is false since it allows $x=y$.

Every natural number greater than 1 has at least one prime factor

$\forall n \in \mathbb{N}, n > 1, \exists p \in \mathbb{N}, p \mbox{ prime} \land p \mid n$
$p \mbox{ prime} \Leftrightarrow \forall q \in \mathbb{N}, q \not= 1, q \not= p, q \nmid p$
So you get
$\forall n \in \mathbb{N}, n > 1, \exists p \in \mathbb{N}, [\forall q \in \mathbb{N}, q \not= 1, q \not= p, q \nmid p] \land p \mid n$
Which you can also write
$\forall n \in \mathbb{N}, n > 1, \exists p \in \mathbb{N}, \forall q \in \mathbb{N}, q \not= 1, q \not= p, q \nmid p \land p \mid n$


TITLE: Bi-annihilator of a subspace of the dual of an infinite-dimensional vector space
QUESTION [3 upvotes]: Let $V$ be an infinite-dimensional vector space and $V^*$ its dual.
For a linear subspace $W\subset V$ define $W^ \circ\subset V^*$ as the subspace of linear forms on $V$ vanishing on $W$.
Dually, for $\Gamma\subset V^*$ define $\Gamma^\diamond \subset V$ as the set of vectors $v\in V$ such that $\gamma(v)=0$ for all linear forms $\gamma\in \Gamma$.
It is slightly surprising but not too difficult to show that that we  have for all subspaces $W\subset V$ the equality $(W^\circ) ^\diamond=W$.
But is it true that for all $\Gamma\subset V^*$ we have $(\Gamma^\diamond)^\circ=\Gamma$ ?
And is there a reference (article, book, lecture notes,...) where this problem is mentioned?

REPLY [5 votes]: No, $(\Gamma^\diamond)^\circ$ need not always equal $\Gamma$.  Let $\mathcal B$ be a basis for $V$, and let $\Gamma$ be the span of the 'dual' set $\{e_b \mathrel: b \in \mathcal B\}$, so $e_b(c)$ is the Iverson bracket $[b = c]$ for all $b, c \in \mathcal B$.  Then $\Gamma^\diamond$ is $0$, so $(\Gamma^\diamond)^\circ$ is all of $V^*$; but $\Gamma$ itself does not contain, for example, the element $\sum_{b \in \mathcal B} e_b$ of $V^*$.

REPLY [3 votes]: The equality is false in general.
Here is a counterexample: fix a basis $v_i, i\in I$ of $V$ and consider the set of coordinate linear forms $v^*_i, i\in I$.
These forms are linearly independant but never form a basis since $V$ is infinite-dimensional.
So complete these forms to a basis $(v^*_j), j\in J$ with $J\setminus I\neq\emptyset$.
Choose $l\in J\setminus I$ and put $J'=J\setminus \{l\}$
If you define $\Gamma \subset V^*$ as the vector space generated by the $v_j^*, j\in J'$, then $\Gamma^\diamond =0$ (since already the subspace of $V^*$ generated by the $v_i^*, i\in I$ kill all vectors in $V$) so that $\Gamma\subsetneq (\Gamma^\diamond)^\circ=\{0\}^\circ=V^*$ yielding the required counterexample.


TITLE: Uniformly most powerful test for two-sided hypothesis
QUESTION [0 upvotes]: I need some help:
Prove that a uniformly most powerful test for a level $\alpha\in(0,1)$ doesn't exist for the test $H_0:\mu=\mu_0$ versus $H_1:\mu\neq\mu_0$, while $\mu,\mu_0\in\mathbb{R}$.

REPLY [0 votes]: It's not true in general. so you can not prove it. according to
Uniformly Most Powerful Test for a Uniform Sample
UMP test for
$\left\{\begin{array}{cc} H_0: & \theta= \theta_0 \\
H_1: & \theta \neq \theta_0
\end{array}  \right.$
exist. ( $\{x_i \}_{i=1}^{n} \overset{i.i.d}{\sim} Uniform(0,\theta)$)


TITLE: the image of normal subgroups
QUESTION [2 upvotes]: I want to find an example of a homomorphism $f:G\to H$ such that $A$ is a normal subgroup of $G,$ but $f(A)$ is not so in $H.$
I know that if $f$ was onto $f(A )$ must be normal , but otherwise i want to find an example!

REPLY [4 votes]: Take $A=G$ as a subgroup of $H$ that is not normal. Then define $f:G\rightarrow H$ by $x\mapsto x$
Notice here that $A=G$ is a normal subgroup of itself.
$f(A)=f(G)=G$ is not normal in $H$

REPLY [2 votes]: Let $H < G$ which is not normal, and the injection $H \to G$. Then $A = H$ is normal in $H$, but not in $G$.


TITLE: Prove that $\sqrt{1+\pi^2}, \hspace{2pt}\pi-\sqrt\pi,\hspace{2pt} \pi^2+\pi+\sqrt{1+2\pi}$ are not algebraic
QUESTION [0 upvotes]: Considering the fact that $\pi$ is not an algebraic number, I need to prove these expressions are not algebraic :
$$\sqrt{1+\pi^2},\hspace{10pt} \pi-\sqrt\pi,\hspace{10pt} \pi^2+\pi+\sqrt{1+2\pi}$$
Thanks
Some explanations will be appreciated ! Thank you all for your time .

REPLY [2 votes]: Set $c=\pi^2+\pi+\sqrt{1+2\pi}$. Then
$$
\sqrt{1+2\pi}=c-\pi^2-\pi
$$
and
$$
1+2\pi=(c-\pi^2-\pi)^2
$$
and, rearranging terms, you get a polynomial with coefficients in $\mathbb{Q}[c]$ that is satisfied by $\pi$. This means that $\pi$ is algebraic over $\mathbb{Q}[c]$. But you know that $\pi$ is transcendental, so also $c$ must be transcendental, otherwise you'd get that
$$
[\mathbb{Q}[\pi]:\mathbb{Q}]=
[\mathbb{Q}[\pi]:\mathbb{Q}[c]][\mathbb{Q}[c]:\mathbb{Q}]
$$
would be finite.


TITLE: Is it possible to test if a computable number is rational or integer?
QUESTION [21 upvotes]: Is it possible to algorithmically test if a computable number is rational or integer?
In other words, would it be possible for a library that implements computable numbers to provide the functions isInteger or isRational?
I am guessing that it is not possible, and that this is somehow related to the fact that it is not possible to test if two numbers are equal, but I don't see how to prove it.
Edit: A computable number $x$ is given by a function $f_x(\epsilon)$ that can return a rational approximation of $x$ with precision $\epsilon$: $|x - f_x(\epsilon)| \leq \epsilon$, for any $\epsilon > 0$. Given such function, is it possible to test if $x \in \mathrm{Q}$ or $x \in \mathrm{Z}$?

REPLY [42 votes]: It is easy to get confused about what it means to "represent" or "implement" a real number. In fact, we are witnessing a discussion in the comments where the representation is contentious. So let me address this first.
How do we know that an implementation is correct?
The theory which explains how to represent things in a computer is realizability. The basic idea is that, given a set $X$, we pick a datatype $\tau$ and to every $x \in X$ a set of values of type $\tau$ which realize it. We write $v \vdash x \in X$ when $v$ is a value that realizes $x$. For example (I shall use Haskell for no good reason), a sensible implementation of $\mathbb{N}$ might be the datatype Integer where $v \vdash k \in \mathbb{N}$ when $v$ evaluates to the numeral $\overline{k}$ (thus in particular -42 does not represent a natural number, and neither does a diverging program). But some joker could walk by and suggest that we use Bool to represent natural numbers with $\mathtt{True} \vdash 42 \in \mathbb{N}$ and $\mathtt{False} \vdash n \in \mathbb{N}$ for $n \neq 42$. Why is this incorrect? We need a criterion.
In the case of "joker numbers" the easy observation is that addition cannot be implemented. Suppose I tell you I have two numbers, both represented by $\mathtt{False}$. Can you give a realizer for their sum? Well, that depends on whether the sum is 42, but you cannot tell. Since addition is an "essential part of what natural numbers are", this is unacceptable. In other words, implementation is not about sets, but about structures, i.e., we have to represent sets in such a way that it is possible to also implement the relevant structure. Let me stress this:

We implement structures, not bare sets. Therefore, we have to be able to implement the entire structure, together with operations and all the axioms, in order for the implementation to be correct.

If you do not abide by this principle, then you have to suggest an alternative mathematical criterion of correctness. I do not know of one.
Example: representation of natural numbers
For natural numbers the relevant structure is described by Peano axioms, and the crucial axiom that has to be implemented is induction (but also $0$, successor, $+$ and $\times$). We can compute, using realizability, what the implementation of induction does. It turns out to be a map (where nat is the yet unknown datatype which represents natural numbers)

induction : 'a -> (nat -> 'a -> 'a) -> 'nat -> 'a

satisfying induction x f zero = x and induction x f (succ n) = f n (induction x f n). All this comes out of realizability. We have a criterion: an implementation of natural numbers is correct when it allows an implementation of Peano axioms. A similar result would be obtained if we used the characterization of numbers as the initial algebra for the functor $X \mapsto 1 + X$.
Correct implementation of real numbers
Let us turn attention to the real numbers and the question at hand. The first question to ask is "what is the relevant structure of the real numbers?" The answer is: Archimedean Cauchy complete ordered field. This is the established meaning of "real numbers". You do not get to change it, it has been fixed by others for you (in our case the alternative Dedekind reals turn out to be isomorphic to the Cauchy reals, which we are considering here.) You cannot take away any part of it, you are not allowed to say "I do not care about implementing addition", or "I do not care about the order". If you do that, you must not call it "real numbers", but something like "real numbers where we forget the linear order".
I am not going to go into all the details, but let me just explain how the various parts of the structure give various operations on reals:

the Archimedean axiom is about computing rational approximations of reals
the field structure gives the usual arithmetical operations
the linear order gives us a semidecidable procedure for testing $x < y$
the Cauchy completeness gives us a function lim : (nat -> real) -> real which takes a (representation of) rapid Cauchy sequence and returns its limit. (A sequence $(x_n)_n$ is rapid if $|x_n - x_m| \leq 2^{-\min(n,m)}$ for all $m, n$.)

What we do not get is a test function for equality. There is nothing in the axioms for reals which asks that $=$ be decidable. (In contrast, the Peano axioms imply that the natural numbers are decidable, and you can prove that by implementing eq : nat -> nat -> Bool using only induction as a fun exercise).
It is a fact that the usual decimal representation of reals that humanity uses is bad because with it we cannot even implement addition. Floating point with infinite mantissa fails as well (exercise: why?). What works, however is signed digit representation, i.e., one in which we allow negative digits as well as positive ones. Or we could use sequences of rationals which satisfy the rapid Cauchy test, as stated above.
The Tsuyoshi representation also implements something, but not $\mathbb{R}$
Let us consider the following representation of reals: a real $x$ is represented by a pair $(q,b)$ where $(q_n)_n$ is a rapid Cauchy sequence converging to $x$ and $b$ is a Boolean indicating whether $x$ is an integer. For this to be a representation of the reals, we would have to implement addition, but as it turns out we cannot compute the Boolean flags. So this is not a representation of the reals. But it still does represent something, namely the subset of the reals $\mathbb{Z} \cup (\mathbb{R} \setminus \mathbb{Z})$. Indeed, according to the realizability interpretation a union is implemented with a flag indicating which part of the union we are in. By the way, $\mathbb{Z} \cup (\mathbb{R} \setminus \mathbb{Z})$ is a not equal to $\mathbb{R}$, unless you believe in excluded middle, which cannot be implemented and is therefore quite irrelevant for this discussion. We are of forced by computers to do things intuitionistically.
We cannot test whether a real is an integer
Finally, let me answer the question that was asked. We now know that an acceptable representation of the reals is one by rapid Cauchy sequences of rationals. (An important theorem states that any two representations of reals which are acceptable are actually computably isomorphic.)
Theorem: Testing whether a real is an integer is not decidable.
Proof. Suppose we could test whether a real is an integer (of course, the real is realized by a rapid Cauchy sequence). The idea, which will allow you to prove a much more general theorem if you want, is to construct a rapid Cauchy sequence $(x_n)_n$ of non-integers which converges to an integer. This is easy, just take $x_n = 2^{-n}$. Next, solve the Halting problem as follows. Given a Turing machine $T$, define a new sequence $(y_n)_n$ by
$$y_n = \begin{cases}
x_n & \text{if $T$ has not stopped within $n$ steps}\\
x_m & \text{if $T$ stopped in step $m$ and $m \leq n$}
\end{cases}$$
That is, the new sequence looks like the sequence $(x_n)_n$ as long as $T$ runs, but then it gets "stuck" at $x_m$ if $T$ halts in step $m$. Very importantly, the new sequence is also a rapid Cauchy sequence (and we can prove this without knowing whether $T$ halts). Therefore, we can compute its limit $z = \lim_n y_n$, because our representation of reals is correct. Test whether $z$ is an integer. If it is, then it must be $0$ and this only happens if $T$ runs forever. Otherwise, $z$ is not an integer, so $T$ must have stopped. QED.
Exercise: adapt the above proof to show that we cannot test for rational numbers. Then adapt it to show we cannot test for anything non-trivial (this is a bit harder).
Sometimes people get confused about all this testing business. They think we have proved that we can never test whether a real is an integer. But surely, 42 is a real and we can tell whether it is an integer. In fact, any particular real we come up with, $\sin 11$, $88 \ln 89$, $e^{\pi \sqrt{163}}$, etc., we can perfectly well tell whether they are integers. Precisely, we can tell because we have extra information: these reals are not given to us as sequences, but rather as symbolic expressions from which we can compute the Tsuyoshi bit. As soon as the only information we have about the real is a sequence of rational approximations converging to it (and I do not mean a symbolic expression describing the sequence, but a black box which outputs the $n$-th term on input $n$) then we will be just as helpless as machines.
The moral of the story
It makes no sense to talk about implementation of a set unless we know what sort of operations we want to perform on it.

REPLY [6 votes]: It is undecidable whether a given computable number is equal to zero.
(So your rational approximation oracle returns 0 for every ε you've tried?  Maybe you just haven't given it a small enough ε.)
Thus, it's undecidable whether a given computable number between -½ and +½ is an integer.

REPLY [2 votes]: A function being computable is a stronger than the function being continuous, i.e. any computable function needs to be continuous in the information topology.
You want to see if the function $F:\mathbb{R} \to \{Yes,No\}$ defined by
$$
F(r) =
\begin{cases}
YES & r\in \mathbb{Q}\\
NO & o.w.
\end{cases}
$$
is computable.
Let's assume that similar to what you wrote, real numbers are given by black-boxes and a real number black-box $r$ can be used to obtain a rational number of the form $\frac{k}{2^{n}}$ such that $r \in [\frac{k}{2^n},\frac{k+1}{2^n}]$ for any natural number $n$.
Then your function is not continuous and is therefore not computable.
Here is a more direct and intuitive adversary argument: Let $M$ be a machine. Assume that the number is $0$. Assume that for $n$ the black-box returns $[-\frac{1}{2^n},\frac{1}{2^n}]$. If $M$ doesn't halt we are done. If it halts, let the largest number that $M$ has asked an approximation for to be $m$. Then information that $M$ has seen up to this point is consistent with the real number being any real number in the interval $[-\frac{1}{2^m},\frac{1}{2^m}]$, i.e. $M$ doesn't have enough information to correctly solve the problem. If $M$ answers $NO$ then the answer is incorrect. If it answers $YES$ then we can consider running $M$ on any irrational number in the interval $[-\frac{1}{2^m},\frac{1}{2^m}]$ and $M$ will incorrectly answer $YES$ since it will get exactly the same information from the black-box. Therefore $M$ cannot solve the problem correctly.
The proof that any computable function needs to be continuous is similar.


TITLE: Inverse of Laplace operator
QUESTION [1 upvotes]: Let $\ L=\frac{-d^2}{dx^2}$ defined on $\ H^2(]0,1[) ∩ H^1_0(]0,1[)$ \ $\ L$ is the laplacian operator in one dimension \ how can we express the inverse of $\ L $

REPLY [1 votes]: The Green's function might need a slight modification.
\begin{align*}
G\left(x,y\right) & =\frac{2}{\pi^{2}}\sum_{n=1}^{\infty}\frac{1}{n^{2}}\sin\left(n\pi x\right)\sin\left(n\pi y\right)\\
 & =\begin{cases}
\left(1-y\right)x & 0\leq x\leq y\\
y\left(1-x\right) & y\leq x\leq1
\end{cases}.
\end{align*}
In details:
Let $L=-\frac{d^{2}}{dx^{2}}$, with domain as specified. Then $L$
is selfadjoint, and it has eigenfunctions (normalized) $\left\{ \varphi_{n}\left(x\right):=\sqrt{2}\sin\left(n\pi x\right)\mathrel{;}n\in\mathbb{N}\right\} $,
so that
\begin{align*}
L^{-1} & =\sum_{n=1}^{\infty}\lambda_{n}^{-1}\varphi_{n}\left(x\right)\varphi_{n}\left(y\right)=\frac{2}{\pi^{2}}\sum_{n=1}^{\infty}\frac{1}{n^{2}}\sin\left(n\pi x\right)\sin\left(n\pi y\right).
\end{align*}
For the Green's function, assume $LG\left(\cdot,y\right)=\delta\left(\cdot-y\right)$,
$y\in\left(0,1\right)$, subject to the condition $G\left(0,y\right)=G\left(0,1\right)=0$.
So
\begin{align*}
G\left(x,y\right) & =\begin{cases}
ax & x\leq y\\
c\left(x-1\right) & x\geq y
\end{cases}
\end{align*}
where $a=\frac{c\left(y-1\right)}{y}$, since $G\left(\cdot,y\right)$
is continuous at $x=y$. But
\begin{align*}
\lim_{x\rightarrow y^{-}}G_{x}\left(x,y\right)-\lim_{x\rightarrow y^{+}}G_{x}\left(x,y\right) & =1,
\end{align*}
so $c=-y$.


text:
\begin{definition}[Definition:Circumcircle of Triangle]
The '''circumcircle of a [[Definition:Triangle (Geometry)|triangle]]''' is the [[Definition:Circle|circle]] that passes through all three of the [[Definition:Triangle (Geometry)|triangle]]'s [[Definition:Vertex of Polygon|vertices]].
This [[Definition:Circle|circle]] is said to [[Definition:Circle Circumscribed around Polygon|circumscribe]] the [[Definition:Triangle (Geometry)|triangle]].
:[[File:Circumcircle.png|320px]]
=== [[Definition:Circumcircle of Triangle/Circumcenter|Circumcenter]] ===
{{:Definition:Circumcircle of Triangle/Circumcenter}}
=== [[Definition:Circumcircle of Triangle/Circumradius|Circumradius]] ===
{{:Definition:Circumcircle of Triangle/Circumradius}}
\end{definition}

config: wiki


TITLE: How to find distance s1?
QUESTION [0 upvotes]: Here is a screenshot of a problem from my textbook

My book says
$$S₁= 2t (\tan(\theta_{2})\cos(\theta_{1}))$$
I know the distance for $\overline{AC}$ is
$$2t (\tan(\theta_{2})) $$
but where did $\cos(\theta_1) $  come from?
How is  $S_{1}=$ distance between $\overline{AC}$ and $(\cos(\theta_{1}))$?

REPLY [1 votes]: Please see the diagram.
From $\triangle ADB$, we can calculate that $DB = t\tan\theta_2$. Now $AC = 2DB = 2t\tan\theta_2$.
Denote point $X$ such that $\triangle AXC$ is a right triangle, with the right angle being at $X$. $\angle ACX = \theta_1$, by the theorem of corresponding angles. Since $AC$ is the hypotenuse of this triangle, and $CX$ is the adjacent side to the angle $ACX$, $s_1 = AC\cos\theta_1 = 2t\tan\theta_2\cos\theta_1$, as required.
Just to add on, the comment on the question bringing up unstated assumptions is a very fair one. But this problem seems to be one on the classical theory of refraction and partial/total internal reflection in the topic of optics (in physics), so the ray diagram and its inherent assumptions are fairly obvious. Nevertheless, it would be good for the asker to confirm this.


TITLE: Forces Problem from Halliday, Resnick and Krane
QUESTION [0 upvotes]: To my knowledge, the net force on $m_1$ is $T - m_1 \cdot g \cdot \sin(\theta_1) = m_1 \cdot a$ assuming that the block accelerates up the incline. However, the solutions manual states that the acceleration is given by $m_1 \cdot a = g\cdot m_1 \cdot \sin(\theta_1) - T$ in that same situation.
I have been at this for half an hour now and can't understand how acceleration is opposite in direction to the tension force when moving up the incline. Thanks in advance for the help.

REPLY [0 votes]: It looks like you're choosing the positive coordinate axis to be up the slope, while the solution is choosing the positive axis down the slope. The motion will be the same regardless of how one chooses coordinates (as it must be).


TITLE: Calculatin differential
QUESTION [0 upvotes]: I have this seemingly simple question that gives me hard time:
Subset $\Omega \subset \mathbb{R}^n$ is open and $f\in C^1(\Omega;\mathbb{R}^n)$. If $x,y\in\Omega$, where $x=(x_1,...,x_n)$, then we have following $$f(x)-f(y)=\int_0 ^1 \frac{\partial}{\partial t}f(y+t(x-y))dt= \int_0 ^1 \sum_{k=1}^{n}\frac{\partial}{\partial x_k}f(y+t(x-y))(x_k-y_k)dt$$ $$=\int_0 ^1 f'(y+t(x-y))(x-y)dt$$
My calculus is a bit rusty, so can you help explaining some of the previous steps?

REPLY [1 votes]: Use the chain rule, let $z = y + t (x-y)$:
$$ \frac{\partial f(y + t (x-y))}{\partial t} = f'(z)\frac{\partial z}{\partial t} = f'(z)(x-y)$$
Since you are working in $\mathbb{R}^n$ then actually
$f'(z) (x-y)$ means $(x-y)\cdot \nabla f(z)$


TITLE: Confusion in notation, problem from Herstein's book.
QUESTION [0 upvotes]: Let $G$ be a nonempty set closed under an associative product, which in
  addition satisfies :
A. There exists an $e$ in G such that $a⋅e=a$ for all $a∈G$.
B. Given $a∈G$, there exists an element $y(a)∈G$ such that $a⋅y(a)=e$.
Prove that $G$ must be a group under this product.

This is a question from I N Herstein. It was already asked before. Here is the solution
https://math.stackexchange.com/a/1200617/581242

Let, $ab=e\land bc=e\tag {1}$ for some $b,c\in G$. And, $ae=a\tag{2}$
  From $(2)$, $$eae=ea\implies(ab)a(bc)=ea\implies
 ((ab)(ab))c=ea\implies ec=ea\tag{3}$$
Similarly, $$ae=a\implies a(bc)=a\implies (ab)c=a\implies
 ec=a\tag{4}$$
$(3)$ and $(4)$ implies, $$ea=a$$
Also from $(3)$ and $(1)$, $$(bab)(bca)=e\implies
 b((ab)(bc)a)=e\implies ba=e$$

I don't understand the notation used in first line. What does $$ab=e\land bc=e$$ mean?

REPLY [1 votes]: In logic, $\land$ means "and" and $\lor$ means "or."


TITLE: I want to know the volume of the enclosed space 0<xy<z<1
QUESTION [0 upvotes]: Here, $x,y$ and $z$ are real numbers and they are satisfied with these all conditions.
$0 \leqq x \leqq 1$
$0 \leqq y \leqq 1$
$0 \leqq z \leqq 1$
$0 \leqq xy \leqq z \leqq 1$
And I would like to know the volume of the space which $x,y$ and $z$ satisfied with these conditions in $\mathbb{R}^3$
In this case, should I calculate below integral?
$$ \int_0^1 dz \int_0^{z/x} dy \int_0^{z/y} dz$$

REPLY [1 votes]: Your integration limits don't make sense, even after correcting the last $\mathrm{d}z$ to $\mathrm{d}x$.  Just read it from left to right.  First we have $0\leq z\leq1$.  That's fine.  Then we have $0\leq y\leq {z\over x}$, but $x$ hasn't been defined yet, so this can't be right.
The way to do it is shown in the comments.  $x$ and $y$ can take any values in $[0,1]$ and $z$ must be $\geq$ their product.


TITLE: Understanding why Lagrange Interpolation works
QUESTION [2 upvotes]: I'm confused as to why dividing $p(x) = q(x) / q(x_1)$ works. I understand that $q(x)$ is designed such that any $x_i$ for $2 \leq i \leq d+1$ will make $q(x)=0$. I know that $q(x_1)$ will not make $q(x) =0$ because there is no $(x-x_1)$. So, you keep generating a $q(x)$ that is non-zero for $x_i$ for all $2 \leq i \leq d+1$, and then divide by their respective $q(x_i)$ ? How does this give us $p(x)$ ?

REPLY [1 votes]: The method described in your question explains how to obtain a solution when $y_1=1$ and all other $y_i$ are zero, but of course it doesn't give you the general solution. But once you have this, there isn't much left to do.
Denote $p_1(x)$ the polynomial obtained by this method. Using the same construction, you can obtain other polynomials $p_i(x)$, which are such that $p_i(x_i)=1$ and $p_i(x_j)=0$ for all $i,j$.
Once you have all these polynomials, you can find the interpolating polynomial for any given $y_1,\dots y_{d+1}$; it is simply equal to
$$p(x)=\sum_{k=1}^{d+1}y_kp_k(x).$$
Indeed, you can check that this polynomial is such that
$$p(x_i)=\sum_{k=1}^{d+1}y_kp_k(x_i)=y_ip(x_i)=y_i$$
(the other terms in the sum being all zero).


TITLE: Proving $f: \mathbb{R}^n \to \mathbb{R}^m$ is continuous if $f^{-1} (O)$ := {${x \in \mathbb{R}^n} : f(x) \in O$} $\subset \mathbb{R}^n$
QUESTION [0 upvotes]: How can I prove that a transfomation $f: \mathbb{R}^n \to \mathbb{R}^m$ is continuous if and only if the inverse image $f^{-1} (O)$ := {${x  \in \mathbb{R}^n} : f(x) \in O$} $\subset \mathbb{R}^m$ of every open set $O \subset \mathbb{R}^n$ is also an open set?
I know that there exists a continuous function, for example $h:\mathbb{R}^n \to \mathbb{R}^m$ and an open set $O \subset \mathbb{R}^n$, so that $h(O) \subset \mathbb{R}^m$ isn't open. I could just take $h \equiv 0$ and $O = \mathbb{R}^n$, then $h(O)$ = {$0$} $\subset \mathbb{R}^m$ which isn't open. This implies that the statement doesn't work for images of open sets.
Further I know that a transformation $f: \mathbb{R}^n \to \mathbb{R}^m$ is  continuous exactly when the inverse image $f^{-1}(A)$ of every closed set $A \subset \mathbb{R}^m$ is a closed set.
But that still doesn't help me much with how proving the statement on top..

REPLY [1 votes]: First of all, usually the definition of a continuous map $f:X \to Y$ is that for every open $U \subset Y$ $f^{-1}(U)\subset X$ is open.
So if you already know the statement for closed sets: Let $f: X \to Y$ be "continuous" (your closed set definition). Consider $U \subset Y$ open, then $Y\setminus U$ is closed. and hence $f^{-1}(Y\setminus U)= X \setminus f^{-1} (U) $ is closed, wich means $f^{-1}(U)$ is open.


TITLE: Reducing trigonometry functions in exponents
QUESTION [0 upvotes]: I need to know a bit about replacing the trig functions in exponents with simpler algebraic equations which can then be reduced with logarithms. For example, an answer to this question said that
$$e^{\sin(x)} = e^{\frac{1}{2i}(e^{ix} - e^{-ix})}$$
applies, which would do the trick. As well as this, using a few basic log rules it could be adjusted to deal with any base as well. I do have a few questions though.
Firstly, how does this work? It looks like it has a fair bit to do with Euler's identity but I don't understand the exact derivation.
Next, how do I generalize it (For example, with the exponent as $4*sin^2(3x)$)? I honestly don't know if it barely involves changing numbers or if it doesn't apply at all anymore.
Finally, is this transferrable to the other trig functions? If it's connected to Euler's identity as I think it is, it should be easy to change it to apply to cosine instead, and then the two together may somehow be used to derive the transformation (if I can call it that) to tan.
Thank you so much for your help!

REPLY [1 votes]: The derivation of the definitions of $\sin x$ and come from subtracting Euler's formulas.
\begin{align}
e^{ix}&=\cos x+i\sin x\\
e^{-ix}&=\cos x-i\sin x\\
\implies\sin x&=\frac1{2i}(e^{ix}-e^{-ix})
\end{align}
It therefore follows that $$\exp(\sin x)=\exp\left(\frac1{2i}(e^{ix}-e^{-ix})\right).$$
For your example, we can write
\begin{align}
4\sin^2(3x)=4(\sin(3x))^2&=4\left(\frac1{2i}(e^{3ix}-e^{-3ix})\right)^2\\
&=-(e^{3ix}-e^{-3ix})^2\\
&=2-e^{6ix}-e^{-6ix}\\
\implies\exp(4\sin^2(3x))&=\exp(2-e^{6ix}-e^{-6ix}).\\
\end{align}