# Outage Performance and Optimal Design of MIMO-NOMA Enhanced Small Cell
Networks With Imperfect Channel-State Information

Zheng Shi  Hong Wang  Yaru Fu  Guanghua Yang  Shaodan Ma  and Xinrong Ye Zheng
Shi and Guanghua Yang are with the School of Intelligent Systems Science and
Engineering, Jinan University, Zhuhai 519070, China, and also with the State
Key Laboratory of Internet of Things for Smart City, University of Macau,
Macao 999078, China (e-mails<EMAIL_ADDRESS>ghyang@jnu.edu.cn).Hong
Wang is with the School of Communication and Information Engineering, Nanjing
University of Posts and Telecommunications, Nanjing 210003, China, and also
with the National Mobile Communications Research Laboratory, Southeast
University, Nanjing 210096, China (e-mail: wanghong@njupt.edu.cn).Yaru Fu is
with the School of Science and Technology, The Open University of Hong Kong,
Hong Kong SAR, China (e-mail: yfu@ouhk.edu.hk).Shaodan Ma is with the State
Key Laboratory of Internet of Things for Smart City and the Department of
Electrical and Computer Engineering, University of Macau, Macao 999078, China
(e-mail: shaodanma@um.edu.mo).Xinrong Ye is with the School of Physics and
Electronic Information, Anhui Normal University, Wuhu 241002, China (e-mail:
shuchong@mail.ahnu.edu.cn).

###### Abstract

This paper focuses on boosting the performance of small cell networks (SCNs)
by integrating multiple-input multiple-output (MIMO) and non-orthogonal
multiple access (NOMA) in consideration of imperfect channel-state information
(CSI). The estimation error and the spatial randomness of base stations (BSs)
are characterized by using Kronecker model and Poisson point process (PPP),
respectively. The outage probabilities of MIMO-NOMA enhanced SCNs are first
derived in closed-form by taking into account two grouping policies, including
random grouping and distance-based grouping. It is revealed that the average
outage probabilities are irrelevant to the intensity of BSs in the
interference-limited regime, while the outage performance deteriorates if the
intensity is sufficiently low. Besides, as the channel uncertainty lessens,
the asymptotic analyses manifest that the target rates must be restricted up
to a bound to achieve an arbitrarily low outage probability in the absence of
the inter-cell interference. Moreover, highly correlated estimation error
ameliorates the outage performance under a low quality of CSI, otherwise it
behaves oppositely. Afterwards, the goodput is maximized by choosing
appropriate precoding matrix, receiver filters and transmission rates. In the
end, the numerical results verify our analysis and corroborate the superiority
of our proposed algorithm.

###### Index Terms:

Imperfect channel-state information (CSI), multiple-input multiple-output
(MIMO), non-orthogonal multiple access (NOMA), small cell networks (SCNs).

##  I Introduction

With the surge of wireless data traffic and number of mobile devices, it has
become increasingly imperative for 5G cellular networks to offer an improved
spectrum efficiency as well as massive connectivity[1, 2]. As per Cisco visual
networking index (VNI)[3], 5G speeds are estimated to be 13 times higher than
the average mobile connection by 2023, and the number of mobile devices is
forecasted to grow from 8.8 billion in 2018 to 13.1 billion by 2023. To
confront these unprecedented challenges, small cell has been recognized as one
of the promising solutions for 5G to deliver ever-increasing network capacity
and fulfill the accommodation of a large number of users[4]. However,
aggressive frequency reuse will incur severe interference for ultra-dense
small cell deployment[5]. To remedy this issue, the fusion of small cell
networks (SCNs) and other key enabling 5G technologies has received enormous
interest[5, 6, 7]. Specifically, the combination of multiple-input multiple-
output (MIMO) and non-orthogonal multiple access (NOMA) is anticipated to
considerably improve spectral efficiency and support more concurrent
connections for interference-infested SCNs[8, 9, 10, 11].

###  I-A Related Works and Motivations

Since MIMO attains additional spatial diversity gain by deploying multiple
antennas and NOMA exploits extra multi-user diversity gain by taking advantage
of the channel discrepancy, it has been vastly reported that MIMO-NOMA
achieves superior performance over its orthogonal counterpart, i.e., MIMO-
orthogonal multiple access (OMA), in terms of both the reliability and the
spectral efficiency [12, 13, 14, 15]. Nevertheless, the application of MIMO-
NOMA to SCNs still remains in its infancy. Notwithstanding, a wide range of
investigations on single-input single-out (SISO)-NOMA enhanced SCNs have been
conducted in the literature [16, 5, 17]. To mention only a few, by assuming
single antenna at both small cell base stations (BSs) and users, the coverage
probability and achievable rate of non-uniform SCNs were examined in [5]. By
considering large-scale networks, it was studied in [16] that the coverage
probability of small cells depends heavily on the prescribed transmit rates
and power sharing coefficients. In addition, a novel energy-efficient
algorithm was proposed for NOMA aided heterogeneous SCNs with energy
harvesting by assuming perfect channel-state information (CSI) in [17].
Unfortunately, prior works seldom touched MIMO-NOMA enhanced SCNs except for
[8, 9]. In particular, the authors in [8] substantiated the superior
performance of MIMO-NOMA over MIMO-OMA in multi-cell systems in terms of the
ergodic sum-rate. Besides, a new interference mitigation and power allocation
scheme was invented for downlink MIMO-NOMA assisted heterogeneous SCNs [9].
However, the existing literature commonly assume perfect CSI at the receivers
without the practical consideration of imperfect CSI in MIMO-NOMA enhanced
SCNs.

Unfortunately, there are only few available methodologies that can be
leveraged to undertake the performance investigations into MIMO-NOMA enhanced
SCNs with imperfect CSI. To be specific, K. He _et al._ developed a novel
compressed sensing-aided MIMO-NOMA scheme by allowing for imperfect CSI in
[18], where an upper bound for the performance of the proposed scheme was
presented by using restricted isometry property. In [19], the worst-case
Gaussian approximation was employed to derive the total sum rate, with which
the significant spectral and energy efficiency gains over the existing multi-
way relay network were demonstrated. A semi-closed expression was provided for
the outage performance evaluation of multiple-input single-output (MISO)-NOMA
in view of imperfect CSI in [20], in which the assumption of correlated CSI
was made. In [21], Gaussian approximation was used to derive the outage
probability of MIMO-NOMA systems under imperfect CSI. For the sake of
analytical tractability, the uncorrelated estimation errors were assumed and
only power allocation was involved into precoding without accounting for
beamforming design in [21]. However, the developed methodologies in [18, 19,
20, 21] are based on either approximations or simulations, and they cannot be
directly applied into the scenario of large-scale SCNs as well. Moreover,
although a large-scale dense network was considered in [22], the proposed
analytical approach is specific to the multiple-input single-output
(MISO)-NOMA system. In contrast, the performance analysis of the MIMO-NOMA
system will introduce the difficulty of handling random matrix variate.

This motivates us to trigger the performance investigations of MIMO-NOMA
enhanced large-scale SCNs with imperfect CSI at both the transmit and receive
sides. More specifically, Poisson point process (PPP) and Kronecker model are
utilized to capture the randomness of small cell BSs’ positions and the impact
of estimation error correlation, respectively. By employing stochastic
geometry and integral transform, the exact outage probabilities are derived in
closed-form and insightful results are extracted by conducting asymptotic
outage analyses. With the analytical results, the optimal design of MIMO-NOMA
enhanced SCNs is empowered. To the best of our knowledge, this is the first
work that copes with MIMO-NOMA enhanced SCNs under imperfect CSI.
Additionally, in contrast to [13, 12, 21], a more flexible MIMO-NOMA
transmission framework with variable number of clusters/groups is established
in this paper. Furthermore, it is worth noting that the Gaussian approximation
used in [21] cannot be applied in this paper because of the nonexistence of
mean and variance of the aggregated inter-cell interference.

###  I-B Contributions

To recapitulate, five main contributions of this paper are enumerated below.

  1. 1.

Exact compact expressions are derived for the outage probabilities of MIMO-
NOMA enhanced SCNs, where two NOMA grouping policies are considered, namely
random NOMA grouping and distanced-based NOMA grouping.

  2. 2.

The exact outage analyses show that the average outage probability is
irrespective of the intensity of BSs in the interference-limited regime.
Whereas, a counterintuitive behavior that the intensity of BSs below a
threshold is harmful to the average outage performance is observed.

  3. 3.

As the channel uncertainty vanishes, the target transmission rates should be
upper bounded to guarantee an arbitrarily low outage probability in the
absence of the inter-cell interference.

  4. 4.

The correlation coefficient has a positive effect on the outage probability
under the high quality of CSI. Conversely, the correlation coefficient impairs
the outage performance if the quality of CSI is low.

  5. 5.

The goodput maximization of the MIMO-NOMA system is considered by jointly
devising the precoding matrix, receiver filters and transmission rates.
Numerical results exhibit that the proposed scheme outperforms other
benchmarking schemes especially under significant difference between channel
gains, including MIMO-OMA with precoding, MIMO-OMA without precoding and MIMO-
NOMA without precoding.

###  I-C Paper Organization and Notation

The remaining of the paper is organized as follows. Section II delineates the
model of MIMO-NOMA enhanced SCNs. The exact and asymptotic outage analyses are
carried out in Sections III and IV, respectively. Section V devotes itself to
the optimal system design. Numerical results are then presented for
validations in Section VI. Section VII finally concludes this paper.

_Notation_ : The following notations shall be adopted throughout this paper.
The uppercase and lowercase boldface letters denote matrices and vectors,
respectively. ${\bf X}^{\mathrm{T}}$, ${\bf X}^{\mathrm{H}}$, ${\bf X}^{-1}$
and ${\bf X}^{1/2}$ stands for the transpose, conjugate transpose, matrix
inverse and Hermitian square root of matrix ${\bf X}$, respectively.
${\rm{Tr}}(\cdot)$, ${\mathrm{det}}(\cdot)$ and ${\mathrm{diag}}(\cdot)$ refer
to trace, determinant and diagonalization operations, respectively.
$\left\|\cdot\right\|$ represents the Euclidean/Frobenius norm of a
vector/matrix. ${{\mathbb{E}}_{A}}\\{\cdot\\}$ denotes the operator of the
expectation taking over random variable $A$, and the subscription $A$ is
sometimes omitted without leading to ambiguity. $\mathbb{C}^{m\times n}$ is
the set of $m\times n$-dimensional complex matrices. ${\cal
CN}({\boldsymbol{\mu},\boldsymbol{\Sigma}})$ denotes the distribution of a
circularly-symmetric complex Gaussian (CSCG) random variables with mean vector
$\boldsymbol{\mu}$ and covariance matrix $\boldsymbol{\Sigma}$.
$\mathbf{0}_{n}$, $\mathbf{1}_{n}$ and $\mathbf{I}_{n}$ represent an $n\times
1$ all-zeros vector, an $n\times 1$ all-ones vector and an $n\times n$
identity matrix, respectively, and the dimensionality $n$ is sometimes dropped
out for convenience. The symbol ${\rm i}=\sqrt{-1}$ refers to the imaginary
unit. ${\rm Re}\\{\cdot\\}$ represents the real part of a complex number.
$\bigcup{(\cdot)}$ stands for the union operation.
$H^{m,n}_{p,q}(\left.x\right|\cdot)$ represents Fox’s H-function. Any other
notations are deferred to define in the place where they appear.

##  II System Model

This paper considers a MIMO-NOMA enhanced SCNs with imperfect CSI. We assume
that the positions of the small cell BSs obey a homogeneous PPP $\Phi_{b}$ of
intensity $\lambda_{b}$ and the users are also distributed according to a
homogenous PPP $\Phi_{u}$ of intensity $\lambda_{u}$. Each user is associated
with its closest BS[23, 24]. The BSs and the users are equipped with $M$ and
$N$ antennas, respectively. By considering heavily loaded networks, i.e.,
$\lambda_{u}\gg\lambda_{b}$, it is reasonable to assume that there are at
least $2K$ users in a typical Voronoi cell and $K\leq\min\\{M,N\\}$[25, 26].
To accommodate $K$ pairs of NOMA users over MIMO channels within the coverage
of each cell, $2K$ users can be randomly chosen in the Voronoi cell.
Accordingly, the NOMA users are uniformly distributed within the Voronoi cell
[23]. In [27, 25, 28], the distance between a user and its serving BS
approximately follows a Rayleigh distribution with probability density
function (PDF)

$f_{d}(x)\approx 2c\lambda_{b}\pi xe^{-c\lambda_{b}\pi x^{2}},\,x\geq 0,$ (1)

where $c=5/4$. Hereinafter, the system model is introduced by splitting into
three subsections, including transmission model, imperfect CSI and received
SINR.

###  II-A Transmission Model

In this paper, the NOMA scheme is implemented in the power-domain. By
following the similar MIMO-NOMA framework to [13], the $K\times 1$
information-bearing vector ${{\bf{s}}}$ sent by a typical BS $z$ is expressed
as

${{\bf{s}}}=\left({{\beta_{1}}{s_{1}}+{\beta_{\tilde{1}}}{s_{\tilde{1}}}},\cdots,{{\beta_{k}}{s_{k}}+{\beta_{\tilde{k}}}{s_{\tilde{k}}}},\cdots,{{\beta_{K}}{s_{K}}+{\beta_{\tilde{K}}}{s_{\tilde{K}}}}\right)^{\rm
T},$ (2)

where ${s_{k}}$ and ${s_{\tilde{k}}}$ are the signals intended for the near
user and the far user of the $k$-th pair, respectively, and they are
normalized to unit power, i.e.,
${\mathbb{E}}\\{|s_{k}|^{2}\\}={\mathbb{E}}\\{|s_{k^{\prime}}|^{2}\\}\triangleq
1$, ${\beta_{k}}$ and ${\beta_{\tilde{k}}}$ correspond to the power allocation
coefficients with ${\beta_{k}}^{2}+{\beta_{\tilde{k}}}^{2}=1$. The received
signal at user $k$ (or $\tilde{k}$) can be expressed as

$\displaystyle{\bf{y}}_{k}$
$\displaystyle={\sqrt{P\ell\left({\left\|{z-{o_{k}}}\right\|}\right)}}{{\bf{H}}}_{zk}{{\bf{V}}}{{\bf{s}}}+{\sum\nolimits_{x\in{\Phi_{b}}\backslash\left\\{z\right\\}}}\sqrt{\rho_{I}\ell\left({\left\|{x-{o_{k}}}\right\|}\right)}{{\bf{1}}_{N}}{w_{x}}+{\bf{n}}_{k},$
(3)

where ${o_{k}}$ stands for the position of user $k$, $P$ denotes the transmit
power of each data stream, ${{\bf{H}}}_{zk}\in{\mathbb{C}}^{N\times M}$
represents the channel response matrix from the BS $z$ to user $k$,
$x\in{\Phi_{b}}\backslash\left\\{z\right\\}$ stands for the interfering BSs
(the co-channel BSs excluding $z$) and
${{{\bf{V}}}}=\left({{{\bf{v}}_{1}},\cdots,{{\bf{v}}_{K}}}\right)$ is the
normalized $M\times K$ precoding matrix, i.e., $\left\|{\bf v}_{k}\right\|=1$
for $k=1,\cdots,K$, the co-channel interference is modeled according to the
classical shot noise model due to the unavailable CSI between interfering BS
$x\in{\Phi_{b}}$ and user $k$ [13], $\ell\left(d\right)={d^{-\alpha}}$
captures path loss and $\alpha$ stands for the path loss exponent,
${{w}_{{x}}}$ denotes the normalized signal sent by the interferer, $\rho_{I}$
corresponds to the interference power, ${\bf n}_{k}$ denotes the complex-
valued additive Gaussian white noise with variance ${\sigma^{2}}{\bf I}_{M}$.

###  II-B Imperfect CSI

We assume that only partial CSI is available at both the transmitter and the
receiver. More precisely, the imperfect channel model for ${\bf H}_{zk}$ is
formulated as [24]

${\bf H}_{zk}=\hat{\bf H}_{zk}+{\bf E}_{zk}.$ (4)

where $\hat{\bf H}_{zk}$ is the channel state known to BS $z$ and user $k$,
and ${\bf E}_{zk}$ is the channel estimation error. Furthermore, a frequently
used error model is adopted for the estimation error ${\bf E}_{zk}$, i.e.,
Kronecker correlation model. In particular, the channel uncertainty matrix
${\bf E}_{zk}$ is given by

${\bf{E}}_{zk}={{\bf{R}}_{rk}}^{1/2}{{\bf{E}}_{w}}{{\bf{R}}_{tk}}^{1/2}.$ (5)

where ${\rm vec}({\bf E}_{w})\sim{\cal CN}({\bf 0},{\sigma_{h}}^{2}{\bf I})$.
It is worth noting that error correlations take place because of the transmit
and receive antenna structure. Accordingly, ${\bf{E}}_{zk}$ is a CSCG random
matrix, i.e., ${\rm vec}({\bf E}_{zk})\sim{\cal CN}({\bf
0},{\sigma_{h}}^{2}({{\bf R}_{tk}^{\rm T}\otimes{\bf R}_{rk}}))$. It is
assumed that the partial channel informations $\hat{\bf H}_{zk}$, ${\bf
R}_{tk}$ and ${\bf R}_{rk}$ are known to BS $z$ and user $k$. To characterize
the quality of CSI, we introduce the channel K factor defined as [24]

$\mathcal{K}_{zk}=\frac{\|{\hat{\bf
H}_{zk}}\|_{F}^{2}}{\mathbb{E}\\{\left\|{{\bf
E}_{zk}}\right\|_{F}^{2}\\}}=\frac{\|{\hat{\bf
H}_{zk}}\|_{F}^{2}}{{\sigma_{h}}^{2}{\rm Tr}{({{\bf R}_{tk}})}{\rm Tr}({\bf
R}_{rk})},$ (6)

where $\mathcal{K}_{zk}$ represents the power ratio of the known channel part
to the unknown part, which quantifies the quality of CSI. In addition, the
similar definition also applies to the imperfect CSI associated with user
$\tilde{k}$.

###  II-C Received SINR

Similarly to [13], the superimposed signals at the receive sides are
successively decoded in a descending order of link distances. Specifically,
with regard to the near user $k$, by applying receiver filter ${\bf u}_{k}$ to
(3), the desired data stream $k$ can be written as

$\displaystyle{{\hat{s}}_{k}}=$
$\displaystyle\sqrt{P\ell\left({{d_{k}}}\right)}{{\bf{u}}_{k}}^{\rm{H}}\left({{{{\bf{\hat{H}}}}_{zk}}+{{\bf{E}}_{zk}}}\right){{\bf{v}}_{k}}\left({{\beta_{k}}{s_{k}}+{\beta_{\tilde{k}}}{s_{\tilde{k}}}}\right)$
$\displaystyle+\underbrace{\sqrt{P\ell\left({{d_{k}}}\right)}\sum\limits_{i=1,i\neq
k}^{K}{{{\bf{u}}_{k}}^{\rm{H}}\left({{{{\bf{\hat{H}}}}_{zk}}+{{\bf{E}}_{zk}}}\right){{\bf{v}}_{i}}\left({{\beta_{i}}{s_{i}}+{\beta_{\tilde{i}}}{s_{\tilde{i}}}}\right)}}_{{\rm{intra-
cell}}\,{\rm{interference}}}$
$\displaystyle+\underbrace{{\sum_{x\in{\Phi_{b}}\backslash\left\\{z\right\\}}}\sqrt{{\rho_{I}}\ell\left({\left\|{x-{o_{k}}}\right\|}\right)}{{\bf{u}}_{k}}^{\rm{H}}{{\bf{1}}_{N}}{w_{x}}}_{{\rm{inter-
cell\,interference}}}+{{\bf{u}}_{k}}^{\rm{H}}{{\bf{n}}_{k}}.$ (7)

where $\left\|{z-{o_{k}}}\right\|=d_{k}$. By following the NOMA principle, the
technique of successive interference cancellation (SIC) is adopted to cancel
out the message intended to user $\tilde{k}$ with SINR given by

$\displaystyle{{\rm
SINR}_{k\to\tilde{k}}}=\frac{{P\ell\left({{d_{k}}}\right){{\left|{{{\bf{u}}_{k}}^{\rm{H}}{{{\bf{\hat{H}}}}_{zk}}{{\bf{v}}_{k}}}\right|}^{2}}{\beta_{\tilde{k}}}^{2}}}{{P\ell\left({{d_{k}}}\right)\left(\begin{array}[]{c}{\left|{{{\bf{u}}_{k}}^{\rm{H}}{{\bf{E}}_{zk}}{{\bf{v}}_{k}}}\right|^{2}}{\beta_{\tilde{k}}}^{2}+{\left|{{{\bf{u}}_{k}}^{\rm{H}}\left({{{{\bf{\hat{H}}}}_{zk}}+{{\bf{E}}_{zk}}}\right){{\bf{v}}_{k}}}\right|^{2}}{\beta_{k}}^{2}\\\
+\sum\nolimits_{i=1,i\neq
k}^{K}{{{\left|{{{\bf{u}}_{k}}^{\rm{H}}\left({{{{\bf{\hat{H}}}}_{zk}}+{{\bf{E}}_{zk}}}\right){{\bf{v}}_{i}}}\right|}^{2}}}\end{array}\right)+{I_{k}}+{\sigma_{k}^{2}}}},$
(10)

where $\sigma_{k}^{2}=\sigma^{2}\left\|{{{\bf{u}}_{k}}}\right\|^{2}$ and
${I_{k}}$ is the power of the inter-cell interference observed by user $k$,
such that

${I_{k}}={\rho_{I}}{\left|{{{\bf{u}}_{k}}^{\rm{H}}{{\bf{1}}_{N}}}\right|^{2}}{\sum\nolimits_{x\in{\Phi_{b}}\backslash\left\\{z\right\\}}}\ell\left({\left\|{x-{o_{k}}}\right\|}\right).$
(11)

By subtracting the message $s_{\tilde{k}}$, the near user decodes its own
message $s_{k}$ with SINR given by

$\displaystyle{{\rm
SINR}_{k}}=\frac{{P\ell\left({{d_{k}}}\right){{\left|{{{\bf{u}}_{k}}^{\rm{H}}{{{\bf{\hat{H}}}}_{zk}}{{\bf{v}}_{k}}}\right|}^{2}}{\beta_{k}}^{2}}}{{P\ell\left({{d_{k}}}\right)\left({{{\left|{{{\bf{u}}_{k}}^{\rm{H}}{{\bf{E}}_{zk}}{{\bf{v}}_{k}}}\right|}^{2}}+\sum\nolimits_{i=1,i\neq
k}^{K}{{{\left|{{{\bf{u}}_{k}}^{\rm{H}}\left({{{{\bf{\hat{H}}}}_{zk}}+{{\bf{E}}_{zk}}}\right){{\bf{v}}_{i}}}\right|}^{2}}}}\right)+{I_{k}}+{\sigma_{k}^{2}}}}.$
(12)

In an analogous way, the far user $\tilde{k}$ decodes its own message with
SINR given by

${{\rm
SINR}_{\tilde{k}}}=\frac{{P\ell\left({{d_{\tilde{k}}}}\right){{\left|{{{\bf{u}}_{\tilde{k}}}^{\rm{H}}{{{\bf{\hat{H}}}}_{z\tilde{k}}}{{\bf{v}}_{k}}}\right|}^{2}}{\beta_{\tilde{k}}}^{2}}}{{P\ell\left({{d_{\tilde{k}}}}\right)\left(\begin{array}[]{c}{\left|{{{\bf{u}}_{\tilde{k}}}^{\rm{H}}{{\bf{E}}_{z\tilde{k}}}{{\bf{v}}_{k}}}\right|^{2}}{\beta_{\tilde{k}}}^{2}+{\left|{{{\bf{u}}_{\tilde{k}}}^{\rm{H}}\left({{{{\bf{\hat{H}}}}_{z\tilde{k}}}+{{\bf{E}}_{z\tilde{k}}}}\right){{\bf{v}}_{k}}}\right|^{2}}{\beta_{k}}^{2}\\\
+\sum\nolimits_{i=1,i\neq
m}^{K}{{{\left|{{{\bf{u}}_{\tilde{k}}}^{\rm{H}}\left({{{{\bf{\hat{H}}}}_{z\tilde{k}}}+{{\bf{E}}_{z\tilde{k}}}}\right){{\bf{v}}_{i}}}\right|}^{2}}}\end{array}\right)+{I_{\tilde{k}}}+{\sigma_{\tilde{k}}^{2}}}},$
(13)

where
${\sigma_{\tilde{k}}^{2}}=\sigma^{2}\left\|{{{\bf{u}}_{\tilde{k}}}}\right\|^{2}$
and ${\bf u}_{\tilde{k}}$ refers to the detection filter applied at user
$\tilde{k}$, $\left\|{z-{o_{\tilde{k}}}}\right\|=d_{\tilde{k}}$ and the
interference term ${I_{\tilde{k}}}$ reads as

${I_{\tilde{k}}}={\rho_{I}}{\left|{{{\bf{u}}_{\tilde{k}}}^{\rm{H}}{{\bf{1}}_{N}}}\right|^{2}}{\sum\nolimits_{x\in{\Phi_{b}}\backslash\left\\{z\right\\}}}\ell\left({\left\|{x-{o_{\tilde{k}}}}\right\|}\right).$
(14)

It is worth noting that the aggregated interferences ${I_{k}}$ and
${I_{\tilde{k}}}$ are correlated due to the same interfering BSs, i.e.,
$x\in\Phi_{b}\backslash\\{z\\}$.

##  III Outage Probability

Since the outage probability is the key performance metric, this section is
dedicated to conducting the outage analysis for MIMO-NOMA enhanced SCNs with
imperfect CSI. The outage probabilities of users $\tilde{k}$ and $k$ are
separately derived as follows.

###  III-A Outage Probability of the Far User $\tilde{k}$

####  III-A1 Conditional Outage Probability

The outage event at the far user $\tilde{k}$ happens if the channel capacity
is less than the preset transmission rate $R_{\tilde{k}}$. Given the
collection of channel information $\mathcal{H}=\\{\hat{\bf H}_{zk},{\bf
R}_{rk},{\bf R}_{tk},\hat{\bf H}_{z\tilde{k}},{\bf R}_{r\tilde{k}},{\bf
R}_{t\tilde{k}},d_{k},d_{\tilde{k}},\alpha,\sigma^{2},{\sigma_{h}}^{2},\rho_{I},\lambda_{b}\\}$,
the conditional outage probability of user $\tilde{k}$ can be expressed as

$\displaystyle{p_{\tilde{k}|\mathcal{H}}}$
$\displaystyle={\Pr\left({{{\log}_{2}}\left({1+{{\rm
SINR}_{\tilde{k}}}}\right)<{R_{\tilde{k}}}}\right)}=1-\underbrace{\Pr\left({{{\log}_{2}}\left({1+{{\rm
SINR}_{\tilde{k}}}}\right)\geq{R_{\tilde{k}}}}\right)}_{{q_{\tilde{k}|{\mathcal{H}}}}},$
(15)

where ${{q_{\tilde{k}|{\mathcal{H}}}}}$ stands for the conditional successful
probability given the statistical channel information $\mathcal{H}$. We thus
proceed to derive ${{q_{\tilde{k}|{\mathcal{H}}}}}$ next. By using (13), it
follows that

$\displaystyle{q_{\tilde{k}|{\mathcal{H}}}}=\Pr\left\\{\begin{array}[]{l}{\beta_{\tilde{k}}}^{2}{\left|{{\chi_{\tilde{k}}}}\right|^{2}}+{\beta_{k}}^{2}{\left|{{\mu_{\tilde{k}}}+{\chi_{\tilde{k}}}}\right|^{2}}+\sum\limits_{i=1,i\neq
m}^{K}{{{\left|{{\mu_{\tilde{i}}}+{\chi_{\tilde{i}}}}\right|}^{2}}}\leq\frac{{{{\left|{{\mu_{\tilde{k}}}}\right|}^{2}}{\beta_{\tilde{k}}}^{2}}}{{{2^{{R_{\tilde{k}}}}}-1}}-\frac{{{{\sigma_{\tilde{k}}^{2}}}}}{{P\ell\left({{d_{\tilde{k}}}}\right)}}-\frac{{{I_{\tilde{k}}}}}{{P\ell\left({{d_{\tilde{k}}}}\right)}}\end{array}\right\\},$
(17)

where
${\mu_{\tilde{i}}}={{\bf{u}}_{\tilde{k}}}^{\rm{H}}{{{\bf{\hat{H}}}}_{z\tilde{k}}}{{\bf{v}}_{i}}$
and
${\chi_{\tilde{i}}}={{\bf{u}}_{\tilde{k}}}^{\rm{H}}{{\bf{E}}_{z\tilde{k}}}{{\bf{v}}_{i}}$.
Then, by using the identity
${\beta_{\tilde{k}}}^{2}{\left|{{\chi_{\tilde{k}}}}\right|^{2}}+{\beta_{k}}^{2}{\left|{{\mu_{\tilde{k}}}+{\chi_{\tilde{k}}}}\right|^{2}}={\left|{{\beta_{k}}^{2}{\mu_{\tilde{k}}}+{\chi_{\tilde{k}}}}\right|^{2}}+{\beta_{k}}^{2}{\beta_{\tilde{k}}}^{2}{\left|{{\mu_{\tilde{k}}}}\right|^{2}}$,
(17) can be simplified as

$\displaystyle{q_{\tilde{k}|{\mathcal{H}}}}$
$\displaystyle=\Pr\left\\{{{{\left({{\boldsymbol{\tilde{\chi}}}+{\boldsymbol{\tilde{\nu}}}}\right)}^{\rm{H}}}\left({{\boldsymbol{\tilde{\chi}}}+{\boldsymbol{\tilde{\nu}}}}\right)\leq{\tau_{\tilde{k}}}-\frac{{{I_{\tilde{k}}}}}{{P\ell\left({{d_{\tilde{k}}}}\right)}}}\right\\},$
(18)

where
${\boldsymbol{\tilde{\chi}}}={\left({{{\chi}_{\tilde{1}}},\cdots,{{\chi}_{\tilde{K}}}}\right)^{\rm{T}}}={\left({{{\bf{u}}_{\tilde{k}}}^{\rm{H}}{{\bf{E}}_{z\tilde{k}}}{\bf{V}}}\right)^{\rm{T}}}$,
${\tau_{\tilde{k}}}=\left({{1}/{({{2^{{R_{\tilde{k}}}}}-1})}-{\beta_{k}}^{2}}\right){\beta_{\tilde{k}}}^{2}{\left|{{\mu_{\tilde{k}}}}\right|^{2}}-{{{\sigma_{\tilde{k}}^{2}}}}/{({P\ell\left({{d_{\tilde{k}}}}\right)})}$
and
${\boldsymbol{\tilde{\nu}}}={\left({{\mu_{\tilde{1}}},\cdots,{\mu_{\tilde{k-1}}},{\beta_{k}}^{2}{\mu_{\tilde{k}}},{\mu_{\tilde{k+1}}},\cdots,{\mu_{\tilde{K}}}}\right)^{\rm{T}}}$.
In order to guarantee ${q_{\tilde{k}|{\mathcal{H}}}}<1$ from (18), the
condition ${\beta_{k}}^{2}\left({{2^{{R_{\tilde{k}}}}}-1}\right)<1$ should be
satisfied. Apparently, ${\boldsymbol{\tilde{\chi}}}$ is a CSCG random vector
with mean zero and covariance matrix
${\bf{\tilde{\Sigma}}}={\sigma_{h}^{2}}{{\bf{u}}_{\tilde{k}}}^{\rm{H}}{{\bf{R}}_{r\tilde{k}}}{{\bf{u}}_{\tilde{k}}}{\left({{{\bf{V}}^{\rm{H}}}{{\bf{R}}_{t\tilde{k}}}{\bf{V}}}\right)^{\rm{T}}}$,
i.e., ${\boldsymbol{\tilde{\chi}}}\sim{\cal CN}(\bf 0,\tilde{\bf\Sigma})$[24,
eq.13]. In analogous to characteristic function used in [29], applying inverse
Laplace transform to (18) yields

$\displaystyle q_{\tilde{k}|{\mathcal{H}}}$
$\displaystyle={{\mathbb{E}}_{{I_{\tilde{k}}}}}\left\\{{\int\limits_{{{\mathbb{C}}^{K}}}{u\left({{\tau_{\tilde{k}}}-\frac{{{I_{\tilde{k}}}}}{{P\ell\left({{d_{\tilde{k}}}}\right)}}}-{{\left({{\bf{x}}+{\boldsymbol{\tilde{\nu}}}}\right)}^{\rm{H}}}\left({{\bf{x}}+{\boldsymbol{\tilde{\nu}}}}\right)\right){f_{{\boldsymbol{\tilde{\chi}}}}}\left({\bf{x}}\right)}d{\bf{x}}}\right\\}$
$\displaystyle={{\mathbb{E}}_{{I_{\tilde{k}}}}}\left\\{{\int\limits_{{{\mathbb{C}}^{K}}}{\frac{1}{{2\pi{\rm{i}}}}\int\limits_{\tilde{c}-{\rm{i}}\infty}^{\tilde{c}+{\rm{i}}\infty}{\frac{1}{s}{e^{s\left({{\tau_{\tilde{k}}}-\frac{{{I_{\tilde{k}}}}}{{P\ell\left({{d_{\tilde{k}}}}\right)}}}-{{\left({{\bf{x}}+{\boldsymbol{\tilde{\nu}}}}\right)}^{\rm{H}}}\left({{\bf{x}}+{\boldsymbol{\tilde{\nu}}}}\right)\right)}}ds}{f_{{\boldsymbol{\tilde{\chi}}}}}\left({\bf{x}}\right)}d{\bf{x}}}\right\\},\,\tilde{c}>0,$
(19)

where the last step holds by using
$u(x)=\frac{1}{{2\pi{\rm{i}}}}\int_{\tilde{c}-{\rm{i}}\infty}^{\tilde{c}+{\rm{i}}\infty}{{1}/{s}{e^{sx}}ds}$,
$u(x)$ stands for the step unit function and
${f_{{\boldsymbol{\tilde{\chi}}}}}\left({\bf{x}}\right)$ represents the joint
PDF of ${\boldsymbol{\tilde{\chi}}}$. Since the joint PDF of
${\boldsymbol{\tilde{\chi}}}$ is
${f_{\boldsymbol{\tilde{\chi}}}}\left({\bf{x}}\right)={\exp({-{{\bf{x}}^{\rm{H}}}{{\bf{\tilde{\Sigma}}}^{-1}}{\bf{x}}})}/{({{\pi^{K}}\det({\bf{\tilde{\Sigma}}})})}$,
(III-A1) can be further rewritten as

$\displaystyle q_{\tilde{k}|{\mathcal{H}}}=$
$\displaystyle\frac{1}{{{\pi^{K}}\det\left({{\bf{\tilde{\Sigma}}}}\right)}}\frac{1}{{2\pi{\rm{i}}}}\int\limits_{\tilde{c}-{\rm{i}}\infty}^{\tilde{c}+{\rm{i}}\infty}{\frac{{{e^{s{\tau_{\tilde{k}}}}}}}{s}{{\mathbb{E}}_{{I_{\tilde{k}}}}}\left\\{{{e^{-s\frac{{{I_{\tilde{k}}}}}{{P\ell\left({{d_{\tilde{k}}}}\right)}}}}}\right\\}ds}\int\limits_{{{\mathbb{C}}^{K}}}{{e^{-s{{\left({{\bf{x}}+{\boldsymbol{\tilde{\nu}}}}\right)}^{\rm{H}}}\left({{\bf{x}}+{\boldsymbol{\tilde{\nu}}}}\right)-{{\bf{x}}^{\rm{H}}}{{{\bf{\tilde{\Sigma}}}}^{-1}}{\bf{x}}}}}d{\bf{x}}.$
(20)

By using the sum of two quadratic forms [30, eq.8.1.7] and the important
integration [29, eq.22], the inner integral in (20) can be computed as

$\displaystyle\int\limits_{{\mathbb{C}^{K}}}{{e^{-s{{\left({{\bf{x}}+{\boldsymbol{\tilde{\nu}}}}\right)}^{\rm{H}}}\left({{\bf{x}}+{\boldsymbol{\tilde{\nu}}}}\right)-{{\bf{x}}^{\rm{H}}}{{{\bf{\tilde{\Sigma}}}}^{-1}}{\bf{x}}}}}d{\bf{x}}$
$\displaystyle=\int\limits_{{\mathbb{C}^{K}}}{{e^{{{-\left({{\bf{x}}+{\bf{\tilde{a}}}}\right)}^{\rm{H}}}\left({s{\bf{I}}+{{{\bf{\tilde{\Sigma}}}}^{-1}}}\right)\left({{\bf{x}}+{\bf{\tilde{a}}}}\right)-{{{\boldsymbol{\tilde{\nu}}}}^{\rm{H}}}{{\left({\frac{1}{s}{\bf{I}}+{\bf{\tilde{\Sigma}}}}\right)}^{-1}}{\boldsymbol{\tilde{\nu}}}}}d{\bf{x}}}$
$\displaystyle=\frac{{{\pi^{K}}}}{{\det\left({s{\bf{I}}+{{{\bf{\tilde{\Sigma}}}}^{-1}}}\right)}}{e^{-{{{\boldsymbol{\tilde{\nu}}}}^{\rm{H}}}{{\left({\frac{1}{s}{\bf{I}}+{\bf{\tilde{\Sigma}}}}\right)}^{-1}}{\boldsymbol{\tilde{\nu}}}}}.$
(21)

Furthermore, by using the Laplace functional of ${I_{\tilde{k}}}$, the
expectation term in the above can be obtained as [31, eq.5.11]

$\displaystyle{{\mathbb{E}}_{{I_{\tilde{k}}}}}\left\\{{{e^{-s\frac{{{I_{\tilde{k}}}}}{{P\ell\left({{d_{\tilde{k}}}}\right)}}}}}\right\\}$
$\displaystyle={e^{-\pi{\lambda_{b}}{\omega_{\tilde{k}}}{d_{\tilde{k}}}^{2}{s^{\frac{2}{\alpha}}}}},$
(22)

where
${\omega_{\tilde{k}}}=\Gamma\left({1-\frac{2}{\alpha}}\right){\left({\frac{{{\rho_{I}}}}{P}{{\left|{{{\bf{u}}_{\tilde{k}}}^{\rm{H}}{{\bf{1}}_{K}}}\right|}^{2}}}\right)^{\frac{2}{\alpha}}}$.
Thus, by substituting (III-A1) and (22) into (20), we obtain

$\displaystyle q_{\tilde{k}|{\mathcal{H}}}$
$\displaystyle=\frac{1}{{\det\left({{\bf{\tilde{\Sigma}}}}\right)}}\frac{1}{{2\pi{\rm{i}}}}\int\limits_{\tilde{c}-{\rm{i}}\infty}^{\tilde{c}+{\rm{i}}\infty}{\frac{{{e^{-{{{\boldsymbol{\tilde{\nu}}}}^{\rm{H}}}{{\left({\frac{1}{s}{\bf{I}}+{\bf{\tilde{\Sigma}}}}\right)}^{-1}}{\boldsymbol{\tilde{\nu}}}+{\tau_{\tilde{k}}}s{-\pi{\lambda_{b}}{\omega_{\tilde{k}}}{d_{\tilde{k}}}^{2}{s^{\frac{2}{\alpha}}}}}}}}{{s\det\left({s{\bf{I}}+{{{\bf{\tilde{\Sigma}}}}^{-1}}}\right)}}ds}.$
(23)

With the eigenvalue decomposition
${\bf{\tilde{\Sigma}}}=\boldsymbol{\tilde{\Psi}}\boldsymbol{\tilde{\Delta}}{\boldsymbol{\tilde{\Psi}}}^{\rm{H}}$,
(23) reduces to

$q_{\tilde{k}|{\mathcal{H}}}=\frac{1}{{2\pi{\rm{i}}}}\int\limits_{\tilde{c}-{\rm{i}}\infty}^{\tilde{c}+{\rm{i}}\infty}{{e^{{\tau_{\tilde{k}}}s}}\underbrace{\frac{{\prod\nolimits_{i=1}^{K}{{e^{-\frac{{s{{\left|{{\zeta_{\tilde{i}}}}\right|}^{2}}}}{{1+s{\delta_{\tilde{i}}}}}}}}{e^{-\pi{\lambda_{b}}{\omega_{\tilde{k}}}{d_{\tilde{k}}}^{2}{s^{\frac{2}{\alpha}}}}}}}{{s\prod\nolimits_{i=1}^{K}{\left({1+s{\delta_{\tilde{i}}}}\right)}}}}_{{F_{\tilde{k}}}\left(s\right)}ds},$
(24)

where
$\boldsymbol{\tilde{\Delta}}={\rm{diag}}\left({{\delta_{\tilde{1}}},\cdots,{\delta_{\tilde{K}}}}\right)$
and
$\boldsymbol{\tilde{\Psi}}^{\rm{H}}\boldsymbol{\tilde{\nu}}={\left({{\zeta_{\tilde{1}}},\cdots,{\zeta_{\tilde{K}}}}\right)^{\rm{H}}}$.
In order to calculate $q_{\tilde{k}|{\mathcal{H}}}$, a numerical inversion of
Laplace transform developed in [32] can be applied herein. More specifically,
by using the method of Abate and Whitt, (24) can be approximated with an
arbitrarily small discretization error as[32]

$\displaystyle q_{\tilde{k}|{\mathcal{H}}}\approx$
$\displaystyle\frac{{{2^{-M}}{e^{A/2}}}}{{{\tau_{\tilde{k}}}}}\sum\limits_{m=0}^{M}{{{M}\choose{m}}}\left(\begin{array}[]{l}\frac{1}{2}{\rm{Re}}\left\\{{{F_{\tilde{k}}}\left({\frac{A}{{2{\tau_{\tilde{k}}}}}}\right)}\right\\}+\sum\limits_{n=1}^{Q+m}{{{\left({-1}\right)}^{n}}{\rm{Re}}\left\\{{{F_{\tilde{k}}}\left({\frac{{A+2n\pi{\rm{i}}}}{{2{\tau_{\tilde{k}}}}}}\right)}\right\\}}\end{array}\right),$
(26)

where $M$ refers to the number of Euler summation terms, $Q$ refers to the
truncation order, the discretization error is bounded by $|\epsilon|\leq
e^{-A}/(1-e^{-A})$ and the truncation error is manageable by properly choosing
$M$ and $Q$. To have a discretization error up to $10^{-10}$, $A$ is set to
$A\approx 23$. As suggested by [32], $M=11$ and $Q=15$ are typically chosen.

By putting (26) into (15), we can get the final expression of the conditional
outage probability ${p_{k|\mathcal{H}}}$.

####  III-A2 Average Outage Probability

Furthermore, the average outage probability by taking the expectation over the
distance $d_{\tilde{k}}$ can be expressed as

${p_{\tilde{k}}}={{\mathbb{E}}_{{d_{\tilde{k}}}}}\left\\{{p_{\tilde{k}|{\mathcal{H}}}}\right\\}=1-{{\mathbb{E}}_{{d_{\tilde{k}}}}}\left\\{{q_{\tilde{k}|{\mathcal{H}}}}\right\\},$
(27)

where

$\displaystyle{\mathbb{E}_{{d_{\tilde{k}}}}}\left\\{{{q_{\tilde{k}|{\cal
H}}}}\right\\}=\frac{1}{{2\pi{\rm{i}}}}\int\limits_{\tilde{c}-{\rm{i}}\infty}^{\tilde{c}+{\rm{i}}\infty}{{e^{{{\bar{\tau}}_{\tilde{k}}}s}}}\underbrace{\frac{{\prod\nolimits_{i=1}^{K}{{e^{-\frac{{s{{\left|{{\zeta_{\tilde{i}}}}\right|}^{2}}}}{{1+s{\delta_{\tilde{i}}}}}}}}\overbrace{{\mathbb{E}_{{d_{\tilde{k}}}}}\left\\{{{e^{-\frac{{{{\sigma_{\tilde{k}}^{2}}}}}{{P\ell({{d_{\tilde{k}}}})}}s-\pi{\lambda_{b}}{\omega_{\tilde{k}}}{d_{\tilde{k}}}^{2}{s^{\frac{2}{\alpha}}}}}}\right\\}}^{{\varphi_{\tilde{k}}}\left(s\right)}}}{{s\prod\nolimits_{i=1}^{K}{\left({1+s{\delta_{\tilde{i}}}}\right)}}}}_{{g_{\tilde{k}}}\left(s\right)}ds,$
(28)

and
${{\bar{\tau}}_{\tilde{k}}}=\left({{1}/{({{2^{{R_{\tilde{k}}}}}-1})}-{\beta_{k}}^{2}}\right){\beta_{\tilde{k}}}^{2}{\left|{{\mu_{\tilde{k}}}}\right|^{2}}$.
Likewise, by invoking the method of Abate and Whitt, ${p_{\tilde{k}}}$ is
obtained as

$\displaystyle{p_{\tilde{k}}}\approx$ $\displaystyle
1-\frac{{{2^{-M}}{e^{A/2}}}}{{{\bar{\tau}_{\tilde{k}}}}}\sum\limits_{m=0}^{M}{{{M}\choose{m}}}\left(\begin{array}[]{l}\frac{1}{2}{\rm{Re}}\left\\{{{g_{\tilde{k}}}\left({\frac{A}{{2{\bar{\tau}_{\tilde{k}}}}}}\right)}\right\\}+\sum\limits_{n=1}^{Q+m}{{{\left({-1}\right)}^{n}}{\rm{Re}}\left\\{{{g_{\tilde{k}}}\left({\frac{{A+2n\pi{\rm{i}}}}{{2{\bar{\tau}_{\tilde{k}}}}}}\right)}\right\\}}\end{array}\right).$
(30)

It is worth noting that the policy of user pairing directly determines the
distribution of $d_{\tilde{k}}$. Under different grouping policies, we have
different forms of ${\varphi_{\tilde{k}}}\left(s\right)$. In the sequel,
${\varphi_{\tilde{k}}}\left(s\right)$ is derived by considering two grouping
policies, i.e., the random grouping and the distance-based grouping.

######  Theorem 1.

Regarding the random NOMA grouping policy, $2K$ users are randomly grouped
into $K$ pairs. Under this grouping policy,
${\varphi_{\tilde{k}}}\left(s\right)$ is given by

$\displaystyle{\varphi_{\tilde{k}}}\left(s\right)=$
$\displaystyle\frac{{2c}}{{c+{\omega_{\tilde{k}}}{{s}^{\frac{2}{\alpha}}}}}H_{1,1}^{1,1}\left({\left.{\frac{{{{\sigma_{\tilde{k}}^{2}}}}}{P}s{{\left({\pi{\lambda_{b}}\left({c+{\omega_{\tilde{k}}}{{s}^{\frac{2}{\alpha}}}}\right)}\right)}^{-\frac{\alpha}{2}}}}\right|\begin{array}[]{*{20}{c}}{\left({0,\frac{\alpha}{2}}\right)}\\\
{\left({0,1}\right)}\end{array}}\right)$ (33)
$\displaystyle-\frac{{2c}}{{2c+{\omega_{\tilde{k}}}{{s}^{\frac{2}{\alpha}}}}}H_{1,1}^{1,1}\left({\left.{\frac{{{{\sigma_{\tilde{k}}^{2}}}}}{P}s{{\left({\pi{\lambda_{b}}\left({2c+{\omega_{\tilde{k}}}{{s}^{\frac{2}{\alpha}}}}\right)}\right)}^{-\frac{\alpha}{2}}}}\right|\begin{array}[]{*{20}{c}}{\left({0,\frac{\alpha}{2}}\right)}\\\
{\left({0,1}\right)}\end{array}}\right),$ (36)

where $H_{p,q}^{m,n}(\cdot)$ represents the Fox H-function [33, Def.1.1].
Moreover, under the interference-limited regime, i.e., $\sigma^{2}/P\to 0$,
(33) reduces to

${\varphi_{\tilde{k}}}\left(s\right)=\frac{{2c}}{{c+{\omega_{\tilde{k}}}{{s}^{\frac{2}{\alpha}}}}}-\frac{{2c}}{{2c+{\omega_{\tilde{k}}}{{s}^{\frac{2}{\alpha}}}}}.$
(37)

###### Proof.

Since the distance between a user and its serving BS follows a Rayleigh
distribution as (1), the PDF of ${{d_{\tilde{k}}}}$ can be obtained by using
order statistics as

${f_{{d_{\tilde{k}}}}}\left(x\right)=2{F_{d}}(x){f_{d}}(x),$ (38)

where ${F_{d}}(x)=1-{e^{-c{\lambda_{b}}\pi{x^{2}}}}$ corresponds to the
cumulative distribution function (CDF) of $d$. Based on (38),
${\varphi_{\tilde{k}}}\left(s\right)$ can be derived after some algebraic
manipulations as

$\displaystyle{\varphi_{\tilde{k}}}\left(s\right)=$ $\displaystyle
2c{\lambda_{b}}\pi\int\limits_{0}^{\infty}{{e^{-\frac{{{{\sigma_{\tilde{k}}^{2}}}}}{P}s{y^{\frac{\alpha}{2}}}-\pi{\lambda_{b}}\left({{\omega_{\tilde{k}}}{s^{\frac{2}{\alpha}}}+c}\right)y}}dy}-2c{\lambda_{b}}\pi\int\limits_{0}^{\infty}{{e^{-\frac{{{{\sigma_{\tilde{k}}^{2}}}}}{P}s{y^{\frac{\alpha}{2}}}-\pi{\lambda_{b}}\left({{\omega_{\tilde{k}}}{s^{\frac{2}{\alpha}}}+2c}\right)y}}dy}.$
(39)

By using the Parseval’s type property of Mellin transform [34, eq.8.3.21]
along with [33, Eq.1.37], (39) can be represented in terms of Fox’s H-function
as (33). Moreover, by using the asymptotic expansion of [33, Theorem.1.2], the
Fox’s H-function in (33) approaches to 1 as $\sigma^{2}/P\to 0$. Thus, we have
(37). ∎

######  Theorem 2.

Regarding the distance-based NOMA grouping policy, $2K$ users are assumed to
rank in an increasing order according to their distances from their serving
BS. Let $r_{\tilde{k}}$ denote the ranking order of user $\tilde{k}$. Under
this grouping policy, ${\varphi_{\tilde{k}}}\left(s\right)$ is given by

$\displaystyle{\varphi_{\tilde{k}}}\left(s\right)=$ $\displaystyle
c{r_{\tilde{k}}}{2K\choose{{r_{\tilde{k}}}}}\sum\limits_{l=0}^{{r_{\tilde{k}}}-1}{\frac{{{{\left({-1}\right)}^{l}}{{r_{\tilde{k}}}-1\choose{l}}}H_{1,1}^{1,1}\left({\left.{\frac{{{{\sigma_{\tilde{k}}^{2}}}}}{P}s{{\left({\pi{\lambda_{b}}\left({\left({2K-{r_{\tilde{k}}}+l+1}\right)c+{\omega_{\tilde{k}}}{s^{\frac{2}{\alpha}}}}\right)}\right)}^{-\frac{\alpha}{2}}}}\right|\begin{array}[]{*{20}{c}}{\left({0,\frac{\alpha}{2}}\right)}\\\
{\left({0,1}\right)}\end{array}}\right)}{{\left({2K-{r_{\tilde{k}}}+l+1}\right)c+{\omega_{\tilde{k}}}{s^{\frac{2}{\alpha}}}}}}.$
(42)

In the interference-limited regime, i.e., $\sigma^{2}/P\to 0$, it follows that

${\varphi_{\tilde{k}}}\left(s\right)=c{r_{\tilde{k}}}{{2K}\choose{r_{\tilde{k}}}}\sum\limits_{l=0}^{{r_{\tilde{k}}}-1}{\frac{{{{\left({-1}\right)}^{l}}{{{r_{\tilde{k}}}-1}\choose{l}}}}{{\left({2K-{r_{\tilde{k}}}+l+1}\right)c+{\omega_{\tilde{k}}}{s^{\frac{2}{\alpha}}}}}}.$
(43)

###### Proof.

By using order statistics, the PDF of ${{d_{\tilde{k}}}}$ is expressed as

${f_{{d_{\tilde{k}}}}}\left(x\right)={r_{\tilde{k}}}{{2K}\choose{r_{\tilde{k}}}}{\left({{F_{d}}(x)}\right)^{{r_{\tilde{k}}}-1}}{\left({1-{F_{d}}(x)}\right)^{2K-{r_{\tilde{k}}}}}{f_{d}}(x).$
(44)

Then, using (44) along with some basic manipulations results in (2). ∎

Surprisingly, from the above two theorems, the average outage probability
${p_{\tilde{k}}}$ is independent of the intensity of BSs in the interference-
limited regime because ${\varphi_{\tilde{k}}}\left(s\right)$ is irrelevant to
${{\lambda_{b}}}$. This is due to the fact that the effect of $\lambda_{b}$
has two-fold. On one hand, the increase of $\lambda_{b}$ yields severe co-
channel interference, which eventually leads to high outage probability. On
the other hand, the size of the Voronoi cell diminishes with the increase of
$\lambda_{b}$, thus it leads to the decay of path loss.

###  III-B Outage Probability of the Near User $k$

####  III-B1 Conditional Outage Probability

By conducting the SIC at the near user, the near user has to subtract the far
user’s signal prior to decoding its own signal. Hence, the conditional outage
probability of user $k$ given $\mathcal{H}$ is given by

$\displaystyle{p_{k|\mathcal{H}}}$
$\displaystyle=\Pr\left({{{\log}_{2}}\left({1+{{\rm
SINR}_{k\to\tilde{k}}}}\right)<{R_{\tilde{k}}}\bigcup{{{\log}_{2}}\left({1+{{\rm
SINR}_{k}}}\right)<{R_{k}}}}\right)$
$\displaystyle=1-\underbrace{\Pr\left({{{\rm
SINR}_{k\to\tilde{k}}}\geq{2^{{R_{\tilde{k}}}}}-1,{{{{\rm
SINR}_{k}}}\geq{2^{{R_{k}}}}-1}}\right)}_{{q_{k|{\mathcal{H}}}}}.$ (45)

Clearly, ${\rm{SIN}}{{\rm{R}}_{k\to\tilde{k}}}$ and ${\rm{SIN}}{{\rm{R}}_{k}}$
are strongly correlated due to the involvement of the same co-channel
interference and channel uncertainty. The existence of the spatial correlation
challenges the later outage analysis. By using (10) and (12),
${q_{k|{\mathcal{H}}}}$ can be rewritten as

${q_{k|{\mathcal{H}}}}=\Pr\left(\begin{array}[]{l}{\beta_{\tilde{k}}}^{2}{\left|{{\chi_{k}}}\right|^{2}}+{\beta_{k}}^{2}{\left|{{\mu_{k}}+{\chi_{k}}}\right|^{2}}+\sum\nolimits_{i=1,i\neq
k}^{K}{{{\left|{{\mu_{i}}+{\chi_{i}}}\right|}^{2}}}\leq{\theta_{\tilde{k}}}-\frac{{{I_{k}}}}{{P\ell\left({{d_{k}}}\right)}},\\\
{\left|{{\chi_{k}}}\right|^{2}}+\sum\nolimits_{i=1,i\neq
k}^{K}{{{\left|{{\mu_{i}}+{\chi_{i}}}\right|}^{2}}}\leq{\theta_{k}}-\frac{{{I_{k}}}}{{P\ell\left({{d_{k}}}\right)}}\end{array}\right),$
(46)

where
${\mu_{i}}={{\bf{u}}_{k}}^{\rm{H}}{{{\bf{\hat{H}}}}_{zk}}{{\bf{v}}_{i}}$,
${\chi_{i}}={{\bf{u}}_{k}}^{\rm{H}}{{\bf{E}}_{zk}}{{\bf{v}}_{i}}$,
${\theta_{k}}={{{{\left|{{\mu_{k}}}\right|}^{2}}{\beta_{k}}^{2}}}/{({{2^{{R_{k}}}}-1})}-{{{{\sigma_{k}^{2}}}}}/{({P\ell\left({{d_{k}}}\right)})}$
and
${\theta_{\tilde{k}}}={{{{\left|{{\mu_{k}}}\right|}^{2}}{\beta_{\tilde{k}}}^{2}}}/{({{2^{{R_{\tilde{k}}}}}-1})}-{{{{\sigma_{k}^{2}}}}}/{({P\ell\left({{d_{k}}}\right)})}$.

Similarly, it can be easily proved that
${\boldsymbol{{\chi}}}=(\chi_{1},\cdots,\chi_{K})^{\rm
T}={\left({{{\bf{u}}_{k}}^{\rm{H}}{{\bf{E}}_{zk}}{\bf{V}}}\right)^{\rm{T}}}$
is a CSCG random vector with mean zero and covariance matrix
${\bf{\Sigma}}={\sigma_{h}^{2}}{{\bf{u}}_{k}}^{\rm{H}}{{\bf{R}}_{rk}}{{\bf{u}}_{k}}{\left({{{\bf{V}}^{\rm{H}}}{{\bf{R}}_{tk}}{\bf{V}}}\right)^{\rm{T}}}$,
i.e., ${\boldsymbol{{\chi}}}\sim{\cal CN}(\bf 0,\bf\Sigma)$. With the help of
the step unit function, ${q_{k|{\mathcal{H}}}}$ can thus be expressed as

$\displaystyle{q_{k|{\mathcal{H}}}}={{\mathbb{E}}_{{I_{k}}}}\left\\{\begin{array}[]{l}\int\limits_{{\mathbb{C}^{K}}}{u\left({{\theta_{\tilde{k}}}-\frac{{{I_{k}}}}{{P\ell\left({{d_{k}}}\right)}}-{\beta_{\tilde{k}}}^{2}{{\left|{{\chi_{k}}}\right|}^{2}}-{\beta_{k}}^{2}{{\left|{{\mu_{k}}+{\chi_{k}}}\right|}^{2}}-\sum\nolimits_{i=1,i\neq
k}^{K}{{{\left|{{\mu_{i}}+{\chi_{i}}}\right|}^{2}}}}\right)}\\\ \times
u\left({{\theta_{k}}-\frac{{{I_{k}}}}{{P\ell\left({{d_{k}}}\right)}}-{{\left|{{\chi_{k}}}\right|}^{2}}-\sum\nolimits_{i=1,i\neq
k}^{K}{{{\left|{{\mu_{i}}+{\chi_{i}}}\right|}^{2}}}}\right){f_{\boldsymbol{\chi}}}\left({\bf{x}}\right)d{\bf{x}}\end{array}\right\\},$
(49)

where ${f_{\boldsymbol{\chi}}}\left({\bf{x}}\right)$ stands for the joint PDF
of ${\boldsymbol{\chi}}$ and
${f_{\boldsymbol{\chi}}}\left({\bf{x}}\right)={\exp({-{{\bf{x}}^{\rm{H}}}{{\bf{\Sigma}}^{-1}}{\bf{x}}})}/{({{\pi^{K}}\det({\bf{\Sigma}})})}$.
In analogous to (III-A1), by using the inverse Laplace transform of the step
unit function along with some rearrangements, ${q_{k|{\mathcal{H}}}}$ can be
obtained as

${q_{k|{\mathcal{H}}}}=\frac{1}{{{{\left({2\pi{\rm{i}}}\right)}^{2}}}}\int\limits_{c_{1}-{\rm{i}}\infty}^{c_{1}+{\rm{i}}\infty}{\int\limits_{c_{2}-{\rm{i}}\infty}^{c_{2}+{\rm{i}}\infty}{\frac{1}{s}\frac{1}{t}{e^{s{\theta_{\tilde{k}}}+t{\theta_{k}}}}{{\mathbb{E}}_{{I_{k}}}}\left\\{{{e^{-\left({s+t}\right)\frac{{{I_{k}}}}{{P\ell\left({{d_{k}}}\right)}}}}}\right\\}\mathcal{K}\left({s,t}\right)dsdt}},\,c_{1},c_{2}>0,$
(50)

where

$\displaystyle\mathcal{K}\left({s,t}\right)=\int\limits_{{\mathbb{C}^{M}}}{{e^{-(s{\beta_{\tilde{k}}}^{2}+t){{\left|{{x_{k}}}\right|}^{2}}-s{\beta_{k}}^{2}{{\left|{{\mu_{k}}+{x_{k}}}\right|}^{2}}-\left({s+t}\right)\sum\nolimits_{i=1,i\neq
k}^{K}{{{\left|{{\mu_{i}}+{x_{i}}}\right|}^{2}}}}}}{f_{\boldsymbol{\chi}}}\left({\bf{x}}\right)d{\bf{x}}.$
(51)

In regard to the expectation term in (50), it can be derived similarly to (22)
as

$\displaystyle{{\mathbb{E}}_{{I_{k}}}}\left\\{{{e^{-s\frac{{{I_{k}}}}{{P\ell\left({{d_{k}}}\right)}}}}}\right\\}$
$\displaystyle={e^{-\pi{\lambda_{b}}{\omega_{k}}{d_{k}}^{2}{(s+t)^{\frac{2}{\alpha}}}}},$
(52)

where
${\omega_{k}}=\Gamma\left({1-\frac{2}{\alpha}}\right){\left({\frac{{{\rho_{I}}}}{P}{{\left|{{{\bf{u}}_{k}}^{\rm{H}}{{\bf{1}}_{K}}}\right|}^{2}}}\right)^{\frac{2}{\alpha}}}$.
Additionally, in order to get $\mathcal{K}\left({s,t}\right)$, we define the
following vectors ${\boldsymbol{\mu}}=\left(\mu_{1},\cdots,\mu_{K}\right)^{\rm
T}$,

${\bf{A}}={\rm{diag}}\left({\overbrace{\underbrace{s+t,\cdots,s+t}_{\left({k-1}\right)-{\rm{entries}}},s{\beta_{k}}^{2},s+t,\cdots,s+t}^{K-{\rm{entries}}}}\right),$
(53)
${\bf{B}}={\rm{diag}}\left({\overbrace{\underbrace{0,\cdots,0}_{\left({k-1}\right)-{\rm{entries}}},s{\beta_{\tilde{k}}}^{2}+t,0,\cdots,0}^{K-{\rm{entries}}}}\right).$
(54)

Hereby, $\mathcal{K}\left({s,t}\right)$ can be rewritten as

$\displaystyle{\cal K}(s,t)$
$\displaystyle=\frac{1}{{{\pi^{K}}\det\left({\bf{\Sigma}}\right)}}\int\nolimits_{{\mathbb{C}^{K}}}{{e^{-{{\left({{\bf{x}}+{\boldsymbol{\mu}}}\right)}^{\rm{H}}}{\bf{A}}\left({{\bf{x}}+{\boldsymbol{\mu}}}\right)-{{\bf{x}}^{\rm{H}}}\left({{\bf{B}}+{{\bf{\Sigma}}^{-1}}}\right){\bf{x}}}}d{\bf{x}}}.$
(55)

By using the sum of two quadratic forms [30, Eq.8.1.7], we obtain

$\displaystyle{\left({{\bf{x}}+{\boldsymbol{\mu}}}\right)^{\rm{H}}}{\bf{A}}\left({{\bf{x}}+{\boldsymbol{\mu}}}\right)+{{\bf{x}}^{\rm{H}}}\left({{\bf{B}}+{{\bf{\Sigma}}^{-1}}}\right){\bf{x}}={\left({{\bf{x}}+{\boldsymbol{\nu}}}\right)^{\rm{H}}}{\bf{\Xi}}\left({{\bf{x}}+{\boldsymbol{\nu}}}\right)+\phi\left({s,t}\right),$
(56)

where ${\bf{\Xi}}={\bf{A}+\bf{B}}+{{\bf{\Sigma}}^{-1}}=(s+t){\bf
I}+{{\bf{\Sigma}}^{-1}}$,
${\boldsymbol{\nu}}={{\bf{\Xi}}^{-1}}{\bf{A}}\boldsymbol{\mu}$ and

$\displaystyle\phi\left({s,t}\right)$
$\displaystyle={{\boldsymbol{\mu}}^{\rm{H}}}{\bf{A}}{\boldsymbol{\mu}}-{{\boldsymbol{\mu}}^{\rm{H}}}{\bf{A}}{{\bf{\Xi}}^{-1}}{\bf
A}{\boldsymbol{\mu}}={{\boldsymbol{\mu}}^{\rm{H}}}{\bf{A}}{{\bf{\Xi}}^{-1}}{(\bf{B}+{{\bf{\Sigma}}^{-1}})}{\boldsymbol{\mu}}.$
(57)

Substituting (56) into (55) yields

$\displaystyle{\cal K}(s,t)$
$\displaystyle=\frac{{{e^{-\phi\left({s,t}\right)}}}}{{{\pi^{K}}\det\left({\bf{\Sigma}}\right)}}\int\limits_{{\mathbb{C}^{K}}}{{e^{-{{\left({{\bf{x}}+{\boldsymbol{\nu}}}\right)}^{\rm{H}}}{\bf{\Xi}}\left({{\bf{x}}+{\boldsymbol{\nu}}}\right)}}d{\bf{x}}}=\frac{{{e^{-\phi\left({s,t}\right)}}}}{{\det\left({\bf{\Sigma}}\right)\det\left({\bf{\Xi}}\right)}},$
(58)

where the last step holds by using [29, Eq. 22]. By plugging (52) and (58)
into (50), it results in

$\displaystyle{q_{k|\mathcal{H}}}=\frac{1}{{\det\left({\bf{\Sigma}}\right)}}\frac{1}{{{{\left({2\pi{\rm{i}}}\right)}^{2}}}}\int\limits_{c_{1}-{\rm{i}}\infty}^{c_{1}+{\rm{i}}\infty}{\int\limits_{c_{2}-{\rm{i}}\infty}^{c_{2}+{\rm{i}}\infty}{\frac{1}{st\det\left({\bf{\Xi}}\right)}{e^{s{\theta_{\tilde{k}}}+t{\theta_{k}}-\pi{\lambda_{b}}{\omega_{k}}{d_{k}}^{2}{{(s+t)}^{\frac{2}{\alpha}}}-\phi\left({s,t}\right)}}dsdt}}.$
(59)

Moreover, by capitalizing on eigenvalue decomposition
${\bf{\Sigma=\Psi\Delta}}{{\bf{\Psi}}^{\bf{H}}}$, defining
${\bf{\Psi}}=\left({{{\boldsymbol{\psi}}_{1}},\cdots,{{\boldsymbol{\psi}}_{K}}}\right)$
and
${\bf{\Delta}}={\rm{diag}}\left({{\delta_{1}},\cdots,{\delta_{K}}}\right)$,
$\phi\left({s,t}\right)$ in (59) can be simplified as

$\phi\left({s,t}\right)=\sum\limits_{i=1}^{K}{\frac{{{{\boldsymbol{\mu}}^{\rm{H}}}{\bf{A}}{{\boldsymbol{\psi}}_{i}}\left({{\delta_{i}}{{\boldsymbol{\psi}}_{i}}^{\rm{H}}{{\bf
B}{\boldsymbol{\mu}}}+{{\boldsymbol{\psi}}_{i}}^{\rm{H}}{\boldsymbol{\mu}}}\right)}}{{1+\left({s+t}\right){\delta_{i}}}}}.$
(60)

As a consequence, ${q_{k|{\mathcal{H}}}}$ is derived as

$\displaystyle{q_{k|{\mathcal{H}}}}$
$\displaystyle=\frac{1}{{{{\left({2\pi{\rm{i}}}\right)}^{2}}}}\int\limits_{c_{1}-{\rm{i}}\infty}^{c_{1}+{\rm{i}}\infty}{\int\limits_{c_{2}-{\rm{i}}\infty}^{c_{2}+{\rm{i}}\infty}{{e^{s{\theta_{\tilde{k}}}+t{\theta_{k}}}}\underbrace{\frac{{{e^{-\pi{\lambda_{b}}{\omega_{k}}{d_{k}}^{2}{{(s+t)}^{\frac{2}{\alpha}}}-\sum\nolimits_{i=1}^{K}{\frac{{{{\boldsymbol{\mu}}^{\rm{H}}}{\bf{A}}{{\boldsymbol{\psi}}_{i}}\left({{\delta_{i}}{{\boldsymbol{\psi}}_{i}}^{\rm{H}}{{\bf
B}{\boldsymbol{\mu}}}+{{\boldsymbol{\psi}}_{i}}^{\rm{H}}{\boldsymbol{\mu}}}\right)}}{{1+\left({s+t}\right){\delta_{i}}}}}}}}}{{st\prod\nolimits_{i=1}^{K}{\left({1+\left({s+t}\right){\delta_{i}}}\right)}}}}_{{F_{k}}\left(s,t\right)}dsdt}}$
$\displaystyle\triangleq f_{k}(\theta_{\tilde{k}},\theta_{k}).$ (61)

The above two-dimensional inverse Laplace transform can be obtained by means
of the Moorthy algorithm [35]. More specifically,
$f_{k}(\theta_{\tilde{k}},\theta_{k})$ is approximated by adopting a
trapezoidal rule of integration as

$f_{k}(\theta_{\tilde{k}},\theta_{k})\approx\frac{e^{c_{1}\theta_{\tilde{k}}+c_{2}\theta_{k}}}{4T^{2}}\left\\{{2{\rm{Re}}\left\\{\begin{array}[]{l}\sum\nolimits_{{l_{1}}=1}^{\infty}{\sum\nolimits_{{l_{2}}=1}^{\infty}{F_{k}^{{l_{1}},{l_{2}}}E_{\theta_{\tilde{k}},\theta_{k}}^{{l_{1}},{l_{2}}}}+{F_{k}^{{l_{1}},-{l_{2}}}E_{\theta_{\tilde{k}},\theta_{k}}^{{l_{1}},-{l_{2}}}}}\\\
+\sum\nolimits_{{l_{1}}=1}^{\infty}{F_{k}^{{l_{1}},0}E_{\theta_{\tilde{k}}}^{{l_{1}}}}+\sum\nolimits_{{l_{2}}=1}^{\infty}{F_{k}^{0,{l_{2}}}E_{\theta_{k}}^{{l_{2}}}}\end{array}\right\\}+F_{k}^{0,0}}\right\\},$
(62)

where the parameter $T$ determines the sampling period, $L$ is the truncation
order,

$F_{k}^{{l_{1}},{l_{2}}}={F_{k}}\left({{c_{1}}+{\rm
i}{l_{1}}{\pi}/T,{c_{2}}+{\rm i}{l_{2}}{\pi}/T}\right),$ (63)
$E_{\theta_{\tilde{k}},\theta_{k}}^{{l_{1}},{l_{2}}}={e^{{\rm{i}}{l_{1}}\pi\theta_{\tilde{k}}/{T}+{\rm{i}}{l_{2}}\pi\theta_{k}/{T}}}=E_{\theta_{\tilde{k}}}^{{l_{1}}}E_{\theta_{k}}^{{l_{2}}},$
(64)

As proved in [35], the discretization error $E_{r}$ can be controlled by
properly setting $c_{i}$ such that $c_{2}=-1/(2T)\ln((E_{r}-a)/(1-\xi))$,
$\xi=e^{-2Tc_{1}}<E_{r}$. According to the Moorthy algorithm, (62) can be
calculated by truncating the infinite series such that $l_{i}\in[0,L]$. The
tests in [35] suggest $L/\max(\theta_{\tilde{k}},\theta_{k})\in[0.5,0.8]$.
Furthermore, the Epsilon algorithm is invoked to extrapolate the remaining
terms above $L$ terms so as to accelerate the convergence and accuracy[36]. To
be specific, $2{\cal P}+1$ partial sums are used for the Epsilon algorithm,
and the experiments in [37] suggest ${\cal P}=2$.

Thereafter, by putting (III-B1) into (III-B1), we can obtain the conditional
outage probability ${p_{k|\mathcal{H}}}$.

####  III-B2 Approximate Expression of ${p_{k|{\cal H}}}$

Unfortunately, the two-dimensional numerical inversion of Laplace transform
entails considerable computation burden, which precludes to the real-time
optimal system design. In what follows, an approximate expression of
${p_{k|\mathcal{H}}}$ is derived to combat this dilemma. By overlooking the
dependence between ${\rm{SIN}}{{\rm{R}}_{k\to\tilde{k}}}$ and
${\rm{SIN}}{{\rm{R}}_{k}}$, ${p_{k|{\cal H}}}$ is approximated by using
(III-B1) as

$\displaystyle{p_{k|{\cal H}}}\approx
1-\underbrace{\Pr\left({{\rm{SIN}}{{\rm{R}}_{k\to\tilde{k}}}\geq{2^{{R_{\tilde{k}}}}}-1}\right)}_{{q_{k\to\tilde{k}|{\cal
H}}}}\underbrace{\Pr\left({{\rm{SIN}}{{\rm{R}}_{k}}\geq{2^{{R_{k}}}}-1}\right)}_{{q_{k\to
k|{\cal H}}}}.$ (65)

Likewise, ${q_{k\to\tilde{k}|{\cal H}}}$ and ${q_{k\to\tilde{k}|{\cal H}}}$
can be obtained by following the same steps as in Section III-A1. To be
specific, by introducing two vectors
${{\boldsymbol{\nu}}_{1}}={\left({{\mu_{1}},\cdots,{\mu_{k-1}},{\beta_{k}}^{2}{\mu_{k}},{\mu_{k+1}},\cdots,{\mu_{K}}}\right)^{\rm{T}}}$
and
${{\boldsymbol{\nu}}_{2}}={\left({{\mu_{1}},\cdots,{\mu_{k-1}},0,{\mu_{k+1}},\cdots,{\mu_{K}}}\right)^{\rm{T}}}$,
${q_{k\to\tilde{k}|{\cal H}}}$ and ${q_{k\to k|{\cal H}}}$ can be respectively
simplified as

$\displaystyle{q_{k\to\tilde{k}|{\cal
H}}}=\Pr\left({{{\left({{\boldsymbol{\chi}}+{{\boldsymbol{\nu}}_{1}}}\right)}^{\rm{H}}}\left({{\boldsymbol{\chi}}+{{\boldsymbol{\nu}}_{1}}}\right)\geq{\theta_{\tilde{k}}}-{\beta_{k}}^{2}{\beta_{\tilde{k}}}^{2}{{\left|{{\mu_{k}}}\right|}^{2}}}-\frac{{{I_{k}}}}{{P\ell\left({{d_{k}}}\right)}}\right),$
(66) $\displaystyle{q_{k\to k|{\cal H}}}=$
$\displaystyle\Pr\left({{{\left({{\boldsymbol{\chi}}+{{\boldsymbol{\nu}}_{2}}}\right)}^{\rm{H}}}\left({{\boldsymbol{\chi}}+{{\boldsymbol{\nu}}_{2}}}\right)\geq{\theta_{k}}}-\frac{{{I_{k}}}}{{P\ell\left({{d_{k}}}\right)}}\right),$
(67)

where (66) is derived in analogous to (18). Apparently, both (66) and (67) can
be computed as (26) by replacing
$({\tau_{\tilde{k}}},{\boldsymbol{\tilde{\nu}}},{\bf{\tilde{\Sigma}}},{\omega_{\tilde{k}}},{d_{\tilde{k}}})$
with
$({\theta_{\tilde{k}}}-{\beta_{k}}^{2}{\beta_{\tilde{k}}}^{2}{{\left|{{\mu_{k}}}\right|}^{2}},{\boldsymbol{\nu}}_{1},{\bf{\Sigma}},{\omega_{k}},{d_{k}})$
and
$({\theta_{k}},{\boldsymbol{\nu}}_{2},{\bf{\Sigma}},{\omega_{k}},{d_{k}})$,
respectively. The explicit expressions of ${q_{k\to\tilde{k}|{\cal H}}}$ and
${q_{k\to k|{\cal H}}}$ are omitted here to conserve space. Moreover, the
approximate expression of (65) will underestimate the outage performance due
to the neglect of positive correlation between successful events of
subtracting the far user’s message and decoding its own message at user $k$.
This assertion is further justified by the simulation results in Section VI.
Hence, (65) in fact serves as an upper bound of ${p_{k|{\cal H}}}$. With (65),
a robust optimal design can then be enabled.

####  III-B3 Average Outage Probability

Moreover, the average outage probability by taking the expectation with
respect to the distance $d_{k}$ is given by

${p_{k}}={{\mathbb{E}}_{{d_{k}}}}\left\\{{{p_{k|{\mathcal{H}}}}}\right\\}=1-{{\mathbb{E}}_{{d_{k}}}}\left\\{{{q_{k|{\mathcal{H}}}}}\right\\},$
(68)

where ${{\mathbb{E}}_{{d_{k}}}}\left\\{{{q_{k|{\mathcal{H}}}}}\right\\}$ is
obtained by

$\displaystyle{{\mathbb{E}}_{{d_{k}}}}\left\\{{{q_{k|{\mathcal{H}}}}}\right\\}=$
$\displaystyle\frac{1}{{{{\left({2\pi{\rm{i}}}\right)}^{2}}}}\int\limits_{{c_{1}}-{\rm{i}}\infty}^{{c_{1}}+{\rm{i}}\infty}{\int\limits_{{c_{2}}-{\rm{i}}\infty}^{{c_{2}}+{\rm{i}}\infty}{{e^{s{{\bar{\theta}}_{\tilde{k}}}+t{{\bar{\theta}}_{k}}}}\underbrace{\frac{{{e^{-\sum\nolimits_{i=1}^{K}{\frac{{{\boldsymbol{\mu}^{\rm{H}}}{\bf{A}}{\boldsymbol{\psi}_{i}}\left({{\delta_{i}}{\boldsymbol{\psi}_{i}}^{\rm{H}}{\bf{B}}\boldsymbol{\mu}+{\boldsymbol{\psi}_{i}}^{\rm{H}}\boldsymbol{\mu}}\right)}}{{1+\left({s+t}\right){\delta_{i}}}}}}}\overbrace{\mathbb{E}\left\\{{{e^{-\frac{{{{\sigma_{k}^{2}}}}}{{P\ell\left({{d_{k}}}\right)}}\left({s+t}\right)-\pi{\lambda_{b}}{\omega_{k}}{d_{k}}^{2}{{(s+t)}^{\frac{2}{\alpha}}}}}}\right\\}}^{{\varphi_{k}}\left({s+t}\right)}}}{{st\prod\nolimits_{i=1}^{K}{\left({1+\left({s+t}\right){\delta_{i}}}\right)}}}}_{{g_{k}}\left({s,t}\right)}dsdt}}.$
(69)

where
${{\bar{\theta}}_{k}}={\left|{{\mu_{k}}}\right|^{2}}{\beta_{k}}^{2}/({2^{{R_{k}}}}-1)$
and
${{\bar{\theta}}_{\tilde{k}}}={\left|{{\mu_{k}}}\right|^{2}}{\beta_{\tilde{k}}}^{2}/({2^{{R_{\tilde{k}}}}}-1)$.
By using the Moorthy algorithm, ${p_{k}}$ can be consequently derived as

$\displaystyle{p_{k}}$ $\displaystyle\approx
1-\frac{e^{c_{1}\bar{\theta}_{\tilde{k}}+c_{2}\bar{\theta}_{k}}}{4T^{2}}\left\\{{2{\rm{Re}}\left\\{{\begin{array}[]{*{20}{l}}{\sum\nolimits_{{l_{1}}=1}^{L}{\sum\nolimits_{{l_{2}}=1}^{L}{g_{k}^{{l_{1}},{l_{2}}}E_{\bar{\theta}_{\tilde{k}},\bar{\theta}_{k}}^{{l_{1}},{l_{2}}}}}+{g_{k}^{{l_{1}},-{l_{2}}}E_{\bar{\theta}_{\tilde{k}},\bar{\theta}_{k}}^{{l_{1}},-{l_{2}}}}}\\\
{+\sum\nolimits_{{l_{1}}=1}^{L}{g_{k}^{{l_{1}},0}E_{\bar{\theta}_{\tilde{k}}}^{{l_{1}}}}+\sum\nolimits_{{l_{2}}=1}^{L}{g_{k}^{0,{l_{2}}}E_{\bar{\theta}_{k}}^{{l_{2}}}}}\end{array}}\right\\}+g_{k}^{0,0}}\right\\},$
(72)

where $g_{k}^{{l_{1}},{l_{2}}}={g_{k}}\left({{c_{1}}+{\rm
i}{l_{1}}{\pi}/T,{c_{2}}+{\rm i}{l_{2}}{\pi}/T}\right)$.

Similarly to Theorem 1 and Theorem 2, ${\varphi_{k}}\left({s+t}\right)$ are
obtained by considering the random and the distance-based policies.

######  Theorem 3.

Under the random NOMA grouping policy, ${\varphi_{k}}\left(s+t\right)$ is
given by

$\displaystyle{\varphi_{k}}\left(s+t\right)=\frac{{2c}}{{2c+{\omega_{k}}{(s+t)^{\frac{2}{\alpha}}}}}H_{1,1}^{1,1}\left({\left.{\frac{{{{\sigma_{k}^{2}}}}}{P}(s+t){{\left({\pi{\lambda_{b}}\left({2c+{\omega_{k}}{(s+t)^{\frac{2}{\alpha}}}}\right)}\right)}^{-\frac{\alpha}{2}}}}\right|\begin{array}[]{*{20}{c}}{\left({0,\frac{\alpha}{2}}\right)}\\\
{\left({0,1}\right)}\end{array}}\right),$ (75)

In the interference-limited regime, i.e., $\sigma^{2}/P\to 0$,
${\varphi_{k}}\left(s\right)$ can be simplified as

${\varphi_{k}}\left(s+t\right)=\frac{{2c}}{{2c+{\omega_{k}}{(s+t)^{\frac{2}{\alpha}}}}}.$
(76)

###### Proof.

Similarly to (38), the PDF of ${{d_{k}}}$ is given by

${f_{{d_{k}}}}\left(x\right)=2\left({1-{F_{d}}(x)}\right){f_{d}}(x),$ (77)

By using (77), (75) can be obtained accordingly. ∎

######  Theorem 4.

Under the distance-based NOMA grouping policy, ${\varphi_{k}}\left(s+t\right)$
can be obtained similarly to (2) by directly replacing $\tilde{k}$,
$\sigma_{\tilde{k}}^{2}$, $r_{\tilde{k}}$ and $s$ with $k$, $\sigma_{k}^{2}$,
$r_{k}$ and $(s+t)$, respectively. Thereon, in the interference-limited
regime, i.e., $\sigma^{2}/P\to 0$, ${\varphi_{k}}\left(s\right)$ collapses to

$\displaystyle{\varphi_{k}}\left(s+t\right)=c{r_{k}}{{2K}\choose{r_{k}}}\sum\limits_{l=0}^{{r_{k}}-1}{\frac{{{{\left({-1}\right)}^{l}}{{{r_{k}}-1}\choose{l}}}}{{\left({2K-{r_{k}}+l+1}\right)c+{\omega_{k}}{(s+t)^{\frac{2}{\alpha}}}}}},$
(78)

where $r_{k}$ denotes the ranking order of user $k$ and $r_{k}<r_{\tilde{k}}$.

###### Proof.

The proof is parallel to that of Theorem 2 and hence omitted. ∎

Clearly, the average outage probability ${p_{k}}$ is independent of the
intensity of BSs because ${\varphi_{k}}\left(s\right)$ is irrelevant to
${{\lambda_{b}}}$.

##  IV Asymptotic Outage Analysis

Unfortunately, the exact outage expressions are too cumbersome to gain extra
helpful insights. To address this issue, this section seeks to investigate the
asymptotic behavior of the outage probability as the quality of CSI improves,
i.e., ${\mathcal{K}}_{zk}\to\infty$ or $\sigma_{h}^{2}\to 0$. For simplicity
and tractability, we specialize our analysis to the case of no inter-cell
interference, i.e., $\lambda_{b}=0$. Similarly to [24], the tool of Chernoff
bound is resorted to offer upper bounds for the outage probabilities. The
asymptotic outage probabilities of the far and near users are individually
discussed by splitting into two subsections.

###  IV-A Asymptotic Outage Probability of the Far User $\tilde{k}$

If $\lambda_{b}=0$, the conditional outage probability of user $\tilde{k}$ is

$\displaystyle{p_{\tilde{k}|\mathcal{H}}}$
$\displaystyle=\Pr\left\\{{{{\left({{\boldsymbol{\tilde{\chi}}}+{\boldsymbol{\tilde{\nu}}}}\right)}^{\rm{H}}}\left({{\boldsymbol{\tilde{\chi}}}+{\boldsymbol{\tilde{\nu}}}}\right)>{\tau_{\tilde{k}}}}\right\\},$
(79)

By using Chernoff bound, ${p_{\tilde{k}|\mathcal{H}}}$ is upper bounded by

${p_{\tilde{k}|{\cal
H}}}\leq{e^{-{\tau_{\tilde{k}}}\tilde{s}}}\mathbb{E}\left\\{{{e^{\tilde{s}{{\left({{\boldsymbol{\tilde{\chi}}}+{\boldsymbol{\tilde{\nu}}}}\right)}^{\rm{H}}}\left({{\boldsymbol{\tilde{\chi}}}+{\boldsymbol{\tilde{\nu}}}}\right)}}}\right\\},$
(80)

for any $\tilde{s}>0$. With the joint PDF of ${\boldsymbol{\tilde{\chi}}}$ and
the important integral (III-A1), (80) can be rewritten as

$\displaystyle{p_{\tilde{k}|{\cal
H}}}\leq\frac{{{e^{-{\tau_{\tilde{k}}}\tilde{s}}}}}{{{\pi^{K}}\det\left({{\bf{\tilde{\Sigma}}}}\right)}}\int\limits_{{\mathbb{C}^{K}}}{{e^{\tilde{s}{{\left({{\bf{x}}+{\boldsymbol{\tilde{\nu}}}}\right)}^{\rm{H}}}\left({{\bf{x}}+{\boldsymbol{\tilde{\nu}}}}\right)-{{\bf{x}}^{\rm{H}}}{{{\boldsymbol{\tilde{\Sigma}}}}^{-1}}{\bf{x}}}}}d{\bf{x}}$
$\displaystyle=\frac{{{e^{-{\tau_{\tilde{k}}}\tilde{s}}}}}{{\det\left({{\bf{I}}-\tilde{s}{\bf{\tilde{\Sigma}}}}\right)}}{e^{{{{\bf{\tilde{\nu}}}}^{\rm{H}}}{{\left({\frac{1}{\tilde{s}}{\bf{I}}-{\bf{\tilde{\Sigma}}}}\right)}^{-1}}{\bf{\tilde{\nu}}}}}$
$\displaystyle={e^{-\left({{\tau_{\tilde{k}}}\tilde{s}+\sum\limits_{i=1}^{K}{\ln\left({1-\tilde{s}{\delta_{\tilde{i}}}}\right)}+\sum\limits_{i=1}^{K}{\frac{{\tilde{s}{{\left|{{\zeta_{\tilde{i}}}}\right|}^{2}}}}{{\tilde{s}{\delta_{\tilde{i}}}-1}}}}\right)}},$
(81)

where $\tilde{s}\in(0,\min\\{{{1/\delta_{\tilde{i}}}},i\in[1,K]\\})$. To
circumvent the effect of ${\sigma_{h}^{2}}$ on $\tilde{s}$, we define
${\delta_{\tilde{i}}}={\sigma_{h}^{2}}{\upsilon_{\tilde{i}}}$, where
${\delta_{\tilde{i}}}$ are the eigenvalues of
${\boldsymbol{\tilde{\Sigma}}}={\sigma_{h}^{2}}{{\bf{u}}_{\tilde{k}}}^{\rm{H}}{{\bf{R}}_{r\tilde{k}}}{{\bf{u}}_{\tilde{k}}}{\left({{{\bf{V}}^{\rm{H}}}{{\bf{R}}_{t\tilde{k}}}{\bf{V}}}\right)^{\rm{T}}}$.
Accordingly, $\\{\upsilon_{\tilde{1}},\cdots,\upsilon_{\tilde{K}}\\}$ are the
eigenvalues of
${{\bf{u}}_{\tilde{k}}}^{\rm{H}}{{\bf{R}}_{r\tilde{k}}}{{\bf{u}}_{\tilde{k}}}{\left({{{\bf{V}}^{\rm{H}}}{{\bf{R}}_{t\tilde{k}}}{\bf{V}}}\right)^{\rm{T}}}$.
Hereby, (IV-A) can be expressed as

$\displaystyle{p_{\tilde{k}|{\cal H}}}$
$\displaystyle\leq{e^{-\frac{{\bar{s}}}{{\sigma_{h}^{2}}}{\left({{\tau_{\tilde{k}}}+\sum\limits_{i=1}^{K}{\frac{{{{\left|{{\zeta_{\tilde{i}}}}\right|}^{2}}}}{{\bar{s}{\upsilon_{\tilde{i}}}-1}}}}\right)}}}{{e^{-\sum\limits_{i=1}^{K}{\ln\left({1-\bar{s}{\upsilon_{\tilde{i}}}}\right)}}}},$
(82)

where
$\bar{s}={\sigma_{h}^{2}}\tilde{s}\in(0,\min\\{1/{\upsilon_{\tilde{i}}},i\in[1,K]\\})$.
Clearly from (82), the following theorem concerning the asymptotic behavior of
${p_{\tilde{k}|{\cal H}}}$ can be found.

######  Theorem 5.

As the channel uncertainty vanishes, i.e., ${\mathcal{K}}_{zk}\to\infty$ or
${{\sigma_{h}^{2}}}\to 0$, the outage probability ${p_{\tilde{k}|{\cal H}}}$
decays to zero if the target transmission rate satisfies the following
sufficient condition

${R_{\tilde{k}}}<{\log_{2}}\left({1+\frac{{{\beta_{\tilde{k}}}^{2}{{\left|{{\mu_{\tilde{k}}}}\right|}^{2}}}}{{\sum\nolimits_{i=1}^{K}{{{\left|{{\zeta_{\tilde{i}}}}\right|}^{2}}}+{\beta_{k}}^{2}{\beta_{\tilde{k}}}^{2}{{\left|{{\mu_{\tilde{k}}}}\right|}^{2}}}+\frac{{\sigma_{\tilde{k}}^{2}}}{{P\ell\left({{d_{\tilde{k}}}}\right)}}}}\right).$
(83)

###### Proof.

As ${{\sigma_{h}^{2}}}\to 0$, the right-hand side (RHS) of the inequality (82)
approaches to zero only if
${{\tau_{\tilde{k}}}+\sum\nolimits_{i=1}^{K}{{{{{\left|{{\zeta_{\tilde{i}}}}\right|}^{2}}}}/{({\bar{s}{\upsilon_{\tilde{i}}}-1})}}}>0$.
By using the definition of ${\tau_{\tilde{k}}}$ together with some
rearrangements, setting $\bar{s}=0$ yields the upper bound of
${R_{\tilde{k}}}$. ∎

It should be highlighted that Theorem 5 differs from the results in [24],
where only a single data stream is considered. Besides, a loose upper bound of
${R_{\tilde{k}}}$ is provided in [24]. Particularly, the simulation results in
Section VI justify the correctness of Theorem 5. However, it is worth
mentioning that the same property does not carry over to the case of
$\lambda_{b}\neq 0$, that is, an arbitrary low outage target cannot be
warranted for any nonzero target rates.

###  IV-B Asymptotic Outage Probability of the Near User $k$

From (III-B1), ${p_{k|{\cal H}}}$ is upper bounded by using inclusion-
exclusion principle as

$\displaystyle{p_{k|{\cal H}}}\leq$
$\displaystyle\Pr\left({{{\log}_{2}}\left({1+{\rm{SIN}}{{\rm{R}}_{k\to\tilde{k}}}}\right)<{R_{\tilde{k}}}}\right)+\Pr\left({{{\log}_{2}}\left({1+{\rm{SIN}}{{\rm{R}}_{k}}}\right)<{R_{k}}}\right).$
(84)

By setting $\lambda_{b}=0$ together with the definitions
${{\boldsymbol{\nu}}_{1}}=({\mu_{1}},\cdots,{\mu_{k-1}},{\beta_{k}}^{2}{\mu_{k}},{\mu_{k+1}},\cdots,{\mu_{K}})^{\rm{T}}$
${{\boldsymbol{\nu}}_{2}}=({\mu_{1}},\cdots,{\mu_{k-1}},0,{\mu_{k+1}},\cdots,{\mu_{K}})^{\rm{T}}$,
(84) can be rewritten as

$\displaystyle{p_{k|{\cal H}}}\leq$
$\displaystyle\Pr\left({{{\left({{\boldsymbol{\chi}}+{{\boldsymbol{\nu}}_{1}}}\right)}^{\rm{H}}}\left({{\boldsymbol{\chi}}+{{\boldsymbol{\nu}}_{1}}}\right)>{\theta_{\tilde{k}}}-{\beta_{k}}^{2}{\beta_{\tilde{k}}}^{2}{{\left|{{\mu_{k}}}\right|}^{2}}}\right)+\Pr\left({{{\left({{\boldsymbol{\chi}}+{{\boldsymbol{\nu}}_{2}}}\right)}^{\rm{H}}}\left({{\boldsymbol{\chi}}+{{\boldsymbol{\nu}}_{2}}}\right)>{\theta_{k}}}\right),$
(85)

where the first term on the RHS of the inequality is obtained similarly to
(18). In analogous to (80) and (IV-A), applying Chernoff bound to (85) gives
rise to

$\displaystyle{p_{k|{\cal H}}}<$
$\displaystyle{e^{-\left({\left({{\theta_{\tilde{k}}}-{\beta_{k}}^{2}{\beta_{\tilde{k}}}^{2}{{\left|{{\mu_{k}}}\right|}^{2}}}\right){s}+\sum\limits_{i=1}^{K}{\ln\left({1-{s}{\delta_{i}}}\right)}+\sum\limits_{i=1}^{K}{\frac{{{s}{{\left|{{{\boldsymbol{\psi}}_{i}}{{\boldsymbol{\nu}}_{1}}}\right|}^{2}}}}{{s{\delta_{i}}-1}}}}\right)}}+{e^{-\left({{\theta_{k}}{s}+\sum\limits_{i=1}^{K}{\ln\left({1-{s}{\delta_{i}}}\right)}+\sum\limits_{i=1}^{K}{\frac{{{s}{{\left|{{{\boldsymbol{\psi}}_{i}}{{\boldsymbol{\nu}}_{2}}}\right|}^{2}}}}{{s{\delta_{i}}-1}}}}\right)}},$
(86)

where $s\in(0,\min\\{{{1/\delta_{i}}},i\in[1,K]\\})$.

Similarly to (82), we denote the eigenvalues of
${{\bf{u}}_{k}}^{\rm{H}}{{\bf{R}}_{rk}}{{\bf{u}}_{k}}{\left({{{\bf{V}}^{\rm{H}}}{{\bf{R}}_{tk}}{\bf{V}}}\right)^{\rm{T}}}$
by $\\{\upsilon_{1},\cdots,\upsilon_{K}\\}$. Thus, we have
${\delta_{i}}={\sigma_{h}^{2}}{\upsilon_{i}}$, where ${\delta_{i}}$ are the
eigenvalues of ${\bf{\Sigma}}$. (86) can thus be rewritten as

$\displaystyle{p_{k|{\cal
H}}}<{e^{-\sum\limits_{i=1}^{K}{\ln\left({1-\hat{s}{\upsilon_{i}}}\right)}}}\left({e^{-\frac{{\hat{s}}}{{\sigma_{h}^{2}}}\left({{\theta_{\tilde{k}}}-{\beta_{k}}^{2}{\beta_{\tilde{k}}}^{2}{{\left|{{\mu_{k}}}\right|}^{2}}+\sum\limits_{i=1}^{K}{\frac{{{{\left|{{{\boldsymbol{\psi}}_{i}}^{\rm
H}{{\boldsymbol{\nu}}_{1}}}\right|}^{2}}}}{{\hat{s}{\upsilon_{i}}-1}}}}\right)}}\right.\left.+{e^{-\frac{{\hat{s}}}{{\sigma_{h}^{2}}}\left({{\theta_{k}}+\sum\limits_{i=1}^{K}{\frac{{{{\left|{{{\boldsymbol{\psi}}_{i}}^{\rm
H}{{\boldsymbol{\nu}}_{2}}}\right|}^{2}}}}{{\hat{s}{\upsilon_{i}}-1}}}}\right)}}\right),$
(87)

where $\hat{s}={\sigma_{h}^{2}}s\in(0,\min\\{1/{\upsilon_{i}},i\in[1,K]\\})$.
On the basis of (87), the following theorem regarding the asymptotic behavior
of ${p_{k|{\cal H}}}$ can be proved.

######  Theorem 6.

As the channel uncertainty vanishes, i.e., ${\mathcal{K}}_{zk}\to\infty$ or
${{\sigma_{h}^{2}}}\to 0$, the outage probability ${p_{k|{\cal H}}}$ converges
to zero if the target transmission rates satisfy the following sufficient
conditions

${R_{\tilde{k}}}<{\log_{2}}\left({1+\frac{{{\beta_{\tilde{k}}}^{2}{{\left|{{\mu_{k}}}\right|}^{2}}}}{{\sum\nolimits_{i=1}^{K}{{{\left|{{\boldsymbol{\psi}_{i}}^{\rm
H}{\boldsymbol{\nu}_{1}}}\right|}^{2}}}+{\beta_{k}}^{2}{\beta_{\tilde{k}}}^{2}{{\left|{{\mu_{k}}}\right|}^{2}}}+\frac{{\sigma_{k}^{2}}}{{P\ell\left({{d_{k}}}\right)}}}}\right),$
(88)
${R_{k}}<{\log_{2}}\left({1+\frac{{{{\left|{{\mu_{k}}}\right|}^{2}}{\beta_{k}}^{2}}}{{\sum\nolimits_{i=1}^{K}{{{\left|{{{\boldsymbol{\psi}}_{i}}^{\rm
H}{{\boldsymbol{\nu}}_{2}}}\right|}^{2}}}+\frac{{\sigma_{k}^{2}}}{{P\ell\left({{d_{k}}}\right)}}}}}\right).$
(89)

###### Proof.

As ${{\sigma_{h}^{2}}}\to 0$, the RHS of the inequality (87) approaches to
zero only if
${\theta_{\tilde{k}}}-{\beta_{k}}^{2}{\beta_{\tilde{k}}}^{2}{\left|{{\mu_{k}}}\right|^{2}}+\sum\nolimits_{i=1}^{K}{{{{{\left|{{{\bf{\psi}}_{i}}{{\bf{\nu}}_{1}}}\right|}^{2}}}}/{({\hat{s}{\upsilon_{i}}-1})}}>0$
and
${\theta_{k}}+\sum\nolimits_{i=1}^{K}{{{{{\left|{{{\bf{\psi}}_{i}}{{\bf{\nu}}_{2}}}\right|}^{2}}}}/{({\hat{s}{\upsilon_{i}}-1})}}>0$.
By using the definitions of ${\theta_{\tilde{k}}}$ and ${\theta_{k}}$, the
upper bounds of ${R_{\tilde{k}}}$ and ${R_{k}}$ can be obtained by setting
$\hat{s}=0$. ∎

##  V Optimal System Design

By making use of the foregoing analytical results, the parameters of MIMO-NOMA
enhanced SCNs can be properly configured to adapt to the imperfect fading
channels while ensuring the quality of service, including the precoding matrix
(i.e., $\bf V$), receiver filters (i.e., ${\bf u}_{k},{\bf
u}_{\tilde{k}},k\in[1,K]$), power allocation coefficients (i.e., $\beta_{k}$,
$\beta_{\tilde{k}},k\in[1,K]$) and the transmission rates (i.e., $R_{k}$,
$R_{\tilde{k}},k\in[1,K]$). To exemplify this, we confine our attention to
maximize the long-term goodput that measures the number of successfully
conveyed information bits per transmission[38, 39, 40]. Given the known
channel information $\mathcal{H}$, the conditional goodput of the MIMO-NOMA
system is expressed as

$\displaystyle{\mathcal{T}_{g|{\mathcal{H}}}}$
$\displaystyle=\sum\limits_{k=1}^{K}{{R_{k}}({1-{p_{k|\mathcal{H}}}})+{R_{\tilde{k}}}({1-{p_{\tilde{k}|\mathcal{H}}}})}.$
(90)

Thereon, the maximization of goodput while guaranteeing low outage constraints
is posed by

$\begin{array}[]{*{20}{c
l}}{\mathop{\max}\limits_{{\bf{V}},\\{{{\bf{u}}_{k}},{{\bf{u}}_{\tilde{k}}},{R_{k}},{R_{\tilde{k}}}:k\in[1,K]\\}}}&{{\mathcal{T}_{g|{\cal
H}}}}\\\ {{\rm{s}}{\rm{.t}}{\rm{.}}}&{{p_{k|{\cal
H}}}\leq\varepsilon,{p_{\tilde{k}|{\cal H}}}\leq\varepsilon,k\in[1,K]}\\\
{}\hfil&{\left\|{{{\bf{v}}_{k}}}\right\|=1,k\in[1,K]}\\\
{}\hfil&{{{\beta_{k}}^{2}}\left({{2^{{R_{\tilde{k}}}}}-1}\right)<1,k\in[1,K]}\end{array},$
(91)

where $\varepsilon$ denotes the maximum endurable outage probability. In
general, the power allocation coefficients $\beta_{k}$ and transmission rates
$R_{k}$, $R_{\tilde{k}}$ cannot be jointly optimized in (91) in order to not
violate the intention of NOMA principle. This is because as disclosed by [41],
the joint optimization of $\beta_{k}$ and $R_{k}$&$R_{\tilde{k}}$ would behave
like waterfilling algorithm if the user fairness is disregarded, that is, all
the power would be allocated to the users under benign channel conditions and
no information bits would be delivered to the users experienced bad channel
conditions. In light of this observation, we fix the power allocation
coefficient to maintain the user fairness in the sequel.

Unfortunately, due to the complex outage expressions, it is virtually
impossible to get the globally optimal solution to (91). To tackle this
difficulty, the concept of the signal/interference alignment proposed by [13]
is applied to choose an appropriate ${\bf{V}}$, ${{\bf{u}}_{k}}$ and
${{\bf{u}}_{\tilde{k}}}$ with low implementation overhead111Intrinsically, the
signal alignment can be regarded as a linear decorrelator [42, Sec. 8.3.1],
which is a combination of the projection operation and matched filter. The
linear decorrelator is frequently used to suppress the inter-stream
interference.. In particular, in order to substantially suppress the inter-
pair interference, it is necessary to impose the constraints on
${\mu_{i}}={\mu_{\tilde{i}}}=0,i\neq k$. This leads to

$\left({\begin{array}[]{*{20}{c}}{{{\bf{u}}_{k}}^{\rm{H}}{{{\bf{\hat{H}}}}_{zk}}}\\\
{{{\bf{u}}_{\tilde{k}}}^{\rm{H}}{{{\bf{\hat{H}}}}_{z\tilde{k}}}}\end{array}}\right){{\bf{v}}_{i}}={\bf
0},\,i\neq k,k\in[1,K],$ (92)

In order to capitalize on the signal alignment, we apply the following
decomposition ${\bf V}={\bf L}{\bf P}$, where ${\bf L}\in{\mathbb{C}}^{M\times
K}$, ${\bf L}^{\rm H}{\bf L}={\bf I}_{K}$ and ${\bf P}\in{\mathbb{C}}^{K\times
K}$. In some sense, the involvement of ${\bf L}$ is equivalent to choose $K$
transmit antennas among $M$ ones. Hence, a similar algorithm as suggested in
[13] can be developed to obtain ${\bf L}$, as shown in Algorithm 1.
Furthermore, by defining
${{\bf{g}}_{k}}^{\rm{H}}={{\bf{u}}_{k}}^{\rm{H}}{{{\bf{\hat{H}}}}_{zk}}{\bf{L}}$
and
${{\bf{g}}_{\tilde{k}}}^{\rm{H}}={{\bf{u}}_{\tilde{k}}}^{\rm{H}}{{{\bf{\hat{H}}}}_{z\tilde{k}}}{\bf{L}}$,
we have
${\left({{{\bf{g}}_{k}},{{\bf{g}}_{\tilde{k}}}}\right)^{\rm{H}}}{{\bf{p}}_{i}}={\bf{0}}$,
where ${\bf P}=({{\bf{p}}_{1}},\cdots,{{\bf{p}}_{K}})$. To ensure the
existence of ${{\bf{p}}_{i}}$, the prerequisite of undertaking the signal
alignment is ${{\bf{g}}_{k}}={{\bf{g}}_{\tilde{k}}}$, which ensures the
existence of both ${{\bf{u}}_{k}}$ and ${{\bf{u}}_{\tilde{k}}}$. As proved by
[13], ${{{\bf{u}}_{k}}}$ and ${{{\bf{u}}_{\tilde{k}}}}$ are thus given by

$\left({\begin{array}[]{*{20}{c}}{{{\bf{u}}_{k}}}\\\
{{{\bf{u}}_{\tilde{k}}}}\end{array}}\right)={\bf U}_{k}{\bf z}_{k},$ (93)

where ${\bf U}_{k}\in{\mathbb{C}}^{2N\times(2N-K)}$ consists of the $(2N-K)$
right singular vectors of
$({{{({{{{\bf{\hat{H}}}}_{zk}}{\bf{L}}})}^{\rm{H}}},-{{({{{{\bf{\hat{H}}}}_{z\tilde{k}}}{\bf{L}}})}^{\rm{H}}}})$
corresponding to its zero singular values, ${\bf
z}_{k}\in{\mathbb{C}}^{(2N-K)\times 1}$ and $|{\bf z}_{k}|^{2}=1$. In [13], an
algorithm with low complexity has been provided to properly choose ${\bf
z}_{k}$. After determining ${{{\bf{u}}_{k}}}$ and ${{{\bf{u}}_{\tilde{k}}}}$,
${{\bf{g}}_{k}}$ and ${{\bf{g}}_{\tilde{k}}}$ are obtained. Accordingly, ${\bf
P}$ is then given by ${\bf{P}}={{\bf{G}}^{-{\rm{H}}}}{\bf{D}}$, where
${\bf{G}}=({\bf g}_{1},\cdots,{\bf g}_{K})$ and ${\bf{D}}={\rm
diag}(1/\sqrt{({{\bf{G}}^{-1}}{{\bf{G}}^{-{\rm{H}}}})_{1,1}},\cdots,1/\sqrt{({{\bf{G}}^{-1}}{{\bf{G}}^{-{\rm{H}}}})_{K,K}})$
ensures $\left\|{{{\bf{v}}_{k}}}\right\|=\left\|{{{\bf{p}}_{k}}}\right\|=1$.
Therefore, ${\bf V}$ can be calculated as

${\bf V}={\bf L}{{\bf{G}}^{-{\rm{H}}}}{\bf{D}}.$ (94)

An algorithm for the design of ${\bf V}$ is listed in Algorithm 1.

Algorithm 1 The selection of the precoding matrix ${\bf{V}}$

1:for $i=1$ to $M!/(M-K)!$ do

2: ${\bf L}\leftarrow{\bf Y}_{i}$, ${\bf Y}_{i}$ is a matrix of size ${M\times
K}$ that contains exactly one entry of 1 in each column, at most one entry of
1 in each row and 0s elsewhere;

3: Find ${{\bf{u}}_{k}},{{\bf{u}}_{\tilde{k}}}$ by using [13, Algorithm 1];

4: Construct ${\bf{G}}=({\bf g}_{1},\cdots,{\bf g}_{K})^{\rm{H}}$, where
${{\bf{g}}_{k}}^{\rm{H}}={{\bf{u}}_{k}}^{\rm{H}}{{{\bf{\hat{H}}}}_{zk}}{\bf{L}}$
and $k\in[1,K]$;

5: Get the effective channel gain for each user pair,
$\gamma_{k}=1/({{\bf{G}}^{-1}}{{\bf{G}}^{-{\rm{H}}}})_{k,k}$;

6: Find the smallest effective channel gain, i.e.,
$\gamma_{\min,i}=\min\\{\gamma_{1},\cdots,\gamma_{K}\\}$;

7:end for

8:Find the index $i$ that maximizes the smallest effective channel gain, i.e.,
$i^{*}=\mathop{\arg}\limits_{i\in[1,M!/(M-K)!]}\max\gamma_{\min,i}$;

9:${\bf V}\leftarrow{\bf L}{{\bf{G}}^{-{\rm{H}}}}{\bf{D}}$, where ${\bf
L}={\bf Y}_{i^{*}}$.

After deciding ${\bf{V}}$, ${{\bf{u}}_{k}}$ and ${{\bf{u}}_{\tilde{k}}}$, the
original problem of (91) collapses to

$\begin{array}[]{*{20}{c
l}}{\mathop{\max}\limits_{\\{{R_{k}},{R_{\tilde{k}}}:k\in[1,K]\\}}}&{\sum\limits_{k=1}^{K}{{R_{k}}({1-{p_{k|\mathcal{H}}}})+{R_{\tilde{k}}}({1-{p_{\tilde{k}|\mathcal{H}}}})}}\\\
{{\rm{s}}{\rm{.t}}{\rm{.}}}&{{p_{k|{\cal
H}}}\leq\varepsilon,{p_{\tilde{k}|{\cal H}}}\leq\varepsilon,k\in[1,K]}\\\
{}\hfil&{{{\beta_{k}}^{2}}\left({{2^{{R_{\tilde{k}}}}}-1}\right)<1,k\in[1,K]}\end{array}.$
(95)

By noticing that the objective function of (95) is a sum of
${{{R_{k}}({1-{p_{k|\mathcal{H}}}})+{R_{\tilde{k}}}({1-{p_{\tilde{k}|\mathcal{H}}}})}}$
for all $k\in[1,K]$, ${p_{k|{\cal H}}}$ and ${p_{\tilde{k}|{\cal H}}}$ only
depend on ${R_{k}}$ and ${R_{\tilde{k}}}$ (independent of ${R_{i}}$ and
${R_{\tilde{i}}}$ for $i\neq k$), (95) can further be decoupled as $K$
independent subproblems, i.e.,

$\begin{array}[]{*{20}{c
l}}{\mathop{\max}\limits_{{R_{k}},{R_{\tilde{k}}}}}&{{{R_{k}}({1-{p_{k|\mathcal{H}}}})+{R_{\tilde{k}}}({1-{p_{\tilde{k}|\mathcal{H}}}})}}\\\
{{\rm{s}}{\rm{.t}}{\rm{.}}}&{{p_{k|{\cal
H}}}\leq\varepsilon,{p_{\tilde{k}|{\cal H}}}\leq\varepsilon}\\\
{}\hfil&{{{\beta_{k}}^{2}}\left({{2^{{R_{\tilde{k}}}}}-1}\right)<1}\end{array},k\in[1,K].$
(96)

(96) can be solved with many popular programming tools, such as sequential
convex programming (SCP), interior-point method, etc.222To examine the
computational complexity of (96), the interior-point method is taken as an
example. Since there are only two variables to be optimized in each
subproblem, the computational complexity of each subproblem is mainly affected
by the number of Newton iterations. According to the complexity analysis of
the interior-point method in [43], the total number of Newton iterations is
linearly proportional to
$\left\lceil{\frac{{\ln\left({m/\left({{t^{\left(0\right)}\epsilon}}\right)}\right)}}{{\ln\mu}}}\right\rceil$,
where $m$ is the number of inequality constraints, $\mu$ is the increasing
factor, $t^{(0)}$ is the initial value of the parameter that sets the accuracy
of the barrier approximation, $\epsilon$ specifies the final duality gap and
$\left\lceil\cdot\right\rceil$ denotes the ceil operation. Accordingly, the
computational complexity of solving the $K$ independent subproblems is
linearly proportional to
$K\left\lceil{\frac{{\ln\left({m/\left({{t^{\left(0\right)}\epsilon}}\right)}\right)}}{{\ln\mu}}}\right\rceil$.
In fact, ${\ln\left({m/\left({{t^{\left(0\right)}\epsilon}}\right)}\right)}$
refers to the logarithm of the ratio of the initial duality gap
$m/t^{\left(0\right)}$ to the final duality gap $\epsilon$. In practice, the
factor $\mu$ is usually assigned with a larger value around 2 to 100, which
leads to the total number of iterations on the order of a few tens. As a
consequence, the computational complexity of solving (96) is very low..

##  VI Numerical Analysis

In this section, the simulation results are presented for verifications and
discussions. Unless otherwise specified, the simulation parameters are
tabulated in Table I333The intention behind the linear relationship between
$R_{k}$ and $R_{\tilde{k}}$ is to plot 2-D figures for the purpose of clear
illustrations., where an exponential correlation profile is used to model the
covariance matrices ${\bf R}_{rk}=(\kappa^{|i-j|})_{1\leq i,j\leq N}$ and
${\bf R}_{tk}=(\kappa^{|i-j|})_{1\leq i,j\leq M}$[24]. For illustration, the
pairing policy in [12, 44] is considered for the distanced-based grouping
policy, where the farthest user is paired with the nearest user, and the rest
of users are paired in the same fashion. Henceforth, the ranking orders with
respect to the distance-based NOMA grouping policy are given by
$r_{1},r_{\tilde{1}},r_{2},r_{\tilde{2}}=1,4,2,3$. Besides, in the succeeding
examples, the technique of the signal alignment is applied to relieve the
intra-cell interference.

TABLE I: Simulation Parameters. Parameters | Values | Parameters | Values
---|---|---|---
The number of transmit antennas $M$ | 3 | Gaussian noise power $\sigma^{2}$ | -$99$ dBm
The number of receive antennas $N$ | 2 | Distances $d_{\tilde{k}}=2.5d_{k}$ | 50 m
The number of user pairs $K$ | 2 | Target rates $R_{k}=2R_{\tilde{k}}$ | 1 bps/Hz
Path loss exponent $\alpha$ | 3.5 | Coverage radius of simulated network | 5 km
Intensity of users $\lambda_{u}$ |  $2\times 10^{-4}$ users/m2 | Exponential correlation coefficient $\kappa$ | 0.9
Intensity of BSs $\lambda_{b}$ |  $10^{-5}$ BSs/m2 | Channel K factor $\mathcal{K}_{zk}$ | 20 dB
Transmit power at BS $P$ |  $20$ dBm | 2-norm of the estimated CSI ${\|{\hat{\bf H}_{zk}}\|_{F}^{2}}$ | $MN$
Interference power $\rho_{I}$ |  $15$ dBm | Power allocation factor $\beta_{k}^{2}$ | 0.3

###  VI-A Verifications

Fig. 1 examines the conditional outage probability $p_{k|\mathcal{H}}$ against
the transmission rate $R_{\tilde{k}}$. In Fig. 1, the approximate results of
$p_{k|\mathcal{H}}$ (labeled as “$k$=1-Approx.”) are obtained by using (65).
Clearly, the analytical results are in perfect agreement with the simulation
results. In addition, the conditional outage probability increases with the
transmission rates. This is due to the tradeoff between the system throughput
and reliability. Besides, a slight gap between the approximate and simulated
outage probabilities of the near users can be observed. Furthermore, it is not
hard to find that the exact outage probability is always upper bounded by the
approximate one. This is due to the fact that the strong correlation between
${\rm{SIN}}{{\rm{R}}_{k\to\tilde{k}}}$ and ${\rm{SIN}}{{\rm{R}}_{k}}$ renders
the near user having a high successful probability to reconstruct its own
message if it successfully subtracts the far user’s message. The same trends
can be observed in Fig. 2.

Figure 1: The conditional outage probability $p_{k|\mathcal{H}}$ versus the
transmission rate $R_{\tilde{k}}$ by setting parameters as
$R_{k}=2R_{\tilde{k}}$.

Fig. 2 plots the conditional outage probability $p_{k|\mathcal{H}}$ against
the channel K factor ${\cal K}_{zk}$. As expected, the analytical results
coincide with the simulation ones. Moreover, as ${\cal K}_{zk}$ increases, the
outage probability decays to an outage floor. It is no doubt that the outage
floor is determined by the co-channel interference from surrounding BSs.

Figure 2: The conditional outage probability $p_{k|\mathcal{H}}$ versus the
channel K factor ${\cal K}_{zk}$ by setting $R_{k}=2R_{\tilde{k}}=1$ bps/Hz.

Fig. 3 investigates the average outage probability $p_{k}$ against the
transmission rate $R_{\tilde{k}}$. In both Figs. 3 and 4, the exact and
approximate outage probabilities of the random grouping policy (labeled as
“Random-Exact” and “Random-Approx.”) for far (or near) NOMA user are
calculated by using (33) (or (75)) and (37) (or (76)), respectively. Moreover,
the simulated outage probability (labeled as “Random-Sim.”) is plotted for
verification. Similarly, the numerical results for the distance-based grouping
policy can be obtained according to Theorems 2 and 4, and the corresponding
exact, approximate and simulated outage probabilities are labeled as
“Distance-Exact”, “Distance-Approx.” and “Distance-Sim.”, respectively. From
Fig. 3, we observe a perfect match between the analytical and simulation
results, which corroborates the validity of our analysis. It can also be seen
that the distance-based grouping policy performs better than the random
grouping policy in terms of the outage probability of the near user, while the
impact on the outage probability of the far user behaves differently. This is
due to the fact that the near/far user under the distance-based grouping
policy is closer to/farther away from its associated BS than that under random
grouping policy. Nevertheless, it is noteworthy that the random grouping
policy is frequently used for benchmarking purpose, because it has been
intensively investigated in [45] and [46] that the random grouping policy
provides the lowest sum rate gain.

Figure 3: The average outage probability $p_{k}$ versus the transmission rate
$R_{\tilde{k}}$ by setting $R_{k}=2R_{\tilde{k}}$.

Fig. 4 depicts the average outage probability $p_{k}$ against the intensity of
BSs $\lambda_{b}$. Aside from the similar findings as Fig. 3, it is evident
that the average outage probability is irrelevant to the intensity of BSs
under high $\lambda_{b}$, i.e., interference-limited regime. This is
consistent with our analytical results. Nonetheless, the average outage
performance declines if $\lambda_{b}$ falls below a threshold (e.g., $10^{-6}$
BSs/m2). Indeed, such observations seem to be counterintuitive since the
inter-cell interference is alleviated as $\lambda_{b}$ decreases. However, the
decrease of $\lambda_{b}$ also signifies the expansion of the size of Voronoi
cell, on average, this dominant effect engenders severe path loss.

Figure 4: The average outage probability $p_{k}$ versus the transmission rate
$R_{\tilde{k}}$ by setting $R_{k}=2R_{\tilde{k}}=1$ bps/Hz.

###  VI-B Impacts of System Parameters

Fig. 5 illustrates the impact of the exponential correlation coefficient
$\kappa$ on the conditional outage probability $p_{k|\mathcal{H}}$. It is
worth mentioning that $\lambda_{b}=10^{-7}$ BSs/m2 is designated in this
example to weaken the impact of interference so as to better reflect the
influence of estimation error correlation. Clearly from Fig. 4,
$\lambda_{b}=10^{-7}$ BSs/m2 can be chosen to get rid of interference-limited
regime. It is shown that the correlation coefficient has a positive effect on
the average outage probability $p_{k|\mathcal{H}}$ under low channel K factor
${\cal K}_{zk}$. On the contrary, the correlation coefficient impairs the
outage performance under high ${\cal K}_{zk}$. Moreover, the conditional
outage probability of the far user, i.e., $p_{k|{\cal H}}$, changes very
slightly with respect to the correlation coefficient $\kappa$ under ${\cal
K}_{zk}=20$dB. This is due to the fact that the channel estimation error
diminishes as $\mathcal{K}_{zk}$ increases. As a consequence, the error
correlation under high $\mathcal{K}_{zk}$ has a negligible impact on the
outage probability compared to that of the interference, especially for the
far user with larger path loss.

Figure 5: The conditional outage probability $p_{k|\mathcal{H}}$ versus the
exponential correlation coefficient $\kappa$ with parameters
$R_{k}=2R_{\tilde{k}}=1$bps/Hz and $\lambda_{b}=10^{-7}$BSs/m2.

As shown in Fig. 6, the conditional outage probability $p_{k|\mathcal{H}}$ of
NOMA pair $k=1$ is plotted against the transmission rate $R_{\tilde{k}}$ under
different channel K factor ${\mathcal{K}}_{zk}$. With the increase of
${\mathcal{K}}_{zk}$, the outage curves become steeper. In particular, the
outage probability significantly increases to 1 for ${\mathcal{K}}_{zk}=50$ dB
as $R_{\tilde{k}}$ grows from 1.73 bps/Hz to 1.74 bps/Hz. However, this is not
beyond our expectation which can be elucidated as follows. As unveiled by the
asymptotic outage analysis in Section IV, the upper bound of the transmission
rate $R_{\tilde{k}}=1.737$ bps/Hz, which is in accordance with our
observations in Fig. 6. In contrast, it is found from Theorem 6 that $R_{1}$
and $R_{2}$ are upper bounded by 21.846 bps/Hz and 23.758 bps/Hz,
respectively. Apparently, the outage performance of NOMA users are
considerably restricted by the target rates of the far users, i.e.,
$R_{\tilde{k}}$.

Figure 6: The conditional outage probability $p_{k|\mathcal{H}}$ versus the
transmission rate $R_{\tilde{k}}$ by setting $R_{k}=2R_{\tilde{k}}$ and
$\lambda_{b}=0$ BSs/m2.

###  VI-C Optimal System Design

Fig. 7 compares the maximum achievable goodput of the proposed MIMO-NOMA
enhanced scheme with those of three benchmarking schemes, including MIMO-OMA
with precoding scheme [24], MIMO-OMA without precoding scheme and MIMO-NOMA
without precoding scheme [12, 21]. For comparative fairness, the coefficient
$\beta_{k}^{2}$ is also used for the allocation of the orthogonal resources
(time/bandwidth) in the OMA scheme. Additionally, it is worthwhile to note
that the approximate expression of ${p_{k|{\cal H}}}$ is used to solve (96) to
realize efficient and robust system design. It can be seen from Fig. 7 that
the proposed scheme performs the best particularly for a great difference
between $d_{k}$ and $d_{\tilde{k}}$. This is due to the fact that the
intrinsic nature of NOMA exploits the difference between channel gains to
boost the spectral efficiency. Moreover, the precoding schemes can achieve a
noticeable goodput gain over those without precoding. In addition, as the
channel K factor ${\cal K}_{zk}$ increases (till no estimation error), the
maximum achievable goodput increases up to an upper bound, which obviously
depends on the severity of co-channel interference, i.e., $\lambda_{b}$.

Figure 7: Performance comparison with the three baseline schemes in terms of
the goodput ${\cal T}_{k|\mathcal{H}}$ versus the channel K factor ${\cal
K}_{zk}$ with $\varepsilon=10^{-2}$.

##  VII Conclusions

The performance of MIMO-NOMA enhanced SCNs has been thoroughly investigated by
taking into account the assumption of imperfect CSI. The channel estimation
error and the spatial randomness of BSs have been characterized by
capitalizing on Kronecker correlation model and homogeneous PPP, respectively.
The outage probabilities of MIMO-NOMA systems have been derived in compact
form by considering two distinct NOMA grouping policies, including the random
grouping and the distance-based grouping. It has been proved that the average
outage probabilities are independent of the intensity of BSs in the
interference-limited regime. Interestingly, it has been found that the outage
performance is degraded under a sufficiently low intensity of BSs. Besides,
the asymptotic analyses have been carried out to gain more insights into the
outage behavior under high quality of CSI. Specifically, it has been shown
that the target rates must be limited up to a bound to achieve an arbitrarily
low outage probability in the absence of inter-cell interference. The outage
expressions have empowered the maximization of the goodput by properly
designing the precoding matrix, receiver filters and transmission rates. At
last, the numerical results have also revealed that the estimation error
correlation contributes to improve the outage performance under low quality of
CSI. Whereas, opposite observations have been made under high quality of CSI.
Furthermore, for intelligent resource allocation of NOMA enhanced SCNs, we
plan to extend the machine-learning based approaches as in [47] and [48] to
our future work.

## References

  * [1] F. Al-Turjman, E. Ever, and H. Zahmatkesh, “Small cells in the forthcoming 5G/IoT: Traffic modelling and deployment overview,” _IEEE Commun. Surveys Tuts._ , vol. 21, no. 1, pp. 28–65, First Quarter 2019.
  * [2] X. Sun and L. Dai, “Towards fair and efficient spectrum sharing between LTE and WiFi in unlicensed bands: Fairness-constrained throughput maximization,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 4, pp. 2713–2727, Apr. 2020.
  * [3] “Cisco visual networking index: Forecast and trends, 2018–2023,” Cisco, San Jose, CA, USA, Tech. Rep., Mar. 2020. [Online]. Available: https://www.cisco.com/c/en/us/solutions/collateral/executive-perspectives/annual-internet-report/white-paper-c11-741490.pdf
  * [4] M. De Ree, G. Mantas, A. Radwan, S. Mumtaz, J. Rodriguez, and I. E. Otung, “Key management for beyond 5G mobile small cells: A survey,” _IEEE Access_ , vol. 7, pp. 59 200–59 236, 2019.
  * [5] T. Han, J. Gong, X. Liu, S. R. Islam, Q. Li, Z. Bai, and K. S. Kwak, “On downlink NOMA in heterogeneous networks with non-uniform small cell deployment,” _IEEE Access_ , vol. 6, pp. 31 099–31 109, 2018.
  * [6] W. Tan, M. Matthaiou, S. Jin, and X. Li, “Spectral efficiency of DFT-based processing hybrid architectures in massive MIMO,” _IEEE Wireless Commun. Lett._ , vol. 6, no. 5, pp. 586–589, Sep. 2017.
  * [7] H. Lei, D. Wang, K.-H. Park, I. S. Ansari, J. Jiang, G. Pan, and M.-S. Alouini, “Safeguarding UAV IoT communication systems against randomly located eavesdroppers,” _IEEE Internet Things J._ , vol. 7, no. 2, pp. 1230–1244, Feb. 2019.
  * [8] Z. Wei, L. Yang, D. W. K. Ng, J. Yuan, and L. Hanzo, “On the performance gain of NOMA over OMA in uplink communication systems,” _IEEE Trans. Commun._ , vol. 68, no. 1, pp. 536–568, Jan. 2020.
  * [9] A. Nasser, O. Muta, M. Elsabrouty, and H. Gacanin, “Interference mitigation and power allocation scheme for downlink MIMO–NOMA hetnet,” _IEEE Trans. Veh. Technol._ , vol. 68, no. 7, pp. 6805–6816, Jul. 2019.
  * [10] Y. Liu, M. Zhao, L. Xiao, and S. Zhou, “Pilot domain NOMA for grant-free massive random access in massive MIMO marine communication system,” _China Commun._ , vol. 17, no. 6, pp. 131–144, Jun. 2020.
  * [11] W. Chen, H. Ding, S. Wang, D. B. da Costa, F. Gong, and P. H. J. Nardelli, “Ambient backscatter communications over NOMA downlink channels,” _China Commun._ , vol. 17, no. 6, pp. 80–100, Jun. 2020.
  * [12] Z. Ding, F. Adachi, and H. V. Poor, “The application of MIMO to non-orthogonal multiple access,” _IEEE Trans. Wireless Commun._ , vol. 15, no. 1, pp. 537–552, Jan. 2016.
  * [13] Z. Ding, R. Schober, and H. V. Poor, “A general MIMO framework for NOMA downlink and uplink transmission based on signal alignment,” _IEEE Trans. Wireless Commun._ , vol. 15, no. 6, pp. 4438–4454, Jun. 2016.
  * [14] H. Wang, Y. Fu, R. Song, Z. Shi, and X. Sun, “Power minimization precoding in uplink multi-antenna NOMA systems with jamming,” _IEEE Trans. Green Commun. Netw._ , vol. 3, no. 3, pp. 591–602, Sep. 2019.
  * [15] Y. Fu, M. Zhang, L. Salaun, C. W. Sung, and C. S. Chen, “Zero-forcing oriented power minimization for multi-cell MISO-NOMA systems: A joint user grouping, beamforming and power control perspective,” _IEEE J. Sel. Areas Commun._ , vol. PP, no. 99, pp. 1–1, 2020.
  * [16] Y. Liu, Z. Qin, M. Elkashlan, A. Nallanathan, and J. A. McCann, “Non-orthogonal multiple access in large-scale heterogeneous networks,” _IEEE J. Sel. Areas Commun._ , vol. 35, no. 12, pp. 2667–2680, Dec. 2017\.
  * [17] H. Zhang, M. Feng, K. Long, G. K. Karagiannidis, V. C. Leung, and H. V. Poor, “Energy efficient resource management in SWIPT enabled heterogeneous networks with NOMA,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 2, pp. 835–845, Feb. 2019.
  * [18] K. He, Y. Li, Y. Zhang, K. Zhang, and C. Yin, “A novel tensor CS-based NOMA MIMO system for the downlink of massive mission-critical MTC in 5G and beyond,” _IEEE Access_ , vol. 7, pp. 97 624–97 640, 2019.
  * [19] S. Silva, G. A. A. Baduge, M. Ardakani, and C. Tellambura, “NOMA-aided multi-way massive MIMO relaying,” _IEEE Trans. Commun._ , vol. PP, no. 99, pp. 1–1, 2020.
  * [20] C. Chen, W. Cai, X. Cheng, L. Yang, and Y. Jin, “Low complexity beamforming and user selection schemes for 5G MIMO-NOMA systems,” _IEEE J. Sel. Areas Commun._ , vol. 35, no. 12, pp. 2708–2722, Dec. 2017.
  * [21] J. Cui, Z. Ding, and P. Fan, “Outage probability constrained MIMO-NOMA designs under imperfect CSI,” _IEEE Trans. Wireless Commun._ , vol. 17, no. 12, pp. 8239–8255, Dec. 2018.
  * [22] G. Chen, L. Qiu, and C. Ren, “On the performance of cluster-based MIMO-NOMA in multi-cell dense networks,” _IEEE Trans. Commun._ , vol. 68, no. 8, pp. 4773–4787, Aug. 2020.
  * [23] M. Salehi, H. Tabassum, and E. Hossain, “Accuracy of distance-based ranking of users in the analysis of NOMA systems,” _IEEE Trans. Commun._ , vol. 67, no. 7, pp. 5069–5083, Jul. 2019.
  * [24] J. Park, Y. Sung, D. Kim, and H. V. Poor, “Outage probability and outage-based robust beamforming for MIMO interference channels with imperfect channel state information,” _IEEE Trans. Wireless Commun._ , vol. 11, no. 10, pp. 3561–3573, Oct. 2012.
  * [25] M. Haenggi, “User point processes in cellular networks,” _IEEE Wireless Commun. Lett._ , vol. 6, no. 2, pp. 258–261, Apr. 2017.
  * [26] S. M. Yu and S.-L. Kim, “Downlink capacity and base station density in cellular networks,” in _Proc. 11th international symposium and workshops on modeling and optimization in mobile, ad hoc and wireless networks (WiOpt’13)_ , Tsukuba Science City, Japan, May 2013, pp. 119–124.
  * [27] B. Yu, S. Mukherjee, H. Ishii, and L. Yang, “Dynamic TDD support in the LTE-B enhanced local area architecture,” in _2012 IEEE Globecom Workshops_ , Anaheim, CA, USA, Dec. 2012, pp. 585–591.
  * [28] Y. Wang, M. Haenggi, and Z. Tan, “The meta distribution of the SIR for cellular networks with power control,” _IEEE Trans. Commun._ , vol. 66, no. 4, pp. 1745–1757, Apr. 2017.
  * [29] T. Y. Al-Naffouri, M. Moinuddin, N. Ajeeb, B. Hassibi, and A. L. Moustakas, “On the distribution of indefinite quadratic forms in Gaussian random variables,” _IEEE Trans. Commun._ , vol. 64, no. 1, pp. 153–165, Jan. 2016\.
  * [30] K. Petersen and M. Pedersen, _The Matrix Cookbook_. Technical University of Denmark, 2008. [Online]. Available: https://www.math.uwaterloo.ca/~hwolkowi/matrixcookbook.pdf
  * [31] M. Haenggi, _Stochastic geometry for wireless networks_. Cambridge, U.K.: Cambridge University Press, 2012.
  * [32] J. Abate and W. Whitt, “Numerical inversion of Laplace transforms of probability distributions,” _ORSA Journal on computing_ , vol. 7, no. 1, pp. 36–43, 1995.
  * [33] A. Mathai, R. K. Saxena, and H. J. Haubold, _The H-function_. New York, NY, USA: Springer, 2009.
  * [34] L. Debnath and D. Bhatta, _Integral transforms and their applications_. Boca Raton, FL, USA: CRC press, 2010.
  * [35] M. Moorthy, “Numerical inversion of two-dimensional Laplace transforms—fourier series representation,” _Applied Numerical Mathematics_ , vol. 17, no. 2, pp. 119–127, 1995.
  * [36] J. R. Macdonald, “Accelerated convergence, divergence, iteration, extrapolation, and curve fitting,” _Journal of Applied Physics_ , vol. 35, no. 10, pp. 3034–3041, 1964.
  * [37] L. Brancik, “Error analysis and optimal parameter evaluation in FFT-based 2D-NILT method,” in _Proc. IEEE Asia-Pacific Conference on Circuits and Systems_ , vol. 2, Tainan, Taiwan, Dec. 2004, pp. 809–812.
  * [38] Z. Shi, H. Wang, Y. Fu, G. Yang, S. Ma, F. Hou, and T. A. Tsiftsis, “Zero-forcing based downlink virtual MIMO-NOMA communications in IoT networks,” _IEEE Internet Things J._ , vol. 7, no. 4, pp. 2716–2737, Apr. 2020.
  * [39] Z. Shi, S. Ma, G. Yang, and M.-S. Alouini, “Energy-efficient optimization for HARQ schemes over time-correlated fading channels,” _IEEE Trans. Veh. Technol._ , vol. 67, no. 6, pp. 4939–4953, Apr. 2018.
  * [40] W. Rui and V. K. Lau, “Combined cross-layer design and HARQ for multiuser systems with outdated channel state information at transmitter (CSIT) in slow fading channels,” _IEEE Trans. Wireless Commun._ , vol. 7, no. 7, pp. 2771–2777, Jul. 2008.
  * [41] Z. Shi, S. Ma, H. ElSawy, G. Yang, and M.-S. Alouini, “Cooperative HARQ-assisted NOMA scheme in large-scale D2D networks,” _IEEE Trans. Commun._ , vol. 66, no. 9, pp. 4286–4302, Sep. 2018.
  * [42] D. Tse, _Fundamentals of wireless communication_. Cambridge, U.K.: Cambridge university press, 2005.
  * [43] S. Boyd, S. P. Boyd, and L. Vandenberghe, _Convex optimization_. Cambridge, U.K.: Cambridge university press, 2004\.
  * [44] M. S. Ali, H. Tabassum, and E. Hossain, “Dynamic user clustering and power allocation for uplink and downlink non-orthogonal multiple access (NOMA) systems,” _IEEE Access_ , vol. 4, pp. 6325–6343, Oct. 2016.
  * [45] L. Zhu, J. Zhang, Z. Xiao, X. Cao, and D. O. Wu, “Optimal user pairing for downlink non-orthogonal multiple access (NOMA),” _IEEE Wireless Commun. Lett._ , vol. 8, no. 2, pp. 328–331, Apr. 2018.
  * [46] S. R. Islam, M. Zeng, O. A. Dobre, and K.-S. Kwak, “Resource allocation for downlink NOMA systems: Key techniques and open issues,” _IEEE Wireless Commun._ , vol. 25, no. 2, pp. 40–47, Apr. 2018.
  * [47] Z. Zhou, H. Liao, X. Zhao, B. Ai, and M. Guizani, “Reliable task offloading for vehicular fog computing under information asymmetry and information uncertainty,” _IEEE Trans. Veh. Technol._ , vol. 68, no. 9, pp. 8322–8335, Sep. 2019.
  * [48] H. Liao, Z. Zhou, X. Zhao, L. Zhang, S. Mumtaz, A. Jolfaei, S. H. Ahmed, and A. K. Bashir, “Learning-based context-aware resource allocation for edge-computing-empowered industrial IoT,” _IEEE Internet Things J._ , vol. 7, no. 5, pp. 4260–4277, May 2019.