\section{Introduction}
\label{sec:introduction}
Cooperative unmanned aerial vehicle (UAV) networks, also known as flying ad-hoc networks (FANET), recently become extremely popular for their ease of set-up and use, low price, and high maneuverability in 3D space \cite{WalidSurvey}. In such networks, nodes work cooperatively to make the communication of two far away nodes possible by acting as intermediate repeaters. Using such networks where implementing the infrastructure is not feasible, or the infrastructure is damaged, is crucially important. Wildfire monitoring \cite{wildFire1}, search and rescue after a disaster \cite{disaster}, ISR \cite{battlefield}, ad hoc UAV cloud service provider \cite{gtss}, and border surveillance \cite{border} are some examples to name. A promising architecture for such networks uses software-defined networking (SDN) as its cornerstone to add flexibility to the network control and  management \cite{178}. SDN separates the control plane from the data plane to pass the decision-making process to a software-based controller. There are excellent efforts in designing SDN-based topology management \cite{178}, resource management \cite{180}, path planning \cite{pathPlanning}, node positioning \cite{248}, and monitoring \cite{214}. The SDN-based routing problem, of how to find the optimal path from the source node to the destination, is also tackled for multi-path routing \cite{multipath}, and low mobility UAV networks \cite{lowMob}. However, the SDN-based routing problem for highly dynamic UAV networks is still a challenging unresolved problem to the best of our knowledge.

In this paper, we consider a cooperative UAV network with a central controller. We assume a low-bandwidth reliable communication channel exists between the controller and the UAV nodes. However, the UAVs need to transfer a high volume of data, such as high-resolution pictures, videos, or thermal pictures. Hence, they need a high bandwidth communication channel for data transfer. In this way, the control plane and the data plane are separated, just like any other SDN-based UAV network. It is worth mentioning that the controller can be a ground or aerial controller with a line of sight or non-line of sight communication with the nodes. Such communication could be borrowed from the existing cellular network, be provided by the line-of-sight communication with a controller placed at a higher elevation, or use a dedicated bandwidth for command and control, as it is the common assumption in many UAV networks \cite{WalidSurvey}. The controller is responsible for decision-making in all network management processes. We assume the routing, as one of the control management tasks, is also one of the controller responsibilities. Fig. (\ref{fig::systemModel}) represents a schematic view of our system model. Since the controller is responsible for all network management tasks, any algorithm required to be performed by the controller is desired to be with a low computation and communication complexity to prevent any potential delay in network management. This fact makes the routing problem more challenging in such a network architecture.

\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{figs/architecture.pdf}
    \caption{System model (solid lines represent a direct link where the dashed lines represent a communication which may or may not be direct)}
    \label{fig::systemModel}
\end{figure}

Cooperative UAV networks typically use the conventional routing algorithms designed for the conventional ad-hoc networks such as Mobile Ad-hoc NETwork (MANET) and Vehicular Ad-hoc NETwork (VANET). Ad-hoc On-demand Distance Vector (AODV) \cite{aodv}, Destination Sequenced Distance Vector (DSDV) \cite{dsdv}, and Optimized Link State Routing (OLSR) \cite{olsr} are some instances. Such algorithms aim at finding the shortest path between the source and the destination and do not consider the other network characteristics. Furthermore, the highly dynamic nature of UAV networks leads to fast disconnection and consequently a need for multiple reroute processes \cite{tcpreorder}. Multiple reroute processes dramatically increase the routing communication complexity as well as overall end-to-end communication delay. In this case, the data packets belonging to the same data file/stream may traverse different paths, and hence they arrive out-of-order. The same sideeffect is tackling the multi-path routing algorithms. TCP deals with some of out-of-order packets as the lost packets and orders retransmission.

We recently published an optimized predictive and adaptive routing (OPAR) \cite{opar} for cooperative UAV networks to consider the path lifetime in addition to the path length in the path selection process. The aim was to reduce the reroute process's overheads. We find that considering the path lifetime can significantly improve UAV networks' performance in terms of throughput, flow completion time (FCT), and flow success rate. We further find that considering only the path lifetime and path length in selecting the paths may cause an unbalanced network load. Since the network uses a shared medium, unbalanced load in the network may cause performance degradation in throughput and consequently FCT, especially in the highly loaded networks. Accordingly, in this paper, we extend the OPAR algorithm to propose a load-balanced OPAR (LB-OPAR) routing algorithm by considering the network load in addition to the path lifetime and path length in selecting the paths. We show that this consideration does not impose any extra routing traffic to the network in comparison with OPAR, in the SDN-based setting. Accordingly, we model the entire problem of the SDN-based routing algorithm in a cooperative UAV network with an analytical optimization model. The proposed model is binary linear programming (BLP) problem. While BLP problems are very well-known to be  NP-complete problems, we show that in the specific case of the modeled BLP, we can find the optimized path using a graph-based algorithmic solution. The algorithmic solution is lightweight with the computational complexity of $O(|E^2|)$ where $|E|$ is the number of network links.

LB-OPAR needs only small-size messages to be sent to the controller periodically. Each message contains three consequent positions of the UAV node. The controller predicts the link lifetime of the UAV nodes and uses the predicted values in solving the optimization problem. The prediction algorithm is lightweight, with a computation complexity of $O(|V|^2)$ where $|V|$ is the number of UAV nodes. The controller uses its own data to calculate the communication load on each node. Hence, no further communication is required. The low communication complexity and the low computational complexity of the LB-OPAR algorithm make it a perfect fit for the SDN-based model.

To evaluate the performance of LB-OPAR, we use the ns-3 network simulator \cite{ns3}. We simulate the network to compare LB-OPAR with the well-known baseline routing algorithms AODV, DSDV, and OLSR as benchmarks. We compare the results with those of OPAR as well. We show in Section (\ref{sec:RelatedWork}) that the state-of-the-art uses these benchmarks to find the routes in the cooperative UAV networks or build new routing algorithms on the basis of these benchmarks. We selected the algorithms for comparison that cover the wide range of reactive and proactive algorithms as well as distance vector and link-state ones. We measure the network throughput, FCT, and flow success rate as the performance metrics. We show that LB-OPAR improves the network throughput by up to four times the benchmarks. It decreases the FCT for average $20\%$ and increases the flow success rate by about $30\%$.

This work's main contribution is to propose a load-balanced routing algorithm for SDN-based cooperative ad-hoc networks. More specifically, $i)$ it analytically models the routing problem in the mentioned systems with linear optimization technique and considers the path length, lifetime, and load altogether; $ii)$ it proposes an algorithmic solution for the proposed problem; $iii)$ it analytically proves the solution's complexity; $iv)$ it proposes a lightweight highly accurate link lifetime prediction method; $v)$ it exhaustively evaluates the performance of the proposed algorithm and compares it with the benchmark solutions.

The rest of this paper is organized as follows. The state-of-the-art is reviewed in Section (\ref{sec:RelatedWork}). We propose the analytical model with its theoretical solution of the LB-OPAR in Section (\ref{sec::proposedAlg}). The performance of LB-OPAR is evaluated and compared to the baseline solutions in Section (\ref{sec:NumericalResults}). Finally, we conclude the paper in Section (\ref{sec::conclusion}).

\section{Related Work}
\label{sec:RelatedWork}
Cooperative UAV networks typically use conventional routing algorithms, which are proposed originally to be used in mobile ad-hoc networks (MANETs), vehicular ad-hoc networks (VANETs), or wireless sensor networks (WSNs). However, there are several routing algorithms dedicatedly proposed for cooperative UAV networks. Thus, in this section, we first review the conventional routing algorithms, then we present the algorithms proposed dedicatedly for UAV networks. Finally, we review the proposals for SDN-based UAV network routing algorithms to cover the state-of-the-art.

Generally, there are two categorizations for conventional ad-hoc routing algorithms \cite{adhocRoutin}. Each of them looks at the problem from a different viewpoint. However, they all are proposed to find the path with the least number of hop counts, and they do not consider other metrics such as path quality, reliability, or lifetime. The first looks at the algorithms based on their nature of being proactive, reactive, or hybrid. Proactive routing protocols periodically keep track of the network topology changes, regardless of the demand for new routes. Optimized link-state routing protocol (OLSR) \cite{olsr} and destination sequence distance vector (DSDV) \cite{dsdv} are two well-known instances of this category. In contrast to proactive protocols, reactive protocols attempt to find a route only when a new route is needed. Ad-hoc on-demand distance vector (AODV) \cite{aodv} is a very well-known instance of this category.
The second categorization considers the method of implementing the routing algorithms to be link state or distance vector. In the link-state approach, each node builds a general view of the entire network topology to find the optimized paths. The OLSR protocol belongs to this category. In contrast, each node needs to keep only the next hop toward the destination and the path cost in the distance-vector approach. Both AODV and DSDV use this approach in their routing.

Most of the routing algorithms proposed dedicatedly for UAV networks are built based on the conventional ad-hoc routing algorithms and extend them to consider another metric in their path selection process. As an instance, Bomio et al. \cite{ON} proposed enhanced AODV (E-AODV) by taking into account a route reliability criterion. They considered a 2D area to check whether two nodes are getting closer or not. If their distance is going to be increased, the corresponding link is discarded. Li et al. \cite{LEPR} proposed link-stability estimation-based preemptive routing (LEPR) based on AODV. LEPR finds multiple disjoint paths using the AODV algorithm, calculates the expected transmission count (ETX) metric of \cite{ETX} for each link of the route, and keeps all routes to use them when the path with the best ETX becomes not available anymore.

Pu et al. \cite{ltaolsr} proposed link quality and traffic-load aware OLSR (LTA-OLSR) to find the path with the highest quality from the set of the shortest paths found by OLSR. Song et al. \cite{olsrpdm} proposed OLSR-based mobility and delay prediction (OLSR-PMD) to choose the more stable next-hop nodes based on a prediction algorithm. Rosati et al. \cite{polsr} proposed predictive OLSR (P-OLSR) to add weight to the ETX metric of each link based on the corresponding nodes' relative speeds. Lyu et al. \cite{qngpsr} proposed  Q-Network Enhanced Geographic Ad-Hoc Routing Protocol based on GPSR (QN-GPSR). Authors reduced the premier forwarding of GPSR \cite{gpsr} algorithm by applying Q-learning, a reinforcement learning algorithm, into the GPSR routing protocol. Their results show that QN-GPSR has a higher packet delivery ratio and a lower end-to-end delay than the basic GPSR protocol. However, applying a Q-network increases the computational complexity of the GPSR algorithm by order of magnitude. Liu et al. \cite{qmr} proposed another Q-learning based Routing protocol for Flying Ad Hoc Networks and named it Q-learning based multi-objective routing (QMR). Authors show that QMR provides a higher packet delivery ratio, lower delay, and energy consumption than existing Q-learning-based routing protocols. However, the vast imposed computational complexity by the Q-learning algorithm is the main drawback of this algorithm. Adaptive and Reliable Routing Protocol with Deep Reinforcement Learning (ARdeep) \cite{ardeep} is also proposed to be applied to UAV networks. Although the authors show that it performs better than the basic GPSR algorithm in terms of average end-to-end delay and packet delivery ratio, it increases the computational complexity by order of magnitude, just like other deep learning-based algorithms. Although most of the computational complexity source of ML algorithms belongs to the training phase, a learned model might not be able to lead to high precision results in a network with a highly dynamic nature \cite{ROVIRASUGRANES2022102790}. To address this shortcoming, the learning-based methods keep the learning phase active during the network operation time to be regularly adjusted with the new live data where such online learning needs high computational complexity.

Several other routing algorithms are proposed dedicatedly for UAV networks, but they are special-purpose algorithms or only applicable to the networks with specific missions. For instance, Sunder et al. \cite{fueling} designed a routing protocol for a UAV-based fueling system. Coelho et al. \cite{redefine} designed a routing algorithm for UAV networks where UAVs are used as WiFi repeaters to maximize the coverage. The work of \cite{packageDelivery} designed for package delivery networks. Cheriguene et al. \cite{sempr} proposed a Swarm energy-efficient multicast routing protocol for UAVs flying in group formations, to facilitate the control and information delivery among the UAVs while minimizing inter-UAV packet loss, packet re-transmission, and end-to-end delay.
It is obvious that such mission-specific routing algorithms cannot be used as general-purpose routing for all cooperative UAV networks.

An SDN architecture is proposed to be used in UAV networks to increase the network control flexibility and improves the network performance. Different architectures are proposed for SDN-based UAV networks \cite{191,qos,infocomwkshp2018,SDNtraffic}. All the proposed architectures have a controller or a set of distributed controllers that has/have direct, reliable communication with the UAV nodes. Using such an architecture, any of the network management tasks could be migrated into the controller. However, the migrated tasks do not impose huge computational or communication complexity into the controller. Routing, as an essential network management task, could be migrated into the controller. To the best of our knowledge, the problem of finding an optimal path in an SDN-based UAV network has not yet been substantially studied. Several recently published works mentioned this problem \cite{lowMob,multipath,SDNtraffic}.

Ramaprasath et al. \cite{lowMob} proposed an SDN-based routing algorithm that utilizes AODV in the initialization phase of the network. Alas, this proposal cannot capture the high dynamicity of the UAV networks due to its multiple long wait times during the routing procedure. Hence, this algorithm's networks of interest are the UAV networks with deficient mobility. Authors of \cite{multipath} proposed an SDN-based routing algorithm that finds all possible paths toward the destination and divides the entire data file among them. Each part of the data file is then sent through one of the paths. As the out-of-order packet arrival is one of the main challenges of the highly dynamic networks, this works increases the number of such packets by using multiple paths. Authors also assume that each UAV is equipped with different network interfaces. Hence, it can send the data file part from different paths concurrently, which poses a hardware limitation for networks aimed at using this algorithm. Qi et al. \cite{SDNtraffic} proposed traffic differentiated routing (TDR), a new SDN architecture for UAV networks that includes two types of controllers, cluster controllers and a coordination controller. In this architecture, the network is clustered, and each cluster needs a controller. To arrange different cluster controllers, a coordination controller is also required. The routing algorithm of this paper is proposed to find the routes inside the clusters. Besides their new network architecture, authors assumed that the controller is aware of the current location and the nodes' speed. Using such information, the authors proposed an optimization problem in which finding the optimal solution is not practically feasible. Hence, they use approximation methods to find a near-optimal path. As we mentioned earlier, the routing problem in SDN-based cooperative UAV networks needs substantial work. Accordingly, we propose LB-OPAR as a general-purpose routing algorithm for SDN-based UAV networks in this paper. Table (\ref{tbl::relatedWork}) represents a comparison between routing algorithms proposed dedicated for UAV networks or could be used in them, including LB-OPAR.

\begin{table}[t]

\caption[]{Routing Algorithms for Cooperative UAV Networks. (The comparison of this table is partially based on the results of  Section (\ref{sec:NumericalResults}) of this paper.)}
\resizebox{.7\columnwidth}{!}{
\begin{minipage}{\textwidth}
\begin{tabular}{ | l | l | l | l | l | l | l | l | }
\hline
Algorithm  &Mobility &Routing  &Commun.&Comput. &SDN-based/ &Load  &Success\\
&Support &  Metric &  Overhead&  Overhead&Fully-dist. & Balancing&Rate\\\hline
AODV\cite{aodv} &Good &Path len.&High &Low & Fully-dist.&$\qquad\times$&Medium \\
DSDV\cite{dsdv} & VLimited &Path len.  &VHigh &Low &Fully-dist. &$\qquad\times$&VLow\\
OLSR\cite{olsr}& VLimited &Path len.&VHigh &Low &Fully-dist. &$\qquad\times$&VLow\\
E-AODV\cite{ON}&Good &Distance&High & Low& Fully-dist.&$\qquad\times$&Low\\
LEPR\cite{LEPR}&Good &ETX&High & Low &Fully-dist. &$\qquad\times$&Medium \\
LTA-OLSR\cite{ltaolsr}&VLimited &RSSI+load  &VHigh &Low & Fully-dist.&$\qquad\checkmark$&VLow\\
OLSR-PMD\cite{olsrpdm}&VLimited &Delay  &VHigh &Low &Fully-dist. &$\qquad\times$&VLow\\
P-OLSR\cite{polsr}&VLimited &Weighted ETX  &VHigh &Low & Fully-dist. &$\qquad\times$&VLow\\
QN-GPSR\cite{qngpsr}&Limited &Link Quality  &VHigh & VHigh& SDN-based&$\qquad\times$&Low\\
ARdeep\cite{ardeep}&Medium & PER\footnote{Packet error rate} &-&VHigh & SDN-based&$\qquad\times$ &-\\
QMR\cite{qmr}&Limited & Delay \& Energy &VHigh & VHigh& SDN-based&$\qquad\times$&Low\\
\cite{multipath}&VLimited &Multipath&VHigh &Low & SDN-based&$\qquad\times$&Low\\
\cite{lowMob}&VLimited &Path Len.  &High &VHigh & SDN-based&$\qquad\times$&VLow\\
TDR\cite{SDNtraffic}&VLimited &Link Availability&- &VHigh & SDN-based&$\qquad\times$&-\\
OPAR\cite{opar} &VGood &Path len.  &Low & VLow&SDN-based &$\qquad\times$&VHigh\\
& & \& lifetime & & &&&\\
LB-OPAR &VGood & Path len.,  &VLow &Low &SDN-based &$\qquad\checkmark$&VHigh\\
& & load \& lifetime & & &&&\\\hline
\end{tabular}
\label{tbl::relatedWork}
\end{minipage}}
\end{table}

\section{Proposed Algorithm}
\label{sec::proposedAlg}

We model the routing problem in SDN-based cooperative UAV networks with a linear programming (LP) optimization model. The optimization model considers path length, path lifetime, and network load altogether. We define the path lifetime as the minimum lifetime of all the links forming the path. Hence, we need to know the link lifetime of each network link. Obviously, the link lifetime depends mostly on the nodes' movement.  Hence, we propose a link lifetime prediction algorithm to support LB-OPAR. We further consider the path load as the maximum node load for the nodes participating in the path. We define the node load as the number of the neighboring nodes engaged with data transmission, affecting the medium access accessible with the corresponding node. In other words, the node load is the number of nodes, in the corresponding node's adjacency, that are sending or forwarding data packets, i.e. they are a source of a data flow or engaged as forwarding nodes in a communication flow. Since the controller manages the routes and knows the nodes' location, it can calculate nodes' load. The proposed optimization problem variables are Binary, which means that the model is a Binary LP problem. The Binary LP problems are well-known to be non-polynomial problems in the NP-complete category \cite{karp}. However, we show that the proposed problem is an especial case in which its solution could be found using the graph-based algorithmic method.

In this section, we propose the LP problem, considering that we know the links' lifetime and nodes' load. We then propose the graph-based algorithmic solution, which can find the solution of the Binary LP problem. Finally, we present the details of how to predict the link lifetime and calculate the node load by the controller. Table (\ref{tbl::notation}) represents the notations used throughout this paper.

\begin{table}[t!]
\caption[]{Table of notations }
\resizebox{1\textwidth}{!}{
\begin{minipage}{\textwidth}
\begin{tabular}{ l | l  }
  $n$ & Number of UAVs\\
  $u_i$ & The $i^{th}$ UAV\\
  $G(V,E)$ & Network Graph\\
  $V$ & Set of graph vertices\\
  $E$& Set of graph edges\\
  $e(i,j)$ & The edge from node $u_i$ to node $u_j$\\
  $\mathcal{T}$& Link lifetime prediction matrix\\
  $\mathcal{L}$& The vector of nods' load.\\
  $\tau_{(i,j)}$& Link lifetime of $e(i,j)$\\
  $l_i$& The load of node $u_i$\\
  $x_{(i,j)}$& Binary variable corresponding to $e(i,j)$\\
  $T$& The inverse of path lifetime\\
  $w_1$ & Path length weight\\
  $w_2$ & Path lifetime weight\\
  $w_3$ & Path load weight\\
  $\varsigma$& The matrix of all lifetime-load costs\\
  $\varsigma_{(i,j)}$& the lifetime-load cost of $e_{(i,j)}$\\
  $\varsigma_{path}$&The maximum lifetime-load cost of a path\\
  $x_i(t_0)$&The position of $u_i$ at time $t_0$ along x-axis in 3D space\\
  $y_i(t_0)$&The position of $u_i$ at time $t_0$ along y-axis in 3D space\\
  $z_i(t_0)$&The position of $u_i$ at time $t_0$ along z-axis in 3D space\\
  $\alpha_i$& The azimuthal degree of node $u_i$ direction\\
  $\theta_i$& The polar degree of node $u_i$ direction\\
  $v_i(t_0)$& The speed of node $i$ in time $t_0$\\
  $a_i(t_0)$& The acceleration of node $i$ in time $t_0$\\
  $d_{(t_0)}^{(i,j)}$& Euclidean distance between $u_i$ and $u_j$ at time $t_0$\\
    $d_{(t_0,t_1)}^{(i)}$& Euclidean distance traversed by node $u_i$ in time $[t_0\quad t_1]$\\
  $\Re$& Communication range
\end{tabular}
\label{tbl::notation}
\end{minipage}}
\end{table}

\subsection{LB-OPAR: SDN-Based load balanced Optimized Routing}
\label{sec::optimization}

On the one hand, it is generally a challenging problem to consider multiple metrics in defining an optimization problem. Multi-objective linear programming (MOLP) is a trivial modeling technique in such cases. MOLP solution, if feasible, leads in a set of all efficient extreme points, i.e. all maximal efficient faces. The optimal solution, then, is based on the weights of importance of different objective functions. On the other hand, finding a path requires an optimization problem with integer or Binary variables, where solving the optimization problems, in most cases, results in non-integer solutions. Rounding the non-integer solutions to the closest integer value leads to a non-optimal solution. Finding the optimized solution of most of the proposed optimization problems in routing and path planning is not feasible \cite{SDNtraffic,fueling,packageDelivery}. Accordingly, the optimization problem's design has to be such that finding the optimal solution becomes feasible, if possible. Hence, we carefully designed the following optimization problem, which is a Binary LP problem. In this problem, we assume that we have a matrix of all links' lifetimes, i.e. $\mathcal{T}$, as well as the vector of all nodes' load, i.e. $\mathcal{L}$.

\begin{eqnarray}
 \min_{} & & \sum\limits_{\substack{(i,j)\in E\\i\in\{1,\cdots,n\}\\j\in\{1,\cdots,n\}}}  w_1 x_{(i,j) } + w_2 T +  w_3  L
 \label{TimeObj} \nonumber\\
   \mbox{Subject To:} & &  \sum\limits_{(s,i)\in E} x_{(s,i)} = 1 \label{Const1} \\
 & & \sum\limits_{(i,s)\in E} x_{(i,s)} = 0 \label{Const2} \\
& &\sum\limits_{(i,d)\in E} x_{(i,d)} = 1 \label{Const3} \\
& & \sum\limits_{(d,i)\in E} x_{(d,i)} = 0 \label{Const4} \\
& & \sum\limits_{\substack{(i,j)\in E\\ i\neq s }} x_{(i,j)} = \sum\limits_{\substack{(j,k)\in E \\ j\neq d}} x_{(j,k)} \label{Const5} \\
 && T \geq \frac{x_{(i,j)}}{\tau_{(i,j)}} \label{Const6}\\
 & & x_i =\sum\limits_{\substack{j=0\\ j\ne i}}^{n} x_{(i,j)} \label{Const7}\\
 && L\geq x_i l_i \label{Const8}\\
 & & x_{(i,j)} \in \{0,1\} \label{Const9}\\
  & & 0 \leq T \leq1 \label{Const10}
\end{eqnarray}

The objective function of this problem is to minimize the path length as well as path load, and maximize the path lifetime. The significance of each of the path length, lifetime, and load is set by setting the weights $w_1$, $w_2$, and $w_3$, respectively, where $\sum_{i=1}^{3}w_i=1$. The Binary variables $x_{(i,j)}$s represent whether the corresponding network link $e_{(i,j)}$ participate in the optimal path or not. If the link $e_{(i,j)}$ is selected to be a part of an optimized path, $x_{(i,j)}$ is set to one, elsewhere it is set to zero. In this objective function, $T$ and $L$ represent the inverse path lifetime, and the path load, respectively. While $T$ has to be minimized according to the objective function, it has to be greater or equal to $\frac{x_{(i,j)}}{\tau_{(i,j)}}$ for all links $x_{(i,j)}$s, according to Constraint (\ref{Const6}). It is the idea to find the path lifetime, i.e. the minimum link lifetime for the links participated in the path, and consider it in the objective function. Similarly, the path load $L$ has to be minimized. To find the path load we minimize $L$ in the objective function where in Constraint (\ref{Const8}) we force $L$ to be greater than or equal to $x_il_i$, where according to Constraint (\ref{Const7}) $x_i =\sum_{\substack{j=0, j\ne i}}^{n} x_{(i,j)}$, and $l_i$ is the network load of node $i$.

We need Constraints (\ref{Const1}) and (\ref{Const2}) to let the source node start the path and prevent it from being in a loop. Similarly, we need Constraints (\ref{Const3}) and (\ref{Const4}) to let the destination node be chosen as the last node in the path and terminate the path selecting operation. Constraint (\ref{Const5}) assures that the intermediate node, which receives an active link, has to have an outgoing active link. As we discussed earlier, the proposed problem variables are Binary, which is shown in Constraint (\ref{Const9}). Finally, we need Constraint (\ref{Const10}) in practice to remove the links with a very short lifetime. How short the link could be to be discarded depends on the metric chosen for link lifetime to be second, millisecond, etc.

\subsection{Optimization Problem Solution}

Generally, Binary LP problems are in the NP-complete category \cite{karp}, and hence, they do not have a polynomial-time solutions. In this section, we show that the LB-OPAR problem is a special case of Binary LP problems where the optimal solution can be found via graph-based algorithmic solution. Needless to say, the proposed algorithmic solution might not be applicable to other Binary LP problems. While the solution is as follows, Algorithm (\ref{alg::LB-OPAR}) represents the pseudocode of the algorithm. We prove through Lemma (\ref{lem::opt}) that the output of Algorithm (\ref{alg::LB-OPAR}) is the optimal solution of the Binary LP problem proposed in section (\ref{sec::optimization}). The computation and space complexity of Algorithm (\ref{alg::LB-OPAR}) is also investigated in Lemma (\ref{lem::complexity}).

To find the solution to the proposed optimization problem, i.e. the optimal path, we consider the network as a graph $G(V,E)$ where $V$ and $E$ are the sets of graph vertices and edges, respectively. We consider any UAV node as a vertex in graph $G$. Any two neighbor UAVs which are located at the communication range of one another are shown by an edge in graph $G$. Thus, the edge $e_{(i,j)}$ in graph $G$ represents a physical link between UAV nodes $u_i$ and $u_j$. If the communication range between two nodes is not identical, we can consider $e_{(i,j)}$ as a directional edge starting from node $u_i$ toward node $u_j$, without loss of generality.

\begin{algorithm}
\caption{Finding the solution of the optimization problem}
 \label{alg::LB-OPAR}
\begin{algorithmic}
\STATE{$objValue=\infty$}, \COMMENT{Set the initial objective value equal to $\infty$.}
\STATE{$path=\emptyset$}, \COMMENT{Set the initial path to empty path.}
\STATE{$\varsigma\leftarrow w_2\frac{1}{\mathcal{T}}+w_3\mathcal{L} $}, \COMMENT{Form $\varsigma$, a $n\times n$ matrix containing the lifetime-load cost of all graph edges. The lifetime-load cost for each edge $e_{(i,j)}$ is calculated as $\varsigma_{(i,j)}=w_2\frac{1}{\tau_{(i,j)}}+w_3l_j$.  }
\STATE{$\varsigma_{Sorted}=descentSort(\varsigma)$}, \COMMENT{Sort the lifetime-load cost matrix in descending order.}
\STATE{$G'(V,E)=G(V,E)$}, \COMMENT{Let $G'(V,E)$ equal to the main graph $G(V,E)$.}
\WHILE{$isPath(G'(V,E),src,dst)$}
\STATE{$newPath=BFS(G'(V,E),src,dst)$}, \COMMENT{To find the shortest path with BFS algorithm.}
\IF{$newObjValue < objValue$}
\STATE{$path=newPath$},\COMMENT{Set the path equal to the new path.}
\STATE{$objValue=newObjValue$},\COMMENT{Set the objective value equal to the objective value of the new path.}
\ENDIF
\STATE{$\varsigma_{path}	\leftarrow lifetimeLoadCost(newPath)$}, \COMMENT{Calculate the maximum path lifetime-load cost for the edges forming the $newPath$.}
\STATE{$G'(V,E) \leftarrow remove(\varsigma_{path}, \varsigma_{Sorted},G'(V,E))$}, \COMMENT{Remove all links with the lifetime-load cost greater than or equal to the $\varsigma_{path}$.}
\ENDWHILE

\end{algorithmic}
\end{algorithm}

To find the optimal path, we first define a lifetime-load cost metric for each edge. We name this metric $	\varsigma_{(i,j)}$ for edge $e_{(i,j)}$ and we can calculate this metric as $	\varsigma_{(i,j)}=w_2\frac{1}{\tau_{(i,j)}}+w_3l_j$. The basic idea of Algorithm (\ref{alg::LB-OPAR}) is to use the breadth-first search (BFS) algorithm \cite{bfs} to find the shortest path in graph $G(V,E)$, regardless of its lifetime or load. BFS is the algorithm with the lowest computational complexity to find the shortest path in a graph. We then calculate the objective value for the path and the lifetime-load cost, i.e. $	\varsigma_{(i,j)}$, of all edges forming the path. Next, we remove all graph edges with less or equal lifetime-load cost than the maximum $\varsigma_{(i,j)}$ of the path. We apply the BFS algorithm again on the new graph. If the new path has a lower objective value, we consider it as the optimal path; elsewhere, we discard the new path. The procedure is repeated until no more paths stay between the source and the destination nodes. Algorithm (\ref{alg::LB-OPAR}) shows the detailed pseudocode of the proposed algorithm. We show in Lemma (\ref{lem::opt}) that the output of Algorithm (\ref{alg::LB-OPAR}) is the optimal path. We then calculate the space and computational complexity of this algorithm in Lemma (\ref{lem::complexity}).

\begin{lemma}
\label{lem::opt}
If there is more than one path from the source node toward the destination, Algorithm (\ref{alg::LB-OPAR}) returns the optimal path.
\end{lemma}
\begin{proof}
We prove this lemma using proof by induction. If there is at least one path, we know that the BFS algorithm will find it at the first iteration of the algorithm, and it is the optimal path. Assume that the algorithm runs $(n+1)$ iterations to return the final answer. Again, assume that after $n^{th}$ iteration, path $p$ is returned as the optimal path, and it is optimal until then. To prove the lemma, we need to show that the output of the last iteration, i.e. $(n+1)^{th}$ iteration, will be optimal, too. In this case, we have to show that the link(s) that we removed from the graph in $n^{th}$ iteration did not lead to discarding the optimal path. We use proof by contradiction to prove it. Hence, we assume that there exists a path $\bar{p}$ with an objective value less than that of $p$, which at least one of its links is discarded in $n^{th}$ iteration. It is worth mentioning that, since Algorithm (\ref{alg::LB-OPAR}) first finds the shortest path and then calculates the lifetime-load cost and discards some links, $\bar{p}$ cannot be shorter than $p$. If one of $\bar{p}$ links is discarded, it means that the link has a lifetime-load cost more than or equal to the maximum lifetime-load cost of $p$'s links. When $\bar{p}$ is not shorter than $p$, all of its links have to have lifetime-load costs less than that of the maximum cost of $p$'s links to become the optimal path. We meet a contradiction; hence, the lemma is proved. If there exists more than one path from the source node toward the destination, Algorithm (\ref{alg::LB-OPAR}) returns the optimal path.
\end{proof}

\begin{lemma}
\label{lem::complexity}
The computational complexity of Algorithm (\ref{alg::LB-OPAR}) is $O(|E|^2)$ for connected networks, i.e. there is at least one path between any pair of nodes. This algorithm's space complexity is also $O(|V|+|E|)$.
\end{lemma}
\begin{proof}
We know that the computational complexity of the BFS algorithm is $O(|V|+|E|)$ where $|V|$ and $|E|$ stand for the number of graph vertices and edges, respectively \cite{bfs}. In Algorithm (\ref{alg::LB-OPAR}), in the worst case, we need to perform BFS $|E|$ times. We further need to sort the lifetimes-load cost matrix. The worst-case complexity of sorting the lifetime-load cost matrix is $O(|E|\log|E|)$. Hence, the entire algorithm complexity is $O(|E|(|V|+|E|)+|E|\log|E|)=O(|E|^2+|E||V|+|E|\log|E|)$. We know that $O(|E|)\ge O(|V|)$ for the connected networks. Hence, the computational complexity of Algorithm (\ref{alg::LB-OPAR}) is $O(|E|^2)$.
We further know that the space complexity of the BFS algorithm is $O(|V|)$. We need further $O(|E|)$ space to store the lifetime-load cost matrix. The space complexity of the sort is $O(1)$. Hence, the space complexity of Algorithm (\ref{alg::LB-OPAR}) is $O(|V|+|E|)$.
\end{proof}

It is worthy to mention that we have two reasons in choosing our LP-based methodology instead of workflow scheduling techniques such as the reference one multi-objective heterogeneous earliest-finish-time  (MOHEFT) \cite{MOHEFT}, and the much improved one FAst  workflow  scheduling  approach based  MOHEFT  using BAcktraking and CHeckpointing  (FAMOBACH) \cite{FAMOBACH}. Workflow scheduling techniques, in general, work with much lower computational complexity than LP algorithms and lead to near-optimal solutions. However, our proposed solution finds the exact optimal answer to the problem as we proved it in Lemma (\ref{lem::opt}),  and the computational complexity of our solution is $O(|E|^2)$ at its worst case, as we shown in Lemma (\ref{lem::complexity}), which is fairly low in comparison with the mentioned alternatives.

\subsection{Link Lifetime Prediction and Node Load Calculation}
\label{sec::LLTprediction}

The optimization problem of Section (\ref{sec::optimization}) considers a $n\times n$ matrix $\mathcal{T}$ which contains link lifetimes as well as vector $\mathcal{L}$ with $n$ elements which contain nodes' load. In this section, we describe how we predict the link lifetimes and calculate node load to form $\mathcal{T}$ and $\mathcal{L}$. We start with $\mathcal{T}$ where we need only three consecutive node locations to fill it. Each node regularly sends a packet to the controller, including its last three consecutive locations in 3D space. The controller, in its turn, calculates the azimuthal and polar directions of each node as well as their speed and acceleration. Using this information, the controller calculates the time for each pair of neighboring nodes that they go out of the communication range of each other. To use simpler notations and without the loss of generality, we consider the communication range of all nodes equal and represent it with $\Re$. It is worthy to mention that our predicted lifetimes stay precise if the corresponding nodes do not change their direction or acceleration.

We represent the location of node $u_i$ at time $t$ by the tuple $(x_i(t),y_i(t),z_i(t))$. We show the distance traversed by the node $u_i$ in the time interval $[t_1 \quad t_2]$ by $d_{(t_1,t_2)}^{(i)}$ and it could be easily calculated by Equation (\ref{eq::distance}).

\begin{equation}
\label{eq::distance}
\scriptstyle
    d_{(t_1,t_2)}^{(i)}=\sqrt{[x_i(t_2)-x_i(t_1)]^2 + [y_i(t_2)-y_i(t_1)]^2 + [z_i(t_2)-z_i(t_1)]^2}
\end{equation}

Considering the three consecutive node locations at times $t_0$, $t_1$, and $t_2$, we calculate the azimuthal angle $\alpha_i$ and polar angle $\theta_i$ of the node's movement direction using Equations (\ref{eq::azimuthal}) and (\ref{eq::polar}), respectively.

\begin{equation}
\label{eq::azimuthal}
\alpha_i=tan^{-1}({\frac{y_i(t_2)-y_i(t_1)}{x_i(t_2)-x_i(t_1)}})\\
\end{equation}

\begin{align}
\label{eq::polar}
\theta_i&=cos^{-1}({\frac{z_i(t_2)-z_i(t_1)}{d_{(t_1,t_2)}^{(i)}}})\\\nonumber
&=tan^{-1}({\frac{\sqrt{[x_i(t_2)-x_i(t_1)]^2 + [y_i(t_2)-y_i(t_1)]^2 }}{z_i(t_2)-z_i(t_1)}})
\end{align}

We can calculate speed $v_i$ and acceleration $a_i$ of node $u_i$ respectively using Equations (\ref{eq::3dSpeed}) and (\ref{eq::acc}).

\begin{equation}
\label{eq::3dSpeed}
v_i(t_1)= \frac{\sqrt{d_{(t_0,t_1)}^{(i)}}}{t_1-t_0}
\end{equation}

\begin{equation}
\label{eq::acc}
a_i(t_2)=\frac{v_i(t_2)-v_i(t_1)}{t_2-t_0}
\end{equation}

Next, we can calculate the predicted location of the node $u_i$ at time $(t_2+\Delta t)$ using Equation (\ref{eq::3dNextPosAcc}). Then, we can calculate the distance between nodes $u_i$ and $u_j$ at time $(t_2+\Delta t)$ using Equation (\ref{eq::nodeDistance}). We represent this distance by $d_{(i,j)}(t_2+\Delta t)$ where $\Delta t$ is unknown. Knowing the communication range of nodes, i.e. $\Re$, we can calculate the time $\Delta t$ using Equation (\ref{eq::root}). This variable is the time required for the nodes to go out of the communication range of one another and represents the lifetime of the link $e_{(i,j)}$. Equation (\ref{eq::root}) is a simple nonlinear equation with only one unknown variable, i.e. $\Delta t$. This equation can be easily solved with any lightweight numerical method such as the Bisection method. It is worth mentioning that the LB-OPAR is a modular routing algorithm. Hence, the prediction algorithm proposed here, while its accuracy is validated in Section (\ref{sec::predEval}), could be replaced with any other prediction algorithm with no need for any modification in the overall routing
solution.

{\small
\begin{align}
\label{eq::3dNextPosAcc}
&\begin{cases}
x_i(t_2+\Delta t)=x_i(t_2)+( v_i(t_2) \Delta t+\frac{1}{2}a_i \Delta t^2 ) sin(\theta_i) cos(\alpha_i)\\
y_i(t_2+\Delta t)=y_i(t_2)+(  v_i(t_2)\Delta t+\frac{1}{2}a_i \Delta t^2 ) sin(\theta_i) sin(\alpha_i)\\
z_i(t_2+\Delta t)=z_i(t_2)+( v_i(t_2)\Delta t+\frac{1}{2}a_i \Delta t^2 ) cos(\theta_i)
\end{cases}
\end{align} }

\begin{equation}
\label{eq::nodeDistance}
d_{(t_2+\Delta t)}^{(i,j)}=\\
\sqrt{ \begin{aligned}
&(x_i(t_2+\Delta t)-x_j(t_2+\Delta t))^2+\\
&(y_i(t_2+\Delta t)-y_j(t_2+\Delta t))^2+\\
&(z_i(t_2+\Delta t)-z_j(t_2+\Delta t))^2
\end{aligned}
}
\end{equation}
\begin{equation}
\label{eq::root}
d_{(t_2+\Delta t)}^{(i,j)}-\Re=0
\end{equation}

In a wireless network, a node can start transmitting only when the channel is empty, which happens when no other node in the communication range of the sending node transmits a packet. If, for any reason, more than one node attempt to send a packet simultaneously, a collision will happen. Generally, wireless networks, via their medium access control (MAC) protocols, try to avoid the collision by preventing one of the nodes from sending its own packets at the same time. This fact leads to a delay in flow completion time. In our SDN-based system, the controller is aware of all transmitting flows and their paths. Hence, it can calculate the number of nodes that may attempt to send a packet in the communication range of a sending node. We refer to this number as the node's load.

To calculate each node's load, the controller starts the network operation with all zeros for the nodes' load. Any time the controller receives a route request, it calculates the optimal path and determines the nodes affected by the new route. Any node included in the new route and all its one hop neighbors are considered as the nodes affected by the new route. Since the controller potentially knows all nodes' locations, it can also calculate the distance between each pair of nodes, and it further knows the communication range of each node. Hence the controller can determine the affected nodes by the new route. The controller then increases the load of the determined nodes by one. Every time that a route becomes unavailable, the controller decreases the affected nodes' load by one. In two cases, a route is considered unavailable. First, the source node does not receive the TCP acknowledgment packets, and consequently, the source node orders a reroute. Second, the source node receives the acknowledgment of the last packet of the file and sends a flow termination message to the controller.

\section{Performance Evaluation}
\label{sec:NumericalResults}
To evaluate the performance of LB-OPAR, we use ns-3 network simulator. Ns-3 is a well-known discrete-event simulator widely used for research and development \cite{ns3}. We simulate the network for different scenarios with a different number of nodes and a different number of concurrent file-transfer flows to exhaustively evaluate the network performance in different situations. Network throughput, flow completion time, and flow success rate are measured as the performance evaluation metrics. To show the significant improvement of LB-OPAR in the network performance, we compare its results with the baseline benchmarks AODV, OLSR, DSDV as well as the OPAR algorithm. We select the baseline algorithms that cover the proactive and reactive, link-state and distance-vector, as well as SDN-based and fully distributed routing algorithms. In this section, we first describe the simulation setting. Next, we evaluate the prediction algorithm proposed in Section (\ref{sec::LLTprediction}). Since LB-OPAR is sensitive to its weights $w_1$, $w_2$, and $w_3$, we then perform an analysis on these weights to find the proper value for them in different network settings. Finally, we represent the results of the comparison among LB-OPAR, AODV, DSDV, OLSR, and OPAR, along with the analysis.

\subsection{Simulation Setting}

We simulate a UAV network with $n$ nodes in a $2000\times 300\times 50$ 3D space. We perform two sets of simulation scenarios. In the first set of the simulation scenarios, we vary the number of nodes from 50 to 100 in an increment step of 5. In this set of scenarios, five source nodes start sending their data to five corresponding destination UAVs. The generated traffic is FTP over TCP NewReno. A 5 MB data file is sent in each flow. In the second set of scenarios, we set the number of nodes equal to fifty, i.e. $n=50$, and vary the number of concurrent flows from one to ten. Each simulation instance is run for 500 seconds. UAV nodes move based on two different random movement patterns, 3D random waypoint, and 3D Gauss-Markov model. While random waypoint chooses its waypoints in a completely random manner, the Gauss-Markov model is closer to pre-planned scenarios. It prevents the sudden significant direction change by applying a limitation on the variability of the movement angle. We choose these two mobility models to show that LB-OPAR shows excellent performance in both random and pre-planned scenarios. The speed range is set to $[0 \quad 50]$ $m/s$. We run each simulation instance 20 times with different random seeds, measure the parameters of interest, and report the average results.

We set the transmission power to 7.5 dBm for each UAV. We use the free-space propagation loss model, also known as Friis. The propagation loss of UAV transmission is very well-known to behave close to the Friis model. The constant speed propagation model is also selected as the propagation delay model. UAV nodes use IEEE 802.11b as their wireless communication standard. Table (\ref{tbl::simSetting}) represents the simulation parameters.

\begin{table}[t]
\caption[]{Simulation Setting }
\resizebox{1\textwidth}{!}{
\begin{minipage}{\textwidth}
\begin{tabular}{ l | l  }
 Simulator version & ns-3 3.30\\
  Number of UAVs & $[50\quad 100]$ \\
  Area size & $300\times 2000\times 50 m$\\
  Transmission power & 7.5 dBm\\
  Number of concurrent flows & $[1 \quad 10]$\\
  File size & 5 MB\\
  Simulation time & 500 sec\\
  Speed range & $[0 \quad 50] m/s$\\
  Mobility models & 3D RWP\\
  Traffic type & TCP NewReno\\
  Wireless communication standard & IEEE 802.11b\\
  Propagation loss model & Free-space propagation loss\\
  Propagation delay model & Constant speed propagation delay
\end{tabular}
\label{tbl::simSetting}
\end{minipage}}
\end{table}

We measure the flow success rate, network throughput, and flow completion time to make the comparison. We define flow success rate as the rate of the flows that successfully deliver their entire 5 MB file to the corresponding destination node in the simulation time. It is worth mentioning that if the simulation time is infinite, all flows with any of the routing algorithms can deliver their files successfully. However, practically it is not feasible to keep the simulation without stop time. Hence, if a flow destination does not receive all of the 5 MB file's packets in the simulation time, we consider that flow as unsuccessful. In the real world, each communication flow may have a deadline. The simulation time is the mapping of the communication deadline concept in the simulation world.

We use the standard definition for throughput, which is the rate of successful packet delivery over the communication bandwidth, and it is measured in Mbps. While we report the average throughput for all successful flows, we find that it might not be fair to report the average of raw throughput values. In some cases, a specific routing algorithm works worse than the others and fails in finding the proper routes. As an instance, it may successfully deliver just one file out of ten. Assume that the source node and the destination be physical neighbors for the successful flow. In this case, it achieves a close to $1$ throughput, where the other routing algorithms, which may have a perfect flow success rate, will show much lower throughput. Hence, to report a fair comparison, we weight the throughput with its corresponding flow success rate.

Finally, we define FCT as the time between the first packet of the 5 MB file sent by the source node and the last packet of the same file received by the destination node. Since the simulation time is 500 seconds, some flows may fail in delivering all of their packets to the corresponding destination in the simulation time. Hence, we consider the FCT equal to 500 for the unsuccessful flows as the lower bound for the FCT of those flows.

\subsection{Prediction Evaluation}
\label{sec::predEval}

To evaluate the prediction algorithm presented in Section (\ref{sec::LLTprediction}), we design a simulation scenario and compare its results with those of the extrapolation algorithm as a standard baseline. We use Newton divided differences extrapolation \cite{newton} to predict node's position in time $t_2+\Delta t$, using three consequent node positions in times $t_0$, $t_1$, and $t_2$. It is worthy to note that all other standard extrapolation algorithms, such as Lagrange extrapolation, result in the same value of the Newton divided differences method \cite{numerical}. We choose this algorithm for its better computation complexity in comparison with the others. Next, we put the predicted position in Equation (\ref{eq::nodeDistance}) to calculate the distance between $u_i$ and $u_j$ at time $t_2+\Delta t$, and form Equation (\ref{eq::root}). Then, we use a lightweight numerical method such as Bisection to find the root and equivalently predict the link lifetime between the corresponding nodes. We measure the prediction error of the link lifetime as the absolute difference between the predicted lifetime and the simulated one.

We consider a network with a different number of nodes where the nodes start moving for a warmup process of 100 seconds. We perform the prediction and then continue the simulation for extra 500 seconds. We predict the link lifetimes using both the algorithms of Section (\ref{sec::LLTprediction}) and the extrapolation. We then watch the links to measure their simulated lifetimes. Next, We calculate the expected value of the error for both methods in each scenario. Obviously, we predict the lifetime for the links that already existed at the time of prediction. We consider the lifetime of the links that last alive until the end of the simulation as the remaining simulation time, i.e. 500 seconds. We further consider the lifetime of the nodes in which the prediction returned a value more than the simulation time, equal to the remaining simulation time. We run each scenario 1000 times and calculate the expected value of the error and its standard deviation. While the expected value of the error shows the accuracy of the prediction algorithm, the standard deviation shows its precision. The precision represents the interval of the error distribution around the expected value. We find that the curves of all network densities behave almost identical. Hence, we report only the results for a network with 50 nodes.

Fig. (\ref{fig::accuracy}) shows the cumulative distribution function (CDF) for the results of the expected value and the standard deviation of the error for both prediction algorithms presented in Section (\ref{sec::LLTprediction}) and the one use extrapolation. Fig. (\ref{fig::accuracy}a) shows that the prediction algorithm proposed here is on average about 70 seconds more accurate than the one use extrapolation. Fig. (\ref{fig::accuracy}b) shows the CDF of the standard deviation around the expected value, which represents the prediction precision. The lower standard deviation value shows higher precision, and consequently, a more reliable prediction algorithm. This figure shows that the standard deviation for the proposed prediction algorithm is averagely around 20 seconds, where this value for the extrapolation algorithm is about 90 seconds. Regardless of the accuracy of the algorithm, the standard deviation shows that the proposed algorithm is more reliable than that of extrapolation since it has fewer fluctuations.

\begin{figure}[t!]
	\centering
	\subfloat[Error expected value ]{\includegraphics[width=.5\linewidth]{figs/predictionAccuracy}}
	\subfloat[Error standard deviasion]{\includegraphics[width=.5\linewidth]{figs/predictionAccuracySD}}\\
	\caption{A comparison of CDF of LB-OPAR prediction algorithm versus extrapolation.}
	\label{fig::accuracy}
\end{figure}

\subsection{Weight Analysis}

In this Section, we analyze the path length, path lifetime, and path load weights of the LB-OPAR optimization problem, i.e. $w_i$s. Accordingly, we designed a simulation scenario with the same setting of Table (\ref{tbl::simSetting}) where we vary the weights to measure the sensitivity of the network performance versus this variation. In this simulation, we first analyze the path load weight, i.e. $w_3$. We vary the value of $w_3$ from 0 to 1 in an increment step of $0.1$ and let $w_1=w_2$. Since $\sum_{i=1}^{3}w_i=1$, the value of $w_1$ and $w_2$ could be calculated easily.

Simulation results show that the highest performance in terms of throughput, success rate, and FCT is achieved for the values of $w_3$ mentioned in Table (\ref{tbl::w3}). The results show that for the network with less load, assigning a value more than zero to $w_3$ might not benefit the network performance. However,  increasing the number of concurrent flows needs a higher weight for $w_3$ to achieve better performance in terms of all mentioned performance metrics. Furthermore, we find that there is no optimal case with $w_3=1$ which means path lifetime and path length always have an obvious effect in choosing the optimal path. We also find that increasing the network density has no significant effect on the value of $w_3$.

\begin{table}[t]
\caption[]{Path load weight analysis}
\resizebox{\columnwidth}{!}{
\begin{minipage}{\textwidth}
\begin{tabular}{ | l || l | l | l | l | l | l | l | l | l | l | l | }
\hline
 no. of Flows & 1& 2 &3 &4 &5 &6&7&8&9&10\\\hline
 $w_3$& 0& 0 &0 &0.1 &0 &0.3&0.4&0.7&0.6&0.7\\\hline\hline
  no. of UAVs & 55& 60 &65 &70 &75 &80&85&90&95&100\\\hline
 $w_3$        & 0.1& 0 &0.2 &0.1 &0.2 &0&0&0&0&0.1\\\hline
\end{tabular}
\label{tbl::w3}
\end{minipage}}
\end{table}

Next, we analyze the effect of varying path-length weight, i.e. $w_1$, versus path lifetime weight, i.e. $w_2$, for the fixed path load weight represented in Table (\ref{tbl::w3}). By fixing the value of $w_3$ and varying the values of $w_1$ and $w_2$ we find that the best achievable performance is for the values represented in Tables (\ref{tbl::wsCon}) and (\ref{tbl::wsNode}) for different network loads and different network densities, respectively. Table (\ref{tbl::wsCon}) shows that for the network with less number of concurrent communication flows, i.e. less traffic, the path lifetime weight, i.e. $w_2$ is more important in comparison with path length. However, results show that the weights of path lifetime and path length are almost equal for the network with higher traffic. Table (\ref{tbl::wsCon})  also shows that the network density has no significant effect on the values of path load and path lifetime weights. Hence, for the network with five concurrent communication flows and different numbers of UAVs, the best performance is achievable when $w_1=w_2$, in most cases.

\begin{table}[t]
\caption[]{Weight Analysis for Different Network Loads}
\resizebox{\columnwidth}{!}{
\begin{minipage}{\textwidth}
\begin{tabular}{ | l || l | l | l | l | l | l | l | l | l | l | l | }
\hline
 no. of Flows & 1& 2 &3 &4 &5 &6&7&8&9&10\\\hline
 $w_3$& 0& 0&0 &0.1 &0 &0.3&0.4&0.7&0.6&0.7\\\specialrule{.1em}{.05em}{.05em}
 $w_1$& 0.3 & 0.3    &0.3 &0.2  &0.5   &0.4   &0.4  &0.15   &0.2  &0.15\\
 $w_2$& 0.7 & 0.7 &0.7 &0.7 &0.5 &0.3&0.2&0.15&0.2&0.15\\\hline
  \end{tabular}
\label{tbl::wsCon}
\end{minipage}
}
\end{table}

\begin{table}[t]
\caption[]{Weight Analysis for Different Network Densities}
\resizebox{\columnwidth}{!}{
\begin{minipage}{\textwidth}
\begin{tabular}{ | l || l | l | l | l | l | l | l | l | l | l | l | l | }
\hline
 no. of UAVs & 50& 55 &60 &65 &70 &75 & 80 &85 &90&95&100\\\hline
 $w_3$& 0& 0.1 &0 &0.2 &0.1 &0.2&0&0&0&0 &0.1\\\specialrule{.1em}{.05em}{.05em}
 $w_1$& 0.5 & 0.45    &0.5 &0.4  &0.6   &0.4   &0.5  &0.4   &0.4  &0.5 &0.45\\
 $w_2$& 0.5 & 0.45 &0.5 &0.4 &0.3 &0.4&0.5&0.6&0.6&0.5 &0.45\\\hline
  \end{tabular}
\label{tbl::wsNode}
\end{minipage}
}
\end{table}

\subsection{Simulation Results}

In this section, we compae the LB-OPAR simulation results with those of OPAR, AODV, DSDV, and OLSR, for both RWP and G-M mobility models. We compare the results for flow success rate, network throughput, and flow completion time for different network densities and different network loads, i.e. concurrent flows.

Fig. (\ref{fig::success}) shows the results for the flow success rate. DSDV and OLSR routing protocols show less than $10\%$ of flow success rate, in most cases, which is mainly because of their inability in catching up with the
high dynamicity of the nodes' movement in 3D space. AODV, however, represents better performance in terms of flow success rate in comparison with DSDV and OLSR. It successfully delivers $30\%$ of the flows, in average. Indeed, both OPAR and LB-OPAR significantly outperform the other algorithms by the average of $60\%$ flow success rate. As it is obvious, increasing the network load leads to more significant performance improvement for LB-OPAR compared to OPAR by the average of $5\%$. For the less number of concurrent flows, since the network load weight is zero or close to zero, no significant difference is shown between LB-OPAR and OPAR, as it is the case for the network with 5 concurrent flows and different densities, i.e. Fig. (\ref{fig::success}c) and (\ref{fig::success}d). Something worth mentioning is that AODV, DSDV, and OLSR show a higher success rate under the RWP mobility model than G-M. It seems that the randomness of the RWP model helps them to increase their performance. A similar fact was previously discovered in different network studies such as \cite{jellyfish}. OPAR and LB-OPAR show no significant difference under RWP and G-M mobility models, representing their stability of performing well in different conditions.

\begin{figure}[t!]
	\centering
	\subfloat[Different number of  flows (RWP)]{\includegraphics[width=.5\linewidth]{figs/successRateConRWP}}
	\subfloat[Different number of  flows (G-M)]{\includegraphics[width=.5\linewidth]{figs/successRateConGM}}\\
	\subfloat[Different number of UAVs (RWP)]{  \includegraphics[width=.49\linewidth]{figs/successRateNodeRWP}}
	\subfloat[Different number of UAVs (G-M)]{  \includegraphics[width=.49\linewidth]{figs/successRateNodeGM}}
	\caption{A comparison of the flow success rate for different loads and densities.}
	\label{fig::success}
\end{figure}

Next, we investigate the network throughput for all the combinations of different network loads, network densities, and mobility models. Fig. (\ref{fig::throughput}) shows the results where the achieved throughput by DSDV and OLSR is close to zero due to their failure in handling the fast topology changes. Interestingly and in contrast to the conventional mobile networks, increasing the number of concurrent flows shows no significant effect on the throughput, which means in highly dynamic networks, the bottleneck of the throughput is the reroute process due to the topology change not the load of the network. For this parameter, AODV
shows much better performance in comparison with DSDV and OLSR, while OPAR and LB-OPAR outperforms AODV by $100\%$ improvements. OPAR and LB-OPAR show competitive behavior. However, LB-OPAR outperforms OPAR in higher network loads and under the RWP mobility model for about $20\%$. For different network densities, since there are only five concurrent flows and the network load weight is zero or close to zero, we see no significant differences between LB-OPAR and OPAR. However, they outperform AODV by the average of $100\%$ improvements.

\begin{figure}[t!]
	\centering
	\subfloat[Different number of  flows (RWP)]{\includegraphics[width=.5\linewidth]{figs/throughputConRWP}}
	\subfloat[Different number of flows (G-M)]{\includegraphics[width=.5\linewidth]{figs/throughputConGM}}\\
	\subfloat[Different number of UAVs (RWP)]{  \includegraphics[width=.49\linewidth]{figs/throughputNodeRWP}}
	\subfloat[Different number of UAVs (G-M)]{  \includegraphics[width=.49\linewidth]{figs/throughputNodeGM}}
	\caption{A comparison of the throughput for different loads and densities.}
	\label{fig::throughput}
\end{figure}

Fig. (\ref{fig::FCT}) shows the results for flow completion time. As we mentioned earlier, we consider the simulation time as the lower bound of the FCT for the failed flows. Hence, the results of Fig. (\ref{fig::FCT}) show the lower bound for the different algorithms' FCT. This figure includes the combination of all network densities, loads, and the simulated mobility models. For all simulation settings,  OPAR and LB-OPAR show the lowest FCT compared to the other routing algorithms. However, LB-OPAR for the network with higher loads shows around $15\%$ less FCT in comparison with OPAR, under RWP mobility model. OPAR and LB-OPAR outperforms AODV by averagely $10\%$.  AODV, in its turn, shows a significantly lower FCT than that of OLSR and DSDV where their average FCT is close to the simulation time, i.e. 500 s, due to theire high failure rate. In many cases, OPAR and LB-OPAR show close FCT results, which is one of two cases. The first case is when there are one or  a few number of concurrent flows which lead to negligible effect of network load in choosing the optimal path. The latter case owing to the fact  that in some cases, there is only one path, or there are several paths, but just one of them has a reasonable objective value, where the others are very long or with a very short lifetime. Hence, both the OPAR and LB-OPAR select the same path, which leads to the close results.

\begin{figure}[t!]
	\centering
	\subfloat[Different number of  flows (RWP)]{\includegraphics[width=.5\linewidth]{figs/FCTConRWP}}
	\subfloat[Different number of flows (G-M)]{\includegraphics[width=.5\linewidth]{figs/FCTConGM}}\\
	\subfloat[Different number of UAVs (RWP)]{  \includegraphics[width=.49\linewidth]{figs/FCTNodeRWP}}
	\subfloat[Different number of UAVs (G-M)]{  \includegraphics[width=.49\linewidth]{figs/FCTNodeGM}}
	\caption{A comparison of FCT for different loads and densities.}
	\label{fig::FCT}
\end{figure}

\section{Acknowledgement}
This material is based upon work supported by the Air Force Office of Scientific Research under award number FA9550-20-1-0090 and the National Science Foundation under Grant Number CNS-2034218. Any opinions, findings and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the US government or AFRL.

\section{Conclusions and Future Work}
\label{sec::conclusion}

The highly dynamic nature of cooperative UAV networks causes the conventional routing algorithms, which aim at finding the shortest path, to fail in reaching the optimal network performance. While we have shown, in our recently published work OPAR, that considering the path lifetime in selecting the routes has a significant impact on the network performance, in this paper, we extend OPAR to consider the network load as well. Accordingly, in this paper, we proposed load-balanced OPAR, an SDN-based routing algorithm for cooperative UAV networks. We modeled the entire problem with an analytical model and proposed a lightweight solution for the proposed model. We implemented LB-OPAR in the ns-3 network simulator and exhaustively evaluated its performance. We show that LB-OPAR outperforms the benchmark routing algorithms AODV, DSDV, and OLSR. It also improves the performance of OPAR when the network works under higher loads. We further found that the performance bottleneck in highly dynamic UAV networks is basically the route breaks down, more than any other reason. Hence, proposing a routing algorithm to cover the multiple reroutes problem has significant importance in improving cooperative UAV networks' performance. Since LB-OPAR has a lightweight computational and space complexity, we aim at distributing the SDN controller tasks among all network nodes for fully distributed networks as future work.

\nocite{*}
\bibliographystyle{IEEEtran}


\section{Introduction}
Rational design of energy conversion materials are indispensable to develop cutting edge technologies to exploit the renewable energy resources for future thermal energy management systems. Thermoelectric (TE) technology offers an environmentally benign route for efficient conversion of waste heat into useful electricity. However, low performance TE materials turns out to be a bottleneck for designing high performance TE devices. Performance of a thermoelectric material is determined by a dimensionless quantity, figure of merit, ZT = S$\sigma^2$T/$\kappa$, where S, $\sigma$, and T are the Seebeck coefficient, electrical conductivity, and absolute temperature, respectively. $\kappa$ = $\kappa_e$+$\kappa_l$; is total thermal conductivity which includes contribution from electronic ($\kappa_e$) and lattice ($\kappa_l$) thermal conductivities. Therefore, discovery of advanced functional materials through various design strategies or fine-tuning the existing TE materials properties by an external parameter such as doping, pressure, temperature, strain, etc., to improve the performance of a material. In addition, lone pair can be used as a tool for designing many classes of functional materials\cite{Seshadri2022} in general and low $\kappa_l$ materials in particular because of the lone pair induces strong anharmonicity in the system. On the other hand, (layered) structure causes strong bonding hierarchy $i.e.,$ having rigid and fluctuating sub-lattices in a material has important implications to enhance complex interdependent thermoelectric properties, for instance, high power factor (PF=S$^2\sigma$) and low lattice thermal conductivity ($\kappa_l$).\cite{Pal2018,pandey2020ultralow}

Bismuth oxychalcogenides are emerging class of materials due to their outstanding thermeoelectric and opto-electronic properties.\cite{Tippireddy2021} In particular, dibismuth dioxychalcogenides, Bi$_2$O$_2$Ch (Ch = S, Se, Te), their 2D counterparts,\cite{Yang2019,Li2021,Liang2021} and heterostructures\cite{Bo2022} gained tremendous research interest due to their potential applications in thermoelectrics,\cite{Luu2015,Wang2021,Tippireddy2021} ferroelectrics,\cite{Wu2017} and optoelectronics.\cite{Yang2022,Tippireddy2021,Li2021,Zhang2015,Guo2019} The ternary Bi$_2$O$_2$Ch systems are synthesized by a solid state reaction of well-known binary thermoelectric bismuth chalcogenides Bi$_2$Ch$_3$ (Ch = S, Se, Te) with Bi$_2$O$_3$. Bi$_2$O$_2$S crystallizes in the primitive orthorhombic structure (SG: $Pnmn$, Z=2 formula units (f.u.) per primitive cell), while Bi$_2$O$_2$(Se/Te) crystallizes in the body centred tetragonal (SG : $I4/mmm$, Z = 2 f.u./cell) anti-ThCr$_2$Si$_2$-type structure by having the Bi$_2$O$_2$ (PbO-type structural units in a square planar unit of oxygen atoms mono-caped with Bi atom) and  alternating layers of Se, which are held together by weak electrostatic forces in contrast to van der Waals (vdW) interactions that are generally observed in layered materials.\cite{Wei2019,Guo2021} The stereo-chemical activity of 6$s^2$ lone-pair of Bi$^{3+}$ cation is the main reason for different crystal symmetry of Bi$_2$O$_2$S over heavy chalcogen based Bi$_2$O$_2$Se and Bi$_2$O$_2$Te systems. Infra-red and/or Raman spectra\cite{Xu2019} and strain effects on bulk and monolayer Bi$_2$O$_2$Ch crystals are systematically investigated to explore their possible applications in nano-electronics.\cite{Cheng2018} Structural stability and topological properties of ternary Bi$_2$X$_2$Y (X, Y = O, S, Se, and Te) compounds are studied using dispersion corrected DFT-D2 method.\cite{Ma2018} A high pressure X-ray diffraction study\cite{Bu2020} revealed that Bi$_2$O$_2$S undergoes a structural transition from $Pnnm$ $\rightarrow$ $I4/mmm$ at 6.4 GPa due to suppression of the Bi$^{3+}$ cation 6$s^2$ lone-pair and then 2D $\rightarrow$ 3D structure up on further compression above 13.2 GPa owing to disappearance the lone-pair. Bi$_2$O$_2$Se is one of the most stable Sillen-type compounds under pressure with no structural transition being observed until 30 GPa\cite{Pereira2018}. Bi$_2$O$_2$Te is found to be stable till 50 GPa \cite{Wang2020-1} and it has promising applications in optical fields as a tuning material due to large birefringence.\cite{Chen2019} The intrinsic ultralow $\kappa_l$ behavior\cite{Wang2018,Song2019,Song2020,Guo2021} and thermoelectric properties\cite{Wang2018,Song2020} of Bi$_2$O$_2$Ch (Ch = S, Se, Te) compounds are studied using first principles calculations and Boltzmann transport theory. Effect of hydrostatic pressure on stereochemically active lone-pair and its implications on lattice dynamics and phonon transport of Bi$_2$O$_2$S is scarce in the literature, which is essential, as it provides a fundamental understanding to enhance the performance of Bi$_2$O$_2$S in thermoelectrics and/or opto-electronic applications.

Here, we systematically investigate the crystal and electronic structure, elastic properties, lattice dynamics, and phonon transport properties for the ambient ($Pnnm$) and high pressure ($I4/mmm$) phases of Bi$_2$O$_2$S. By analyzing pressure-dependent static enthalpy, lattice, bond parameters, and elastic constants, we show that Bi$_2$O$_2$S undergoes a structural phase transition from low symmetry ($Pnmn$) to a high symmetry ($I4/mmm$) structure around 4 GPa.  By analysis of electron localization function (ELF) and electronic structure, we further show that under pressure the lone-pair electrons of Bi$^{3+}$ cation are suppressed.  This suppression of lone-pair electrons reduces anharmonicity at high pressure and therefore $\kappa_l$ is enhanced.

\section{Computational details and Methodology}
The ab-initio calculations have been performed using Vienna Ab-initio Simulation Package (VASP)\cite{Kresse1996}. The electron-electron interactions $i.e.$ exchange-correlation are treated by Generalized Gradient Approximation (GGA) within the parametrization of Perdew-Burke-Ernzerhof (PBE). The electron-ion interactions are treated with pseudopotential approach. We have used DFT-D3 method\cite{Grimme2011} to capture the weak vdWs interactions. The following plane wave basis orbitals were considered as valence electrons; Bi: 5d$^{10}$,6s$^2$,6p$^3$; O: 2s$^2$,2p$^4$; and S: 3s$^2$,3p$^4$. A kinetic energy cutoff of 600 eV was used for plane wave basis set expansion and a spacing of 2$\pi$ $\times$ 0.025 $\AA^{-1}$ for k-mesh in the irreducible Brillouin zone. The plane wave cutoff energy as well as k-spacing were varied first to ensure the total energy convergence. Self consistency criteria is chosen in such way that the total energy is converged to 1e$^{-8}$ eV/atom and the maximal force on each atom is less than 1e$^{-3}$ eV/$\AA$. Electronic structure is calculated using Tran Blaha modified Becke Johnson (TB-mBJ) potential\cite{Tran2009} as implemented in WIEN2k package\cite{Blaha2002} including spin-orbit coupling. Crystal structures and electron localization functions are visualized using the VESTA software\cite{VESTA2008}.

Harmonic and anharmonic interatomic force constants (IFCs) are needed to calculate $\kappa_l$ from first principle calculations. For both the systems, harmonic IFCs were obtained using the finite displacement method within the Phonopy package~\cite{togo2008first,togo2015first} with 5$\times$5$\times$2 (500 atom) supercell and $\Gamma$-point only k-grid. To obtain well converged phonon frequencies a high plane wave cutoff of 600 eV and strict energy convergence criteria of 1e$^{-8}$ were used. Long range corrections to phonon frequencies were included by calculation of Born effective charges and dielectric constants within the method proposed by Wang et. al.~\cite{wang2010mixed}. Anharmonic IFCs were calculated using 3$\times$3$\times$2 supercell. For the displaced supercells configurations phonon-phonon interactions were considered out to the 11$^{th}$ nearest neighbors. This requires 972 and 500 DFT calculations for $Pnnm$ and $I4/mmm$ phases, respectively. To reduce the computational cost here, only $\Gamma$-point k-mesh along with plane wave energy cutoff of 520 eV and energy convergence criteria of 1e$^{-8}$ eV was used for Brillouin zone integration. The calculated harmonic and anharmonic force constants were used to solve the phonon Boltzmann transport equation iteratively as implemented in the ShengBTE code~\cite{li2012thermal1,li2012thermal,li2014shengbte}. A Gaussian smearing width of 0.1 was used to approximate the Dirac  distribution, and a q-point mesh of 21$\times$21$\times$9 was used for Brillozin zone integration. In these calculations three-phonon scattering and the isotopic disorder scattering~\cite{tamura1983isotope} from natural isotope mass variation are considered.

\section{Results and Discussion}

\subsection{Crystal structure and structural phase transition}
As discussed in section I, at ambient conditions, Bi$_2$O$_2$S crystallizes in lower symmetry ($Pnmn$) structure (see Figures \ref{fig:str}(a-b)) than heavy chalcogen based systems namely Bi$_2$O$_2$Se and Bi$_2$O$_2$Te. This is due to the presence of stereochemically active $6s^2$ lone-pair of Bi$^{3+}$ cation in Bi$_2$O$_2$S and its activity might be suppressed by chemical pre-compression of heavy chalcogens (Se/Te) in Bi$_2$O$_2$Se and Bi$_2$O$_2$Te systems, which drives them to crystallize in relatively high symmetry ($I4/mmm$) structure (see Figures \ref{fig:str}(c-d)) and the $Pnmn$ structure is distorted analogue of $I4/mmm$ crystal symmetry. To investigate the bonding differences between these two phases further, we examine the ELF~\cite{becke1990simple, seshadri2001visualizing} as depicted in Figures~\ref{fig:str}(e)-(g). The ELF is projected onto the [110] plane, which contains atoms of Bi and S. The red color represents complete electron localization, while the blue color represents almost no localization. The ELF of the two phases appear similar at first glance, but significant differences can be seen around Bi-S bonds. In the $Pnmn$ phase, the Bi lone-pair can be seen by yellow regions of electron localization in Figure~\ref{fig:str}(e). When pressure is applied (Figures ~\ref{fig:str}(f)-(g)), the intensity of the yellow lobes surrounding the Bi atoms is significantly reduced, indicating that the lone-pair of Bi$^{3+}$ cation is suppressed. Due to presence of lone-pair in $Pnmn$ phase, the coordination environment around Bi$^{3+}$ cation differs greatly between these two phases. In $Pnmn$ phase, Bi and S have two in-equivalent bond lengths and these are listed as Bi-S1 (3.00 $\AA$) and Bi-S2 (3.51 $\AA$) in Table \ref{table2} and Figure~\ref{fig:str}(a). While in the $I4/mmm$ phase, Bi and S bonds have four-fold coordination with only one in-equivalent bond with a length of 3.14 $\AA$. Similarly, the bond length between Bi and O atoms in the $Pnmn$ phase ranges from 2.24 to 2.43 $\AA$, whereas in the $I4/mmm$ phase, all Bi-O bonds are identical, with a bond length of 2.19 $\AA$. Overall, the bonding is highly heterogeneous in the $Pnmn$ phase over $I4/mmm$ phase.

A recent HP-XRD study\cite{Bu2020} revealed that Bi$_2$O$_2$S undergoes a structural phase transition from orthorhombic ($Pnmn$) to tetragonal ($I4/mmm$) structure at 6.4 GPa due to suppression of the $6s^2$ lone-pair and then it transforms from 2D to 3D structure upon further compression above 13.2 GPa. To get further insights on structural phase transition and anharmonicity, we have systematically investigated the high pressure behavior of Bi$_2$O$_2$S up to 10 GPa. As shown in Figure~\ref{fig:lattice}(a), the enthalpy difference between ambient and high pressure phases is gradually decreasing with pressure up to 3 GPa and then after both phases exhibit iso-enthalpic nature above 4 GPa. This result suggests a continuous phase transition from $Pnmn$ $\rightarrow$ $I4/mmm$ in Bi$_2$O$_2$S. To get further insights on the structural phase transition, we have analyzed the pressure dependent lattice constants, bond parameters as a function of pressure and equation of state (EOS). The calculated lattice constants are decreasing monotonically with pressure.  Interestingly, the lattice constants 'a' and 'b' are merging into a single lattice constant above 4 GPa and lattice constant 'c' continues to be 'c' for high pressure tetragonal ($I4/mmm$) phase (see Figure \ref{fig:lattice}(b)) and this result is consistent with the variation of lattice constants as a function of pressure as observed in the HP-XRD measurements.\cite{Bu2020} The calculated pressure coefficients by fitting the pressure-lattice parameter data to a quadratic expression show that the $Pnmn$ lattice is highly compressible along a-axis.  As shown in Figure \ref{fig:lattice}(c), the calculated volumes for both $Pnmn$ and $I4/mmm$ phases are converging and become almost equal above 3 GPa until the maximum studied pressure range (0-10 GPa) in this work and it is in excellent agreement with the HP-XRD measurements.\cite{Bu2020}  As illustrated in Figures \ref{fig:BL}(a-b), the calculated in-equivalent bond lengths (Bi-S $\&$ Bi-O) and angles (S-Bi-S $\&$ O-Bi-O) of lower symmetry ambient $Pnmn$ phase are converging and merging into a single in-equivalent bond length/angle for higher symmetry $I4/mmm$ phase above 4 GPa. The static enthalpy calculations and pressure evolution of structural properties of ambient and high pressure phases clearly demonstrate that Bi$_2$O$_2$S undergoes a continuous phase transition from $Pnmn$ $\rightarrow$ $I4/mmm$ under high pressure.

\subsection{Cross-band-gap hybridization and enhanced Born effective charges}
It is well-known fact that the standard DFT functionals underestimate the band gaps\cite{Heyd2005,Verma2017} by 30-40 $\%$ for semiconductors and insulators due to self-interaction, to overcome this problem, we have used TB-mBJ\cite{Tran2009} potential to calculate electronic structure of Bi$_2$O$_2$S. The calculated band gaps are 0.82 and 0.67 eV without and with inclusion of spin-orbit coupling (SOC) within PBE-GGA, respectively for $Pnmn$ phase and the corresponding band gaps using TB-mBJ potential are 1.52 and 1.33 eV. The obtained PBE-GGA band gap 0.82 eV without SOC is consistent with previously obtained band gap of 1.0 eV,\cite{Zhang2015} and the TB-mBJ band gap with SOC 1.33 eV is consistent with HSE06 band gap of 1.25 eV.\cite{Wu2017} As shown by Shi et al.~\cite{Shi2016}, materials with fully occupied lone-pair ($ns^2$, n = 4, 5, 6) show distinct features in their electronic structure by having an additional band (dominated by fully occupied cation $s$-states) below the valence band compared to the materials without lone-pair cation. The conduction band is mainly dominated by extended $p$-states of lone-pair cation in contrast to the cation $s$-states in materials without lone-pair cation. The calculated projected electronic density of states (PDOS) for both ambient orthorhombic and high pressure tetragonal (at 5 GPa) phases of Bi$_2$O$_2$S are presented in Figures \ref{fig:DOS}(a) and (b), respectively. As discussed above, the bottom of the valence band is mainly derived from $s$-states of Bi$^{3+}$ cation for both of these phases ($Pnmn$ and $I4/mmm$). As shown in Figure \ref{fig:DOS}, the top of the valence band is mainly dominated by the $p$-states of S$^{2-}$ and O$^{2-}$ anions and also minor contribution from $s$ $\&$ $p$-states of the Bi$^{3+}$ cation, while the conduction band is mainly dominated by extended $p$-states of the Bi$^{3+}$ cation, which hybridize with anion (S$^{2-}$ and O$^{2-}$) $p$-states of the valence band, thus, results in a significant cross-band-gap hybridization in Bi$_2$O$_2$S. The highly dispersive nature of $p$-states of S$^{2-}$ and O$^{2-}$ anions and Bi$^{3+}$ cation at the top of the valence band and minimum of the conduction band, respectively provides high mobility for light-induced e-h separation and their transport in photocatalysis.\cite{Zhang2015} Moreover, the electronic structure also gives a fundamental insight on stereochemical activity of a lone-pair.\cite{Waghmare2003} The cation $s$-states are broader (see Figure \ref{fig:DOS}(a)) and are inter-mediated by unoccupied cation $p$-states with anion $p$-states. When the activity of lone-pair is diminishing, the cation $s$-states (below the valence band) are narrow and have sharp peaks in the PDOS (see Figure \ref{fig:DOS}(b)). The observed cross-band-gap hybridization leads to a mixed ionic-covalent character, which causes enhanced lattice polarization in lone-pair cation containing compounds.\cite{Du2010,Du2010InX} This could be clearly seen from Born effective charges (BECs), as presented in Table \ref{Born}, the calculated BECs are significantly enhanced for both phases of Bi$_2$O$_2$S compared to the formal charges of cation (+3) and anions (-2). Enhanced BECs bring the lattice in the proximity of ferroelectric instability\cite{Du2010} thereby lowering the lattice thermal conductivity ($\kappa_l$)~\cite{pandey2020lattice}. Therefore, it is very intriguing to investigate elastic and phonon transport properties of lone-pair containing materials in general and in particular for Bi$_2$O$_2$S, especially to understand what is the effect of lone-pair on lattice dynamics and phonon transport at high pressure?

\subsection{Elastic properties, lattice stability and phonon transport}
To explore the mechanical and dynamical stability, we calculated second order elastic constants (SOEC) and phonon dispersion curves, respectively for $Pnmn$ and $I4/mmm$ phases of Bi$_2$O$_2$S at ambient as well as at high pressure. The calculated SOEC at 0 and 5 GPa for both the phases are given in Table \ref{elastic}. At 0 GPa, the $Pnmn$ phase obeys necessary and sufficient conditions of Born stability criteria\cite{Mouhat2014} and no imaginary frequencies are observed, which indicate that the $Pnmn$ phase is mechanically and dynamically stable at ambient pressure (see Figure \ref{fig:PHDISP}(a)). The $I4/mmm$ phase doesn't satisfy the Born stability criteria (C$_{44}$ $\textgreater$ 0) and also it possesses imaginary frequencies along $\Gamma$-M direction (see Figure \ref{fig:PHDISP}(b)) at 0 GPa, which clearly demonstrate that $I4/mmm$ phase is mechanically and dynamically unstable at ambient pressure for Bi$_2$O$_2$S.

When a non-zero hydrostatic pressure is applied to a particular crystal system, the Born stability criteria must be revised with modified elastic constants of a crystal system at a given hydrostatic pressure (P) conditions. According to Sinko and Smirnov\cite{Sinko2002,Sinko2004} the modified elastic constants at a given pressure are $\widetilde{C}_{ii}$ = $C_{ii}$ - P, (for $i$ = 1-6) and $\widetilde{C}_{12}$ = $C_{12}$ + P,  $\widetilde{C}_{13}$ = $C_{13}$ + P and $\widetilde{C}_{23}$ = $C_{23}$ + P. Therefore, the modified Born stability criteria under hydrostatic pressure for orthorhombic crystal symmetry can be found in our previous work,\cite{Yedu2015} whereas for tetragonal system the same are given as follows:
\begin{equation}
\begin{split}
C_{11} - C_{12} - 2P  >  0, \\ P^2 + P (C_{11} + C_{12} + 4C_{13}) + 2C^2_{13} - C_{33}(C_{11} + C_{12}) > 0, \\  C_{44} - P  >  0, \hspace{0.2in} C_{66} - P > 0
\end{split}
\end{equation}
Under high pressure ($\textgreater$ 4 GPa), for instance at 5 GPa (see Table \ref{elastic}), the calculated SOEC for both $Pnmn$ and $I4/mmm$ phases have the following relationships: C$^{Pnmn}_{11}$ $\approx$ C$^{Pnmn}_{22}$ $\approx$ C$^{I4/mmm}_{11}$, C$^{Pnmn}_{44}$ = C$^{Pnmn}_{55}$ = C$^{I4/mmm}_{44}$ and C$^{Pnmn}_{12}$ $\approx$ C$^{Pnmn}_{13}$ = C$^{I4/mmm}_{12}$. These relationships strongly suggest a continuous phase transition from $Pnmn$ to $I4/mmm$ is occurring in Bi$_2$O$_2$S under high pressure and is in good accord with the pressure dependent structural properties such as pressure dependent lattice constants and bond parameters (see Figures \ref{fig:lattice} and \ref{fig:BL}). In addition, no imaginary frequencies are found for the computed phonon dispersion curves at 5 GPa, which strongly suggest that $I4/mmm$ phase becomes dynamically stable at 5 GPa.

As discussed above, under the application of pressure the lone-pair of Bi$^{3+}$ cation is suppressed, which can have important consequences on phonon transport. To investigate this further we analyze the lattice dynamics and phonon transport in both $Pnmn$ and $I4/mmm$ phases of Bi$_2$O$_2$S. The calculated phonon dispersion and phonon density of states (PHDOS) are shown in Figures~\ref{fig:PHDISP}(a) and (b) for $Pnmn$ and $I4/mmm$ phases, respectively. In general, for both $Pnmn$ and $I4/mmm$ phases, the three acoustic modes are comprised of very low frequencies ($\leq$ 1.5 THz) and they overlap with optical branches. In both phases, the optical phonons show very low frequency modes due to presence of heavy metal Bi and they hybridize with the acoustic phonon modes. The optical phonon branches in high-frequency regions ($>$ 4THz) are very dispersive and exhibit large phonon group velocity, indicating the possibility of contribution to $\kappa_l$. The PHDOS in these compounds can be divided into three main regions: Bi dominated region below 1 THz, mixed Bi and S region above 1 THz, and O dominated region above 6 THz. Despite the huge mass contrast between Bi and S (m$_{Bi}$/m$_S$ = 6.5), there is a significant overlap between Bi and S in the PHDOS. This overlap is larger for the $Pnmn$ phase (1 THz - 5.5 THz) than in the $I4/mmm$ phase (3 THz - 5.5 THz). The sulfur PHDOS in the $Pnmn$ phase shows several small peaks in 1 THz - 5.5 THz and a large peak around 5 THz. The $I4/mmm$ phase also shows sulfur peaks in the PHDOS, but these are less pronounced, reflecting more dispersive optical branches. The large overlap between acoustic and optical phonons can lead to enhanced scattering of the heat-carrying acoustic phonons --- suggesting the possibility of relatively lower $\kappa_l$ for the $Pnmn$ phase.

The temperature-dependent $\kappa_l$ of the two phases are compared in Figure~\ref{fig:Kappa}(a) along in-plane and cross-plane directions. As seen, $\kappa_l$ of the ambient $Pnmn$ phase is much lower than the high pressure $I4/mmm$ phase. At room temperature, the average $\kappa_l$ values are 1.47 and 3.63 W/m-K for $Pnmn$ and $I4/mmm$ phase, respectively. Since both $Pnmn$ and $I4/mmm$ phases have an anisotropic crystal structure, the same is reflected in the calculated $\kappa_l$. Due to the orthorhombic structure, the $Pnmn$ phase exhibits different $\kappa_l$ along in-plane x and y directions and are become comparable at higher temperatures. As can be seen from Figure~\ref{fig:Kappa}(a) both phases show anisotropy in $\kappa_l$ with $\kappa^{\parallel}_l$/$\kappa^{\perp}_l$ ratio as 2.14 and 2.58 in $Pnmn$ and $I4/mmm$ phase, respectively. For instance, the calculated thermal anisotropy 2.14 of ambient phase of Bi$_2$O$_2$S is comparable with thermal anisotropy 2.0 for Bi$_2$O$_2$Se in $I4/mmm$ symmetry.\cite{Guo2021} As expected for layered materials with weak electrostatic interactions,\cite{Guo2021} the $\kappa_l$ along the cross-plane direction is smaller than the in-plane direction. Interestingly the calculated thermal transport anisotropy is much smaller than what is typically observed in vdW layered materials such as, MoS$_2$ where $\kappa^{\parallel}_l$/$\kappa^{\perp}_l$  of nearly 30 was found~\cite{meng2019thermal}. This clearly distinguishes between the role of interlayer weak electrostatic and vdW bonding in determining the thermal anisotrpy in layered materials.

To understand the thermal transport behavior further we calculate the cumulative $\kappa_l$ ($\kappa_{ac}$) as a function of phonon frequency, which is shown in Figure~\ref{fig:Kappa}(b) and (c) for $Pnmn$ and $I4/mmm$ phases.  As can be seen there, in both the phases, the optical phonons contribute significantly to $\kappa_l$ particularly along the in-plane direction. The phonon frequency range contributing to cross-plane and in-plane $\kappa_l$ is also quite different. The cross-plane $\kappa_l$ is mostly contributed phonons with frequency $<$ 5 THz, whereas phonons with frequencies up to 12 THz contribute to in-plane $\kappa_l$ in both the phases. This can be understood by analyzing phonon group velocities which are shown in the Figure~\ref{fig:Kappa}(d) and (e) for $Pnmn$ and $I4/mmm$ phases, respectively. In both the phases for the majority of the frequency range $v_g$ is higher along in-plane direction than the cross-plane direction. This is more prominent for optical phonons in the frequency window of 6-12 THz, where in-plane $v_g$ is more than two times larger than the cross-plane $v_g$. These higher in-plane $v_g$ are the main reason for higher $\kappa_l$ along the in-plane direction. For the $Pnmn$ phase, the $v_g$ along in-plane x and y directions are comparable which is consistent with the comparable $\kappa_l$ along these directions.

We then focus on exploring the origin of the relatively low $\kappa_l$ in the $Pnmn$ phase. To explain the calculated differences in $\kappa_l$, we look into three factors, phonon group velocity, specific heat, and three phonon scattering rates, that determine $\kappa_l$. As explained above the phonon group velocities are comparable in both phases --- therefore phonon velocities are not the reason for the lower $\kappa_l$ in the $Pnmn$ phase. As for the second-factor specific heat, the $Pnmn$ phase has only slightly higher ($\sim$ 1.3 \%) specific heat than the $I4/mmm$ phase. This indicates that the difference in phonon anharmonicity should be the main reason for different $\kappa_l$. The phonon anharmonicity governs  the strength of phonon-phonon interactions which determines the overall intrinsic thermal resistance. Greater anharmonicity results in shorter phonon lifetimes and lower $\kappa_l$. Indeed, as shown in Figure~\ref{fig:Kappa}(f), the three phonon scattering rates of the $Pnmn$ phase are significantly larger than the $I4/mmm$ phase particularly in the frequency range 6-15 THz. The optical phonons in $Pnmn$ phase are more dispersive than in $I4/mmm$ phase, which enables it to satisfy the conservation laws required for three phonon scattering resulting in higher scattering phase space and therefore higher scattering rates.

As introduced before the main difference between orthorhombic $Pnmn$ and tetragonal $I4/mmm$ phases, is the absence of Bi 6$s^2$ lone-pair electrons in the $I4/mmm$ phase, which may explain its relatively less anharmonicity.  To gain further insights into the role of lone-pairs electrons in bonding and anharmonicity, we next calculate the potential energy by off-centering the Bi, O, and S atoms along Cartesian x and z directions. The resulting potential energy plots are presented in Figures~\ref{fig:PES}(a) and (b) for $Pnmn$ and $I4/mmm$ phases, respectively. As can be seen there, the displacement potential energy of Bi and S atoms are significantly different along x and z directions in both $Pnmn$ and $I4/mmm$ phases. We find that for Bi and S atomic displacements in the x or y directions, the potential energy curve is relatively flat, whereas displacement along z-direction requires more energy. In both phases the energy required for O displacement is comparable. The potential energy analysis stipulates that Bi and S vibrational modes could be easily excited by low energy. This is consistent with the calculated PHDOS, where despite the huge mass difference, PHDOS of Bi and S atoms show significant overlap.  Furthermore, for Bi vibration in $Pnmn$ phase potential energy curve shows asymmetric features and deviates from parabolic (harmonic) behavior. The polynomial fit of this potential energy well gives anharmonic coefficient as 2.99 and 1.68 along x and z directions, respectively. Such large deviation from typical parabolic behavior corroborates the strong anharmonicity of this phase. For comparison,  in $I4/mmm$ phase the anharmonic coefficients for Bi atomic vibrations along z direction is 1.56, and for Bi displacement along x direction the potential is harmonic (coefficient of anharmonic fit is zero). The above displacement potential energy analysis confirms that the $Pnmn$ phase is more anharmonic than $I4/mmm$ phase.  This is further confirmed by plotting the mean square displacements (MSD) for Bi and S atoms which are shown in Figure~\ref{fig:PES}(c) for both $Pnmn$ and $I4/mmm$ phases. The MSD for Bi and S atoms are higher in $Pnmn$ phase  (Bi = 0.012 $\AA^2$, and S= 0.024 $\AA^2$ at 300 K), than in the $I4/mmm$ phase (Bi = 0.009 $\AA^2$, and S= 0.011 $\AA^2$ at 300 K). This is expected because the longer Bi-S2 bond lengths in $Pnmn$ phase make it easier to move the atoms from their equilibrium positions during the vibration than in the $I4/mmm$ phase. Overall, the anharmonicity of the $6s^2$ lone-pair of Bi$^{3+}$ cation plays a pivotal role in determining the phonon transport in Bi$_2$O$_2$S and it is suppressed (in $I4/mmm$ phase) upon compression as can be seen from chemical pre-compression of Se/Te atoms in Bi$_2$O$_2$Se and Bi$_2$O$_2$Te compounds.

\section{Conclusions}
In summary, we have systematically investigated the effect of hydrostatic pressure on stereo-chemically active 6$s^2$ lone-pair containing Bi$^{3+}$ cation with weak electrostatic interlayer bonding compound, Bi$_2$O$_2$S and their implications on crystal structure, lattice dynamics and phonon transport properties using first principles calculations and Boltzmann transport theory. We predict that Bi$_2$O$_2$S undergoes a continuous structural phase transition from low symmetry distorted ($Pnmn$) structure to a high symmetry ordered ($I4/mmm$) structure at around 4 GPa, which is comparable to the transition pressure 6.4 GPa determined from HP-XRD measurements\cite{Bu2020}. The mechanical and dynamical stability of both $Pnmn$ and $I4/mmm$ phases are confirmed from the calculated elastic constants and phonon dispersion curves, respectively. The obtained enhanced Born effective charges (BECs) clearly demonstrate that both the phases show significant cross-band-gap hybridization from their electronic structure suggesting near ferroelectric instability and is favorable to achieve low $\kappa_l$. The obtained low $\kappa_l$ value 1.47 W-m/K for $Pnmn$ phase over 3.63 W-m/K for $I4/mmm$ phase is due to suppression of lone-pair of Bi$^{3+}$ cation at high pressure which reduces anharmonicity in $I4/mmm$ phase. The highly anharmonic behavior of Bi and S atoms of $Pnmn$ phase demonstrated from PES and MSD at ambient pressure is mainly responsible for low $\kappa_l$. Overall, the present work provides an insight on how hydrostatic pressure effects the phonon transport in  stereo-chemically active 6$s^2$ lone-pair containing Bi$_2$O$_2$S. The present study could stimulate exploring pressure dependent phonon transport properties in materials especially with lone-pair cation, which will provide a new avenue of pressure modulated low $\kappa_l$ materials for future thermal energy applications.

\section{Acknowledgments}
NYK and TP contributed equally to this manuscript. NYK would like to thank Institute for Advanced Computational Science, Stony Brook University for providing computational resources (Seawulf cluster). TP is supported by Research Foundation Flanders (FWO-Vl). SCRR would like thank RGUKT Basar for providing computational facilities.