Uplink Non-Orthogonal Multiple Access over Mixed RF-FSO Systems

TL;DR

This paper analyzes a relay-assisted uplink NOMA system with mixed RF and FSO links, deriving performance metrics and demonstrating the advantages of FSO backhaul for high data rates and reliability.

Contribution

It provides the first comprehensive analysis of uplink NOMA over mixed RF-FSO systems, including closed-form outage probability and ergodic capacity expressions.

Findings

FSO backhaul outperforms RF backhaul in high-throughput scenarios.

Derived closed-form outage probability expressions for mixed RF-FSO NOMA.

Showed the superiority of FSO links for high-reliability NOMA communications.

Abstract

In this paper, we consider a relay-assisted uplink non-orthogonal multiple access (NOMA) system. In this system, two radio frequency (RF) users are grouped for simultaneous transmissions, over each resource block, to an intermediate relay. The relay then forwards the amplified version of the users' aggregated signals, in the presence of multiuser interference, to a relatively far destination. In order to cope with the users' ever-increasing desire for higher data rates, a high-throughput free-space optics (FSO) link is employed as the relay-destination backhaul link. It is assumed that the FSO backhaul link is subject to Gamma-Gamma turbulence with pointing error. Also, a Rayleigh fading model is considered for the user-relay access links. Under these assumptions, we derive closed-form expressions for the outage probability and tractable forms, involving only one-dimensional integrals, for the ergodic capacity. Moreover, the outage probability and ergodic capacity analysis are extended to the conventional RF-backhauled systems in the presence of multiuser interference to both relay and destination nodes, and Rician fading for the relay-destination RF link. Our results reveal the superiority of FSO backhauling for high-throughput and high-reliability NOMA systems compared to RF backhauling. This work can be considered as a general analysis of dual-hop uplink NOMA systems as well as the first attempt to incorporate power-domain NOMA in mixed RF-FSO systems.

Tables

Table 1. Table I: Some of the important parameters used for simulations [ 36 ] .

Coefficient	Value
Responsivity of the photodetector, $\rho$	$0.5$ $V^{-1}$
Electrical-to-optical conversion coefficient, $\eta$	$1$
Optical receiver aperture radius, $r$	$10$ $\mathrm{cm}$
Laser beam divergence angle, $\varphi$	$2$ $\mathrm{mrad}$
Noise power at the relay and destination RF receivers, ${\sigma }_{ℛ}^2$ and $N_0$	$-80$ $\mathrm{dBm}$
Noise variance at the destination FSO receiver, ${\sigma }_{𝒟}^2$	${10}^{-14}$ $A^2$
Reference distance of the RF link, $d_{\mathrm{ref}}^{\mathrm{RF}}$	$80$ $m$
Number of interfering users to the relay, $K$	$10$
Number of interfering users to the destination, $K_d$	$10$
Transmitter and receiver antenna gains of the user-relay RF links, $(G_{t,i}^{\mathrm{RF}},G_{r,ℛ}^{\mathrm{RF}})$	$(5,8)$ $\mathrm{dBi}$
Transmitter and receiver antenna gains of the relay-destination RF backhaul link, $(G_{t,ℛ}^{\mathrm{RF}},G_{r,𝒟}^{\mathrm{RF}})$	$(10,15)$ $\mathrm{dBi}$
Number of iterations for numerical simulations, $N_t$	$5\times {10}^6$
Length of the FSO backhaul link, $d_{ℛ𝒟}$	$1200$ $m$
Wavelength of FSO signal, ${\lambda }^{\mathrm{FSO}}$	$1550$ $\mathrm{nm}$
Frequency of RF signal, $f^{\mathrm{RF}}$	$3$ $\mathrm{GHz}$
Parameter of the Rician distribution for the RF backhaul link, $\Omega$	$6$ $\mathrm{dB}$
Path-loss exponent of RF links, $\nu$	$3.5$

Equations

$$y_{\mathcal{R}}=\sum _{i=1}^2{x}_i{\tilde{h}}_i\sqrt{a_i{L}_{i}P}+\sum _{k=1}^K{x}_k^{\prime }{\tilde{h}}_k^{\prime }\sqrt{L_k^{\prime }{p}_k^{\prime }}+n_{\mathcal{R}},$$

$$\displaystyle\!f_{\tilde{g}}(\tilde{g})=\frac{\alpha\beta\xi^{2}}{A_{0}\Gamma(\alpha)\Gamma(\beta)}{\rm G}_{1,3}^{3,0}\left[\frac{\alpha\beta}{A_{0}}\tilde{g}{\bigg{|}\begin{matrix}\xi^{2}\\ \xi^{2}-1,\alpha\!-\!1,\beta\!-\!1\end{matrix}}\right],$$

$$y_{\mathcal{D}}=$$

$$+n_{\mathcal{D}},$$

$${\gamma }_{{\pi }_1}^{(1)}=\frac{a_1{L}_{1}P{\tilde{g}}^2|{\tilde{h}}_1{|}^2}{a_2{L}_{2}P{\tilde{g}}^2|{\tilde{h}}_2{|}^2+{\sum }_{k=1}^K{L}_k^{\prime }{p}_k^{\prime }{\tilde{g}}^2|{\tilde{h}}_k^{\prime }{|}^2+{\tilde{g}}^2{\sigma }_{\mathcal{R}}^2+C_{\mathcal{D}}},$$

$${\gamma }_{{\pi }_1}^{(2)}=\frac{a_2{L}_{2}P{\tilde{g}}^2|{\tilde{h}}_2{|}^2}{{\sum }_{k=1}^K{L}_k^{\prime }{p}_k^{\prime }{\tilde{g}}^2|{\tilde{h}}_k^{\prime }{|}^2+{\tilde{g}}^2{\sigma }_{\mathcal{R}}^2+C_{\mathcal{D}}},$$

$$\displaystyle\Pr(\gamma^{(1)}_{\pi_{1}}<\gamma^{(1)}_{\rm th},\pi_{1}){\Big{|}}_{\gamma^{(1)}_{\rm th}<1}$$

$$\displaystyle\stackrel{{\scriptstyle(b)}}{{=}}\operatorname{\mathbb{E}}_{|\tilde{h}_{2}|^{2}\!<\!{J^{(1)}_{\rm th}}\big{(}\boldsymbol{\mathcal{I}}_{1}\!+\!C_{\mathcal{D}}/(a_{1}L_{1}P\tilde{g}^{2})\big{)}}\!\!\left[\exp\left(-|\tilde{h}_{2}|^{2}\!\times\!10^{-s/10}\right)-\exp\left(-\gamma^{(1)}_{\rm th}\!\Big{[}|\tilde{h}_{2}|^{2}\!\times\!10^{-s/10}+{\boldsymbol{\mathcal{I}}_{1}}+C_{\mathcal{D}}/(a_{1}L_{1}P\tilde{g}^{2})\Big{]}\right)\right]$$

$$\displaystyle=\operatorname{\mathbb{E}}_{\boldsymbol{\mathcal{I}}_{1},\tilde{g}}\left[({1+10^{-s/10}})^{-1}\times\left[1-\exp\left(-\left[{1+10^{-s/10}}\right]\!{J^{(1)}_{\rm th}}\big{(}\boldsymbol{\mathcal{I}}_{1}\!+\!C_{\mathcal{D}}/(a_{1}L_{1}P\tilde{g}^{2})\big{)}\right)\right]\right]-\!\operatorname{\mathbb{E}}_{\boldsymbol{\mathcal{I}}_{1},\tilde{g}}\!\bigg{[}({1\!+\!\gamma^{(1)}_{\rm th}\!\!\times\!10^{-s/10}})^{-1}$$

$$\displaystyle\times\!\left[1\!-\!\exp\!\left(\!-\!\left[{1\!+\!\gamma^{(1)}_{\rm th}\!\!\times\!\!10^{-s/10}}\right]\!{J^{(1)}_{\rm th}}\big{(}\boldsymbol{\mathcal{I}}_{1}\!+\!C_{\mathcal{D}}/(a_{1}\!L_{1}\!P\tilde{g}^{2})\big{)}\right)\right]\!\exp\!\left(\!-\gamma^{(1)}_{\rm th}\!\Big{[}{\boldsymbol{\mathcal{I}}_{1}}\!+\!C_{\mathcal{D}}/(a_{1}\!L_{1}\!P\tilde{g}^{2})\Big{]}\!\right)\!\!\bigg{]}$$

$$\displaystyle\stackrel{{\scriptstyle(c)}}{{=}}\frac{10^{s/10}}{1+10^{s/10}}\left(1\!-\!\exp\!\bigg{(}\frac{-J^{(1)}_{{\rm th},1}\sigma^{2}_{\mathcal{R}}}{a_{1}L_{1}P}\bigg{)}\operatorname{\mathbb{E}}_{\tilde{g}}\!\!\left[\exp\!\bigg{(}\frac{-J^{(1)}_{{\rm th},1}C_{\mathcal{D}}}{a_{1}L_{1}P\tilde{g}^{2}}\bigg{)}\right]\prod_{k=1}^{K}\frac{a_{1}L_{1}P}{a_{1}L_{1}P+{J^{(1)}_{{\rm th},1}}L^{\prime}_{k}p^{\prime}_{k}}\right)\!-\!\frac{10^{s/10}}{\gamma^{(1)}_{\rm th}+10^{s/10}}\Bigg{(}\!\!\exp\!\bigg{(}\frac{-\gamma^{(1)}_{\rm th}\sigma^{2}_{\mathcal{R}}}{a_{1}L_{1}P}\bigg{)}$$

$$\displaystyle\times\!\operatorname{\mathbb{E}}_{\tilde{g}}\!\!\left[\exp\!\bigg{(}\frac{-\gamma^{(1)}_{\rm th}C_{\mathcal{D}}}{a_{1}L_{1}P\tilde{g}^{2}}\bigg{)}\right]\prod_{k=1}^{K}\frac{a_{1}L_{1}P}{a_{1}L_{1}P+{\gamma^{(1)}_{\rm th}}L^{\prime}_{k}p^{\prime}_{k}}-\exp\bigg{(}\frac{-\sigma^{2}_{\mathcal{R}}J^{(1)}_{{\rm th},2}}{a_{1}L_{1}P}\bigg{)}\operatorname{\mathbb{E}}_{\tilde{g}}\!\!\left[\exp\!\bigg{(}\frac{-C_{\mathcal{D}}J^{(1)}_{{\rm th},2}}{a_{1}L_{1}P\tilde{g}^{2}}\bigg{)}\right]\prod_{k=1}^{K}\frac{a_{1}L_{1}P}{a_{1}L_{1}P+L^{\prime}_{k}p^{\prime}_{k}J^{(1)}_{{\rm th},2}}\Bigg{)}.$$

$$P_{\mathrm{out}}^{(1)}$$

$$\displaystyle\stackrel{{\scriptstyle(b)}}{{=}}1\!-\!\Big{[}\Pr(\gamma^{(1)}_{\pi_{1}}>\gamma^{(1)}_{\rm th},\pi_{1})+$$

$$\displaystyle\Pr(\gamma^{(2)}_{\pi_{2}}>\gamma^{(2)}_{\rm th},\pi_{2})\times\Pr(\gamma^{(1)}_{\pi_{2}}>\gamma^{(1)}_{\rm th},\pi_{2})/P(\pi_{2})\Big{]},\!$$

$$P({\pi }_1)$$

$$=E_{|{\tilde{h}}_2{|}^2}[\exp (-|{\tilde{h}}_2{|}^2\times 10^{-s/10})]$$

$$=(1+10^{-s/10}{)}^{-1},$$

$$P_{\mathrm{cov}|{\pi }_1}^{(1)}$$

$$P_{\mathrm{cov}|{\pi }_2}^{(1)}$$

$$\displaystyle\operatorname{\mathbb{E}}_{\boldsymbol{\mathcal{I}}_{1},\tilde{g}}\left[\exp\left(-C\big{[}\boldsymbol{\mathcal{I}}_{1}+C_{\mathcal{D}}/(a_{1}L_{1}P\tilde{g}^{2})\big{]}\right)\right]=$$

$$\displaystyle\!\exp\!\bigg{(}\!\frac{-C\sigma^{2}_{\mathcal{R}}}{a_{1}L_{1}P}\!\bigg{)}\operatorname{\mathbb{E}}_{\tilde{g}}\!\!\left[\exp\!\bigg{(}\!\frac{-CC_{\mathcal{D}}}{a_{1}\!L_{1}\!P\tilde{g}^{2}}\!\!\bigg{)}\!\right]\prod_{k=1}^{K}\!\frac{a_{1}L_{1}P}{a_{1}\!L_{1}\!P\!+\!C\!L^{\prime}_{k}p^{\prime}_{k}},\!$$

$$\displaystyle\Pr(\gamma^{(1)}_{\pi_{1}}<\gamma^{(1)}_{\rm th},\pi_{1}){\Big{|}}_{\gamma^{(1)}_{\rm th}\geq 1}\!\!{=}\frac{10^{s/10}}{1\!+\!10^{s/10}}-\frac{10^{s/10}}{\gamma^{(1)}_{\rm th}\!+\!10^{s/10}}\exp\!\bigg{(}\!\frac{-\gamma^{(1)}_{\rm th}\sigma^{2}_{\mathcal{R}}}{a_{1}L_{1}P}\!\bigg{)}\!\operatorname{\mathbb{E}}_{\tilde{g}}\!\!\left[\exp\!\bigg{(}\frac{-\gamma^{(1)}_{\rm th}C_{\mathcal{D}}}{a_{1}L_{1}P\tilde{g}^{2}}\bigg{)}\right]\prod_{k=1}^{K}\frac{a_{1}L_{1}P}{a_{1}L_{1}P\!+\!{\gamma^{(1)}_{\rm th}}L^{\prime}_{k}p^{\prime}_{k}}.$$

$$\displaystyle\mathcal{G}(A)\!\triangleq\!\operatorname{\mathbb{E}}_{\tilde{g}}\!\left[\exp\!\left(\!-\frac{A}{\tilde{g}^{2}}\right)\!\right]\!\!=\!\frac{\xi^{2}\!\times\!2^{\alpha+\beta-2}}{2\pi\Gamma(\alpha)\Gamma(\beta)}{\rm G}_{6,1}^{0,6}\left[\frac{16A_{0}^{2}}{A(\alpha\beta)^{2}}{\bigg{|}\begin{matrix}1,(2-\xi^{2})/2,(1-\alpha)/2,(2-\alpha)/2,(1-\beta)/2,(2-\beta)/2\\ -\xi^{2}/2\end{matrix}}\right].$$

$$\displaystyle\Pr(\gamma^{(2)}_{\pi_{2}}<\gamma^{(2)}_{\rm th},\pi_{2}){\Big{|}}_{\gamma^{(2)}_{\rm th}<1}\!$$

$$\displaystyle\times\Bigg{[}\exp\!\bigg{(}\frac{-\gamma^{(2)}_{\rm th}\sigma^{2}_{\mathcal{R}}}{a_{2}L_{2}P}\bigg{)}\mathcal{G}\!\bigg{(}\frac{\gamma^{(2)}_{\rm th}C_{\mathcal{D}}}{a_{2}L_{2}P}\bigg{)}\prod_{k=1}^{K}\frac{a_{2}L_{2}P}{a_{2}L_{2}P+{\gamma^{(2)}_{\rm th}}L^{\prime}_{k}p^{\prime}_{k}}-\exp\!\bigg{(}\frac{-\sigma^{2}_{\mathcal{R}}J^{(2)}_{{\rm th},2}}{a_{2}L_{2}P}\bigg{)}\mathcal{G}\!\bigg{(}\frac{C_{\mathcal{D}}J^{(2)}_{{\rm th},2}}{a_{2}L_{2}P}\bigg{)}\prod_{k=1}^{K}\frac{a_{2}L_{2}P}{a_{2}L_{2}P+L^{\prime}_{k}p^{\prime}_{k}J^{(2)}_{{\rm th},2}}\Bigg{]}.$$

$$\displaystyle\Pr(\gamma^{(2)}_{\pi_{2}}<\gamma^{(2)}_{\rm th},\pi_{2}){\Big{|}}_{\gamma^{(2)}_{\rm th}\geq 1}\!\!{=}\frac{1}{1\!+\!10^{s/10}}-\frac{10^{-s/10}}{\gamma^{(2)}_{\rm th}\!+\!10^{-s/10}}\exp\!\bigg{(}\!\frac{-\gamma^{(2)}_{\rm th}\sigma^{2}_{\mathcal{R}}}{a_{2}L_{2}P}\!\bigg{)}\mathcal{G}\!\bigg{(}\frac{\gamma^{(2)}_{\rm th}C_{\mathcal{D}}}{a_{2}L_{2}P}\bigg{)}\prod_{k=1}^{K}\frac{a_{2}L_{2}P}{a_{2}L_{2}P\!+\!{\gamma^{(2)}_{\rm th}}L^{\prime}_{k}p^{\prime}_{k}}.$$

$$\displaystyle\!\!\!\!\!\Pr(\gamma^{(1)}_{\pi_{2}}\!>\!\gamma^{(1)}_{\rm th}\!,\pi_{2})\!=\!\frac{1}{1\!+\!10^{s/10}}\!\times\!\exp\!\bigg{(}\!\!\frac{-\gamma^{(1)}_{\rm th}\sigma^{2}_{\mathcal{R}}(1\!+\!10^{s/10})}{a_{1}L_{1}P}\bigg{)}\mathcal{G}\!\bigg{(}\frac{\gamma^{(1)}_{\rm th}C_{\mathcal{D}}(1\!+\!10^{s/10})}{a_{1}L_{1}P}\bigg{)}\!\!\prod_{k=1}^{K}\!\frac{a_{1}L_{1}P}{a_{1}L_{1}P+{\gamma^{(1)}_{\rm th}}L^{\prime}_{k}p^{\prime}_{k}(1\!+\!10^{s/10})}.$$

$$P_{\mathrm{out}}^{(2)}$$

$$\displaystyle\Pr(\gamma^{(1)}_{\pi_{1}}>\gamma^{(1)}_{\rm th},\pi_{1})\times\Pr(\gamma^{(2)}_{\pi_{1}}>\gamma^{(2)}_{\rm th},\pi_{1})/P(\pi_{1})\Big{]},$$

$$\mathrm{Pr}({\gamma }_{{\pi }_1}^{(2)}>{\gamma }_{\mathrm{th}}^{(2)},{\pi }_1)=(1+10^{-s/10}{)}^{-1}\times$$

$$\displaystyle\exp\!\bigg{(}\frac{-\gamma^{(2)}_{\rm th}\sigma^{2}_{\mathcal{R}}(1\!+\!10^{-s/10})}{a_{2}L_{2}P}\bigg{)}\mathcal{G}\!\bigg{(}\frac{\gamma^{(2)}_{\rm th}C_{\mathcal{D}}(1\!+\!10^{-s/10})}{a_{2}L_{2}P}\bigg{)}$$

$$\times \prod _{k=1}^K\frac{a_2{L}_{2}P}{a_2{L}_{2}P+{\gamma }_{\mathrm{th}}^{(2)}{L}_k^{\prime }{p}_k^{\prime }(1+10^{-s/10})}.$$

$$R_{\Sigma }=$$

$$=$$

$$P_{\mathrm{out}}^{\Sigma }=\mathrm{Pr}\left(|{\tilde{h}}_1{|}^2<{\gamma }_{\mathrm{th}}^{\Sigma }\left[{\mathcal{I}}_1+\frac{C_{\mathcal{D}}}{a_1{L}_{1}P{\tilde{g}}^2}\right]-|{\tilde{h}}_2{|}^2\times 10^{-s/10}\right).$$

$$P_{\mathrm{out}}^{\Sigma }=$$

$$=$$

$$\displaystyle{\hskip 5.69046pt}-\frac{10^{s/10}}{10^{s/10}-1}\times\exp\!\bigg{(}\frac{-\gamma^{\Sigma}_{\rm th}\sigma^{2}_{\mathcal{R}}}{a_{1}L_{1}P}\bigg{)}{\mathcal{G}}\!\bigg{(}\frac{\gamma^{\Sigma}_{\rm th}C_{\mathcal{D}}}{a_{1}L_{1}P}\bigg{)}\prod_{k=1}^{K}\frac{a_{1}L_{1}P}{a_{1}L_{1}P+\gamma^{\Sigma}_{\rm th}L^{\prime}_{k}p^{\prime}_{k}}.$$

$$f_{{\kappa }_b}({\kappa }_b)=\frac{1+\Omega }{{\mathrm{e}}^{\Omega }}\exp \left(-(1+\Omega ){\kappa }_b\right)I_0\left(2\sqrt{\Omega (1+\Omega ){\kappa }_b}\right),$$

Full text

I Introduction

Non-orthogonal multiple access (NOMA) is widely considered as one of the enabling technologies for the fifth generation (5G) wireless networks. With its two general power- and code-domain forms, NOMA can potentially pave the way toward higher throughput, lower latency, improved fairness, higher reliability, and massive connectivity [2]. Motivated by these fascinating advantages, extensive research activities have been carried out in the past few years to advance NOMA in diverse directions [3, 4].

Excavating the rich literature on NOMA, some research activities have focused on the code-domain NOMA which, generally speaking, attempts to serve a set of users in a smaller set of orthogonal resource blocks using a pattern matrix. In this context, variety of multiple access mechanisms have been proposed such as sparse code multiple access (SCMA) [5], lattice partition multiple access (LPMA) [6], interleave-grid multiple access (IGMA) [7], and pattern division multiple access (PDMA) [8]. However, code-domain NOMA usually suffers from a high detection complexity due to the need for complex multiuser detection methods such as maximum likelihood (ML) detection, massage passing algorithm (MPA), and maximum a posteriori (MAP) detection. To this end, some of the recent work in code-domain NOMA have focused on lowering the detection complexity, e.g., using list sphere decoding [9], and recursive detection approaches enabled by sophisticated designs of the overall pattern matrix [10].

Power-domain NOMA, on the other hand, has attracted more attention from the research and industrial communities because of its relative simplicity. For example, a two-user downlink power-domain NOMA system, also referred to as multiuser superposition transmission, has been proposed for inclusion in the Third Generation Partnership Project (3GPP) long-term evolution advanced (LTE-A) standard [11]. Motivated by this, power-domain NOMA has been explored in various directions, including multiple-input multiple-output systems [12], cooperative transmission [13], simultaneous wireless information and power transfer [14, 15, 16], ultra-reliable and low-latency communications [17], visible-light communications [18], and millimeter-wave (mmWave) communications [19], (see, e.g., [3] for a comprehensive survey).

In a variety of applications, there is a need to transmit the users’ data to a central unit or a wired base station (BS); however, given the limited power of the users, it is not feasible for the users to directly communicate with the relatively far destination. To this end, several recent works have considered the relaying problem in downlink and uplink NOMA communications. In particular, capacity analysis of a simple cooperative relaying system, consisting of a source, a relay, and a destination node is provided in [20]. The outage probabilities and ergodic sum rate of a downlink two-user NOMA system, with a full-duplex relay helping one of the users, are characterized in [21]. Performance of downlink NOMA transmission with an intermediate amplify-and-forward (AF) relay for multiple-antenna systems, and over Nakagami-$m$ fading channels is investigated in [22] and [23], respectively. The performance of coordinated direct and relay transmission for two-user downlink and uplink NOMA systems is investigated in [24] and [25], respectively. Hybrid decode-and-forward (DF) and AF relaying in NOMA systems is proposed in [26], and forwarding strategy selection (either AF or DF) problem is explored in [27]. Moreover, a comprehensive performance evaluation of a two-user downlink cooperative NOMA system is provided in [28], where one of the users acts as a relay switching between half-duplex and full-duplex modes. In addition, a dynamic DF-based cooperative scheme has been proposed in [29] for downlink NOMA transmission with spatially random users.

The aforementioned prior works often assume communications in the absence of external multiuser interference to the NOMA users and the conventional sub-6 GHz radio frequency (RF) for the backhaul links. However, the available bandwidth in the sub-6 GHz band is scarce and falls short of supporting the users’ aggressive demand for the higher data rates, especially when NOMA is employed in the user-relay access links to provide higher throughputs. In this case, the relay-destination backhaul link can pose a severe bottleneck on the end-to-end performance and substantially negate the NOMA advantages through reducing the users’ achievable throughput and reliability which can in turn even increase their latency.

A potential approach to overcome the aforementioned drawback is to utilize higher frequency bands, e.g., through the deployment of mmWave and/or free-space optics (FSO) backhaul links [30, 31]. MmWave communication is usually preferred for relatively short communication lengths due to the severe propagation conditions at millimeter frequencies [32]. FSO links, on the other hand, can provide much more available bandwidth and support ranges in the order of several kilometers [33]. To this end, in this paper, we investigate the performance of uplink NOMA transmission over mixed RF-FSO systems. In particular, an AF relay is employed to forward the amplified received signal from the Rayleigh fading access links to the destination through an ultra high-throughput directive interference-free FSO link subject to Gamma-Gamma (GG) fading with beam misalignment error. This paper can be considered as a general analysis of dual-hop uplink NOMA systems, and also an initial attempt to incorporate power-domain NOMA in mixed RF-FSO systems.

Our main contributions can be summarized as follows.

• •

We consider a general dual-hop uplink NOMA transmission subject to the presence of multiuser interference from some independent users. Such interference can be induced, e.g., due the co-channel interference from nearby users aiming to communicate to other relays or destinations. The inclusion of external interference can be also helpful in the analysis of mmWave NOMA, where side-lobes of nearby mmWave beams cause inter-beam interference to the power-domain NOMA users grouped over a given mmWave beam [34, 35].

• •

We apply dynamic-order decoding to determine the detection order of the NOMA users at the destination.

• •

We derive the closed-form expressions for the individual- and sum-rate outage probabilities of the mixed RF-FSO uplink NOMA system with respect to dynamic-order decoding at the destination, AF scheme at the relay, Rayleigh fading for the user-relay access links, and GG turbulence with the inclusion of pointing error for the relay-destination backhaul FSO link.

• •

We further derive the outage probability closed forms (in terms of an infinite series that can effectively be approximated by some finite number of terms) for the RF-backhauled system when the relay-destination backhaul link is subject to Rician fading, and both relay and destination nodes are subject to external multiuser interference.

• •

Average individual- and sum-rate formulas are characterized for both FSO- and RF-backhauled systems up to only one-dimensional integrals over the fading coefficient of the backhaul link. That is equivalent to say that the ergodic capacity closed-form expressions are obtained for the single-hop uplink NOMA subject to some exterior multiuser interference, or better to say, for the aforementioned dual-hop system model given each realization of the backhaul fading coefficient.

• •

Extensive numerical results are provided to validate the accuracy of the derived formulas and also ascertain the system performance over different channel conditions.

The rest of the paper is organized as follows. In Section II, we describe the system model. In Section III, we derive the individual- and sum-rate outage probability closed-form formulas for both FSO- and RF-backhauled dual-hop uplink NOMA systems. Section IV is devoted to the ergodic capacity analysis for the same system model, Section V provides the numerical results, and Section VI concludes the paper.

II System Model

As shown in Fig. 1, consider two RF users $\mathcal{U}_{1}$ and $\mathcal{U}_{2}$ grouped together for uplink NOMA transmission to an AF relay $\mathcal{R}$. Denote the composite RF channel gain of the $\mathcal{U}_{i}-\mathcal{R}$ link by $h_{i}=\sqrt{L_{i}}\tilde{h}_{i}$, $i=1,2$, where $L_{i}$ and $\tilde{h}_{i}$ are the path-loss gain and the independent-and-identically-distributed (i.i.d.) Rayleigh fading coefficient of the $\mathcal{U}_{i}-\mathcal{R}$ RF link, respectively. The path-loss gain is given by $L_{i}={G^{\mathrm{RF}}_{t,i}G^{\mathrm{RF}}_{r,\mathcal{R}}}\times\left[{\lambda^{\mathrm{RF}}}/{(4\pi d^{\mathrm{RF}}_{\mathrm{ref}})}\right]^{2}\times\left[{d^{\mathrm{RF}}_{\mathrm{ref}}}/{d_{i}^{\mathrm{RF}}}\right]^{\nu}$ [36, Eq. (2)] in which $\lambda^{\mathrm{RF}}$ is the wavelength of the RF signal, $G^{\mathrm{RF}}_{t,i}$ and $G^{\mathrm{RF}}_{r,\mathcal{R}}$ are the RF transmit and receive antenna gains of the $\mathcal{U}_{i}-\mathcal{R}$ link, respectively, $d^{\mathrm{RF}}_{\mathrm{ref}}$ is a reference distance for the antenna far-field, $d^{\mathrm{RF}}_{i}$ is the $\mathcal{U}_{i}-\mathcal{R}$ link length, and $\nu$ is the RF path-loss exponent. Furthermore, assume that the uplink transmission to the relay is affected by undesired multiuser interference from $K$ interfering users $\mathcal{I}_{k}$, $k=1,2,...,K$, each with the transmit power $p^{\prime}_{k}$, path-loss gain $L^{\prime}_{k}$, and i.i.d. Rayleigh fading coefficient $\tilde{h}^{\prime}_{k}$. As explained in Section I, this interference can be from the users scheduled for the concurrent transmission to some other relays in the cellular network or any other non-vanishing interference during the desired transmission block. The received signal by the relay can then be expressed as

[TABLE]

where $x_{i}$ and $x^{\prime}_{k}$ are the transmit symbols by $\mathcal{U}_{i}$ and $\mathcal{I}_{k}$, respectively, $a_{1}$ and $a_{2}=1-a_{1}$ are the power coefficients, and $n_{\mathcal{R}}$ is the additive white Gaussian noise (AWGN) of the relay receiver with mean zero and variance $\sigma^{2}_{\mathcal{R}}$. Note that for the users with independent Rayleigh fading, all fading gains $|\tilde{h}_{i}|^{2}$’s and $|\tilde{h}^{\prime}_{k}|^{2}$’s have an exponential distribution with mean one (to ensure that fading neither amplifies nor attenuates the received power) as $f_{|\tilde{h}_{i}|^{2}}(x)=f_{|\tilde{h}^{\prime}_{k}|^{2}}(x)=\exp(-x)$, $x\geq 0$.

The received signal $y_{\mathcal{R}}$ at the relay is then converted to optical signal using intensity-modulation direct-detection (IM/DD), and is amplified with a constant gain $G$ to keep the disparity between the power levels of different NOMA users for successive interference cancellation (SIC) detection at the destination. In this case, the transmitted optical signal by the relay toward the destination $\mathcal{D}$ can be expressed as $S_{\mathcal{R}}=G(1+\eta y_{\mathcal{R}})$, where $\eta$ is the electrical-to-optical conversion coefficient [37]. The transmitted signal then undergoes the FSO channel with the composite gain $g=g_{l}\tilde{g}$ where $g_{l}$ is the path-loss gain of the $\mathcal{R-D}$ FSO backhaul link, with the length $d_{\mathcal{RD}}$, defined as $g_{l}=\rho\times 10^{-\kappa d_{\mathcal{RD}}/10}$ where $\rho$ is the responsivity of the photodetector, and $\kappa$ is the weather-dependent attenuation coefficient [36]. Moreover, $\tilde{g}=g_{p}g_{f}$ is the total fading coefficient due to pointing error $g_{p}$ and optical turbulence $g_{f}$. In the case of GG optical turbulence with beam misalignment, the distribution of $\tilde{g}$ can be expressed as [38]

[TABLE]

where $\alpha$ and $\beta$ are the fading parameters of the GG distribution, $\xi$ is the ratio of the equivalent beam radius and the pointing error displacement standard deviation (jitter) measured at the receiver, $\Gamma(\cdot)$ is the gamma function [39, Eq. (8.310)], and ${\rm G}[\cdot]$ is the Meijer’s G-function [39, Eq. (9.301)]. Furthermore, $A_{0}$ is the geometric loss in the case of perfect beam alignment (zero radial displacement) defined as $A_{0}=[{\rm erf}(\sqrt{\pi}r/(\sqrt{2}\phi d_{\mathcal{RD}}))]^{2}$ in which ${\rm erf}(\cdot)$ is the error function, $r$ is the receiver aperture radius, and $\phi$ is the transmitter beam divergence angle.

The destination then filters out the direct current (DC) component of $g_{l}\tilde{g}G$ from $g_{l}\tilde{g}S_{\mathcal{R}}+n_{\mathcal{D}}$ to obtain the received signal as

[TABLE]

where $n_{\mathcal{D}}$ is the destination AWGN with mean zero and variance $\sigma^{2}_{\mathcal{D}}$.

We assume that the NOMA users are indexed based on their path-loss gains, i.e., $L_{1}\geq L_{2}$, and the power allocation strategy proposed in [40] is adopted to determine $a_{1}$ and $a_{2}$ as $a_{1}L_{1}=a_{2}L_{2}\times 10^{s/10}$ where $s\geq 0$ is the power back-off step; hence, $a_{1}=L_{2}\times 10^{s/10}/(L_{1}+L_{2}\times 10^{s/10})$ and $a_{2}=L_{1}/(L_{1}+L_{2}\times 10^{s/10})$. We further consider dynamic-order decoding at the destination [41, 42], that is expected to achieve a higher performance compared to fixed-order decoding assuming that the BS has perfect knowledge about the channel state information (CSI) and sorts the NOMA users based on their instantaneous received power. In fact, based on the principles of uplink power-domain NOMA [40, 43], the BS sorts the users based on their channel conditions from the best to the worst. Therefore, depending on the fading coefficients $\tilde{h}_{1}$ and $\tilde{h}_{2}$, the detection order is $\pi_{1}=(1,2)$, meaning that the first user is decoded first, if $a_{1}L_{1}|\tilde{h}_{1}|^{2}\geq a_{2}L_{2}|\tilde{h}_{2}|^{2}$; otherwise, the detection order is $\pi_{2}=(2,1)$.

III Outage Probability Analysis

In this section, we first characterize the individual- and sum-rate outage probabilities for mixed RF-FSO NOMA systems and then extend the results to the RF-backhauled systems.

III-A Individual-Rate Outage Analysis

Note that if the detection order is $\pi_{1}$, the SIC receiver first treats the signal from the second NOMA user as noise to decode $x_{1}$ with the signal-to-interference-plus-noise ratio (SINR) given by

[TABLE]

and then, after removing the received power from the first user, decodes $x_{2}$ with the SINR given by

[TABLE]

where $C_{\mathcal{D}}\triangleq\sigma^{2}_{\mathcal{D}}/(\eta^{2}g_{l}^{2}G^{2})$. Similarly, when the detection order is $\pi_{2}$ the SINR values $\gamma^{(1)}_{\pi_{2}}$ and $\gamma^{(2)}_{\pi_{2}}$ can be obtained by properly changing the indices in (4) and (5).

Let $\gamma^{(i)}_{\rm th}=2^{R^{(i)}_{\rm th}}-1$ denote the threshold SINR for an IM/DD FSO link to achieve the desired data rate $R^{(i)}_{\rm th}$, $i=1,2$. Then the outage probability of the first user $\mathcal{U}_{1}$ in achieving an individual rate of $R^{(1)}_{\rm th}$ can be characterized as

[TABLE]

where step $(a)$ follows from the law of total probability by defining $P_{{\rm out}|\pi_{i}}^{(1)}$, $i=1,2$, as the conditional outage probability of the first NOMA user given the decoding order $\pi_{i}$. Moreover,

[TABLE]

and $P(\pi_{2})=1-P(\pi_{1})=(1+10^{s/10})^{-1}$ are the probabilities of having decoding orders $\pi_{1}$ and $\pi_{2}$, respectively. Furthermore, step $(b)$ follows, first, by defining $P_{{\rm cov}|\pi_{i}}^{(1)}\triangleq 1-P_{{\rm out}|\pi_{i}}^{(1)}$, $i=1,2$, as the probability of successfully achieving $R^{(1)}_{\rm th}$ for $\mathcal{U}_{1}$ conditioned on the decoding order $\pi_{i}$, and then noting that the correct detection of $x_{1}$ for the decoding order $\pi_{2}$ also requires successful decoding of the preceding symbol $x_{2}$, i.e.,111More precisely, $P_{{\rm cov}|\pi_{2}}^{(1)}$ in (III-A) should be written in the form of the joint probability $P_{{\rm cov}|\pi_{2}}^{(1)}=\Pr(\gamma^{(2)}_{\pi_{2}}>\gamma^{(2)}_{\rm th},\gamma^{(1)}_{\pi_{2}}>\gamma^{(1)}_{\rm th}|\pi_{2})=\Pr(\gamma^{(2)}_{\pi_{2}}>\gamma^{(2)}_{\rm th}|\pi_{2})\times\Pr(\gamma^{(1)}_{\pi_{2}}>\gamma^{(1)}_{\rm th}|\gamma^{(2)}_{\pi_{2}}>\gamma^{(2)}_{\rm th},\pi_{2})$ since the events $\{\gamma^{(1)}_{\pi_{2}}>\gamma^{(1)}_{\rm th}\}$ and $\{\gamma^{(2)}_{\pi_{2}}>\gamma^{(2)}_{\rm th}\}$ are not independent due to the presence of multiuser interference and the backhaul link imposing common random variables $\boldsymbol{\mathcal{I}}_{1}$ (defined in (III-A)) and $\tilde{g}$ on both $\gamma^{(1)}_{\pi_{2}}$ and $\gamma^{(2)}_{\pi_{2}}$. However, these two events are independent conditioned on $\boldsymbol{\mathcal{I}}_{1}$ and $\tilde{g}$. Therefore, we can calculate $\Pr(\gamma^{(1)}_{\pi_{2}}>\gamma^{(1)}_{\rm th}|\gamma^{(2)}_{\pi_{2}}>\gamma^{(2)}_{\rm th},\pi_{2})=\operatorname{\mathbb{E}}_{\boldsymbol{\mathcal{I}}_{1},\tilde{g}}\left[\Pr(\gamma^{(1)}_{\pi_{2}}>\gamma^{(1)}_{\rm th}|\gamma^{(2)}_{\pi_{2}}>\gamma^{(2)}_{\rm th},\pi_{2},\boldsymbol{\mathcal{I}}_{1},\tilde{g})\right]=\operatorname{\mathbb{E}}_{\boldsymbol{\mathcal{I}}_{1},\tilde{g}}\left[\Pr(\gamma^{(1)}_{\pi_{2}}>\gamma^{(1)}_{\rm th}|\pi_{2},\boldsymbol{\mathcal{I}}_{1},\tilde{g})\right]=\Pr(\gamma^{(1)}_{\pi_{2}}>\gamma^{(1)}_{\rm th}|\pi_{2})$.

[TABLE]

In the following, we calculate the three joint probabilities in (III-A) to ascertain the outage probability of the first user $\mathcal{U}_{1}$.

In order to calculate $\Pr(\gamma^{(1)}_{\pi_{1}}>\gamma^{(1)}_{\rm th},\pi_{1})$ we first note that $\Pr(\gamma^{(1)}_{\pi_{1}}>\gamma^{(1)}_{\rm th},\pi_{1})=P(\pi_{1})\Pr(\gamma^{(1)}_{\pi_{1}}>\gamma^{(1)}_{\rm th}|\pi_{1})=P(\pi_{1})[1-\Pr(\gamma^{(1)}_{\pi_{1}}<\gamma^{(1)}_{\rm th}|\pi_{1})]=P(\pi_{1})-\Pr(\gamma^{(1)}_{\pi_{1}}<\gamma^{(1)}_{\rm th},\pi_{1})$. Then using (4), $\Pr(\gamma^{(1)}_{\pi_{1}}<\gamma^{(1)}_{\rm th},\pi_{1})$ for $\gamma^{(1)}_{\rm th}<1$ can be calculated as (III-A) shown at the top of this page where, in step $(a)$, ${\boldsymbol{\mathcal{I}}_{1}}\triangleq(\sum_{k=1}^{K}{L^{\prime}_{k}p^{\prime}_{k}}|\tilde{h}^{\prime}_{k}|^{2}+\sigma^{2}_{\mathcal{R}})/(a_{1}L_{1}P)$ is the sum of the power of multiuser interference and noise, at the relay, normalized to the average power of the first NOMA user. Moreover, step $(b)$ follows, first, by defining the constant $J^{(1)}_{\rm th}\triangleq 10^{s/10}\times\gamma^{(1)}_{\rm th}/(1-\gamma^{(1)}_{\rm th})>0$ for $\gamma^{(1)}_{\rm th}<1$, and then noting that $\Pr(X<Y,X\geq Z)$ for three random variables (RVs) $X$, $Y$, and $Z$ can be calculated using the law of total probability as $\Pr(X<Y,X\geq Z)=\Pr(Z\leq X<Y,Z<Y)$ since $\Pr(Z\leq X<Y,Z\geq Y)=0$. Finally, step $(c)$ of (III-A) follows by noting that for any constant $C$

[TABLE]

due to the independence of ${\boldsymbol{\mathcal{I}}}_{1}$ and $\tilde{g}$, and then applying the independency among $|\tilde{h}^{\prime}_{k}|^{2}$’s to get \operatorname{\mathbb{E}}_{\boldsymbol{\mathcal{I}}_{1}}\!\big{[}\exp(-C\boldsymbol{\mathcal{I}}_{1})\big{]}=\exp\big{(}\frac{-C\sigma^{2}_{\mathcal{R}}}{a_{1}L_{1}P}\big{)}\prod_{k=1}^{K}\operatorname{\mathbb{E}}_{|{\tilde{h}^{\prime}}_{k}|^{2}}\left[\exp\left(-CL^{\prime}_{k}p^{\prime}_{k}|\tilde{h}^{\prime}_{k}|^{2}/(a_{1}L_{1}P)\right)\right]. Furthermore, in step $(c)$ of (III-A), $J^{(1)}_{{\rm th},1}\triangleq J^{(1)}_{\rm th}({1+10^{-s/10}})$ and $J^{(1)}_{{\rm th},2}\triangleq\gamma^{(1)}_{\rm th}+J^{(1)}_{\rm th}({1+\gamma^{(1)}_{\rm th}\times 10^{-s/10}})$.

We should further emphasize that (III-A) is obtained for $\gamma^{(1)}_{\rm th}<1$. If $\gamma^{(1)}_{\rm th}\geq 1$, the upper limit of $|\tilde{h}_{1}|^{2}$ in the equality $(a)$ of (III-A) is always greater than its lower limit meaning that the condition $(1-\gamma^{(1)}_{\rm th})|\tilde{h}_{2}|^{2}\times 10^{-s/10}<\gamma^{(1)}_{\rm th}\left[{\boldsymbol{\mathcal{I}}_{1}}+C_{\mathcal{D}}/(a_{1}L_{1}P\tilde{g}^{2})\right]$ holds for all values of $|\tilde{h}_{2}|^{2}$ and there is no need to impose such an extra condition on the derivation of the corresponding probability. Consequently, by averaging over $|\tilde{h}_{2}|^{2}$, $\boldsymbol{\mathcal{I}}_{1}$, and $\tilde{g}$, $\Pr(\gamma^{(1)}_{\pi_{1}}<\gamma^{(1)}_{\rm th},\pi_{1})$ for $\gamma^{(1)}_{\rm th}\geq 1$ can be derived as (11) shown at the top of the next page.

Finally, using (III-A) for $\gamma^{(1)}_{\rm th}<1$ or (11) for $\gamma^{(1)}_{\rm th}\geq 1$, one can obtain the coverage probability of the first NOMA user for the decoding order $\pi_{1}$ as $P_{\rm cov}^{(1)}(\pi_{1})=\Pr(\gamma^{(1)}_{\pi_{1}}>\gamma^{(1)}_{\rm th},\pi_{1})=(1+10^{-s/10})^{-1}-\Pr(\gamma^{(1)}_{\pi_{1}}<\gamma^{(1)}_{\rm th},\pi_{1})$. However, the closed-form characterization of $P_{\rm cov}^{(1)}(\pi_{1})$ still requires the calculation of expressions of the form $\operatorname{\mathbb{E}}_{\tilde{g}}\!\left[\exp\left(-A/\tilde{g}^{2}\right)\right]$, where $A$ is a constant and $\tilde{g}$ is distributed according to (2). To do so, we first apply [44, Eq. (11)] and [39, Eq. (9.31.2)] to write $\exp\left(-A/\tilde{g}^{2}\right)$ in the form of a Meijer’s G-function as \exp\left(-A/\tilde{g}^{2}\right)={\rm G}_{1,0}^{0,1}\left[\tilde{g}^{2}/A\big{|}^{1}_{-}\right]. Then we can apply [44, Eq. (21)] to calculate the infinite integral of product of Meijer’s G-functions involved in $\operatorname{\mathbb{E}}_{\tilde{g}}\!\!\left[\exp\left(-A/\tilde{g}^{2}\right)\right]=\int_{0}^{\infty}\exp\left(-A/\tilde{g}^{2}\right)f_{\tilde{g}}(\tilde{g})d\tilde{g}$ as (12) shown at the top of the next page. Note that the order of G-function in (12) is reduced using [39, Eq. (9.31.1)]. For the ease of notation, hereafter, we denote $\operatorname{\mathbb{E}}_{\tilde{g}}\!\!\left[\exp\left(-A/\tilde{g}^{2}\right)\right]$ by $\mathcal{G}(A)$ for any constant $A$.

Similarly, the second term in (III-A) can be obtained, first, by writing $\Pr(\gamma^{(2)}_{\pi_{2}}>\gamma^{(2)}_{\rm th},\pi_{2})=P(\pi_{2})-\Pr(\gamma^{(2)}_{\pi_{2}}<\gamma^{(2)}_{\rm th},\pi_{2})$. Then using the symmetry of the problem, it can be shown that $\Pr(\gamma^{(2)}_{\pi_{2}}<\gamma^{(2)}_{\rm th},\pi_{2})$ for $\gamma^{(2)}_{\rm th}<1$ and $\gamma^{(2)}_{\rm th}\geq 1$ can be obtained as (III-A) and (14), respectively, shown at the top of the next page, where $\mathcal{G}(\cdot)$ is given in (12), and $J^{(2)}_{\rm th}\triangleq 10^{-s/10}\times\gamma^{(2)}_{\rm th}/(1-\gamma^{(2)}_{\rm th})>0$ is defined for $\gamma^{(2)}_{\rm th}<1$. Also, $J^{(2)}_{{\rm th},1}\triangleq J^{(2)}_{\rm th}({1+10^{s/10}})$ and $J^{(2)}_{{\rm th},2}\triangleq\gamma^{(2)}_{\rm th}+J^{(2)}_{\rm th}({1+\gamma^{(2)}_{\rm th}\times 10^{s/10}})$.

Moreover, the last probability term in (III-A) can be calculated by first writing $\Pr(\gamma^{(1)}_{\pi_{2}}>\gamma^{(1)}_{\rm th},\pi_{2})=\Pr\left(\gamma^{(1)}_{\rm th}\left[{\boldsymbol{\mathcal{I}}_{1}}+C_{\mathcal{D}}/(a_{1}L_{1}P\tilde{g}^{2})\right]<|\tilde{h}_{1}|^{2}\!\!<|\tilde{h}_{2}|^{2}\times 10^{-s/10}\right)$, where $\gamma^{(1)}_{\pi_{2}}$ can be expressed similar to (5). Then using a similar approach to (III-A), the closed-form expression for all values of $\gamma^{(1)}_{\rm th}$ can be expressed as (15) shown at the top of this page. This completes the closed-form characterization of the outage probability of the first NOMA user $\mathcal{U}_{1}$.

Finally, the outage probability of the second NOMA user $\mathcal{U}_{2}$ can be characterized as

[TABLE]

where $\Pr(\gamma^{(2)}_{\pi_{2}}>\gamma^{(2)}_{\rm th},\pi_{2})$ and $\Pr(\gamma^{(1)}_{\pi_{1}}>\gamma^{(1)}_{\rm th},\pi_{1})$ have already been calculated, and, using the symmetry of the problem, $\Pr(\gamma^{(2)}_{\pi_{1}}>\gamma^{(2)}_{\rm th},\pi_{1})$ can be obtained as

[TABLE]

Remark 1: The latter analysis suggests that the outage probability of the second NOMA user can be characterized using the preceding analysis by substituting $-s$ for $s$ and appropriate change of indexing $1\leftrightarrow 2$. This is because the only difference between $\mathcal{U}_{1}$ and $\mathcal{U}_{2}$ is that the user with a lower average gain is labeled as the second user, i.e., $a_{2}L_{2}=a_{1}L_{1}\times 10^{-s/10}$.

III-B Sum-Rate Outage Analysis

For the uplink NOMA transmission with logarithmic functions for the rates as $R^{(i)}_{\pi_{j}}=\log_{2}(1+\gamma^{(i)}_{\pi_{j}})$, $i,j\in\{1,2\}$, it can be verified that the sum of the NOMA users, regardless of their decoding order, can be expressed as

[TABLE]

Denoting the fractional term of the logarithm argument in (III-B) by $\gamma_{\Sigma}$, the sum-rate outage probability defined as $P_{\rm out}^{\Sigma}=\Pr(\gamma_{\Sigma}<\gamma^{\Sigma}_{\rm th})$, where $\gamma^{\Sigma}_{\rm th}=2^{R^{\Sigma}_{\rm th}}-1$ is the threshold equivalent SINR to achieve the desired sum rate of $R^{\Sigma}_{\rm th}$, can be expressed as

[TABLE]

Let $\mathcal{B}$ represent the event $\{|\tilde{h}_{2}|^{2}<\gamma^{\Sigma}_{\rm th}\!\times\!10^{s/10}[\boldsymbol{\mathcal{I}}_{1}\!+\!C_{\mathcal{D}}\!/\!(a_{1}\!L_{1}\!P\tilde{g}^{2})]\}$, and $\mathcal{O_{S}}\triangleq\{|\tilde{h}_{1}|^{2}\!<\!\gamma^{\Sigma}_{\rm th}\!\left[\boldsymbol{\mathcal{I}}_{1}\!\!+\!C_{\mathcal{D}}\!/\!(a_{1}\!L_{1}\!P\tilde{g}^{2})\right]\!\!-\!|\tilde{h}_{2}|^{2}\!\times\!10^{-s/10}\}$ denote the sum-rate outage event defined in (19). Clearly, $\Pr(\mathcal{O_{S}},\mathcal{B}^{c})=0$ where $\mathcal{B}^{c}$ is the complementary event of $\mathcal{B}$. Therefore, using the law of total probability, $P_{\rm out}^{\Sigma}$ can be expressed as $P_{\rm out}^{\Sigma}=\Pr(\mathcal{O_{S}},\mathcal{B})$, which is calculated in a closed form as (III-B) at the top of the next page.

Remark 2: In the special case of the absence of multiuser interference (except the NOMA users themselves), one can obtain the outage probability closed-form expressions by substituting $L^{\prime}_{k}p^{\prime}_{k}=0$, $\forall k=1,2,...,K$, which replaces all the product terms of the form $\prod_{k=1}^{K}[\cdot]$ by $1$ and summations of the form $\sum_{k=1}^{K}[\cdot]$ by [math].

III-C RF Backhaul Analysis

In this part, for the sake of completeness and comparison, we extend our preceding analysis to the case of conventional RF-backhauled systems. In particular, we investigate the performance of dual-hop uplink NOMA where the $\mathcal{R-D}$ backhaul link forwards the amplified version of $y_{\mathcal{R}}$ in (1) by a gain $G_{b}$ through an RF link with the path-loss gain $L_{b}$ and fading coefficient $\tilde{h}_{b}$, i.e., with the composite channel gain $h_{b}=\sqrt{L_{b}}\tilde{h}_{b}$. For the backhaul link, a line-of-sight (LOS) path should be available from the relay to destination for the applicability of the directive FSO link. Therefore, for the $\mathcal{R-D}$ RF backhaul link we assume Rician fading, with the shape parameter $\Omega$ and scale parameter $\Psi$, which takes into account the effects of both LOS and scattered paths. The scale parameter $\Psi$ is the total average power of fading and hence $\Psi=1$, and $\Omega$ represents the ratio of the power contributions by the LOS path to the remaining scattered paths. In this case, $\kappa_{b}\triangleq|\tilde{h}_{b}|^{2}$ can be characterized according to a non-central chi-square distribution given by

[TABLE]

where $I_{0}(\cdot)$ is the zeroth-order modified Bessel function of the first kind [45].

Given that the second hop is employing the RF band, it is practically possible for the destination receiver to be exposed to the presence of some external interfering users. In what follows we characterize the outage probability of the uplink NOMA transmission when the destination, in addition to the directive signal from the $\mathcal{R-D}$ Rician channel, receives the superposition of the signals of $K_{\rm d}$ interfering users $\mathcal{I}^{\prime\prime}_{j}$, $j=1,2,...,K_{\rm d}$, each with the transmit power $p^{\prime\prime}_{j}$, path-loss gain $L^{\prime\prime}_{j}$, and i.i.d. Rayleigh fading coefficient $\tilde{h}^{\prime\prime}_{j}$. In this case, the received signal can be expressed as $y^{\rm RF}_{\mathcal{D}}=\sqrt{L_{b}}\tilde{h}_{b}G_{b}y_{\mathcal{R}}+\sum_{j=1}^{K_{\rm d}}x^{\prime\prime}_{j}\tilde{h}^{\prime\prime}_{j}\sqrt{L^{\prime\prime}_{j}p^{\prime\prime}_{j}}+n^{\rm RF}_{\mathcal{D}}$, where $x^{\prime\prime}_{j}$ is the transmitted signal by the $j$-th interfering users and $n^{\rm RF}_{\mathcal{D}}$ is the RF noise at the destination with mean zero and variance $N_{0}$. Therefore, it is easy to verify that all the SINR equations (e.g., (4) and (5)) remain the same except that $\tilde{g}^{2}$ and $C_{\mathcal{D}}$ will be replaced by $\kappa_{b}\triangleq|\tilde{h}_{b}|^{2}$ and $(N_{0}+\boldsymbol{\mathcal{I}}^{\prime\prime})/(L_{b}G_{b}^{2})=C^{\rm RF}_{\mathcal{D}}+\boldsymbol{\mathcal{I}}^{\prime\prime}/(L_{b}G_{b}^{2})$, respectively, where $\boldsymbol{\mathcal{I}}^{\prime\prime}\triangleq\sum_{j=1}^{K_{\rm d}}{L^{\prime\prime}_{j}p^{\prime\prime}_{j}}|\tilde{h}^{\prime\prime}_{j}|^{2}$ is the total external interference to the destination and $C^{\rm RF}_{\mathcal{D}}\triangleq N_{0}/(L_{b}G_{b}^{2})$. Therefore, in order to characterize the outage probability formulas, instead of calculating $\operatorname{\mathbb{E}}_{\tilde{g}}[\exp(-BC_{\mathcal{D}}/\tilde{g}^{2})]$ or equivalently ${\mathcal{G}}(BC_{\mathcal{D}})$, we need to evaluate $\operatorname{\mathbb{E}}_{\kappa_{b},\boldsymbol{\mathcal{I}}^{\prime\prime}}[\exp(-B(\boldsymbol{\mathcal{I}}^{\prime\prime}+N_{0})/(L_{b}G_{b}^{2}\kappa_{b}))]$. Note that the constant $B$ represents the factor of $-C_{\mathcal{D}}/\tilde{g}^{2}$ or $C_{\mathcal{D}}$ in the arguments of $\operatorname{\mathbb{E}}_{\tilde{g}}[\exp(\cdot)]$ or ${\mathcal{G}}(\cdot)$, respectively, in the previously-derived formulas. The required expectation can be characterized as

[TABLE]

where $(a)$ is obtained by the partial fraction decomposition of the involved rational function assuming $L^{\prime\prime}_{j}p^{\prime\prime}_{j}\neq L^{\prime\prime}_{l}p^{\prime\prime}_{l}\neq 0$, $\forall j\neq l$, which is a reasonable assumption given that $L^{\prime\prime}_{j}p^{\prime\prime}_{j}$’s are continuous quantities222From the stochastic geometry perspective, the interfering users are located in stochastic locations leading to path-losses that are continuous RVs. Hence, the probability of having exactly equal $L^{\prime\prime}_{k}p^{\prime\prime}_{k}$’s in a real-world scenario is zero. Accordingly, our derived closed-form expressions hold with probability one in real-world networks. However, if, hypothetically, some users have equal distances and hence equal $L^{\prime\prime}_{k}p^{\prime\prime}_{k}$’s, then some of our derived formulas cannot further be simplified, e.g., we cannot write the equality in step $(a)$ of (III-C).. Moreover, step $(b)$ is obtained according to Appendix A, where $\mathsf{G}_{.,.}^{.,.}$ denotes the upper incomplete Meijer’s G-function which is a special case of the generalized upper incomplete Fox’s H-function $\mathcal{H}_{.,.}^{.,.}$ introduced in [46].333A simple implementation of $\mathcal{H}_{.,.}^{.,.}$ in Mathematica is presented in [46, Appendix B].

Remark 3: As a special case, in the absence of external RF interference to the destination, the received signal can be expressed as $y^{\rm RF}_{\mathcal{D}}=\sqrt{L_{b}}\tilde{h}_{b}G_{b}y_{\mathcal{R}}+n^{\rm RF}_{\mathcal{D}}$; hence, $\boldsymbol{\mathcal{I}}^{\prime\prime}=0$. Therefore, all the preceding results in Sections III-A and III-B are valid for the RF-backhauled NOMA relaying system when $C_{\mathcal{D}}^{\rm RF}$ and $\operatorname{\mathbb{E}}_{\kappa_{b}}[\exp(-A/\kappa_{b})]$ are substituted for $C_{\mathcal{D}}$ and $\operatorname{\mathbb{E}}_{\tilde{g}}[\exp(-A/\tilde{g}^{2})]$, respectively, where $A$ is a constant, i.e., the argument of $\mathcal{G}(\cdot)$ in the formulas derived earlier. A particularly interesting form for $\operatorname{\mathbb{E}}_{\kappa_{b}}[\exp(-A/\kappa_{b})]$ can be obtained by applying [39, Eq. (8.447.1)] to expand $I_{0}\!\left(2\sqrt{{\Omega(1\!+\!\Omega)\kappa_{b}}}\right)$ in the form of an infinite series and then applying [39, Eq. (3.471.9)] to express $\operatorname{\mathbb{E}}_{\kappa_{b}}[\exp(-A/\kappa_{b})]$ as

[TABLE]

where $K_{n}(\cdot)$ is the $n$-th order modified Bessel function of the second kind. As verified in [47], the summation in (III-C) can effectively be calculated using only few terms, resulting in closed-form expressions for the outage probabilities.

Finally, it is worth mentioning that given any threshold/target data rate, the individual data rates of $R^{(i)}_{\rm th}(1-P_{\rm out}^{(i)})$, $i=1,2$, and sum rate of $R^{\Sigma}_{\rm th}(1-P_{\rm out}^{\Sigma})$ are achievable where the outage probabilities are calculated with respect to the threshold rates $R^{(1)}_{\rm th}$, $R^{(2)}_{\rm th}$, and $R^{\Sigma}_{\rm th}$ according to the analysis in this section. In addition to this characterization of achievable rate regions, we are also interested to determine the average data rates of the users, known as their ergodic capacity, which is obtained using the assumption that the users’ threshold/target data rates are adjusted by their channel conditions, i.e., $\gamma^{(i)}_{\rm th}=\gamma^{(i)}$ and $\gamma^{\Sigma}_{\rm th}=\gamma_{\Sigma}$, meaning that they can always decode their signals without outage. Such analysis will be proceeded in the next section in a same sequential order as this section.

IV Ergodic Capacity Analysis

IV-A Average Individual Rates

Given that the instantaneous individual rate $R_{i}$ of the $i$-th NOMA user, $i=1,2$, is related to its instantaneous SINR $\gamma^{(i)}$ as $R_{i}=\log_{2}(1+\gamma^{(i)})$, and that $\gamma^{(i)}$ can take two different forms of $\gamma^{(i)}_{\pi_{1}}$ and $\gamma^{(i)}_{\pi_{2}}$ each with the probabilities $P(\pi_{1})$ and $P(\pi_{2})$, respectively, the average individual rate $\overline{R}_{i}$ of the $i$-th NOMA user can be characterized as

[TABLE]

In this subsection, we first characterize the first user’s average rate $\overline{R}_{1}$ (equivalently, calculate the two expectation terms in (IV-A) for $i=1$), and then relate $\overline{R}_{2}$ to $\overline{R}_{1}$.

The first expectation term in (IV-A), for the first user, can be evaluated as (IV-A) shown at the top of the next page, where step $(a)$ is written using (4), and step $(b)$ follows from the cumulative density function (CDF) of the exponential distribution and applying the part-by-part integration equality

[TABLE]

with $c$ being a constant, and $f_{Y}(\cdot)$ and $F_{Y}(\cdot)$ being the probability density function (PDF) and CDF of the RV $Y$, respectively.

The first expectation term of the step $(b)$ in (IV-A), denoted by $V_{1}$, can be calculated as

[TABLE]

where step $(a)$ is obtained using the integral formula [39, Eq. (4.337.1)] with ${\rm Ei}(\cdot)$ being the exponential integral function defined as [39, Eq. (8.211.1)] for negative arguments, and $v_{1,i_{1}}=(1+10^{s/10})/i_{1}$, $i_{1}=1,2$. As it is proven in Appendix C, the expectation term involved in step $(a)$ of (IV-A) can be calculated as (IV-A) shown at the top of the next page where the coefficients for $k=2,3,...,K$ can, recursively, be obtained as

[TABLE]

with the initial values $\beta^{(1)}_{\alpha_{1}}=-\beta^{(1)}_{v_{i_{1}}}=a_{1}L_{1}P/(L^{\prime}_{1}p^{\prime}_{1}v_{1,i_{1}}\!-a_{1}L_{1}P)$. We should remark that this result is obtained with the reasonable assumption $L^{\prime}_{i}p^{\prime}_{i}\neq L^{\prime}_{j}p^{\prime}_{j}\neq 0$, $\forall i\neq j$, and also $a_{1}L_{1}P\neq L^{\prime}_{l}p^{\prime}_{l}v_{1,i_{1}}$, $\forall l=1,2,...,K$.

The second expectation term in (IV-A), denoted by $V_{2}$, can be expressed as

[TABLE]

where step $(a)$ follows by calculating the integration over $|\tilde{h}_{1}|^{2}$ using [39, Eq. (3.352.2)], and step $(b)$ is derived using (42) in Appendix B by defining $v_{2,1}=(1+10^{s/10})/2$ and $v_{2,2}=1$. The expectation term involved in step $(b)$ can be calculated according to (IV-A) and (IV-A) by substituting $v_{2,i_{2}}$ for $v_{1,i_{1}}$.

On the other hand, the second expectation term in (IV-A) can be characterized as (IV-A), shown at the top of the next page, where step $(a)$ is from the definition of $\gamma^{(1)}_{\pi_{2}}$ as the dual of (5). Step $(b)$ is obtained by first evaluating the finite integral over $x$ in step $(a)$ using the part-by-part integration equality in (IV-A) and then applying [39, Eq. (3.352.1)]. Moreover, step $(c)$ follows from applying [39, Eq. (4.337.2)] to evaluate the expectation over $|\tilde{h}_{2}|^{2}$. Finally, step $(d)$ is written after calculating the integral involved in step $(c)$ using (42) in Appendix B. Note that ${\rm Ei}(\cdot)$ is a negative quantity for the negative arguments (see, e.g., [48, Table 1]); therefore, the expression obtained in step $(d)$ of (IV-A) takes positive values. Similarly, the expectation term involved in step $(c)$ can be calculated according to (IV-A) and (IV-A) by substituting $1+10^{s/10}$ for $v_{1,i_{1}}$. Finally, the expectation over $\tilde{g}$ in (IV-A) can be performed through a one-dimensional integral weighted by $f_{\tilde{g}}(\tilde{g})$ given in (2) which, to the best of our knowledge, cannot be calculated in a closed form.

The above analysis complete the characterization of the first NOMA user’s average individual rate. A similar characterization can be obtained for $\mathcal{U}_{2}$ by substituting $-s$ for $s$ and appropriate change of indexing $1\leftrightarrow 2$.

IV-B Average Sum Rate

As discussed before, the sum rate of NOMA users can always be expressed as (III-B), regardless of their decoding order. Therefore, the average sum rate of NOMA users can be calculated as (IV-B) shown at the top of the next page, where step $(a)$ follows by the definition of instantaneous sum rate in (III-B), applying [39, Eq. (4.337.2)], and defining ${\rm eEi}(t)\triangleq{\rm e}^{t}{\rm Ei}(-t)$, $\forall t>0$. Moreover, step $(b)$ is derived by applying [39, Eq. (4.337.2)] and Eq. (42) in Appendix B, and then defining $v_{s,1}=10^{s/10}$, $v_{s,2}=1$, $A_{s,1}=(10^{s/10}-1)^{-1}$, and $A_{s,2}=-(1-10^{-s/10})^{-1}$. The expectation term involved in step $(b)$ can be calculated by first applying (IV-A) and (IV-A) with substitution of $v_{s,i_{s}}$ for $v_{1,i_{1}}$, and then taking a one-dimensional integral over $\tilde{g}$ weighted by $f_{\tilde{g}}(\tilde{g})$.

IV-C RF-Backhauled System

In this subsection, we explain how the preceding ergodic capacity analysis can be extended to RF-backhauled systems. As discussed in Section III-C, in the absence of external interference to the destination, the results for the FSO-backhauled system can be extended to that of RF-backhauled one with substituting $C_{\mathcal{D}}^{\rm RF}=N_{0}/(L_{b}G_{b}^{2})$ and $\kappa_{b}=|\tilde{h}_{b}|^{2}$ for $C_{\mathcal{D}}$ and $\tilde{g}^{2}$, respectively. Therefore, the only difference will be the calculation of expectation terms of the form $\operatorname{\mathbb{E}}_{\kappa_{b}}\!\left[\exp\!\left(c^{\prime}_{b}[{\sigma^{2}_{\mathcal{R}}\!+\!C_{\mathcal{D}}^{\rm RF}/\kappa_{b}}]\right)\!{\rm Ei}\!\left(-c^{\prime}_{b}[{\sigma^{2}_{\mathcal{R}}\!+\!C_{\mathcal{D}}^{\rm RF}/\kappa_{b}}]\right)\right]$ instead of $\operatorname{\mathbb{E}}_{\tilde{g}}\!\left[\exp\!\left(c^{\prime}_{b}[{\sigma^{2}_{\mathcal{R}}\!+\!C_{\mathcal{D}}/\tilde{g}^{2}}]\right)\!{\rm Ei}\!\left(-c^{\prime}_{b}[{\sigma^{2}_{\mathcal{R}}\!+\!C_{\mathcal{D}}/\tilde{g}^{2}}]\right)\right]$, in all of the preceding ergodic capacity analysis, where $\kappa_{b}$ is distributed according to (21), and $c^{\prime}_{b}>0$ is a constant, e.g., $c^{\prime}_{b}=v_{1,i_{1}}/(a_{1}L_{1}P)$ in the first expectation term of (IV-A). As we know, such expectations cannot be characterized in closed forms and need to be calculated through one-dimensional integrals weighted by $f_{\kappa_{b}}(\kappa_{b})$.

Moreover, when the destination is subject to the presence of $K_{\rm d}$ interfering users, $C_{\mathcal{D}}/\tilde{g}^{2}$ in all of the previously-derived ergodic capacity formulas should be replaced by $(\boldsymbol{\mathcal{I}}^{\prime\prime}+N_{0})/(L_{b}G_{b}^{2}\kappa_{b})=C_{\mathcal{D}}^{\rm RF}/\kappa_{b}+\boldsymbol{\mathcal{I}}^{\prime\prime}/(L_{b}G_{b}^{2}\kappa_{b})$, where $\boldsymbol{\mathcal{I}}^{\prime\prime}\triangleq\sum_{j=1}^{K_{\rm d}}{L^{\prime\prime}_{j}p^{\prime\prime}_{j}}|\tilde{h}^{\prime\prime}_{j}|^{2}$ is the total external interference to the destination. Therefore, in order to characterize the ergodic capacity, instead of calculating $\operatorname{\mathbb{E}}_{\tilde{g}}\!\left[\exp\!\left(c^{\prime}_{b}[{\sigma^{2}_{\mathcal{R}}\!+\!C_{\mathcal{D}}/\tilde{g}^{2}}]\right)\!{\rm Ei}\!\left(-c^{\prime}_{b}[{\sigma^{2}_{\mathcal{R}}\!+\!C_{\mathcal{D}}/\tilde{g}^{2}}]\right)\right]$, we have to evaluate

[TABLE]

By defining $B_{\kappa_{b}}\triangleq\sigma^{2}_{\mathcal{R}}\!+C_{\mathcal{D}}^{\rm RF}/\kappa_{b}$ and $\alpha^{\prime}_{j}(\kappa_{b})\triangleq L^{\prime\prime}_{j}p^{\prime\prime}_{j}/(L_{b}G_{b}^{2}\kappa_{b})$, $j=1,2,...,K_{\rm d}$, and then applying a similar approach to Appendix C, one can show that $E_{b}$ in (IV-C) can be calculated as (34) shown at the top of the next page

where the coefficients for $k=2,3,...,K_{\rm d}$ are, recursively, defined as

[TABLE]

with the initial values $\beta^{(1)}_{\alpha^{\prime}_{1}}=-\beta^{(1)}_{c^{\prime}_{b}}=L_{b}G_{b}^{2}\kappa_{b}/(L^{\prime\prime}_{1}p^{\prime\prime}_{1}c^{\prime}_{b}-L_{b}G_{b}^{2}\kappa_{b})$. Finally, the ergodic capacity formulas for the RF-backhauled system can be obtained after one-dimensional integrals over $\kappa_{b}$ to perform the expectations of the form expressed in (34). We should again emphasize that this result is valid under the reasonable assumption $L^{\prime\prime}_{i}p^{\prime\prime}_{i}\neq L^{\prime\prime}_{j}p^{\prime\prime}_{j}\neq 0$, $\forall i\neq j$, and also $L_{b}G_{b}^{2}\kappa_{b}\neq L^{\prime\prime}_{l}p^{\prime\prime}_{l}c^{\prime}_{b}$, $\forall l=1,2,...,K_{\rm d}$.

Remark 4: Our characterization of the average individual- and sum-rate formulas for both FSO- and RF-backhauled systems are up to only one-dimensional integrals over the fading coefficient of the backhaul link. That is equivalent to say that the ergodic capacity closed-form expressions are obtained for the single-hop uplink NOMA subject to some exterior multiuser interference, or better to say, for the considered dual-hop system model given each realization of the backhaul fading coefficient.

V Numerical Results

In this section, we present the numerical results to evaluate the performance of uplink NOMA over mixed RF-FSO and dual-hop RF/RF systems, and corroborate the correctness of the derived outage probability and ergodic capacity formulas. Some of the parameters considered for simulations (unless explicitly specified) are listed in Table I. For the multiuser interference, we consider the product of $L^{\prime}_{k}p^{\prime}_{k}$, $k=1,2,...,K$, to be the $k$-th element of the vector $K_{\mathcal{I},R}P_{0}L_{2}{\boldsymbol{u}^{\prime}_{10}}$ where $P_{0}=1$ $\mathrm{m}\mathrm{W}$, $K_{\mathcal{I,R}}\geq 0$ is a constant to define the upper bound of the received power from each interfering user to the relay as a factor of $P_{0}L_{2}$, and ${\boldsymbol{u}^{\prime}_{10}}=(0.6957,\!0.6279,\!0.4504,\!0.4736,\!0.9497,\!0.0835,0.2798,0.4470,\\ 0.5876,0.8776)$ is a length-10 vector of uniformly generated numbers over the interval $(0,1)$. Similarly, we consider the external interference to the RF receiver of the destination in the case of the RF-backhauled system to be of the form $K_{\mathcal{I},D}P_{0}L_{2}{\boldsymbol{u}^{\prime\prime}_{10}}$ with the constant $K_{\mathcal{I},D}\geq 0$ defining the upper bound of the received power from each interfering user to the destination as a factor of $P_{0}L_{2}$ and ${\boldsymbol{u}^{\prime\prime}_{10}}=(0.5259,\!0.9635,\!0.5688,\!0.2584,\!0.2959,\!0.7439,0.9797,0.3491,\\ 0.8371,0.5587)$ being another length-10 vector of uniformly generated numbers over the interval $(0,1)$.

Fig. 2 shows the individual-rate outage performance of the uplink mixed RF-FSO NOMA system for three different values of the power back-off step $s=0$, $10$, and $25$ $\mathrm{d}\mathrm{B}$. We assume that the first and second users are in distances of $d^{\mathrm{RF}}_{1}=100$ $\mathrm{m}$ and $d^{\mathrm{RF}}_{2}=200$ $\mathrm{m}$ from the relay, respectively; this together with the values in Table I implies $L_{1}=9.04\times 10^{-8}$ and $L_{2}=8\times 10^{-9}$. The other adopted parameters are listed in the caption of Fig. 2. For $s=0$ we will have $a_{1}L_{1}=a_{2}L_{2}$; therefore, one should expect lower outage probabilities for the second NOMA user given its lower threshold SINR $\gamma_{\rm th}^{(2)}=0.4$ (equivalently, lower target rate $R_{\rm th}^{(2)}=0.485$ $\mathrm{b}\mathrm{i}\mathrm{t}\mathrm{s}\mathrm{/}\mathrm{s}\mathrm{e}\mathrm{c}\mathrm{/}\mathrm{H}\mathrm{z}$). However, by increasing $s$ a larger fraction of power will be assigned to the first NOMA user, and $\mathcal{U}_{1}$ achieves lower outage probabilities even if it demands a larger SINR threshold. As a consequence, increasing $s$ will decrease the outage probability of $\mathcal{U}_{1}$ and increase the outage probability of $\mathcal{U}_{2}$. Moreover, the match between the analytical results and Monte-Carlo simulations corroborates the correctness of the derived closed-form expressions for the individual-rate outage probabilities.

The comparison between NOMA and orthogonal multiple access (OMA) is also depicted in Fig. 2. We assume, for the OMA transmission, that the total transmission time is equally divided between the two users and each user employs the entire transmission power $P$ during its corresponding time slot. Therefore, denoting the SINR of the $i$-th OMA user by $\gamma^{(i)}_{\rm OMA}$, $i\in\{1,2\}$, each OMA user will have the rate $R^{(i)}_{\rm OMA}=0.5\log_{2}(1+\gamma^{(i)}_{\rm OMA})$. Then it is easy to verify that, in order to achieve the target data rate $R_{\rm th}^{(i)}=\log_{2}(1+\gamma_{\rm th}^{(i)})$, each $i$-th OMA user has to satisfy the threshold SINR of $\gamma^{(i)}_{{\rm th},\rm OMA}=(1+\gamma_{\rm th}^{(i)})^{2}-1$. It is observed that NOMA transmission is in favor of the first user except for very small values of $s$ while the second user experiences an opposite situation. In particular for $s=10$ $\mathrm{d}\mathrm{B}$ where $a_{1}L_{1}=10a_{2}L_{2}$, i.e., when the roughly an order of magnitude difference between the channel gains $L_{1}$ and $L_{2}$ is reflected in the NOMA power allocation policy, both of the NOMA users (especially the second user) achieve better performance than the OMA transmission. This suggests the possible existence of some optimal values for $s$ (see also Fig. 5) given some specific figures of merit.

Fig. 3 illustrates the individual-rate outage performance of the mixed RF-FSO NOMA system for $s=5$ $\mathrm{d}\mathrm{B}$ and different values of the threshold SINRs. As expected, outage performance degrades with increasing the threshold SINRs. More importantly, the induced interference between NOMA users due to the non-orthogonal operation limits the outage performance for the large values of the threshold SINRs and prevents achieving small enough outage probabilities even for the large values of the transmitted power. As a result, the system outage performance saturates while the saturation limits are higher for the larger values of the threshold SINRs. Note that the channel is stochastic and the decoding order of the NOMA users changes dynamically. Therefore, the first-decoded user dynamically changes between $\mathcal{U}_{1}$ and $\mathcal{U}_{2}$; hence, both NOMA users experience interference and cannot meet large values of the threshold SINRs even if $P$ is increased.

Sum-rate outage probability results of the mixed RF-FSO NOMA system are depicted in Fig. 4 for $\gamma^{\Sigma}_{\rm th}=1.2$, three different power back-off steps $s=0$, $10$, and $100$ $\mathrm{d}\mathrm{B}$, and two different distances of the second user $d^{\mathrm{RF}}_{2}=200$ and $400$ $\mathrm{m}$. The other parameters are the same as Fig. 2. For the sake of comparison, the simulation results are also provided for the OMA sum-rate outage probabilities. Similar to Fig. 2, we assume that the transmission time is equally divided between the two users while each user transmits with the full power $P$ during its transmission time. Then it is easy to verify that the sum-rate outage probability of the OMA transmission in satisfying the target rate $R^{\Sigma}_{\rm th}=\log_{2}(1+\gamma^{\Sigma}_{\rm th})$ can be obtained as

[TABLE]

For the case of $d^{\mathrm{RF}}_{2}=2d^{\mathrm{RF}}_{1}=200$ $\mathrm{m}$, which corresponds to $L_{1}=9.04\times 10^{-8}$ and $L_{2}=8\times 10^{-9}$, one can see that by increasing $s$ the sum-rate outage probability decreases first and then increases, at least for higher transmit powers. Specifically, for $s=10$ $\mathrm{d}\mathrm{B}$ where $a_{1}L_{1}=10a_{2}L_{2}$, i.e., when the roughly an order of magnitude difference between the channel gains $L_{1}$ and $L_{2}$ is reflected in the NOMA power allocation policy, NOMA outperforms OMA in terms of the sum-rate outage probability. One can easily verify that in the special case of second user having a very poor channel condition compared to the first user, NOMA, with assigning the whole power to the first user (equivalently, $s\to\infty$), can achieve twice the OMA rate. This can be observed also from the case of $d^{\mathrm{RF}}_{2}=4d^{\mathrm{RF}}_{1}=400$ $\mathrm{m}$, corresponding to $L_{1}=9.04\times 10^{-8}$ and $L_{2}=7.06\times 10^{-10}$, with the choice of $s=100$ $\mathrm{d}\mathrm{B}$ (at least at reasonably low transmit powers where $R^{(2)}_{\rm OMA}=0.5\log_{2}(1+\gamma^{(2)}_{\rm OMA})\approx 0$).

To better observe the trend of sum-rate outage probability as a function of $s$, we have characterized in Fig. 5 the sum-rate outage performance of the system for $P=25$ $\mathrm{d}\mathrm{B}\mathrm{m}$, and different values of $d^{\mathrm{RF}}_{2}$ and $\gamma^{\Sigma}_{\rm th}$. It is observed that the outage performance decreases first and then increases with the increase of $s$. Therefore, given any set of system parameters, there is a power back-off step $s^{*}$ minimizing the sum-rate outage probability. However, such a $s^{*}$ is not necessarily the best operation point as such an operation region may depend to the individual outage probabilities and achievable rates and not only to the sum-rate outage probability. By comparing the four plots corresponding to $\gamma_{\rm th}^{\Sigma}=1.2$ we can observe that increasing $d^{\mathrm{RF}}_{2}$, i.e., differentiating more between the channel qualities of the NOMA users, shifts $s^{*}$ to the larger values. Moreover, the outage probability is higher for the larger values of $d^{\mathrm{RF}}_{2}$ before the minimum of $s^{*}$’s and is lower after their maximum. Furthermore, we have observed that increasing $\gamma_{\rm th}^{\Sigma}$, while $d^{\mathrm{RF}}_{2}$ is kept fixed, increases the outage probability and does not change the outage-minimizing $s$. Finally, it can be observed that as $s$ increases, we have $a_{1}\to 1$ and $a_{2}\to 0$. Therefore, further increasing $s$, beyond a certain large threshold, does not noticeably change the power allocation coefficients. Consequently, the outage probability saturates while increasing $s$.

The individual- and sum-rate outage probability results of the dual-hop RF/RF NOMA system with $d_{\mathcal{RD}}=500$ $\mathrm{m}$ are characterized in Fig. 6 for three different values of $K_{\mathcal{I,D}}=0$, $0.1$ and $1$. As observed, even small values of $K_{\mathcal{I,D}}$ result in significant performance degradation, mainly due to the weak desired signal at the destination. This figure, in addition to confirming the accuracy of the results derived in Section III-C, highlights the superiority of FSO backhauling in terms of outage performance (compare the case of $K_{\mathcal{I,D}}=0$ in Fig. 6 with the corresponding plots in Figs. 2 and 4, and note that the length of backhaul link here is much less than that of FSO backhaul), especially for longer backhaul ranges, noticing that FSO links with even lower SINRs may be preferred given their much more available bandwidths compared to RF backhaul links. We should emphasize that increasing the number of interfering users increases the outage probability in a similar fashion. However, the results are not included in this paper due to the space limitation.

Average individual and sum rates of the mixed RF-FSO NOMA system are characterized in Fig. 7 for $d^{\mathrm{RF}}_{1}=100$ $\mathrm{m}$, $d^{\mathrm{RF}}_{2}=200$ $\mathrm{m}$, $d_{\mathcal{RD}}=1200$ $\mathrm{m}$, and three values of the power back-off step $s=0$, $10$, and $30$ $\mathrm{d}\mathrm{B}$. It is observed that, at relatively low transmit powers, increasing $s$ (within the range of $s$ in this figure) increases the rate of the first NOMA user and decreases the rate of the second user while the average sum rate is also increased. This is because $a_{1}L_{1}=a_{2}L_{2}\times 10^{s/10}$ and $s\geq 0$; therefore, more power is assigned to the first NOMA user for the larger values of $s$ which increases the first user’s rate. On the other hand, for the larger values of $s$ with higher probabilities the decoding order is $\pi_{1}$, i.e., the first user is decoded first in the presence of interference signal from the second NOMA user; hence, by increasing $P$ the interference power to the first user increases limiting its performance and saturating its average rate at larger $P$’s. This is while the second NOMA user is usually decoded after removing the interference signal of the first user, especially for larger values of $s$; hence, the average rate of the second NOMA user can increase without saturation. Given the increasing trend of the rate of the second NOMA user and the saturating trend of the first NOMA user, the average rate curves of these two users intersect at some power $P$ such that the first user has a higher rate before that crossing point and the second user outperforms after that power. Also as $s$ increases, this crossing point shifts to the right meaning that the first NOMA user has higher data rates over wider ranges of the transmit power for the larger values of $s$. We have also included the simulation results of the OMA average rates. The simulation results confirm the superiority of NOMA over OMA in terms of individual and sum rates (see, e.g., the results for $s=10$ $\mathrm{d}\mathrm{B}$).

The impact of the FSO link parameters on the average sum rate of the NOMA and OMA mixed RF-FSO systems is investigated in Fig. 8 for $s=25$, $d^{\mathrm{RF}}_{1}=100$ $\mathrm{m}$, $d^{\mathrm{RF}}_{2}=200$ $\mathrm{m}$, $d_{\mathcal{RD}}=1200$ $\mathrm{m}$, and different values of $\kappa$, $\alpha$, $\beta$, and $\zeta$. It is observed that the average sum rate decreases by increasing $\kappa$ (equivalent to having larger path-losses) and decreasing $\alpha$, $\beta$, and $\zeta$ (equivalent to stronger atmospheric turbulence). As a result, in the extreme atmospheric conditions the performance of FSO backhauling might be inferior to that of RF backhauling; this necessitates adaptive switching between RF and FSO infrastructure to simultaneously take the advantage of the potentials of both FSO (e.g., higher data rates and longer ranges) and RF (e.g., resilience to atmospheric adverse conditions) backhauling. Furthermore, note that for all of the considered scenarios NOMA average sum rate outperforms that of the OMA transmission.

Average individual and sum rates of the dual-hop RF/RF system are characterized in Fig. 9 for $s=8$ $\mathrm{d}\mathrm{B}$, $d^{\mathrm{RF}}_{1}=100$ $\mathrm{m}$, $d^{\mathrm{RF}}_{2}=200$ $\mathrm{m}$, and two different values of $d_{\mathcal{RD}}=500$ and $1200$ $\mathrm{m}$. It is observed that as the backhaul link length increases the point for which the individual rates intersect (see, e.g., Fig. 7) moves to the right hand side, i.e., the first NOMA user will have higher average rates over a wider range of the transmit powers. This figure, in addition to confirming the accuracy of the results derived in Section IV-C, further highlights the superiority of FSO backhauling in terms of average achievable rates, especially for longer backhaul ranges, given much more available bandwidth of optical links.

VI Conclusions and Future Directions

We studied the performance of uplink NOMA over mixed RF-FSO systems where the NOMA users employ RF signals for concurrent transmission to an intermediate relay which forwards the amplified version of the received signal in the presence of multiuser interference to the destination through an ultra high-speed FSO link. We adopted dynamic-order decoding to determine the priority of the users at the destination based on their instantaneous CSI. The analysis of the paper carried out in a general model of the channel and system, such as including different aspects of the RF and FSO channels and taking into account the presence of multiuser interference, e.g., due the co-channel interference from some nearby users aiming to communicate to some other relays or destinations. We also extended the results to the case of RF-backhauled systems, i.e., the conventional dual-hop RF/RF systems in the presence of multiuser interference. In terms of the outage probability performance, we derived the closed-form expressions for both mixed RF-FSO and dual-hop RF/RF NOMA systems. Moreover, we characterized the average individual and sum rates of both FSO- and RF-backhauled systems up to only one-dimensional integrals over the fading coefficient of the backhaul link. We presented extensive numerical results confirming the accuracy of the derived outage probability and average rate formulas and characterizing the impact of various design and channel parameters on the system performance.

Our results revealed the superiority of FSO backhauling compared to RF backhauling in terms of outage probability and ergodic capacity. This advantage can further be highlighted, as a future work, by analyzing the throughput of the dual-hop RF-FSO and RF/RF NOMA systems (e.g., when the users of the access links have enough data for transmission and the backhaul link is imposing a bottleneck) by taking into account the much higher available bandwidth of optical links. Moreover, given the sensitivity of FSO links to the atmospheric conditions (see, e.g., Fig. 8), it is inevitable to design hybrid RF/FSO backhaul links to simultaneously take the advantages of both FSO and RF systems; this necessitates further studies to characterize the performance of relay-assisted NOMA systems when the backhaul link dynamically switches between FSO and RF systems given properly-defined switching policies. More importantly, the results of this paper provide novel expressions for the outage probabilities and ergodic capacities of dual-hop RF-FSO and RF/RF systems subject to some external interference to the RF receivers. This is particularly relevant to mmWave NOMA systems which recently have received extensive attention. In fact, in mmWave NOMA systems, inter-beam interference due to the side-lobes of nearby mmWave beams can adversely affect the power-domain NOMA users grouped over a given mmWave beam. As a result, the analysis in this paper can significantly pave the way toward future research on mmWave NOMA systems.

Appendix A Closed-Form Expression for Eq. (III-C)

In order to derive the closed form of Eq. (III-C), we need to calculate integrals of the form

[TABLE]

for some positive constants $a,b,c,d$, and integer $n$.444The use of symbols in this appendix should not be confused with the rest of the paper as we only reuse some symbols during the proof steps. Using [39, Eq. (8.447.1)] to expand $I_{0}\!\left(\sqrt{dx}\right)$ as an infinite series, we can rewrite (37) as

[TABLE]

Hence, the problem boils down to the closed-form calculation of $\mathcal{J}_{2}$ defined in (38). To do so, we first use [49, Eq. (07.34.03.0228.01)] together with [39, Eq. (9.301)] to write $\exp(-b/x)$ in the integral form

[TABLE]

where $\mathcal{C}$ is a complex contour of integration ensuring the convergence of the above Mellin-Barnes type integral, e.g., a vertical line starting from the point $-\epsilon-j\infty$ and terminating at the point $-\epsilon+j\infty$ for any $\epsilon>0$.555Based on [50, Sec. 1.1], this is a proper choice for $\mathcal{C}$ that separates all the poles (i.e., all nonnegative integers) to the right, as required for the convergence of (39). Now, by substituting (39) in $\mathcal{J}_{2}$ we have

[TABLE]

where $(a)$ is calculated using [39, Eq. (3.383.10)] when noting that $\operatorname{Re}(n+k+1-s)=n+k+\epsilon+1>0$, and $(b)$ is obtained using the definition of the generalized upper incomplete Fox’s H-function $\mathcal{H}_{.,.}^{.,.}$ [46]. Note that the upper incomplete gamma function $\Gamma(s\!-\!n\!-\!k,\!ac)$ does not have any poles, and the poles of the other two gamma functions in $(a)$, i.e., $\Gamma(-s)$ and $\Gamma(n+k+1-s)$ are at $l$ and $k+n+1+l$, respectively, for any nonnegative integer $l$. Therefore, all the poles of $\Gamma(-s)$ and $\Gamma(n+k+1-s)$ lie on the right side of the vertical line considered above for $\mathcal{C}$; thus, according to [46], the aforementioned vertical line for the contour $\mathcal{C}$ also guarantees the convergence of the Mellin-Barnes type integral in $(a)$. We should further mention that step $(c)$ is obtained using [46, Eq. (A.9)], and step $(d)$ is from the definition of the upper incomplete Meijer’s G-function $\mathsf{G}_{.,.}^{.,.}$, which is a special case of $\mathcal{H}_{.,.}^{.,.}$ obtained when the second element of all triples corresponding to the coefficients of $s$ or $-s$ in the arguments of (incomplete) gamma functions are equal to one. Finally, by inserting $n=K_{\rm d}$, $a=1+\Omega$, $b=BC_{\mathcal{D}}^{\rm RF}$, $c=BL^{\prime\prime}_{l}p^{\prime\prime}_{l}/(L_{b}G_{b}^{2})$, and $d=4\Omega(1\!+\!\Omega)$ in (A) and then (38), we get the closed-form expression in (III-C).

Appendix B Useful Integral Equations over ${\rm Ei}(\cdot)$

Lemma 1, Corollary 1, and Corollary 2 can be found in the literature with slight variations. However, independent proofs are included for the sake of completeness.

Lemma 1: For any $c_{1},c_{2}>0$, $a<0$, and $b\in\mathbb{R}$ such that (s.t.) $a+b<0$, we have [49, Eq. (06.35.21.0013.01)]

[TABLE]

where $f(t)|_{c_{1}}^{c_{2}}=f(c_{2})-f(c_{1})$ for the function $f(t)$.

Proof: Note that for $c_{1},c_{2}>0$, $a<0$, and $a+b<0$, the exponential integral functions involved in both right- and left-hand sides (RHS and LHS) of (41) have negative arguments. In that case, based on the definition of ${\rm Ei}(x)$ for $x<0$ [39, Eq. (8.211.1)] we have $\frac{d}{dt}{\rm Ei}(lt)=-\frac{d}{dt}\int_{-lt}^{\infty}\frac{{\rm e}^{-x}}{x}dx={\rm e}^{lt}/t$, $\forall\leavevmode\nobreak\ \!lt\!<\!0$. Then it is easy to verify that the derivative of ${\rm e}^{bt}{\rm Ei}(at)-{\rm Ei}([a+b]t)$ is equal to $b{\rm e}^{bt}{\rm Ei}(at)$; this completes the proof.

Corollary 1: $\forall c_{1}>0$, $a<0$, and $b\in\mathbb{R}$ s.t. $a+b<0$

[TABLE]

Proof: Note first based on [48, Eq. (2)] that ${\rm Ei}(-\infty)=\lim_{x\to-\infty}{\rm Ei}(x)=\lim_{x\to-\infty}{\rm e}^{x}/x=0$. Then by substituting $c_{2}=\infty$ in (41), we have ${\rm Ei}([a+b]c_{2})=0$ since $a+b<0$. Moreover, ${\rm e}^{bc_{2}}{\rm Ei}(ac_{2})|_{c_{2}=\infty}=\lim_{x\to\infty}{\rm e}^{(a+b)x}/ax=0$. Therefore, (42) can be inferred from (41) by inserting $c_{2}=\infty$.

Corollary 2: $\forall c_{2}>0$, $a<0$, and $b\in\mathbb{R}$ s.t. $a+b<0$ [39, Eq. (5.231.2)]

[TABLE]

Proof: By substituting $c_{1}=0$ in (41), we have ${\rm e}^{bc_{1}}{\rm Ei}(ac_{1})-{\rm Ei}([a+b]c_{1})=\lim_{x\to 0}{\rm Ei}(ax)-{\rm Ei}([a+b]x)$. Then using [48, Eq. (1)], $\lim_{x\to 0}{\rm Ei}(x)=\gamma+\ln|x|$, where $\gamma=0.57721$ is the Euler’s constant. Therefore, for $c_{1}=0$ we have ${\rm e}^{bc_{1}}{\rm Ei}(ac_{1})-{\rm Ei}([a+b]c_{1})=-\ln\left(|a+b|/|a|\right)=-\ln\left(1+b/a\right)$. This completes the proof using (41). Furthermore, note that by substituting $c_{2}=\infty$ in (43) and using the discussions along with the proof of Corollary 1 one can obtain

[TABLE]

which is the same result reported in [39, Eq. (6.224.1)].

Appendix C Proof of (IV-A) and (IV-A)

In this appendix, we show how the expectation term involved in $V_{1}$ in (IV-A) can be calculated in a closed form over $\boldsymbol{\mathcal{I}}_{1}$. To begin with, let us, for the ease of notation, define $\alpha_{k}\triangleq L^{\prime}_{k}p^{\prime}_{k}/(a_{1}L_{1}P)$, $k=1,2,...,K$, $B_{\tilde{g}}\triangleq[\sigma^{2}_{\mathcal{R}}+C_{\mathcal{D}}/\tilde{g}^{2}]/(a_{1}L_{1}P)$, $\overline{\boldsymbol{\mathcal{I}}}_{1,l}\triangleq\sum_{k=l+1}^{K}\alpha_{k}|\tilde{h}^{\prime}_{k}|^{2}+B_{\tilde{g}}$, $l=1,2,...,K-1$, and $\overline{\boldsymbol{\mathcal{I}}}_{1,K}\triangleq B_{\tilde{g}}$. In order to take the expectation over $\boldsymbol{\mathcal{I}}_{1}$, we first average over $|\tilde{h}^{\prime}_{1}|^{2}$ given that the other $|\tilde{h}^{\prime}_{k}|^{2}$’s are some constant quantities, then average over $|\tilde{h}^{\prime}_{2}|^{2}$, and so on. The expectation of the term involved in (IV-A) over $|\tilde{h}^{\prime}_{1}|^{2}$ can be calculated as

[TABLE]

where step $(a)$ is derived by letting $f_{|\tilde{h}^{\prime}_{1}|^{2}}(x_{1})={\rm e}^{-x_{1}}$ and then $t_{1}=\alpha_{1}x_{1}+\overline{\boldsymbol{\mathcal{I}}}_{1,1}$, and step $(b)$ follows from (42). Note that in step $(a)$, $c_{1}=\overline{\boldsymbol{\mathcal{I}}}_{1,1}>0$, $a=-v_{1,i_{1}}<0$, and $a+b=-1/\alpha_{1}<0$, allowing to use Corollary 1 in (42). Using a similar procedure, the expectation of $V_{1,1}$ over $|\tilde{h}^{\prime}_{2}|^{2}$ can be expressed as

[TABLE]

For the ease of notation, let us further define ${\rm eEi}(x)\triangleq{\rm e}^{x}{\rm Ei}(-x)$, $\forall x>0$, and denote by $\beta^{(k)}_{\alpha_{i}}$, $i=1,2,...,k$, the coefficient of ${\rm eEi}(\overline{\boldsymbol{\mathcal{I}}}_{1,k}/\alpha_{i})$ in $V_{1,k}$ after taking the $k$-th expectation, i.e., averaging over $|\tilde{h}^{\prime}_{k}|^{2}$, $k=1,2,...,K$, and $\beta^{(k)}_{v_{i_{1}}}$ the coefficient of ${\rm eEi}(v_{1,i_{1}}\overline{\boldsymbol{\mathcal{I}}}_{1,k})$ in $V_{1,k}$. Using this notation, it is easy to observe that $V_{1,3}\triangleq\operatorname{\mathbb{E}}_{|\tilde{h}^{\prime}_{3}|^{2}}[V_{1,2}]$ can be expressed as

[TABLE]

where, based on (C), $\beta^{(2)}_{\alpha_{1}}=-\frac{1}{\alpha_{1}v_{1,i_{1}}\!-\!1}\frac{1}{\alpha_{2}/\alpha_{1}\!-\!1}$, $\beta^{(2)}_{\alpha_{2}}=\frac{1}{\alpha_{1}v_{1,i_{1}}\!-\!1}\frac{1}{\alpha_{2}/\alpha_{1}\!-\!1}-\frac{1}{\alpha_{1}v_{1,i_{1}}\!-\!1}\frac{1}{\alpha_{2}v_{1,i_{1}}\!-\!1}$, and $\beta^{(2)}_{v_{i_{1}}}=\frac{1}{\alpha_{1}v_{1,i_{1}}\!-\!1}\frac{1}{\alpha_{2}v_{1,i_{1}}\!-\!1}$.

Using some inductive arguments and with the help of (C) it can be shown that in the $k$-th expectation step, i.e., calculating $V_{1,k}$, we will have $k+1$ distinct terms of ${\rm eEi}(\overline{\boldsymbol{\mathcal{I}}}_{1,k}/\alpha_{1}),{\rm eEi}(\overline{\boldsymbol{\mathcal{I}}}_{1,k}/\alpha_{2}),...,{\rm eEi}(\overline{\boldsymbol{\mathcal{I}}}_{1,k}/\alpha_{k})$, and ${\rm eEi}(v_{1,i_{1}}\overline{\boldsymbol{\mathcal{I}}}_{1,k})$. Therefore,

[TABLE]

where the coefficients for $k=2,3,...,K$ can, recursively, be obtained as

[TABLE]

with the initial values $\beta^{(1)}_{\alpha_{1}}=-\beta^{(1)}_{v_{i_{1}}}=(\alpha_{1}v_{1,i_{1}}\!-1)^{-1}$. Finally, we need to prove (48) with the coefficients in (C) by recurrence. To this end, the base case is true by inserting $k=1$ in (48) and (C) and then recalling $V_{1,1}$ from (C). Moreover, assuming that the induction hypothesis (i.e., (48) and (C) for $k$) is true, we need to show that the induction step holds for $k+1$. Using the definition of $V_{1,{k+1}}\triangleq\operatorname{\mathbb{E}}_{|\tilde{h}^{\prime}_{k+1}|^{2}}[V_{1,k}]$ we have

[TABLE]

where step $(a)$ follows from the induction hypothesis in (48), step $(b)$ is by recalling $\overline{\boldsymbol{\mathcal{I}}}_{1,k}\triangleq\alpha_{k+1}|\tilde{h}^{\prime}_{k+1}|^{2}+\overline{\boldsymbol{\mathcal{I}}}_{1,k+1}$, and step $(c)$ is obtained using (42) in a similar approach to (C). Hence, using step $(c)$ in (C), we can characterize $V_{1,{k+1}}$ as

[TABLE]

which is equal to $\beta^{(k+1)}_{v_{i_{1}}}{\rm eEi}(v_{1,i_{1}}\overline{\boldsymbol{\mathcal{I}}}_{1,{k+1}})+\sum_{i=1}^{k+1}\beta^{(k+1)}_{\alpha_{i}}{\rm eEi}(\overline{\boldsymbol{\mathcal{I}}}_{1,{k+1}}/\alpha_{i})$ with the coefficients defined similar to (C) by substituting $k+1$ for $k$. This suggests that the induction step holds and completes the proof by recurrence. At the end, by substituting $k=K$ in (48) and (C) and noting that $\overline{\boldsymbol{\mathcal{I}}}_{1,K}=B_{\tilde{g}}$, one can obtain (IV-A) and (IV-A).

Acknowledgment

The authors would like to thank Mr. S. M. Azimi-Abarghouyi for the helpful discussions on an earlier draft version of this paper.