Multiantenna Multisource Multirelay Spectrum Sharing Networks: Buffer-Aware Direct Link Activation for Higher Reliability and Lower Overhead

Cooperative relaying can be implemented in spectrum sharing networks to extend the range of reliable communication. In this paper, we incorporate multiple-antenna technology and buffer-aided relaying to guarantee a highly reliable connectivity for the secondary network with several source nodes. The multi-antenna configuration for the sources and destination nodes suggests that there are several potential source-to- destination channels which can be auspicious for the network in two ways. First: to balance up the buffer states. Second: to expand the link selection opportunity at each time step. Motivated by this rationale, we propose a buffer-aware communication protocol to incorporate the direct transmissions along the relaying links without incurring excess overhead for circulation of channel-state-information (CSI). Considering Nakagami-m fading, we derive closed-form expressions for the end-to-end (ete) outage probability and average packet delay of the secondary network under the proposed protocol. Through Monte-Carlo simulations, we investigate the influential network parameters and evaluate the proposed technique in comparison to two benchmark schemes, one with buffer-based link-prioritization and one without prioritization. Findings demonstrate that the proposed protocol outperforms both benchmarks in terms of outage probability and ete delay, especially as the number of nodes scales up. However, the superior performance comes at the cost of more CSI circulation. Furthermore, it is shown that if the global CSI cannot be collected accurately, then the dependency of the link selection on the global CSI should be relaxed to mitigate the performance loss. The theoretical and Monte-Carlo results coincide in several simulation examples, verifying the presented theoretical analysis.


I. INTRODUCTION
Over the course of last decades, the steadily growing demand for fast and ubiquitous connectivity has caused an astronomical utilization of electromagnetic spectrum. As a limited resource, the spectrum is prone to scarcity and in turn, it is crucial to enhance the efficiency of spectrum utilization [1]. In this context, enabled by cognitive radio technology, the The associate editor coordinating the review of this manuscript and approving it for publication was Chuan Heng Foh . spectrum sharing networks can effectively mitigate the spectrum congestion through a careful secondary access to the licensed frequency bands [2]. However, the range of reliable communication in secondary networks is constrained because the transmit power is stringently regulated to protect the primary network from interference. To resolve this issue and extend the communication to a far-located secondary user, an intermediate relay node can be employed [3], [4], [5]. Spectrum sharing relay networks are considered one of the most promising technologies to cope with the VOLUME 10, 2022 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ spectrum scarcity problem. However, providing a highly reliable connectivity for a large secondary network is a challenging problem, especially for the higher transmission rates. In this context, enabling technologies such as buffer-aided relaying [6], [7], [8], [9], [10], [11], [12], [13], [14] and multiple-antenna technology [15], [16], [17], [18], [19] can be incorporated to substantially improve the reliability within the secondary network.

A. RELATED LITERATURE
It has been shown that data-buffering at the relaying terminal would bring about a marked improvement in terms of reliability and throughput [20], [21], [22]. Such a superior performance over the conventional buffer-less schemes stems form the flexibility that data buffering brings to the link activation: any relaying link can be activated at any time for reception from the source or transmission towards destination as long as the buffer is not completely full or empty [23]. 1 The advantage of buffer-aided relaying becomes more pronounced in a multi-relay scenario because the effect of flexibility in link selection scales up as the number of relay nodes grows [24], [25], [26]. Evidently, the link selection plays a key role in any communication system with data-queuing [27], [28], [29], [30], [31] including buffer-aided relaying networks. One of the earliest selection techniques for the multi-relay configuration was coined max-max relay selection (MMRS); in this technique the odd time slots are allotted to the links between the source and relays and even time slots to ones between relays and the destination [24]. The max link relay selection (MLRS) was introduced shortly later; in this technique, each time slot is opportunistically assigned to strongest link among the available relaying links [25]. In comparison, the MLRS offers a lower outage probability while the MMRS can provide a lower transmission delay because it operates based on a prefixed temporal schedule. However, since a full (empty) buffer would not be available for reception from the source (transmission towards the destination), in the follow up contributions, researchers proposed a number of buffer-aware link selection protocols to mitigate the performance loss caused by such phenomena [32], [33], [34], [35], [36], [37], [38], [39], [40], [41], [42].
Motivated by the advantages of data-buffering, several research works incorporated the buffer-aided relaying techniques into the conventional spectrum sharing networks [6], [7], [8], [9], [10], [11], [12], [13], [14]. In detail, Chen et al. investigated the best relay selection policy with several halfduplex decode-and-forward (DF) relays and derived the outage probability [6]. A collaborative communication scenario between the primary and secondary users was explored in [7] where the secondary source is equipped with a buffer and occasionally assists the primary network as relaying terminal. A rate-optimal adaptive link selection technique similar to [21] with a single DF relay was studied in [9] and [10] where transmit power of the secondary nodes are assumed prefixed in the former while in the latter, the transmit power is a random parameter governed by the interference channel. Zhang et. al explored a generalized buffer-aware relay selection scheme [12]. The effect of direct link transmissions was explored in [13] in terms of throughput maximization and in [14] in terms of outage minimization. Recently, this concept was applied to a mixed RF/FSO spectrum sharing model [11].
On the other hand, the multi-antenna (MA) technology is yet another promising solution that can be employed to improve the reliability without incurring further spectrum consumption [43]. The concept of multi-antenna spectrum sharing relay networks has been extensively investigated in the past years [15], [16], [17], [18], [19]. Particularly, Mana et. al considered a multiple-inputmultiple-output (MIMO) secondary base station assisting the primary network as an amplify-and-forward (AF) relay while transmitting its own message to the secondary destination [15]. Deng et. al investigated the outage probability and bit error rate for a cognitive network composed of a MIMO source, a half-duplex DF relay, and a destination [16]. Yeoh et. al presented a similar contribution, but in the presence of interference from the primary transmitter [18]. To develop upon previous works, Elsaadany and Hamouda considered a multi-relay transmission scenario [19]; however, in their model, only secondary destination is a MA terminal.
Moreover, for a wide range of applications, the wireless networks are deployed to provide connectivity for several users [44], [45], [46], [47]. Thus, it is of practical importance to design spectrum sharing networks based on the multi-user requirements [48], [49], [50]. Recently, some researchers studied the network scheduling problem for the multi-user buffer-aided relaying networks [51], [52], [53], [54], [55], [56], [57], [58]. Particularly, an optimization problem with fixed transmission power was solved in [51] to maximize the average throughput for a unicast scenario, and in [52] for multicast scenario. A maximization problem similar to [21] was developed in [53] with a single half-duplex relay and multiple source nodes while Li et al. addressed the average throughput maximization with a full duplex relay [59]. Zhang et al. combined the simultaneous transmission of several source nodes and maximum likelihood detection at the relay and found the optimal power allocation to minimize the outage probability [55]. Jamali et al. found the optimal link activation policy as well as the optimal transmit power allocation to maximize the average throughput of the bidirectional relaying system [56]. An optimization problem similar to [56] was solved later in [57] for a bidirectional transmission model with multiple relay nodes. Nomikos et al. proposed NOMA transmission for the uplink communication where the interference cancellation at the relay nodes is performed to decode the simultaneously transmitted signals of the users [58].

B. MOTIVATION AND CONTRIBUTIONS
Interestingly, despite the great potentials of an all-multiantenna configuration for the spectrum sharing networks, the resource allocation problem for such networks has not been investigated. Moreover, it is envisaged that the jointly incorporation of multiple-antenna technology and bufferaided relaying can markedly enhance the reliability of the secondary networks. In this context, the multi-antenna deployment and multi-source configuration of the network implies that there are several potential source-to-destination channels which can be utilized in two ways. First: to balance out the buffer states. Second: to expand the link selection opportunity at each time step. However, incorporating the direct transmissions would require further CSI circulation in the network. Thus, a communication protocol need to be developed to effectively leverage the available resources towards achieving a highly reliable secondary connectivity. Motivated by these points, we aim to develop a buffer-aware communication protocol for the all-MA secondary network which smartly incorporates the direct transmissions along the relaying links to harness the benefits of potential direct channels while mitigating the overhead required for CSI circulation.
Although the MA buffer-aided relaying model has been recently investigated for the conventional networks in few works [60], [61], [62], [63], [64], to the best of our knowledge, applying such a combination to the spectrum sharing networks has not been explored yet. Furthermore, as highlighted in Table 1, there are other gaps in the related literature in terms of antenna selection, fading analysis, network configuration, and direct transmission opportunity. This work also bridges the aforementioned gaps; our major novelties and contributions can be listed as follows: • This work leverages the combination of multi-antenna technology and buffer-aided relaying to boost the reliability of the spectrum sharing networks considering multi-source multi-relay topology and multi-antenna configuration at each node.
• A novel communication protocol is introduced under which the secondary transmissions are prioritized based on the buffer status of the relays. In contrast to the state of the art [60], [61], [62], [63], [64], our proposed protocol dynamically incorporates direct communications along the relaying transmissions to balance out the buffer dynamics and to increase the number of employable communication links. Furthermore, to mitigate the additional overhead caused due to CSI circulation of the direct channels, we propose a novel buffer-aware CSI collection technique to be used at the network coordinator node.
• Since the system model as well as the proposed scheduling technique are new and subject to spectrum sharing constraint, a new theoretical framework is developed to provide a fast and accurate tool for the performance study. In this vein, a block fading Nakagami-m fading environment is considered for all the network channels and the antenna selection technique is incorporated to reduce the implementation complexity. It is worth mentioning that the presented works in [61], [62], and [63] lack fading analysis and [60] considered a less general Rayleigh model. In the following sections, we describe the system model, propose scheduling protocols, develop a mathematical framework to evaluate reliability and delay performance of the proposed protocol, we create two benchmarks to present a fair performance comparison, and finally we evaluate our proposed secondary network through computer simulations.

II. SYSTEM MODEL
As illustrated in Fig. 1, the spectrum sharing relay network that we consider includes N secondary sources (SS), a secondary destination (SD), K decode-and-forward (DF) secondary relays (SR), and a primary network in the vicinity. A centralized network structure is considered and the global channel state information (CSI) of the network is assumed to VOLUME 10, 2022 be available at the network's coordinator node. The SR nodes operate in the half-duplex mode and over equal-length time slots. Each relay is equipped with a data buffer that can store a maximum number of M packets. For the sake of tractability, the sources, the relays, and the corresponding number of packets stored in the buffers are respectively distinguished by indices as s n , r k , and Q k where n ∈ {1, . . . , N }, k ∈ {1, . . . , K }, and Q k = {0, . . . , M }, i.e., Q k shows the number of packets currently stored in the buffer of k-th relay. The packets are assumed to be numbered so that the destination can rectify the disarrangement of received packets due to dynamic link activation.
All the nodes in the network are equipped with multiple antennas where L s , L r , and L d respectively represent the number of antennas at the source, relay, and destination nodes. To facilitate implementation while harnessing the benefits of multi-antenna technology, transmit antenna selection (TAS) and maximum ratio combining (MRC) schemes are performed at the transmitter and receiver sides of the communications [65]. To distinguish the manifold wireless channels, we let s l n represent the l-th antenna of n-th secondary source, r j k the j-th antenna of k-th secondary relay, d j the j-th antenna of secondary destination, and p denote the primary network. Adopting these terminologies throughout the paper, we refer to the channel coefficient between the transmitter a and receiver b as h ab where a ∈ {s l n , r j k } and b ∈ {r j k , d j , p}. Also, the associated channel gains are denoted g ab , i.e., g ab = |h ab | 2 . For instance, g s 2 1 r 1 2 is the channel gain between the second antenna of the first source and the first antenna of the second relay. The SS as well as the SR are considered as two centralized clusters such that the clustered nodes are located roughly in the same distance to the other network nodes. Therefore, the channels from the nodes within a cluster to an another node in the network are considered independent and identically (i.i.d) distributed random variables (RV). All links in the network are assumed to undergo Nakagami-m block fading and therefore, the channel coefficients are constant within each block of data while changing independently from one block to the other; the probability density function (PDF) as well as the cumulative distribution function (CDF) for the typical channel gain g ab are given by [66] where, ab is the average power level, m ab is the fading severity, (.) is the Gamma function, and (., .) is the upper incomplete Gamma function [67], [68], [69]. Due to the spectrum sharing regulation, the transmit power at the secondary nodes must be controlled such that the interference seen at the primary network remains below the interference temperature, denoted I th . Also, a maximum transmittable power limit, denoted P M , is considered and signals at the receivers are perturbed by the AWGN that has a variance of σ 2 .
Throughout this manuscript, we refer to the secondary links between the sources and relays as S-R links; similarly, the terms S-D and R-D are used to refer to the direct links between the sources and destination, and the relaying links between the relays and destination, respectively.

III. PROPOSED COMMUNICATION PROTOCOL
Since the number available links for communication directly affects the reliability of network, it is critically important to develop a buffer-aware communication protocol to avoid performance loss caused by the unavailability of buffers [23]. To this end, since the CSI circulation and link selection protocol are intertwined, we can jointly design both mechanisms to mitigate the overhead of CSI circulation. In the following, we first clarify our proposed link selection regardless of the overhead required for CSI circulation and then, we describe our buffer-aware CSI circulation mechanism to mitigate the overhead.

A. PROPOSED LINK SELECTION TECHNIQUE
It is clear that maintaining buffers as close as possible to half-full can effectively preserve the long term availability of the relays for transmission and/or reception [37]. Assuming that the information of buffers is collected, the coordinator updates two decision variables, namely N lg and N ld , where N lg is the number of relays whose buffer is less than half-full and N ld the ones whose buffer is more than half-full. Mathematically, these two parameters can be described as where F is a function which takes in a set and returns the number of non-empty elements, and . is the floor operation. 2 Next, the buffer-status of the relaying network B is determined based on N lg and N ld and Q k as bellow Once information regarding the channel and buffer states are collected, the S-R, R-D, and S-D transmissions are prioritized depending on B: 1) if B = ce, then only S-R transmission is considered upon which the strongest S-R link is selected and then communication starts if the selected link is able to support the network's target rate. Otherwise, there will be no transmission.
2) if B = cf, then only R-D transmission is considered and subsequently, the strongest R-D link will be activated for transmission if it can support the network's target rate. Otherwise, there will be no transmission. 3) if B = lg, then the transmission window is prioritized as follows: first S-R channels are considered and the best S-R link is selected to start communication. If the selected link is not realizable for the network's target rate, the time slot is given to S-D and the best S-D link is selected for transmission. However, if the selected S-D link is also unrealizable, then the R-D channels will be considered and the best link is selected to start communication. If all the links fail, there will be no transmission. 4) if B = ld, then similar to the lg scenario, the S-R, S-D, and R-D are prioritized but with the following pattern: first priority is given to R-D, second to S-D, and third to the S-R channels. 5) if B = ow, then similar to the lg and ld scenarios, the S-R, S-D, and R-D are prioritized but with the following pattern: first priority is given to R-D, second to S-R, and third to the S-D channels. In the following, we propose a novel buffer-aware CSI collection technique to mitigate the overhead required for link selection.

B. PROPOSED CSI COLLECTION TECHNIQUE
The CSI of the channels are estimated at the relays and destination through listening to the pilot signals that source and relays broadcast periodically. Prior to each transmission slot, a sequential three-stage CSI collection algorithm is executed at the network coordinator unit. In the first stage, the information of buffer states is collected; in the second stage, the coordinator collects the CSI of the direct links from the destination according to the buffer states; in the third stage, the CSI of the relaying links will be collected based on the information gathered in the first and second stages. The three stages are clearly expounded in the following: In the first stage, the buffer-states are collected and B is determined as (3). In the second stage, it is decided whether or not collect the CSI of the direct channels based on B: if B ∈ {ce, cf}, then the CSI will not be collected; However, if B ∈ {lg, ld, ow}, then the CSI will be collected. In the third stage, the information gathered in the first and second stages are utilized to collect the information of the relaying links: The reason behind this approach is that in the {ce, cf} scenarios only S-R and R-D channels can be activated and therefore, there is no need to collect the CSI of other links. On the other hand, in the lg and ld cases, we aim to activate the S-D link instead of the one relaying link whose activation would aggravate the off-balance buffers. As the result if, the S-D is realizable for the target rate, then there is no need to collect the CSI of the S-R and R-D links in the ld and lg cases, respectively. Lastly, in the ow case, since the S-D link would be considered for activation when both of the S-R and R-D links are not realizable, its CSI should be collected. Table 2 summarizes several possible scenarios for the network transmissions and buffer dynamics considering K = 4 and M = 3. The third column shows the availability of the relays where N-A means ''not-applicable''. For instance, the queues at the relays are given by [0, 2, 1, 2] for the state S c ; since the half-full threshold is 1.5 = 1, there is one buffer less than half-full and two buffers more than halffull, yielding N lg = 1, N ld = 2, and ld status as N ld > N lg . Similarly, S 1 represents the ce status; S a as well as S d represent the lg status; S c , S e , S g , and S h belong to the ld status, S b and S f to the ow, and S 256 to the cf status. Note that for K = 4 and M = 3 we have Q k ∈ {0, 1, 2, 3} where k ∈ {1, 2, 3, 4}. Therefore, the last state is numbered 4 4 = 256 since the total number of unique states can be found as (M + 1) K where the additional one represents the completely-empty state of the buffers.

C. AN ILLUSTRATIVE EXAMPLE
It is worth mentioning that the last column of the table highlights the transmission priority given the buffer status where xµy means that the transmission priority is given to y if x is not realizable. As can be seen, when buffers tend to be fuller, e.g., S e through S h , the R-D and S-D transmissions are prioritized over S-R links. Similarly, when the buffers tend to be emptier, S-R and S-D transmissions are prioritized over the R-D links.

D. PRELIMINARIES
Due to data-buffering at relays, we develop a queuing model for our network as depicted in Fig. 2(a). For better tractability, we simplify our queuing model as the equivalent model in Fig. 2(b) with a source, a queuing relay, and a destination where the equivalent transmission rates on the S-D, S-R, and VOLUME 10, 2022 R-D links denoted ρ sd , ρ sr , and ρ rd , respectively. We let S i denote the i-th unique combination of buffer states as Labeling the buffer states based on the queue-length, the queue status of the equivalent model is characterized as Q e ∈ {S 1 , . . . , S (M +1) K }, i.e., S 1 and S (M +1) K imply the ce and cf cases, respectively. Furthermore, we let sr i and rd i respectively represent the number of available relays for S-R and R-D transmissions given the buffer status S i ; mathematically speaking, where Q k (S i ) signifies the queue length at the k-th relay given S i . An example: we get sr 1 = 3, sr 256 = 0 in Table 2. Finally, as mentioned earlier, the secondary nodes perform MRC and TAS to reduce the implementation costs while harnessing the benefits of multi-antenna technology. In this context, letting γ ab denote the SNR of the link between transmitter a and receiver b, the SNR of the strongest links for the R-D, S-R, and S-D channels can be organized as where γ r l k d j is the SNR of the channel between l-th antenna of k-th relay and j-th antenna of destination; similarly, γ s l n r j k is the SNR of the channel between l-th antenna of n-th source and j-th antenna of k-th relay; also, γ s l n d j is the SNR of the channel between l-th antenna of n-th source and j-th antenna of destination. Considering the underlay spectrum sharing constraint, it is straightforward to show that

IV. OUTAGE ANALYSIS
The outage probability is an important metric for qualifying the reliability of a communication network. A communication outage occurs when the SNR of the channel is less than a minimum required value depending on the target bit rate. If we let W denote the transmission bandwidth, the required SNR threshold associated with a fixed target rate ρ is given by = 2 ρ W − 1. As illustrated in Fig. 3, the outage occurs if none of the available links are realizable in a given buffer status. Thus, the ete SNR of the network can be expressed as Adopting the Markov chain analysis framework, the secondary network's outage probability can be obtained as [14] where π i is the steady state probability for the i-th buffer state S i , and i is the associated conditional outage probability.
Here, from (5) and (6), it is inferred that γ s l n d j and γ s l n r j k are not independent due to the common RV g s l n p and as the result, it is not straightforward to directly derive the closed-form expression of i . However, it can be shown that i can be indirectly derived as (See Appendix A) where, In the following, we shall derive i and π i in closed-form.
A. i DERIVATION As (10) suggests, in order to derive i , it is required to find the closed-form expressions of A i , B i , and C i . In the following, we shall derive each term sequentially.

1) DERIVATION OF A i
By applying the concept of total probability to handle the cumbersome RV g s l n p and some modifications, A i can be rewritten as where f g s l n p (x) is the PDF of g s l n p . Substituting the SNR expressions from (6), (11) can be written as the sum of two terms as A i = (A 1 + A 2 ) L s N (See Appendix B): where α sd = m sd L d , α sr = L r m sr , β sd = m sd / sd , β sr = m sr / sr , β sp = m sp / sp , and γ (., .) represents the lower incomplete Gamma function. The first term can be readily obtained as However, solving the integral involved in A 2 is highly tedious, especially if m sd or m sr are not integers. Thus, it might be needed to solve it numerically. However, if m sd and m sr are integers, then a closed-form expression can be derived as (See Appendix C) where, V is a family set whose subsets include the unique summation indexes that add up to sr i , i.e., {v 1 , v 2 , . . . , v α sr +1 }| α sr +1 i=1 v i = sr i . Also, V means summation over all the subsets of V where for each subset we let

2) DERIVATION OF B i
This term is associated with R-D transmissions and since R-D channels are independent, it is straightforward to show that Similar to A i , placing the mathematical forms of SNRs into (15) and applying the concept total probability to unravel the cumbersome minimum function in γ r l k d j , i.e., min P M , I th /g r l k p , it is concluded that B i can be derived It is clear that B 2 , in general, has to be calculated using a numerical integration technique; however, a closed form Here, similar to A i and B i , after placing the mathematical expressions of the γ s l n r j k into (18) and applying the total probability concept, the closed form expressions of C i can be obtained as C i = (C 1 + C 2 ) NL s where C 1 and C 2 are given by (See Appendix D)    Table 2. Specifically, p +ce , p −cf , p +lg , p +ld , p −lg , p −ld , p +ow , and p −ow represent the transition probabilities associated with the ce, cf, lg, ld, and ow cases. Note: as it will be clarified in Section IV B, self-transitions can be found indirectly based on statistical properties of Markov transition matrix.
Similar to A 2 and B 2 , the integration involved in C 2 needs to be calculated numerically in general; however, in the case of integer m sr , a closed-form expression can be derived as where the definitions of V, C V , δ v , η v , and V are identical to ones given in (14).

B. π i DERIVATION
In order to obtain the steady state probabilities, it is required to find the transition matrix characterizing the Markov model. The transitions between two states are unique and predicate on the selection of a certain relay such that where P r (SR s ) = 1/ sr i and P r (SR s ) = 1/ rd i are respectively the probabilities that s-th secondary relay is selected for S-R and R-D transmissions. 3 Moreover, since there are three possible scenarios for transition between two states, we classify the transition probabilities in three major groups to facilitate the analysis. Fig. 4 illustrates an example for the possible transitions where S 1 , S a , . . . , S 256 are the same states described in Table 2.

1) FORWARD TRANSITIONS
this group represents the increase in the number of stored packets which occurs due to the successful S-R communication; Defining p +lg i , p +ld i , and p +ow i to distinguish forward transitions in lg, ld, and ow cases, it is straightforward to show that where sr i and sr i at the denominators account for the relay selection probability given in (22). In the following, we shall derive p +lg i ,p +ld i , and p +ow i subsequently. In the case of p +lg i , a careful observation on (10) reveals P r (γ SR > ) = 1 − C i . Therefore, the closed-form expression for p +lg i is given by In terms of p +ld i , because γ RD is uncorrelated with γ SR and γ SD , p +ld i can be reworded as 3 Note that the set of i.i.d S-R channels is not related to the set of i.i.d R-D channels and that the two sets do not necessarily have the same fading profile, i.e., m sr and sr can be different from m rd and rd .
Following the same approach we adopted in (11), a more tractable expression for (25) can be obtained as where B i is derived as explained in (16). After some algebraic modification, it is straightforward to show that where A i is derived as explained in (12) and D i can be obtained using the same derivation steps used for A i as D i = (D 1 + D 2 ) NL s where D 1 and D 2 are given by Similar to A 2 , B 2 , and C 2 , in general, the integral in (29) needs to be calculated numerically. However, for the integer m sd , a closed-form theoretical expression can be found as Regarding p +ow i , since γ RD and γ SR are uncorrelated, a closed-form expression can be readily derived as 2) REVERSE TRANSITIONS this group renders the decrease in the number of stored packets which occurs upon a successful R-D communication; we use a minus sign to distinguish this group. Similar to the forward transitions, defining p −lg i , p −ld i , and p −ow i to distinguish the lg, ld, and ow cases, the associated transition probabilities can be written as In the following we shall derive the closed form expressions of p −ld i , p −lg i , and p −ow i subsequently. From (32) and (10), it is clear P r (γ RD > ) = 1 − B i and therefore, On the other hand, since γ RD is independent from γ SR and γ SD , p −lg i can be derived as where the last step stems from the definition of A i in (10).

3) SELF-TRANSITIONS
this group encapsulates the instances in which the number of stored packets remains unchanged; however, there is no need to directly derive such probabilities thanks to the stochastic properties of the transition matrix, i.e., the sum of the elements on each column should be equal to one. Here, knowing the forward and reverse transition probabilities, we firstly form a transition matrix as T (i, j) = t i,j where i stands for the matrix row, j matrix column such that i = j, and for all the i, j ∈ 1, . . . , (M + 1) K , we have t ij = P r (S j → S i ). Next, we fill the diagonal probabilities such that the summation of each column is equal to one. Finally, the steady state probabilities can be obtained by solving the following equation: where, π = π 1 , . . . , π (M +1) K is the state probability column-wise vector; I M ,K is an (M +1) K by (M +1) K identity matrix; T is the modified transition matrix that incorporates the normalization equation i.e.,

V. DELAY ANALYSIS
The communication delay can be defined as the time a packet spends in the network from the instance of departure from the source until it is decoded at the destination. Evidently, the ete delay is a random variable and in turn, we concentrate on the average ete packet delay as the performance metric.
whereτ is the average ete packet delay,τ sd the average packet delay if direct transmission is chosen, andτ rd is the average delay if relaying transmission on the R-D link is chosen. Thus, for a large number of packets, it is implied that N sd N total = P sd , N sr N total = P sr , and N rd N total = P rd where P sd , P sr , and P rd are the probabilities of the successful transmission on the S-D, S-R, and R-D links given by where the arrow expression A→B signifies that the link between A and B is selected and the link is not in outage. It follows that the ete packet delay is given bȳ Finally, since direct transmissions is not a queuing system, we haveτ sd = 1 whereas for the relaying transmission, the associated delay can be found applying the Little's Law: Regarding the activation probability of the S-R and R-D links, a careful examination of (38) reveals that the P r (S→R|S i ) and P r (R→D|S i ) are linearly proportional to the transition probabilities derived earlier in Section IV. B. Particularly, for the forward transitions we have P r (S → R|S i ∈ lg) = sr i p +lg i , P r (S → R|S i ∈ ld) = sr i p +ld i , and P r (S → R|S i ∈ ow) = sr i p +ow i ; on the other hand, for the reverse transitions we have: . Furthermore, regarding P r (S→D|S i ), it is required to derive the activation probability which depends on the prioritization mechanism in lg, ld, and ow cases and can be summarized as Note that in (41), P sd = 0 in the ce and cf because the direct transmission is not considered. Here, a careful investigation of (41) reveals that the activation probabilities in lg and ld cases can be readily obtained using the expressions we derived earlier for A i , B i , C i , and D i . Specifically, (41) can be written as

VI. BUFFER OCCUPANCY PERCENTAGE
It is insightful to investigate the overall buffer occupancy at the relays. Evidently, the average queue length of the buffers can be used to determine the occupancy percentage:

VII. SIMULATION EXAMPLES
In this section, we provide several simulation examples to evaluate the performance of proposed secondary network in comparison to benchmarks MLRS [25] and EFPbRS [39]. However, since the aforementioned techniques do not consider several users and multiple antennas, we first extend them to our model in the following.

A. BENCHMARK SCHEMES 1) MLRS EXTENSION
Based on the original MLRS technique, the strongest relaying link among all available links is selected at each time step [25]. Hence, if χ sr(t) and χ rd(t) respectively show the set of available S-R and R-D links at the t-th time slot, the ete SNR under the extended MLRS technique in our model is given by where γ * S−R (χ sr(t) ) = max

2) EFPbLS EXTENSION
To account for the impact of buffer-availability, Manoj et al. developed a buffer-aware link selection policy where the transmissions are prioritized based on the number of relays with completely full and empty buffer states. Hence, we call this technique the Full-Empty Prioritization-based Link Selection (FEPbLS). 1) if one or more buffers are completely full, select the best R-D link. If the best R-D link fails to support the target rate, select the best S-R link among the relays with completely empty buffer. 2) if no buffer is completely full, select the best S-R link among the relays with completely empty buffer. 3) if neither of the previous cases is true, select the best relaying link among all the available links. Therefore, the ete SNR of the network under EFPbLS can be formulated as where χ sr(t) and χ rd(t) are respectively, the set of S-R and R-D links associated with the relays whose buffers are completely empty and full. Also, χ sr(t) , χ rd(t) , γ * S−R (χ sr(t) ) and γ * R−D (χ rd(t) ) have the same definitions given for (44).

B. SIMULATION RESULTS
All the Monte-Carlo results are obtained through averaging over 10 6 independent packet trials where the error bound between the analytical and Monte-Carlo values are no more than 0.1%. It is important to mention that the average packet delay is normalized to the length of the packet and therefore, it is given as the number of packets. Also, the network parameters used for each simulation example would be given in the caption.
The curves in Fig. 5 illustrate the effect of direct transmissions on the performance of the network. More specifically, Fig. 5(a) depicts the outage probability versus the strength of direct channels for different number of antennas at the destination. As can be seen, as the direct link strength improves in comparison to the relaying channels, the positive influence of additional antenna at SD becomes more pronounced to the extent that for L d = 6 and sd ≈ 45%( sr = rd ) the outage probability reduces more than 15%, i.e., from 0.03 to 0.025. Similarly, the packet delay curves in Fig 5(b) shows that such an increase in the number of direct channels and/or the links strength would result in a lower ete delay. To underline the effect of direct transmissions, in the  following simulation examples we consider a moderately weak direct links by setting sd ≈ 0.45% sr = rd .
The simulation examples illustrated in Fig. 6 depict the outage and delay performance of the network versus interference temperature in comparison to the MLRS and EFP-bLS techniques and for different numbers of users, relays, and buffering capacities. In specific, the curves in Fig. 6(a) clearly shows how increasing N , K , and M results in higher reliability. Furthermore, a careful comparison between the outage behavior of the proposed protocol and MLRS technique reveals that as the N , K , and M grows, the proposed technique substantially outperforms both of the MLRS and EFPbLS techniques; this can be explained by the fact that the proposed protocol was designed to preserve the availability of the relays and in turn, the higher packet arrival rate due to the proliferation of S-R channels is balanced out by more often activation of the R-D links. Although EFPbLS in contrast to MLRS is also aimed to address this issue, it is not able to preserve the same level of buffer-availability compared to the proposed protocol; this can be explained by the fact that EFPbLS is focused only on the completely full and empty states. Also, as number of users and relays increases the importance of the buffer-aware links scheduling becomes more conspicuous and a large performance gap occurs.
On the other hand, Fig. 6(b) depicts the average ete packet delay versus interference temperature in comparison to the benchmarks for the same different K ,N , M scenarios given in Fig. 6(a). As can be seen, the proposed protocol outperforms both of the MLRS and EFPbLS techniques. However, the overall impact of increasing K ,N , and M is a higher ete packet delay.
Furthermore, to shed light on the large performance gap between the proposed protocol and benchmarks in higher SNRs and bigger values of N , M , and K , in Fig. 7, we illustrate the average queue length of the buffers for the examples presented in Fig. 6 with N = M = K = 4. Evidently, as the transmit power budget increases in Fig. 7, i.e., I th → 0, the buffer occupancy percentage under MLRS technique does not improve conspicuously and lowers down slightly around 80% range. In the case of EFPbLS, the buffer occupancy decreases into the 70% range. However, the buffer-occupancy of the buffers under the proposed protocol quickly plunges into 50% range which clearly explains why it provides a superior outage and delay performance. More specifically, since each source and relay are equipped with several antennas the ramifications of an unavailable relay would be highly consequential as can be seen in the gap in the associated outage and delay curves in Fig. 6.
The curves in Fig. 8 reflect the achievable effective capacity for different target rates. As can be seen, the maximum achievable effective capacity under the proposed protocol is larger than that under the MLRS and EFPbLS techniques. This can be justified by the superior outage performance of the proposed technique because increasing the target rate would translate into the higher outage probability, lowering the achievable effective capacity. (Note: the effective capacity C eff can be readily evaluated as C eff = ρP o ete (ρ) where P o ete (ρ) is the outage probability if the target rate is ρ.) The simulation examples in Fig. 9 demonstrate the effect different fading profiles on the ete performance of the network where m represents the fading severity parameter for all the network channels. As the fading channels become stronger ( sr = rd → 1), it can be seen that the outage and delay performance of network improves under different severity profiles in Fig. 9(a) and Fig. 9(b), respectively. Moreover, as the m grows bigger, the outage performance degrades which is due to the fact that the behavior of outage probability over different values of m depends on the target rate of the network [70]. It is important to mention that, for the non-integer fading severity parameters, the outage probability as well as the packet delay can only be obtained using a numerical integration. In this context, the fading scenarios with m = 1.5 and m = 2.5 were considered in Fig. 6 to highlight this point. Furthermore, this figure clearly shows that the proposed protocol outperforms the benchmarks in different fading scenarios.
The simulation examples in Fig. 10 depicts the average required CSI circulation for link selection in terms of number of channels and under MLRS, EFPbLS, proposed protocol with buffer-aware CSI circulation, and the proposed protocol without buffer-aware CSI circulation (labeled as ''Proposed*''). This figure reveals two interesting facts. First: as can be seen, the proposed buffer-aware CSI circulation results in a lowered overhead compared to Proposed*. Second: in the high SNR regime the proposed protocol and the EFP-bLS would circulate almost the CSI of twelve and ten links respectively whereas the MLRS circulation slightly fluctuates around seven links. This is due to the fact that the MLRS results in a higher ce and cf states. Hence, the superior performance of the proposed protocol in terms of delay and reliability comes at the cost of higher overhead for CSI circulation. Furthermore, since all the three techniques would acquire the buffer-state information prior to each transmission slot, it is clear that the MLRS and EFPbLS would require 2K bits as the number of relays increase while the proposed technique would require 3K bits. These requirements stem from the fact that the buffer dynamics in MLRS and EFPbLS are quantized in three cases (ce, cf, ow) wheres the the buffer dynamics under the proposed technique would be quantized in five cases (ce, cf, lg, ld, ow).
Evidently, the amount of information that a certain communication protocol requires for operation reflects a measure of implementation complexity. (a larger overhead requires more resources, e.g., more bandwidth for feedback channels). Since the instantaneous required information can vary over time depending on the buffer status, it can be concluded that the proposed protocol is more complex than EFPbLS and MLRS in the average sense as illustrated in Fig. 10. Overall, since the proposed protocol can markedly outperform the MLRS and EFPbLS techniques in term of reliability and delay, its higher complexity would pale in significance if the prime goal is achieving a higher reliability and a lower delay.
The simulation examples in Fig. 11 illustrates the sensitivity of the proposed communication protocol with respect to the accuracy of the collected CSI in comparison to MLRS and EFPbLS benchmarks. Particularly, adopting the outdated CSI model for the interference channels, we investigated the outage and delay performance of the three protocols in Fig. 11(a) and Fig. 11(b). Mathematically, the outdated CSI for the interference channels can be modeled as [71] g xp = ηĝ xp + (1 − η) g xp (48) whereg xp is the outdated CSI of the channel between x and primary network such that x ∈ {s l n , r j k }, η is the correlation parameter determining the severity of the outdatedness,ĝ xp is random variable with the same statistical properties of the actual CSI g xp . (The higher η, the lower resemblance of the available CSIg xp to the actual CSI g xp .) The figures reveal three important facts. First: the proposed protocol has a higher sensitivity to the CSI accuracy as the gap between the curves is larger compared to the gaps in the benchmark schemes. Second, for the exact and slightly outdated CSI scenarios, the proposed protocol outperforms the benchmark schemes over the entire range of interference temperatures. Third, for the highly outdated CSI scenario the proposed protocol is outperformed by EFPbLS; Such a trend can be justified by the following fact: since EFbLS relies less on the CSI for link selection, it results in a lower performance loss in highly outdated CSI. This is an interesting observation because it illustrates that if the accurate CSI cannot be collected, then the dependency of the link selection protocol on the CSI must be relaxed. This can be investigated as an interesting future research direction.

VIII. CONCLUSION
In this paper, we jointly incorporated the multiple-antenna technology and buffer aided relaying to provide a highly reliable connectivity for the secondary network with several source users. To harness the potential benefits of direct channels, which inheres in such network model, we proposed a novel communication protocol which incorporates the direct transmissions along the relaying link to improve the reliability and average delay. In this regard, to mitigate the additional overhead that is caused due to CSI circulation of the direct channels, the link selection and CSI collection are performed based on the buffer-states of the relays. Considering an all-multi-antenna configuration and Nakagami-m fading channels, we derived closed-form expressions for the outage probability and average packet delay of the secondary network under the proposed protocol. Through several simulation examples, we demonstrated that a highly reliable connectivity for the secondary network can be achieved by increasing the number of antennas and buffering capacity. In this context, the proposed communication protocol can significantly outperform two benchmark schemes, namely, MLRS and EFPbLS in terms of reliability and delay. Furthermore, we demonstrate that if the global CSI of the network cannot be collected accurately, then the dependency of the link selection mechanism on the global CSI should be relaxed to mitigate performance loss. Also, the theoretical and Monte-Carlo results in the simulation examples are in excellent agreement which verifies the correctness of the derived closed-form expressions.
Here, according to the proposed protocol, it is clear that sr i = 0 when all the buffers are full and similarly rd i = 0 when all the buffers are empty. Also, in the rest of circumstances the number of available relays depends on the buffer status. In other words, letting i = 1 and i = (M + 1) K respectively represent the ce and cf cases, it is clear that for the ow case we have i = 1, (M + 1) K . Therefore, if i = 1, then we would only have S-R links available and in turn, i = P r γ SR < | sr i = K . On the other hand, if i = (M + 1) K , then we would only have R-D links available and therefore, i = P r γ rd < | rd i = K . However, if i = 1, (M + 1) K , then both of the S-R and R-D channels would be available; also, since γ SD and γ SR are correlated due to the common transmit antenna, i for this case can be simplified as i = A i B i where A i = P r max (γ SR , γ SD ) < | sr i and B i = P r γ RD < | rd i . Organizing such cases, (9) is obtained.
Appendix B: Placing the SNR expressions into (11) using (6) and considering the fact that the S-R channels are i.i.d RVs, A i can be written as x   f g s l n p (x) dx, (50) where W = L r j=1 g s l n r j k and Z = L d j=1 g s l n d j . Here, by applying the concept of total probability, we can rewrite expressions inside A i as the product of two independent probabilities conditioned on g s l n p , i.e., A i = E X {A i |x} where X = g s l n p and E X {.} represents the expectation operation. Next, breaking the integration region into two parts to decouple the minimum function and using (2), A i can be expressed as A i = (A 1 + A 2 ) L s N where A 1 and A 2 are given as (12).

Appendix C:
Rewording the two Gamma incomplete functions as a series representation, then applying the multinomial expansion, the exponential part of the integrand in A 2 can be express as where the definition of the new parameters at the right hand side of (51) are given in the paragraph after (14). Here, placing (51) and utilizing the series representation for the other Gamma incomplete function in (12) and then solving the resulting integration task, (14) is obtained. Appendix D: In this case, since the S-R channels are a set of i.i.d random variables, (18) can be rewritten as Here, placing the SNR expressions into (52) using (6) and then breaking the integration region into two parts to decouple the minimum function similar to (50), the expressions (19) and (20) are obtained.