# Exploiting Underlay Spectrum Sharing in Cell-Free Massive MIMO Systems

Diluka Loku Galappaththige, and Gayan Amarasuriya Aruma Baduge Authors are
with the School of Electrical, Computer, and Biomedical Engineering, Southern
Illinois University, Carbondale, IL, USA 62901, Email: {diluka.lg,
<EMAIL_ADDRESS>This work in part has been presented at IEEE
International Conference on Communications (ICC), May 2019, Shanghai, China,
[1] and also in IEEE ICC May 2020, Dublin, Ireland [2].

###### Abstract

We investigate the coexistence of underlay spectrum sharing in cell-free
massive multiple-input multiple-output (MIMO) systems. A primary system with
geographically distributed primary access points (P-APs) serves a multitude of
primary users (PUs), while a secondary system serves a large number of
secondary users (SUs) in the same primary/licensed spectrum by exploiting the
underlay spectrum sharing. To mitigate the secondary co-channel interference
inflected at PUs, stringent secondary transmit power constraints are defined
for the secondary access points (S-APs). A generalized pilots sharing scheme
is used to locally estimate the uplink channels at P-APs/S-APs, and thereby,
conjugate precoders are adopted to serve PUs/SUs in the same time-frequency
resource element. Moreover, the effect of a user-centric AP clustering scheme
is investigated by assigning a suitable set of APs to a particular user. The
impact of estimated downlink (DL) channels at PUs/SUs via DL pilots beamformed
by P-APs/S-APs is investigated. The achievable primary/secondary rates at
PUs/SUs are derived for the statistical DL and estimated DL CSI cases. User-
fairness for PUs/SUs is achieved by designing efficient transmit power control
policies based on a multi-objective optimization problem formulation of joint
underlay spectrum sharing and max-min criteria. The proposed orthogonal
multiple-access based analytical framework is also extended to facilitate non-
orthogonal multiple-access. Our analysis and numerical results manifest that
the primary/secondary performance of underlay spectrum sharing can be boosted
by virtue of the average reduction of transmit powers/path-losses, uniform
coverage/service, and macro-diversity gains, which are inherent to distributed
transmissions/receptions of cell-free massive MIMO.

##  I Introduction

Massive multiple-input multiple-output (MIMO) operating in sub-6 GHz can
simultaneously serve many users in the same time-frequency resource element by
virtue of aggressive spatial multiplexing gains rendered by large base-station
(BS) antenna arrays [1, 2, 3, 4]. Co-located massive MIMO in which all BS
antennas are packed into the same array is currently being deployed in the
United States [5]. Thus, the co-located massive MIMO enabled with fully-
digital beamforming has already become a reality [6].

Recently, a distributed/cell-free massive MIMO architecture, which deploys a
large number of distributed access points (APs), is proposed to enhance the
coverage probability and hence to provide a uniformly better service to users
within a much larger geographical area [7, 8, 9, 10, 11, 12, 13, 14].
Coordinated multi-point (CoMP) and network MIMO [15] are two related
technologies, which exploit the notion of cooperation among BSs located within
dedicated cells. On the contrary to these existing techniques, all APs are
deployed in a cell-free architecture aiming to jointly serve all users for a
given/large geographical area with no cell-boundaries. The uplink (UL) channel
state information (CSI) is locally estimated at each AP via pilots send by
users, and these APs are connected to a central processing unit (CPU) via a
fronthaul/backhaul network [8]. By exploiting the channel reciprocity of time-
division duplexing (TDD) mode of operation, APs acquire the downlink (DL) CSI
via UL channel estimates and thereby design precoders for DL data
transmission. The requirement for CSI exchange among APs via the CPU strictly
depends on the AP precoder design. Owing to the fact that the distributed APs
enable a user-centric architecture, the average transmission distances of a
cell-free massive MIMO are inherently smaller than that of the co-located
counterpart [16]. This benefit translates into transmit power and path-loss
reductions, which in turn lead to boosted energy efficiencies. The distributed
APs circumvent the impediments caused by spatially correlated fading and
shadowing due to large obstacles. The underlying macro-diversity gains can be
exploited to boost the achievable rates [17, 18]. Joint beamfoming at a
massive number of APs also enables unprecedented spatial multiplexing gains.
Thus, cell-free massive MIMO enables a large number of concurrent connections
with a guaranteed uniform service throughout a given/large geographical area.

Cognitive radio techniques based on the spectrum sharing are extensively
explored to mitigate spectrum scarcity and underutilization/holes for the
next-generation wireless communication systems [19, 20, 21]. In particular,
many spectrum sharing techniques are evolved through three main paradigms,
namely underlay, overlay, and interweave [19]. The overlay spectrum sharing
involves sophisticated signal processing techniques and requires codebook
knowledge of the non-cognitive users. The interweave spectrum sharing adopts
opportunistic frequency reuse over the available spectrum holes and hence
requires stringent activity information of the non-cognitive users. The
underlay spectrum sharing allows cognitive/secondary users (SUs) to
simultaneously operate within the licensed spectrum if the secondary
interference inflicted at the non-cognitive/primary users (PUs) is below a
certain threshold. The main reasons for adopting the underlay spectrum sharing
in this paper are its implementation simplicity and the achievable high
spectrum utilization [19, 20, 21], compared to the significantly sophisticated
overlay and interweave counterparts. Thus, in the underlay spectrum sharing,
an underlaid secondary system can be simultaneously operated in the same
licensed primary spectrum by defining stringent secondary transmit power
constraints such that the secondary co-channel interference (CCI) caused to
the PUs always falls below a predefined primary interference threshold (PIT)
[19].

###  I-A Related prior research on underlay spectrum sharing with massive MIMO

In [21], an initial foundation for investigating the feasibility of underlay
spectrum sharing in massive MIMO with co-located antenna arrays at BSs is
established. The BSs in [21] use the maximum ratio transmission (MRT) for
signal transmission, and the DL achievable rate of the secondary system which
is underlaid in a primary massive MIMO system is investigated. Reference [20]
investigates the impact of inherent intra/inter-cell pilot contamination in
multi-cell multi-user massive MIMO system with underlay spectrum sharing. In
[22], a dual-hop enabled spectrum sharing system is analyzed, and the
achievable rates are derived. Thereby, the detrimental effects of inter/intra-
cell pilot contamination are investigated. The fundamental performance limits
for relay selection strategies in massive MIMO two-way relaying are explored
for perfect CSI in [23]. Moreover, in [23], the asymptotic achievable rates
are derived for the best relay selection by deriving the asymptotic signal-to-
interference-plus-noise ratio (SINR). In [24], the achievable rates of
reserve-TDD based underlay spectrum sharing are presented. In [13], the
achievable rates of underlay spectrum sharing co-located massive MIMO non-
orthogonal multiple-access (NOMA) are presented. In [25], a low-complexity
sub-optimal user-clustering technique for NOMA based underlay spectrum sharing
in cell-free massive MIMO is proposed, and thereby, the achievable rates are
derived for fixed transmit power allocation. In [26, 27], the performance
bounds of spectrum sharing for massive MIMO with stochastic BS/user locations
are derived. Moreover, by investigating pilot contamination, path-loss
inversion power control, and receiver association policies, the secondary
interference for a random cognitive massive MIMO system is characterized in
[26, 27]. In [28], the quality-of-service aware power allocation and user
selection schemes are studied for cognitive massive MIMO systems. Pilot
decontamination techniques are proposed to asymptotically mitigate the
residual interference in an underlaid single user massive MIMO cognitive radio
system in [29]. In [17], the macro-multiplexing gain achieved from
optimization of antenna locations is characterized in terms of the ambient
dimension of the cell and the path-loss exponent. Reference [18] derives the
upper and lower bounds of the achievable rate with the perfect/imperfect CSI
for cell-free massive MIMO systems. Thereby, [18] shows that the bounds of the
achievable rate converge to a common lower bound owing to the extra distance-
diversity or macro-diversity gain offered by distributed antennas in cell-free
massive MIMO. In [30], the effects of finite capacity of fronthauls in the
presence of residual hardware impairments at the users and APs are
investigated by deriving the achievable rates for the compress-forward-
estimate, estimate-compress-forward, and estimate-multiply-compress-forward
strategies. Reference [31] proposes a low-complexity power allocation scheme
to maximize the energy efficiency for a cell-free massive MIMO system with
user centric approach operating at millimeter-wave (mmWave) frequencies. In
[32], two pilot assignments, namely the user-group and interference-based
$K$-means schemes are proposed for the structured massive access.

###  I-B Our motivation

The aforementioned related prior references [20, 22, 24, 23, 26, 29, 21] have
investigated the coexistence of massive MIMO and underlay spectrum sharing
with only co-located antenna arrays at the BSs. The closely related references
[13] and [25], respectively, investigate the rate performance of co-located
and cell-free massive MIMO NOMA with underlay spectrum sharing by only
considering fixed power allocation. To the best of our knowledge, multi-
objective max-min fairness-based transmit power allocation, impact of
beamformed DL pilots by APs, and the achievable rates with DL estimated CSI at
users have not yet been investigated. Thus, in our paper, we fill this gap by
exploring multi-objective power allocation, impact of DL pilots, and
practically realizable performance bounds of underlay spectrum sharing in
cell-free massive MIMO with imperfectly estimated UL/DL CSI. Since the APs are
spatially distributed over a given geographical area, the distributed
transmissions of cell-free massive MIMO architecture can be beneficial in
providing uniformly better average rate to users by virtue of mitigating the
near-far effect via max-min power control than that of the co-located massive
MIMO counterparts. The average amount of secondary CCI inflicted at a
particular PU from cell-free/distributed massive MIMO transmissions can be
better constrained as specified by the PIT with low-power distributed APs.
Moreover, cell-free massive MIMO has more robustness against the detrimental
effects of correlated small/large scale fading than the co-located counterpart
[8]. Both primary and secondary systems can provide a higher coverage
probability because there are no cell boundaries, and users are much closer to
the APs in cell-free setting. Thus, the cell-free massive MIMO can
significantly boost the performance of underlay spectrum sharing.

###  I-C Our contribution and its difference relative to the existing
literature

Our main contribution is to investigate different UL/DL CSI cases at
APs/users, and transmit power control and their effects on the achievable
rates of underlay spectrum sharing in orthogonal multiple-access (OMA)/NOMA-
based cell-free massive MIMO. Specifically, we derive the performance metrics
by proposing max-min based multi-objective transmit power control algorithms
and by exploring the deleterious impact of imperfectly estimated UL CSI at
APs, availability of estimated DL CSI at the users, and effects of using
statistical DL CSI with imperfect successive interference cancellation (SIC)
signal decoding. Both the primary system and the secondary system, which is
underlaid within the primary licensed spectrum, adopt a generalized pilot
sharing scheme to minimize the training overhead. Thus, the pilot sequences
sent by PUs and SUs are used to locally estimate the UL channels at primary
access points (P-APs)/secondary access points (S-APs), respectively. Thereby,
the impact of imperfectly estimated CSI is considered for our analysis. Our
performance metrics for cell-free underlay spectrum sharing are categorized by
taking into account the availability of long-term/statistical CSI and
estimated DL CSI at PUs/SUs. When only UL channel estimation is used, the
users rely on statistical CSI for signal decoding. Nonetheless, when APs
beamform DL pilots, the users adopt estimated DL CSI to decode signals. Thus,
the performance bounds are established for these two user CSI cases. Moreover,
the effect of user-centric achieved through clustering the APs that serve for
a particular user in both primary an secondary systems is investigated for
cell-free underlay spectrum sharing.

To mitigate the secondary CCI inflected at PUs due to simultaneous
transmission, the stringent secondary transmit power constraints are
introduced at S-APs. Thus, by defining a PIT at PUs, the secondary transmit
power is constrained such that the total secondary CCI at any PU falls bellow
the predefined PIT. Then, in the presence of imperfectly estimated UL CSI at
APs and DL CSI at users with intra-system pilot contamination and estimation
errors, the achievable rates for PUs/SUs are derived. A multi-objective
transmit power control algorithm is designed based on the max-min fairness
criterion to guarantee a uniform quality-of-service among all the users.
Moreover, the above OMA-based system model is extended to facilitate NOMA
transmissions, and the corresponding performance bounds are established. The
practical viability of underlay spectrum sharing in cell-free massive MIMO is
explored by using numerical results through our analysis and Monte-Carlo
simulations.

This paper goes well beyond our related conference papers [1, 2] by presenting
a multi-objective optimization of max-min fairness-based transmit power
allocation for the OMA-based primary system, user-centric AP clustering
aspects, beamforming of DL pilots, impact of estimated DL channel estimates,
adoption of estimated DL CSI for signal decoding, and the corresponding
achievable rates with estimated DL CSI for underlay spectrum sharing within
cell-free massive MIMO. All numerical results/comparisons, except for Fig. 8,
and their descriptions are distinctive from those of [1, 2], and the
corresponding figures have been regenerated for different system parameters
with respect to [1, 2].

Notation: $\mathbf{z}^{T}$ denotes the transpose of $\mathbf{z}$ . The
conjugate of $z$ is denoted by $z^{*}$. The notation
$z\sim\mathcal{CN}(\cdot,\cdot)$ denotes that $z$ is a complex-valued
circularly symmetric Gaussian distributed random variable. The operators
$\mathbb{E}\\!\left[{\cdot}\right]$ and
$\mathbb{V}\mathrm{ar}\\!\left[{\cdot}\right]$ are the expectation and
variance, respectively.

##  II System, channel and signal models

###  II-A System and channel model for OMA/NOMA

$\mathrm{P\mbox{-}\\!AP_{m}}$$\mathrm{S\mbox{-}\\!AP_{n}}$$U_{P}(k)$$U_{S}(l)$$v_{ml}$$f_{mk}$$g_{nl}$$u_{nk}$$\mathrm{P\mbox{-}CPU}$$\mathrm{S\mbox{-}CPU}$$\mathrm{P\mbox{-}\\!AP_{m}}$$\mathrm{S\mbox{-}\\!AP_{n}}$${v_{mbj}}$${f_{mai}}$${u_{nai}}$${g_{nbj}}$$\mathrm{P\mbox{-}CPU}$$\mathrm{S\mbox{-}CPU}$$\mathrm{PC_{a}}$$U_{S}(bL)$$U_{S}(b1)$$U_{S}(bj)$$\mathrm{SC_{b}}$$U_{P}(aK)$$U_{P}(a1)$$U_{P}(ai)$a.
OMA system modelb. NOMA system model Figure 1: A cell-free massive MIMO system
models with underlay spectrum sharing.

We consider a TDD cell-free massive MIMO system with underlay spectrum sharing
(see Fig. 1a). A secondary system is underlaid within a primary licensed
spectrum in order to enhance the overall spectrum efficiency by eliminating
spectrum holes. The primary system having $M$ single-antenna P-APs serves $K$
single-antenna PUs, while the secondary system with $N$ single-antenna S-APs
uses the same time-frequency spectrum to serve $L$ single-antenna SUs. We
introduce the secondary transmit power constraints for all S-APs to ensure
that the performance of the primary system is not hindered by the simultaneous
secondary transmission in the same primary spectrum. Thus, the CCI caused by
the secondary system at PUs falls below a predefined PIT, which defines a
upper limit for the CCI endurance ability of PUs. The PIT constraint mitigates
excessive secondary CCI at PUs. A synchronized operation between the primary
and the secondary systems is assumed such that all P-APs and S-APs
simultaneously serve all PUs and SUs by adopting spatial multiplexing rendered
by cell-free massive MIMO [8]. Moreover, P-APs/S-APs are connected to their
respective primary/secondary central processing units (P-CPU/S-CPU). The $k$th
PU and the $l$th SU are denoted by $U_{P}(k)$ and $U_{S}(l)$, respectively.

In Fig. 1a, $f_{mk}$, $g_{nl}$, $v_{ml}$, and $u_{nk}$ are the channel
coefficients between the $m$th P-AP and $U_{P}(k)$, the $n$th S-AP and
$U_{S}(l)$, the $m$th P-AP and $U_{S}(l)$, and the $n$th S-AP and $U_{P}(k)$,
respectively, where $m\in\\{1,\cdots,M\\}$, $k\in\\{1,\cdots,K\\}$,
$n\in\\{1,\cdots,N\\}$, and $l\in\\{1,\cdots,L\\}$. The above channels can be
modeled in a unified manner as

$\displaystyle{h}_{ab}=\tilde{h}_{ab}{{\zeta}^{1/2}_{{h}_{ab}}},$ (1)

where $h\in\\{f,g,u,v\\}$, $a\in\\{m,n\\}$, and $b\in\\{k,l\\}$. Here,
${\tilde{h}_{ab}\sim\mathcal{CN}(0,1)}$ captures the independent quasi-static
Rayleigh fading and stays fixed during the coherence interval, while
$\zeta_{{h}_{ab}}$ accounts for the large-scale fading, including path-loss
and shadow fading. Since the large-scale coefficients stay fixed for several
coherence intervals, it is assumed that the large-scale coefficients are known
a-prior at both P-APs and S-APs [8]. Thus, the estimation of these large-scale
fading coefficients can be done in once about tens/hundreds of coherence
intervals [33].

Next, we extend our cell-free massive MIMO underlay spectrum sharing
techniques to facilitate NOMA. To this end, we consider a system setup with
$AK$ number of PUs and $BL$ number of SUs, where $A$ and $B$ are the numbers
of primary clusters (PCs) and secondary clusters (SCs), respectively (see Fig.
1b). Based on the spatial directions of the users [24, 34], $K$ and $L$ number
of PUs and SUs are assigned to each PC and SC, respectively. The $k$th PU in
the $a$th PC and the $l$th SU in the $b$th SC are denoted by $U_{P}(a,k)$ and
$U_{S}(b,l)$, respectively. The channel between the $m$th P-AP and
$U_{P}(a,k)$ is denoted by $f_{mak}$, where $m\in\\{1,\cdots,M\\}$,
$a\in\\{1,\cdots,A\\}$, and $k\in\\{1,\cdots,K\\}$. The channel between the
$n$th S-AP and $U_{S}(b,l)$ is represented by $g_{nbl}$, where
$n\in\\{1,\cdots,N\\}$, $b\in\\{1,\cdots,B\\}$, and $l\in\\{1,\cdots,L\\}$.
Further, $v_{mbl}$ and $u_{nak}$ denote the interference channels between the
$m$th P-AP and $U_{S}(b,l)$, and the $n$th S-AP and $U_{P}(a,k)$,
respectively. These channels are modeled similar to (1) as

$\displaystyle{h}_{qrs}=\tilde{h}_{qrs}{{\zeta}^{1/2}_{{h}_{qrs}}},$ (2)

where $h\in\\{f,g,u,v\\}$, $q\in\\{m,n\\}$, $r\in\\{a,b\\}$, and
$s\in\\{k,l\\}$. Moreover, $\tilde{h}_{qrs}\sim\mathcal{CN}(0,1)$ captures the
small-scale fading, while ${{\zeta}_{{h}_{qrs}}}$ captures the large-scale
fading, including path-loss and shadowing.

###  II-B The UL channel estimation for OMA/NOMA

For OMA, the UL channels are estimated locally at P-APs and S-APs by using
respective user pilots [8]. During UL channel estimation period, $\tau_{p}$
symbols out of the coherence interval having $\tau_{c}$ symbols are used to
transmit the UL pilot. Then, by using these pilot sequences, the channels
$f_{mk}$ and $g_{nl}$ are estimated at P-APs and S-APs, respectively. In
practice, the number of orthogonal pilot sequences is limited as it is defined
by the channel coherence interval [33]. To reduce the pilot overhead and to
increase the number of served PUs/SUs, in this work, the pilot sequences are
shared among PUs and SUs by adopting the following pilot sharing strategy. It
is assumed that $Q\\!\leq\\!\min\\!\left(\\!K,L\\!\right)$ number of pilot
sequences having a length of $\tau_{p}$ symbol duration is shared among PUs
and SUs. Then, we define the complete pilot sets used by PUs
$(\mathbf{\Phi}_{P})$ and SUs $(\mathbf{\Phi}_{S})$ as

$\displaystyle\mathbf{\Phi}_{P}=\left[\mathbf{\Phi};\tilde{\mathbf{\Phi}}_{P}\right]\qquad\text{and}\qquad\mathbf{\Phi}_{S}=\left[\mathbf{\Phi};\tilde{\mathbf{\Phi}}_{S}\right],$
(3)

where $\mathbf{\Phi}\in\mathbb{C}^{Q\times\tau_{p}}$ denotes the pilots shared
by $Q$ PUs/SUs having a length of $\tau_{p}$ symbol duration.
$\tilde{\mathbf{\Phi}}_{P}\in\mathbb{C}^{\left(K-Q\right)\times\tau_{p}}$, and
$\tilde{\mathbf{\Phi}}_{S}\in\mathbb{C}^{\left(L-Q\right)\times\tau_{p}}$ are
the pilots assigned for the remaining $\left(K-Q\right)$ PUs or
$\left(L-Q\right)$ SUs, respectively. We can define the orthogonal properties
among these pilots as $\mathbf{\Phi}^{H}\tilde{\mathbf{\Phi}}_{P}=0$,
$\mathbf{\Phi}^{H}\tilde{\mathbf{\Phi}}_{S}=0$, and
$\tilde{\mathbf{\Phi}}_{P}^{H}\tilde{\mathbf{\Phi}}_{S}=0$. We define
$\mathbf{\Phi}_{P}=\left[\boldsymbol{\phi}_{P_{1}}^{T},\cdots,\boldsymbol{\phi}_{P_{k}}^{T},\cdots,\boldsymbol{\phi}_{P_{K}}^{T}\right]^{T}$
and
$\mathbf{\Phi}_{S}=\left[\boldsymbol{\phi}_{S_{1}}^{T},\cdots,\boldsymbol{\phi}_{S_{l}}^{T},\cdots,\boldsymbol{\phi}_{S_{L}}^{T}\right]^{T}$,
where $\boldsymbol{\phi}_{P_{k}}\in\mathbb{C}^{1\times\tau_{p}}$ and
$\boldsymbol{\phi}_{S_{l}}\in\mathbb{C}^{1\times\tau_{p}}$ are the pilot
sequences sent by $U_{P}(k)$ and $U_{S}(l)$, respectively, and
$\|\boldsymbol{\phi}_{P_{k}}\|^{2}=1$ and
$\|\boldsymbol{\phi}_{S_{l}}\|^{2}=1$ for $k\in\\{1,\cdots,K\\}$ and
$l\in\\{1,\cdots,L\\}$. Then, we can write the pilot signal received at the
$m$th P-AP and the $n$th S-AP as

$\displaystyle\mathbf{y}^{\prime}_{P_{m}}$ $\displaystyle=$
$\displaystyle\sqrt{P_{p}}\sum\nolimits_{k=1}^{K}f_{mk}\boldsymbol{\phi}_{P_{k}}+\sqrt{P_{p}}\sum\nolimits_{l=1}^{L}v_{ml}\boldsymbol{\phi}_{S_{l}}+\mathbf{n}_{P_{m}}^{\prime},$
(4a) $\displaystyle\mathbf{y}^{\prime}_{S_{n}}$ $\displaystyle=$
$\displaystyle\sqrt{P_{p}}\sum\nolimits_{l=1}^{L}g_{nl}\boldsymbol{\phi}_{S_{l}}+\sqrt{P_{p}}\sum\nolimits_{k=1}^{K}u_{nk}\boldsymbol{\phi}_{P_{k}}+\mathbf{n}_{S_{n}}^{\prime},$
(4b)

where $P_{p}=\tau_{p}P$ and $P$ denotes the average transmitted pilot power at
each PU/SU. Moreover, $\mathbf{n}_{P_{m}}^{\prime}$ and
$\mathbf{n}_{S_{n}}^{\prime}$ are additive white Gaussian noise (AWGN)
vectors, having independent and identically distributed (i.i.d.)
$\mathcal{CN}(0,1)$ elements at the $m$th P-AP and the $n$th S-AP,
respectively. The sufficient statistics for estimating $f_{mk}$ and $g_{nl}$
can be obtained by projecting $\boldsymbol{\phi}_{P_{k}}^{H}$ and
$\boldsymbol{\phi}_{S_{l}}^{H}$ onto (4a) and (4b), respectively, as

$\displaystyle y_{P_{mk}}$ $\displaystyle=$
$\displaystyle\boldsymbol{\phi}_{P_{k}}^{H}\mathbf{y}^{\prime}_{P_{m}}=\sqrt{P_{p}}f_{mk}+\sqrt{P_{p}}v_{mk}+n_{P_{m}},$
(5a) $\displaystyle y_{S_{nl}}$ $\displaystyle=$
$\displaystyle\boldsymbol{\phi}_{S_{l}}^{H}\mathbf{y}^{\prime}_{S_{n}}=\sqrt{P_{p}}g_{nl}+\sqrt{P_{p}}u_{nl}+n_{S_{n}},$
(5b)

where $\\{\phi_{P_{k}},\phi_{S_{l}}\\}\\!\in\\!{\bf\Phi}$,
$n_{P_{m}}\\!\\!=\\!\boldsymbol{\phi}_{P_{k}}^{H}\mathbf{n}_{P_{m}}^{\prime}\\!\sim\\!\mathcal{CN}(0,1)$,
and
$n_{S_{n}}\\!\\!=\\!\boldsymbol{\phi}_{S_{l}}^{H}\mathbf{n}_{S_{n}}^{\prime}\\!\sim\\!\mathcal{CN}(0,1)$
as $\boldsymbol{\phi}_{P_{k}}$ and $\boldsymbol{\phi}_{S_{l}}$ are unitary
vectors.

Proposition 1: The minimum mean square error (MMSE) estimates of $f_{mk}$ and
$g_{nl}$ are given by

$\displaystyle\hat{f}_{mk}$ $\displaystyle=$
$\displaystyle\left({\mathbb{E}\\!\left[{y_{P_{mk}}^{*}f_{mk}}\right]}\big{/}{\mathbb{E}\\!\left[{|y_{P_{mk}}|^{2}}\right]}\right)y_{P_{mk}}=c_{P_{mk}}y_{P_{mk}},$
(6a) $\displaystyle\hat{g}_{nl}$ $\displaystyle=$
$\displaystyle\left({\mathbb{E}\\!\left[{y_{S_{nl}}^{*}g_{nl}}\right]}\big{/}{\mathbb{E}\\!\left[{|y_{S_{nl}}|^{2}}\right]}\right)y_{S_{nl}}=c_{S_{nl}}y_{S_{nl}},$
(6b)

where $c_{P_{mk}}$ and $c_{S_{nl}}$ are given by

$\displaystyle c_{P_{mk}}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{\tau_{p}P}{\zeta}_{{f}_{mk}}}{\tau_{p}P\left({\zeta}_{{f}_{mk}}+{\zeta}_{{v}_{mk}}\right)+1}\qquad\text{and}\qquad
c_{S_{nl}}=\frac{\sqrt{\tau_{p}P}{\zeta}_{{g}_{nl}}}{\tau_{p}P\left({\zeta}_{{g}_{nl}}+{\zeta}_{{u}_{nl}}\right)+1}.$
(7)

###### Proof.

Appendix A-A. ∎

Remark 1: The MMSE channel estimates in (6a) and (6b) are valid for the
coexistence of PUs/SUs with shared pilot sequences defined by (3). The MMSE
channel estimation with the conventional orthogonal pilots for cell-free
massive MIMO is reported in [8].

Due to TDD channel reciprocity, P-APs and S-APs utilize locally estimated
$\hat{f}_{mk}$ and $\hat{g}_{nl}$ as DL CSI to construct their precoders [33].
Furthermore, the actual channels can be written as

$\displaystyle{f}_{mk}$ $\displaystyle=$
$\displaystyle\hat{f}_{mk}+\epsilon_{f_{mk}}\qquad\text{and}\qquad{g}_{nl}=\hat{g}_{nl}+\epsilon_{g_{nl}},$
(8)

where $\epsilon_{f_{mk}}$ and $\epsilon_{g_{nl}}$ are the estimation errors,
which are independent of the corresponding channel estimates yielded from
orthogonality property of MMSE criterion [35].

For NOMA, the UL channel $(f_{mak})$ is also estimated locally at P-APs from
the pilots sent by PUs within PCs. Again, we assume that
$Q\\!\leq\\!\min(A,B)$ number of pilots is shared among PCs and SCs as

$\displaystyle\mathbf{\Phi}_{PA}=\left[\mathbf{\Phi};\tilde{\mathbf{\Phi}}_{PA}\right]\qquad\text{and}\qquad\mathbf{\Phi}_{SB}=\left[\mathbf{\Phi};\tilde{\mathbf{\Phi}}_{SB}\right],$
(9)

where $\mathbf{\Phi}\in\mathbb{C}^{Q\times\tau_{p}}$ denotes the shared
$Q$-pilot sequence among PCs and SCs. Then, the remaining $(A-Q)$ PCs or
$(B-Q)$ SCs are assigned with the pilot sequences
$\tilde{\mathbf{\Phi}}_{PA}\in\mathbb{C}^{(Q-A)\times\tau_{p}}$ and
$\tilde{\mathbf{\Phi}}_{SB}\in\mathbb{C}^{(Q-B)\times\tau_{p}}$, respectively.
Furthermore,
$\mathbf{\Phi}_{PA}=\left[\boldsymbol{\phi}_{P_{1}}^{T},\cdots,\boldsymbol{\phi}_{P_{a}}^{T},\cdots,\boldsymbol{\phi}_{P_{A}}^{T}\right]^{T}$
and
$\mathbf{\Phi}_{SB}=\left[\boldsymbol{\phi}_{S_{1}}^{T},\cdots,\boldsymbol{\phi}_{S_{b}}^{T},\cdots,\boldsymbol{\phi}_{S_{B}}^{T}\right]^{T}$,
where $\boldsymbol{\phi}_{P_{a}}\in\mathbb{C}^{1\times\tau_{p}}$ and
$\boldsymbol{\phi}_{S_{b}}\in\mathbb{C}^{1\times\tau_{p}}$ are the pilot
sequences assigned to the $a$th PC and the $b$th SC, satisfying
$\|\boldsymbol{\phi}_{P_{a}}\|^{2}=1$ and
$\|\boldsymbol{\phi}_{S_{b}}\|^{2}=1$ for $a\in\\{1,\cdots,A\\}$ and
$b\in\\{1,\cdots,B\\}$. The received pilot vector at the $m$th P-AP can be
written as

$\displaystyle\mathbf{y}_{P_{m}}$ $\displaystyle=$
$\displaystyle\sqrt{P_{p}}\sum\nolimits_{a=1}^{A}\sum\nolimits_{k=1}^{K}f_{mak}\boldsymbol{\phi}_{P_{a}}+\sqrt{P_{p}}\sum\nolimits_{b=1}^{B}\sum\nolimits_{l=1}^{L}v_{mbl}\boldsymbol{\phi}_{S_{b}}+\mathbf{n}_{P_{m}},$
(10)

where $\mathbf{n}_{P_{m}}\sim\mathcal{CN}(0,1)$ is the AWGN vector at the
$m$th P-AP. Then, by projecting $\boldsymbol{\phi}_{P_{a}}^{H}\in{\bf\Phi}$
into (10), we obtain a sufficient statistic to estimate $f_{mak}$ as

$\displaystyle y_{P_{ma}}$ $\displaystyle=$
$\displaystyle\boldsymbol{\phi}_{P_{a}}^{H}\mathbf{y}_{P_{m}}=\sqrt{P_{p}}\sum\nolimits_{i=1}^{K}f_{mak}+\sqrt{P_{p}}\sum\nolimits_{j=1}^{L}v_{maj}+n_{P_{m}}.$
(11)

By following steps similar to those in Appendix A-A, the MMSE estimate of
$f_{mak}$ can be derived as

$\displaystyle\hat{f}_{mak}$ $\displaystyle=$
$\displaystyle\frac{\mathbb{E}\\!\left[{y_{P_{ma}}^{*}f_{mak}}\right]}{\mathbb{E}\\!\left[{|y_{P_{ma}}|^{2}}\right]}y_{P_{ma}}=\frac{\sqrt{\tau_{p}P}{\zeta}_{{f}_{mak}}}{\tau_{p}P\left(\sum_{i=1}^{K}{\zeta}_{{f}_{mai}}+\sum_{j=1}^{L}{\zeta}_{{v}_{maj}}\right)+1}y_{P_{ma}}.$
(12)

Since $y_{P_{ma}}$ is Gaussian distributed, we have
$\hat{f}_{mak}\sim\mathcal{CN}(0,\alpha_{f_{mak}})$, where $\alpha_{f_{mak}}$
is given as

$\displaystyle\alpha_{f_{mak}}$ $\displaystyle=$
$\displaystyle\mathbb{E}\\!\left[{|\hat{f}_{mak}|^{2}}\right]=\frac{\tau_{p}P\zeta_{{f}_{mak}}^{2}}{\tau_{p}P\left(\sum_{i=1}^{K}\zeta_{{f}_{mai}}+\sum_{j=1}^{L}\zeta_{{v}_{maj}}\right)+1}.$
(13)

Then, the channel estimation error of $f_{mak}$ is defined as
$\epsilon_{f_{mak}}=f_{mak}-\hat{f}_{mak}\sim\mathcal{CN}(0,\zeta_{{f}_{mak}}-\alpha_{f_{mak}})$.

Remark 2: The parameters $y_{S_{nb}}$, $\hat{g}_{nbl}$, and $\alpha_{g_{nbl}}$
corresponding to the secondary system can be obtained by replacing subscripts
$\\{P,m,a,K,k,i\\}$ of (11), (12), and (13), respectively, by
$\\{S,n,b,L,l,j\\}$. Hence, the explicit presentation of UL channel estimates
at S-APs is omitted for the sake of brevity.

###  II-C AP Clustering for OMA

The P-APs/S-APs are clustered to serve a particular set of users based on the
locally estimated channel gains. We assume that only
$\mathcal{M}_{P}=\\{1,\cdots,M_{P}\\}$ and
$\mathcal{N}_{S}=\\{1,\cdots,N_{S}\\}$ sets of P-APs/S-APs are assigned to
$U_{P}(k)$ and $U_{S}(l)$, respectively [36]. At each P-AP/S-AP, a set of
transmission determination coefficients can be introduced based on the
estimated channels such that $\delta_{P_{mk}}=\delta_{S_{nl}}=1$ for
$m\in\mathcal{M}_{P}$, $n\in\mathcal{N}_{S}$, and
$\delta_{P_{mk}}=\delta_{S_{nl}}=0$, otherwise.

Remark 3: When the aforementioned AP clustering is adopted for the proposed
cell-free underlay spectrum sharing, the primary/secondary CCI can be reduced,
and hence, the achievable rates of both primary/secondary systems can be
boosted with respect to a conventional unclustered system as depicted in Fig.
5 in our numerical results in Section VIII.

###  II-D Signal model for OMA

Due to implementation simplicity and near optimal performance in large AP
regime [33], conjugate precoding is used at P-APs and S-APs to transmit
signals towards their respective users via the channel estimates in (6a) and
(6b). We can write the transmitted signals at the $m$th P-AP and the $n$th
S-AP as

$\displaystyle{x}_{P_{m}}$ $\displaystyle=$
$\displaystyle\sqrt{P_{P}}\sum\nolimits_{k=1}^{K}\delta_{P_{mk}}{\eta_{P_{mk}}^{1/2}\hat{f}_{mk}^{*}q_{P_{k}}}\quad\text{and}\quad{x}_{S_{n}}=\sqrt{P_{S}}\sum\nolimits_{l=1}^{L}\delta_{S_{nl}}{\eta_{S_{nl}}^{1/2}\hat{g}_{nl}^{*}q_{S_{l}}},$
(14)

where $\eta_{P_{mk}}$ and $\eta_{S_{nl}}$ are the power allocation
coefficients at the $m$th P-AP and the $n$th S-AP, respectively. Moreover,
$P_{P}$ and $P_{S}$ denote the maximum allowable transmit powers at each P-AP
and S-AP, respectively. Here, $\eta_{S_{nl}}$ is selected to satisfy the total
transmit power constraints given in (28). The signals intended to $U_{P}(k)$
and $U_{S}(l)$ are denoted by $q_{P_{k}}$ and $q_{S_{l}}$, respectively, and
they satisfy $\mathbb{E}\\!\left[{|q_{P_{k}}|^{2}}\right]=1$ and
$\mathbb{E}\\!\left[{|q_{S_{l}}|^{2}}\right]=1$. Then, the received signals at
$U_{P}(k)$ and $U_{S}(l)$ can be written as

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!{r}_{P_{k}}$ $\displaystyle=$
$\displaystyle\sum\nolimits_{m=1}^{M}{f_{mk}{x}_{P_{m}}}+\sum\nolimits_{n=1}^{N}{u_{nk}{x}_{S_{n}}}+n_{P_{k}}\quad\text{and}\quad{r}_{S_{l}}=\sum\nolimits_{n=1}^{N}{g_{nl}{x}_{S_{n}}}+\sum\nolimits_{m=1}^{M}{v_{ml}{x}_{P_{m}}}+n_{S_{l}},$
(15)

where $n_{P_{k}}\sim\mathcal{CN}(0,1)$ and $n_{S_{l}}\sim\mathcal{CN}(0,1)$
are AWGN at $U_{P}(k)$ and $U_{S}(l)$, respectively. Here, ${x}_{P_{m}}$ and
${x}_{S_{n}}$ are given in (14). We can rearrange the received signal at
$U_{P}(k)$ in (15) as

$\displaystyle{r}_{P_{k}}$ $\displaystyle=$
$\displaystyle\sqrt{P_{P}}\sum\nolimits_{m=1}^{M}\delta_{P_{mk}}{\eta_{P_{mk}}^{1/2}f_{mk}\hat{f}_{mk}^{*}q_{P_{k}}}+\sqrt{P_{P}}\sum\nolimits_{m=1}^{M}\sum\nolimits_{i\neq
k}^{K}\delta_{P_{mi}}{\eta_{P_{mi}}^{1/2}f_{mk}\hat{f}_{mi}^{*}q_{P_{i}}}$
(16)
$\displaystyle+\sqrt{P_{S}}\sum\nolimits_{n=1}^{N}\sum\nolimits_{j=1}^{L}\delta_{S_{nj}}{\eta_{S_{nj}}^{1/2}u_{nk}\hat{g}_{nj}^{*}q_{S_{j}}}+n_{P_{k}},$

where the first term represents the desired signal component at $U_{P}(k)$,
while the inter-system interference caused by beamforming uncertainty of
conjugate precoding is captured by the second term. The third term accounts
for the intra-system interference yielded from the secondary CCI. Similarly,
the received signal at $U_{S}(l)$ in (15) can be rewritten as

$\displaystyle{r}_{S_{l}}$ $\displaystyle=$
$\displaystyle\sqrt{P_{S}}\sum\nolimits_{n=1}^{N}\delta_{S_{nl}}{\eta_{S_{nl}}^{1/2}g_{nl}\hat{g}_{nl}^{*}q_{S_{l}}}+\sqrt{P_{S}}\sum\nolimits_{n=1}^{N}\sum\nolimits_{j\neq
l}^{L}\delta_{S_{nj}}{\eta_{S_{nj}}^{1/2}g_{nl}\hat{g}_{nj}^{*}q_{S_{j}}}$
(17)
$\displaystyle+\sqrt{P_{P}}\sum\nolimits_{m=1}^{M}\sum\nolimits_{i=1}^{K}\delta_{P_{mi}}{\eta_{P_{mi}}^{1/2}v_{ml}\hat{f}_{mi}^{*}q_{P_{i}}}+n_{S_{l}}.$

###  II-E Signal model for NOMA

Again, P-APs employ conjugate precoding by using the channel estimates in
(12). The transmitted signal at the $m$th P-AP can be written as

$\displaystyle{x}_{P_{m}}$ $\displaystyle=$
$\displaystyle\sqrt{P_{P}}\sum\nolimits_{a=1}^{A}\sum\nolimits_{i=1}^{K}{\eta_{P_{mai}}^{1/2}\hat{f}_{mai}^{*}q_{P_{ai}}},$
(18)

where $\eta_{P_{mai}}$ is the transmit power control coefficient at the $m$th
P-AP. Here, $\mathbb{E}\\!\left[{|q_{P_{ai}}|^{2}}\right]=1$, where
$q_{P_{ai}}$ is the signal intended for $U_{P}(a,i)$. Thus, the received
signal at $U_{P}(a,k)$ can be written as

$\displaystyle{r}_{P_{ak}}$ $\displaystyle=$
$\displaystyle\sqrt{P_{P}}\sum\nolimits_{m=1}^{M}{\eta_{P_{mak}}^{1/2}f_{mak}\hat{f}_{mak}^{*}q_{P_{ak}}}+\sqrt{P_{P}}\sum\nolimits_{m=1}^{M}\sum\nolimits_{i\neq
k}^{K}{\eta_{P_{mai}}^{1/2}f_{mak}\hat{f}_{mai}^{*}q_{P_{ai}}}$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!+\\!\sqrt{P_{P}}\sum\nolimits_{m=1}^{M}\\!\sum\nolimits_{a^{\prime}\neq
a}^{A}\\!\sum\nolimits_{i=1}^{K}\\!\eta_{P_{ma^{\prime}i}}^{1/2}f_{mak}\hat{f}_{ma^{\prime}i}^{*}q_{P_{a^{\prime}i}}\\!+\\!\sqrt{P_{S}}\sum\nolimits_{n=1}^{N}\\!\sum\nolimits_{b=1}^{B}\\!\sum\nolimits_{j=1}^{L}\\!{\eta_{S_{nbj}}^{1/2}u_{nak}\hat{g}_{nbj}^{*}q_{S_{bj}}}\\!+\\!n_{P_{ak}},$

where $n_{P_{ak}}\sim\mathcal{CN}(0,1)$ is the AWGN at $U_{P}(a,k)$. To apply
the power-domain NOMA, we assume that the users in the $a$th PC are ordered
based on the effective channel strength as [24, 12, 34]

$\displaystyle\mathbb{E}\\!\left[{\left|\sum\nolimits_{m=1}^{M}\hat{f}_{ma1}\right|^{2}}\right]\geq\cdots\geq\mathbb{E}\\!\left[{\left|\sum\nolimits_{m=1}^{M}\hat{f}_{mak}\right|^{2}}\right]\geq\cdots\geq\mathbb{E}\\!\left[{\left|\sum\nolimits_{m=1}^{M}\hat{f}_{maK}\right|^{2}}\right].$
(20)

In power-domain NOMA, higher transmit powers are allocated for the users with
weaker channel conditions. Thus, the transmit powers are ordered as [24, 12,
34]

$\displaystyle P_{P_{a1}}\leq\cdots\leq P_{P_{ak}}\leq\cdots\leq P_{P_{aK}},$
(21)

where $P_{P_{ak}}=P_{P}\eta_{P_{mak}}$. Consequently, $U_{P}(a,k)$ aims to
decode the signal intended for $U_{P}(a,i)$ for $\forall i\geq k$ provided
that $U_{P}(a,k)$ can decode its own signal. Hence, $U_{P}(a,k)$ may
successively cancel the intra-cluster interference from $U_{P}(a,i)$ before
decoding its own signal for $\forall i\geq k$, and the residual interference
due to SIC error propagation must also captured. Moreover, $U_{P}(a,k)$ treats
the signals for users $\forall i<k$ as interference [24, 12, 34]. To this end,
the received signal at $U_{P}(a,k)$ upon imperfect SIC with error propagation
can be written as

$\displaystyle{r}_{P_{ak}}$ $\displaystyle=$
$\displaystyle\underbrace{\sqrt{P_{P}}\sum\nolimits_{m=1}^{M}\eta_{P_{mak}}^{1/2}f_{mak}\hat{f}_{mak}^{*}q_{P_{ak}}}_{\text{Desired
signal}}+\underbrace{\sqrt{P_{P}}\sum\nolimits_{m=1}^{M}\sum\nolimits_{i=1}^{k-1}\eta_{P_{mai}}^{1/2}f_{mai}\hat{f}_{mai}^{*}q_{P_{ai}}}_{\text{Intra-
cluster interference after SIC}}$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!+\underbrace{\sqrt{P_{P}}\sum\nolimits_{m=1}^{M}\sum\nolimits_{i=k+1}^{K}\eta_{P_{mai}}\left(f_{mak}\hat{f}_{mai}^{*}q_{P_{ai}}-\mathbb{E}\\!\left[{f_{mak}\hat{f}_{mai}^{*}}\right]\hat{q}_{P_{ai}}\right)}_{\text{Error
propagation due to imperfect SIC}}$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!+\\!\underbrace{\sqrt{P_{P}}\\!\sum\nolimits_{m=1}^{M}\\!\sum\nolimits_{a^{\prime}\neq
a}^{A}\\!\sum\nolimits_{i=1}^{K}\\!\eta_{P_{ma^{\prime}i}}^{1/2}f_{mak}\hat{f}_{ma^{\prime}i}^{*}q_{P_{a^{\prime}i}}}_{\text{Intra-
system
interference}}\\!+\\!\underbrace{\sqrt{P_{S}}\\!\sum\nolimits_{n=1}^{N}\\!\sum\nolimits_{b=1}^{B}\\!\sum\nolimits_{j=1}^{L}\\!\eta_{S_{nbj}}^{1/2}u_{nak}\hat{g}_{nbj}^{*}q_{S_{bj}}}_{\text{Inter-
system interference}}\\!+\\!\underbrace{n_{P_{ak}}}_{\text{AWGN}}.$

Remark 4: In our proposed cell-free NOMA-aided underlay spectrum sharing, the
perfect SIC is not feasible due to intra-cluster pilot contamination, intra-
system interference, inter-system interference, channel estimation errors, and
statistical CSI knowledge at the users. Thus, the residual interference caused
by imperfect SIC needs to be modeled. The third term in (II-E) captures the
error propagation due to imperfect SIC in which $\hat{q}_{P_{ai}}$ is the
estimate of ${q}_{P_{ai}}$. Since ${q}_{P_{ai}}$ is Gaussian distributed,
$\hat{q}_{P_{ai}}$ and ${q}_{P_{ai}}$ are assumed to be jointly Gaussian
distributed with a normalized correlation coefficient $\vartheta_{P_{ai}}$ as
[12]

$\displaystyle{q}_{P_{ai}}=\vartheta_{P_{ai}}\hat{q}_{P_{ai}}+e_{P_{ai}},$
(23)

where $\hat{q}_{P_{ai}}\sim\mathcal{CN}(0,1)$,
$e_{P_{ai}}\sim\mathcal{CN}(0,\sigma_{e_{P_{ai}}}^{2}/(1+\sigma_{e_{P_{ai}}}^{2}))$,
and $\vartheta_{P_{ai}}=1/\sqrt{1+\sigma_{e_{P_{ai}}}^{2}}$. Furthermore,
$\hat{q}_{P_{ai}}$ and $e_{P_{ai}}$ are statistically independent. Thus, the
third term in (II-E) can be used to capture the residual interference caused
by error propagation of imperfect SIC when evaluating the SINR and achievable
rate.

Remark 5: The received signal at $U_{S}(b,l)$ for the secondary system
$({r}_{S_{bl}})$ and the signal intended for $(U_{S}(b,j)$, ${q}_{S_{bj}})$
can be obtained by replacing the subscripts $\\{P,M,m,A,a,K,k,i\\}$ of (II-E)
and (23), respectively, by $\\{S,N,n,B,b,L,l,j\\}$.

###  II-F Secondary Transmit Power Constraints for OMA

We constrain the transmit power at each S-AP to guarantee that the secondary
CCI inflected at PUs falls below the PIT of each PU. Thus, we define the total
transmit power at the $n$th S-AP as

$\displaystyle
P_{S_{n}}\triangleq\sum\nolimits_{l=1}^{L}P_{S}\delta_{S_{nl}}\eta_{S_{nl}}\quad\text{for}\quad
n\in\\{1,\cdots,N\\}\quad\text{and}\quad l\in\\{1,\cdots,L\\},$ (24)

where $\sum_{l=1}^{L}\delta_{S_{nl}}\eta_{S_{nl}}\leq 1$. Moreover, the
transmit power allocation coefficient at the $n$th S-AP for $U_{S}(l)$ is
represented by $\eta_{S_{nl}}$. Thus, the CCI received at $U_{P}(k)$ from all
S-APs can be written as

$\displaystyle y_{k}$ $\displaystyle=$
$\displaystyle\sum\nolimits_{n=1}^{N}u_{nk}x_{S_{n}}=\sqrt{P_{S}}\sum\nolimits_{n=1}^{N}\sum\nolimits_{l=1}^{L}\delta_{S_{nl}}\eta_{S_{nl}}^{1/2}u_{nk}\hat{g}_{nl}^{*}q_{S_{l}},$
(25)

where $x_{S_{n}}$ is the transmitted signal at the $n$th S-AP and defined in
(14).

Proposition 2: The total average secondary CCI power ($P_{I_{k}}$) inflicted
at $U_{P}(k)$ is given by

$\displaystyle P_{I_{k}}$ $\displaystyle=$
$\displaystyle\mathbb{E}\\!\left[{\left|y_{k}\right|^{2}}\right]={P_{S}}\underbrace{\sum\nolimits_{n=1}^{N}\sum\nolimits_{l=1}^{L}\delta_{S_{nl}}\eta_{S_{nl}}\mathbb{E}\\!\left[{\left|{u_{nk}\hat{g}_{nl}^{*}}\right|^{2}}\right]}_{Z_{k}},$
(26)

where $Z_{k}$ is defined as

$\displaystyle Z_{k}$ $\displaystyle=$
$\displaystyle\sum\nolimits_{n=1}^{N}\sum\nolimits_{l=1}^{L}\delta_{S_{nl}}\eta_{S_{nl}}\rho_{g_{nl}}\zeta_{{u}_{nk}}+\sum\nolimits_{n=1}^{N}\delta_{S_{nk}}\eta_{S_{nk}}\rho_{u_{nk}}^{2},$
(27)

where $\rho_{g_{nl}}\triangleq\sqrt{\tau_{p}P}c_{S_{nl}}\zeta_{{g}_{nl}}$ and
$\rho_{u_{nk}}\triangleq\sqrt{\tau_{p}P}c_{S_{nk}}\zeta_{{u}_{nk}}$.

###### Proof.

Appendix A-B. ∎

Thus, we can give the secondary transmit power constraint at the $n$th S-AP as
follows:

$\displaystyle
P_{S_{n}}=\min{\left(P_{S},{I_{T_{1}}}\big{/}{Z_{1}},\cdots,{I_{T_{k}}}\big{/}{Z_{k}},\cdots,{I_{T_{K}}}\big{/}{Z_{K}}\right)},$
(28)

where $I_{T_{k}}$ is the interference threshold of $U_{P}(k)$.

###  II-G Secondary transmit power control for NOMA

As per Section II-F, to ensure that the performance of primary system is not
hindered by the secondary system, we constrain the transmit power of S-APs.
The secondary CCI received at $U_{P}(a,k)$ is given by

$\displaystyle y_{ak}$ $\displaystyle=$
$\displaystyle\sum\nolimits_{n=1}^{N}u_{nak}x_{S_{nb}}=\sqrt{P_{S}}\sum\nolimits_{n=1}^{N}\sum\nolimits_{b=1}^{B}\sum\nolimits_{l=1}^{L}\eta_{S_{nbl}}^{1/2}u_{nak}\hat{g}_{nbl}^{*}q_{S_{bl}},$
(29)

where $x_{S_{nb}}$ is the transmit signal intended for $L$ users in the $b$th
SC. Then, the total average secondary CCI ($P_{I_{ak}}$) inflicted at
$U_{P}(a,k)$ can be derived as

$\displaystyle P_{I_{ak}}$ $\displaystyle=$
$\displaystyle\mathbb{E}\\!\left[{\left|y_{ak}\right|^{2}}\right]={P_{S}}\underbrace{\sum\nolimits_{n=1}^{N}\sum\nolimits_{b=1}^{B}\sum\nolimits_{l=1}^{L}\eta_{S_{nl}}\mathbb{E}\\!\left[{\left|{u_{nak}\hat{g}_{nbl}^{*}}\right|^{2}}\right]}_{Z_{ak}},$
(30)

where $Z_{ak}$ can be derived by following steps similar to those in Appendix
A-B as

$\displaystyle Z_{ak}$ $\displaystyle=$
$\displaystyle\sum\nolimits_{n=1}^{N}\sum\nolimits_{b=1}^{B}\sum\nolimits_{l=1}^{L}\eta_{S_{nbl}}\alpha_{g_{nbl}}\zeta_{{u}_{nak}}+\sum\nolimits_{n=1}^{N}\sum\nolimits_{l=1}^{L}\eta_{S_{nal}}\alpha_{f_{mak}}^{2}\left(\frac{\zeta_{{u}_{uak}}\zeta_{{g}_{nal}}}{\zeta_{{f}_{mak}}}\right)^{2}.$
(31)

Thus, we derive the secondary transmit power constraint at the $n$th S-AP as

$\displaystyle
P_{S_{n}}=\min{\left(P_{S},{I_{T_{11}}}\big{/}{Z_{11}},\cdots,{I_{T_{ak}}}\big{/}{Z_{ak}},\cdots,{I_{T_{AK}}}\big{/}{Z_{AK}}\right)},$
(32)

where $I_{T_{ak}}$ is the interference threshold of $U_{P}(a,k)$.

##  III Achievable Rate analysis of OMA-Aided Underlay Spectrum Sharing

###  III-A Achievable rate analysis for the primary system

When the P-APs/S-APs do not beamform DL pilots for the acquisition of DL CSI,
the PUs/SUs are unaware of the instantaneous DL channel coefficients. Thus,
the PUs/SUs must relay on the long-term/statistical DL channel coefficients
for signal detection [33]. This is a typical scenario in TDD-based co-located
massive MIMO in which only UL pilots are used to estimate UL channels at the
P-APs/S-APs, and the instantaneous DL channel coefficients can be approximated
by their statistical counterparts thanks to channel hardening property. To
this end, the signal received at $U_{P}(k)$ can be rearranged as

$\displaystyle\\!\\!\\!\\!r_{P_{k}}\\!$ $\displaystyle=$
$\displaystyle\\!\underbrace{\sqrt{P_{P}}\mathbb{E}\\!\left[{\sum\nolimits_{m=1}^{M}\\!\delta_{P_{mk}}{\eta_{P_{mk}}^{1/2}f_{mk}\hat{f}_{mk}^{*}}}\right]q_{P_{k}}}_{\text{Desired
signal}}$ (33)
$\displaystyle\\!\\!\\!\\!+\underbrace{\sqrt{P_{P}}\left({\sum\nolimits_{m=1}^{M}\\!\delta_{P_{mk}}{\eta_{P_{mk}}^{1/2}f_{mk}\hat{f}_{mk}^{*}}}\\!-\\!{\mathbb{E}\\!\left[{\sum\nolimits_{m=1}^{M}\\!\delta_{P_{mk}}{\eta_{P_{mk}}^{1/2}f_{mk}\hat{f}_{mk}^{*}}}\right]}\right)q_{P_{k}}}_{\text{Detection
uncertainty}}$
$\displaystyle\\!\\!\\!\\!+\\!\underbrace{\sqrt{P_{P}}\sum\nolimits_{i\neq
k}^{K}{\sum\nolimits_{m=1}^{M}\delta_{P_{mi}}{\eta_{P_{mi}}^{1/2}f_{mk}\hat{f}_{mi}^{*}}}q_{P_{i}}}_{\text{Inter-
system interference caused by beamforminmg
uncertainty}}\\!+\\!\underbrace{\sqrt{P_{S}}\sum\nolimits_{j=1}^{L}{\sum\nolimits_{n=1}^{N}\delta_{S_{nj}}{\eta_{S_{nj}}^{1/2}u_{nk}\hat{g}_{nj}^{*}}}q_{S_{j}}}_{\text{Inter-
system interference caused by secondary
CCI}}\\!+\\!\underbrace{n_{P_{k}}}_{\text{AWGN}}\\!.$

In (33), the effective noise consists of the sum of (i) interference caused by
detection uncertainty, (ii) intra-system interference due to beamforming
uncertainty, (iii) inter-system interference due to secondary CCI, and (iv)
AWGN. The desired signal term and effective noise are uncorrelated. Due to the
law of large number, the latter can be treated as worst-case independently
distributed Gaussian noise [37, 33]. Theorem 1: The effective SINR at
$U_{P}(k)$ is given by

$\displaystyle\gamma_{P_{k}}=\frac{{P_{P}}\left|\mathbb{E}\\!\left[{\sum\nolimits_{m=1}^{M}\delta_{P_{mik}}{\eta_{P_{mk}}^{1/2}f_{mk}\hat{f}_{mk}^{*}}}\right]\right|^{2}}{P_{P}\mathbb{V}\mathrm{ar}\\!\left[{\sum\nolimits_{m=1}^{M}\delta_{P_{mk}}{\eta_{P_{mk}}^{1/2}f_{mk}\hat{f}_{mk}^{*}}}\right]+\sum\nolimits_{j=1}^{2}{I_{P}}_{j}+\mathbb{E}\\!\left[{\left|n_{P_{k}}\right|^{2}}\right]},$
(34)

where ${I_{P}}_{j}$ for $j\in\\{1,2\\}$ can be defined as

$\displaystyle\\!\\!\\!\\!\\!{I_{P}}_{1}\\!=\\!{P_{P}}\mathbb{E}\\!\left[{\left|\sum\nolimits_{i\neq
k}^{K}\\!{\sum\nolimits_{m=1}^{M}\\!\delta_{P_{mi}}{\eta_{P_{mi}}^{1/2}f_{mk}\hat{f}_{mi}^{*}}}\right|^{2}}\right]\\!\quad\text{and}\quad{I_{P}}_{2}\\!=\\!{P_{S}}\mathbb{E}\\!\left[{\left|\sum\nolimits_{j=1}^{L}\\!{\sum\nolimits_{n=1}^{N}\\!\delta_{S_{nj}}{\eta_{S_{nj}}^{1/2}u_{nk}\hat{g}_{nj}^{*}}}\right|^{2}}\right]\\!.{}$
(35)

By evaluating the expectation and variance terms in (34) and (35), the SINR is
given by

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\gamma_{P_{k}}\\!=\\!\frac{P_{P}\left(\sum\nolimits_{m=1}^{M}\delta_{P_{mk}}\eta_{P_{mk}}^{1/2}\rho_{f_{mk}}\right)^{2}}{P_{P}\\!\sum\nolimits_{m=1}^{M}\sum\nolimits_{i=1}^{K}\delta_{P_{mi}}\eta_{P_{mi}}\rho_{f_{mi}}\zeta_{f_{mk}}\\!+\\!{P_{S}}\\!\sum\nolimits_{n=1}^{N}\sum\nolimits_{j=1}^{L}\delta_{S_{nj}}\eta_{S_{nj}}\rho_{g_{nj}}\zeta_{{u}_{nk}}\\!+\\!{P_{S}}\\!\sum\nolimits_{n=1}^{N}\delta_{S_{nk}}\eta_{S_{nk}}\rho_{{u}_{nk}}^{2}\\!+\\!1}\\!,$
(36)

where $\rho_{f_{mk}}\triangleq\sqrt{\tau_{p}P}c_{P_{mk}}\zeta_{{f}_{mk}}$.

###### Proof.

Appendix B. ∎

Then, we define the achievable rate of $U_{P}(k)$ as follows:

$\displaystyle
R_{P_{k}}=\left({\left(\tau_{c}-\tau_{p}\right)}/{\tau_{c}}\right)\mathrm{log}_{2}\left(1+\gamma_{P_{k}}\right),$
(37)

where the effective portion of coherence interval for payload data
transmission is captured by the pre-log factor
$({\left(\tau_{c}-\tau_{p}\right)}/{\tau_{c}})$ and $\gamma_{P_{k}}$ is
defined in (36).

###  III-B Achievable rate definition for the secondary system

By following Theorem 1, we derive the achievable rate of $U_{S}(l)$ as
follows:

$\displaystyle\\!\\!\\!\\!\\!\\!R_{S_{l}}\\!=\\!\left({\left(\tau_{c}\\!-\\!\tau_{p}\right)}/{\tau_{c}}\right)\mathrm{log}_{2}\left(1+\gamma_{S_{l}}\right)\,\,\;\text{and}\,\,\;\gamma_{S_{l}}\\!=\\!\frac{{P_{S}}\left|\mathbb{E}\\!\left[{\sum\nolimits_{n=1}^{N}\delta_{S_{nl}}{\eta_{S_{nl}}^{1/2}g_{nl}\hat{g}_{nl}^{*}}}\right]\right|^{2}}{P_{S}\mathbb{V}\mathrm{ar}\\!\left[{\sum\nolimits_{n=1}^{N}\delta_{S_{nl}}{\eta_{S_{nl}}^{1/2}g_{nl}\hat{g}_{nl}^{*}}}\right]\\!+\\!\sum\nolimits_{j=1}^{2}{I_{S}}_{j}\\!+\\!\mathbb{E}\\!\left[{\left|w_{S_{l}}\right|^{2}}\right]},$
(38)

where ${I_{S}}_{j}$ for $j\in\\{1,2\\}$ can be defined as

$\displaystyle\\!\\!\\!\\!\\!\\!I_{S_{1}}\\!=\\!{P_{S}}\mathbb{E}\\!\left[{\left|\sum\nolimits_{j\neq
l}^{L}\\!{\sum\nolimits_{n=1}^{N}\\!\delta_{S_{nj}}{\eta_{S_{nj}}^{1/2}g_{nl}\hat{g}_{nj}^{*}}}\right|^{2}}\right]\quad\text{and}\quad
I_{S_{2}}\\!=\\!P_{P}\mathbb{E}\\!\left[{\left|\sum\nolimits_{i=1}^{K}\\!{\sum\nolimits_{m=1}^{M}\delta_{P_{mi}}\\!\eta_{P_{mi}}^{1/2}v_{ml}\hat{f}_{mi}^{*}}\right|^{2}}\right].$
(39)

Thus, by deriving the expectation and variance term in (38) in similar manner
to (34), the SINR of $U_{S}(l)$ can be computed as

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\gamma_{S_{l}}\\!=\\!\frac{{P_{S}}\left(\sum\nolimits_{n=1}^{N}\delta_{S_{nl}}\eta_{S_{nl}}^{1/2}\rho_{g_{nl}}\right)^{2}}{P_{S}\\!\sum\nolimits_{n=1}^{N}\sum\nolimits_{j=1}^{L}\delta_{S_{nj}}\eta_{S_{nj}}\rho_{g_{nj}}\zeta_{g_{nl}}\\!+\\!P_{P}\\!\sum\nolimits_{m=1}^{M}\sum\nolimits_{i=1}^{K}\delta_{P_{mi}}\eta_{P_{mi}}\rho_{f_{mi}}\zeta_{{v}_{ml}}\\!+\\!P_{P}\\!\sum\nolimits_{m=1}^{M}\delta_{P_{ml}}\eta_{P_{ml}}\rho_{{v}_{ml}}^{2}\\!+\\!1}\\!.$
(40)

##  IV Transmit Power Control

In cell-free massive MIMO, the user-fairness must be guaranteed in terms of
the achievable rate in order to provide a uniform quality-of-service to all
users. To this end, max-min power control algorithms have been shown to be
optimal in the sense of user-fairness in presence of near-far effects [33, 38,
39, 40]. For our proposed system model, a multi-objective optimization problem
(MOOP) [41] is most suited as both primary and secondary systems operate
simultaneously. Furthermore, the secondary transmit power constraints in (28)
must also be considered when formulating this MOOP.

We compute the optimal power allocation coefficients of P-APs/S-APs to
maximize the minimum achievable DL rate among all PUs/SUs by invoking the max-
min optimization criterion [33]. Since the rates in (37) and (38) are
monotonically increasing functions of their arguments, we can equivalently
replace $R_{P_{k}}$ and $R_{S_{l}}$ by $\gamma_{P_{k}}$ and $\gamma_{S_{l}}$,
which are given in (36) and (40), respectively. Thus, we formulate a max-min
transmit power control problem by introducing a common SINR
$\left(\lambda_{\gamma}\right)$ for primary/secondary systems and by defining
slack variables $\beta_{P_{mk}}\triangleq\eta_{P_{mk}}^{1/2}$ and
$\beta_{S_{nl}}\triangleq\eta_{S_{nl}}^{1/2}$ as

$\displaystyle\underset{\beta_{P_{mk}},\beta_{S_{nl}}\forall{m,k,n,l}}{\text{maximize}}\quad$
$\displaystyle\left(\gamma_{P_{k}}\right)^{w_{P}}\left(\gamma_{S_{l}}\right)^{w_{S}}=\lambda_{\gamma},$
(41a) subject to $\displaystyle
C_{1}:\gamma_{P_{k}}\geq\lambda_{\gamma}\quad\text{and}\quad
C_{2}:\gamma_{S_{l}}\geq\lambda_{\gamma},$ (41e) $\displaystyle
C_{3}:\sum\nolimits_{k=1}^{K}\delta_{P_{mk}}\beta_{P_{mk}}^{2}\rho_{f_{mk}}\leq
1\quad\text{and}\quad\sum\nolimits_{l=1}^{L}\delta_{S_{nl}}\beta_{S_{nl}}^{2}\rho_{g_{nl}}\leq
1,$ $\displaystyle C_{4}:P_{S_{n}}\leq{I_{T_{k}}}/{Z_{k}},$ $\displaystyle
C_{5}:0\leq\beta_{P_{mk}}\quad\text{and}\quad 0\leq\beta_{S_{nl}},$

where $w_{P}$ and $w_{S}$ are the priorities assigned for primary and
secondary achievable SINR, respectively. Moreover, $C_{3}$ is obtained by
invoking the maximum allowable transmit power constraints at the $m$th P-AP
and the $l$th S-AP as follows:

$\displaystyle\mathbb{E}\\!\left[{|x_{P_{m}}|^{2}}\right]\leq
P_{P}\;\Rightarrow\;\mathbb{E}\\!\left[{\left|\sqrt{P_{P}}\sum\nolimits_{k=1}^{K}\delta_{P_{mk}}{\eta_{P_{mk}}^{1/2}\hat{f}_{mk}^{*}q_{P_{k}}}\right|^{2}}\right]\leq
P_{P}$ (42)
$\displaystyle\sum\nolimits_{k=1}^{K}\delta_{P_{mk}}\eta_{P_{mk}}\mathbb{E}\\!\left[{|\hat{f}_{mk}^{*}q_{P_{k}}|^{2}}\right]\leq
1\;\;\Rightarrow\;\;\sum\nolimits_{k=1}^{K}\delta_{P_{mk}}\eta_{P_{mk}}\rho_{f_{mk}}\leq
1.$

It can be shown that
$\sum_{l=1}^{L}\delta_{S_{nl}}\eta_{S_{nl}}\rho_{g_{nl}}\leq 1$. We adopt the
secondary transmit power constraint in (28) to obtain $C_{4}$. Since the
objective functions in (41) are quasi-concave functions, we can show that the
underlaying optimization problem is also quasi-concave [8]. Thus, an optimal
solution can be found by using the Bisection method as shown in Algorithm 1.

Algorithm 1 Bisection Algorithm

0: Path-losses between all P-APs/S-APs and PUs/SUs, and the average transmit
powers of PUs/SUs.

0: The power control coefficients $\eta_{P_{mk}}$ for $m\in\\{1,\cdots,M\\}$,
$k\in\\{1,\cdots,K\\}$, and $\eta_{S_{nl}}$ for $n\in\\{1,\cdots,N\\}$,
$l\in\\{1,\cdots,L\\}$. Initialization: Define an initial region for the
objective functions by choosing appropriate values for $\lambda_{min}$ and
$\lambda_{max}$. Choose a tolerance $\epsilon>0$.

1: while $\lambda_{max}-\lambda_{min}>\epsilon$ do

2: Calculate,
$\lambda_{\gamma}={\left(\lambda_{max}+\lambda_{min}\right)}/{2}$.

3: Solve the convex feasibility problem, which can be formulated as
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\|\boldsymbol{V}_{\gamma_{P_{k}}}\|\leq\frac{1}{\sqrt{\lambda}}\left(\sum\nolimits_{m=1}^{M}\delta_{P_{mk}}\beta_{P_{mk}}\rho_{f_{mk}}\right)\qquad\text{and}\qquad\|\boldsymbol{V}_{\gamma_{S_{l}}}\|\leq\frac{1}{\sqrt{\lambda}}\left(\sum\nolimits_{n=1}^{N}\delta_{S_{nl}}\beta_{S_{nl}}\rho_{g_{nl}}\right),$
(43) which is subjected to $C_{2}$, $C_{3}$, and $C_{4}$ given in (41e),
(41e), and (41e), respectively. Moreover,
$\boldsymbol{V}_{\gamma_{P_{k}}}\triangleq\left[\mathbf{v}_{P_{1}}^{T},\;\frac{P_{S}}{P_{P}}\mathbf{v}_{P_{2}}^{T},\;\frac{P_{S}}{P_{P}}\mathbf{v}_{P_{3}}^{T},\;\frac{1}{\sqrt{P_{P}}}\right]$,
where
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\mathbf{v}_{P_{1}}$
$\displaystyle\triangleq$
$\displaystyle\left[\delta_{P_{11}}\beta_{P_{11}}\sqrt{\rho_{f_{11}}\zeta_{{f}_{1k}}},\cdots,\delta_{P_{MK}}\beta_{P_{MK}}\sqrt{\rho_{f_{MK}}\zeta_{{f}_{Mk}}}\,\right]^{T}\\!\\!\\!,$
(44a)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\mathbf{v}_{P_{2}}$
$\displaystyle\triangleq$
$\displaystyle\left[\delta_{S_{11}}\beta_{S_{11}}\sqrt{\rho_{g_{11}}\zeta_{{u}_{1k}}},\cdots,\delta_{S_{NL}}\beta_{S_{NL}}\sqrt{\rho_{g_{NL}}\zeta_{{u}_{Nk}}}\,\right]^{T}\\!\\!\\!,$
(44b)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\mathbf{v}_{P_{3}}$
$\displaystyle\triangleq$
$\displaystyle\left[\delta_{S_{1k}}\beta_{S_{1k}}c_{P_{1k}}\zeta_{{u}_{1k}},\cdots,\delta_{S_{Nk}}\beta_{S_{Nk}}c_{P_{Nk}}\zeta_{{u}_{Nk}}\,\right]^{T}\\!\\!\\!,$
(44c) and
$\boldsymbol{V}_{\gamma_{S_{l}}}\triangleq\left[\mathbf{v}_{S_{1}}^{T},\;\frac{P_{P}}{P_{S}}\mathbf{v}_{S_{2}}^{T},\;\frac{P_{P}}{P_{S}}\mathbf{v}_{S_{3}}^{T},\;\frac{1}{\sqrt{P_{S}}}\right]$,
where
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\mathbf{v}_{S_{1}}$
$\displaystyle\triangleq$
$\displaystyle\left[\delta_{S_{11}}\beta_{S_{11}}\sqrt{\rho_{g_{11}}\zeta_{{g}_{1l}}},\cdots,\delta_{S_{NL}}\beta_{S_{NL}}\sqrt{\rho_{g_{NL}}\zeta_{{g}_{Nl}}}\,\right]^{T}\\!\\!\\!,$
(45a)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\mathbf{v}_{S_{2}}$
$\displaystyle\triangleq$
$\displaystyle\left[\delta_{P_{11}}\beta_{P_{11}}\sqrt{\rho_{f_{11}}\zeta_{{v}_{1l}}},\cdots,\delta_{P_{MK}}\beta_{P_{MK}}\sqrt{\rho_{f_{MK}}\zeta_{{v}_{Ml}}}\,\right]^{T}\\!\\!\\!,$
(45b)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\mathbf{v}_{S_{3}}$
$\displaystyle\triangleq$
$\displaystyle\left[\delta_{P_{1l}}\beta_{P_{1l}}c_{S_{1l}}\zeta_{{v}_{1l}},\cdots,\delta_{P_{Ml}}\beta_{P_{Ml}}c_{S_{Ml}}\zeta_{{v}_{Ml}}\,\right]^{T}\\!\\!\\!.$
(45c)

4: If the status of the problem is feasible, then set
$\lambda_{min}=\lambda_{\gamma}$, otherwise set
$\lambda_{max}=\lambda_{\gamma}$.

5: Stop if $\lambda_{max}-\lambda_{min}<\epsilon$. Otherwise go to Step 2.

6: end while

7: return $\eta_{P_{mk}}=\beta_{P_{mk}}^{2}$ and
$\eta_{S_{nl}}=\beta_{S_{nl}}^{2}$for $m\in\\{1,\cdots,M\\}$,
$k\in\\{1,\cdots,K\\}$, $n\in\\{1,\cdots,N\\}$, and $l\in\\{1,\cdots,L\\}$.

##  V The implication of DL pilot transmission

In our achievable rate analysis in Section III-A, the PUs/SUs are assumed to
be unaware of DL channel estimates, and this is a typical assumption in co-
located massive MIMO literature [33]. It is aimed at minimizing the pilot
overhead to preserve system scalability. Thus, the users adopt long-term DL
channel statistics to decode the received signals as the DL channel
coefficients can be tightly approximated by their average counterparts by
virtue of channel hardening [37]. Nevertheless, it has been shown in [42] that
channel hardening in cell-free massive MIMO occurs only when a large number of
APs is distributed in close-vicinity, and hence, the adopting statistical DL
CSI for signal decoding at users may hinder the system performance. To
circumvent this, DL pilots can be beamformed to estimate DL channels at the
users, and this approach ensures that the DL pilot sequence length does not
scale with the number of APs. Next, we investigate the impact of DL pilots for
the proposed cell-free massive MIMO with underlay spectrum sharing.

###  V-A DL channel estimation

We define the effective DL desired and interference channel coefficients of
the primary system based on the signal received at $U_{P}(k)$ in (16) as
follows:

$\displaystyle\mu_{P_{ki}}\triangleq\sum\nolimits_{m=1}^{M}\delta_{P_{mi}}\eta_{P_{mi}}^{1/2}f_{mk}\hat{f}_{mi}^{*}\qquad\text{and}\qquad\lambda_{P_{kj}}\triangleq\sum\nolimits_{n=1}^{N}\delta_{S_{nj}}\eta_{S_{nj}}^{1/2}u_{nk}\hat{g}_{nj}^{*}.$
(46)

To estimate DL channels, we need an additional $\tau_{p,d}$ symbol duration
for transmitting DL pilots towards PUs/SUs. We consider the same pilot sharing
technique that was used for UL pilot transmission. The same pilot sequences
will be used between P-APs/S-APs and PUs/SUs. For the sake of exposition, we
denote the primary and secondary DL pilot sequences by
$\boldsymbol{\phi}_{P_{k},d}$, respectively, for $k\in\\{1,\cdots,K\\}$ and
$\boldsymbol{\phi}_{S_{l},d}$ for $l\in\\{1,\cdots,L\\}$. The pilot signal
sent by the $m$th P-AP can be written as

$\displaystyle\mathbf{x}_{P_{m},d}$ $\displaystyle=$
$\displaystyle\sqrt{P_{p,d}}\sum\nolimits_{i=1}^{K}\delta_{P_{mi}}\eta_{P_{mi}}^{1/2}\hat{f}_{mi}^{*}\boldsymbol{\phi}_{P_{i},d},$
(47)

where $P_{p,d}\triangleq\tau_{p,d}P_{d}$ and $P_{d}$ is the average DL pilot
transmit power at each P-AP/S-AP. Moreover, $\hat{f}_{mi}$ is defined in (6a).
Next, the pilot vector received at $U_{P}(k)$ can be written as

$\displaystyle\mathbf{y}^{\prime}_{P_{k},d}$ $\displaystyle=$
$\displaystyle\sum\nolimits_{m=1}^{M}f_{mk}\mathbf{x}_{P_{m},d}+\sum\nolimits_{n=1}^{N}u_{nk}\mathbf{x}_{S_{n},d}+\mathbf{n}^{\prime}_{P_{k},d},$
(48)

where $\mathbf{n}^{\prime}_{P_{k},d}$ is the AWGN vector with i.i.d.
$\mathcal{CN}(0,1)$ elements, at $U_{P}(k)$. Then, we rewrite the received
pilot vector at $U_{P}(k)$ by substituting (47) into (48) as

$\displaystyle\mathbf{y}^{\prime}_{P_{k},d}$ $\displaystyle=$
$\displaystyle\sqrt{P_{p,d}}\left(\sum\nolimits_{i=1}^{K}\mu_{P_{ki}}\boldsymbol{\phi}_{P_{i},d}+\sum\nolimits_{j=1}^{L}\lambda_{P_{kj}}\boldsymbol{\phi}_{S_{j},d}\right)+\mathbf{n}^{\prime}_{P_{k},d},$
(49)

where $\mu_{P_{ki}}$ and $\lambda_{P_{kj}}$ are the effective desired and
interference DL channels (46). To estimate the effective primary desired DL
channel, the sufficient statistics can be obtained by projecting
$\boldsymbol{\phi}_{P_{k},d}^{H}$ into (49) as

$\displaystyle y_{P_{k},d}$ $\displaystyle=$
$\displaystyle\boldsymbol{\phi}_{P_{k},d}^{H}\mathbf{y}^{\prime}_{P_{k},d}=\sqrt{P_{p,d}}\left(\mu_{P_{kk}}+\lambda_{P_{kk}}\right)+n_{P_{k},d},$
(50)

where
$n_{P_{k},d}=\boldsymbol{\phi}_{P_{k},d}^{H}\mathbf{n}^{\prime}_{P_{k},d}\sim\mathcal{CN}(0,1)$
is the AWGN at $U_{P}(k)$.

Proposition 3: With beamformed DL pilots, the MMSE estimate of $\mu_{P_{kk}}$
(46) is given by

$\displaystyle\\!\\!\\!\\!\\!\hat{\mu}_{P_{kk}}\\!$ $\displaystyle=$
$\displaystyle\\!\mathbb{E}\\!\left[{\mu_{P_{kk}}}\right]\\!+\\!\frac{\mathbb{C}\mathrm{ov}\\!\left({\mu_{P_{kk}}y_{P_{k},d}^{*}}\right)}{\mathbb{C}\mathrm{ov}\\!\left({y_{P_{k},d}y_{P_{k},d}^{*}}\right)}\left(y_{P_{k},d}-\mathbb{E}\\!\left[{y_{P_{k},d}}\right]\right)$
(51) $\displaystyle=$
$\displaystyle\mathbb{E}\\!\left[{\mu_{P_{kk}}}\right]+\frac{\sqrt{P_{p,d}}\mathbb{V}\mathrm{ar}\\!\left[{\mu_{P_{kk}}}\right]}{{P_{p,d}}\left(\mathbb{V}\mathrm{ar}\\!\left[{\mu_{P_{kk}}}\right]+\mathbb{V}\mathrm{ar}\\!\left[{\lambda_{P_{kk}}}\right]\right)+1}\times\left(y_{P_{k},d}-\sqrt{P_{p,d}}\left(\mathbb{E}\\!\left[{\mu_{P_{kk}}}\right]+\mathbb{E}\\!\left[{\lambda_{P_{kk}}}\right]\right)\right).$

Then, the MMSE estimate of ${\mu}_{P_{kk}}$ is given by evaluating (51) as

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\hat{\mu}_{P_{kk}}\\!$
$\displaystyle=$
$\displaystyle\\!\frac{\\!\sqrt{\\!P_{p,d}}v_{P_{kk}}y_{P_{k},d}\\!+\\!\\!\sum\nolimits_{m=1}^{M}\\!\delta_{P_{mk}}\eta_{P_{mk}}^{1/2}\rho_{f_{mk}}\\!\\!+\\!P_{p,d}\\!\left(\\!\\!u_{P_{kk}}\\!\\!\sum\nolimits_{m=1}^{M}\\!\delta_{P_{mk}}\eta_{P_{mk}}^{1/2}\rho_{f_{mk}}\\!\\!-\\!v_{P_{kk}}\\!\sum\nolimits_{n=1}^{N}\\!\delta_{S_{nk}}\eta_{S_{nk}}^{1/2}\rho_{u_{nk}}\\!\right)\\!}{P_{p,d}\left(v_{P_{kk}}+u_{P_{kk}}\right)+1}.$
(52)

In (52), $v_{P_{kk}}$ and $u_{P_{kk}}$ are defined as

$\displaystyle
v_{P_{kk}}\triangleq\sum\nolimits_{m=1}^{M}\delta_{P_{mk}}\eta_{P_{mk}}\zeta_{{f}_{mk}}\rho_{f_{mk}}\qquad\text{and}\qquad
u_{P_{kk}}\triangleq\sum\nolimits_{n=1}^{N}\delta_{S_{nk}}\eta_{S_{nk}}\zeta_{{u}_{nk}}\rho_{u_{nk}}.$
(53)

The actual effective/desired DL channel coefficient is given by
$\mu_{P_{kk}}=\hat{\mu}_{P_{kk}}+\epsilon^{\mu}_{P_{kk}}$, where
$\epsilon^{\mu}_{P_{kk}}$ is estimation error, which is independent of
respective channel estimate.

###### Proof.

Appendix A-C. ∎

Remark 6: It is worth noting that the parameters $\mu_{S_{lj}}$,
$\lambda_{S_{li}}$, $\hat{\mu}_{S_{ll}}$, $v_{S_{ll}}$, and $u_{S_{ll}}$
corresponding to the secondary system can be obtained by replacing subscripts
$\\{P,m,M,k,i\\}$ in (46), (52), and (53), respectively, by $\\{S,n,N,l,j\\}$.

###  V-B Primary achievable DL rate with DL pilots

We can rewrite the received signal at $U_{P}(k)$ in (16) via the effective
desired and interference DL channel coefficients as follows:

$\displaystyle{r}_{P_{k},d}$ $\displaystyle=$
$\displaystyle\sqrt{P_{P}}\mu_{P_{kk}}q_{P_{k}}+\sqrt{P_{P}}\sum\nolimits_{i\neq
k}^{K}\mu_{P_{ki}}q_{P_{i}}+\sqrt{P_{S}}\sum\nolimits_{j=1}^{L}\lambda_{P_{kj}}q_{S_{j}}+n_{P_{k},d},$
(54)

where $n_{P_{k},d}\sim\mathcal{CN}(0,1)$ is the AWGN at $U_{P}(k)$. Then, by
using the DL channel estimate at $U_{P}(k)$, (54) can be rearranged to
facilitate decoding the desired signal as

$\displaystyle{r}_{P_{k},d}$ $\displaystyle=$
$\displaystyle\underbrace{\sqrt{P_{P}}\mathbb{E}\\!\left[{\left(\mu_{P_{kk}}|\hat{\mu}_{P_{kk}}\right)}\right]q_{P_{k}}}_{\text{Desired
signal}}+\underbrace{\sqrt{P_{P}}\left(\left(\mu_{P_{kk}}|\hat{\mu}_{P_{kk}}\right)q_{P_{k}}-\mathbb{E}\\!\left[{\left(\mu_{P_{kk}}|\hat{\mu}_{P_{kk}}\right)}\right]q_{P_{k}}\right)}_{\text{Detection
uncertainty}}$ (55)
$\displaystyle+\underbrace{\sqrt{P_{P}}\sum\nolimits_{i\neq
k}^{K}\left(\mu_{P_{ki}}|\hat{\mu}_{P_{kk}}\right)q_{P_{i}}}_{\text{Inter-
system
interference}}+\underbrace{\sqrt{P_{S}}\sum\nolimits_{j=1}^{L}\left(\lambda_{P_{kj}}|\hat{\mu}_{P_{kk}}\right)q_{S_{j}}}_{\text{Intra-
system interference}}+\underbrace{n_{P_{k},d}}_{\text{AWGN}}.$

By using (55), the SINR at $U_{P}(k)$ can be given as

$\displaystyle\gamma_{P_{k},d}=\frac{P_{P}|\mathbb{E}\\!\left[{\mu_{P_{kk}}|\hat{\mu}_{P_{kk}}}\right]|^{2}}{P_{P}\sum\nolimits_{i=1}^{K}\mathbb{E}\\!\left[{|\mu_{P_{ki}}|^{2}|\hat{\mu}_{P_{kk}}}\right]-P_{P}|\mathbb{E}\\!\left[{\mu_{P_{kk}}|\hat{\mu}_{P_{kk}}}\right]|^{2}+P_{S}\sum\nolimits_{j=1}^{L}\mathbb{E}\\!\left[{|\lambda_{P_{kj}}|^{2}|\hat{\mu}_{P_{kk}}}\right]+1}.$
(56)

From the facts that (i) $\mu_{P_{kk}}$ is Gaussian distributed, (ii)
$\hat{\mu}_{P_{kk}}$ and $\epsilon^{\mu}_{P_{kk}}$ are independent, and (iii)
$\mu_{P_{kk}}$, $\mu_{P_{ki}}$, and $\lambda_{P_{kj}}$ are independent for
$i,j\neq k$ [42], we rewrite the SINR in (56) as follows:

$\displaystyle\gamma_{P_{k},d}$ $\displaystyle=$
$\displaystyle\frac{P_{P}|\hat{\mu}_{P_{kk}}|^{2}}{P_{P}\left(\sum\nolimits_{i\neq
k}^{K}\mathbb{E}\\!\left[{|\mu_{P_{ki}}|^{2}}\right]+\mathbb{E}\\!\left[{|\epsilon^{\mu}_{P_{kk}}|^{2}}\right]\right)+P_{S}\sum\nolimits_{j=1}^{L}\mathbb{E}\\!\left[{|\lambda_{P_{kj}}|^{2}}\right]+1}.$
(57)

Thereby, we define the achievable DL rate as

$\displaystyle
R_{P_{k},d}=\left({\tau_{d,d}}/{\tau_{c}}\right)\mathbb{E}_{\hat{\mu}_{P_{kk}}}\\!\left[{\mathrm{log}_{2}\left(1+\gamma_{P_{k},d}\right)}\right],$
(58)

where
$\tau_{d,d}\\!=\\!\tau_{c}\\!-\\!\left(\\!\tau_{p}\\!+\\!\tau_{p,d}\\!\right)$.
From Jensen’s inequality, an upper bound for DL rate at $U_{P}(k)$ is derived
as

$\displaystyle
R^{ub}_{P_{k},d}=\left({\tau_{d,d}}/{\tau_{c}}\right)\mathrm{log}_{2}\left(1+\mathbb{E}_{\hat{\mu}_{P_{kk}}}\\!\left[{\gamma_{P_{k},d}}\right]\right),$
(59)

where $\mathbb{E}_{\hat{\mu}_{kk}}\\!\left[{\gamma_{P_{k},d}}\right]$ is given
by

$\displaystyle\mathbb{E}_{\hat{\mu}_{P_{kk}}}\\!\left[{\gamma_{P_{k},d}}\right]=\frac{P_{P}\mathbb{E}\\!\left[{|\hat{\mu}_{P_{kk}}|^{2}}\right]}{P_{P}\left(\sum\nolimits_{i\neq
k}^{K}\mathbb{E}\\!\left[{|\mu_{P_{ki}}|^{2}}\right]+\mathbb{E}\\!\left[{|\epsilon^{\mu}_{P_{kk}}|^{2}}\right]\right)+P_{S}\sum\nolimits_{j=1}^{L}\mathbb{E}\\!\left[{|\lambda_{P_{kj}}|^{2}}\right]+1}.$
(60)

Proposition 4: By evaluating expectation terms in (60), the effective SINR at
$U_{P}(k)$ in the case of beamformed DL pilots by the P-APs can be derived as
follows:

$\displaystyle\\!\\!\\!\\!\mathbb{E}_{\hat{\mu}_{P_{kk}}}\\!\left[{\gamma_{P_{k},d}}\right]\\!=\\!\frac{P_{P}\left(\sum\nolimits_{m=1}^{M}\sum\nolimits_{m^{\prime}=1}^{M}\delta_{P_{mk}}\eta_{P_{mk}}^{1/2}\delta_{P_{m^{\prime}k}}\eta_{P_{m^{\prime}k}}^{1/2}\rho_{f_{mk}}\rho_{f_{m^{\prime}k}}+v_{P_{kk}}-\kappa^{\epsilon}_{P_{kk}}\right)}{\\!P_{P}\\!\left(\sum\nolimits_{i\neq
k}^{K}v_{P_{ki}}\\!+\\!\kappa^{\epsilon}_{P_{kk}}\right)\\!+\\!P_{S}\\!\left(\sum\nolimits_{n=1}^{N}\sum\nolimits_{j=1}^{L}\delta_{S_{nj}}\eta_{S_{nj}}\rho_{g_{nj}}\zeta_{{u}_{nk}}\\!+\\!\sum\nolimits_{n=1}^{N}\delta_{S_{nk}}\eta_{S_{nk}}\rho_{u_{nk}}^{2}\right)\\!+\\!1},$
(61)

where $\kappa^{\epsilon}_{P_{kk}}$ and $v_{P_{ki}}$ are defined as

$\displaystyle\\!\\!\\!\\!\\!\kappa^{\epsilon}_{P_{kk}}\\!\triangleq\\!\frac{\left(1+P_{p,d}u_{P_{kk}}\right)^{2}v_{P_{kk}}+\left(P_{p,d}v_{P_{kk}}\right)^{2}u_{P_{kk}}+P_{p,d}v_{P_{kk}}^{2}}{\left(P_{p,d}\left(v_{P_{kk}}+u_{P_{kk}}\right)+1\right)^{2}}\quad\text{and}\quad
v_{P_{ki}}\\!\triangleq\\!\sum\nolimits_{m=1}^{M}\delta_{P_{mi}}\eta_{P_{mi}}\rho_{f_{mi}}\zeta_{{f}_{mk}},$
(62)

where $P_{p,d}\triangleq\tau_{p,d}P_{d}$.

###### Proof.

Appendix C. ∎

###  V-C Secondary achievable DL rate with DL pilots

We define an upper bound of the DL achievable rate at $U_{S}(l)$ via steps
similar to (54)-(59) as

$\displaystyle
R^{ub}_{S_{l},d}=\left({\tau_{d,d}}/{\tau_{c}}\right)\mathrm{log}_{2}\left(1+\mathbb{E}_{\hat{\mu}_{S_{ll}}}\\!\left[{\gamma_{S_{l},d}}\right]\right),$
(63)

where $\mathbb{E}_{\hat{\mu}_{S_{ll}}}\\!\left[{\gamma_{S_{l},d}}\right]$ is
given by

$\displaystyle\mathbb{E}_{\hat{\mu}_{S_{ll}}}\\!\left[{\gamma_{S_{l},d}}\right]$
$\displaystyle=$
$\displaystyle\frac{P_{S}\mathbb{E}\\!\left[{|\hat{\mu}_{S_{ll}}|^{2}}\right]}{P_{S}\left(\sum\nolimits_{j\neq
l}^{L}\mathbb{E}\\!\left[{|\mu_{S_{lj}}|^{2}}\right]+\mathbb{E}\\!\left[{|\epsilon^{\mu}_{S_{ll}}|^{2}}\right]\right)+P_{P}\sum\nolimits_{i=1}^{K}\mathbb{E}\\!\left[{|\lambda_{S_{li}}|^{2}}\right]+1}.$
(64)

Then, by following steps similar to (61) to evaluate the expectation terms in
(64), the effective SINR at $U_{S}(l)$ can be written as

$\displaystyle\\!\\!\\!\\!\\!\\!\
\\!\\!\\!\\!\\!\mathbb{E}_{\hat{\mu}_{S_{ll}}}\\!\left[{\gamma_{S_{l},d}}\right]\\!=\\!\frac{P_{S}\left(\sum\nolimits_{n=1}^{N}\sum\nolimits_{n^{\prime}=1}^{N}\delta_{S_{nl}}\eta_{S_{nl}}^{1/2}\delta_{S_{n^{\prime}l}}\eta_{S_{n^{\prime}l}}^{1/2}\rho_{g_{nl}}\rho_{g_{n^{\prime}l}}+v_{S_{ll}}-\kappa^{\epsilon}_{S_{ll}}\right)}{P_{S}\left(\sum\nolimits_{j\neq
l}^{L}v_{S_{lj}}\\!+\\!\kappa^{\epsilon}_{S_{ll}}\right)\\!+\\!P_{P}\left(\sum\nolimits_{m=1}^{M}\sum\nolimits_{i=1}^{K}\delta_{P_{mi}}\eta_{P_{mi}}\rho_{f_{mi}}\zeta_{{v}_{ml}}\\!+\\!\sum\nolimits_{m=1}^{M}\delta_{P_{ml}}\eta_{P_{ml}}\rho_{v_{ml}}^{2}\right)\\!+\\!1},$
(65)

where $\kappa^{\epsilon}_{S_{ll}}$ and $v_{S_{lj}}$ are defined as

$\displaystyle\kappa^{\epsilon}_{S_{ll}}\triangleq\frac{\left(1+P_{p,d}u_{S_{ll}}\right)^{2}v_{S_{ll}}+\left(P_{p,d}v_{S_{ll}}\right)^{2}u_{S_{ll}}+P_{p,d}v_{S_{ll}}^{2}}{\left(P_{p,d}\left(v_{S_{ll}}+u_{S_{ll}}\right)+1\right)^{2}}\quad\text{and}\quad
v_{S_{lj}}\triangleq\sum\nolimits_{n=1}^{N}\delta_{S_{nj}}\eta_{S_{nj}}\rho_{g_{nj}}\zeta_{{g}_{nl}}.$
(66)

Remark 7: We reveal through our numerical results in Section VIII that the
adoption of estimated DL CSI at the PUs/SUs for signal decoding can be
exploited to boost the achievable rates by circumventing the less prevalent
channel hardening property in cell-free massive MIMO compared to that of co-
located counterpart. The underlying implication is that the assumption of
statistical DL channels are approximately equal to the instantaneous
counterparts may not be accurate for cell-free massive MIMO.

##  VI Achievable Rate analysis of NOMA-Aided Underlay Spectrum Sharing

###  VI-A Primary achievable rate with NOMA

We rearrange the received signal at $U_{P}(a,k)$ (II-E) to decode the desired
signal as

$\displaystyle{r}_{P_{ak}}$ $\displaystyle=$
$\displaystyle\underbrace{\sqrt{P_{P}}\mathbb{E}\\!\left[{\sum\nolimits_{m=1}^{M}\eta_{P_{mak}}^{1/2}f_{mak}\hat{f}_{mak}^{*}}\right]q_{P_{ak}}}_{\text{Desired
signal}}$ (67)
$\displaystyle+\underbrace{\sqrt{P_{P}}\sum\nolimits_{m=1}^{M}\eta_{P_{mak}}^{1/2}f_{mak}\hat{f}_{mak}^{*}q_{P_{ak}}-\sqrt{P_{P}}\mathbb{E}\\!\left[{\sum\nolimits_{m=1}^{M}\eta_{P_{mak}}^{1/2}f_{mak}\hat{f}_{mak}^{*}}\right]q_{P_{ak}}}_{\text{Detection
uncertinity}}+w_{{P_{ak}}},$

where $w_{{P_{ak}}}$ denotes the effective noise at $U_{P}(a,k)$ containing
intra-cluster interference after SIC, error propagation due to imperfect SIC,
intra-system interference, inter-system interference, and AWGN given in
(II-E). From (67), we write the SINR at $U_{P}(a,k)$ as

$\displaystyle\gamma_{P_{ak}}=\frac{{P_{P}}\left|\mathbb{E}\\!\left[{\sum\nolimits_{m=1}^{M}{\eta_{P_{mak}}^{1/2}f_{mak}\hat{f}_{mak}^{*}}}\right]\right|^{2}}{P_{P}\mathbb{V}\mathrm{ar}\\!\left[{\sum\nolimits_{m=1}^{N}{\eta_{P_{mak}}^{1/2}f_{mak}\hat{f}_{mak}^{*}}}\right]+{P_{P}}\sum\nolimits_{i=1}^{4}{I_{P_{ai}}}+\mathbb{E}\\!\left[{\left|n_{P_{ak}}\right|^{2}}\right]},$
(68)

where ${I_{P_{ai}}}$ for $i\in\\{1,2,3,4\\}$ can be defined as

$\displaystyle{I_{P_{a1}}}$ $\displaystyle=$
$\displaystyle\mathbb{E}\\!\left[{\\!\left|\sum\nolimits_{m=1}^{M}\sum\nolimits_{i=1}^{k-1}\eta_{P_{mai}}^{1/2}f_{mak}\hat{f}_{mak}^{*}\right|^{2}}\right],{}$
(69a) $\displaystyle{I_{P_{a2}}}$ $\displaystyle=$
$\displaystyle\mathbb{E}\\!\left[{\left|\sum\nolimits_{m=1}^{M}\sum\nolimits_{i=k+1}^{K}\left(f_{mak}\hat{f}_{mai}^{*}q_{P_{ai}}-\mathbb{E}\\!\left[{f_{mak}\hat{f}_{mai}^{*}}\right]\hat{q}_{ai}\right)\right|^{2}}\right],{}$
(69b) $\displaystyle{I_{P_{a3}}}$ $\displaystyle=$
$\displaystyle\mathbb{E}\\!\left[{\left|\sum\nolimits_{m=1}^{M}\sum\nolimits_{a^{\prime}\neq
a}^{A}\sum\nolimits_{i=1}^{K}\eta_{P_{ma^{\prime}i}}^{1/2}f_{mak}\hat{f}_{ma^{\prime}i}^{*}\right|^{2}}\right],$
$\displaystyle{I_{P_{a4}}}$ $\displaystyle=$
$\displaystyle\frac{{P_{S}}}{{P_{P}}}\mathbb{E}\\!\left[{\left|\sum\nolimits_{n=1}^{N}\sum\nolimits_{b=1}^{B}\sum\nolimits_{j=1}^{L}\eta_{S_{nbj}}^{1/2}u_{nak}\hat{g}_{nbj}^{*}\right|^{2}}\right].{}$
(69c)

Then, we compute the SINR by evaluating the expectation and variance terms in
(68) as

$\displaystyle\\!\\!\\!\\!\\!\\!\gamma_{P_{ak}}\\!=\\!\frac{P_{P}\left(\sum\nolimits_{m=1}^{M}\eta_{P_{mak}}^{1/2}\alpha_{f_{mak}}\right)^{2}}{P_{P}\\!\sum\nolimits_{m=1}^{M}\sum\nolimits_{a=1}^{A}\sum\nolimits_{i=1}^{K}\eta_{P_{mai}}\alpha_{f_{mai}}\zeta_{f_{mak}}\\!+\\!{P_{P}}\\!\sum\nolimits_{i=1}^{3}I^{\prime}_{P_{ai}}\\!+\\!\sum\nolimits_{n=1}^{N}\sum\nolimits_{b=1}^{B}\sum\nolimits_{j=1}^{L}\eta_{S_{nbj}}\alpha_{g_{nbj}}\zeta_{{u}_{nak}}\\!+\\!1},$
(70)

where $I^{\prime}_{P_{ai}}$ for $i\in\\{1,2,3\\}$ is defined as

$\displaystyle{I^{\prime}_{P_{a1}}}$ $\displaystyle=$
$\displaystyle\sum\nolimits_{m=1}^{M}\sum\nolimits_{i=1}^{k-1}\eta_{P_{mai}}\alpha_{f_{mak}}^{2}\left(\frac{\zeta_{{f}_{mai}}}{\zeta_{{f}_{mak}}}\right)^{2},{}$
(71a) $\displaystyle{I^{\prime}_{P_{a2}}}$ $\displaystyle=$ $\displaystyle
2\sum\nolimits_{m=1}^{M}\sum\nolimits_{i=k+1}^{K}(1-\vartheta_{P_{ai}})\eta_{P_{mai}}\alpha_{f_{mak}}^{2}\left(\frac{\zeta_{{f}_{mai}}}{\zeta_{{f}_{mak}}}\right)^{2},{}$
(71b) $\displaystyle{I^{\prime}_{P_{a3}}}$ $\displaystyle=$
$\displaystyle\frac{{P_{S}}}{{P_{P}}}\sum\nolimits_{n=1}^{N}\sum\nolimits_{j=1}^{L}\eta_{S_{naj}}\alpha_{f_{mak}}^{2}\left(\frac{\zeta_{{u}_{nak}}\zeta_{{g}_{naj}}}{\zeta_{{f}_{mak}}}\right)^{2}.{}$
(71c)

The derivation of (70) follows steps similar to those in Appendix B, and
hence, it is omitted for the sake of brevity. Next, the achievable rate of
$U_{P}(a,k)$ and the sum rate of primary system are given by

$\displaystyle
R_{P_{ak}}=\left({\left(\tau_{c}-\tau_{p}\right)}/{\tau_{c}}\right)\mathrm{log}_{2}\left(1+\gamma_{P_{ak}}\right)\quad\text{and}\quad
R_{P}=\sum\nolimits_{a=1}^{A}\sum\nolimits_{k=1}^{K}R_{P_{ak}},$ (72)

where $\gamma_{P_{ak}}$ is defined in (70).

###  VI-B Secondary achievable rate with NOMA

We follow a similar analysis to Section VI-A for deriving the sum rate of the
secondary system as

$\displaystyle
R_{S}=\sum\nolimits_{b=1}^{B}\sum\nolimits_{l=1}^{L}R_{S_{bl}},$ (73)

where $R_{S_{bl}}$ is the achievable rate of $U_{S}(b,l)$ and given by

$\displaystyle
R_{S_{bl}}=\left({\left(\tau_{c}-\tau_{p}\right)}/{\tau_{c}}\right)\mathrm{log}_{2}\left(1+\gamma_{S_{bl}}\right).$
(74)

In (74), we obtain the SINR at $U_{S}(b,l)$ denoted by $\gamma_{S_{bl}}$ by
following (70) and replacing the primary system variables with respective
secondary system variables in (70) as

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\gamma_{S_{bl}}\\!=\\!\frac{P_{S}\left(\sum\nolimits_{n=1}^{N}\eta_{S_{nbl}}^{1/2}\alpha_{g_{nbl}}\right)^{2}}{P_{S}\sum\nolimits_{n=1}^{N}\sum\nolimits_{b=1}^{B}\sum\nolimits_{j=1}^{L}\eta_{S_{nbj}}\alpha_{g_{nbj}}\zeta_{g_{nbl}}\\!+\\!P_{S}\sum\nolimits_{j=1}^{3}I^{\prime}_{S_{bj}}\\!+\\!\sum\nolimits_{m=1}^{M}\sum\nolimits_{a=1}^{A}\sum\nolimits_{j=1}^{K}\eta_{P_{mai}}\alpha_{f_{mai}}\zeta_{{v}_{mbl}}\\!+\\!1},$
(75)

where $I^{\prime}_{S_{bi}}$ for $i\in\\{1,2,3\\}$ is given by

$\displaystyle{I^{\prime}_{S_{b1}}}$ $\displaystyle=$
$\displaystyle\sum\nolimits_{n=1}^{N}\sum\nolimits_{j=1}^{l-1}\eta_{S_{nbj}}\alpha_{g_{nbl}}^{2}\left(\frac{\zeta_{{g}_{nbj}}}{\zeta_{{g}_{nbl}}}\right)^{2},{}$
(76a) $\displaystyle{I^{\prime}_{S_{b2}}}$ $\displaystyle=$ $\displaystyle
2\sum\nolimits_{n=1}^{N}\sum\nolimits_{j=l+1}^{L}(1-\vartheta_{S_{bj}})\eta_{S_{nbj}}\alpha_{g_{nbl}}^{2}\left(\frac{\zeta_{{g}_{nbj}}}{\zeta_{{g}_{nbl}}}\right)^{2},{}$
(76b) $\displaystyle{I^{\prime}_{S_{b3}}}$ $\displaystyle=$
$\displaystyle\frac{{P_{P}}}{P_{S}}\sum\nolimits_{m=1}^{M}\sum\nolimits_{i=1}^{K}\eta_{P_{mbi}}\alpha_{g_{nbl}}^{2}\left(\frac{\zeta_{{v}_{mbl}}\zeta_{{f}_{mbi}}}{\zeta_{{g}_{nbl}}}\right)^{2}.{}$
(76c)

##  VII The implication of DL pilot transmission with NOMA

The effective DL desired and interference channel coefficients at $U_{P}(a,k)$
are defined from (II-E) as

$\displaystyle\mu_{P_{a^{\prime}i}^{ak}}\triangleq\sum\nolimits_{m=1}^{M}\eta_{P_{ma^{\prime}i}}^{1/2}f_{mak}\hat{f}_{ma^{\prime}i}^{*}\qquad\text{and}\qquad\lambda_{P_{bj}^{ak}}\triangleq\sum\nolimits_{n=1}^{N}\eta_{S_{nbj}}^{1/2}u_{nak}\hat{g}_{nbj}^{*}.$
(77)

The P-APs/S-APs again use the same pilot sequences in (9) to beamform DL
pilots toward PUs/SUs. Then, a sufficient statistic to estimate the desired
effective DL channel at $U_{P}(a,k)$ is given by

$\displaystyle
y_{P_{ak},d}=\sqrt{P_{p,d}}\left(\sum\nolimits_{i=1}^{K}\mu_{P_{ai}^{ak}}+\sum\nolimits_{j=1}^{L}\lambda_{P_{aj}^{ak}}\right)+n_{P_{ak},d},$
(78)

where $n_{P_{ak},d}\sim\mathcal{CN}(0,1)$ is the AWGN at $U_{P}(a,k)$.
Thereby, the MMSE estimate of $\mu_{P_{at}^{ak}}$ can be derived as [42, 35]

$\displaystyle\\!\\!\\!\\!\\!\hat{\mu}_{P_{at}^{ak}}$ $\displaystyle=$
$\displaystyle\mathbb{E}\\!\left[{\mu_{P_{at}^{ak}}}\right]\\!+\\!\frac{\mathbb{C}\mathrm{ov}\\!\left({\mu_{P_{at}^{ak}}y_{P_{ak},d}^{*}}\right)}{\mathbb{C}\mathrm{ov}\\!\left({y_{P_{ak},d}y_{P_{ak},d}^{*}}\right)}\left(y_{P_{ak},d}-\mathbb{E}\\!\left[{y_{P_{ak},d}}\right]\right)$
(79) $\displaystyle=$
$\displaystyle\mathbb{E}\\!\left[{\mu_{P_{at}^{ak}}}\right]+\frac{\sqrt{P_{p,d}}\mathbb{V}\mathrm{ar}\\!\left[{\mu_{P_{at}^{ak}}}\right]}{{P_{p,d}}\left(\sum\nolimits_{i=1}^{K}\mathbb{V}\mathrm{ar}\\!\left[{\mu_{P_{ai}^{ak}}}\right]+\sum\nolimits_{j=1}^{L}\mathbb{V}\mathrm{ar}\\!\left[{\lambda_{P_{aj}^{ak}}}\right]\right)+1}\times\left(y_{P_{k},d}-\mathbb{E}\\!\left[{y_{P_{ak},d}}\right]\right).$

By evaluating (79) via steps similar to those in Appendix A-C, the MMSE
estimate of $\mu_{P_{at}^{ak}}$ can be derived as

$\displaystyle\hat{\mu}_{P_{at}^{ak}}=\theta_{P_{at}^{ak}}+\Omega_{P_{at}^{ak}}\left(y_{P_{k},d}-\sqrt{P_{p,d}}\left(\sum\nolimits_{i=1}^{K}\theta_{P_{ai}^{ak}}+\sum\nolimits_{j=1}^{L}\psi_{P_{aj}^{ak}}\right)\right),$
(80)

where
$\theta_{P_{ai}^{ak}}=\sum\nolimits_{m=1}^{M}\eta_{P_{mai}}^{1/2}\alpha_{f_{mak}}\zeta_{{f}_{mai}}/\zeta_{{f}_{mak}}$,
$\psi_{P_{aj}^{ak}}=\sum\nolimits_{n=1}^{N}\eta_{S_{naj}}^{1/2}\alpha_{u_{nak}}\zeta_{{g}_{nak}}/\zeta_{{g}_{naj}}$,
and

$\displaystyle\\!\\!\\!\\!\Omega_{P_{at}^{ak}}$ $\displaystyle=$
$\displaystyle{\sqrt{P_{p,d}}\left(\varrho_{P_{at}^{ak}}^{\mu}-\theta_{P_{at}^{ak}}^{2}\right)}\Big{/}\left({P_{p,d}\left(\sum\nolimits_{i=1}^{K}\left(\varrho_{P_{ai}^{ak}}^{\mu}-\theta_{P_{ai}^{ak}}^{2}\right)+\sum\nolimits_{j=1}^{L}\left(\varrho_{P_{bj}^{ak}}^{\lambda}-\psi_{P_{aj}^{ak}}^{2}\right)\right)+1}\right),$
(81a) $\displaystyle\\!\\!\\!\\!\varrho_{P_{ai}^{ak}}^{\mu}$ $\displaystyle=$
$\displaystyle\sum\nolimits_{m=1}^{M}\eta_{P_{mai}}\alpha_{f_{mai}}\left(\zeta_{{f}_{mak}}+\alpha_{f_{mak}}\right),\quad\text{and}\quad\varrho_{P_{bj}^{ak}}^{\lambda}=\sum\nolimits_{n=1}^{N}\eta_{S_{nbj}}\left(\zeta_{{u}_{nak}}\alpha_{g_{nbj}}+\right).$
(81b)

The actual effective DL channel gain is given by
${\mu}_{P_{at}^{ak}}=\hat{\mu}_{P_{at}^{ak}}+\epsilon_{P_{at}^{ak}}^{\mu}$,
where $\epsilon_{P_{at}^{ak}}^{\mu}$ is an estimation error, which is
independent of the channel estimate $\hat{\mu}_{P_{at}^{ak}}$.

###  VII-A Primary achievable DL rate with DL pilots in NOMA

It is assumed that the users in same cluster are ordered based on the
effective channel gains as per (20), and hence, SIC can be adopted to decode
the power-domain NOMA signals at the users [34]. Thereby, the post-processed
signal after an imperfect SIC operation at $U_{P}(a,k)$ can be written as

$\displaystyle{r}_{P_{ak},d}$ $\displaystyle=$
$\displaystyle\underbrace{\sqrt{P_{P}}\mu_{P_{ak}^{ak}}q_{P_{ak}}}_{\text{Desired
signal}}+\underbrace{\sqrt{P_{P}}\sum\nolimits_{i=1}^{k-1}\mu_{P_{ai}^{ak}}q_{P_{ai}}}_{\text{Intra-
cluster interference after
SIC}}+\underbrace{\sqrt{P_{P}}\sum\nolimits_{i=k+1}^{K}\epsilon_{P_{ai}^{ak}}^{\mu}q_{P_{ai}}}_{\text{Error
propagation due to imperfect SIC}}$ (82)
$\displaystyle+\underbrace{\sqrt{P_{P}}\sum\nolimits_{a^{\prime}\neq
a}^{A}\sum\nolimits_{i=1}^{K}\mu_{P_{a^{\prime}i}^{ak}}q_{P_{a^{\prime}i}}}_{\text{Intra-
system
interference}}+\underbrace{\sqrt{P_{S}}\sum\nolimits_{b=1}^{B}\sum\nolimits_{j=1}^{L}\lambda_{P_{bj}^{ak}}q_{S_{bj}}}_{\text{Inter-
system interference}}+\underbrace{n_{P_{ak}}}_{\text{AWGN}}.$

By using (82), the corresponding SINR at $U_{P}(a,k)$ can be derived as

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\gamma_{P_{ak},d}\\!=\\!\frac{P_{P}\left|\mathbb{E}\\!\left[{\mu_{P_{ak}^{ak}}|\hat{\mu}_{P_{ak}^{ak}}}\right]\right|^{2}}{\\!P_{P}\\!\\!\sum\nolimits_{a^{\prime}=1}^{A}\\!\sum\nolimits_{i=1}^{K}\\!\mathbb{E}\\!\left[{\left|\mu_{P_{a^{\prime}i}^{ak}}|\hat{\mu}_{P_{kk}}\right|^{\\!2}}\right]\\!\\!-\\!P_{P}\\!\\!\sum\nolimits_{i=k}^{K}\\!\left|\mathbb{E}\\!\left[{\mu_{P_{ak}^{ak}}|\hat{\mu}_{P_{ak}^{ak}}}\right]\\!\right|^{\\!2}\\!\\!+\\!P_{S}\\!\\!\sum\nolimits_{b=1}^{B}\\!\sum\nolimits_{j=1}^{L}\\!\mathbb{E}\\!\left[{\left|\lambda_{P_{aj}^{ak}}|\hat{\mu}_{P_{ak}^{ak}}\right|^{\\!2}}\right]\\!\\!+\\!1\\!}\\!.$
(83)

By using techniques similar to those used in (57), the expected value of
$\gamma_{P_{ak},d}$ can be written as

$\displaystyle\mathbb{E}_{\hat{\mu}_{P_{ak}^{ak}}}\\!\left[{\gamma_{P_{ak},d}}\right]$
$\displaystyle=$
$\displaystyle\frac{P_{P}\mathbb{E}\\!\left[{\left|\hat{\mu}_{P_{ak}^{ak}}\right|^{2}}\right]}{P_{P}\sum\nolimits_{a^{\prime}\neq
a}^{A}\sum\nolimits_{i=k}^{K}\mathbb{E}\\!\left[{\left|\mu_{P_{ai}^{ak}}\right|^{2}}\right]+\sum_{j=1}^{2}I^{\hat{\mu}}_{j}+P_{S}\sum\nolimits_{b=1}^{B}\sum\nolimits_{j=1}^{L}\mathbb{E}\\!\left[{\left|\lambda_{P_{bj}^{ak}}\right|^{2}}\right]+1},$
(84)

where $I^{\hat{\mu}}_{j}$ for $j\in\\{1,2\\}$ is given as

$\displaystyle I^{\hat{\mu}}_{1}$ $\displaystyle=$ $\displaystyle
P_{P}\sum\nolimits_{i=k+1}^{K}\mathbb{E}\\!\left[{\left|\hat{\mu}_{P_{ak}^{ak}}\right|^{2}}\right]\quad\text{and}\quad
I^{\hat{\mu}}_{2}=P_{P}\sum\nolimits_{i=1}^{K}\mathbb{E}\\!\left[{\left|\epsilon^{\mu}_{P_{ai}^{ak}}\right|^{2}}\right].$
(85a)

By evaluating the expectation terms in (84), the average SINR at $U_{P}(a,k)$
can be derived as

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\mathbb{E}_{\hat{\mu}_{P_{ak}^{ak}}}\\!\left[{\gamma_{P_{ak},d}}\right]\\!$
$\displaystyle=$
$\displaystyle\\!\frac{P_{P}\Phi_{P_{ak}^{ak}}^{\mu}}{\\!P_{P}\\!\sum\nolimits_{a^{\prime}\neq
a}^{A}\\!\sum\nolimits_{i=1}^{K}\\!\varrho_{P_{ai}^{ak}}^{\mu}\\!+\\!P_{P}\\!\sum\nolimits_{i=k+1}^{K}\\!\Phi_{P_{ak}^{ak}}^{\mu}\\!+\\!P_{P}\\!\sum\nolimits_{i=1}^{K}\\!\\!\left(\\!\varrho_{P_{ai}^{ak}}^{\mu}\\!-\\!\Phi_{P_{ai}^{ak}}^{\mu}\\!\right)\\!\\!+\\!P_{S}\\!\sum\nolimits_{b=1}^{B}\\!\sum\nolimits_{j=1}^{L}\\!\varrho_{P_{bj}^{ak}}^{\lambda}\\!\\!+\\!1\\!},$
(86)

where $\Phi_{P_{ai}^{ak}}^{\mu}$ is given by

$\displaystyle\\!\\!\\!\\!\\!\Phi_{P_{ai}^{ak}}^{\mu}$
$\displaystyle\triangleq$
$\displaystyle\theta_{P_{at}^{ak}}^{2}\\!+\\!\Omega_{P_{at}^{ak}}^{2}\left({P_{p,d}}\left(\sum\nolimits_{i=1}^{K}\varrho_{P_{ai}^{ak}}^{\mu}\\!+\\!\sum\nolimits_{j=1}^{L}\varrho_{P_{bj}^{ak}}^{\lambda}\right)\\!+\\!1\\!-\\!{P_{p,d}}\left(\sum\nolimits_{i=1}^{K}\theta_{P_{ai}^{ak}}\\!+\\!\sum\nolimits_{j=1}^{L}\psi_{P_{aj}^{ak}}\right)^{2}\right).$
(87)

Then, an upper bound for the DL rate at $U_{P}(a,k)$ with estimated DL
channels via the beamformed pilots can be derived as
$R^{ub}_{P_{ak},d}\\!=\\!\left({\tau_{d,d}}/{\tau_{c}}\right)\mathrm{log}_{2}\left(\\!1\\!+\\!\mathbb{E}_{\hat{\mu}_{P_{ak}^{ak}}}\\!\left[{\gamma_{P_{ak},d}}\right]\right)$,
where
$\mathbb{E}_{\hat{\mu}_{P_{ak}^{ak}}}\\!\left[{\gamma_{P_{ak},d}}\right]$ is
defined in (86).

###  VII-B Secondary achievable DL rate with DL pilots in NOMA

By following an analysis similar to Section VII-A, the achievable rate at
$U_{S}(b,l)$ with estimated DL channels can be derived as
$R^{ub}_{S_{bl},d}=\left({\tau_{d,d}}/{\tau_{c}}\right)\mathrm{log}_{2}\left(1+\mathbb{E}_{\hat{\mu}_{P_{bl}^{bl}}}\\!\left[{\gamma_{S_{bl},d}}\right]\right)$,
where
$\mathbb{E}_{\hat{\mu}_{P_{bl}^{bl}}}\\!\left[{\gamma_{S_{bl},d}}\right]$ is
defined as

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\mathbb{E}_{\hat{\mu}_{S_{bl}^{bl}}}\\!\left[{\gamma_{S_{bl},d}}\right]\\!$
$\displaystyle=$
$\displaystyle\\!\frac{P_{S}\Phi_{S_{bl}^{bl}}^{\mu}}{\\!P_{S}\\!\sum\nolimits_{b^{\prime}\neq
b}^{B}\sum\nolimits_{j=1}^{L}\varrho_{S_{bj}^{bl}}^{\mu}\\!+\\!P_{S}\\!\sum\nolimits_{j=l+1}^{L}\Phi_{S_{bl}^{bl}}^{\mu}\\!+\\!P_{S}\\!\sum\nolimits_{j=1}^{L}\\!\left(\varrho_{S_{bj}^{bl}}^{\mu}\\!-\\!\Phi_{S_{bj}^{bl}}^{\mu}\right)\\!+\\!P_{P}\\!\sum\nolimits_{a=1}^{A}\sum\nolimits_{i=1}^{K}\varrho_{S_{ai}^{bl}}^{\lambda}\\!+\\!1},$
(88)

where $\Phi_{S_{bl}^{bl}}^{\mu}$, $\varrho_{S_{bj}^{bl}}^{\mu}$, and
$\varrho_{S_{ai}^{bl}}^{\lambda}$ can be obtained by replacing the primary
subscripts $\\{P,M,m,A,a,K,k,i\\}$ with respective secondary subscripts
$\\{S,N,n,B,b,L,l,j\\}$ in (87), (81b), and (81b), respectively.

##  VIII Numerical Results

In this section, our numerical results are presented to obtained useful
insights. The simulation parameters are as follows: $\tau_{c}=196$,
$\tau_{p}=\tau_{p,d}=\max(K,L)$, and
$\zeta_{{h}_{ab}}=(d_{0}/d_{ab})^{\nu}\times 10^{\varphi_{ab/10}}$, where
$d_{ab}$ is transmission distance between the $a$th P-AP/S-AP and the $b$th
PU/SU, $d_{0}$ is the reference distance, and $\nu$ is the path-loss exponent.
Here, $10^{\varphi_{ab/10}}$ captures the shadow fading with
$\varphi_{ab}\sim\mathcal{N}(0,8)$. In an area of $800\times 800\,m^{2}$, the
P-APs/S-APs are uniformly distributed, while PUs/SUs are randomly placed.

Figure 2: User rate versus the primary transmit power for $K=L=10$, $M=N=32$,
$I_{T}=0$ dB, $\sigma_{P}^{2}=\sigma_{S}^{2}=1$, $\nu=2.4$, and $d_{0}=1$ m.
Moreover, $P_{S}=P_{P}/2$ and
$P_{S_{n}}=\min{\left(P_{S},{I_{T_{1}}}\big{/}{Z_{1}},\cdots,{I_{T_{k}}}\big{/}{Z_{k}},\cdots,{I_{T_{K}}}\big{/}{Z_{K}}\right)}$.

In Fig. 2, the implications of our max-min based multi-objective transmit
power allocation are investigated. To this end, the achievable user rates of
the primary and secondary systems are plotted against the primary transmit
power per P-AP ($P_{P}$). The maximum allowable secondary transmit power
($P_{S}$) at each S-AP is kept at $P_{P}/2$. Then, the primary and secondary
transmit power control coefficients are computed by using the proposed max-min
algorithm in Section IV. The rates with max-min power allocation are compared
with those with uniform power allocation. The pair {$U_{P}(1),U_{S}(1)$} is
the users with strongest channels, while the pair {$U_{P}(10),U_{S}(10)$}
represents the users with weakest channel gains. Fig. 2 reveals that PUs/SUs
experience distinct achievable rates when the uniform power allocation is
adopted. Thus, the achievable rates are dependent on the detrimental near-far
effects. When the proposed max-min power control is employed, all PUs/SUs
achieve their respective common rates regardless of the near-far effects. For
instance, at $P_{P}=0$ dBW, the weaker user $U_{P}(10)$ achieves a rate gain
of $241.3\%$ from our max-min power allocation over the uniform power
allocation.

Figure 3: The primary sum rate versus the primary transmit power for $K=L=10$,
$M=N=32$, $I_{T}=0\,dB$, $\sigma_{P}^{2}=\sigma_{S}^{2}=1$, $\nu=2.4$, and
$d_{0}=1$ m.

In Fig. 3 and Fig. 4, an achievable rate comparison is presented for cell-
free/co-located, statistical/estimated DL CSI with max-min/uniform power
allocations. In this context, the achievable sum rates of the primary and
secondary systems, respectively, are plotted in Fig. 3 and Fig. 4 as a
function of the primary transmit power per P-AP ($P_{P}$). In Fig. 4, $P_{S}$
is set to $P_{P}/2$ at each S-AP. When the uniform power allocation is
adopted, the achievable sum rate of the primary system increases monotonically
with $P_{P}$. However, for the secondary system, the achievable sum rate
gradually increases up to a maximum in the low $P_{P}$ regime since $P_{S}$ is
proportional to $P_{P}$, and then it decreases as $P_{P}$ grows without bound.
The reason for this behavior is that when $P_{P}$ increases, the respective
secondary maximum allowable transmit power also increases since
$P_{S}=P_{P}/2$. Thus, the secondary sum rate grows gradually until the
secondary transmit power constraints in (28) are met. At this point S-APs
transmit signals with their maximum allowable transmit power $P_{S}$, and the
secondary sum rate reaches a maximum. Simultaneously, the primary transmit
powers at the P-APs keep increasing, and this causes a high level of primary
CCI at SUs. Consequently, the secondary sum rate decreases as $P_{P}$ grows
without bound. When max-min transmit power allocation algorithm is adopted for
both systems, the primary and secondary achievable sum rates increase with the
primary transmit power $P_{P}$. Furthermore, Fig. 3 and Fig. 4 reveal that the
both systems achieve higher sum rates when DL CSI is adopted at the users for
signal decoding over the statistical CSI case. In particular, the achieve rate
performance of the proposed underlay spectrum sharing in cell-free massive
MIMO is compared with that of the co-located counterpart in Fig. 3 and Fig. 4.
This comparison shows that the cell-free version outperforms the co-located
case in terms of the sum rate of the underlay spectrum sharing. For instance,
in Fig. 4, the cell-free based secondary system achieves the rate gains of
$75.6\%$, $119.6\%$, and $130.0\%$ for the uniform power allocation, max-min
power allocation, and DL CSI cases, respectively, compared to that of the co-
located counterpart at $P_{P}=5$ dBW.

Figure 4: The secondary sum rate versus the primary transmit power for
$K=L=10$, $M=N=32$, $I_{T}=0\,dB$, $\sigma_{P}^{2}=\sigma_{S}^{2}=1$,
$\nu=2.4$, and $d_{0}=1$ m. Moreover, $P_{S}=P_{P}/2$ and
$P_{S_{n}}=\min{\left(P_{S},{I_{T_{1}}}\big{/}{Z_{1}},\cdots,{I_{T_{k}}}\big{/}{Z_{k}},\cdots,{I_{T_{K}}}\big{/}{Z_{K}}\right)}$.
Figure 5: The impact of user-centric AP clustering. The achievable sum rate
versus the primary transmit power with/without AP clustering for $K=L=10$,
$M=N=32$, $I_{T}=0$ dB, $\sigma_{P}^{2}=\sigma_{S}^{2}=1$, $\nu=2.4$, and
$d_{0}=1$ m.

In Fig. 5, the effect of AP clustering is investigated. To this end the
achievable sum rates of primary and secondary systems are plotted against the
primary transmit power per P-AP with/without AP clustering by adopting max-min
power allocation. Fig. 5 reveals that AP clustering boosts the achievable
rates of both primary and secondary systems. For example, at $P_{P}$ of 10
dBW, the primary and secondary systems achieve sum rate gains of $39.5\%$ and
$48.7\%$, respectively, when the user-centric AP clustering is adopted over
the case of uniform AP deployment without a predefined AP clustering scheme.
The reason for this behavior is that, when a certain number of APs is
allocated for a particular user, the CCI from the remaining users of the own
system and the other system is reduced. This reduced level of CCI translates
into achievable rate gains.

Figure 6: The secondary user rate versus PIT ($I_{T}$) for $K=L=4$, $M=N=32$,
$P_{P}=0$ dBW, $\sigma_{P}^{2}=\sigma_{S}^{2}=1$, $\nu=2.4$, and $d_{0}=1$ m.

In Fig. 6, the effects of secondary transmit power constraints on the
achievable secondary user rates are explored. Two sets of rate curves are
plotted as a function of PIT ($I_{T}$) for $M=N=32$ and with/without DL CSI at
SUs by keeping the maximum allowable transmit power $P_{S}=P_{P}/2$ at each
S-AP. Fig. 6 clearly reveals that the secondary user rates grow exponentially
in the low $I_{T}$ regime in the both cases. In high $I_{T}$ regime, the SU
rates for the both CSI cases saturate to a maximum. The reason for this
behavior is that in low $I_{T}$ regime, a high amount of secondary CCI at PUs
is allowed, whereas in high $I_{T}$ regime, the SU rates saturate when the
secondary transmit power meets the transmit power constraints in (28). Fig. 6
also reveals that the SU rates can be boosted when the estimated DL CSI from
beamformed DL pilots is used over the statistical/long-term DL CSI counterpart
for signal decoding at SUs. For example, $U_{S}(3)$ and $U_{S}(4)$ achieve
rate gains of 20.9 % and 18.6 %, respectively, by using the DL CSI at
$I_{T}=10$ dB compared to those of without DL CSI.

Figure 7: The secondary sum rate versus PIT ($I_{T}$) for $K=L=4$, $P_{P}=0$
dBW, $\sigma_{P}^{2}=\sigma_{S}^{2}=1$, $\nu=2.4$, and $d_{0}=1$ m.

In Fig. 7, the implication of the number of P-APs/S-APs is investigated by
plotting the achievable secondary sum rate as a function of the PIT ($I_{T}$)
with statistical DL CSI at the SUs. By varying the number of P-APs and S-APs
as $M=N=32$, $M=N=72$, $M=N=128$, $M=N=200$, and $M=N=288$, five sets of rate
curves are plotted. Fig. 7 shows that the maximum saturation of the secondary
sum rate is heavily depend on $N/M$. The reason for this behavior can be
described as follows: In the low $I_{T}$ regime, the secondary only can
transmit smaller powers because PUs can withstand only to very small level of
secondary CCI. Thus, those low secondary transmit powers result in smaller
secondary achievable rates. Moreover, when $I_{T}$ grows large, it allows to
inflect a high amount of secondary CCI at the PUs, and thus, the S-APs can
transmit high power and it exponentially increases the secondary sum rate.
Once the secondary power constraints in (28) are met, the secondary sum rate
saturates to a maximum.

Figure 8: A comparison between the achievable sum rate of NOMA and OMA for
$M=N=32$, $P_{P}=0$ dBW, $\sigma_{P}^{2}=\sigma_{S}^{2}=1$, $\nu=2.4$,
$\vartheta_{P_{ai}}=\vartheta_{S_{bj}}=0.2$, and $d_{0}=1$ m.

In Fig. 8, a comparison of OMA versus NOMA in terms of the achievable sum
rates of the primary and secondary systems is presented. To this end, the sum
rates are plotted as a function of the number of users that can be served
simultaneously. The analytical curves are plotted by using (72) and (73) for
primary and secondary systems, respectively. Since the coherence interval
$\tau_{c}=196$, the maximum number of users that can be served by OMA is
limited to 196. On the other hand, NOMA can serve more number of users than
that of OMA because of user clustering. However, Fig. 8 reveals that OMA
outperforms NOMA in the regime of low number of users. The reason for this
behavior is that the intra-cluster pilot contamination due to the shared
pilots among NOMA clusters and the residual interference caused by error
propagation from imperfect SIC hinder the achievable rates of NOMA in low user
regime. However, NOMA-aided cell-free underlay spectrum sharing is beneficial
in boosting the number of simultaneous served SUs and the achievable sum rates
at high user rate regime.

Figure 9: The impact of DL pilots in NOMA. The achievable sum rate versus the
primary transmit power for $M=N=32$, $A=B=2$, $K=L=2$, $I_{T}=0\,dB$,
$\sigma_{P}^{2}=\sigma_{S}^{2}=1$, $\nu=2.4$, and $d_{0}=1$ m. Moreover,
$P_{S}=P_{P}/2$ and
$P_{S_{n}}=\min{\left(P_{S},{I_{T_{1}}}\big{/}{Z_{1}},\cdots,{I_{T_{k}}}\big{/}{Z_{k}},\cdots,{I_{T_{K}}}\big{/}{Z_{K}}\right)}.$

In Fig. 9, the impact of DL pilots on the achievable sum rate of NOMA-aided
underlay spectrum sharing for cell-free massive MIMO is studied. The
analytical rate curves for NOMA with estimated DL CSI are plotted via our
analysis in (86) and (88). Fig. 9 shows that the achievable sum rate of the
primary system with/without DL pilots grows monotonically with $P_{P}$.
Nevertheless, the sum rate of the secondary system with/without DL channel
training gradually grows until a maximum in the low power regime, and then, it
starts to decrease as $P_{P}$ continues to grow large. This is because the
secondary transmit power reaches its maximum limit and also due to higher
levels of primary CCI at SUs as per the description of Fig. 4. Moreover, Fig.
9 reveals that the achievable rates can be boosted when the users adopt the
estimated DL CSI via beamformed pilots for signal decoding at PUs/SUs. For
instance, at a primary transmit power of 0 dBW, the primary system with DL
pilots achieves a sum rate gain of 52.4% over a system that only relies on
statistical DL CSI (without DL pilots) for signal decoding at the PUs.

##  IX Conclusion

The practical feasibility of deploying underlay spectrum sharing in cell-free
massive MIMO has been investigated by adopting UL/DL pilot-based channel
estimations, max-min based MOOP, OMA/NOMA and the corresponding achievable
rates for both primary and secondary systems. User fairness has been
guaranteed by adopting an MOOP-based max-min fairness algorithm for
P-APs/S-APs, while satisfying the secondary transmit power constraints, which
are subjected to PIT. The achievable rates for both systems have been derived
for locally estimated UL CSI at P-APs/S-APs and DL CSI at PUs/SUs. The impact
of DL channel estimation at PUs/SUs via beamforming of DL pilots to boost the
achievable rates has been studied. We reveal that our proposed MOOP-based max-
min transmit power control algorithm can significantly boost the achievable
rate of the both systems over the uniform power allocation by using carefully
designed secondary transmit power constraints at each S-AP. The effect of
user-centric cell-free massive MIMO in which the P-APs/S-APs are clustered to
serve a particular PU/SU has been investigated. The trade-offs between OMA and
NOMA in terms of the number of concurrently served PUs/SUs and the achievable
rates have been explored, and thereby, it has been revealed that NOMA-aided
underlay spectrum sharing can be beneficial in satisfying future massive
access demands in a cell-free set-up. It has been shown that the achievable
rates can be boosted when the PUs/SUs estimate DL channels from the beamformed
DL pilots and adopt this estimated DL CSI to decode signals instead of solely
relying on statistical DL CSI. Thus, our performance analysis establishes that
a secondary system can be operated within the same primary/licensed spectrum
without hindering the primary system performance in a cell-free massive MIMO
set-up.

## Appendix A

###  A-A Derivation of the MMSE estimates in (6a) and (6b)

By substituting (5a) into (6a), $\hat{f}_{mk}$ can be derived as [35]

$\displaystyle\hat{f}_{mk}$ $\displaystyle=$
$\displaystyle\frac{\mathbb{E}\\!\left[{\left(\sqrt{\tau_{p}P}\left.\\!{f}\\!\right.\\!^{*}\\!_{mk}+\sqrt{\tau_{p}P}\left.\\!{v}\\!\right.\\!^{*}\\!_{mk}+{n}_{P_{m}}^{*}\right)f_{mk}}\right]}{\mathbb{E}\\!\left[{|\sqrt{\tau_{p}P}f_{mk}+\sqrt{\tau_{p}P}v_{mk}+{n}_{P_{m}}|^{2}}\right]}y_{P_{mk}}$
(89) $\displaystyle=$
$\displaystyle\frac{\sqrt{\tau_{p}P}\mathbb{E}\\!\left[{\left.\\!{f}\\!\right.\\!^{*}\\!_{mk}{f}_{mk}}\right]}{\tau_{p}P\mathbb{E}\\!\left[{|f_{mk}|^{2}}\right]+\tau_{p}P\mathbb{E}\\!\left[{|v_{mk}|^{2}}\right]+\mathbb{E}\\!\left[{|{n}_{P_{m}}|^{2}}\right]}y_{P_{mk}}\stackrel{{\scriptstyle(a)}}{{=}}c_{P_{mk}}y_{P_{mk}},$

where $c_{P_{mk}}$ is defined in (7) and
${n}_{P_{m}}=\mathbf{n}_{P_{m}}\boldsymbol{\phi}_{P_{k}}^{H}$. The step $(a)$
is written by using the fact that $f_{mk}$, $v_{mk}$, and $n_{P_{m}}$ are
Gaussian random variables with zero means and then, by evaluating the
expectations of their squared norms [43]. Similarly, the MMSE estimate of
$g_{nl}$ can be derived as given in (6b).

###  A-B Derivation of $Z_{k}$ in (27)

By substituting (6b) and (5b) into (26), $Z_{k}$ is computed as

$\displaystyle Z_{k}$ $\displaystyle=$
$\displaystyle\sum\nolimits_{n=1}^{N}\delta_{S_{nk}}\eta_{S_{nk}}c_{S_{nk}}^{2}\mathbb{E}\\!\left[{\left|{u_{nk}\left(\sqrt{P_{p}}g_{nk}+\sqrt{P_{p}}u_{nk}+n_{S_{n}}\right)^{*}}\right|^{2}}\right]$
(90) $\displaystyle+\sum\nolimits_{n=1}^{N}\sum\nolimits_{l\neq
k}^{L}\delta_{S_{nl}}\eta_{S_{nl}}c_{S_{nl}}^{2}\mathbb{E}\\!\left[{\left|{u_{nk}\left(\sqrt{P_{p}}g_{nl}+\sqrt{P_{p}}u_{nl}+n_{S_{n}}\right)^{*}}\right|^{2}}\right]$
$\displaystyle=$
$\displaystyle\sum\nolimits_{n=1}^{N}\\!\delta_{S_{nk}}\eta_{S_{nk}}c_{S_{nk}}^{2}\zeta_{{u}_{nk}}\\!\left({P_{p}}\\!\left(\zeta_{{g}_{nk}}\\!+\\!2\zeta_{{u}_{nk}}\right)\\!+\\!1\right)\\!+\\!\sum\nolimits_{n=1}^{N}\sum\nolimits_{l\neq
k}^{L}\\!\delta_{S_{nl}}\eta_{S_{nl}}c_{S_{nl}}^{2}\zeta_{{u}_{nk}}\\!\left({P_{p}}\left(\zeta_{{g}_{nl}}\\!+\\!\zeta_{{u}_{nl}}\right)\\!+\\!1\right)$
$\displaystyle=$
$\displaystyle\sum\nolimits_{n=1}^{N}\sum\nolimits_{l=1}^{L}\delta_{S_{nl}}\eta_{S_{nl}}\rho_{g_{nl}}\zeta_{{u}_{nk}}+\sum\nolimits_{n=1}^{N}\delta_{S_{nk}}\eta_{S_{nk}}\rho_{u_{nk}}^{2}.$

###  A-C Derivation of the MMSE estimate of $\mu_{P_{kk}}$ in (52)

By noting that $f_{mk}=\hat{f}_{mk}+\epsilon_{f_{mk}}$,
$\mathbb{E}\\!\left[{\mu_{P_{kk}}}\right]$ can be calculated as

$\displaystyle\\!\\!\\!\\!\mathbb{E}\\!\left[{\mu_{P_{kk}}}\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}\\!\left[{\sum\nolimits_{m=1}^{M}\delta_{P_{mk}}\eta_{P_{mk}}^{1/2}\hat{f}_{mk}\hat{f}_{mk}^{*}}\right]+\mathbb{E}\\!\left[{\sum\nolimits_{m=1}^{M}\delta_{P_{mk}}\eta_{P_{mk}}^{1/2}\epsilon_{f_{mk}}\hat{f}_{mk}^{*}}\right]=\sum\nolimits_{m=1}^{M}\delta_{P_{mk}}\eta_{P_{mk}}^{1/2}\rho_{f_{mk}},$
(91)

where $\rho_{f_{mk}}\triangleq\sqrt{\tau_{p}P}c_{P_{mk}}\zeta_{{f}_{mk}}$.
Then, $\mathbb{E}\\!\left[{\lambda_{P_{kk}}}\right]$ can be evaluated as

$\displaystyle\mathbb{E}\\!\left[{\lambda_{P_{kk}}}\right]$ $\displaystyle=$
$\displaystyle\mathbb{E}\\!\left[{\sum\nolimits_{n=1}^{N}\delta_{S_{nk}}\eta_{S_{nk}}^{1/2}u_{nk}\hat{g}_{nk}^{*}}\right]\stackrel{{\scriptstyle(b)}}{{=}}\mathbb{E}\\!\left[{\sum\nolimits_{n=1}^{N}\delta_{S_{nk}}\eta_{S_{nk}}^{1/2}c_{S_{nk}}u_{nk}\left(\sqrt{\tau_{p}P}\left(\hat{g}_{nk}^{*}+\hat{u}_{nk}^{*}\right)+n_{S_{n}}\right)}\right]$
(92) $\displaystyle=$
$\displaystyle\sum\nolimits_{n=1}^{N}\delta_{S_{nk}}\eta_{S_{nk}}^{1/2}\rho_{u_{nk}},$

where the step $(b)$ is written by using (6b) and (5b). The variance term of
$\mu_{P_{kk}}$ in (51) is given by

$\displaystyle\mathbb{V}\mathrm{ar}\\!\left[{\mu_{P_{kk}}}\right]=\mathbb{E}\\!\left[{|\mu_{P_{kk}}|^{2}}\right]-|\mathbb{E}\\!\left[{\mu_{P_{kk}}}\right]|^{2},$
(93)

where the first expectation term in (93) can be derived as

$\displaystyle\mathbb{E}\\!\left[{|\mu_{P_{kk}}|^{2}}\right]=\mathbb{E}\\!\left[{\left|\sum\nolimits_{m=1}^{M}\delta_{P_{mk}}\eta_{P_{mk}}^{1/2}\hat{f}_{mk}\hat{f}_{mk}^{*}\right|^{2}}\right]+\mathbb{E}\\!\left[{\left|\sum\nolimits_{m=1}^{M}\delta_{P_{mk}}\eta_{P_{mk}}^{1/2}\epsilon_{f_{mk}}\hat{f}_{mk}^{*}\right|^{2}}\right]$
$\displaystyle=\\!\\!\sum\nolimits_{m=1}^{M}\\!\\!\delta_{P_{mk}}\eta_{P_{mk}}\mathbb{E}\\!\left[{\\!|\hat{f}_{mk}|^{4}\\!}\right]\\!\\!+\\!\\!\sum\nolimits_{m=1}^{M}\\!\sum\nolimits_{m^{\prime}\neq
m}^{M}\\!\\!\delta_{P_{mk}}\eta_{P_{mk}}^{1/2}\delta_{P_{m^{\prime}k}}\eta_{P_{m^{\prime}k}}^{1/2}\mathbb{E}\\!\left[{\\!|\hat{f}_{mk}|^{2}|\hat{f}_{m^{\prime}k}|^{2}\\!}\right]$
$\displaystyle\quad+\sum\nolimits_{m=1}^{M}\\!\\!\delta_{P_{mk}}\eta_{P_{mk}}\\!\\!\left(\zeta_{{f}_{mk}}\\!-\\!\rho_{f_{mk}}\right)\rho_{f_{mk}}$
$\displaystyle=\\!2\\!\\!\sum\nolimits_{m=1}^{M}\\!\\!\\!\delta_{P_{mk}}\eta_{P_{mk}}\rho_{f_{mk}}^{2}\\!\\!+\\!\\!\sum\nolimits_{m=1}^{M}\\!\sum\nolimits_{m^{\prime}\neq
m}^{M}\\!\\!\\!\delta_{P_{mk}}\eta_{P_{mk}}^{1/2}\delta_{P_{m^{\prime}k}}\eta_{P_{m^{\prime}k}}^{1/2}\rho_{f_{mk}}\rho_{f_{m^{\prime}k}}\\!\\!+\\!\\!\sum\nolimits_{m=1}^{M}\\!\\!\\!\delta_{P_{mk}}\eta_{P_{mk}}\\!\left(\zeta_{{f}_{mk}}\\!\\!-\\!\rho_{f_{mk}}\right)\\!\rho_{f_{mk}}$
$\displaystyle=\sum\nolimits_{m=1}^{M}\sum\nolimits_{m^{\prime}=1}^{M}\delta_{P_{mk}}\eta_{P_{mk}}^{1/2}\delta_{P_{m^{\prime}k}}\eta_{P_{m^{\prime}k}}^{1/2}\rho_{f_{mk}}\rho_{f_{m^{\prime}k}}+\sum\nolimits_{m=1}^{M}\delta_{P_{mk}}\eta_{P_{mk}}\zeta_{{f}_{mk}}\rho_{f_{mk}}.$
(94)

By substituting (91) and (A-C) into (93), we have

$\displaystyle\mathbb{V}\mathrm{ar}\\!\left[{\mu_{P_{kk}}}\right]=\sum\nolimits_{m=1}^{M}\delta_{P_{mk}}\eta_{P_{mk}}\zeta_{{f}_{mk}}\rho_{f_{mk}}\triangleq
v_{P_{kk}}.$ (95)

The variance of $\lambda_{P_{kk}}$ can be given as

$\displaystyle\mathbb{V}\mathrm{ar}\\!\left[{\lambda_{P_{kk}}}\right]=\mathbb{E}\\!\left[{|\lambda_{P_{kk}}|^{2}}\right]-|\mathbb{E}\\!\left[{\lambda_{P_{kk}}}\right]|^{2}.$
(96)

Thus, the first expectation term in (96) can be derived as

$\displaystyle\mathbb{E}\\!\left[{|\lambda_{P_{kk}}|^{2}}\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}\\!\left[{\left|\sum\nolimits_{n=1}^{N}\delta_{S_{nk}}\eta_{S_{nk}}^{1/2}{u}_{nk}\hat{g}_{nk}^{*}\right|^{2}}\right]\stackrel{{\scriptstyle(c)}}{{=}}\sum\nolimits_{n=1}^{N}\delta_{S_{nk}}\eta_{S_{nk}}c_{S_{nk}}^{2}\zeta_{{u}_{nk}}\left(\tau_{p}P\left(\zeta_{{g}_{nk}}+2\zeta_{{u}_{nk}}\right)+1\right)$
(97) $\displaystyle=$
$\displaystyle\sum\nolimits_{n=1}^{N}\delta_{S_{nk}}\eta_{S_{nk}}\rho_{g_{nk}}\zeta_{{u}_{nk}}+\sum\nolimits_{n=1}^{N}\delta_{S_{nk}}\eta_{S_{nk}}\rho_{u_{nk}}^{2},$

where the step $(c)$ is written by following steps similar to those used in
the step $(b)$ in (92). Therefore, by using the results in (92) and (97), we
have

$\displaystyle\mathbb{V}\mathrm{ar}\\!\left[{\lambda_{P_{kk}}}\right]=\sum\nolimits_{n=1}^{N}\delta_{S_{nk}}\eta_{S_{nk}}\zeta_{{u}_{nk}}\rho_{g_{nk}}\triangleq
u_{P_{kk}}.$ (98)

Then, by using (91), (92) , (95), and (98), the effective channel estimate can
be derived as (52) [42].

##  Appendix B Derivation of SINR in (34)

The expectation in numerator of (34) can be derived as

$\displaystyle\mathbb{E}\\!\left[{\sum\nolimits_{m=1}^{M}\delta_{P_{mk}}{\eta_{P_{mk}}^{1/2}f_{mk}\hat{f}_{mk}^{*}}}\right]=\sum\nolimits_{m=1}^{M}\eta_{P_{mk}}^{1/2}\mathbb{E}\\!\left[{\left(\hat{f}_{mk}+\epsilon_{f_{mk}}\right)\left.\\!{f}\\!\right.\\!^{*}\\!_{mk}}\right]\stackrel{{\scriptstyle(d)}}{{=}}\sum\nolimits_{m=1}^{M}\delta_{P_{mk}}\eta_{P_{mk}}^{1/2}\rho_{f_{mk}},$
(99)

where $\epsilon_{f_{mk}}$ is an estimation error of $f_{mk}$ such that
$f_{mk}=\hat{f}_{mk}+\epsilon_{f_{mk}}$, satisfying
$\mathbb{E}\\!\left[{\epsilon_{f_{mk}}\left.\\!{f}\\!\right.\\!^{*}\\!_{mk}}\right]=0$.
In (99), the step $(d)$ is written by substituting (6a) and then evaluating
expectation term as

$\displaystyle\mathbb{E}\\!\left[{|\hat{f}_{mk}|^{2}}\right]=c_{P_{mk}}^{2}\mathbb{E}\\!\left[{|y_{P_{mk}}|^{2}}\right]=\sqrt{\tau_{p}P}\zeta_{f_{mk}}c_{P_{mk}}=\rho_{f_{mk}}.$
(100)

Then, the variance term in (34) can be derived as

$\displaystyle\mathbb{V}\mathrm{ar}\\!\left[{\sum\nolimits_{m=1}^{M}\delta_{P_{mk}}{\eta_{P_{mk}}^{1/2}f_{mk}\hat{f}_{mk}^{*}}}\right]$
$\displaystyle=\sum\nolimits_{m=1}^{M}\delta_{P_{mk}}\eta_{P_{mk}}\left(\mathbb{E}\\!\left[{\left|\left(\hat{f}_{mk}+\epsilon_{f_{mk}}\right)\left.\\!{\hat{f}}\\!\right.\\!^{*}\\!_{mk}\right|^{2}}\right]-\rho_{f_{mk}}^{2}\right)$
(101)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!=\sum\nolimits_{m=1}^{M}\delta_{P_{mk}}\eta_{P_{mk}}\left(2\rho_{f_{mk}}^{2}+\rho_{f_{mk}}\left(\zeta_{f_{mk}}-\rho_{f_{mk}}\right)-\rho_{f_{mk}}^{2}\right)=\sum\nolimits_{m=1}^{M}\delta_{P_{mk}}\eta_{P_{mk}}\rho_{f_{mk}}\zeta_{f_{mk}}.$

The expectation of the first term in (35) can be computed as

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\mathbb{E}\\!\left[{\left|\sum\nolimits_{i\neq
k}^{K}{\sum\nolimits_{m=1}^{M}\delta_{P_{mi}}{\eta_{P_{mi}}^{1/2}f_{mk}\hat{f}_{mi}^{*}}}\right|^{2}}\right]$
$\displaystyle\stackrel{{\scriptstyle(e)}}{{=}}\sum\nolimits_{i\neq
k}^{K}\sum\nolimits_{m=1}^{M}\mathbb{E}\\!\left[{\left|\delta_{P_{mi}}\eta_{P_{mi}}^{1/2}f_{mk}c_{P_{mi}}y_{P_{mi}}\right|^{2}}\right]$
(102)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\stackrel{{\scriptstyle(f)}}{{=}}\sum\nolimits_{i\neq
k}^{K}\sum\nolimits_{m=1}^{M}\delta_{P_{mi}}\eta_{P_{mi}}c_{P_{mi}}^{2}\zeta_{f_{mk}}\left(\tau_{p}P\left(\zeta_{f_{mi}}+\zeta_{v_{mi}}\right)+1\right)=\sum\nolimits_{i\neq
k}^{K}\sum\nolimits_{m=1}^{M}\delta_{P_{mi}}\eta_{P_{mi}}\rho_{f_{mi}}\zeta_{{f}_{mk}},$

where the steps $(e)$ and $(f)$ are written by using (6a) and (5a),
respectively. The expectation in the second term in (35) can be derived by
following steps similar to those used in (90) as follows:

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\mathbb{E}\\!\left[{\left|\sum\nolimits_{j=1}^{L}{\sum\nolimits_{n=1}^{N}\delta_{S_{nj}}{\eta_{S_{nj}}^{1/2}u_{nk}\hat{g}_{nj}^{*}}}\right|^{2}}\right]$
$\displaystyle=$
$\displaystyle\sum_{n=1}^{N}\delta_{S_{nk}}\eta_{S_{nk}}c_{S_{nk}}^{2}\zeta_{u_{nk}}\left(\tau_{p}P\left(\zeta_{g_{nk}}+2\zeta_{u_{nk}}\right)+1\right)$
(103) $\displaystyle+\sum\nolimits_{j\neq
k}^{L}\sum\nolimits_{n=1}^{N}\delta_{S_{nj}}\eta_{S_{nj}}c_{S_{nj}}^{2}\zeta_{u_{nk}}\left(\tau_{p}P\left(\zeta_{g_{nj}}+\zeta_{u_{nj}}\right)+1\right)$
$\displaystyle=$
$\displaystyle\sum\nolimits_{j=1}^{L}\sum\nolimits_{n=1}^{N}\delta_{S_{nj}}\eta_{S_{nj}}\rho_{g_{nj}}\zeta_{{u}_{nk}}+\sum\nolimits_{n=1}^{N}\delta_{S_{nk}}\eta_{S_{nk}}\rho_{u_{nk}}^{2}.$

By substituting (99), (101), (102), and (103) into (34), the desired SINR can
be derived as in (36).

##  Appendix C Derivation of SINR in (61)

The expectation of magnitude squared error $(\epsilon^{\mu}_{P_{kk}})$ in (60)
can be calculated as

$\displaystyle\mathbb{E}\\!\left[{|\epsilon^{\mu}_{P_{kk}}|^{2}}\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}\\!\left[{|\mu_{P_{kk}}-\hat{\mu}_{P_{kk}}|^{2}}\right]\stackrel{{\scriptstyle(g)}}{{=}}\frac{\left(1+P_{p,d}u_{P_{kk}}\right)^{2}\mathbb{V}\mathrm{ar}\\!\left[{\mu_{P_{kk}}}\right]+\left(P_{p,d}v_{P_{kk}}\right)^{2}\mathbb{V}\mathrm{ar}\\!\left[{\lambda_{P_{kk}}}\right]+P_{p,d}v_{P_{kk}}^{2}}{\left(P_{p,d}\left(v_{P_{kk}}+u_{P_{kk}}\right)+1\right)^{2}}$
(104) $\displaystyle=$
$\displaystyle\frac{\left(1+P_{p,d}u_{P_{kk}}\right)^{2}v_{P_{kk}}+\left(P_{p,d}v_{P_{kk}}\right)^{2}u_{P_{kk}}+P_{p,d}v_{P_{kk}}^{2}}{\left(P_{p,d}\left(v_{P_{kk}}+u_{P_{kk}}\right)+1\right)^{2}}\kappa^{\epsilon}_{P_{kk}},$

where $\kappa^{\epsilon}_{P_{kk}}$ is defined in the first term of (62) and
the step $(g)$ is written by substituting (52) and then, evaluating the
expectation. The expectation term of $\mu_{P_{ki}}$ in the denominator of (60)
can be derived as

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\mathbb{E}\\!\left[{|\mu_{P_{ki}}|^{2}}\right]\\!$
$\displaystyle=$
$\displaystyle\\!\mathbb{E}\\!\left[{\left|\sum\nolimits_{m=1}^{M}\\!\delta_{P_{mi}}\eta_{P_{mi}}^{1/2}\hat{f}_{mk}\hat{f}_{mi}^{*}\right|^{2}}\right]\\!\\!+\\!\mathbb{E}\\!\left[{\left|\sum\nolimits_{m=1}^{M}\\!\delta_{P_{mi}}\eta_{P_{mi}}^{1/2}\epsilon_{f_{mk}}\hat{f}_{mi}^{*}\right|^{2}}\right]$
(105) $\displaystyle=$
$\displaystyle\sum\nolimits_{m=1}^{M}\\!\delta_{P_{mi}}\eta_{P_{mi}}\zeta_{f_{mk}}\rho_{f_{mi}}\\!\triangleq\\!v_{P_{ki}}.$

Then, the expectation term in (60) with $\lambda_{P_{kj}}$ is evaluated
similar to the steps those used in (103) as

$\displaystyle\\!\\!\\!\\!\\!\\!\sum\nolimits_{j=1}^{L}\mathbb{E}\\!\left[{|\lambda_{P_{kj}}|^{2}}\right]\\!=\\!\sum\nolimits_{j=1}^{L}\sum\nolimits_{n=1}^{N}\delta_{S_{nj}}\eta_{S_{nj}}\rho_{g_{nj}}\zeta_{{u}_{nk}}\\!+\\!\sum\nolimits_{n=1}^{N}\delta_{S_{nk}}\eta_{S_{nk}}\rho_{u_{nk}}^{2}.$
(106)

The expectation term in the numerator of (60) is given by

$\displaystyle\mathbb{E}\\!\left[{|\hat{\mu}_{P_{kk}}|^{2}}\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}\\!\left[{|\mu_{P_{kk}}-\epsilon^{\mu}_{P_{kk}}|^{2}}\right]=\mathbb{E}\\!\left[{|\mu_{P_{kk}}|^{2}}\right]-\mathbb{E}\\!\left[{|\epsilon^{\mu}_{P_{kk}}|^{2}}\right].$
(107)

By substituting (A-C) and (104) into (107), we have

$\displaystyle\mathbb{E}\\!\left[{|\hat{\mu}_{P_{kk}}|^{2}}\right]=\sum\nolimits_{m=1}^{M}\sum\nolimits_{m^{\prime}=1}^{M}\delta_{P_{mk}}\eta_{P_{mk}}^{1/2}\delta_{P_{m^{\prime}k}}\eta_{P_{m^{\prime}k}}^{1/2}\rho_{f_{mk}}\rho_{f_{m^{\prime}k}}+v_{P_{kk}}-\kappa^{\epsilon}_{P_{kk}}.$
(108)

By substituting (104), (105), (106), and (108) into (60), the effective DL
SINR at $U_{P}(k)$ is evaluated as (61). By following step similar to
(104)-(108), the DL SINR at $U_{S}(l)$ can be computed as (65).

## References

  * [1] D. L. Galappaththige and G. Amarasuriya, “Cell-Free Massive MIMO with Underlay Spectrum-Sharing,” in _IEEE Int. Conf. on Commun. (ICC)_ , 2019, pp. 1–7.
  * [2] ——, “NOMA-Aided Cell-Free Massive MIMO with Underlay Spectrum-Sharing,” in _IEEE Int. Conf. on Commun. (ICC)_ , 2020, pp. 1–6.
  * [3] T. L. Marzetta, “Noncooperative Cellular Wireless with Unlimited Numbers of Base Station Antennas,” _IEEE Trans. Wireless Commun._ , vol. 9, no. 11, pp. 3590–3600, 2010.
  * [4] E. G. Larsson, O. Edfors, F. Tufvesson, and T. L. Marzetta, “Massive MIMO for Next Generation Wireless Systems,” _IEEE Commun. Mag._ , vol. 52, no. 2, pp. 186–195, 2014.
  * [5] “Sprint Unveils Six 5G-Ready Cities; Significant Milestone Toward Launching First 5G Mobile Network in the U.S.” https://newsroom.sprint.com/sprint-unveils-5G-ready-massive-MIMO-markets.htm, accessed: 2020-05-08.
  * [6] E. Björnson _et al._ , “Massive MIMO is a Reality? What is Next?: Five Promising Research Directions for Antenna Arrays,” _Digital Signal Process._ , vol. 94, pp. 3 – 20, 2019.
  * [7] H. Q. Ngo _et al._ , “Cell-Free Massive MIMO: Uniformly Great Service for Everyone,” in _IEEE 16th Int. Workshop on Signal Process. Adv. in Wireless Commun. (SPAWC)_ , June 2015, pp. 201–205.
  * [8] ——, “Cell-Free Massive MIMO versus Small Cells,” _IEEE Trans. Wireless Commun._ , vol. 16, no. 3, pp. 1834–1850, Mar. 2017.
  * [9] E. Nayebi _et al._ , “Precoding and Power Optimization in Cell-Free Massive MIMO Systems,” _IEEE Trans. Wireless Commun._ , vol. 16, no. 7, pp. 4445–4459, July 2017.
  * [10] T. C. Mai, H. Q. Ngo, M. Egan, and T. Q. Duong, “Pilot Power Control for Cell-Free Massive MIMO,” _IEEE Trans. Veh. Technol._ , vol. 67, no. 11, pp. 11 264–11 268, 2018.
  * [11] T. X. Doan, H. Q. Ngo, T. Q. Duong, and K. Tourki, “On the Performance of Multigroup Multicast Cell-Free Massive MIMO,” _IEEE Commun. Lett._ , vol. 21, no. 12, pp. 2642–2645, Dec 2017.
  * [12] Y. Li and G. A. Aruma Baduge, “NOMA-Aided Cell-Free Massive MIMO Systems,” _IEEE Wireless Commun. Lett._ , vol. 7, no. 6, pp. 950–953, 2018\.
  * [13] ——, “Underlay Spectrum-Sharing Massive MIMO NOMA,” _IEEE Commun. Lett._ , vol. 23, no. 1, pp. 116–119, Jan. 2019.
  * [14] J. Zhang _et al._ , “Prospective Multiple Antenna Technologies for Beyond 5G,” _IEEE J. Sel. Areas Commun._ , vol. 38, no. 8, pp. 1637–1660, 2020.
  * [15] K. Hosseini, W. Yu, and R. S. Adve, “Large-Scale MIMO versus Network MIMO for Multicell Interference Mitigation,” in _IEEE 15th Int. Workshop on Signal Process. Adv. in Wireless Commun. (SPAWC)_ , 2014, pp. 70–74.
  * [16] L. Dai, “A Comparative Study on Uplink Sum Capacity with Co-Located and Distributed Antennas,” _IEEE J. Sel. Areas Commun._ , vol. 29, no. 6, pp. 1200–1213, 2011.
  * [17] E. Koyuncu, “Performance Gains of Optimal Antenna Deployment in Massive MIMO Systems,” _IEEE Trans. Wireless Commun._ , vol. 17, no. 4, pp. 2633–2644, 2018.
  * [18] P. Liu, K. Luo, D. Chen, and T. Jiang, “Spectral Efficiency Analysis of Cell-Free Massive MIMO Systems with Zero-Forcing Detector,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 2, pp. 795–807, 2020.
  * [19] A. Goldsmith, S. A. Jafar, I. Maric, and S. Srinivasa, “Breaking Spectrum Gridlock with Cognitive Radios: An Information Theoretic Perspective,” _Proc. IEEE_ , vol. 97, no. 5, pp. 894–914, May 2009.
  * [20] H. Al-Hraishawi, G. A. Aruma Baduge, H. Q. Ngo, and E. G. Larsson, “Multi-Cell Massive MIMO Uplink with Underlay Spectrum Sharing,” _IEEE Trans. on Cogn. Commun. Netw._ , vol. 5, no. 1, pp. 119–137, 2019.
  * [21] L. Wang _et al._ , “Massive MIMO in Spectrum Sharing Networks: Achievable Rate and Power Efficiency,” _IEEE Syst. J._ , vol. 11, no. 1, pp. 20–31, March 2017.
  * [22] Y. Li, D. Kudathanthirige, and G. A. A. Baduge, “Massive MIMO Relay Networks with Underlay Spectrum Sharing,” _IEEE Trans. on Cogn. Commun. Netw._ , vol. 4, no. 4, pp. 677–691, Dec. 2018.
  * [23] S. Silva, M. Ardakani, and C. Tellambura, “Interference Suppression and Energy Efficiency Improvement with Massive MIMO and Relay Selection in Cognitive Two-Way Relay Networks,” _IEEE Trans. Green Commun. Netw._ , vol. 4, no. 2, pp. 326–339, 2020.
  * [24] D. Kudathanthirige and G. A. A. Baduge, “NOMA-Aided Multicell Downlink Massive MIMO,” _IEEE J. Sel. Areas Commun._ , vol. 13, no. 3, pp. 612–627, 2019.
  * [25] F. Rezaei, A. R. Heidarpour, C. Tellambura, and A. Tadaion, “Underlaid Spectrum Sharing for Cell-Free Massive MIMO-NOMA,” _IEEE Commun. Lett._ , vol. 24, no. 4, pp. 907–911, 2020.
  * [26] S. Kusaladharma and C. Tellambura, “Secondary User Interference Characterization for Spatially Random Underlay Networks with Massive MIMO and Power Control,” _IEEE Trans. Veh. Technol._ , vol. 66, no. 9, pp. 7897–7912, Sep. 2017.
  * [27] ——, “Massive MIMO based Underlay Networks with Power Control,” in _IEEE Int. Conf. on Commun. (ICC)_ , May 2016, pp. 1–6.
  * [28] S. Chaudhari and D. Cabric, “QoS Aware Power Allocation and User Selection in Massive MIMO Underlay Cognitive Radio Networks,” _IEEE Trans. on Cogn. Commun. and Netw._ , vol. 4, no. 2, pp. 220–231, 2018.
  * [29] M. Filippou, D. Gesbert, and H. Yin, “Decontaminating Pilots in Cognitive Massive MIMO Networks,” in _IEEE Int. Symp. Wireless Commun. Systems_ , Paris, France, Aug. 2012, pp. 816–820.
  * [30] H. Masoumi and M. J. Emadi, “Performance Analysis of Cell-Free Massive MIMO System with Limited Fronthaul Capacity and Hardware Impairments,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 2, pp. 1038–1053, 2020.
  * [31] M. Alonzo, S. Buzzi, A. Zappone, and C. D’Elia, “Energy-Efficient Power Control in Cell-Free and User-Centric Massive MIMO at Millimeter Wave,” _IEEE Trans. Green Commun. and Netw._ , vol. 3, no. 3, pp. 651–663, 2019.
  * [32] S. Chen _et al._ , “Structured Massive Access for Scalable Cell-Free Massive MIMO Systems,” _IEEE J. Sel. Areas Commun._ , pp. 1–1, 2020.
  * [33] T. L. Marzetta, E. G. Larsson, H. Yang, and H. Q. Ngo, _Fundamentals of Massive MIMO_. Cambridge University Press, Cambridge, UK, 2016.
  * [34] 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.
  * [35] S. M. Kay, _Fundamentals of Statistical Signal Processing: Estimation Theory_. Englewood Cliffs, NJ: Prentice Hall, 1993.
  * [36] E. Björnson and L. Sanguinetti, “Scalable Cell-Free Massive MIMO Systems,” _IEEE Trans. Commun._ , vol. 68, no. 7, pp. 4247–4261, 2020\.
  * [37] H. Q. Ngo and E. G. Larsson, “No Downlink Pilots are Needed in TDD Massive MIMO,” _IEEE Trans. Wireless Commun._ , vol. 16, no. 5, pp. 2921–2935, May 2017.
  * [38] P. Marbach, “Priority Service and Max-Min Fairness,” _IEEE/ACM Trans. Netw._ , no. 5, pp. 733–746, Oct 2003.
  * [39] B. Radunovic and J. Le Boudec, “A Unified Framework for Max-Min and Min-Max Fairness with Applications,” _IEEE/ACM Trans. Netw._ , vol. 15, no. 5, pp. 1073–1083, Oct 2007.
  * [40] L. Zheng, D. W. H. Cai, and C. W. Tan, “Max-Min Fairness Rate Control in Wireless Networks: Optimality and Algorithms by Perron-Frobenius Theory,” _IEEE Trans. Mobile Comput._ , vol. 17, no. 1, pp. 127–140, Jan 2018.
  * [41] E. Björnson, E. A. Jorswieck, M. Debbah, and B. Ottersten, “Multiobjective Signal Processing Optimization: The way to Balance Conflicting Metrics in 5G Systems,” _Signal Process. Mag._ , vol. 31, no. 6, pp. 14–23, Nov 2014.
  * [42] G. Interdonato, H. Q. Ngo, E. G. Larsson, and P. Frenger, “How Much Do Downlink Pilots Improve Cell-Free Massive MIMO?” in _IEEE Glob. Commun. Conf. (GLOBECOM)_ , Dec 2016, pp. 1–7.
  * [43] H. Cramer, _Random Variables and Probability Distributions_. Cambridge University Press, 1970.