# Zero-Forcing Based Downlink Virtual MIMO-NOMA Communications in IoT Networks

Zheng Shi, Hong Wang, Yaru Fu, Guanghua Yang, Shaodan Ma, Fen Hou, and
Theodoros A. Tsiftsis Copyright (c) 20xx IEEE. Personal use of this material
is permitted. However, permission to use this material for any other purposes
must be obtained from the IEEE by sending a request to pubs-
<EMAIL_ADDRESS>received July 6, 2019; revised October 1,
2019; revised November 11, 2019; accepted November 27, 2019. This work was
supported in part by the National Key Research and Development Program of
China under Grant 2017YFE0120600, in part by National Natural Science
Foundation of China under Grants 61801192, 61801246 and 61601524, in part by
Guangdong Basic and Applied Basic Research Foundation under Grant
2019A1515012136, in part by the Science and Technology Planning Project of
Guangdong Province under Grants 2018B010114002 and 2019B010137006, in part by
the Science and Technology Development Fund, Macau SAR (File no. 0036/2019/A1,
File no. 037/2017/AMJ and File no. SKL-IOTSC2018-2020), in part by Open
Research Foundation of National Mobile Communications Research Laboratory of
Southeast University under Grant 2018D09, and in part by the Research
Committee of University of Macau under Grants MYRG2018-00156-FST and
MYRG2016-00171-FST. The associate editor coordinating the review of this paper
and approving it for publication was Jun Wu. (Corresponding author: Guanghua
Yang.)Zheng Shi, Guanghua Yang and Theodoros A. Tsiftsis are with the
Institute of Physical Internet and the School of Intelligent Systems Science
and Engineering, Jinan University, Zhuhai 519070, China (e-mails:
<EMAIL_ADDRESS><EMAIL_ADDRESS>theo_tsiftsis@jnu.edu.cn).H. 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 Department of Information Systems Technology and Design, Singapore
University of Technology and Design, Singapore 487372
(email:yaru_fu@sutd.edu.sg).Shaodan Ma and Fen Hou are with the State Key
Laboratory of Internet of Things for Smart City and the Department of
Electrical and Computer Engineering, University of Macau, Macao S.A.R. 999078,
China (e-mails<EMAIL_ADDRESS>fenhou@um.edu.mo).

###### Abstract

To support massive connectivity and boost spectral efficiency for internet of
things (IoT), a downlink scheme combining virtual multiple-input multiple-
output (MIMO) and non-orthogonal multiple access (NOMA) is proposed. All the
single-antenna IoT devices in each cluster cooperate with each other to
establish a virtual MIMO entity, and multiple independent data streams are
requested by each cluster. NOMA is employed to superimpose all the requested
data streams, and each cluster leverages zero-forcing detection to de-
multiplex the input data streams. Only statistical channel state information
(CSI) is available at base station to avoid the waste of the energy and
bandwidth on frequent CSI estimations. The outage probability and goodput of
the virtual MIMO-NOMA system are thoroughly investigated by considering
Kronecker model, which embraces both the transmit and receive correlations.
Furthermore, the asymptotic results facilitate not only the exploration of
physical insights but also the goodput maximization. In particular, the
asymptotic outage expressions provide quantitative impacts of various system
parameters and enable the investigation of diversity-multiplexing tradeoff
(DMT). Moreover, power allocation coefficients and/or transmission rates can
be properly chosen to achieve the maximal goodput. By favor of Karush-Kuhn-
Tucker conditions, the goodput maximization problems can be solved in closed-
form, with which the joint power and rate selection is realized by using
alternately iterating optimization. Besides, the optimization algorithms tend
to allocate more power to clusters under unfavorable channel conditions and
support clusters with higher transmission rate under benign channel
conditions.

###### Index Terms:

Internet of things (IoT), Kronecker model, non-orthogonal multiple access
(NOMA), virtual multiple-input multiple-output (MIMO), zero-forcing (ZF).

##  I Introduction

With the phenomenal growth of the internet of things (IoT) devices in 5G
scenarios (including smart city, connected health, industrial internet,
vehicle network), a huge number of wireless-enabled sensors are broadly
deployed. A great deal of sensing data is collected sporadically, and then
forwarded to the cloud via connected IoT devices, cellular networks and
internet [1, 2]. It has been predicted in Ericsson’s report that the number of
short-range and cellular IoT devices will reach 17.8 billion and 4.1 billion
by 2024, respectively [3]. Moreover, the visual networking index (VNI)
released by Cisco has forecasted that machine to machine (M2M) connections in
IoT networks will grow from 6.1 billion in 2017 to 14.6 billion by 2022 [4].
The proliferation of IoT devices will pose the necessity of massive wireless
connections to confront the ongoing paradigm of 5G and beyond. In addition,
the non-orthogonal multiple access (NOMA) has been recognized as a promising
technology enabler for IoT to offer massive connectivity over limited wireless
resources [5, 6, 7, 8]. The essence of NOMA is to serve multiple users over
the same frequency/time/code resource by means of superposition coding and
successive interference cancellation (SIC). In particular, NOMA is capable of
attaining an improved spectral efficiency over its conventional counterpart,
orthogonal multiple access (OMA). In addition, NOMA can also strike a balanced
tradeoff between the system throughput and users’ fairness. All these
intrinsic characteristics of NOMA are beneficial to the development and
deployment of IoT. Moreover, in order to get rid of the high complexity and
severe error propagation of SIC, the combination of NOMA and OMA technologies
has been recognized as a feasible solution for resource-limited IoT networks
[9, 7]. Particularly, the hybrid multiple access (MA) scheme first groups
users into pairs to perform NOMA, then multiplexes different pairs over
orthogonal resources.

Despite the prominent advantages of NOMA, other enhanced technologies are also
needed to further boost the throughput and reliability of IoT networks. Note
that multiple-input-multiple-output (MIMO) can exploit extra degrees of
freedom [10, 11], and thus the integration of NOMA with MIMO has became an
appealing solution for both massive connectivity support and capacity
provision recently [12].

###  I-A Related Works

The maximization of the ergodic capacity was studied for MIMO-NOMA system with
a pair of users in [13]. To accommodate as many pairs of users as possible,
the concept of signal alignment was devised to cancel out inter-pair
interference within MIMO-NOMA framework in [14, 15]. By using the signal
alignment, the signal receptions between different pairs of users can be
handled independently. Hence, the investigations of MIMO-NOMA in [13, 15, 14]
were essentially specific to single/multiple pairs of NOMA users. To support
more than two users in each NOMA group, multi-user MIMO-NOMA transmission
strategies were developed in [16, 17] by leveraging group interference
cancellation in conjunction with minimum mean square error (MMSE) detector. To
meet dynamic QoS requirements in IoT communications, novel precoding and
detection strategies were proposed for MIMO-NOMA systems to differ users’
effective channel gains in [18]. A coordinated multipoint (CoMP)-enabled NOMA
scheme was proposed to attain an improved spectral efficiency in [19], where
distributed single-antenna base stations (BSs) collaborate together to form a
virtual MIMO. Nonetheless, most of the existing works assumed frequent
feedback of channel state information (CSI) at the transmitter. This
assumption leads to high energy consumption, signaling and computational
overhead, which is apparently impractical for bandwidth-limited, energy-
constrained and latency-sensitive IoT networks [20, 21]. In [22, 23], more
effective CSI feedback was considered into NOMA-MIMO schemes. For further
suppression of system overhead, the sender was assumed to have only
statistical knowledge of CSI [24, 25, 26, 27, 28]. To be specific, the rate
and outage performance was analyzed for single-stream MIMO-NOMA systems in
[24]. In regard to the similar systems, a power allocation scheme was further
proposed for maximizing the ergodic sum rate in [25]. Furthermore, a linear
MIMO-NOMA model with complex-valued power coefficients was derived by
considering multi-stream transmissions in [26]. Unlike [24, 25, 26] that
assume independent fading channels, the joint design of power allocation and
receive beamforming was developed for single-stream MIMO-NOMA systems with
receive correlation in [27]. Whereas, the outage performance of multi-stream
MIMO-NOMA systems was examined in [28] by considering transmit correlation.

IoT devices usually equip with a single antenna for the sake of limited size,
low cost and energy saving, which in turn hampers the exploitation of spatial
multiplexing gain from MIMO technology. To combat this issue, virtual MIMO
technology has emerged as an attractive candidate by coordinating multiple
nodes to form one virtual antenna array entity [29, 30]. The realization of
virtual MIMO brings new challenges of data sharing, time/frequency
synchronization, user grouping/clustering, and CSI estimation [31, 32].
Nevertheless, it was disclosed in [33, 34] that the virtual MIMO system
surpasses the simple single-input single-output (SISO) system in terms of the
consumed energy and delay. In [35], the virtual MIMO was proved to be a viable
framework to improve the data rate and energy efficiency. Furthermore,
experimental measurements were carried out for the virtual MIMO to confirm its
throughput enhancement in indoor-to-outdoor scenarios in [36]. In order to
prolong the lifetime of IoT networks as well as improve its capacity, the
application of virtual MIMO to IoT has received an ever-increasing research
attention recently[37, 38, 39, 40]. To name a few, in [37], a multi-antenna
unmanned aerial vehicle (UAV) serves as an on-demand flexible platform, which
communicates with a cluster of single-antenna IoT devices through a virtual
MIMO link. In [38], a virtual MIMO approach was invented for the delivery of
sensor data to a distant sink node in wireless body area networks (WBAN).

###  I-B Scope of the Paper

To offer massive number of connections and implement MIMO transmissions in IoT
scenarios, this paper proposes a novel virtual MIMO-NOMA transmission scheme
for IoT communications. More specifically, downlink communications between a
BS and multiple clusters are taken into account, wherein all the single-
antenna IoT devices in each cluster collaborates with each other to create a
virtual MIMO entity, and a number of independent data streams are requested by
each cluster. On one hand, NOMA technique is employed at the BS to
simultaneously serve all the clusters in power-domain. On the other hand, each
cluster capitalizes on zero-forcing (ZF) detection to de-multiplex all the
desired input data streams. To avoid the waste of the energy and bandwidth on
frequent channel estimations, the BS is assumed to have only statistical
knowledge of CSI. An in-depth analysis of the outage and goodput performance
of the virtual MIMO-NOMA system is carried out, where a more practical channel
model, i.e., Kronecker model, is adopted to account for both the transmit and
receive correlations. More importantly, the asymptotic outage probability is
derived in a compact and simple form, which offers not only insightful
conclusions but also a simple analytical treatment for the optimal system
design. With the asymptotic results, the maximization of the goodput is
enabled with closed-form solution. Moreover, a modified version of the
proposed scheme, i.e., virtual MIMO-hybrid MA, is devised to balance the
tradeoff between two conflicting goals of enhancing spectral efficiency and
saving computation/energy resources for IoT networks. Numerical verifications
are finally presented.

###  I-C Contributions

The main contributions of this paper can be summarized as follows:

  1. 1.

The outage probability and goodput of the downlink virtual MIMO-NOMA-aided IoT
communication system are scrutinized in this paper. To characterize the impact
of spatial correlation between antennas in realistic propagation environments,
Kronecker model is employed to capture both the transmit and receive
correlations. However, Kronecker channel model poses significant challenges to
the performance evaluation of the proposed scheme. As the most fundamental
performance metric, the outage probability can be accurately approximated in a
compact form by invoking the methodologies of moment generating function (MGF)
and multivariate analysis.

  2. 2.

The asymptotic outage probability is then derived in closed-form to ease the
extraction of insightful results. In particular, the asymptotic outage
probability quantifies the impacts of spatial correlation, transmission rates
and power allocation coefficients. In addition, the diversity order is
$\mathpzc d=N_{r}-M+1$, where $M$ and $N_{r}$ stand for the number of data
streams requested by each cluster and the number of IoT devices within each
cluster, respectively.

  3. 3.

Furthermore, diversity-multiplexing tradeoff (DMT) is studied to figure out
the multiplexing benefit of MIMO-NOMA channels.

  4. 4.

The asymptotic outage probability also enables the maximization of goodput in
closed-form by favor of Karush-Kuhn-Tucker (KKT) conditions. Three
optimization algorithms are thereon developed for the proper selection of the
transmission powers and/or rates. It is revealed that the developed algorithms
incline to allocate more power to clusters under unfavorable channel
conditions and support clusters with higher transmission rate under benign
channel conditions.

###  I-D Outline and Notations

The remainder of the paper is outlined as follows. The virtual MIMO-NOMA
system model is presented together with two fundamental performance metrics in
Section II. Section III conducts the exact and asymptotic analyses of the
outage probability. With the analytical results, the goodput maximization is
then enabled in Section IV. The numerical results are shown for verifications
in Section V.

_Notations:_ We shall use the following notations throughout the paper. Bold
uppercase and lowercase letters are used to denote matrices and vectors,
respectively. ${\bf X}^{\mathrm{T}}$, ${\bf X}^{\mathrm{H}}$, ${\bf X}^{-1}$
and ${\bf X}^{1/2}$ denote the transpose, conjugate transpose, matrix inverse
and Hermitian square root of matrix ${\bf X}$, respectively. ${\bf
X}^{\dagger}$ represents the Moore-Penrose pseudo-inverse of $\bf X$ and ${\bf
X}^{\dagger}={\left({{{\bf{X}}}^{\rm{H}}{{\bf{X}}}}\right)^{-1}}{{\bf{X}}}^{\rm{H}}$.
$\mathrm{vec}$, $\rm{tr}$, $\mathrm{det}$ and $\mathrm{diag}$ are the
operators of vectorization, trace, determinant and diagonalization,
respectively. $[{\bf X}]_{ij}$ is the $(i,j)$-th element of $\bf X$.
$\mathbf{0}_{n}$ and $\mathbf{I}_{n}$ stand for $1\times n$ all-zero vector
and $n\times n$ identity matrix, respectively. $\mathbb{C}^{m\times n}$
denotes the set of $m\times n$-dimensional complex matrices. The symbol ${\rm
i}=\sqrt{-1}$ is the imaginary unit. ${{\mathrm{E}}_{A}}$ denotes the operator
of the expectation taking over random variable $A$. $o(\cdot)$ denotes
little-O notation. $(\cdot)_{n}$ represents Pochhammer symbol.
$\bigcup{(\cdot)}$ stands for the union operation. The symbol $\simeq$ means
“asymptotically equal to”. $\propto$ stands for “directly proportional to”.
$\Gamma(x)$ denotes Gamma function. $G^{m,n}_{p,q}(\left.\cdot\right|x)$
stands for Meijer G-function. $\gamma(s,x)$ denotes the lower incomplete gamma
function. Any other notations will be defined in the place where they take
place.

##  II System Model and Performance Metrics

Figure 1: A downlink virtual MIMO-NOMA system model.

###  II-A System Model

As shown in Fig. 1, this paper considers a downlink virtual MIMO-NOMA system
with one BS/central sink node serving $K$ clusters in IoT networks. Herein, we
assume that the BS is equipped with $N_{t}$ antennas. Besides, the $K$
clusters are predefined to ease management, and each cluster consists of
$N_{r}$ neighbouring single-antenna IoT devices111In general, the
pairing/clustering process could be either predefined or dynamically
scheduled. In particular, similarly to [37], the predefined clustering divides
all the IoT devices into multiple clusters, each with $N_{r}$ adjacent
devices. It is noteworthy that our subsequent analysis is also applicable to
the case of clusters with unequal sizes. Whereas, with regard to the dynamic
clustering, the BS first sends a message to notify all the IoT devices within
the cell so as to trigger the clustering process via a self-organizing
protocol. Upon receiving the notification message, each device periodically
broadcasts a life messages (LifeMsg) containing its unique identity, and
meanwhile searches its neighbouring devices within a certain range by
overhearing their LifeMsgs. After the completion of the search, all the
neighbouring devices collaborate and share their antennas with each other to
form a cluster. Each cluster then reports its relative position (or path loss)
and statistical knowledge of CSI (e.g., channel correlations, average channel
gains) to the BS. Once the BS gathers all the clusters’ reports, it determines
the decoding order of SIC for these clusters according to their relative
positions, and then the transmission rates and powers are properly devised
based on the channel statistics [41, 9, 42]. The necessary information (e.g.,
decoding order) for detections is thereafter fed back from the BS to the
clusters. Nonetheless, the further investigation of implementing the
pairing/clustering is out of the scope of this paper. For more details,
interested readers are referred to [23, 7, 43].[37]. These devices share their
antennas with each other through an error-free intra-cluster interface, in
order to create a virtual MIMO system[44, 45]. We can easily extend our design
to the general case with more than one antennas mounted on the IoT devices.
This is based on the fact that an IoT device with multiple antennas can be
viewed as multiple single-antenna IoT devices. Moreover, we assume that
$M(\leq\min{(N_{t},N_{r})})$ independent data streams are requested by each
cluster. In doing so, a power-domain NOMA is used to multiplex these clusters
over the same time-frequency resource with the help of superposition coding at
the BS. Thus, the received signal at the $k$-th nearest cluster to the BS is
given by

${{\bf{y}}_{k}}={\sqrt{\ell(d_{k})}}{{\bf{H}}_{k}}{{\bf{V}}}{{\bf{s}}}+{{\bf{n}}_{k}},$
(1)

where ${{\bf{H}}_{k}}\in{\mathbb{C}}^{N_{r}\times N_{t}}$ denotes the channel
response matrix from the BS to cluster $k$; ${\bf V}=({\bf v}_{1},\cdots,{\bf
v}_{M})\in{\mathbb{C}}^{N_{t}\times M}$ represents the transmit beamforming
matrix at the BS, i.e., $\left\|{{{\bf{v}}_{m}}}\right\|=1$ for
$m=1,\cdots,M$; ${{\bf{s}}}=({s}_{1},\cdots,{s}_{M})^{\rm T}$ is the
superimposed symbol vector at the BS, and the entries of the transmit symbol
vector are drawn independently from a complex symmetric circularly normal
distributions with zero-mean and variance of $P$, i.e., ${{\bf{s}}}\sim{\cal
CN}({\bf 0},P{\bf I}_{M})$, and $P$ denotes the total transmit power of each
data stream; ${{\bf{n_{k}}}}$ is the additive Gaussian white noise (AWGN)
vector with variance $\sigma^{2}$, i.e., ${{\bf{n_{k}}}}\sim{\cal CN}({\bf
0},\sigma^{2}{\bf I}_{N_{r}})$. A composite channel model is applied to
encompass both small-scale Rayleigh fading and large-scale path loss. In
particular, the path loss is modelled by using Friis equation, such that
$\ell(d_{k})=\mathcal{K}{d_{k}}^{-\alpha}$ with path loss exponent $\alpha\geq
2$, where $d_{k}$ denotes the distance between the BS and cluster $k$,
$\mathcal{K}$ represents the free space power received at the reference
distance $d_{0}=1$m. It is noticeable that all the IoT devices within each
cluster are assumed to suffer from almost identical large-scale fadings
because of a sufficiently long distance between the BS and clusters relative
to the cluster size [46]. Without loss of generality, the distances between
the BS and clusters are sorted in an increasing order such that $d_{1}\leq
d_{2}\leq\cdots\leq d_{K}$.

Due to the close spacing between adjacent antenna elements no matter at the BS
or within each cluster, the resultant mutual antenna coupling brings about the
spatial correlation between elements in realistic propagation
environments[47]. Particularly in downlink IoT communications, high
correlation between the transmit antennas may occur in the absence of local
scattering [48]. Furthermore, in IoT network, devices or sensors located
indoors (e.g., airports, malls, house and buildings) usually suffer from
shadowing effects. The shadowing induced by the same obstacle (tree/building)
will incur high correlation between propagation paths[46]. Overlooking the
existence of spatial correlation will overestimate the system performance and
offer misleading design guide. In order to capture the negative impact of the
antenna correlation at transceivers, Kronecker channel model is widely adopted
by decomposing the spatial correlation into two independent components, i.e.,
transmit and receive correlations[49]. According to the Kronecker model, the
channel matrix between the BS and cluster $k$, ${{\bf{H}}_{k}}$, is
represented by [50, 51]

${{\bf{H}}_{k}}={{\bf R}_{r}}^{1/2}{\bf H}_{w}{{\bf R}_{t}}^{1/2},$ (2)

where ${\bf R}_{t}$ and ${\bf R}_{r}$ denote transmit- and receive-side
covariance matrices, respectively, the elements of matrix ${\bf H}_{w}$ are
independent and identically distributed (i.i.d.) complex zero-mean unit-
variance Gaussian random variables. Thus, the distribution of ${\bf H}_{k}$
follows as ${\rm vec}\left({\bf H}_{k}\right)\sim{\cal CN}({\bf
0}_{N_{t}N_{r}},{{\bf R}_{t}}^{\rm T}\otimes{\bf R}_{{{r}}})$. It is
worthwhile to point out that ${\bf R}_{t}$ and ${\bf R}_{r}$ are positive
semi-definite Hermitian matrices. The validity of Kronecker model has been
corroborated regardless of antenna configurations and intra-array spacing,
provided that the transmit and receive angular power profiles are separable
[52]. Besides, in order to conserve energy, relieve signaling overhead and
avoid frequent CSI feedback in large-scale energy- and bandwidth-limited IoT
networks, we assume that only statistical CSI is available to the BS [20, 21].

Since NOMA is employed to simultaneously accommodate multiple clusters by
sharing the same physical resources, the information bearing vector of $M$
data streams can be constructed by using superposition coding as follows

${\bf{s}}=\sqrt{P}\sum\limits_{k=1}^{K}{{\rm
diag}\left\\{\sqrt{\zeta_{1,k}},\cdots,\sqrt{\zeta_{M,k}}\right\\}{{\bf{x}}_{k}}},$
(3)

where ${{\bf{x}}_{k}}=(x_{1,k},\cdots,x_{M,k})^{\rm T}$ is the signal vector
intended for cluster $k$, $\zeta_{m,k}$ stands for the power allocation
coefficient of signal $x_{m,k}$, which is subject to the normalization
constraint $\sum\nolimits_{k=1}^{K}{{\zeta_{m,k}}}={1}/{M}$ for all
$m\in[1,M]$, to ensure the fairness among different data streams. On the other
hand, the application of NOMA can guarantee clusters’ fairness by taking
advantage of the difference in channel qualities between the BS and
clusters[23].

The linear ZF detection is a vital technology to offer low complexity and
implementation cost for IoT communications[53, 18]. Hence, we assume that the
received signal is reconstructed by using a ZF detector. Specifically, the ZF
estimate of the transmitted symbol streams at cluster $k$ can be obtained by
multiplying ${\bf y}_{k}$ with
$({\sqrt{\ell(d_{k})}}{{{{\bf{H}}_{k}}{\bf{V}}}})^{\dagger}$ as

$\displaystyle{\bf{\hat{s}}}$
$\displaystyle={\bf{s}}+\frac{1}{\sqrt{\ell(d_{k})}}\left({{{{\bf{H}}_{k}}{\bf{V}}}}\right)^{\dagger}{{\bf{n}}_{k}}.$
(4)

With the linear ZF receiver, the decodings for $M$ data streams can be
independently performed. This eventually induces a very low decoding
complexity.

After ZF detection, the SIC as one of the key components of NOMA technique is
applied to separate superimposed signals sequentially. To simplify the design
of power allocation coefficient and better utilize NOMA scheme, the decoding
order at receivers is determined according to the decreasing order of
transmission distances, i.e., $d_{K}\geq\cdots\geq d_{2}\geq d_{1}$[54, 9]. In
particular, cluster $k$ will detect the message of cluster $i$ prior to
decoding its own message ${\bf x}_{k}$ such that $i>k$, and then subtracting
the signal ${\bf x}_{i}$ from its observation successively. In the meantime,
the signal ${\bf x}_{i}$ intended for cluster $i$ ($<k$) will be regarded as
interference at cluster $k$. Therefore, for the $k$-th nearest cluster to the
BS, the instantaneous received signal-to-interference-plus-noise ratio (SINR)
to detect the signal $x_{m,i}$ that is extracted from data stream $m$ is given
by

$\displaystyle{\gamma_{m,k\to i}}$
$\displaystyle=\frac{{P{\zeta_{m,i}}}}{{P\sum\limits_{l=1}^{i-1}{{\zeta_{m,l}}}+\frac{\sigma^{2}}{\ell({d_{k}})}{{\left[{{{\left({{{\bf{H}}_{k}}{\bf{V}}}\right)}^{\dagger}}{{\left({{{\left({{{\bf{H}}_{k}}{\bf{V}}}\right)}^{\dagger}}}\right)}^{\rm{H}}}}\right]}_{mm}}}}$
$\displaystyle=\frac{{{{\bar{\gamma}}}{\ell({d_{k}})}{\zeta_{m,i}}}}{{{{\bar{\gamma}}}{\ell({d_{k}})}\sum\limits_{l=1}^{i-1}{{\zeta_{m,l}}}+{{\left[{{{\left({{{\left({{{\bf{H}}_{k}}{\bf{V}}}\right)}^{\rm{H}}}{{\bf{H}}_{k}}{\bf{V}}}\right)}^{-1}}}\right]}_{mm}}}},$
(5)

where $i\geq k$ and ${{\bar{\gamma}}}={{P}}/{{{\sigma^{2}}}}$ denotes the
average transmit SNR.

###  II-B Performance Metrics

The absence of perfect instantaneous CSI at the BS undoubtedly results in
decoding errors at receivers. In this paper, we are concerned with the
error/outage probability. More specifically, the event that any cluster fails
to decode any of its own data streams can be deemed as an error/outage. By
following the proposed scheme, the outage event takes place if the cluster
(e.g., cluster $k$) is unable to either subtract the interfering streams
(e.g., $x_{m,i}$, $i>k$) or recover its own data stream (e.g., $x_{m,k}$).
From an information-theoretical perspective, the outage probability of cluster
$k$ for decoding its own data stream $m$ is therefore expressed as

$\displaystyle
p_{m,k}^{out}=\Pr\left[{\bigcup\limits_{i=k}^{K}{{{\log}_{2}}\left({1+{\gamma_{m,k\to
i}}}\right)<{R_{m,i}}}}\right],$ (6)

where ${R_{m,i}}$ is the target transmission rate for data stream $x_{m,k}$.

Since the occurrence of decoding failures will yield the deterioration of the
target achievable rate, the long-term average throughput should be precisely
estimated. Towards this goal, the concept of goodput is commonly used to
measure the number of information bits successfully delivered per
transmission[55]. The goodput of the virtual MIMO-NOMA system can be expressed
in terms of outage probabilities as

${T_{g}}=\sum\limits_{m=1}^{M}{\sum\limits_{k=1}^{K}{\left({1-p_{m,k}^{out}}\right){R_{m,k}}}}.$
(7)

Unambiguously, the outage probability plays a crucial role in performance
evaluation. For the sake of thorough investigation, we have to derive a
closed-form and tractable expression for the outage probability first. With
the analytical result, we are concerned with the maximization of the goodput
for the virtual MIMO-NOMA system in this paper.

##  III Analysis of Outage Probability

Although the outage probability of the MIMO-NOMA system under semi-correlated
MIMO fading channels has been analyzed in [28], the analytical results are
inapplicable to that under full-correlated MIMO fading channels, as shown in
(2). The Kronecker correlation model in (2) significantly challenges the exact
analysis of the outage probability. Moreover, the asymptotic analysis of the
outage probability is also conducted to gain more helpful insights.

###  III-A Exact Analysis

By putting (II-A) into (6), the probability of the event that cluster $k$
fails to decode its own data stream $m$, $p_{m,k}^{out}$, can be simplified as

$\displaystyle p_{m,k}^{out}$ $\displaystyle=$
$\displaystyle\Pr\left[{{{{{\left[{{{\left({{{\left({{{\bf{H}}_{k}}{\bf{V}}}\right)}^{\rm{H}}}{{\bf{H}}_{k}}{\bf{V}}}\right)}^{-1}}}\right]}_{mm}}}}>{{{{\bar{\gamma}}}{\theta_{m,k}}\ell({d_{k}})}}}\right],$
(8)

where
${\theta_{m,k}}=\min\left\\{{\left.{{{{\zeta_{m,i}}}}/{({{2^{{R_{m,i}}}}-1})}-\sum\nolimits_{l=1}^{i-1}{{\zeta_{m,l}}}}\right|k\leq
i\leq K}\right\\}$. It is worth noting that the condition of
${\theta_{m,k}}>0$ should be always satisfied to successfully eliminate inter-
cluster interference while performing SIC for NOMA. Hereby, the constraints
$\left\\{{\theta_{m,k}}>0,k\in[1,K]\right\\}$ lead to
$\left\\{{R_{m,i}}<{\log_{2}}\left({1+{{{\zeta_{m,i}}}}/{{\sum\nolimits_{l=1}^{i-1}{{\zeta_{m,l}}}}}}\right),\,i\in[2,K]\right\\}$.

By defining ${\bf{Z}}={{{\bf{H}}_{k}}{\bf{V}}}$, the elements of ${\bf{Z}}$
are circularly symmetric complex Gaussian distributed as ${\rm vec}\left({\bf
Z}\right)\sim{\cal CN}({\bf 0},{{{\bf{R}}_{t^{\prime}}}}^{\rm T}\otimes{\bf
R}_{{{r}}})$, where
${{\bf{R}}_{t^{\prime}}}={{\bf{V}}^{\rm{H}}}{{\bf{R}}_{t}}{\bf{V}}$. Then,
${\bf{Z}}$ can be rewritten as
${\bf{Z}}={{\bf{R}}_{r}}^{1/2}{{\bf{Z}}_{w}}{{\bf{R}}_{t^{\prime}}}^{1/2}$,
where ${{\bf{Z}}_{w}}$ denotes a $N_{r}\times M$ complex Gaussian matrix with
i.i.d. ${\cal CN}(0,1)$ entries. With the above definitions, $p_{m{},k}^{out}$
can be rewritten as

$p_{m{},k}^{out}=\Pr\left[{{{{{\left[{{{\left({{{\bf{Z}}^{\rm{H}}}{\bf{Z}}}\right)}^{-1}}}\right]}_{mm}}}}>{{{{\bar{\gamma}}}{\theta_{m,k}}\ell({d_{k}})}}}\right].$
(9)

Unfortunately, there are only a few relevant research results concerning (9)
so far. More specifically, by solely considering the transmit correlation,
[56, Theorem 3.2.10] can be employed to get (9) in closed-form [57, 22, 28].
However, the introduction of the receive correlation will exceedingly
complicate the performance analysis due to the intractability of tackling the
integral involving hypergeometric functions of matrix argument. By employing
the vector-matrix partition, M. Kiessling _et al._ in [58] investigated the
symbol error rate of a $2\times 2$ MIMO case with correlations occurring at
both transmit- and receive-sides. However, there is still no general result
available for deriving (9). Motivated by this, the outage probability for
arbitrary data stream $m$, $p_{m,k}^{out}$, can be accurately approximated as

$\displaystyle
p_{m,k}^{out}\approx\tilde{F}\left(\frac{{{{\left[{{{\bf{R}}_{t^{\prime}}}^{-1}}\right]}_{mm}}}}{{{{\bar{\gamma}}}{\theta_{m,k}}\ell({d_{k}})}}\right),$
(10)

where $\tilde{F}\left(x\right)$ is given by (11) at the top of the next page,
${\lambda_{1}},\cdots,{\lambda_{N_{r}}}$ denote the $N_{r}$ eigenvalues of
${{\bf{R}}_{r}}$and the proof is detailed in Appendix A.

$\tilde{F}\left(x\right)=\frac{{\Gamma\left({{N_{r}}}\right)\det\left({{{\left\\{{{\lambda_{i}}^{{N_{r}}-1}G_{2,3}^{2,1}\left({\left.{\begin{array}[]{*{20}{c}}{1,{N_{r}}}\\\
{1,{N_{r}}-M+1,0}\end{array}}\right|\frac{x}{{{\lambda_{i}}}}}\right)}\right\\}}_{\scriptstyle
1\leq i\leq{N_{r}}\hfill\atop\scriptstyle
j=1\hfill}},{{\left\\{{{\lambda_{i}}^{{N_{r}}-j}}\right\\}}_{\scriptstyle
1\leq i\leq{N_{r}}\hfill\atop\scriptstyle 2\leq
j\leq{N_{r}}\hfill}}}\right)}}{{\Gamma\left({{N_{r}}-M+1}\right)\det\left({\left\\{{{\lambda_{i}}^{{N_{r}}-j}}\right\\}_{1\leq
i,j\leq{N_{r}}}}\right)}}.$ (11)

It it noticeable that the approximation in (10) will become an equality
whenever ${{{\bf{R}}_{r}}}=\lambda_{0}{\bf I}_{N_{r}}$. Clearly, the
approximation in (10) would become accurate if fading channels experience a
low correlation at cluster’s side. In practice, this assumption makes sense
because IoT devices are placed sufficiently far apart from each other relative
to the close-spacing between transmit antennas at the base station. Moreover,
(10) can be further simplified by considering the following two special cases.

####  III-A1 $M=1$

Suppose that only a single data stream is requested by each cluster in
downlink IoT communications, i.e., $M=1$, the Meijer-G function in (11)
reduces to

$\displaystyle
G_{2,3}^{2,1}\left({\left.{\begin{array}[]{*{20}{c}}{1,{N_{r}}}\\\
{1,{N_{r}}-M+1,0}\end{array}}\right|x}\right)$
$\displaystyle=G_{1,2}^{1,1}\left({\left.{\begin{array}[]{*{20}{c}}1\\\
{1,0}\end{array}}\right|x}\right)$ (16)
$\displaystyle=\gamma\left({1,x}\right)=1-{e^{-x}},$ (17)

where the first equality holds by using the definition of Meijer G-function
[59], and the second equality holds by using [60, eq.06.07.26.0006.01].
Hereby, the outage probability can be rewritten as (18), shown at the top of
the following page.

$p_{m,k}^{out}=\frac{{\det\left({{{\left\\{{{\lambda_{i}}^{{N_{r}}-1}\left({1-{e^{-\frac{{{{\left[{{{\bf{R}}_{t^{\prime}}}^{-1}}\right]}_{mm}}}}{{{\lambda_{i}}\bar{\gamma}{\theta_{m,k}}\ell({d_{k}})}}}}}\right)}\right\\}}_{\scriptstyle
1\leq i\leq{N_{r}}\hfill\atop\scriptstyle
j=1\hfill}},{{\left\\{{{\lambda_{i}}^{{N_{r}}-j}}\right\\}}_{\scriptstyle
1\leq i\leq{N_{r}}\hfill\atop\scriptstyle 2\leq
j\leq{N_{r}}\hfill}}}\right)}}{{\det\left({\left\\{{{\lambda_{i}}^{{N_{r}}-j}}\right\\}_{1\leq
i,j\leq{N_{r}}}}\right)}},$ (18)

####  III-A2 $\lambda_{1}=\cdots=\lambda_{N_{r}}=\lambda_{0}$

If all the eigenvalues of ${{\bf{R}}_{r}}$ are identical, i.e.,
$\lambda_{1}=\cdots=\lambda_{N_{r}}=\lambda_{0}$, the approximation in (10)
turns out to be an equality. By substituting (11) into (10), we discover that
both the denominator and numerator of (11) equal to zero because of identical
rows in determinant. Hence, as proved in Appendix C, $p_{m,k}^{out}$ can be
simplified as

$\displaystyle p_{m,k}^{out}$
$\displaystyle=\frac{{\gamma\left({{N_{r}}-M+1,\frac{{{{\left[{{{\bf{R}}_{t^{\prime}}}^{-1}}\right]}_{mm}}}}{{{\lambda_{0}}{{\bar{\gamma}}}{\theta_{m,k}}\ell({d_{k}})}}}\right)}}{{\Gamma\left({{N_{r}}-M+1}\right)}},$
(19)

which is consistent with the result in [28].

###  III-B Asymptotic Analysis

To extract more meaningful results from the exact analysis and facilitate the
optimal system design as well, this paper focuses on the asymptotic behaviour
of the outage probability under high SNR, i.e., $\bar{\gamma}\to\infty$. To
this end, it is imperative to carry out the asymptotic analysis of
$\tilde{F}\left(x\right)$ as $x\to 0$. As proved in Appendix D,
$\tilde{F}\left(x\right)$ can be expanded as

$\displaystyle\tilde{F}$
$\displaystyle\left(x\right)=\frac{{\Gamma\left({{N_{r}}}\right)}}{{\Gamma\left({{N_{r}}-M+1}\right)\det\left({\left\\{{{\lambda_{i}}^{{N_{r}}-j}}\right\\}}\right)}}$
$\displaystyle\times\sum\limits_{\vartheta={N_{r}}-M+1}^{{N_{r}}-1}{\frac{{{{\left({-1}\right)}^{2\vartheta-{N_{r}}+M}}{{x}^{\vartheta}}}}{{\vartheta!\left({\vartheta-\left({{N_{r}}-M+1}\right)}\right)!\left({{N_{r}}-1-\vartheta}\right)!}}}$
$\displaystyle\times\det\left(\left\\{{{{{\lambda_{i}}}}^{{N_{r}}-\vartheta-1}}\ln{{\lambda_{i}}},\left.{{{{{\lambda_{i}}}}^{{N_{r}}-j}}}\right|_{j=2}^{{N_{r}}}\right\\}\right)$
$\displaystyle+\frac{{\Gamma\left({{N_{r}}}\right)}}{{\Gamma\left({{N_{r}}-M+1}\right)\det\left({\left\\{{{\lambda_{i}}^{{N_{r}}-j}}\right\\}}\right)}}$
$\displaystyle\times\sum\limits_{\vartheta={N_{r}}}^{\infty}{\frac{{{{\left({-1}\right)}^{\vartheta+M-2}}\left({\vartheta-{N_{r}}}\right)!}}{{\vartheta!\left({\vartheta-{N_{r}}+M-1}\right)!}}{x^{\vartheta}}}$
$\displaystyle\times\det\left(\left\\{{{{{\lambda_{i}}}}^{{N_{r}}-\vartheta-1}},\left.{{{{{\lambda_{i}}}}^{{N_{r}}-j}}}\right|_{j=2}^{{N_{r}}}\right\\}\right).$
(20)

With the series expansion of $\tilde{F}\left(x\right)$ in (III-B), the
asymptotic analysis of the outage probability is thus empowered. The first
summation of (III-B) in the right-hand side of equation exists if and only if
$M>1$, thus the asymptotic analysis of $p_{m,k}^{out}$ should be separated
into two cases, i.e., $M=1$ and $M>1$.

####  III-B1 $M=1$

If only a single data stream in each cluster is assumed, i.e., $M=1$, the
asymptotic expression of $p_{m,k}^{out}$ as $\bar{\gamma}\to\infty$ can be
obtained by putting (III-B) into (10) as

$\displaystyle p_{m,k}^{out}$
$\displaystyle\simeq\frac{1}{{\det\left({\left\\{{{\lambda_{i}}^{{N_{r}}-j}}\right\\}}\right)}}\sum\limits_{\vartheta={N_{r}}}^{\infty}{\frac{{{{\left({-1}\right)}^{\vartheta-1}}}}{{\vartheta!}}{{\left({\frac{{{{\left[{{{\bf{R}}_{t^{\prime}}}^{-1}}\right]}_{mm}}}}{{\bar{\gamma}{\theta_{m,k}}\ell({d_{k}})}}}\right)}^{\vartheta}}}$
$\displaystyle\quad\times\det\left(\left\\{{{{{\lambda_{i}}}}^{{N_{r}}-\vartheta-1}},\left.{{{{{\lambda_{i}}}}^{{N_{r}}-j}}}\right|_{j=2}^{{N_{r}}}\right\\}\right)$
$\displaystyle=\frac{1}{{{N_{r}}!\det\left({{{\bf{R}}_{r}}}\right)}}{\left({\frac{{{{\left[{{{\bf{R}}_{t^{\prime}}}^{-1}}\right]}_{mm}}}}{{\bar{\gamma}{\theta_{m,k}}\ell({d_{k}})}}}\right)^{{N_{r}}}}+o\left({{{\bar{\gamma}}^{-{N_{r}}}}}\right),$
(21)

where the last step holds by using

$\det\left(\left\\{{{{{\lambda_{i}}}}^{-1}},\left.{{{{{\lambda_{i}}}}^{{N_{r}}-j}}}\right|_{j=2}^{{N_{r}}}\right\\}\right)=\\\
{\left({-1}\right)^{{N_{r}}-1}}\det{\left({{{\bf{R}}_{r}}}\right)^{-1}}\det\left({\left\\{{{\lambda_{i}}^{{N_{r}}-j}}\right\\}}\right).$
(22)

####  III-B2 $M>1$

Likewise, if $M>1$, by ignoring high order terms
$o\left({{{\bar{\gamma}}^{-\left({{N_{r}}-M+1}\right)}}}\right)$ in the right-
hand side of (III-B), $p_{m,k}^{out}$ is asymptotic to

$\displaystyle p_{m,k}^{out}$
$\displaystyle\simeq\frac{{{{\left({-1}\right)}^{{N_{r}}-M}}\left({{N_{r}}-1}\right)!}}{{\left({{N_{r}}-M}\right)!\left({{N_{r}}-M+1}\right)!\left({M-2}\right)!}}$
$\displaystyle\times\frac{\det\left(\left\\{{{{{\lambda_{i}}}}^{M-2}}\ln{{\lambda_{i}}},\left.{{{{{\lambda_{i}}}}^{{N_{r}}-j}}}\right|_{j=2}^{{N_{r}}}\right\\}\right)}{\det\left(\left\\{{{\lambda_{i}}^{{N_{r}}-j}}\right\\}\right)}$
$\displaystyle\times{\left({\frac{{{{\left[{{{\bf{R}}_{t^{\prime}}}^{-1}}\right]}_{mm}}}}{{\bar{\gamma}{\theta_{m,k}}\ell({d_{k}})}}}\right)^{{N_{r}}-M+1}}+o\left({{{\bar{\gamma}}^{-\left({{N_{r}}-M+1}\right)}}}\right).$
(23)

Clearly from both (III-B1) and (III-B2), the impacts of spatial correlation,
transmission rates, power allocation coefficients and the numbers of IoT
devices and data streams on the outage performance can be quantified.
Moreover, it is revealed that the diversity order achieved by the virtual
MIMO-NOMA system at fixed rates is $\mathpzc d=N_{r}-M+1$, where the diversity
order is defined as $\mathpzc
d=\mathop{\lim}\nolimits_{\bar{\gamma}\to\infty}{{\log{p_{m,k}^{out}}}}/{{\log\bar{\gamma}}}$.
The diversity order reflects the decreasing slope of the outage probability
relative to the transmit SNR. It is apparent that increasing the number of IoT
devices within the cluster is beneficial to the improvement in the reception
reliability. However, massive IoT devices in the cluster entails prohibitively
high energy consumption and operating expenditure of channel estimation and
signaling synchronization, which is unfavorable for IoT networks.

###  III-C Diversity-Multiplexing Tradeoff (DMT)

To study the DMT in MIMO-NOMA channels, we assume that the power allocation
coefficients scale as
${\zeta_{m,k}}\simeq{c_{m,k}}{{\bar{\gamma}}^{-{\upsilon_{m,k}}}}$, where
${c_{m,k}}$ is a constant and ${\upsilon_{m,k}}\geq 0$. In particular, the
power allocation scheme can be treated as a fixed one if ${\upsilon_{m,k}}=0$
for all $m\in[1,M]$ and $k\in[1,K]$. To reap the multiplexing benefit of MIMO-
NOMA channels, let the transmission rate for the data stream of $x_{m,k}$
scale as ${R_{m,k}}={r_{m,k}}{\log_{2}}\bar{\gamma}$, where ${r_{m,k}}\geq 0$
corresponds to multiplexing gain [61]. On the scale of interest, by following
the above settings together with ${\theta_{m,k}}>0$ and
$\sum\nolimits_{k=1}^{K}{{\zeta_{m,k}}}={1}/{M}$, the conditions of
${\upsilon_{m,1}}\geq\cdots\geq{\upsilon_{m,K}}=0$ and
${r_{m,i}}\leq{\upsilon_{m,i-1}}-{\upsilon_{m,i}}$ for $2\leq i\leq K$ must be
satisfied. Accordingly, on the basis of (III-B1) and (III-B2), the outage
probability $p_{m,k}^{out}$ at target transmission rates
${R_{m,1}},\cdots,{R_{m,K}}$ is proportional to

$\displaystyle p_{m,k}^{out}$
$\displaystyle\propto{\left({{{\bar{\gamma}{\theta_{m,k}}}}}\right)^{-({N_{r}}-M+1)}}$
$\displaystyle\propto{\left({{{\bar{\gamma}\min\left\\{{{c_{m,i}}{{\bar{\gamma}}^{-{\upsilon_{m,i}}-{r_{m,i}}}},k\leq
i\leq K}\right\\}}}}\right)^{-({N_{r}}-M+1)}}$
$\displaystyle\propto{{\bar{\gamma}}^{-\left({{N_{r}}-M+1}\right)\left({1-{\upsilon_{m,k}}-{r_{m,k}}}\right)}},$
(24)

where the last step holds by using the inequality
${r_{m,i}}+{\upsilon_{m,i}}\leq{\upsilon_{m,i-1}}\leq{r_{m,i-1}}+{\upsilon_{m,i-1}}$.
Therefore, (III-C) gives rise to a DMT for the data stream of $x_{m,k}$ as

$\displaystyle
d_{m,k}^{*}(\\{r_{m,k}\\}_{k\in[1,K]})=-\lim\limits_{\bar{\gamma}\to\infty}\frac{\log
p_{m,k}^{out}(\\{{r_{m,k}}{\log_{2}}\bar{\gamma}\\}_{k\in[1,K]})}{\log\bar{\gamma}}$
$\displaystyle=\left({{N_{r}}-M+1}\right)\left({1-{\upsilon_{m,k}}-{r_{m,k}}}\right),$
(25)

where $d_{m,k}^{*}(\\{r_{m,k}\\}_{k\in[1,K]})$ represents the diversity gain
for $x_{m,k}$, and the diversity gain at fixed rates typically refers to the
diversity order. By recalling
${r_{m,k}}+{\upsilon_{m,k}}\leq{r_{m,k-1}}+{\upsilon_{m,k-1}}$ and
${\upsilon_{m,K}}=0$, the diversity gains obey the following relationship

$\displaystyle d_{m,1}^{*}(\\{r_{m,k}\\}_{k\in[1,K]})\leq
d_{m,2}^{*}(\\{r_{m,k}\\}_{k\in[1,K]})\leq\cdots$ $\displaystyle\leq
d_{m,K}^{*}(\\{r_{m,k}\\}_{k\in[1,K]})=({{N_{r}}-M+1})(1-{r_{m,K}}).$ (26)

By taking the fixed power allocation scheme (i.e., ${\upsilon_{m,k}}=0$ for
all $k\in[1,K]$) as an example, the condition of
${r_{m,i}}\leq{\upsilon_{m,i-1}}-{\upsilon_{m,i}}$ for $2\leq i\leq K$ leads
to vanishing multiplexing gains, i.e., ${r_{m,2}}=\cdots={r_{m,K}}=0$. This
consequently yields the DMT of
$d_{m,1}^{*}(\\{r_{m,k}\\}_{k\in[1,K]})=\left({{N_{r}}-M+1}\right)\left({1-{r_{m,1}}}\right)$
and
$d_{m,2}^{*}(\\{r_{m,k}\\}_{k\in[1,K]})=\cdots=d_{m,K}^{*}(\\{r_{m,k}\\}_{k\in[1,K]})={{N_{r}}-M+1}$.

##  IV Goodput Maximization

The optimal design of MIMO-NOMA system has been intensively studied by
assuming perfect instantaneous CSI at the transmitter [15, 23]. However, this
assumption limits their widespread applications into resource-constrained IoT
scenarios, where only the statistical knowledge of CSI is usually known to the
transmitter. Nevertheless, the optimal system settings were seldom discussed
in the circumstances, because the optimal design using the exact outage
probability would be unmanageable. To overcome this drawback, we resort to the
tractable expression of the asymptotic outage probability to ease the system
optimization. To exemplify this, the goodput maximization is considered next.

Since reception failures are inevitable in the absence of perfect CSI in
practice, the goodput is a pivotal performance metric to evaluate the
achievable long term average throughput. It is clear that the goodput depends
on the preset system configurations, such as transmission rates and power
allocation coefficients. Accordingly, these adjustable parameters should be
optimally devised by utilizing the channel statistical information
systematically. In the sequel, the power allocation coefficients and/or
transmission rates are properly chosen to maximize the goodput for the
proposed scheme.

To enable the optimization with low computational complexity and close-form
solution, the manageable expression of the asymptotic outage probability is
adopted herein. For notational convenience, the asymptotic expression of the
outage probability is unified by combining (III-B1) and (III-B2) as

$p_{m,k}^{out}\simeq{\phi_{m,k}}{\theta_{m,k}}^{-\left({{N_{r}}-M+1}\right)},$
(27)

where ${\phi_{m,k}}$ is given by (28), shown at the top of the next page.

${\phi_{m,k}}=\left\\{{\begin{array}[]{*{20}{c}}{\frac{1}{{{N_{r}}!\det\left({{{\bf{R}}_{r}}}\right)}}{{\left({\frac{{{{\left[{{{\bf{R}}_{t^{\prime}}}^{-1}}\right]}_{mm}}}}{{\bar{\gamma}\ell({d_{k}})}}}\right)}^{{N_{r}}}},}&{M=1}\\\
{\frac{{{{\left({-1}\right)}^{{N_{r}}-M}}\left({{N_{r}}-1}\right)!\det\left(\left\\{{{{{\lambda_{i}}}}^{M-2}}\ln{{\lambda_{i}}},\left.{{{{{\lambda_{i}}}}^{{N_{r}}-j}}}\right|_{j=2}^{{N_{r}}}\right\\}\right)}}{{\left({{N_{r}}-M}\right)!\left({{N_{r}}-M+1}\right)!\left({M-2}\right)!\det\left({\left\\{{{\lambda_{i}}^{{N_{r}}-j}}\right\\}}\right)}}{{\left({\frac{{{{\left[{{{\bf{R}}_{t^{\prime}}}^{-1}}\right]}_{mm}}}}{{\bar{\gamma}\ell({d_{k}})}}}\right)}^{{N_{r}}-M+1}},}&{M>1}\end{array}}\right.,$
(28)

###  IV-A Optimal Power Allocation

Given the transmission rates, the problem of goodput maximization via power
allocation is casted as

$\begin{array}[]{*{20}{c
l}}{\mathop{\max}\limits_{{\left\\{{{\zeta_{m,k}}}\right\\}_{m\in[1,M]\hfill\atop
k\in\left[{1,K}\right]\hfill}}}}&{{T_{g}}}\\\
{{\rm{s}}{\rm{.t}}{\rm{.}}}&{\sum\nolimits_{k=1}^{K}{{\zeta_{m,k}}}=\frac{1}{M},\,m\in[1,M],}\\\
{}\hfil&{{\zeta_{m,k}}>0,m\in[1,M],k\in[1,K],}\\\
{}\hfil&{{{\log}_{2}}\left({1+\frac{{{\zeta_{m,k}}}}{{\sum\nolimits_{l=1}^{k-1}{{\zeta_{m,l}}}}}}\right)>{R_{m,k}}},\\\
{}\hfil&{\quad\quad\quad\quad m\in[1,M],k\in\left[{2,K}\right]},\\\
\end{array}$ (29)

where the first and last constraints are used to warrant the data streams’
fairness and NOMA transmissions, respectively.

By introducing the auxiliary variables, the optimization problem can be
reformulated as (30) at the top of the next page.

$\begin{array}[]{*{20}{c
l}}{\mathop{\max}\limits_{{\left\\{{{\zeta_{m,k}}}\right\\}_{m\in[1,M]\hfill\atop
k\in\left[{1,K}\right]\hfill}},{\left\\{{{\theta_{m,k}}}\right\\}_{m\in\left[{1,M}\right]\hfill\atop
k\in\left[{1,K}\right]\hfill}}}}&{\sum\nolimits_{m=1}^{M}{\sum\nolimits_{k=1}^{K}{\left({1-{\phi_{m,k}}{\theta_{m,k}}^{-\left({{N_{r}}-M+1}\right)}}\right){R_{m,k}}}}}\\\
{{\rm{s}}{\rm{.t}}{\rm{.}}}&{\sum\nolimits_{k=1}^{K}{{\zeta_{m,k}}}=\frac{1}{M},\,m\in[1,M],}\\\
{}\hfil&{{\zeta_{m,k}}>0,m\in[1,M],k\in[1,K],}\\\
{}\hfil&{{{\log}_{2}}\left({1+\frac{{{\zeta_{m,k}}}}{{\sum\nolimits_{l=1}^{k-1}{{\zeta_{m,l}}}}}}\right)>{R_{m,k}},{\mkern
1.0mu}k\in\left[{2,K}\right],}\\\
{}\hfil&{{\theta_{m,k}}=\min\left\\{{\left.{\frac{{{\zeta_{m,i}}}}{{{2^{{R_{m,i}}}}-1}}-\sum\nolimits_{l=1}^{i-1}{{\zeta_{m,l}}}}\right|k\leq
i\leq K}\right\\}}.\end{array}$ (30)

Notice that the objective function is an increasing function of
$\theta_{m,k}$, the maximum goodput is achieved at the upper bound of
$\theta_{m,k}$. The last constraint can thus be rewritten as
${{\theta_{m,k}}\leq{{{\zeta_{m,i}}}}/{({{2^{{R_{m,i}}}}-1})}-\sum\nolimits_{l=1}^{i-1}{{\zeta_{m,l}}}}$
equivalently. As a consequence, (30) can be converted into a convex
optimization problem with linear constraints, as shown in (31) at the top of
the following page.

$\begin{array}[]{*{20}{c
l}}{\mathop{\min}\limits_{{\left\\{{{\zeta_{m,k}}}\right\\}_{m\in[1,M]\hfill\atop
k\in\left[{1,K}\right]\hfill}},{\left\\{{{\theta_{m,k}}}\right\\}_{m\in\left[{1,M}\right]\hfill\atop
k\in\left[{1,K}\right]\hfill}}}}&{\sum\nolimits_{m=1}^{M}{\sum\nolimits_{k=1}^{K}{{\phi_{m,k}}{R_{m,k}}{\theta_{m,k}}^{-\left({{N_{r}}-M+1}\right)}}}}\\\
{{\rm{s}}{\rm{.t}}{\rm{.}}}&{\sum\nolimits_{k=1}^{K}{{\zeta_{m,k}}}=\frac{1}{M},\,m\in[1,M],}\\\
{}\hfil&{{\zeta_{m,k}}>0,m\in[1,M],k\in[1,K],}\\\
{}\hfil&{{\theta_{m,k}}>0,m\in[1,M],k\in[1,K],}\\\
{}\hfil&{{\theta_{m,k}}\leq\frac{{{\zeta_{m,i}}}}{{{2^{{R_{m,i}}}}-1}}-\sum\nolimits_{l=1}^{i-1}{{\zeta_{m,l}}},m\in[1,M],k\in\left[{1,K}\right],i\in\left[{k,K}\right]}.\end{array}$
(31)

Alternatively, the optimization problem can be solved by using KKT conditions.
The optimal solution to (31) is given by the following theorem.

######  Theorem 1.

The optimal power allocation coefficients
${{{{\boldsymbol{\zeta}}}}_{m}^{*}}={\left({{{\zeta}_{m,1}^{*}},\cdots,{{\zeta}_{m,K}^{*}}}\right)^{\rm{T}}}$
are given by

${{{\boldsymbol{\zeta}}}_{m}^{*}}=\frac{{{{\bf{L}}_{m}}^{-1}{{\bf{b}}_{m}}}}{{M{{\bf{1}}_{K}}^{\rm{T}}{{\bf{L}}_{m}}^{-1}{{\bf{b}}_{m}}}},$
(32)

where ${\bf 1}_{K}$ is an $K\times 1$ vector of ones, and

${{\bf{b}}_{m}}=\left({{{\left({\frac{{\left({{N_{r}}-M+1}\right){\phi_{m,1}}{R_{m,1}}}}{{{\bf
e}_{1}}^{\rm
T}{{\bf{U}}_{m}}^{-1}{{\bf{1}}_{K}}}}\right)}^{\frac{1}{{{N_{r}}-M+2}}}},\cdots,}\right.\\\
{\left.{{{\left({\frac{{\left({{N_{r}}-M+1}\right){\phi_{m,K}}{R_{m,K}}}}{{{\bf
e}_{K}}^{\rm
T}{{\bf{U}}_{m}}^{-1}{{\bf{1}}_{K}}}}\right)}^{\frac{1}{{{N_{r}}-M+2}}}}}\right)^{\rm{T}}},$
(33) $\displaystyle{{\bf{U}}_{m}}$ $\displaystyle={{\bf{L}}_{m}}^{\rm{T}}$
$\displaystyle=\left({\begin{array}[]{*{20}{c}}{\frac{1}{{{2^{{R_{m,1}}}}-1}}}&{-1}&\cdots&{-1}\\\
0&{\frac{1}{{{2^{{R_{m,2}}}}-1}}}&\cdots&{-1}\\\
\vdots&\vdots&\ddots&\vdots\\\
0&0&\cdots&{\frac{1}{{{2^{{R_{m,K}}}}-1}}}\end{array}}\right),$ (38)

and ${\bf e}_{k}$ denotes a column vector with one as its $k$-th element and
zeros elsewhere. Moreover, ${\theta_{m,k}^{*}}$ is given by

$\displaystyle{\theta_{m,k}^{*}}$
$\displaystyle=\frac{{\zeta_{m,k}^{*}}}{{{2^{{R_{m,k}}}}-1}}-\sum\limits_{l=1}^{k-1}{\zeta_{m,l}^{*}}=\frac{{{{\bf{a}}_{m,k}}^{\rm{T}}{{\bf{L}}_{m}}^{-1}{{\bf{b}}_{m}}}}{{M{{\bf{1}}_{K}}^{\rm{T}}{{\bf{L}}_{m}}^{-1}{{\bf{b}}_{m}}}},$
(39)

where
${{\bf{a}}_{m,k}}={\left({\underbrace{-1,\cdots,-1}_{k-1},\frac{1}{{{2^{{R_{m,k}}}}-1}},\underbrace{0,\cdots,0}_{K-k}}\right)}^{\rm{T}}$.

###### Proof.

Please refer to Appendix E. ∎

By putting (39) into (27), $p_{m,k}^{out}$ can be obtained as

$p_{m,k}^{out}\simeq{M^{{N_{r}}-M+1}}{\phi_{m,k}}{\left({\frac{{{{\bf{1}}_{K}}^{\rm{T}}{{\bf{L}}_{m}}^{-1}{{\bf{b}}_{m}}}}{{{{\bf{a}}_{m,k}}^{\rm{T}}{{\bf{L}}_{m}}^{-1}{{\bf{b}}_{m}}}}}\right)^{{N_{r}}-M+1}}.$
(40)

Then plugging (40) into (7) yields the asymptotic goodput.

Based on Theorem 1, we arrive at the following remark regarding the power
allocation coefficients.

######  Remark 1.

The coefficients $\theta_{m,1}^{*},\cdots,\theta_{m,K}^{*}$ follow the
ascending order of their indices such that
$\theta_{m,1}^{*}<\cdots<\theta_{m,K}^{*}$, which results in

$\frac{{\zeta_{m,1}^{*}}}{{{2^{{R_{m,1}}}}-1}}<\frac{{\zeta_{m,2}^{*}}}{{{2^{{R_{m,2}}}}-1}}<\cdots<\frac{{\zeta_{m,K}^{*}}}{{{2^{{R_{m,K}}}}-1}},$
(41)

where ${{\zeta_{m,k}^{*}}}/{({{2^{{R_{m,k}}}}-1})}$ is the ratio of the power
allocation coefficient to the received SNR threshold, which characterizes the
cluster fairness. Accordingly, the result in (41) indicates that more
transmission power is allocated to cluster under worse channel condition. This
demonstrates that NOMA scheme can maintain the fairness among clusters while
achieving the maximum goodput. On the other hand, this is consistent with the
result of the DMT studied in Section III-C if we substitute the settings of
${\zeta_{m,k}}\simeq{c_{m,k}}{{\bar{\gamma}}^{-{\upsilon_{m,k}}}}$ and
${R_{m,k}}={r_{m,k}}{\log_{2}}\bar{\gamma}$ into (41).

###### Proof.

Please refer to Appendix F. ∎

######  Remark 2.

The asymptotic outage probability is lower bounded as

$p_{m,k}^{out\\_asy}\geq{M^{{N_{r}}-M+1}}{\phi_{m,k}}{\left({{2^{{R_{m,k}}}}-1}\right)^{{N_{r}}-M+1}}.$
(42)

The right hand side of the inequality (42) is equal to the asymptotic outage
probability of cluster $k$ under OMA transmissions, in which all the power is
allocated to cluster $k$’s data stream $m$.

###### Proof.

Please refer to Appendix G. ∎

Roughly speaking from Remark 2, the outage probability is an increasing
function of transmission rate $R_{m,k}$. With the definition of the goodput,
the transmission rate has two opposite effects on the goodput. Specifically,
on one hand, the increase of transmission rate indicates more information bits
delivered to the cluster. On the other hand, it results in high outage
probability. Accordingly, the goodput approaches to zero no matter when
$R_{m,k}$ tends to zero or infinity. All in all, the transmission rate should
be well designed to attain the maximum goodput. However, by putting (40) into
(7), the complex expression of the goodput makes it fairly impossible to
obtain the optimal transmission rate in closed-form, and even intractable to
numerically solve the problem. Hence, a suboptimal algorithm is developed for
the joint optimization of transmission powers and rates, which is deferred to
Section IV-C.

###  IV-B Optimal Rate Selection

Furthermore, for fixed values of the power allocation coefficients
${\zeta_{m,k}}$, the problem of goodput maximization through optimal rate
selection while implementing NOMA transmissions is formulated as

$\begin{array}[]{*{20}{c
l}}{\mathop{\max}\limits_{{\left\\{{{R_{m,k}}}\right\\}_{m\in[1,M]\hfill\atop
k\in\left[{1,K}\right]\hfill}}}}&{{T_{g}}}\\\
{{\rm{s}}{\rm{.t}}{\rm{.}}}&{{R_{m,k}}>0,m\in[1,M],k\in[1,K],}\\\
{}\hfil&{{{\log}_{2}}\left({1+\frac{{{\zeta_{m,k}}}}{{\sum\nolimits_{l=1}^{k-1}{{\zeta_{m,l}}}}}}\right)>{R_{m,k}}},\\\
{}\hfil&{\quad\quad\quad\quad m\in[1,M],k\in\left[{2,K}\right]}.\end{array}$
(43)

Analogously, by introducing the auxiliary variables $\theta_{m,k}$, (43) can
be further reformulated as (44), shown at the top of the next page.

$\begin{array}[]{*{20}{c}}{\mathop{\max}\limits_{{{\left\\{{{R_{m,k}}}\right\\}}_{m\in[1,M]\hfill\atop
k\in\left[{1,K}\right]\hfill}},{{\left\\{{{\theta_{m,k}}}\right\\}}_{m\in\left[{1,M}\right]\hfill\atop
k\in\left[{1,K}\right]\hfill}}}}&{\sum\nolimits_{m\in\left[{1,M}\right]\hfill\atop
k\in\left[{1,K}\right]\hfill}{\left({1-{\phi_{m,k}}{\theta_{m,k}}^{-\left({{N_{r}}-M+1}\right)}}\right){R_{m,k}}}}\\\
{{\rm{s}}.{\rm{t}}.}&{{R_{m,k}}>0,m\in[1,M],k\in[1,K],}\\\
{}\hfil&{{\theta_{m,k}}>0,m\in[1,M],k\in[1,K],}\\\
{}\hfil&{{R_{m,i}}\leq{{\log}_{2}}\left({1+\frac{{{\zeta_{m,i}}}}{{{\theta_{m,k}}+\sum\nolimits_{l=1}^{i-1}{{\zeta_{m,l}}}}}}\right),m\in[1,M],k\in\left[{1,K}\right],i\in\left[{k,K}\right].}\end{array}$
(44)

By means of KKT conditions, the optimal solution to (44) is obtained in the
following theorem.

######  Theorem 2.

The optimal transmission rate is given by

$R_{m,k}^{*}={\log_{2}}\left({1+\frac{{{\zeta_{m,k}}}}{{\theta_{m,k}^{*}+\sum\nolimits_{l=1}^{k-1}{{\zeta_{m,l}}}}}}\right),$
(45)

where $\theta_{m,k}^{*}$ is the zero point of $\varrho(x)-\ell(x)=0$, and

$\varrho(x)=\frac{{{\zeta_{m,k}}x^{{N_{r}}-M+2}\left({1-\phi_{m,k}x^{-\left({{N_{r}}-M+1}\right)}}\right)}}{{\left({{N_{r}}-M+1}\right)\phi_{m,k}{{\ln}}\left({1+\frac{{{\zeta_{m,k}}}}{{x+\sum\nolimits_{l=1}^{k-1}{{\zeta_{m,l}}}}}}\right)}},$
(46)
$\ell(x)=\left({x+\sum\nolimits_{l=1}^{k-1}{{\zeta_{m,l}}}}\right)\left({x+\sum\nolimits_{l=1}^{k}{{\zeta_{m,l}}}}\right).$
(47)

Besides, the optimal rate satisfies the following relationship

$\frac{{\zeta_{m,1}}}{{{2^{{R_{m,1}^{*}}}}-1}}<\frac{{\zeta_{m,2}}}{{{2^{{R_{m,2}^{*}}}}-1}}<\cdots<\frac{{\zeta_{m,K}}}{{{2^{{R_{m,K}^{*}}}}-1}},$
(48)

which is in perfect accordance with the analysis of the DMT in Section III-C.

###### Proof.

Please see Appendix H. ∎

###  IV-C Joint Power Allocation and Rate Selection

At last, we consider the maximization of goodput by jointly optimize power
allocation coefficients and transmission rates. Unfortunately, it is virtually
impossible to derive closed-form solution for this joint optimization problem.
It has been elaborated in Section IV-A that the non-linear fractional form of
the outage probability in (40) hinders further optimization of transmission
rate. Nevertheless, the power allocation coefficients and transmission rates
can be sub-optimally designed by using alternately iterating optimization,
which is based on both the results in Sections IV-A and IV-B [62]. More
specifically, by fixing the transmission rates, the power allocation
coefficients are optimized with (32). Then by fixing the optimal power
allocation, the optimal transmission rate can be obtained with (45). We
successively conduct the above two steps, until a negligible difference
between the optimal objective values of two consecutive iterations occurs. The
pseudocode of the joint power and rate optimization is provided in Algorithm
1. The local convergence of the proposed suboptimal algorithm can be
elucidated as follows. For ease of exposition, we denote by
${T_{g}}(\left(R_{m,k}\right)_{M\times
K},\left({{{\zeta_{m,k}}}}\right)_{M\times K})$ the goodput corresponding to
the fixed transmissions rates $\left(R_{m,k}\right)_{M\times K}$ and powers
$\left({{{\zeta_{m,k}}}}\right)_{M\times K}$. According to Algorithm 1, the
optimal values of the goodput computed by any two consecutive iterative steps
are in an increasing order as

$\displaystyle{T_{g}}\left(\left(R_{m,k}^{(\imath)}\right)_{M\times
K},\left({{{\zeta_{m,k}^{(\imath)}}}}\right)_{M\times K}\right)$
$\displaystyle\leq{T_{g}}\left(\left(R_{m,k}^{(\imath)}\right)_{M\times
K},\left({{{\zeta_{m,k}^{(\imath+1)}}}}\right)_{M\times K}\right)$
$\displaystyle\leq{T_{g}}\left(\left(R_{m,k}^{(\imath+1)}\right)_{M\times
K},\left({{{\zeta_{m,k}^{(\imath+1)}}}}\right)_{M\times K}\right),\,\imath\geq
0,$ (49)

where $\imath$ denotes the iteration number, and the inequalities hold thanks
to the global optimality of both the power allocation and rate selection
algorithms developed in Section IV-A and IV-B, respectively. Moreover, by
putting (42) into (7), $T_{g}$ is upper bounded, because

$\displaystyle{T_{g}}\leq$
$\displaystyle\sum\limits_{m=1}^{M}{\sum\limits_{k=1}^{K}{\left({1-{M^{{N_{r}}-M+1}}{\phi_{m,k}}{{\left({{2^{{R_{m,k}}}}-1}\right)}^{{N_{r}}-M+1}}}\right){R_{m,k}}}},$
(50)

where the right-hand side of (IV-C) is obviously upper bounded owing to its
concavity with respect to $({R_{m,k}})_{M\times K}$. Accordingly, the
monotonic and bounded sequence of $\\{{T_{g}}((R_{m,k}^{(\imath)})_{M\times
K},({{{\zeta_{m,k}^{(\imath)}}}})_{M\times K})\\}$ is convergent. This
consequently substantiates the convergence of Algorithm 1.

Algorithm 1 Suboptimal Joint Power Allocation and Rate Selection Algorithm

1:Generate the initial rates $\left(R_{m,k}^{(0)}\right)_{M\times K}$

2:Set $\imath=0$ and ${T_{g}^{*}}^{(\imath)}=-1$

3:repeat

4: $\imath=\imath+1$

5: Given $\left(R_{m,k}^{(\imath-1)}\right)_{M\times K}$, update
$\left({{{\zeta_{m,k}^{(\imath)}}}}\right)_{M\times K}$ with (32)

6: Given $\left({{{\zeta_{m,k}^{(\imath)}}}}\right)_{M\times K}$, update
$\left(R_{m,k}^{(\imath)}\right)_{M\times K}$ with (45)

7: ${T_{g}^{*}}^{(\imath)}={T_{g}}((R_{m,k}^{(\imath)})_{M\times
K},({{{\zeta_{m,k}^{(\imath)}}}})_{M\times K})$

8:until satisfying termination criterion
${T_{g}^{*}}^{(\imath)}-{T_{g}^{*}}^{(\imath-1)}<\epsilon$

9:$\left(\zeta_{m,k}^{*}\right)_{M\times
K}=\left(\zeta_{m,k}^{(\imath)}\right)_{M\times K}$

10:$\left(R_{m,k}^{*}\right)_{M\times
K}=\left(R_{m,k}^{(\imath+1)}\right)_{M\times K}$

From both (41) and (48), the proposed sub-optimal algorithm would prefer to
allocate more transmission power to clusters under bad channel conditions and
support them with higher transmission rate under good channel conditions.

##  V Numerical Results and Discussions

In this section, numerical results are presented for verifications and
discussions. For illustration, the parameters of the virtual MIMO-NOMA system
are set as follows unless otherwise specified. All the data streams are
assumed to have the same transmission rate $R_{m,i}$, i.e., $R_{m,i}=R$.
According to the prerequisite stated in (III-A), the power allocation
coefficients can be devised as [28] to implement NOMA such that

$\displaystyle{\zeta_{m,k}}=\left({\frac{1}{M}-\sum\limits_{l=k+1}^{K}{{\zeta_{m,l}}}}\right)\left({1-\varepsilon{2^{-R}}}\right),k=2,3,...,K,$
$\displaystyle{\zeta_{m,1}}=\frac{1}{M}-\sum\limits_{l=2}^{K}{{\zeta_{m,l}}},$
(51)

where $\varepsilon$ lies within $\left[{0,1}\right]$. The transmit beamforming
matrix $\bf V$ is set with ones on its principal diagonal and zeros elsewhere.
The antenna correlation is modelled by using the exponential correlation
structure as ${\bf R}_{t}=\left([\rho^{|i-j|}]_{1\leq i,j\leq N_{t}}\right)$
and ${\bf R}_{r}=\left([\rho^{|i-j|}]_{1\leq i,j\leq N_{r}}\right)$ [30]. The
system parameters $N_{t}$, $N_{r}$, $M$, $\varepsilon$, $\alpha$,
$\mathcal{K}$, $\rho$, $\bar{\gamma}$ and $R$ are assumed to be $3$, $3$, $3$,
$0.7$, 3, 1, $0.5$, $70$dB and $2$bps/Hz, respectively. We assume that the
base station is placed at the centre of the cell and $K=4$ clusters are
served. The cluster centers are located at positions $(10{\rm m},0)$,
$(0,20{\rm m})$, $(0,-30{\rm m})$ and $(-40{\rm m},0)$.

###  V-A Numerical Verifications

Fig. 2 depicts the outage probability $p_{1,1}^{out}$ against the transmit SNR
$\bar{\gamma}$. As expected from Fig. 2, the Monte Carlo simulation results
confirm the validity of the analytical results obtained from (10), and there
is a perfect match between the analytical and asymptotic results in high SNR
regime. Besides, it is also shown in Fig. 2 that the decreasing slope of the
outage curve increases with the number of data streams $M$, but is independent
of the path loss exponent $\alpha$. This is consistent with the analysis of
the diversity order, i.e., $\mathpzc d=N_{r}-M+1$. Whereas, since $\alpha$
characterizes the severity of path loss, the decrease of $\alpha$ is favorable
for decreasing the outage probability. Moreover, the outage probability
increases with $M$ given a fixed value of $\bar{\gamma}$, because the total
power constraint results in the reduction of the transmit power allocated to
each data stream as $M$ increases.

Figure 2: The outage probability $p_{1,1}^{out}$ versus the transmit SNR
$\bar{\gamma}$.

Fig. 3 plots the goodput $T_{g}$ versus the transmit SNR $\bar{\gamma}$. It is
observed that the theoretical results coincide well with the simulation ones.
It is clear from Fig. 3 that less severe fading, i.e., small $\alpha$,
improves the goodput. However, unlike Fig. 2, the virtual MIMO-NOMA system may
benefit from the increase of $M$ in terms of the goodput, albeit degrading the
outage performance. Furthermore, it can be seen from Fig. 3 that the goodput
is upper bounded as SNR increases because of the prescribed target
transmission rates $R_{m,k}$. Hence, the observations imply the necessity of
optimally designing the transmission rate to maximize the goodput, which will
be illustrated in Section V-D.

Figure 3: The goodput $T_{g}$ versus the transmit SNR $\bar{\gamma}$.

###  V-B Impact of Spatial Correlation

The effects of the spatial correlation on the outage probability and goodput
are investigated in Figs. 4 and 5, respectively. From both figures, the
theoretical and simulation results are in perfect agreement under low spatial
correlation due to the assumption made in (A). It can also be observed from
Figs. 4 and 5 that the spatial correlation adversely impacts the outage
probability and goodput. In particular, the spatial correlation would severely
impair the performance of the virtual MIMO-NOMA system if $\rho>0.5$, which
conforms to the conclusion drawn in [63, Fig. 3.5].

Figure 4: The impact of the spatial correlation on the outage probability
$p_{1,1}^{out}$. Figure 5: The impact of the spatial correlation on the
goodput $T_{g}$.

###  V-C Comparison with Five Baseline Schemes

To illustrate the advantages of the proposed scheme, we choose five reference
schemes for benchmarking, including the NOMA scheme (labelled as “NOMA”) [64],
the conventional OMA scheme (labelled as “OMA”) that serves as a baseline in
[23], the virtual MIMO scheme (labelled as “Virtual MIMO”) [29], the hybrid MA
scheme (labelled as “Hybrid MA”) in [9], the virtual MIMO-hybrid MA scheme
(labelled as “Virtual MIMO-Hybrid MA”). Herein, the virtual MIMO-hybrid MA
scheme is a modified version of the proposed scheme, which groups virtual MIMO
entities into pairs to conduct NOMA and then accommodates different pairs with
OMA. More specifically, we assume for the NOMA, OMA and hybrid MA schemes that
each IoT device requests one data stream from the BS. With regard to both the
proposed and virtual MIMO-based schemes, this amounts to the case of three
data streams requested by each cluster, i.e., $M=3$. For fair comparison,
time-division multiple access (TDMA) is considered for the OMA, virtual MIMO,
hybrid MA and virtual MIMO-hybrid MA schemes to serve all the IoT
devices/clusters/pairs with equal-length time slots. In addition, for both the
NOMA and hybrid MA schemes, each IoT device can be treated as a special
cluster with only a single IoT device. It is worthwhile to note that our
analytical results also apply to the virtual MIMO-hybrid MA scheme as long as
user pairing is determined. This is because each NOMA pair can be tackled
independently. For illustration, an extreme pairing policy in [23] is
considered for both hybrid MA-based schemes, where the device/virtual MIMO
entity having the best channel condition is paired with the one having the
worst channel condition, and the rest of devices/entities are paired in the
same manner.

Fig. 6 compares the goodputs of the proposed scheme and five baseline schemes.
It is obviously observed that the proposed scheme performs much better than
the five baseline schemes in high SNR regime. However, as opposed to both the
OMA and virtual MIMO schemes, the superior performance of the proposed scheme
cannot be perfectly exhibited at low SNR. This is due to the fact that the
significant advantages of NOMA are based on the premise of high SNR[43].
Whereas, this result no longer holds at low SNR. It is worth noting that the
condition of high SNR is commonly satisfied to achieve a very low outage
(e.g., $10^{-5}$) in 5G systems [65]. For instance, by considering a carrier
bandwidth of 100MHz, thermal noise level of -174dBm/Hz and maximum output
power budget of 24dBm, the transmit SNR is calculated as $118$dB [66]. It is
also observed from Fig. 6 that the proposed scheme generally outperforms the
NOMA scheme in terms of the goodput. For example, by aiming at
$T_{g}=6$bps/Hz, the required transmit SNRs for the proposed scheme and NOMA
scheme must be set as 65dB and 77dB, respectively. Henceforth, the
introduction of the virtual MIMO yields a reduction of the transmit SNR by
12dB, which remarkably conserves energy with guaranteed QoS for IoT networks.
In addition, it is shown that the proposed scheme is a little bit inferior to
the NOMA scheme under low SNR because of the negative effect of antenna
correlations. The similar observations can be found coincidentally between the
OMA and virtual MIMO schemes. Besides, it can be seen that both hybrid MA-
based schemes achieve a higher goodput than the proposed scheme at low SNR.
Nevertheless, the maximum achievable goodput of the hybrid MA scheme in [9] is
tremendously restricted. In contrast, the virtual MIMO-hybrid MA scheme is
able to strike a balance between boosting spectral efficiency and reducing
system overhead for resource-limited IoT networks.

Figure 6: The comparison of the goodput of the proposed scheme and five
baseline schemes with $N_{t}=N_{r}=M=3$, $\varepsilon=0.7$, $\rho=0.5$ and
$R=2$bps/Hz.

###  V-D Goodput Maximization

In this section, Figs. 7 and 8 are plotted to demonstrate the necessity of the
goodput maximization through power allocation and/or rate selection. For
illustrative purposes, we assume that the base station only serves two
clusters with coordinates $(10{\rm m},0)$ and $(0,20{\rm m})$, and a single
data stream is delivered to each cluster, i.e., $M=1$. Accordingly, if the
values of the transmission rates $R_{1,1}$ and $R_{1,2}$ are fixed, the power
allocation coefficients $\zeta_{1,1}$ and $\zeta_{1,2}$ can be optimally
chosen as studied in Section IV-A. With the optimal power allocation, a three-
dimensional figure is displayed in Fig. 7 to describe the relationship between
the maximum goodput and the transmission rates. It is clearly observed that
the transmission rates drastically influences the maximum achievable goodput.
Moreover, in order to reach the maximum goodput, Fig. 7 exhibits that the near
cluster is inclined to have a higher transmission rate than the distant
cluster. It therefore justifies the result in Section IV-C that the proposed
joint optimization tends to provide the cluster under better channel condition
with a higher data rate. Whereas, the difference between the optimal
transmission rates diminishes with the increase of the transmit SNR.

Figure 7: The optimal goodput versus the transmission rates.

With the optimal rate selection, Fig. 8 plots the maximum goodput versus
$\zeta_{1,1}$, where $\zeta_{1,2}=1-\zeta_{1,1}$. It can seen from Fig. 8 that
the power allocation coefficients should be properly selected to attain the
maximum achievable goodput. In particular, if $\zeta_{1,1}=0$ and
$\zeta_{1,2}=1$, NOMA scheme degenerates to its counterpart, i.e., OMA.
Apparently, neither $\zeta_{1,1}=0$ nor $\zeta_{1,2}=1$ will yield the maximum
goodput. Conversely, the optimal power allocation coefficient often lies in
between $[0,1]$. In addition, it is shown that more transmission power is
allocated to the farther user. It further substantiates the result in Section
IV-C that the joint power and rate optimization prefers to assign more
transmission power to the cluster under worse channel condition. The
difference between the optimal values of the power allocation coefficients
becomes conspicuous in high SNR regime. This is because more transmission
power is demanded to support the increasing data rate for the cluster under
poor channel condition at high SNR, as elucidated in Fig. 7.

Figure 8: The optimal goodput versus the power allocation coefficient.

Fig. 9 testifies the convergence behavior of the suboptimal joint power
allocation and rate selection algorithm, namely Algorithm 1. As seen from Fig.
9, the monotonic convergence of Algorithm 1 validates the correctness of the
convergence analysis in Section IV-C. Moreover, it is shown that the goodput
performance almost converges to the upper bound within 15 iterations, which
demonstrates the effectiveness of Algorithm 1. In addition, it is as expected
that the optimal goodput grows with the transmit SNR $\bar{\gamma}$.

Figure 9: The convergence behavior of Algorithm 1 with the maximum allowable
stopping criteria $\epsilon=10^{-5}$.

##  VI Conclusion

This paper has proposed a downlink virtual MIMO-NOMA scheme in IoT networks.
We have assumed only statistical CSI available at the base station to avoid
frequent instantaneous SINR reporting and signaling overhead so as to save the
energy and bandwidth for IoT networks. The outage probability and goodput of
the virtual MIMO-NOMA system have been thoroughly examined by using Kronecker
correlation model. The asymptotic results not only have provided physical
insights but also have paved a way for the goodput maximization. In
particular, the asymptotic outage expressions quantify the impacts of various
system parameters and make it possible to examine the DMT of MIMO-NOMA
channels. Additionally, the power allocation coefficients and/or transmission
rates have been properly chosen to reach the maximum goodput. By virtue of KKT
conditions, the optimal solutions to the goodput maximization problems have
been derived in closed-form. Besides, it has been found that the proposed
algorithms tend to allocate more power to clusters under unfavorable channel
conditions and support clusters with higher transmission rate under benign
channel conditions.

##  Appendix A Proof of (10)

In what follows we specialize to the derivation of $p_{m=1,k}^{out}$ first.
The outage probability for arbitrary data stream $m$ can be got in an
analogous way. By partitioning ${\bf{Z}}$ as
${\bf{Z}}=[{{{\bf{z}}_{1}}},{{\bf{\tilde{Z}}}}]$ along with matrix inversion
in block form in [67, Exercise 5.16], $p_{m{{=}}1,k}^{out}$ can be further
expressed as

$\displaystyle p_{m=1,k}^{out}=$
$\displaystyle\quad\Pr\left[{{{\bf{z}}_{1}}^{\rm{H}}\left({{\bf{I}}-{\bf{\tilde{Z}}}{{\left({{{{\bf{\tilde{Z}}}}^{\rm{H}}}{\bf{\tilde{Z}}}}\right)}^{-1}}{{{\bf{\tilde{Z}}}}^{\rm{H}}}}\right){{\bf{z}}_{1}}<\frac{1}{{{{\bar{\gamma}}}{\theta_{1,k}}\ell({d_{k}})}}}\right].$
(52)

In analogous to [58, eq. 15], $p_{m{{=}}1,k}^{out}$ can be derived as (53),
shown at the top of the next page,

$p_{m{{=}}1,k}^{out}=\Pr\left[{{{\bf{u}}^{\rm{H}}}{{\bf{R}}_{r}}^{1/2}\underbrace{\left({{\bf{I}}-{{\bf{R}}_{r}}^{1/2}{{{\bf{\tilde{Z}}}}_{w}}{{\left({{{{\bf{\tilde{Z}}}}_{w}}^{\rm{H}}{{\bf{R}}_{r}}{{{\bf{\tilde{Z}}}}_{w}}}\right)}^{-1}}{{{\bf{\tilde{Z}}}}_{w}}^{\rm{H}}{{\bf{R}}_{r}}^{1/2}}\right)}_{\bf{Q}}{{\bf{R}}_{r}}^{1/2}{\bf{u}}<\frac{{{{\left[{{{\bf{R}}_{t^{\prime}}}^{-1}}\right]}_{11}}}}{{{{\bar{\gamma}}}{\theta_{1,k}}\ell({d_{k}})}}}\right],$
(53)

where ${{{\bf{\tilde{Z}}}}_{w}}\in{\mathbb{C}}^{N_{r}\times(M-1)}$ and ${\bf
u}\in{\mathbb{C}}^{N_{r}\times 1}$ have i.i.d. zero-mean complex Gaussian
elements with unity variance. By defining
${X_{k}}={{\bf{u}}^{\rm{H}}}{{\bf{R}}_{r}}^{1/2}{\bf{Q}}{{\bf{R}}_{r}}^{1/2}{\bf{u}}$,
it follows that

$p_{m{=}1,k}^{out}=F_{{{X_{k}}}}\left({{{{\left[{{{\bf{R}}_{t^{\prime}}}^{-1}}\right]}_{11}}}}/{({{{\bar{\gamma}}}{\theta_{1,k}}\ell({d_{k}})})}\right),$
(54)

where $F_{{{X_{k}}}}\left(x\right)$ is the cumulative distribution function
(CDF) of $X_{k}$. Conditioned on
${{\bf{R}}_{r}}^{1/2}{\bf{Q}}{{\bf{R}}_{r}}^{1/2}$, we arrive at

$\displaystyle
F_{{{X_{k}}}}\left(x\right)={{\mathbb{E}}_{{{{\bf{\tilde{Z}}}}_{w}}}}\left\\{{\Pr\left[{\left.{{X_{k}}<x}\right|{{\bf{R}}_{r}}^{1/2}{\bf{Q}}{{\bf{R}}_{r}}^{1/2}}\right]}\right\\}$
$\displaystyle={{\mathbb{E}}_{{{{\bf{\tilde{Z}}}}_{w}}}}\left\\{{{F_{\left.{{X_{k}}}\right|{{\bf{R}}_{r}}^{1/2}{\bf{Q}}{{\bf{R}}_{r}}^{1/2}}}\left({\left.x\right|{{\bf{R}}_{r}}^{1/2}{\bf{Q}}{{\bf{R}}_{r}}^{1/2}}\right)}\right\\},$
(55)

where
${F_{\left.{{X_{k}}}\right|{{\bf{R}}_{r}}^{1/2}{\bf{Q}}{{\bf{R}}_{r}}^{1/2}}}\left(x\right)$
denotes the conditional CDF of $X_{k}$ given
${{\bf{R}}_{r}}^{1/2}{\bf{Q}}{{\bf{R}}_{r}}^{1/2}$. To proceed, we recourse to
the method of MGF to derive
${F_{\left.{{X_{k}}}\right|{{\bf{R}}_{r}}^{1/2}{\bf{Q}}{{\bf{R}}_{r}}^{1/2}}}\left(x\right)$.

Given ${{{\bf{R}}_{r}}^{1/2}{\bf{Q}}{{\bf{R}}_{r}}^{1/2}}$, $X_{k}$ is a
Hermitian quadratic form in complex normal variables, and the MGF of $X_{k}$
conditioned on ${{{\bf{R}}_{r}}^{1/2}{\bf{Q}}{{\bf{R}}_{r}}^{1/2}}$ can be
obtained by using Turin’s result from [68] as

$\displaystyle{\mathcal{M}_{\left.{{X_{k}}}\right|{{\bf{R}}_{r}}^{1/2}{\bf{Q}}{{\bf{R}}_{r}}^{1/2}}}\left(s\right)={{\det\left({{\bf{I}}-s{{\bf{R}}_{r}}^{1/2}{\bf{Q}}{{\bf{R}}_{r}}^{1/2}}\right)}}^{-1}.$
(56)

Then, the CDF
${F_{\left.{{X_{k}}}\right|{{\bf{R}}_{r}}^{1/2}{\bf{Q}}{{\bf{R}}_{r}}^{1/2}}}\left(y\right)$
can be derived by invoking inverse Laplace transform as[69, eq.3.6.8]

${F_{\left.{{X_{k}}}\right|{{\bf{R}}_{r}}^{1/2}{\bf{Q}}{{\bf{R}}_{r}}^{1/2}}}\left(x\right)\\\
={{\cal
L}^{-1}}\left\\{\frac{{{\mathcal{M}_{\left.{{X_{k}}}\right|{{\bf{R}}_{r}}^{1/2}{\bf{Q}}{{\bf{R}}_{r}}^{1/2}}}}\left({-s}\right)}{s}\right\\}\left(x\right).$
(57)

By substituting (56) into (57), it follows that

$\displaystyle{F_{\left.{{X_{k}}}\right|{{\bf{R}}_{r}}^{1/2}{\bf{Q}}{{\bf{R}}_{r}}^{1/2}}}\left(x\right)$
$\displaystyle\quad=\frac{1}{{2\pi{\rm i}}}\int\nolimits_{a-{\rm
i}\infty}^{a+{\rm
i}\infty}{\frac{{{e^{sx}}}}{{s\det\left({{\bf{I}}+s{{\bf{R}}_{r}}^{1/2}{\bf{Q}}{{\bf{R}}_{r}}^{1/2}}\right)}}ds}$
$\displaystyle\quad=\frac{1}{{2\pi{\rm i}}}\int\nolimits_{a-{\rm
i}\infty}^{a+{\rm
i}\infty}{\frac{{{e^{sx}}}}{s}{}_{1}{F_{0}}\left({1;-s{{\bf{R}}_{r}}^{1/2}{\bf{Q}}{{\bf{R}}_{r}}^{1/2}}\right)ds},$
(58)

where ${\rm i}=\sqrt{-1}$, the second step holds by using the identity
$\det{\left({{\bf{I}}-{\bf{X}}}\right)^{-a}}={}_{1}{F_{0}}\left({a;{\bf{X}}}\right)$
from [70], and ${}_{p}{F_{q}}\left({{\bf{X}}}\right)$ denotes the
hypergeometric function of one Hermitian matrix [56]. Hence, the CDF of
$X_{k}$, i.e., $F_{{{X_{k}}}}\left(x\right)$, can be obtained by taking
expectation of (A) over the distribution of ${{{\bf{\tilde{Z}}}}_{w}}$ as

$\displaystyle
F_{{{X_{k}}}}\left(x\right)={{\mathbb{E}}_{{{{\bf{\tilde{Z}}}}_{w}}}}\left\\{{{F_{\left.{{X_{k}}}\right|{{\bf{R}}_{r}}^{1/2}{\bf{Q}}{{\bf{R}}_{r}}^{1/2}}}\left({\left.x\right|{{\bf{R}}_{r}}^{1/2}{\bf{Q}}{{\bf{R}}_{r}}^{1/2}}\right)}\right\\}$
$\displaystyle=\frac{1}{{2\pi{\rm i}}}\int\nolimits_{a-{\rm i}\infty}^{a+{\rm
i}\infty}{\frac{{{e^{sx}}}}{s}{{\rm{E}}_{{{{\bf{\tilde{Z}}}}_{w}}}}\left\\{{{}_{1}{F_{0}}\left({1;-s{{\bf{R}}_{r}}^{1/2}{\bf{Q}}{{\bf{R}}_{r}}^{1/2}}\right)}\right\\}ds}$
$\displaystyle=\frac{1}{{2\pi{\rm i}}}\int\nolimits_{a-{\rm i}\infty}^{a+{\rm
i}\infty}{\frac{{{e^{sx}}}}{s}{{\rm{E}}_{{{{\bf{\tilde{Z}}}}_{w}}}}\left\\{{{}_{1}{F_{0}}\left({1;-s{{\bf{R}}_{r}}{\bf{Q}}}\right)}\right\\}ds}.$
(59)

where the second step holds by exchanging the operations of the integration
and the expectation, and the last step holds by using
${}_{p}{F_{q}}\left({{\bf{AB}}}\right)={}_{p}{F_{q}}\left({{\bf{A}}^{1/2}{\bf{B}}{\bf{A}}^{1/2}}\right)$
for any positive definite matrix $\bf A$ and symmetric matrix $\bf B$ [56,
p237, eq. 17]. Noticing that $\bf Q$ is an idempotent matrix with $N_{r}-M+1$
eigenvalues of value 1 and all other eigenvalues of value 0. More
specifically, by defining ${\bf
Y}={{\bf{R}}_{r}}^{1/2}{{{\bf{\tilde{Z}}}}_{w}}$, $\bf Y$ can be written as
${\bf Y}={\bf{U\Sigma}}{{\bf{\Lambda}}^{\rm{H}}}$ based on singular value
decomposition, where ${\bf{\Sigma}}={\rm
diag}\left\\{\sqrt{{{\delta_{i}}}}\right\\}\in{\mathbb{R}}_{+}^{N_{r}\times(M-1)}$,
${\bf U}\in{{\cal U}\left({{N_{r}}}\right)}$, ${\bf\Lambda}\in{{\cal
U}\left({{M-1}}\right)}$ and ${{\cal U}\left({{N}}\right)}$ stands for the
group of unitary $N\times N$ matrices. Accordingly, $\bf Q$ can be rewritten
as

${\bf{Q}}={\bf{U}}\underbrace{\left({{\bf{I}}-{\bf{\Sigma}}{{\left({{{\bf{\Sigma}}^{\rm{H}}}{\bf{\Sigma}}}\right)}^{-1}}{{\bf{\Sigma}}^{\rm{H}}}}\right)}_{\bf{V}}{{\bf{U}}^{\rm{H}}},$
(60)

where ${\bf{V}}$ is a diagonal matrix as

${\bf{V}}={\rm{diag}}\left\\{{\underbrace{1,\cdots,1}_{{N_{r}}-M+1},\underbrace{0,\cdots,0}_{M-1}}\right\\}.$
(61)

Thus, (A) can be further obtained as

$\displaystyle F_{{{X_{k}}}}\left(x\right)=$
$\displaystyle\quad\frac{1}{{2\pi{\rm i}}}\int\nolimits_{a-{\rm
i}\infty}^{a+{\rm
i}\infty}{\frac{{{e^{sx}}}}{s}{{\rm{E}}_{\bf{U}}}\left\\{{{}_{1}{F_{0}}\left({1;-s{{\bf{R}}_{r}}{{\bf{U}}}{\bf{V}}{\bf
U}^{\rm{H}}}\right)}\right\\}ds}.$ (62)

To get (A), we have to obtain the distribution of random unitary matrix $\bf
U$. With the definition of $\bf Y$ and the singular value decomposition, we
arrive at (A), as shown at the top of the following page,

$\displaystyle\left({d{\bf{Y}}}\right)$
$\displaystyle=\det{\left({{{\bf{R}}_{r}}}\right)^{M-1}}\left({d{{{\bf{\tilde{Z}}}}_{w}}}\right)$
$\displaystyle=\left({M-1}\right)!\prod\limits_{j=1}^{M-1}{\left({{N_{r}}-j}\right)!\left({M-1-j}\right)!}\prod\limits_{1\leq
i<j\leq
M-1}{{{\left({{\delta_{j}}-{\delta_{i}}}\right)}^{2}}}\prod\limits_{i=1}^{M-1}{{\delta_{i}}^{{N_{r}}-M+1}}d{\delta_{1}}\cdots
d{\delta_{M-1}}\left({d{\bf{U}}}\right)\left({d{\bf{\Lambda}}}\right),$ (63)

where the symbol $(d{\bf X})$ stands for the exterior product of the elements
of $d{\bf X}$[71, 56, 72]. In particular, $\left({d{\bf{U}}}\right)$ and
$\left({d{\bf{\Lambda}}}\right)$ denote the standard Haar measures of ${{\cal
U}\left({{N_{r}}}\right)}$ and ${{\cal U}\left({{M-1}}\right)}$, respectively.
Hence, by combining (A) with the probability density function (PDF) of
$\tilde{\bf Z}_{w}$, i.e.,
${f_{{{{\bf{\tilde{Z}}}}_{w}}}}({{{{\bf{\tilde{Z}}}}_{w}}})={{{\pi^{-{N_{r}}\left({M-1}\right)}}}}{\rm{etr}}({-{{{\bf{\tilde{Z}}}}_{w}}{{{\bf{\tilde{Z}}}}_{w}}^{\rm{H}}})$,
the PDF of $\bf U$ is derived as

${f_{\bf{U}}}\left({\bf{U}}\right)=\frac{{\left({M-1}\right)!\prod\nolimits_{j=1}^{M-1}{\left({{N_{r}}-j}\right)!\left({M-1-j}\right)!}}}{{\det{{\left({{{\bf{R}}_{r}}}\right)}^{M-1}}{\pi^{{N_{r}}\left({M-1}\right)}}}}\times\\\
\int\nolimits_{{\bf\Delta}\in{\mathbb{R}}_{+}^{{N_{r}}\times\left({M-1}\right)}}{{\rm{etr}}\left({-{{\bf{R}}_{r}}^{-1}{\bf{U}}{\bf\Delta}{{\bf{U}}^{\rm{H}}}}\right)\times}\\\
\prod\limits_{1\leq i<j\leq
M-1}{{{\left({{\delta_{j}}-{\delta_{i}}}\right)}^{2}}}\prod\limits_{j=1}^{M-1}{{\delta_{i}}^{{N_{r}}-M+1}}d{\delta_{1}}\cdots
d{\delta_{M-1}},$ (64)

where ${\bf{\Delta}}={\rm
diag}\left\\{{{{\delta_{i}}}}\right\\}\in{\mathbb{R}}_{+}^{N_{r}\times(M-1)}$
and ${\rm{etr}}({\bf X})={\rm{exp}}({\rm tr}\\{\bf X\\})$. Unfortunately, the
complex form of the distribution of $\bf U$ given by (64) will render the
derivation of (A) cumbersome. Motivated by the finding that
${f_{\bf{U}}}\left({\bf{U}}\right)=1$ if the receive antennas undergo
uncorrelated fadings, i.e., ${{{\bf{R}}_{r}}}={\bf I}_{N_{r}}$, (A) can be
approximated as

$\displaystyle F_{{{X_{k}}}}\left(x\right)$
$\displaystyle\approx\frac{1}{{2\pi{\rm i}}}\int\nolimits_{a-{\rm
i}\infty}^{a+{\rm i}\infty}\frac{{{e^{sx}}}}{s}$
$\displaystyle\times{\int\nolimits_{{\bf{U}}\in{\cal
U}\left({{N_{r}}}\right)}}{{}_{1}{F_{0}}\left({1;-s{{\bf{R}}_{r}}{{\bf{U}}}{\bf{V}}{\bf
U}^{\rm{H}}}\right)(d{\bf{U}})}ds$
$\displaystyle\triangleq\tilde{F}\left(x\right).$ (65)

Afterwards, by identifying the integral in (A) with [56, Theorem 7.3.3],
$\tilde{F}\left(x\right)$ can be further rewritten as

$\displaystyle\tilde{F}\left(x\right)$ $\displaystyle=\frac{1}{{2\pi{\rm
i}}}\int\nolimits_{a-{\rm i}\infty}^{a+{\rm i}\infty}{\frac{{{e^{sx}}}}{s}}$
$\displaystyle\times\int\nolimits_{{\bf{U}}\in{\cal
U}\left({{N_{r}}}\right)}{{}_{1}{F_{0}}\left({1;-s{{\bf{R}}_{r}}{{\bf{U}}}{\bf{V}}{\bf
U}^{\rm{H}}}\right)(d{\bf{U}})}ds$ $\displaystyle=\frac{1}{{2\pi{\rm
i}}}\int\nolimits_{a-{\rm i}\infty}^{a+{\rm
i}\infty}{\frac{{{e^{sx}}}}{s}{}_{1}{F_{0}}^{\left({{N_{r}}}\right)}\left({1;-s{{\bf{R}}_{r}},{\bf{V}}}\right)ds},$
(66)

where ${}_{p}{F_{q}}^{(m)}\left({{\bf{X}}},{\bf Y}\right)$ denotes the
hypergeometric function of two $m\times m$ Hermitian matrices [56]. To derive
(A) in closed-form,
${}_{1}{F_{0}}^{\left({{N_{r}}}\right)}\left({1;-s{{\bf{R}}_{r}},{\bf{V}}}\right)$
is simplified as (67), shown at the top of the next page,

$\displaystyle{}_{1}{F_{0}}^{\left({{N_{r}}}\right)}\left({1;-s{{\bf{R}}_{r}},{\bf{V}}}\right)$
$\displaystyle=\frac{\det\left({{{\left\\{{{{{{\lambda_{i}}}}^{{N_{r}}-1}}{{}_{2}{F_{1}}}\left({1,{N_{r}}-M+1;{N_{r}};-s{\lambda_{i}}}\right)}\right\\}}_{\scriptstyle
1\leq i\leq{N_{r}}\hfill\atop\scriptstyle
j=1\hfill}},{{\left\\{{{\lambda_{i}}^{{N_{r}}-j}}\right\\}}_{\scriptstyle
1\leq i\leq{N_{r}}\hfill\atop\scriptstyle 2\leq
j\leq{N_{r}}\hfill}}}\right)}{{\det\left({{{\left\\{{{\lambda_{i}}^{{N_{r}}-j}}\right\\}}_{1\leq
i,j\leq{N_{r}}}}}\right)}},$ (67)

where ${}_{p}{F_{q}}(a_{1},\cdots,a_{p};b_{1},\cdots,b_{q};x)$ denotes the
hypergeometric function, ${\lambda_{1}},\cdots,{\lambda_{N_{r}}}$ denote the
$N_{r}$ eigenvalues of ${{\bf{R}}_{r}}$, and the proof of (67) is relegated to
Appendix B.

Substituting (67) into (A) leads to (68), shown at the top of the next page.

$\displaystyle\tilde{F}\left(x\right)$
$\displaystyle=\frac{\det\left({{{\left\\{{{{{{\lambda_{i}}}}^{{N_{r}}-1}}\frac{1}{{2\pi{\rm{i}}}}\int\nolimits_{a-{\rm{i}}\infty}^{a+{\rm{i}}\infty}{\frac{{{e^{sx}}}}{s}{}_{2}{F_{1}}\left({1,{N_{r}}-M+1;{N_{r}};-s{\lambda_{i}}}\right)ds}}\right\\}}_{\scriptstyle
1\leq i\leq{N_{r}}\hfill\atop\scriptstyle
j=1\hfill}},{{\left\\{{{\lambda_{i}}^{{N_{r}}-j}}\right\\}}_{\scriptstyle
1\leq i\leq{N_{r}}\hfill\atop\scriptstyle 2\leq
j\leq{N_{r}}\hfill}}}\right)}{{\det\left({\left\\{{{\lambda_{i}}^{{N_{r}}-j}}\right\\}_{1\leq
i,j\leq{N_{r}}}}\right)}}.$ (68)

By using [59, eq.9.34.8], the Gauss hypergeometric function can be expressed
in terms of Meijer G-function as

${{}_{2}{F_{1}}\left({1,{N_{r}}-M+1;{N_{r}};-s{\lambda_{1}}}\right)}=\\\
\frac{{\Gamma\left({{N_{r}}}\right)}}{{\Gamma\left({{N_{r}}-M+1}\right)}}G_{2,2}^{1,2}\left({\left.{\begin{array}[]{*{20}{c}}{0,M-{N_{r}}}\\\
{0,1-{N_{r}}}\end{array}}\right|s{\lambda_{1}}}\right).$ (69)

where $\Gamma(x)$ denotes Gamma function. By inverting Laplace transform [60,
eq.07.34.22.0003.01], we reach

$\displaystyle\frac{1}{{2\pi{\rm i}}}\int\nolimits_{a-{\rm i}\infty}^{a+{\rm
i}\infty}{\frac{{{e^{sx}}}}{s}{}_{2}{F_{1}}\left({1,{N_{r}}-M+1;{N_{r}};-s{\lambda_{i}}}\right)ds}$
$\displaystyle=\frac{{\Gamma\left({{N_{r}}}\right)}}{{\Gamma\left({{N_{r}}-M+1}\right)}}G_{3,2}^{1,2}\left({\left.{\begin{array}[]{*{20}{c}}{0,M-{N_{r}},1}\\\
{0,1-{N_{r}}}\end{array}}\right|\frac{{{\lambda_{i}}}}{x}}\right)$ (72)
$\displaystyle=\frac{{\Gamma\left({{N_{r}}}\right)}}{{\Gamma\left({{N_{r}}-M+1}\right)}}G_{2,3}^{2,1}\left({\left.{\begin{array}[]{*{20}{c}}{1,{N_{r}}}\\\
{1,{N_{r}}-M+1,0}\end{array}}\right|\frac{x}{{{\lambda_{i}}}}}\right),$ (75)

where the last step holds by using [59, eq.9.31.2]. Plugging (A) into (68),
$\tilde{F}\left(x\right)$ can be finally derived as (11). Substituting (A) and
(11) into (54), the approximate expression of $p_{m=1,k}^{out}$ can be
derived. In the similar fashion, the outage probability for arbitrary data
stream $m$, $p_{m,k}^{out}$, can be obtained as (10).

##  Appendix B Proof of (67)

With [56, Definition 7.3.2],
${}_{1}{F_{0}}^{\left({{N_{r}}}\right)}\left({1;-s{{\bf{R}}_{r}},{\bf{V}}}\right)$
can be expanded as

${}_{1}{F_{0}}^{\left({{N_{r}}}\right)}\left({1;-s{{\bf{R}}_{r}},{\bf{V}}}\right)=\sum\limits_{k=0}^{\infty}{\sum\limits_{\boldsymbol{\kappa}}{{{[1]}_{\boldsymbol{\kappa}}}\frac{{{C_{\boldsymbol{\kappa}}}\left({-s{{\bf{R}}_{r}}}\right){C_{\boldsymbol{\kappa}}}\left({\bf{V}}\right)}}{{k!{C_{\boldsymbol{\kappa}}}\left({\bf{I}}\right)}}}},$
(76)

where ${\boldsymbol{\kappa}}=(k_{1},\cdots,k_{N_{r}})$ is a partition of an
integer $k$ with $k_{1}\geq k_{2}\geq\cdots k_{N_{r}}\geq 0$ and
$k_{1}+\cdots+k_{N_{r}}=k$, ${\left[a\right]_{\boldsymbol{\kappa}}}$ is the
complex multivariate hypergeometric coefficient and
${\left[a\right]_{\boldsymbol{\kappa}}}=\prod\limits_{i=1}^{m}{{{\left({a-i+1}\right)}_{{k_{i}}}}}$,
${\left(a\right)_{k}}=a\left({a+1}\right)\cdots\left({a+k-1}\right)$ and
${\left(a\right)_{0}}=1$. ${C_{\boldsymbol{\kappa}}}\left({\bf{X}}\right)$ is
defined as the complex zonal polynomial of ${\bf
X}\in{\mathbb{C}}^{N_{r}\times N_{r}}$ associated with
${\boldsymbol{\kappa}}$, and is given by [70, 73]

${C_{\boldsymbol{\kappa}}}\left({\bf{X}}\right)={\chi_{\left[{\boldsymbol{\kappa}}\right]}}\left(1\right){\chi_{\left[{\boldsymbol{\kappa}}\right]}}\left({\bf{X}}\right),$
(77)

where ${\chi_{\left[{\boldsymbol{\kappa}}\right]}}\left(1\right)$ is the
dimension of the representation $[{\boldsymbol{\kappa}}]$ of the symmetric
group on $k$ symbols given by

${\chi_{\left[{\boldsymbol{\kappa}}\right]}}\left(1\right)=k!\frac{{\prod\nolimits_{i<j}^{{N_{r}}}{\left({{k_{i}}-{k_{j}}-i+j}\right)}}}{{\prod\nolimits_{i=1}^{{N_{r}}}{\left({{k_{i}}+m-i}\right)!}}},$
(78)

and ${\chi_{\left[{\boldsymbol{\kappa}}\right]}}\left(\bf X\right)$ is the
character of the representation $\left[{\boldsymbol{\kappa}}\right]$ of the
linear group given as a symmetric function of the eigenvalues
$\mu_{1},\cdots,\mu_{N_{r}}$ of $\bf X$ by

${\chi_{\left[{\boldsymbol{\kappa}}\right]}}\left({\bf{X}}\right)=\frac{{\det\left({\left\\{{{\mu_{i}}^{{k_{j}}+{N_{r}}-j}}\right\\}}\right)}}{{\det\left({\left\\{{{\mu_{i}}^{{N_{r}}-j}}\right\\}}\right)}}.$
(79)

Particularly, ${C_{\boldsymbol{\kappa}}}\left({\bf{I}}\right)$ is given by[73]

${C_{\boldsymbol{\kappa}}}\left({\bf{I}}\right)=k!\frac{{{{\left[{\prod\nolimits_{i<j}^{{N_{r}}}{\left({{k_{i}}-{k_{j}}-i+j}\right)}}\right]}^{2}}}}{{\prod\nolimits_{i=1}^{{N_{r}}}{\left({{k_{i}}+{N_{r}}-i}\right)!}\prod\nolimits_{i=1}^{{N_{r}}}{\left({{N_{r}}-i}\right)!}}}.$
(80)

By using [74, eq. 5.11]

${[1]_{\boldsymbol{\kappa}}}=\left\\{{\begin{array}[]{*{20}{c}}{{k_{1}}!,}&{{k_{2}}=\cdots={k_{{N_{r}}}}=0}\\\
{0,}&{\rm else}\end{array}}\right.,$ (81)

(76) can be rewritten as

${}_{1}{F_{0}}^{\left({{N_{r}}}\right)}\left({1;-s{{\bf{R}}_{r}},{\bf{V}}}\right)=\sum\limits_{k=0}^{\infty}{\frac{{{C_{{{\boldsymbol{\kappa}}_{k}}}}\left({-s{{\bf{R}}_{r}}}\right){C_{{{\boldsymbol{\kappa}}_{k}}}}\left({\bf{V}}\right)}}{{{C_{{{\boldsymbol{\kappa}}_{k}}}}\left({\bf{I}}\right)}}},$
(82)

where ${{\boldsymbol{\kappa}}_{k}}=\left[{k,0,\cdots,0}\right]$. Moreover,
with (78), ${\chi_{\left[{{{\boldsymbol{\kappa}}_{k}}}\right]}}\left(1\right)$
can be obtained as
${\chi_{\left[{{{\boldsymbol{\kappa}}_{k}}}\right]}}\left(1\right)=1$. From
(80), we have
${C_{{{\boldsymbol{\kappa}}_{k}}}}\left({\bf{I}}\right)={\left({{N_{r}}}\right)_{k}}/{k!}$.
Accordingly, (82) collapses to

${}_{1}{F_{0}}^{\left({{N_{r}}}\right)}\left({1;-s{{\bf{R}}_{r}},{\bf{V}}}\right)=\sum\limits_{k=0}^{\infty}{\frac{{{\chi_{\left[{{{\boldsymbol{\kappa}}_{k}}}\right]}}\left({-s{{\bf{R}}_{r}}}\right){\chi_{\left[{{{\boldsymbol{\kappa}}_{k}}}\right]}}\left({\bf{V}}\right)}}{{{C_{{{\boldsymbol{\kappa}}_{k}}}}\left({\bf{I}}\right)}}}\\\
=\sum\limits_{k=0}^{\infty}{{{\left({-s}\right)}^{k}}\frac{{k!}{{\chi_{\left[{{{\boldsymbol{\kappa}}_{k}}}\right]}}\left({{{\bf{R}}_{r}}}\right){\chi_{\left[{{{\boldsymbol{\kappa}}_{k}}}\right]}}\left({\bf{V}}\right)}}{{\left({{N_{r}}}\right)_{k}}}},$
(83)

where the last step holds by using
${\chi_{\left[{\boldsymbol{\kappa}}\right]}}\left({a{\bf{X}}}\right)={a^{k}}{\chi_{\left[{\boldsymbol{\kappa}}\right]}}\left({\bf{X}}\right)$.
By using [71, Proposition 3],
${\chi_{\left[{{{\boldsymbol{\kappa}}_{k}}}\right]}}\left({\bf{V}}\right)$ can
be simplified as

$\displaystyle{\chi_{\left[{{{\boldsymbol{\kappa}}_{k}}}\right]}}\left({\bf{V}}\right)$
$\displaystyle={\chi_{\left[{{\boldsymbol{\kappa}}_{k}^{\prime}}\right]}}\left({{{\bf{V}}_{\left({1:{N_{r}}-M+1}\right)\times\left({1:{N_{r}}-M+1}\right)}}}\right)$
$\displaystyle=\frac{{{C_{{{\boldsymbol{\kappa}}_{k}^{\prime}}}}\left({\bf{I}}\right)}}{{{\chi_{\left[{{{\boldsymbol{\kappa}}_{k}^{\prime}}}\right]}}\left(1\right)}}=\frac{\left({{{N_{r}}-M+1}}\right)_{k}}{k!}.$
(84)

where the sub-matrix
${{{\bf{V}}_{\left({1:{N_{r}}-M+1}\right)\times\left({1:{N_{r}}-M+1}\right)}}}$
is formed by deleting the last $M-1$ rows and the last $M-1$ columns of
$\bf{V}$ and ${{{\boldsymbol{\kappa}}_{k}^{\prime}}}$ is an $(N_{r}-M+1)\times
1$ vector of ones. Substituting (B) into (83) gives

${}_{1}{F_{0}}^{\left({{N_{r}}}\right)}\left({1;-s{{\bf{R}}_{r}},{\bf{V}}}\right)=\sum\limits_{k=0}^{\infty}{{{\left({-s}\right)}^{k}}}\\\
\times\frac{{\left({{N_{r}}-M+1}\right)_{k}}}{{\left({{N_{r}}}\right)_{k}}}{\chi_{\left[{{{\boldsymbol{\kappa}}_{k}}}\right]}}\left({{{\bf{R}}_{r}}}\right)\\\
=\sum\limits_{k=0}^{\infty}{{{\left({-s}\right)}^{k}}\frac{{\left({{N_{r}}-M+1}\right)_{k}}}{{\left({{N_{r}}}\right)_{k}}}\frac{{\det\left({\left\\{{{\lambda_{i}}^{{k_{j}}+{N_{r}}-j}}\right\\}}\right)}}{{\det\left({\left\\{{{\lambda_{i}}^{{N_{r}}-j}}\right\\}}\right)}}},$
(85)

where $k_{1}=k$ and $k_{2}=\cdots=k_{N_{r}}=0$. After some algebraic
manipulations, (85) can be simplified as (86) at the top of the next page.

${}_{1}{F_{0}}^{\left({{N_{r}}}\right)}\left({1;-s{{\bf{R}}_{r}},{\bf{V}}}\right)=\frac{1}{{\det\left({\left\\{{{\lambda_{i}}^{{N_{r}}-j}}\right\\}}\right)}}\sum\limits_{k=0}^{\infty}{{{\left({-s}\right)}^{k}}\frac{{{{\left({{N_{r}}-M+1}\right)}_{k}}}}{{{{\left({{N_{r}}}\right)}_{k}}}}\det\left({{{\left\\{{{{\left({{\lambda_{i}}}\right)}^{k+{N_{r}}-1}}}\right\\}}_{\scriptstyle
1\leq i\leq{N_{r}}\hfill\atop\scriptstyle
j=1\hfill}},{{\left\\{{{\lambda_{i}}^{{N_{r}}-j}}\right\\}}_{\scriptstyle
1\leq i\leq{N_{r}}\hfill\atop\scriptstyle 2\leq
j\leq{N_{r}}\hfill}}}\right)}\\\
=\frac{1}{{\det\left({\left\\{{{\lambda_{i}}^{{N_{r}}-j}}\right\\}}\right)}}\det\left({{{\left\\{{{\lambda_{i}}^{{N_{r}}-1}\sum\limits_{k=0}^{\infty}{\frac{{{{\left({{N_{r}}-M+1}\right)}_{k}}}}{{{{\left({{N_{r}}}\right)}_{k}}}}}{{\left({-s{\lambda_{i}}}\right)}^{k}}}\right\\}}_{\scriptstyle
1\leq i\leq{N_{r}}\hfill\atop\scriptstyle
j=1\hfill}},{{\left\\{{{\lambda_{i}}^{{N_{r}}-j}}\right\\}}_{\scriptstyle
1\leq i\leq{N_{r}}\hfill\atop\scriptstyle 2\leq j\leq{N_{r}}\hfill}}}\right),$
(86)

With the definition of hypergeometric function, (86) can be further expressed
as (67).

##  Appendix C Proof of (19)

By virtue of the generalized L’Hôpital’s rule [75, Theorem 1.2.4],
$p_{m,k}^{out}$ can be obtained as

$\displaystyle
p_{m,k}^{out}=\frac{1}{{\Gamma\left({{N_{r}}-M+1}\right)}}{\frac{{{\partial^{{N_{r}}-1}}}}{{\partial{\lambda}^{{N_{r}}-1}}}\Biggl{[}{{\lambda}^{{N_{r}}-1}}\times}$
$\displaystyle\,\,\,\,\,{{G_{2,3}^{2,1}\left({\left.{\begin{array}[]{*{20}{c}}{1,{N_{r}}}\\\
{1,{N_{r}}-M+1,0}\end{array}}\right|\frac{{{{\left[{{{\bf{R}}_{t^{\prime}}}^{-1}}\right]}_{mm}}}}{{{\lambda}\bar{\gamma}{\theta_{m,k}}\ell({d_{k}})}}}\right)}\Biggr{]}_{{\lambda}={\lambda_{0}}}}.$
(89)

The $({{N_{r}}-1})$-order partial derivative of Meijer G-function in (C) can
be derived as (92), as shown at the top of the following page, wherein the
first and last steps hold by using [59, eq.9.31.5], and the second equality
holds by using [60, eq.07.34.20.0013.02].

$\displaystyle\eqref{eqn:deri_meijerg}$
$\displaystyle\coloneqq{\left({\frac{{{{\left[{{{\bf{R}}_{t^{\prime}}}^{-1}}\right]}_{mm}}}}{{\bar{\gamma}{\theta_{m,k}}\ell({d_{k}})}}}\right)^{{N_{r}}-1}}\frac{{{\partial^{{N_{r}}-1}}}}{{\partial{\lambda_{1}}^{{N_{r}}-1}}}\left[{G_{2,3}^{2,1}\left({\left.{\begin{array}[]{*{20}{c}}{2-{N_{r}},1}\\\
{2-{N_{r}},2-M,1-{N_{r}}}\end{array}}\right|\frac{{{{\left[{{{\bf{R}}_{t^{\prime}}}^{-1}}\right]}_{mm}}}}{{{\lambda_{1}}\bar{\gamma}{\theta_{m,k}}\ell({d_{k}})}}}\right)}\right]$
(92)
$\displaystyle=\frac{{{\partial^{{N_{r}}-1}}}}{{\partial{z^{{N_{r}}-1}}}}{\left[{G_{2,3}^{2,1}\left({\left.{\begin{array}[]{*{20}{c}}{2-{N_{r}},1}\\\
{2-{N_{r}},2-M,1-{N_{r}}}\end{array}}\right|{z^{-1}}}\right)}\right]_{z=\frac{{{\lambda_{1}}{{\bar{\gamma}}}{\theta_{m,k}}\ell({d_{k}})}}{{{{\left[{{{\bf{R}}_{t^{\prime}}}^{-1}}\right]}_{mm}}}}}}$
(95)
$\displaystyle={\left({\frac{{{{\left[{{{\bf{R}}_{t^{\prime}}}^{-1}}\right]}_{mm}}}}{{{\lambda_{1}}{{\bar{\gamma}}}{\theta_{m,k}}\ell({d_{k}})}}}\right)^{{N_{r}}-1}}G_{1,2}^{1,1}\left({\left.{\begin{array}[]{*{20}{c}}{2-{N_{r}}}\\\
{2-M,1-{N_{r}}}\end{array}}\right|\frac{{{{\left[{{{\bf{R}}_{t^{\prime}}}^{-1}}\right]}_{mm}}}}{{{\lambda_{1}}{{\bar{\gamma}}}{\theta_{m,k}}\ell({d_{k}})}}}\right)$
(98) $\displaystyle=G_{1,2}^{1,1}\left({\left.{\begin{array}[]{*{20}{c}}1\\\
{{N_{r}}-M+1,0}\end{array}}\right|\frac{{{{\left[{{{\bf{R}}_{t^{\prime}}}^{-1}}\right]}_{mm}}}}{{{\lambda_{1}}{{\bar{\gamma}}}{\theta_{m,k}}\ell({d_{k}})}}}\right),$
(101)

Applying [60, eq.06.07.26.0006.01] to (92) then leads to

$\displaystyle\frac{{{\partial^{{N_{r}}-1}}{{{{\lambda}}}^{{N_{r}}-1}}G_{2,3}^{2,1}\left({\left.{\begin{array}[]{*{20}{c}}{1,{N_{r}}}\\\
{1,{N_{r}}-M+1,0}\end{array}}\right|\frac{{{{\left[{{{\bf{R}}_{t^{\prime}}}^{-1}}\right]}_{mm}}}}{{{\lambda}\bar{\gamma}{\theta_{m,k}}\ell({d_{k}})}}}\right)}}{{\partial{\lambda}^{{N_{r}}-1}}}$
(104)
$\displaystyle={{\gamma\left({{N_{r}}-M+1,\frac{{{{\left[{{{\bf{R}}_{t^{\prime}}}^{-1}}\right]}_{mm}}}}{{{\lambda}{{\bar{\gamma}}}{\theta_{m,k}}\ell({d_{k}})}}}\right)}},$
(105)

Substituting (104) into (C) eventually gives (19).

##  Appendix D Proof of (III-B)

To start, the asymptotic expression of Meijer G-function
$G_{2,3}^{2,1}\left({\cdot|z}\right)$ for $z\to 0$ is derived first. By using
[60, eq.07.34.06.0045.01], the Meijer G-function
$G_{2,3}^{2,1}\left({\cdot|z}\right)$ can be written as

$\displaystyle
G_{2,3}^{2,1}\left({\left.{\begin{array}[]{*{20}{c}}{1,{N_{r}}}\\\
{1,{N_{r}}-M+1,0}\end{array}}\right|z}\right)=$ (108)
$\displaystyle-\sum\limits_{\vartheta=1}^{\infty}{{\rm{Res}}\left\\{{\frac{{\Gamma\left({1-\theta}\right)\Gamma\left({{N_{r}}-M+1-\theta}\right)\Gamma\left(\theta\right)}}{{\Gamma\left({{N_{r}}-\theta}\right)\Gamma\left({1+\theta}\right)}}{z^{\theta}},\theta=\vartheta}\right\\}},$
(109)

where ${\rm Res}(\textperiodcentered;v)$ stands for the residue at $s=v$.
Excluding the second order poles at $\\{N_{r}-M+1,\cdots,N_{r}-1\\}$, all the
other poles are simple. By using residue theorem, (108) can be further
expanded as

$\displaystyle
G_{2,3}^{2,1}\left({\left.{\begin{array}[]{*{20}{c}}{1,{N_{r}}}\\\
{1,{N_{r}}-M+1,0}\end{array}}\right|z}\right)=$ (112)
$\displaystyle\sum\limits_{\vartheta=1}^{{N_{r}}-M}{\frac{{{{\left({-1}\right)}^{\vartheta-1}}\left({{N_{r}}-M-\vartheta}\right)!}}{{\vartheta!\left({{N_{r}}-1-\vartheta}\right)!}}{z^{\vartheta}}}+\sum\limits_{\vartheta={N_{r}}-M+1}^{{N_{r}}-1}{\varphi\left(\vartheta\right){z^{\vartheta}}}$
$\displaystyle+\sum\limits_{\vartheta={N_{r}}-M+1}^{{N_{r}}-1}{\frac{{{{\left({-1}\right)}^{2\vartheta-{N_{r}}+M-1}}{z^{\vartheta}}\ln
z}}{{\vartheta!\left({\vartheta-\left({{N_{r}}-M+1}\right)}\right)!\left({{N_{r}}-1-\vartheta}\right)!}}}$
$\displaystyle+\sum\limits_{\vartheta={N_{r}}}^{\infty}{\frac{{{{\left({-1}\right)}^{\vartheta+M-2}}\left({\vartheta-{N_{r}}}\right)!}}{{\vartheta!\left({\vartheta-{N_{r}}+M-1}\right)!}}{z^{\vartheta}}},$
(113)

where $\varphi\left(\vartheta\right)$ is given by

$\displaystyle\varphi\left(\vartheta\right)={\left({-1}\right)^{2\vartheta-{N_{r}}+M-1}}\times$
$\displaystyle{\left[{\frac{d}{{d\theta}}\frac{{{{\left({\Gamma\left({\vartheta+1-\theta}\right)}\right)}^{2}}}}{{\prod\limits_{j=0}^{\vartheta-1}{\left({\theta-j}\right)}\prod\limits_{j={N_{r}}-M+1}^{\vartheta-1}{\left({\theta-j}\right)}\Gamma\left({{N_{r}}-\theta}\right)}}}\right]_{\theta=\vartheta}}.$
(114)

Substituting (112) into (11) yields (D) at the top of the next page.

$\displaystyle\tilde{F}\left(x\right)$
$\displaystyle=\frac{{\Gamma\left({{N_{r}}}\right)}}{{\Gamma\left({{N_{r}}-M+1}\right)\det\left({\left\\{{{\lambda_{i}}^{{N_{r}}-j}}\right\\}}\right)}}\sum\limits_{\vartheta=1}^{{N_{r}}-M}{\frac{{{{\left({-1}\right)}^{\vartheta-1}}\left({{N_{r}}-M-\vartheta}\right)!}}{{\vartheta!\left({{N_{r}}-1-\vartheta}\right)!}}{{x}^{\vartheta}}}\det\left(\left\\{{{{{\lambda_{i}}}}^{{N_{r}}-\vartheta-1}},\left.{{{{{\lambda_{i}}}}^{{N_{r}}-j}}}\right|_{j=2}^{{N_{r}}}\right\\}\right)$
$\displaystyle+\frac{{\Gamma\left({{N_{r}}}\right)}}{{\Gamma\left({{N_{r}}-M+1}\right)\det\left({\left\\{{{\lambda_{i}}^{{N_{r}}-j}}\right\\}}\right)}}\sum\limits_{\vartheta={N_{r}}-M+1}^{{N_{r}}-1}{\varphi\left(\vartheta\right){{x}^{\vartheta}}}\det\left(\left\\{{{{{\lambda_{i}}}}^{{N_{r}}-\vartheta-1}},\left.{{{{{\lambda_{i}}}}^{{N_{r}}-j}}}\right|_{j=2}^{{N_{r}}}\right\\}\right)$
$\displaystyle+\frac{{\Gamma\left({{N_{r}}}\right)}}{{\Gamma\left({{N_{r}}-M+1}\right)\det\left({\left\\{{{\lambda_{i}}^{{N_{r}}-j}}\right\\}}\right)}}\sum\limits_{\vartheta={N_{r}}-M+1}^{{N_{r}}-1}{\frac{{{{\left({-1}\right)}^{2\vartheta-{N_{r}}+M-1}}}{{x}^{\vartheta}}}{{\vartheta!\left({\vartheta-\left({{N_{r}}-M+1}\right)}\right)!\left({{N_{r}}-1-\vartheta}\right)!}}}$
$\displaystyle\times\det\left(\left\\{{{{{\lambda_{i}}}}^{{N_{r}}-\vartheta-1}}\ln\left({\frac{x}{{{\lambda_{i}}}}}\right),\left.{{{{{\lambda_{i}}}}^{{N_{r}}-j}}}\right|_{j=2}^{{N_{r}}}\right\\}\right)$
$\displaystyle+\frac{{\Gamma\left({{N_{r}}}\right)}}{{\Gamma\left({{N_{r}}-M+1}\right)\det\left({\left\\{{{\lambda_{i}}^{{N_{r}}-j}}\right\\}}\right)}}\sum\limits_{\vartheta={N_{r}}}^{\infty}{\frac{{{{\left({-1}\right)}^{\vartheta+M-2}}\left({\vartheta-{N_{r}}}\right)!}}{{\vartheta!\left({\vartheta-{N_{r}}+M-1}\right)!}}{x^{\vartheta}}}\det\left(\left\\{{{{{\lambda_{i}}}}^{{N_{r}}-\vartheta-1}},\left.{{{{{\lambda_{i}}}}^{{N_{r}}-j}}}\right|_{j=2}^{{N_{r}}}\right\\}\right).$
(115)

Since
$\det\left(\left\\{{{{{\lambda_{i}}}}^{\nu}},\left.{{{{{\lambda_{i}}}}^{{N_{r}}-j}}}\right|_{j=2}^{{N_{r}}}\right\\}\right)=0$
for $0\leq\nu\leq{N_{r}}-2$, (D) can be further rearranged as (III-B).

##  Appendix E Proof of Theorem 1

With regard to (31), the associated Lagrangian is written as (E) at the top of
the following page,

$\displaystyle{\cal
L}\left({{{\left\\{{{\zeta_{m,k}}}\right\\}}_{m\in\left[{1,M}\right]\hfill\atop
k\in\left[{1,K}\right]\hfill}},{{\left\\{{{\theta_{m,k}}}\right\\}}_{m\in\left[{1,M}\right]\hfill\atop
k\in\left[{1,K}\right]\hfill}},{{\left\\{{{\eta_{m}}}\right\\}}_{m\in\left[{1,M}\right]}},{{\left\\{{{\xi_{m,k}}}\right\\}}_{m\in\left[{1,M}\right]\hfill\atop
k\in\left[{1,K}\right]\hfill}},{{\left\\{{{\varsigma_{m,k}}}\right\\}}_{m\in[1,M]\hfill\atop
k\in[1,K]\hfill}},{{\left\\{{{\omega_{m,k,i}}}\right\\}}_{m\in[1,M]\hfill\atop{k\in\left[{1,K}\right]\hfill\atop
i\in\left[{k,K}\right]\hfill}}}}\right)$
$\displaystyle=\sum\nolimits_{m\in\left[{1,M}\right]\hfill\atop
k\in\left[{1,K}\right]\hfill}{{\phi_{m,k}}{R_{m,k}}{\theta_{m,k}}^{-\left({{N_{r}}-M+1}\right)}}+\sum\nolimits_{m\in\left[{1,M}\right]}{{\eta_{m}}\left({\sum\nolimits_{k=1}^{K}{{\zeta_{m,k}}}-\frac{1}{M}}\right)}-\sum\nolimits_{m\in\left[{1,M}\right]\hfill\atop
k\in\left[{1,K}\right]\hfill}{{\xi_{m,k}}{\zeta_{m,k}}}$
$\displaystyle\quad-\sum\nolimits_{m\in\left[{1,M}\right]\hfill\atop
k\in\left[{1,K}\right]\hfill}{{\varsigma_{m,k}}{\theta_{m,k}}}+\sum\nolimits_{m\in[1,M]\hfill\atop{k\in\left[{1,K}\right]\hfill\atop
i\in\left[{k,K}\right]\hfill}}{{\omega_{m,k,i}}\left({{\theta_{m,k}}-\frac{{{\zeta_{m,i}}}}{{{2^{{R_{m,i}}}}-1}}+\sum\nolimits_{l=1}^{i-1}{{\zeta_{m,l}}}}\right)},$
(116)

where ${{\eta_{m}}}$, ${{\xi_{m,k}}}$, ${{\varsigma_{m,k}}}$ and
${{\omega_{m,k,i}}}$ refer to the KKT multipliers. On the basis of (E), the
optimal solution should satisfy the necessary KKT conditions (117)-(122) as
follows

${\xi_{m,k}^{*}}{\zeta_{m,k}^{*}}=0,$ (117)
${\varsigma_{m,k}^{*}}{\theta_{m,k}^{*}}=0,$ (118)
${\omega_{m,k,i}^{*}}\left({{\theta_{m,k}^{*}}-\frac{{{\zeta_{m,i}^{*}}}}{{{2^{{R_{m,i}}}}-1}}+\sum\nolimits_{l=1}^{i-1}{{\zeta_{m,l}^{*}}}}\right)=0,$
(119) $\sum\nolimits_{k=1}^{K}{{\zeta_{m,k}^{*}}}=\frac{1}{M},$ (120)

where (121) and (122) are shown at the top of the next page,
${\xi_{m,k}^{*}}\geq 0$, ${\varsigma_{m,k}^{*}}\geq 0$ and
${\omega_{m,k,i}^{*}}\geq 0$.

${\left.{\frac{{\partial{\cal
L}}}{{\partial{\zeta_{m,k}}}}}\right|_{{{\left\\{{\zeta_{m,k}^{*}}\right\\}}_{m\in\left[{1,M}\right]\hfill\atop
k\in\left[{1,K}\right]\hfill}},{{\left\\{{\theta_{m,k}^{*}}\right\\}}_{m\in\left[{1,M}\right]\hfill\atop
k\in\left[{1,K}\right]\hfill}},{{\left\\{{\eta_{m}^{*}}\right\\}}_{m\in\left[{1,M}\right]}},{{\left\\{{\xi_{m,k}^{*}}\right\\}}_{m\in\left[{1,M}\right]\hfill\atop
k\in\left[{1,K}\right]\hfill}},{{\left\\{{\varsigma_{m,k}^{*}}\right\\}}_{m\in[1,M]\hfill\atop
k\in[1,K]\hfill}},{{\left\\{{\omega_{m,k,i}^{*}}\right\\}}_{m\in[1,M]\hfill\atop{k\in\left[{1,K}\right]\hfill\atop
i\in\left[{k,K}\right]\hfill}}}}}=0,$ (121) ${\left.{\frac{{\partial{\cal
L}}}{{\partial{\theta_{m,k}}}}}\right|_{{{\left\\{{\zeta_{m,k}^{*}}\right\\}}_{m\in\left[{1,M}\right]\hfill\atop
k\in\left[{1,K}\right]\hfill}},{{\left\\{{\theta_{m,k}^{*}}\right\\}}_{m\in\left[{1,M}\right]\hfill\atop
k\in\left[{1,K}\right]\hfill}},{{\left\\{{\eta_{m}^{*}}\right\\}}_{m\in\left[{1,M}\right]}},{{\left\\{{\xi_{m,k}^{*}}\right\\}}_{m\in\left[{1,M}\right]\hfill\atop
k\in\left[{1,K}\right]\hfill}},{{\left\\{{\varsigma_{m,k}^{*}}\right\\}}_{m\in[1,M]\hfill\atop
k\in[1,K]\hfill}},{{\left\\{{\omega_{m,k,i}^{*}}\right\\}}_{m\in[1,M]\hfill\atop{k\in\left[{1,K}\right]\hfill\atop
i\in\left[{k,K}\right]\hfill}}}}}=0,$ (122)

By putting (E) into (121) and (122), respectively, we reach

$\eta_{m}^{*}-\xi_{m,k}^{*}-\frac{1}{{{2^{{R_{m,k}}}}-1}}\sum\nolimits_{t\in\left[{1,k}\right]}{\omega_{m,t,k}^{*}}\\\
+\sum\nolimits_{t\in\left[{k+1,K}\right]\hfill\atop
i\in\left[{t,K}\right]\hfill}{\omega_{m,t,i}^{*}}=0,$ (123)
$-\left({{N_{r}}-M+1}\right){\phi_{m,k}}{R_{m,k}}{\theta_{m,k}^{*}}^{-\left({{N_{r}}-M+2}\right)}\\\
-\varsigma_{m,k}^{*}+\sum\nolimits_{i\in\left[{k,K}\right]}{\omega_{m,k,i}^{*}}=0.$
(124)

Notice ${{\zeta_{m,k}}>0}$ and ${{\theta_{m,k}}>0}$, (117) and (118) suggest
$\xi_{m,k}^{*}=0$ and $\varsigma_{m,k}^{*}=0$, respectively. Thus, (123) and
(124) can be rewritten as

$\eta_{m}^{*}-\frac{1}{{{2^{{R_{m,k}}}}-1}}\sum\limits_{t\in\left[{1,k}\right]}{\omega_{m,t,k}^{*}}+\sum\limits_{t\in\left[{k+1,K}\right]\hfill\atop
i\in\left[{t,K}\right]\hfill}{\omega_{m,t,i}^{*}}=0,$ (125)
$-\left({{N_{r}}-M+1}\right){\phi_{m,k}}{R_{m,k}}{\theta_{m,k}^{*}}^{-\left({{N_{r}}-M+2}\right)}\\\
+\sum\nolimits_{i\in\left[{k,K}\right]}{\omega_{m,k,i}^{*}}=0.$ (126)

To proceed, we assume that
${i_{m,k}}=\mathop{\arg}\limits_{i}\min\left\\{{{{\zeta_{m,i}^{*}}}/{({{2^{{R_{m,i}}}}-1})}-\sum\nolimits_{l=1}^{i-1}{\zeta_{m,l}^{*}},i\in\left[{k,K}\right]}\right\\}$,
and the index ${i_{m,k}}$ uniquely exists for simplicity. This assumption will
be further demonastrated in Appendix F. It is readily proved that
${i_{m,1}}\leq\cdots\leq{i_{m,K}}$, and the uniqueness of ${i_{m,k}}$ means

$\theta_{m,k}^{*}<\frac{{\zeta_{m,i}^{*}}}{{{2^{{R_{m,i}}}}-1}}+\sum\nolimits_{l=1}^{i-1}{\zeta_{m,l}^{*}},\,i\in[k,K]\wedge
i\neq{i_{m,k}}.$ (127)

Accordingly, (119) indicates

$\omega_{m,k,i\neq{i_{m,k}}}^{*}=0.$ (128)

Besides, by setting $k=K$ in (126), we have

$\omega_{m,K,K}^{*}=\left({{N_{r}}-M+1}\right){\phi_{m,K}}{R_{m,K}}{\theta_{m,K}^{*}}^{-\left({{N_{r}}-M+2}\right)}.$
(129)

Since $\theta_{m,k}^{*}$ is upper bounded from (127), it follows from (129)
that $\omega_{m,K,K}^{*}>0$. Then setting $k=K$ in (125) leads to

$\displaystyle\eta_{m}^{*}$
$\displaystyle=\frac{1}{{{2^{{R_{m,K}}}}-1}}\sum\limits_{t\in\left[{1,K}\right]}{\omega_{m,t,K}^{*}}\geq\frac{1}{{{2^{{R_{m,K}}}}-1}}\omega_{m,K,K}^{*}>0.$
(130)

Combining (125) and (128) results in

$\eta_{m}^{*}=\frac{1}{{{2^{{R_{m,k}}}}-1}}\sum\limits_{t\in\left[{1,k}\right]}{\omega_{m,t,k}^{*}}-\sum\limits_{t\in\left[{k+1,K}\right]}{\omega_{m,t,{i_{m,t}}}^{*}}.$
(131)

By setting $k=1$, it follows from (131) that

$\eta_{m}^{*}=\frac{1}{{{2^{{R_{m,1}}}}-1}}\omega_{m,1,1}^{*}-\sum\limits_{t\in\left[{2,K}\right]}{\omega_{m,t,{i_{m,t}}}^{*}}.$
(132)

With (130), (132) indicates that $\omega_{m,1,1}^{*}>0$. According to (128),
we deduce that $i_{m,1}=1$, i.e., $\omega_{m,1,i\neq 1}^{*}=0$. Then setting
$k=2$ in (131), we have

$\eta_{m}^{*}=\frac{1}{{{2^{{R_{m,2}}}}-1}}\left({\omega_{m,1,2}^{*}+\omega_{m,2,2}^{*}}\right)-\sum\limits_{t\in\left[{3,K}\right]}{\omega_{m,t,{i_{m,t}}}^{*}}.$
(133)

Since $\omega_{m,1,2}^{*}=0$, (133) suggests that $\omega_{m,2,2}^{*}>0$.
Accordingly, $i_{m,2}=2$ and $\omega_{m,2,i\neq 2}^{*}=0$ follow from (128).
Similarly, by successively setting $k=3,\cdots,K$ in (131), we have
${i_{m,k}}=k$. Therefore, from (119), it follows that

$\theta_{m,k}^{*}=\frac{{\zeta_{m,k}^{*}}}{{{2^{{R_{m,k}}}}-1}}-\sum\limits_{l=1}^{k-1}{\zeta_{m,l}^{*}}.$
(134)

By combining (127) with (134), we arrive at
$\theta_{m,1}^{*}<\cdots<\theta_{m,K}^{*}$. Meanwhile, (131) and (126)
respectively reduce to

$\eta_{m}^{*}=\frac{1}{{{2^{{R_{m,k}}}}-1}}\omega_{m,k,k}^{*}-\sum\limits_{t\in\left[{k+1,K}\right]}{\omega_{m,t,t}^{*}},$
(135)
$-\left({{N_{r}}-M+1}\right){\phi_{m,k}}{R_{m,k}}{\theta_{m,k}^{*}}^{-\left({{N_{r}}-M+2}\right)}\\\
+\omega_{m,k,k}^{*}=0.$ (136)

By solving the system of equations constituted by (135) for $k=1,\cdots,K$,
${{{{\boldsymbol{\omega}}}}_{m}}={\left({{{\omega}_{m,1,1}^{*}},\cdots,{{\omega}_{m,K,K}^{*}}}\right)^{\rm{T}}}$
can be obtained as

${{{{\boldsymbol{\omega}}}}_{m}}=\eta_{m}^{*}{{\bf{U}}_{m}}^{-1}{{\bf{1}}_{K}},$
(137)

where ${\bf 1}_{K}$ is $K\times 1$ vector of ones and ${{\bf{U}}_{m}}$ are
defined in (1). We define ${\bf e}_{k}$ as a column vector with one as the
$k$-th element and zeros elsewhere, and ${{\omega}_{m,k,k}^{*}}$ can then be
represented by

${{\omega}_{m,k,k}^{*}}=\eta_{m}^{*}{{\bf e}_{k}}^{\rm
T}{{\bf{U}}_{m}}^{-1}{{\bf{1}}_{K}}.$ (138)

Substituting (134) and (138) into (136) gives rise to

$-\left({{N_{r}}-M+1}\right){\phi_{m,k}}{R_{m,k}}\times\\\
{\left({\frac{{\zeta_{m,k}^{*}}}{{{2^{{R_{m,k}}}}-1}}-\sum\limits_{l=1}^{k-1}{\zeta_{m,l}^{*}}}\right)^{-\left({{N_{r}}-M+2}\right)}}\\\
+\eta_{m}^{*}{{\bf e}_{k}}^{\rm T}{{\bf{U}}_{m}}^{-1}{{\bf{1}}_{K}}=0.$ (139)

By defining
${\tilde{\zeta}_{m,k}^{*}}={{\zeta}_{m,k}^{*}}{\eta_{m}^{*}}^{\frac{1}{{{N_{r}}-M+2}}}$,
(139) can be rewritten as

$\frac{{\tilde{\zeta}_{m,k}^{*}}}{{{2^{{R_{m,k}}}}-1}}-\sum\limits_{l=1}^{k-1}{\tilde{\zeta}_{m,l}^{*}}\\\
={\left({\frac{{\left({{N_{r}}-M+1}\right){\phi_{m,k}}{R_{m,k}}}}{{{\bf
e}_{k}}^{\rm
T}{{\bf{U}}_{m}}^{-1}{{\bf{1}}_{K}}}}\right)^{\frac{1}{{{N_{r}}-M+2}}}}.$
(140)

By solving the system of equations in (140),
${{{\tilde{\boldsymbol{\zeta}}}}_{m}}={\left({{{\tilde{\zeta}}_{m,1}^{*}},\cdots,{{\tilde{\zeta}}_{m,K}^{*}}}\right)^{\rm{T}}}$
can be derived as

${{{\tilde{\boldsymbol{\zeta}}}}_{m}}={{\bf{L}}_{m}}^{-1}{{\bf{b}}_{m}},$
(141)

where ${{\bf{b}}_{m}}$ and ${{\bf{L}}_{m}}$ are defined in (33) and (1),
respectively.

Substituting (141) into (120) along with the definition of
${{\tilde{\zeta}}_{m,k}^{*}}$, ${\eta_{m}^{*}}^{\frac{1}{{{N_{r}}-M+2}}}$ is
given by

${\eta_{m}^{*}}^{\frac{1}{{{N_{r}}-M+2}}}={{M{{\bf{1}}_{K}}^{\rm{T}}{{\bf{L}}_{m}}^{-1}{{\bf{b}}_{m}}}},$
(142)

Defining by
${{{{\boldsymbol{\zeta}}}}_{m}^{*}}={\left({{{\zeta}_{m,1}^{*}},\cdots,{{\zeta}_{m,K}^{*}}}\right)^{\rm{T}}}$
the vector of power allocation coefficients,
${{{{\boldsymbol{\zeta}}}}_{m}^{*}}$ is finally expressed as (32). With (134),
${\theta_{m,k}^{*}}$ can be expressed as (39).

##  Appendix F Proof of Remark 1

${{\bf{L}}_{m}}$ is a lower triangular matrix and its matrix inverse is still
a lower triangular matrix. Specifically, ${{\bf{L}}_{m}}^{-1}$ can be derived
as

${{\bf{L}}_{m}}^{-1}=\left({{{\left[{\prod\limits_{t\in\left\\{{i,j}\right\\}}{\left({{2^{{R_{m,t}}}}-1}\right)}\prod\limits_{t=j+1}^{i-1}{{2^{{R_{m,t}}}}}}\right]}_{K\geq
i\geq j\geq 1}}}\right),$ (143)

where all entries of ${{\bf{L}}_{m}}^{-1}$ with $i<j$ are zero. By using
(143), ${{\bf{L}}_{m}}^{-1}{{\bf{b}}_{m}}$ can be explicitly expressed as
(144) at the top of the next page,

${{\bf{L}}_{m}}^{-1}{{\bf{b}}_{m}}=\left({\begin{array}[]{*{20}{c}}\vdots\\\
{\underbrace{\sum\nolimits_{j=1}^{i}{\prod\nolimits_{t\in\left\\{{i,j}\right\\}}{\left({{2^{{R_{m,t}}}}-1}\right)}\prod\nolimits_{t=j+1}^{i-1}{{2^{{R_{m,t}}}}}}{{\left({\frac{{\left({{N_{r}}-M+1}\right){\phi_{m,j}}{R_{m,j}}}}{{{{\bf
e}_{j}}^{\rm
T}{{\bf{U}}_{m}}^{-1}{{\bf{1}}_{K}}}}}\right)}^{\frac{1}{{{N_{r}}-M+2}}}}}_{{\beta_{m,i}},1\leq
i\leq K}}\\\ \vdots\end{array}}\right),$ (144)

where $\beta_{m,i}>0$. With (144),
${{{\bf{1}}_{K}}^{\rm{T}}{{\bf{L}}_{m}}^{-1}{{\bf{b}}_{m}}}$ and
${{{\bf{a}}_{m,k}}^{\rm{T}}{{\bf{L}}_{m}}^{-1}{{\bf{b}}_{m}}}$ can be written
as

${{{\bf{1}}_{K}}^{\rm{T}}{{\bf{L}}_{m}}^{-1}{{\bf{b}}_{m}}}=\sum\limits_{i=1}^{K}{{\beta_{m,i}}},$
(145)
${{{\bf{a}}_{m,k}}^{\rm{T}}{{\bf{L}}_{m}}^{-1}{{\bf{b}}_{m}}}=\frac{{{\beta_{m,k}}}}{{{2^{{R_{m,k}}}}-1}}-\sum\limits_{i=1}^{k-1}{{\beta_{m,i}}}.$
(146)

In order to prove that the sequence $\theta_{m,1}^{*},\cdots,\theta_{m,K}^{*}$
is in an ascending order of $k$, i.e., $\theta_{m,k+1}^{*}>\theta_{m,k}^{*}$.
From (39), it suffices to show that

$\frac{{{\beta_{m,k+1}}}}{{{2^{{R_{m,k+1}}}}-1}}-\sum\limits_{i=1}^{k}{{\beta_{m,i}}}>\frac{{{\beta_{m,k}}}}{{{2^{{R_{m,k}}}}-1}}-\sum\limits_{i=1}^{k-1}{{\beta_{m,i}}}.$
(147)

This amounts to proving
${{{\beta_{m,k+1}}}}/{({{2^{{R_{m,k+1}}}}-1})}>{{{2^{{R_{m,k}}}}{\beta_{m,k}}}}/{({{2^{{R_{m,k}}}}-1})}$,
which can be justified by substituting the definition of $\beta_{m,k}$ into
(144). This result is consistent with the assumption of (127). By using (134)
together with $\theta_{m,1}^{*}<\cdots<\theta_{m,K}^{*}$, it follows that

$\frac{{\zeta_{m,k}^{*}}}{{{2^{{R_{m,k}}}}-1}}<\frac{{\zeta_{m,k+1}^{*}}}{{{2^{{R_{m,k+1}}}}-1}}-\zeta_{m,k}^{*}<\frac{{\zeta_{m,k+1}^{*}}}{{{2^{{R_{m,k+1}}}}-1}}.$
(148)

This leads to (41).

##  Appendix G Proof of Remark 2

By using (145) and (146), we have

$\displaystyle\frac{{{{\bf{1}}_{K}}^{\rm{T}}{{\bf{L}}_{m}}^{-1}{{\bf{b}}_{m}}}}{{{{\bf{a}}_{m,k}}^{\rm{T}}{{\bf{L}}_{m}}^{-1}{{\bf{b}}_{m}}}}$
$\displaystyle=\frac{{\sum\nolimits_{i=1}^{K}{{d_{m,i}}}}}{{\frac{{{d_{m,k}}}}{{{2^{{R_{m,k}}}}-1}}-\sum\nolimits_{i=1}^{k-1}{{d_{m,i}}}}}\geq\frac{{\sum\nolimits_{i=1}^{K}{{d_{m,i}}}}}{{\frac{{{d_{m,k}}}}{{{2^{{R_{m,k}}}}-1}}}}$
$\displaystyle\geq\frac{{{d_{m,k}}}}{{\frac{{{d_{m,k}}}}{{{2^{{R_{m,k}}}}-1}}}}={2^{{R_{m,k}}}}-1.$
(149)

Hence, the asymptotic outage probability in (40) is lower bounded by (42).

##  Appendix H Proof of Theorem 2

The associated Lagrangian is given by (150) at the top of this page.

$\bar{\mathcal{L}}\left({{{\left\\{{{R_{m,k}}}\right\\}}_{m\in\left[{1,M}\right]\hfill\atop
k\in\left[{1,K}\right]\hfill}},{{\left\\{{{\theta_{m,k}}}\right\\}}_{m\in\left[{1,M}\right]\hfill\atop
k\in\left[{1,K}\right]\hfill}},{{\left\\{{{\mathchar
22\relax\mkern-10.0mu\lambda_{m,k}}}\right\\}}_{m\in\left[{1,M}\right]\hfill\atop
k\in\left[{1,K}\right]\hfill}},{{\left\\{{{\iota_{m,k}}}\right\\}}_{m\in\left[{1,M}\right]\hfill\atop
k\in\left[{1,K}\right]\hfill}},{{\left\\{{{\varpi_{m,k,i}}}\right\\}}_{m\in\left[{1,M}\right]\hfill\atop{k\in\left[{1,K}\right]\hfill\atop
i\in\left[{k,K}\right]\hfill}}}}\right)\\\
=-\sum\nolimits_{m\in\left[{1,M}\right]\hfill\atop
k\in\left[{1,K}\right]\hfill}{\left({1-{\phi_{m,k}}{\theta_{m,k}}^{-\left({{N_{r}}-M+1}\right)}}\right){R_{m,k}}}-\sum\nolimits_{m\in\left[{1,M}\right]\hfill\atop
k\in\left[{1,K}\right]\hfill}{{\mathchar
22\relax\mkern-10.0mu\lambda_{m,k}}{R_{m,k}}}-\sum\nolimits_{m\in\left[{1,M}\right]\hfill\atop
k\in\left[{1,K}\right]\hfill}{{\iota_{m,k}}{\theta_{m,k}}}\\\
+\sum\nolimits_{m\in\left[{1,M}\right]\hfill\atop{k\in\left[{1,K}\right]\hfill\atop
i\in\left[{k,K}\right]\hfill}}{{\varpi_{m,k,i}}\left({{R_{m,i}}-{{\log}_{2}}\left({1+\frac{{{\zeta_{m,i}}}}{{{\theta_{m,k}}+\sum\nolimits_{l=1}^{i-1}{{\zeta_{m,l}}}}}}\right)}\right)},$
(150)

The necessary KKT conditions are given as

$\mathchar 22\relax\mkern-10.0mu\lambda_{m,k}^{*}R_{m,k}^{*}=0,$ (151)
$\iota_{m,k}^{*}\theta_{m,k}^{*}=0,$ (152)
$\varpi_{m,k,i}^{*}\left({R_{m,i}^{*}-{{\log}_{2}}\left({1+\frac{{{\zeta_{m,i}}}}{{\theta_{m,k}^{*}+\sum\limits_{l=1}^{i-1}{{\zeta_{m,l}}}}}}\right)}\right)=0,$
(153)

(154) and (155) at the top of the next page.

${\left.{\frac{{\partial\bar{\mathcal{L}}}}{{\partial{R_{m,k}}}}}\right|_{{{\left\\{{R_{m,k}^{*}}\right\\}}_{m\in\left[{1,M}\right]\hfill\atop
k\in\left[{1,K}\right]\hfill}},{{\left\\{{\theta_{m,k}^{*}}\right\\}}_{m\in\left[{1,M}\right]\hfill\atop
k\in\left[{1,K}\right]\hfill}},{{\left\\{{\mathchar
22\relax\mkern-10.0mu\lambda_{m,k}^{*}}\right\\}}_{m\in\left[{1,M}\right]\hfill\atop
k\in\left[{1,K}\right]\hfill}},{{\left\\{{\iota_{m,k}^{*}}\right\\}}_{m\in\left[{1,M}\right]\hfill\atop
k\in\left[{1,K}\right]\hfill}},{{\left\\{{\varpi_{m,k,i}^{*}}\right\\}}_{m\in\left[{1,M}\right]\hfill\atop{k\in\left[{1,K}\right]\hfill\atop
i\in\left[{k,K}\right]\hfill}}}}}=0,$ (154)
${\left.{\frac{{\partial\bar{\mathcal{L}}}}{{\partial{\theta_{m,k}}}}}\right|_{{{\left\\{{R_{m,k}^{*}}\right\\}}_{m\in\left[{1,M}\right]\hfill\atop
k\in\left[{1,K}\right]\hfill}},{{\left\\{{\theta_{m,k}^{*}}\right\\}}_{m\in\left[{1,M}\right]\hfill\atop
k\in\left[{1,K}\right]\hfill}},{{\left\\{{\mathchar
22\relax\mkern-10.0mu\lambda_{m,k}^{*}}\right\\}}_{m\in\left[{1,M}\right]\hfill\atop
k\in\left[{1,K}\right]\hfill}},{{\left\\{{\iota_{m,k}^{*}}\right\\}}_{m\in\left[{1,M}\right]\hfill\atop
k\in\left[{1,K}\right]\hfill}},{{\left\\{{\varpi_{m,k,i}^{*}}\right\\}}_{m\in\left[{1,M}\right]\hfill\atop{k\in\left[{1,K}\right]\hfill\atop
i\in\left[{k,K}\right]\hfill}}}}}=0,$ (155)

Substituting (150) into (154) and (155) respectively leads to

$-\left({1-\phi_{m,k}\theta{{{}_{m,k}^{*}}^{-\left({{N_{r}}-M+1}\right)}}}\right)-\mathchar
22\relax\mkern-10.0mu\lambda_{m,k}^{*}\\\
+\sum\nolimits_{t\in\left[{1,k}\right]}{\varpi_{m,t,k}^{*}}=0,$ (156)
$\displaystyle-\left({{N_{r}}-M+1}\right)\phi_{m,k}{\theta_{m,k}^{*}}^{-\left({{N_{r}}-M+2}\right)}R_{m,k}^{*}-\iota_{m,k}^{*}+$
$\displaystyle\sum\limits_{i\in\left[{k,K}\right]}{\frac{{{\varpi_{m,k,i}^{*}\zeta_{m,i}}}}{{\left({\theta_{m,k}^{*}+\sum\limits_{l=1}^{i-1}{{\zeta_{m,l}}}}\right)\left({\theta_{m,k}^{*}+\sum\limits_{l=1}^{i}{{\zeta_{m,l}}}}\right)\ln
2}}}$ $\displaystyle=0.$ (157)

Note that $R_{m,k}^{*},\theta_{m,k}^{*}>0$, it follows from (151) and (152)
that $\mathchar 22\relax\mkern-10.0mu\lambda_{m,k}^{*}$ and
$\iota_{m,k}^{*}=0$. Thus, (156) and (H) can be rewritten as

${1-\phi_{m,k}\theta{{{}_{m,k}^{*}}^{-\left({{N_{r}}-M+1}\right)}}}=\sum\limits_{t\in\left[{1,k}\right]}{\varpi_{m,t,k}^{*}},$
(158)
$\displaystyle\left({{N_{r}}-M+1}\right)\phi_{m,k}{\theta_{m,k}^{*}}^{-\left({{N_{r}}-M+2}\right)}R_{m,k}^{*}=$
$\displaystyle\sum\limits_{i\in\left[{k,K}\right]}{\frac{{{\varpi_{m,k,i}^{*}\zeta_{m,i}}}}{{\left({\theta_{m,k}^{*}+\sum\limits_{l=1}^{i-1}{{\zeta_{m,l}}}}\right)\left({\theta_{m,k}^{*}+\sum\limits_{l=1}^{i}{{\zeta_{m,l}}}}\right)\ln
2}}}.$ (159)

For further simplification, we define
${{\bar{i}}_{m,k}}=\mathop{\arg}\limits_{i}\min\left\\{{{{\log}_{2}}\left({1+{{{\zeta_{m,i}}}}/{({\theta_{m,k}^{*}+\sum\nolimits_{l=1}^{i-1}{{\zeta_{m,l}}}})}}\right),i\in\left[{k,K}\right]}\right\\}$
and we assume that ${{\bar{i}}_{m,k}}$ uniquely exists. Hence, we have

${R_{m,i}^{*}<{{\log}_{2}}\left({1+\frac{{{\zeta_{m,i}}}}{{\theta_{m,k}^{*}+\sum\limits_{l=1}^{i-1}{{\zeta_{m,l}}}}}}\right)},\,i\in[k,K]\wedge
i\neq{{\bar{i}}_{m,k}}.$ (160)

According to (153), we deduce that
$\varpi_{m,k,i\neq{{\bar{i}}_{m,k}}}^{*}=0$. Moreover, by setting $k=1$ in
(158), we have

${1-\phi_{m,1}\theta{{{}_{m,1}^{*}}^{-\left({{N_{r}}-M+1}\right)}}}={\varpi_{m,1,1}^{*}},$
(161)

It is well known that the outage probability without perfect CSI at the BS is
nonzero unless no information is conveyed. Hence, it is reasonable to assume
that
$p_{m,1}^{out\\_asy}={1-\phi_{m,1}\theta{{{}_{m,1}^{*}}^{-\left({{N_{r}}-M+1}\right)}}}>0$.
Accordingly, it follows from (161) that ${\varpi_{m,1,1}^{*}}>0$,
${{\bar{i}}_{m,1}}=1$ and $\varpi_{m,1,i\neq 1}^{*}=0$. In an analogous way,
it is readily proved that ${{\bar{i}}_{m,k}}=k$ by successively setting
$k=2,3,\cdots$ in (161), i.e., $\varpi_{m,k,i\neq k}^{*}=0$ and
$\varpi_{m,k,k}^{*}>0$. With the results, (158) and (H) respectively
degenerates to

${1-\phi_{m,k}\theta{{{}_{m,k}^{*}}^{-\left({{N_{r}}-M+1}\right)}}}={\varpi_{m,k,k}^{*}},$
(162)
$\displaystyle\left({{N_{r}}-M+1}\right)\phi_{m,k}{\theta_{m,k}^{*}}^{-\left({{N_{r}}-M+2}\right)}R_{m,k}^{*}=$
$\displaystyle{\frac{{{\varpi_{m,k,k}^{*}\zeta_{m,k}}}}{{\left({\theta_{m,k}^{*}+\sum\limits_{l=1}^{k-1}{{\zeta_{m,l}}}}\right)\left({\theta_{m,k}^{*}+\sum\limits_{l=1}^{k}{{\zeta_{m,l}}}}\right)\ln
2}}}.$ (163)

Moreover, combining (153) with $\varpi_{m,k,k}^{*}>0$ yields (45). Then
putting (45) and (162) into (H) together with some algebraic manipulations
produces

$\varrho(\theta{{}_{m,k}^{*}})=\ell(\theta{{}_{m,k}^{*}}),$ (164)

where $\varrho(x)$ and $\ell(x)$ are defined in (46) and (47), respectively.

Clearly from (164), the intersection of two curves ${\varrho(x)}$ and
${\ell(x)}$ gives the optimal value $\theta_{m,k}^{*}$. The curves of
${\varrho(x)}$ and ${\ell(x)}$ are plotted in Fig. 10. Clearly, ${\varrho(x)}$
and ${\ell(x)}$ are both increasing functions. Additionally, it is not hard to
prove that ${\varrho(x)}$ increases much faster than ${\ell(x)}$ as
$x\to\infty$. Hence, the zero point of (164), i.e., $\theta_{m,k}^{*}$,
exists, and $\theta_{m,k}^{*}$ can be numerically solved by using bisection
method. After determining $\theta_{m,k}^{*}$, $R_{m,k}^{*}$ can be calculated
by using (45).

$x$$y$$\ell(x)$$\varrho(x)$$\theta_{m,k}^{*}$$\phi_{m,k}^{\frac{1}{NR-M+1}}$$\sum\nolimits_{l=1}^{k-1}{{\zeta_{m,l}}}\sum\nolimits_{l=1}^{k}{{\zeta_{m,l}}}$
Figure 10: The increasing monotonicity of ${\varrho(x)}$ and ${\ell(x)}$.

By combining (45) and (160), the same conclusion can be drawn as Remark 1,
that is, the optimal transmission rates satisfy the inequalities in (48). The
proof of Theorem 2 is thus accomplished.

## References

  * [1] S. Han, X. Xu, S. Fang, Y. Sun, Y. Cao, X. Tao, and P. Zhang, “Energy efficient secure computation offloading in NOMA-based mMTC networks for IoT,” _IEEE Internet Things J._ , vol. 6, no. 3, pp. 5674–5690, Jun. 2019\.
  * [2] W. Hao, Z. Chu, F. Zhou, S. Yang, G. Sun, and K.-K. Wong, “Green communication for NOMA-based CRAN,” _IEEE Internet Things J._ , vol. 6, no. 1, pp. 666–678, Jul. 2018.
  * [3] A. Ericsson, “Ericsson mobility report,” _Recuperado de: https://www. ericsson. com/assets/local/mobility-report/documents/2018/ericsson-mobility-report-june-2018. pdf (Accedido el 09/09/2018)_ , pp. 3–13, Jun. 2019.
  * [4] V. Cisco, “Cisco visual networking index: Forecast and trends, 2017–2022,” _White Paper_ , 2019.
  * [5] R. Duan, J. Wang, C. Jiang, H. Yao, Y. Ren, and Y. Qian, “Resource allocation for multi-UAV aided IoT NOMA uplink transmission systems,” _IEEE Internet Things J._ , vol. 6, no. 4, pp. 7025–7037, Aug. 2019.
  * [6] M. Alkhawatrah, Y. Gong, G. Chen, S. Lambotharan, and J. A. Chambers, “Buffer-aided relay selection for cooperative NOMA in the internet of things,” _IEEE Internet Things J._ , vol. 13, no. 3, pp. 524–537, Jun. 2019\.
  * [7] Z. Ding, Y. Liu, J. Choi, Q. Sun, M. Elkashlan, C. L. I, and H. V. Poor, “Application of non-orthogonal multiple access in LTE and 5G networks,” _IEEE Commun. Mag._ , vol. 55, no. 2, pp. 185–191, Feb. 2017.
  * [8] W. A. Al-Hussaibi and F. H. Ali, “Efficient user clustering, receive antenna selection, and power allocation algorithms for massive MIMO-NOMA systems,” _IEEE Access_ , vol. 7, pp. 31 865–31 882, Feb. 2019.
  * [9] Y. Liu, Z. Ding, M. Elkashlan, and H. V. Poor, “Cooperative non-orthogonal multiple access with simultaneous wireless information and power transfer,” _IEEE J. Sel. Areas Commun._ , vol. 34, no. 4, pp. 938–953, Apr. 2016.
  * [10] W. Lv, R. Wang, J. Wu, J. Xu, P. Li, and J. Dou, “Degrees of freedom of the circular multirelay MIMO interference channel in IoT networks,” _IEEE Internet Things J._ , vol. 5, no. 3, pp. 1957–1966, Jun. 2018.
  * [11] J. Wu, B. Tan, J. Wu, and R. Wang, “Efficient soft video MIMO design to combine diversity and spatial multiplexing gain,” _IEEE Internet Things J._ , vol. 6, no. 3, pp. 5461–5472, Jun. 2019.
  * [12] L. Dai, B. Wang, Y. Yuan, S. Han, I. Chih-Lin, and Z. Wang, “Non-orthogonal multiple access for 5G: solutions, challenges, opportunities, and future research trends,” _IEEE Commun. Mag._ , vol. 53, no. 9, pp. 74–81, Sep. 2015\.
  * [13] Q. Sun, S. Han, I. Chin-Lin, and Z. Pan, “On the ergodic capacity of MIMO NOMA systems,” _IEEE Wireless Commun. Lett._ , vol. 4, no. 4, pp. 405–408, Aug. 2015.
  * [14] Y. Liu, G. Pan, H. Zhang, and M. Song, “On the capacity comparison between MIMO-NOMA and MIMO-OMA,” _IEEE Access_ , vol. 4, pp. 2123–2129, May 2016\.
  * [15] 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.
  * [16] H. Wang, R. Zhang, R. Song, and S.-H. Leung, “A novel power minimization precoding scheme for MIMO-NOMA uplink systems,” _IEEE Commun. Lett._ , vol. 22, no. 5, pp. 1106–1109, May 2018.
  * [17] 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.
  * [18] Z. Ding, L. Dai, and H. V. Poor, “MIMO-NOMA design for small packet transmission in the internet of things,” _IEEE access_ , vol. 4, pp. 1393–1405, Mar. 2016.
  * [19] Y. Al-Eryani, E. Hossain, and D. I. Kim, “Generalized coordinated multipoint (GCoMP)-enabled NOMA: Outage, capacity, and power allocation,” _IEEE Trans. Commun._ , vol. PP, no. 99, pp. 1–1, Jul. 2019\.
  * [20] Z. Ding, Z. Yang, P. Fan, and H. V. Poor, “On the performance of non-orthogonal multiple access in 5G systems with randomly deployed users,” _IEEE Signal Process. Lett._ , vol. 21, no. 12, pp. 1501–1505, Dec. 2014.
  * [21] Z. Yang, Z. Ding, P. Fan, and G. K. Karagiannidis, “On the performance of non-orthogonal multiple access systems with partial channel information,” _IEEE Trans. Commun._ , vol. 64, no. 2, pp. 654–667, Feb. 2016.
  * [22] Z. Ding and H. V. Poor, “Design of massive-MIMO-NOMA with limited feedback,” _IEEE Signal Process. Lett._ , vol. 23, no. 5, pp. 629–633, May 2016.
  * [23] 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.
  * [24] M.-Y. Gong, Z. Yang, and B. Lyu, “Antenna diversity for downlink MIMO-NOMA systems with partial channel state information,” _IEEE Commun. Lett._ , vol. 22, no. 10, pp. 2172–2175, Oct. 2018.
  * [25] M.-y. Gong and Z. Yang, “The application of antenna diversity to NOMA with statistical channel state information,” _IEEE Trans. Veh. Technol._ , vol. 68, no. 4, pp. 3755–3765, Apr. 2019.
  * [26] D. Tong, Y. Ding, Y. Liu, and Y. Wang, “A MIMO-NOMA framework with complex-valued power coefficients,” _IEEE Trans. Veh. Technol._ , vol. 68, no. 3, pp. 2244–2259, Mar. 2019.
  * [27] 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.
  * [28] Z. Shi, G. Yang, Y. Fu, H. Wang, and S. Ma, “Performance analysis of MIMO-NOMA systems with randomly deployed users,” in _Proc. IEEE Global Commun. Conf. (GLOBECOM’18)_ , Abu Dhabi, United Arab Emirates, Dec. 2018, pp. 1–7.
  * [29] M. Dohler, E. Lefranc, and H. Aghvami, “Space-time block codes for virtual antenna arrays,” in _Proc. IEEE Int. Symp. Personal Indoor and Mobile Radio Commun. (PIMRC’02)_ , Papeete, Tahiti, French, Mar. 2002, pp. 414–417.
  * [30] X. Li, X. Yang, L. Li, J. Jin, N. Zhao, and C. Zhang, “Performance analysis of distributed MIMO with ZF receivers over semi-correlated $\mathcal{K}$ fading channels,” _IEEE Access_ , vol. 5, pp. 9291–9303, Jun. 2017.
  * [31] H.-W. Chang and L.-C. Wang, “A low-complexity uplink multiuser scheduling for virtual MIMO systems,” _IEEE Trans. Veh. Technol._ , vol. 65, no. 1, pp. 463–466, Jan. 2016.
  * [32] A. Alexiou, “Wireless world 2020: Radio interface challenges and technology enablers,” _IEEE Veh. Technol. Mag._ , vol. 9, no. 1, pp. 46–53, Feb. 2014\.
  * [33] S. Cui, A. J. Goldsmith, and A. Bahai, “Energy-efficiency of MIMO and cooperative MIMO techniques in sensor networks,” _IEEE J. Sel. Areas Commun._ , vol. 22, no. 6, pp. 1089–1098, Aug. 2004.
  * [34] S. K. Jayaweera, “Virtual MIMO-based cooperative communication for energy-constrained wireless sensor networks,” _IEEE Trans. Wireless Commun._ , vol. 5, no. 5, pp. 984–989, Jun. 2006.
  * [35] Z. Rafique, B.-C. Seet, and A. Al-Anbuky, “Performance analysis of cooperative virtual MIMO systems for wireless sensor networks,” _Sensors_ , vol. 13, no. 6, pp. 7033–7052, May 2013.
  * [36] I. Dey, M. M. Butt, and N. Marchetti, “Throughput analysis for virtual MIMO WSNs over measured MIMO channels,” _IEEE Trans. Instrum. Meas._ , vol. 68, no. 1, pp. 297–299, Jan. 2019.
  * [37] W. Feng, J. Wang, Y. Chen, X. Wang, N. Ge, and J. Lu, “UAV-aided MIMO communications for 5G internet of things,” _IEEE Internet Things J._ , vol. 6, no. 2, pp. 1731–1740, Oct. 2018.
  * [38] S. Kisseleff, I. F. Akyildiz, and W. Gerstacker, “Distributed beamforming for magnetic induction based body area sensor networks,” in _Proc. IEEE Global Commun. Conf. (GLOBECOM’16)_ , Washington, DC, USA, Dec. 2016, pp. 1–7.