Code-Domain NOMA in Massive MIMO: When Is It Needed?

1 Code-domain NOMA in Massive MIMO: When is it Needed? Mai T. P. Le, Luca Sanguinetti, Senior Member, IEEE, Emil Björnson, Senior arXiv:2003.01281v1 [cs.IT] 3 Mar 2020 Member, IEEE, Maria-Gabriella Di Benedetto, Fellow, IEEE Abstract In overloaded Massive MIMO (mMIMO) systems, wherein the number K of user equipments (UEs) exceeds the number of base station antennas M , it has recently been shown that non-orthogonal multiple access (NOMA) can increase the sum spectral efficiency. This paper aims at identifying cases where code-domain NOMA can improve the spectral efficiency of mMIMO in the classical regime in which K < M . Novel spectral efficiency expressions are provided for the uplink and downlink with arbitrary spreading signatures and spatial correlation matrices. Particular attention is devoted to the planar arrays that are currently being deployed in pre-5G and 5G networks and which are characterized by limited spatial resolution. Numerical results show that mMIMO with such planar arrays can benefit from NOMA in scenarios where the UEs are spatially close to each other. A UE grouping scheme is proposed for NOMA-aided mMIMO systems that is applicable to the spatial correlation matrices of the UEs that are currently active in each cell. Numerical results are used to investigate the performance of the algorithm under different operating conditions and types of spreading signatures (orthogonal, sparse and random sets). The analysis reveals that orthogonal signatures provide the highest average spectral efficiency. Index Terms Massive MIMO, uniform linear array, planar rectangular array, spatial correlation matrices, codedomain NOMA, spectral efficiency, channel estimation, arbitrary spreading signatures. Part of this paper was presented at IEEE PIMRC 2019 [1]. Mai T. P. Le is with The University of Danang – University of Science and Technology, Da Nang, Viet Nam, and with the University of Rome “La Sapienza”, 00184 Rome, Italy ([email protected]). L. Sanguinetti is with the University of Pisa, Dipartimento di Ingegneria dell’Informazione, 56122 Pisa, Italy ([email protected]). E. Björnson is with the Department of Electrical Engineering (ISY), Linköping University, 58183 Linköping, Sweden ([email protected]). M.-G. Di Benedetto is with the University of Rome “La Sapienza”, 00184 Rome, Italy ([email protected]). 2 I. I NTRODUCTION Massive MIMO (mMIMO) [2] and Non-Orthogonal Multiple Access (NOMA) [3], [4] are two physical layer technologies that have received large attention in recent years. While mMIMO has already made it into the 5G standard [5], the NOMA functionality remains to be standardized. Since mMIMO will likely be a mainstream feature in 5G networks, it is important to determine if and how NOMA can improve its performance. This is the main topic of this paper. A. Related Work and Motivation Conventional multiple access schemes assign orthogonal resources to each user equipment (UE). This provides restricted/dedicated resources per UE but eliminates inter-UE interference. It is well-known that this approach is inefficient if the interference can be controlled in some other domain [3], [6]–[8]; the power and code domains are typically used for interference suppression in NOMA, while the spatial domain is used for mMIMO. While prior investigations addressed only one of these three domains, some recent works consider systems that combine NOMA and mMIMO. The vast majority of the state-of-the-art contributions in this direction investigate the performance of power-domain NOMA when combined with mMIMO (see [9]–[14] and references therein). The gains are, however, generally limited since, to be efficient, power-domain NOMA requires UEs channels to be non-orthogonal, while a core feature of mMIMO is to make UE channels nearly orthogonal [9]. On the other hand, the combination of code-domain NOMA with mMIMO has received limited attention so far. The investigation in [15] addresses the pilot transmission phase and analyzes two pilot structures, namely, orthogonal and superimposed deterministic pilots. It was shown that the superimposed approach achieves better performance in a high mobility environment with a large number of UEs. The uplink (UL) spectral efficiency and bit error rate performance of mMIMO with a code-domain NOMA scheme, called interleaved division multiple-access, were evaluated in [16] with a low-complexity iterative data-aided channel estimation scheme and different suboptimal detection schemes, such as maximal ratio (MR) and zero-forcing (ZF) combining. In [17], the authors considered the UL of an overloaded setting without any channel state information (CSI). Low density spreading signatures were applied and a blind belief propagation detector was proposed. In [18], the mean squared error of code-domain NOMA was considered as the performance metric of an overloaded mMIMO system. 3 In contrast to prior work, this paper provides a theoretical analysis of the combination of code-domain NOMA and mMIMO in the classical mMIMO regime, that is for an underloaded system. B. Contributions The spectral efficiency (SE) of a classical mMIMO system grows without bound as M → ∞ when the spatial correlation properties of the interfering UEs’ channels are sufficiently different [19], [20]. Nevertheless, the SE that is achieved at any finite M can potentially be improved. In particular, there might be use cases where the UEs are located close to each other, such as in public hubs like stadiums, offices in high-rise buildings, train stations, and public outdoor events, wherein the UEs’ spatial channel correlation properties may be very similar and, thus, a very large number of antennas is needed to deliver acceptable performance when relying solely on the spatial processing provided by classical mMIMO. Orthogonal time-frequency scheduling algorithms that deal with this situation are described in [21], [22], but can these potentially be improved using NOMA? The main objective of this paper is to answer a simple question: What are (if any) the potential benefits of code-domain NOMA with mMIMO in those use cases? To provide some intuition about the role that NOMA can play, Section II first considers the UL of a case study setup with a single cell, K = 2 active UEs and perfectly known line-of-sight (LoS) propagation channels. The base station (BS) is equipped with M = 64 antennas deployed on a uniform linear array (ULA) with half-wavelength spacing. The analysis is carried out for maximum ratio (MR) and minimum mean square error (MMSE) combining schemes for UEs that are located spatially close to each other such that the array cannot resolve the UE angles. This is known as an unfavorable propagation scenario in the mMIMO literature [2], [22]. The analysis is then extended in Sections III and IV to both the UL and DL of a general multicell mMIMO system with NOMA using arbitrary spreading signatures for pilot and data transmission. Novel general SE expressions are provided with arbitrary spatial correlation matrices, that are used to design combining and precoding schemes, and to evaluate system performance for two configurations of antenna arrays and channel models; that is, the 2D one-ring channel model for a ULA and the 3D one-ring channel model for a planar array. In Section V, these SEs are used to confirm the preliminary analysis of Section II for the case study setup with M = 64 and K = 2. To fully take advantage of NOMA in a general setup with multiple UE, in Section VI we propose a per-cell UE grouping algorithm based on the k-means algorithm and using the 4 chordal distance between spatial correlation matrices as a similarity score metric. The proposed algorithm is applicable to both overloaded and underloaded scenarios. C. Outline and notation The paper is organized as follows. Section II provides some intuition on why code-domain NOMA can be useful with mMIMO: a case study setup with a single-cell network, two UEs and deterministic LoS channels. Section III introduces a general signal model for NOMA-aided mMIMO with multicell operation, arbitrary spreading signatures and spatial correlation matrices. The achievable SEs in the UL and DL are derived in Section IV, and used to select the optimal combining and precoding schemes. Numerical results are used to quantify the SEs in the case study setup and to validate the intuition provided in Section II. A UE grouping algorithm is developed in Section VI. The performance of NOMA-aided mMIMO is evaluated in Section VII under different operating conditions. Conclusions are drawn in Section VIII. Notation: Lower-case boldface and upper-case boldface letters are used to denote column vectors (e.g., x, y) and matrices (e.g. X, Y), respectively, while scalars are denoted by lower/uppercase italic letters (e.g. x, y, X, Y ). We denote [xi ] and [X]i,j the ith element of the vector x and (i, j)th element of the matrix X, respectively. kxk2 denotes the L2 -norm of vector x, i.e. kxk2 = qP pP 2 , whereas the Frobenius norm of matrix X is denoted by kX k = 2 |[x] | i F i i,j |[Xi,j ]| . XT , X∗ , XH , trX, E{X} are the transpose, the complex conjugate, the conjugate transpose, the trace and the expectation of the matrix X, respectively. The operator ⊗ stands for the Kronecker product. CM ×N denotes the set of complex-valued N × M matrices. The circularly symmetric complex Gaussian distribution with zero mean and correlation matrix R is denoted by NC (0, R). II. A G ENTLE S TART: S INGLE - CELL DEPLOYMENT WITH TWO UE S AND LOS CHANNELS To showcase what benefits code-domain NOMA can bring in a multi-antenna system, we consider the UL of a single-cell network where the BS is equipped with a uniform linear array of M antennas with half-wavelength spacing, and receives signals simultaneously from K = 2 single-antenna UEs. We denote by hk ∈ CM for k = 1, 2 the channel between UE k and the BS. We further assume free-space LoS channels, leading to the following deterministic channel T √  response [2, Sec. 1.3.2]: hk = βk 1, ejπ sin(φk ) , . . . , ejπ(M −1) sin(φk ) where βk is the large-scale fading attenuation and φk ∈ [0, 2π) is the angle-of-arrival (AoA) from UE k, measured from the 5 broadside of the BS array. We assume that UEs use N-length spreading signatures for UL data transmission, where N is a positive integer. We call uk ∈ CN the spreading signature randomly assigned to UE k and assume that kuk k2 = N. The N × 2 matrix U = [u1 , u2 ] ∈ CN ×2 denotes the signature matrix. The received signal Y ∈ CM ×N for the duration of the spreading signatures is Y = s1 h1 uT1 + s2 h2 uT2 + N, (1) where si ∼ NC (0, p) is the data signal from UE i and N ∈ CM ×N is thermal noise with i.i.d. 2 elements distributed as NC (0, σul ). We assume that β1 = β2 = β and define the average received 2 signal-to-noise ratio (SNR) as SNRul = βp/σul . Note that, in the absence of spreading signatures, (1) reduces to the classical mMIMO signal model for the UL. To detect s1 from Y in (1), the BS uses the combining vector v1 ∈ CM N , multiplied by the vectorized version of Y, to obtain v1H vec (Y) = s1 v1H g1 + s2 v1H g2 + v1H vec (N) , (2)
where gk = vec hk uTk = uk ⊗ hk ∈ CM N for k = 1, 2 is the effective channel vector. By treating the interference as noise, the achievable SE for UE 1 is SE1 = 1 EU {log2 (1 + γ1 )} , N (3) where γ1 is the signal-to-interference-and-noise ratio (SINR) γ1 = |v1H g1 |2 1 |v1H g2 |2 + SNR vH v1 ul 1 (4) and the expectation is taken with respect to the random assignment of signatures. The pre-log factor 1 N accounts for the fraction of samples used for transmitting the spreading signatures and it is smaller than 1 as it would be the case with classical mMIMO. However, if the signatures are properly associated with the UEs, the SE can be higher. To better understand this, we now design the combiner v1 in (2), which must be selected as a function of {g1 , g2 }, rather than {h1 , h2 } as would be the case in classical mMIMO. We begin by considering the popular MR combining with perfect channel knowledge, defined as v1 = g1 , leading to γ1MR = 1 | M1 hH1 h2 |2 | N1 uH1 u2 |2 + 1 M N SNRul , (5) 6 given that1 g1H g1 = MN and |g1H g2 |2 = |hH1 h2 |2 |uH1 u2 |2 . We note that [2, Sec. 1.3.2]    sin(M Ω12 ) if sin(φ1 ) 6= sin(φ2 ) 1 H M sin(Ω12 ) h h2 =  M 1 1 if sin(φ1 ) = sin(φ2 ) (6) with Ω12 = π(sin(φ1 ) − sin(φ2 ))/2. The term | M1 hH1 h2 |2 | N1 uH1 u2 |2 accounts for the interference generated by UE 2 and MNSNRul represents the received SNR in the absence of interference. From (6), it follows that the interference is stronger when the AoAs are similar to each other. However, if the UEs are associated to orthogonal codes/signatures (i.e., uH1 u2 = 0), the interference vanishes irrespective of the similarity of the AoAs, and the SE grows without limit as SNRul → ∞. On the contrary, it saturates to log2 (1 + 1/| M1 hH1 h2 |2 ) with mMIMO, due to the residual interference. Instead of using the suboptimal MR combining, we note that γ1 in (4) is a generalized Rayleigh quotient with respect to v1 and thus is maximized by the minimum mean square error (MMSE) combining vector [2, Sec. 1.3.3]: v1 = 2 X gi giH + i=1 leading to γ1MMSE = g1H 1 IM N g2 g2H + SNRul !−1 1 IM N SNRul (a) g1 = MNSNRul !−1 g1 , | M1 hH1 h2 |2 | N1 uH1 u2 |2 1− 1 1 + M N SNR ul (7) ! (8) where (a) follows from the matrix inversion lemma. The above SINR contains the same terms as (5), but has a different structure. In (5), | M1 hH1 h2 |2 | N1 uH1 u2 |2 must be interpreted as the perfomance loss due to the cancellation of the interference generated by UE 2. Similar to MR combining, this performance loss increases as the signals arrive from similar angles, but can be controlled (or even reduced to zero) by using spreading signatures. To quantitatively compare the different schemes, Fig. 1 shows the SE of UE 1 when M = 64 and SNR = 0 dB with MR (Fig. 1a) and MMSE combining (Fig. 1b) combining schemes. The nominal angle of UE 1 is fixed at φ1 = 30◦ while the angle of UE 2 varies from −60◦ to −60◦ . NOMA is employed with spreading signatures of length N = 2, which are either taken from an orthogonal set or randomly picked up from an assemble of ±1. Irrespective of the combining scheme and type of spreading signatures, mMIMO-NOMA outperforms mMIMO when the UEs 1 (A ⊗ B)H = AH ⊗ BH and (A ⊗ B)(C ⊗ D) = AC ⊗ BD 7 7 SE of UE 1 [bit/s/Hz] 6 5 4 3 2 1 0 -60 mMIMO mMIMO-NOMA w/ orth. codes mMIMO-NOMA w/ rand. codes -30 0 30 60 Angle of interfering UE [degree] (a) MR combining 7 SE of UE 1 [bit/s/Hz] 6 5 4 3 2 1 0 -60 mMIMO mMIMO-NOMA w/ orth. codes mMIMO-NOMA w/ rand. codes -30 0 30 60 Angle of interfering UE [degree] (b) MMSE combining Fig. 1: SE with two code-domain NOMA approaches and mMIMO for M = 64 under LoS propagation with φ1 = 30◦. MR (Fig.1a) vs. MMSE combining (Fig. 1b) with perfect CSI are considered. are closely located, meaning in this case |φ2 − φ1 | ≤ 5◦ . The reason is that mMIMO is unable to spatially separate the UEs in this case. However, mMIMO achieves higher SE with both combining schemes already for |φ2 − φ1 | ≥ 8◦ , which is a relatively small angular difference. The bottom line message of Fig. 1 is that there exist specific cases where NOMA can provide benefits if utilized with BSs equipped with many antennas M, even when M ≫ K. However, several strong assumptions were made in this example; that is, single-cell operation with only 2 8 UEs and LoS propagation with perfect CSI. Moreover, the 64 antennas were deployed on a large uniform linear array with half-wavelength spacing, which is unlikely to be the case in practice [20]. The question thus is: What happens in the UL and DL of practical mMIMO networks where these assumptions are not met? III. S YSTEM M ODEL We consider an mMIMO network composed of L cells. The BS in each cell is equipped with M antennas and simultaneously serves K single-antenna UEs. We assume that the BSs and UEs operate according to a TDD protocol with a data transmission phase and a pilot phase for channel estimation. We consider the standard block fading TDD protocol [2, Sec. 2.1] in which each coherence block consists of τc channel uses, whereof τp are used for UL pilots, τu for UL data, and τd for DL data, with τc = τp + τu + τd . We denote by hjlk ∈ CM the channel between UE k in cell l and BS j. In each coherence block, an independent correlated Rayleigh fading channel realization is drawn:
hjlk ∼ NC 0M , Rjlk , (9) j where Rjlk ∈ CM ×M is the spatial correlation matrix. The normalized trace βlk = tr(Rjlk )/M is the average channel gain from BS j to UE k in cell l. A. Channel Modeling The spatial correlation matrix Rjlk describes both the array geometry and the multipath propagation environment. Models for generation of Rjlk with arbitrary array geometries and environments can be found in [2, Sec. 7.3]. In this paper, we consider the following two physically motivated models: 1) 2D one-ring channel model: This model considers a ULA with half-wavelength spacing j and average path loss βlk [21], [2, Sec. 2.6]. The antennas and UEs are located in the same horizontal plane, thus the azimuth angle is sufficient to determine the directivity. It is assumed that the scatterers are uniformly distributed in the angular interval [ϕjlk −∆, ϕjlk +∆], where ϕjlk is the nominal geographical angle-of-arrival (AoA) and ∆ is the angular spread. This makes the (m1 , m2 )th element of Rjlk equal to j Z ∆  j j βlk ejπ(m1 −m2 ) sin(ϕlk +ϕ) dϕ. Rlk m1 ,m2 = 2∆ −∆ (10) 9 2) 3D one-ring channel model: This model considers a uniform planar array with the halfwavelength horizontal and vertical antenna spacing [2, Sec. 7.3]. We consider a quadratic √ √ array consisting of M horizontal rows with M antennas each, which restricts M to be the square of an integer. In this case, the (m1 , m2 )th element of Rjlk is given by ZZ  j j 2 ) sin(θ) jπ(m1 −m2 ) cos(θ) sin(ϕ) e|jπ(m1 −m Rlk m1 ,m2 = βlk {z } |e {z } f (ϕ, θ)dϕdθ, Vertical correlation (11) Horizontal correlation where f (ϕ, θ) is the joint probability density function of the azimuth ϕ and elevation θ angles. Following [2, Sec. 7.3.2], the 3D model is implemented by assuming that the BS height is 25 m, the UE height is 1.5 m, and a uniform angular distribution is used. We adopt a relative small value azimuth ϕ = 2◦ thorough the paper. The elevation θ of each UE is defined based on its distance to the BS of interest [2, Sec. 7.3.2]. With fixed ϕ = 2◦ , θ in this model ranges from about 3◦ to about 43◦ . Although the 2D model has been commonly used in the mMIMO literature (cf. [21], [23]), the 3D model definitely better reflects the typical pre-5G and 5G mMIMO array configurations [24]. While a 64-antenna ULA can have a high angular resolution in the azimuth domain and no resolution in the elevation domain, an 8 × 8 planar array has a mediocre resolution in both domains. This might have an important impact on the spatial multiplexing capabilities, depending on where the UEs are located. B. Channel Estimation The UL pilot signature of UE k in cell j is denoted by the vector φjk ∈ Cτp and satisfies √ kφjk k2 = τp . The elements of φjk are scaled by the square-root of the pilot power pjk and transmitted over τp channel uses, giving the received signal Yjp ∈ CM ×τp at BS j: Yjp = K X √ |i=1 pji hjji φTji + {z Desired pilots } L K X X √ pli hjli φTli + Npj , l=1,l6=j i=1 | {z Inter-cell pilots } (12) |{z} Noise 2 where Npj ∈ CM ×τp is noise with i.i.d. elements distributed as NC (0, σul ). Note that we are not assuming mutually orthogonal pilot signatures, but arbitrary spreading signatures. Hence, the MMSE estimator of hjjk takes a more complicated form than in prior works, as described in [2, Sec. 3.2]. 10 Lemma 1. The MMSE estimate of hjli is b j = √pli φH ⊗ Rj h li li li where Yjp is given in (12) and Qjli L X K X = l′ =1 i′ =1
Qjli
−1
vec Yjp , 2 IM τp . pl′ i′ (φl′ i′ φHl′ i′ ) ⊗ Rjl′ i′ + σul (13) (14) b j is independent of h b j and has correlation matrix Cj = The estimation error h̃jli = hjli − h li li li E{h̃jli (h̃jli )H } = Rjli − Φjli with

−1 φli ⊗ Rjli . Φjli = pli φHli ⊗ Rjli (Qjli ) (15) Proof: The proof follows from standard results and is provided in Appendix A. Note that the MMSE estimate in (13) holds for any choice of pilot signatures {φli }, that can be arbitrarily taken from orthogonal, non-orthogonal, random, or sparse sets. In classical mMIMO, orthogonal pilot signatures are usually employed, leading to the following simplified MMSE estimation expression [2, Sec. 3.2]: b j = √pli Rj Qj h li li li where Qjli = X
−1
Yjp φli , (16) 2 IM , pl′ i′ τp Rjl′ i′ + σul (17) (l′ ,i′ )∈Pli and Pli is the set collecting the indices of UEs that utilize the same pilot as UE i in cell l. C. UL and DL data transmissions While classical mMIMO only uses spreading signatures for UL pilot transmission, mMIMO with NOMA utilizes N-length spreading signatures also for UL data transmission, N being a positive integer. We denote by ujk ∈ CN the spreading signature assigned to UE k in cell j and assume that kujk k2 = N. As for pilot transmission, the spreading signatures {ujk } are also selected from an arbitrary set and different options will be compared below. The received signal Yj ∈ CM ×N at BS j for the duration of a spreading signature is given by Yj = K X i=1 | sji hjji uTji {z } Intra-cell signals + L K X X sli hjli uTli + Nj , l=1,l6=j i=1 | {z Inter-cell interference } |{z} Noise (18) 11 where sli ∼ NC (0, pli ) is the data signal from UE i in cell l with pli being the transmit power 2 ). and Nj ∈ CM ×N is thermal noise with i.i.d. elements distributed as NC (0, σul In the DL, the transmitted signal Xj ∈ CM ×N is given by Xj = K X ςji Wji , (19) i=1 where ςjk ∼ NC (0, ρjk ) is the data signal intended for UE k in cell j and Wji ∈ CM ×N is the corresponding precoding matrix that determines the spatial directivity of the signal. The received signal yjk ∈ CN ×1 at UE k in cell j, during the transmission of a spreading signature, is H yjk = K X ςji (hjjk )H Wji + L K X X ςli (hjjk )H Wli + nHjk , (20) l=1,l6=j i=1 i=1 2 where njk ∈ CN ×1 is thermal noise with i.i.d. elements distributed as NC (0, σdl ). No a priori assumption is made on the structure of the precoding matrices {Wji }. In Section IV-B, those will be designed on the basis of the channel estimates and of spreading signatures used at the UEs for detection. IV. S PECTRAL E FFICIENCY In this section, we will compute the SEs that are achieved in the UL and DL when arbitrary spreading signatures are used and we will design the combining/precoding vectors. A. UL Spectral Efficiency To detect the data signal sjk from Yj in (18), BS j selects the combining vector vjk ∈ CM N , which is multiplied with the vectorized version of Yj to obtain j H H vjk vec (Yj ) = sjk vjk gjk + K X j H + sji vjk gji i=1,i6=k | {z } Intra-cell interference
where glij = vec hjli uHli ∈ CM N or, equivalently, L K X X H H sli vjk glij + vjk vec (Nj ), l=1,l6=j i=1 | {z Inter-cell interference } | {z Noise (21) } glij = uli ⊗ hjli = (uli ⊗ IM ) hjli , (22) is the effective channel vector with correlation matrix E{glij (glij )H } = (uli ⊗ IM ) Rjli (uHli ⊗ IM ) = b j = (uli ⊗ IM ) h bj . bj = uli ⊗ h (uli uH ) ⊗ Rj . The MMSE estimate of gj is obtained as g li li li li li li 12 Note that (21) is mathematically equivalent to the signal model of a classical mMIMO system where the effective channel vectors are distributed as j glk ∼ NC 0M , (ulk uHlk ) ⊗ Rjli and the effective channel estimates are distributed as j blk g ∼ NC 0M , (ulk uHlk ) ⊗ Φjlk
(23)
(24) with Φjlk given by (15). The key difference is the presence of the spreading signatures (used for UL pilot and data transmissions) in the distributions. The ergodic capacity in UL can thus be evaluated by using the well-established lower bounds developed in the mMIMO literature [2]. Lemma 2. If the MMSE estimator is used, an UL SE of UE k in cell j is SEul jk =
1 τu  ul E log2 1 + γjk N τc [bit/s/Hz] , (25) ul where the effective instantaneous SINR γjk is given in ul γjk = H vjk with L K P P l=1,l6=j i=1 Zj = j 2 H bjk pjk |vjk g | H blij (b glij ) + pli g L X K X l=1 i=1 K P i=1,i6=k H j j bji pji g (b gji ) + Zj ! 2 pli (uli uHli ) ⊗ Cjli + σul IM N . vjk (26) The expectation is with respect to the realizations of the effective channels. Proof: The proof follows the same steps as that of [2, Th. 4.1] and is therefore omitted. Unlike the case study example of Section II where perfect CSI was assumed, the pre-log factor 1 τu N τc in (25) accounts for the fraction of samples used for transmitting pilot and data signatures. Whenever N > 1, it is still smaller than τu , τc which would be the case with classical mMIMO. The SE expression in (25) holds for any combining vector and choice of spreading signatures j bjk in the data transmission. MR combining with vjk = g is a possible choice. Similar to (4), the expression in (26) has also the form of a generalized Rayleigh quotient. Thus, the vector that maximizes the SINR can be obtained as stated by the following lemma. Lemma 3. The SINR in (26) is maximized by vjk = pjk L X K X l=1 i=1 H blij (b pli g glij ) + Zj !−1 j bjk g , (27) 13 leading to  X j H ul γjk = pjk (b gjk ) −1 H j blij (b bjk pli g glij ) + Zj  g . (l,i)6=(j,k) (28) Proof: This result follows from [2, Lemma B.10]. j The combining vector vjk in (27) is a function of the effective MMSE estimates {b gjk }, b j } as would be the case in classical mMIMO. We call it NOMA MMSE (Nrather than {h jk MMSE) combining since it also minimizes the mean-squared error (MSE) MSEk = E{|sjk − H vjk vec (Yj ) |2 {b glij }}, that represents the conditional MSE between the data signal sjk and the H received signal vjk vec (Yj ), after receive combining. So far, we have not taken into account the structure of the spreading signatures {ujk }, thus the SE expressions hold for any set of signatures. We will now consider the special case when the signatures are selected from a set of mutually orthogonal vectors. In this case, the estimate of sjk at BS j is obtained by first correlating Yj with the spreading signature ujk and then by multiplying the processed data signal2 Yj ujk ∈ CM by the combining vector v̄jk ∈ CM . We let Cjk denote the set of the indices of all UEs that utilize the same spreading signature as UE k in cell j. It can be easily shown that the SINR is maximized by v̄jk = pjk X (l,i)∈Cjk with Z̄jk = P (l,i)∈Cjk pli Cjli + 2 σul I . N M b j (h b j )H pli h li li + Z̄jk !−1 bj , h jk (29) This leads to the maximum SINR value  ul b j )H γjk = pjk (h jk X −1 b j (h b j )H + Z̄jk  h bj . pli h li li jk (l,i)∈Cjk (30) B. DL Spectral Efficiency We assume that, to detect the data signal ςji from yjk in (20), UE k in cell j correlates yjk with its associated spreading signature ujk to obtain H zjk = yjk ujk = (hjjk )H Wjk ujk ςjk + K X i=1,i6=k 2 (hjjk )H Wji ujk ςji + K L X X (hjjk )H Wli ujk ςli + nHjk ujk . l=1,l6=j i=1 (31) The processed signal Yj ujk is a sufficient statistic for estimating sjk when the signatures are selected from a set of mutually orthogonal vectors, since then there is no loss in useful information as compared to using the originally received signal Yj ; see e.g. [2, App. C.2.1]. 14 We denote the vectorized version of Wli as wli = vec(Wli ) ∈ CM N and observe that
H j H ) wli . (hjjk )H Wli ujk = ujk ⊗ hjjk vec(Wli ) = (gjk (32) Hence, zjk reduces to j H zjk = (gjk ) wjk ςjk + K X j H (gjk ) wji ςji + i=1,i6=k L K X X j H (gjk ) wli ςli + nHjk ujk . (33) l=1,l6=j i=1 As in the UL, (33) is mathematically equivalent to the signal model of classical mMIMO. Characterizing the capacity is harder in the DL than in the UL since it is unclear how the UE j H should best estimate the effective precoded channel (gjk ) wjk needed for decoding. However, an achievable SE can be computed using the so-called hardening capacity bound, which has received great attention in the mMIMO literature [2, Sec. 4.3] and will be adopted here as well. Lemma 4. The DL ergodic channel capacity of UE k in cell j in mMIMO-NOMA is lower bounded by SEdl jk =
1 τd dl log2 1 + γjk N τc [bit/s/Hz], (34) dl where the effective SINR γjk is given as dl γjk = L P K P l=1 i=1 j H ρjk |E{wjk gjk }|2 l 2 |} ρli E{|wli gjk H − . j }|2 ρjk |E{wjk gjk H + (35) 2 σdl The expectations are with respect to the realizations of the effective channels. Proof: The proof follows the same steps as that of [2, Th. 4.6] and is hence omitted. As in the UL, the DL SE in (34) holds for any choice of precoding vectors and spreading signatures. Moreover, the pre-log factor is reduced by a factor N compared to what it would be in classical mMIMO (i.e., τd /τc ). Unlike the UL, optimal precoding design is a challenge since (34) depends on the precoding vectors {wli } of all UEs. A common heuristic approach relies on the UL-DL duality [2, Th. 4.8], which motivates to select the precoding vectors as scaled versions of the combining vectors: wjk = p vjk , E{||vjk ||2 } (36) where the scaling factor is chosen to satisfy the precoding normalization constraint E{||wjk ||2 } = 1. By selecting vjk according to one of the UL combining schemes described earlier, the corresponding precoding scheme is obtained. 15 The expectations in (35) can be computed for any arbitrary precoding scheme by using Monte Carlo simulations. However, similar to [2, Cor. 4.5], we can obtain the closed-form expressions when using MR precoding, as described in the following corollary. Corollary 1. If MR precoding is used with wjk = √ bjk g E{||b gjk ||2 } , the expectations in (35) become

j H |E{wjk gjk }|2 = pjk tr ujk uHjk ⊗ Φjjk and l 2 E{|wli gjk |} H = tr (37)



ujk uHjk ⊗ Rljk (uli uHli ) ⊗ Φjli
. tr (uli uHli ) ⊗ Φlli (38) If the spreading signatures {ujk } are selected from a set of mutually orthogonal vectors, then we can choose Wjk = w̄jk uHjk where w̄jk ∈ CM is the precoding vector associated to UE k in cell j. Therefore, (31) reduces to H zjk = yjk ujk = Nςjk (hjjk )H w̄jk + X Nςli (hjjk )H w̄li + nHjk ujk , (39) (l,i)∈Cjk from which the effective SINR in (35) reads as dl γjk = P (l,i)∈Cjk H ρli E{|w¯li H hljk |2 } − ρjk |E{w̄jk hjjk }|2 + If MR precoding is used with w̄jk form: dl γjk H ρjk |E{wjk hjjk }|2 bjk / =h q 2 σdl N . (40) b jk ||2}, then (40) can be computed in closed E{||h   −1 ρjk pjk τp tr Rjjk (Qjjk ) Rjjk = . 2    l l l −1 l l −1 l l tr Rjk Rli (Qli ) Rli ρli pjk τp tr Rjk (Qli ) Rli X X 2 σdl  +    ρli + N −1 −1 tr Rlli (Qlli ) Rlli tr Rlli (Qlli ) Rlli (l,i)∈Cjk (l,i)∈{Pjk ∩Cjk \(j,k)} | {z } | {z } Non-coherent interference Coherent interference (41) Unlike with mMIMO (e.g., [2, Cor. 4.7]), the strength of coherent and non-coherent interference terms is determined by how similar the spatial correlation matrices Rlli with (l, i) ∈ Cjk and (l, i) ∈ {Pjk ∩ Cjk \ (j, k)} are to Rljk . By assigning orthogonal spreading signatures to the UEs with similar channel conditions, the SE can be higher than with mMIMO. 16 TABLE I: Network parameters Parameter Value Cell size 250 m × 250 m UL noise power σ 2 = −94 UL and DL transmit powers pjk = ρjk = 20 dBm Samples per coherence block τc = 200 Distance between UE i in cell l and BS j dlij Large-scale fading coefficient for the channel between UE i in cell l and BS j j βli Shadow fading between UE i in cell l and BS j V. N UMERICAL ANALYSIS FOR THE CASE STUDY: = −148.1 − 37.6 log10  j dli 1 km  + Flij dB Flij ∼ N (0, 10) S INGLE - CELL DEPLOYMENT WITH TWO UE S To quantify the potential benefits of code-domain NOMA in mMIMO, we begin by considering the simple case study of Section II with L = 1, M = 64, and K = 2, and numerically evaluate the SE for the practical setup described in Table I. For brevity, the analysis is carried out in the UL and MR and MMSE combining using MMSE channel estimation are considered. When NOMA is employed, we assume that orthogonal codes of length N = 2 are assigned to the two UEs. The two practical channel models described in Section III-A are used. A. Is NOMA needed? Similar to Fig. 1, we investigate the SE behavior with respect to the UEs’ locations. We fix the nominal azimuth angle of one UE at 30◦ while we let the nominal azimuth angle of the second one vary from −90◦ to 90◦ . Following the setup in Fig. 1, we impose that the average 1 1 channel gain per antenna stays the same, i.e., β11 = β12 . Figure 2 shows the UL SE of UE 1 with classical mMIMO and mMIMO-NOMA for the 2D and 3D models. With the NOMA scheme, N-MMSE and N-MR perform exactly the same since N = 2 and thus no interference is present—this is why only the N-MMSE curve is reported. Both channel models are considered with a relatively small ASD of ∆ = 2◦ . We observe that classical mMIMO gives higher SE than NOMA in both 2D and 3D models for most of the angles of the interfering UE. Different results are obtained for the case in which the two UEs have very similar angles. This is a challenging setup characterized by unfavorable propagation, wherein NOMA can bring some benefit. 17 9 SE of UE 1 [bit/s/Hz] 8 7 MMSE MR N-MMSE 6 5 4 3 2 1 -90 -60 -30 0 30 60 90 60 90 Angle of interfering UE [degree] (a) UL with the 2D channel model 9 SE of UE 1 [bit/s/Hz] 8 7 MMSE MR N-MMSE 6 5 4 3 2 1 -90 -60 -30 0 30 Angle of interfering UE [degree] (b) UL with the 3D channel model Fig. 2: SE of UE 1 in a single-cell two-user setup with ∆ = 2◦ and M = 64 with mMIMO and mMIMO-NOMA for N = 4, assuming the nominal azimuth angle of the desired UE is fixed at 30◦ , as a function of the azimuth angle of the interfering UE ranging from −90◦ to 90◦ . The 2D (Fig. 2a) and 3D (Fig. 2b) channel models described in Section III-A are considered. For the 2D model, Fig. 2a shows that MMSE largely outperforms NOMA even in this poor favorable propagation condition. This is because MMSE is a sufficiently powerful scheme to reject the interference even when the UEs are very close in space. However, we notice that this is achieved at the cost of a higher computational complexity than with MR [2] since the complexity scales as M 3 . Figure 2 also shows that NOMA can provide some gain compared to 18 1 Variance in (42) 0.8 2D, Gaussian = 2o 3D, Gaussian Uncorrelated = 2o 0.6 0.4 0.2 0 -90 -60 -30 0 30 60 90 Angle of interfering UE [degree] Fig. 3: Behaviour of the variance defined in (42) for the same setup of Fig. 2. Uncorrelated fading is also reported for comparison with 2D and 3D channel models, described in Section III-A. MR, without any increase in complexity. For the 3D model, Fig. 2b reveals that, when the UEs are close in space, NOMA provides the highest SE irrespective of the combining scheme used with mMIMO. This is because the planar array has a smaller spatial resolution, that reduces the spatial interference rejection capabilities of mMIMO and opens the door for complementing it with NOMA. B. A look at the favorable propagation conditions To better understand the above results, Fig. 3 shows the variance ( ) 1 H 1 tr (R111 R112 ) (h ) h 11 12 1 = δ1,12 =V p 1 1 M 2 β11 β12 E{kh111 k2 }E{kh112 k2 } (42) of the two UEs for 2D and 3D models in the same setup of Fig. 2. The variance is quantitatively measuring the level of favorable propagation [2, Eq. (2.19)]. It takes values in the interval 1 δ1,12 ∈ [0, 1], where smaller values represent a higher level of favorable propagation. Specifically, 1 δ1,12 = 1 if R111 and R112 are rank one and have the same dominant eigenvector. In contrast, 1 δ1,12 = 0 if the correlation matrices R111 and R112 are orthogonal, i.e., tr (R111 R112 ) = 0, which is a special case of linearly independent correlation matrices. Note that full orthogonality is unlikely to appear in practice [20]. The variance in (42) equals 1/M 2 for uncorrelated fading channels. However, Fig. 3 shows that the values of (42) changes with angles when considering the 2D and 3D channel models. It 19 achieves its maximum value at 30◦ for both models, which coincides with the angle giving the lowest SE values in Fig. 2. With the 2D model, the peak variance is relatively small (≈ 0.25), leading to comparatively good favorable propagation conditions. This justifies why classical mMIMO performs fairly well in the setup of Fig. 2. On the other hand, the variance is substantially larger (≈ 0.95) with the 3D model. This is because both horizontal and vertical spatial resolutions of the 8 × 8 array is only given by 8 antennas. Therefore, separation of the UEs in any of the two domains cannot be achieved. Hence, the two UEs cause much interference to each other, and thus the SE of mMIMO deteriorates, especially with MR. As shown in Fig. 2, this issue can be solved with NOMA by assigning orthogonal spreading signatures to the UEs with similar channel conditions. A natural question is thus how to group the UEs in a cell into groups that offers favorable propagation conditions. This problem is addressed in the next section. VI. UE G ROUPING A LGORITHM The concept of grouping UEs in mMIMO based on their spatial correlation matrices was introduced in [21], but for the purpose of orthogonal time-frequency scheduling when the UEs in each group have identical low-rank spatial correlation matrices. Channel measurements for mMIMO systems have recently shown that the spatial correlation matrices may have high rank, with a mix of several weak and a few strong eigendirections [25], [26], and vary even between closely spaced UEs. This implies that one cannot separate UEs into groups with orthogonal spatial correlation matrices to guarantee favorable propagation conditions, or expect UEs in the same group to have identical statistics. In other words, the grouping of UEs is highly non-trivial and will be addressed in this section. To this end, we first define the notion of dominant eigenspaces to capture the eigenspace that contains most of the energy of each correlation matrix. Definition 1 (p-Dominant eigenspace). Let A ∈ CM ×M be a Hermitian matrix with eigenvalue decomposition A = UDUH . The p-dominant eigenspace eigp (A) = [u1 . . . up ] is the (tall) unitary matrix composed of the p eigenvectors belonging to its p largest eigenvalues. The problem is how to group the UEs in a cell such that the p−dominating eigenspaces of the (possibly full-rank) correlation matrices of the UEs in each group are similar and different from the correlation matrices of other groups. A similarity score metric for measuring the difference between two eigenspaces is needed. A possible choice is given by the chordal distance. 20 Group 1 Group 2 Group 3 Group 4 Group 5 Group 6 Group 7 Group 8 y [m] 100 50 0 -100 -50 0 50 100 150 x [m] Fig. 4: Resulting association of UE positions to groups from the k-means algorithm with G = 8, p = 6 and a 3D one-ring channel model with a planar quadratic 8 × 8-antenna array. Definition 2 (Chordal distance). The chordal distance dC (A, B) between two matrices A and B is defined as dC (A, B) = kAAH − BBH k2F . (43) Several solutions exist in the literature to form groups on the basis of similarity scores. Among those, we adopt the k-means algorithm, which is widely used and operates as follows. For any cell j, k-means takes as inputs the set of intra-cell spatial correlation matrices {Rjjk ; k = 1, . . . , K}, the desired number of groups G, and the desired number of dominant eigenspace dimensions per group p. The output is a set of G tall unitary matrices {Ūg ∈ CN ×p : g = 1, . . . , G}, representing the center (or mean) of each group, and the sets {Cg : g = 1, . . . , G}, where Cg denotes the index set of UEs belonging to group g. The pseudo-code provided in Algorithm 1 summaries how the algorithm works. The k-means algorithm allows us to partition a cell into geographical regions, which are characterized by correlation matrices spanning almost orthogonal dominant eigenspaces. The algorithm can be applied ’offline’ to a very larger number of correlation matrices, which have been recorded over time to find static, but environment dependent, group spaces. Only the association of UEs to groups need to be computed at run-time. An example of offline grouping is provided in Fig. 4, which shows the resulting association of 1000 UE positions to G = 8 21 Algorithm 1: k-means algorithm for cell j. Input: {Rjjk ; k = 1, . . . , K}, G, and p Output: {Ū1 , . . . , Ūg }, {U1 , . . . , UK }, and {C1 , . . . , CG } /* Initialization */ t←0 for k = 1, . . . , K do Uk ← eigp (Rjjk ) end for g = 1, . . . , G do t Pick ig randomly out of the set {1, . . . , K} \ ∪g−1 g=1 Cg Cgt ← {ig } Ūg ← Uig end /* Iteratively group the p − dominating eigenspaces */ while {C1 , . . . , CG } 6= {C1t−1 , . . . , CGt−1 } do t←t+1 for g = 1 : . . . , G do C1i ← ∅ end for k = 1 : . . . , K do gk⋆ ← arg ming dC (Ūg , Uk ) Cgtk⋆ ← Cgtk⋆ ∪ k end for g = 1 : . . . , G do P  H U U Ūg ← eigp t k k k∈Cg Cg ← Cgt end end groups of p = 6 dimensions under the 3D one-ring model for a planar quadratic 8 × 8 antenna array with half-wavelength-spacing. The UEs are uniformly distributed over a 120◦ sector with 125 m radius. 22 If the number of active UEs is not very large, Algorithm 1 may provide some groups that are empty while others are overloaded (depending on UE positions). In this case, a further step in the k-means algorithm is needed, which assigns exactly N UEs to each group while minimizing the sum of the chordal distance pairs. This can be achieved by using the Hungarian method [27], which yields the modified algorithm described in Algorithm 2. Algorithm 2: UE grouping algorithm for cell j. Input: {Rjjk ; k = 1, . . . , K}, G, and p Output: {C1′ , . . . , CG′ } /* Step 1 : Run k-means algorithm */ {Ū1 , . . . , Ūg }; {U1 , . . . , UK }; /* Step 2 : Assignment problem */ for k = 1 : . . . , K do for g = 1 : . . . , G do dgk ← dC (Ūg , Uk ) end end /* Form Hungarian matrix by replicate each row of matrix dgk for N times */ DH [K, K] ← dgk [G × N, K] for k = 1 : . . . , K do Find the row r with least value dC for r = 1 : . . . , K do rk ← arg minr DH [r, k] ∪ (rk 6= rk′ ) Update the row index corresponding to group slot Cg′ ← Crk end end VII. P ERFORMANCE EVALUATION In this section, we compare the performance of mMIMO with and without NOMA and validate the benefits of the grouping algorithm. We consider a network with L = 4 cells and numerically 23 Average sum SE [bit/s/Hz/cell] 25 mMIMO mMIMO-NOMA w/o grouping mMIMO-NOMA w/ grouping 20 ☎ 15 10 MR ❅ ❘✄ ❅ ✂ 5 ✆ ■ ❅ ❅ MMSE 0 2 4 6 8 Number of groups Fig. 5: Average sum UL SE with mMIMO and mMIMO-NOMA when K = 8 UEs are uniformly distributed over a 30◦ sector. With mMIMO-NOMA, the UE groups are formed either in a random way (i.e., without grouping) or through Algorithm 2. Orthogonal spreading codes are used. evaluate the average sum SE in the UL and DL for the network setup defined in Table I. Each BS is located in the center of its cell, has M = 64 antennas and serves K UEs. The analysis is carried out with both MR and MMSE combining schemes, using MMSE channel estimation. Motivated by the results of Section V, only the 3D channel model with a 8 × 8 planar array and a relative small ∆ = 2◦ is considered. If not otherwise specified, we assume that τp = K orthogonal pilot sequences are used for channel estimation. A. Is UE grouping needed? We begin by assessing the benefits of properly grouping the UEs with mMIMO-NOMA. We consider the UL and assume that the number of UEs per cell is K = 8. From Fig. 2b, it follows that the SE is largely reduced when the UEs are located within a 30◦ sector. Therefore, we assume that the K = 8 UEs are uniformly and independently distributed over a 30◦ sector (oriented as in Fig. 4) at a distance of 100 m from the BS. Fig. 5 illustrates the average sum SE per cell with classical mMIMO and mMIMO-NOMA. The latter is operated without any grouping algorithm and with the grouping algorithm described in Algorithm 2, following Fig. 4 with p = 6. The number of formed groups is G = 2, 4 and 8, which implies that spreading signatures of length N = K/G = 4, 2 and 1 are assigned to the UEs in each group. The sequences are randomly 24 Average sum SE [bit/s/Hz/cell] 40 mMIMO mMIMO-NOMA w/o grouping mMIMO-NOMA w/ grouping 30 ☞ 20 ✌ ❅ ■ ❅ MMSE MR ❅ ❘✄ ❅ ✂ 10 0 8 16 24 32 Number of UEs (K) (a) UL transmission Average sum SE [bit/s/Hz/cell] 40 mMIMO mMIMO-NOMA w/o grouping mMIMO-NOMA w/ grouping 30 MMSE ☞ ✠ 20 ✌ MR ❅ ❘✄ ❅ ✂ 10 0 8 16 24 32 Number of UEs (K) (b) DL transmission Fig. 6: Average sum SE as a function of K with mMIMO and mMIMO-NOMA with random and grouping-based assignment. UL (Fig. 6a) and DL (Fig. 6b) transmission are considered. Orthogonal spreading codes are used. assigned to the active UEs when no grouping is used. Orthogonal spreading signatures are adopted. The results of Fig. 5 show that mMIMO-NOMA with Algorithm 2 achieves better performance of random grouping with both MR and MMSE combining. Compared to mMIMO, both approaches provide some gain, which maintains constant with MMSE when when passing from G = 2 to 4 while it reduces with MR. This is because MMSE combining has better 25 interference capabilities. With G = 8, we have that N = 1 and thus mMIMO-NOMA reduces to mMIMO. In summary, NOMA can bring some benefits compared to mMIMO even in the case that the spreading signatures are randomly assigned. Better performance can be achieved if assigned according to the spatial correlation matrices. The results are in agreement with those of the case study considered in Fig. 2, and confirm that there exist specific cases where NOMA can provide benefits even when M ≫ K. Similar results can be obtained for the DL. B. Varying number of UEs We now consider the case in which the number of active UEs, K, in each cell increases. Similarly to Fig. 5, we assume that the UEs are located close to each other. Unlike Fig. 5, however, we assume that they are equally distributed in four distinct circle clusters with radius r = 20 m, which have K/4 UEs each and are randomly deployed in each cell. This implies that the UEs are already grouped into G = 4 groups per cell. Spreading signatures of length N = K/4 are assigned to the K/4 UEs in each group. Orthogonal spreading signatures are adopted. Similar to Fig. 5, this might be a quite challenging setup for conventional mMIMO due to the insufficient spatial resolution of a planar BS array with M = 64. We compare classical mMIMO and mMIMO-NOMA with and without grouping-based signature assignment. The average sum SE in the UL and DL is shown in Fig. 6. With mMIMO-NOMA without grouping, the spreading sequences are randomly assigned to the UEs in the cell; this means that UEs in the same group can be assigned to the same spreading sequence. The result show that mMIMO-NOMA with proper assignment of sequences performs well in both UL and DL, particularly when using MMSE combining/precoding. mMIMO-NOMA achieves higher SE than classical mMIMO already with K = 8, and the gap slightly increases as K gets larger. With K = 32, the SE gain is 30% in the UL and 50% in the DL. The reason is that mMIMO-NOMA achieves a roughly constant sum SE as K increases, while the SE reduces with K for classical mMIMO due to the lack of favorable propagation conditions. The SE reduction is larger in the DL than in the UL, which might be due to the suboptimality of MMSE precoding and equal DL power allocation. C. Which spreading signatures are more favorable? We will now compare the achievable SE with spreading signatures taken from either orthogonal, random, or sparse sets. In the random case, the N-length signatures are picked up from 26 Average sum SE [bit/s/Hz/cell] 40 mMIMO mMIMO-NOMA w/ orth. codes mMIMO-NOMA w/ rand. codes mMIMO-NOMA w/ sparse codes 30 MMSE ☞ ✠ 20 ✌ MR ❅ ✄ ❘ ❅ 10 ✂ 0 8 16 24 32 Number of users, K Fig. 7: Average sum UL SE as a function of K with mMIMO-NOMA for different spreading signatures of length N = 4. an assemble of {±1}, whereas in the sparse case low-density signatures are used, which have only one non-zero value randomly distributed within the N-length signature. Fig. 7 shows the sum UL SE in the same setup of Fig. 6. We notice that orthogonal signatures give the highest performance with both MR and MMSE combining. While mMIMO-NOMA with orthogonal codes has superior performance for K ≥ 8, mMIMO-NOMA with random codes might provide some gains compared to mMIMO for K ≥ 32. This is because the probability that a given group of UEs is closely located in space increases as K becomes larger. Interestingly, mMIMO-NOMA with MR outperforms mMIMO only when orthogonal codes are used; this is because MR cannot deal with the extra interference coming from the non-orthogonality of random and sparse codes. Similar results are obtained for the DL, and thus are omitted due to space limitations. D. Impact of channel estimation quality The spatial interference rejection capabilities of mMIMO depends on the quality of channel estimates. So far, we have assumed that τp = K orthogonal pilot sequences are used for channel estimation. This is the common approach in mMIMO since it allows each BS to allocate orthogonal pilot sequences among its UEs, which are those originating the strongest interference. However, there might be use cases with stringent latency requirements in which only few samples τp can be dedicated to channel estimation. In these cases, τp will likely be smaller than K and thus UEs within the same cell can be assigned to the same pilot sequence. This gives rise to intra- 27 Average sum SE [bit/s/Hz/cell] 30 mMIMO mMIMO-NOMA 25 MMSE ☞ ✠ 20 ✌ 15 MR ❅ ✄ ❅ ❘ 10 5 ✂ 0 0 8 16 24 32 Number of pilots p Fig. 8: Average sum SE in UL with K = 32 UEs as a function of number pilot signatures, τp , with mMIMO and mMIMO-NOMA. Orthogonal spreading codes with N = 8 are adopted. cell pilot contamination, which inevitably deteriorates the SE of mMIMO. We now investigate if NOMA can bring some benefits in these cases. We consider the SE in the UL in the same setup of Fig. 6 where K = 32 UEs are equally distributed in four circle-areas of radius r = 20 m, which are randomly deployed in the cell area. Orthogonal spreading codes with length N = 8 are used for transmission and properly assigned to the different groups with mMIMO-NOMA. Fig. 8 shows that SE starts reducing when τp < 16 with both mMIMO and mMIMO-NOMA. However, the decrease in performance is slightly lower with mMIMO-NOMA because it does not rely only on the quality of channel estimates for dealing with interference. Particularly, a large gain is observed with NOMA when MMSE is used with only one channel use (i.e., τp = 1) for channel estimation. This is because MMSE is affected much from not having good channel estimates. VIII. C ONCLUSIONS This paper investigated the potential benefits of code-domain NOMA in mMIMO systems with limited spatial resolution. Novel general SE expressions for arbitrary spreading signatures and combining/precoding schemes were provided. We used these expressions to show, by means of simulations, that the SE can be improved by NOMA in cases when poor favorable propagation conditions are experienced by the UEs. This may happen when the UEs are located close to 28 each other and/or when planar arrays with insufficient resolution in the azimuth domain are considered. A per-cell grouping algorithm was developed and used to group the UEs with similar spatial correlation matrices. Numerical results showed that mMIMO NOMA may provide some gains if spreading sequences are assigned to the UEs within the same group. This is valid even if M ≫ K. The analysis was carried out with orthogonal, random, and sparse spreading signatures, revealing that orthogonal spreading sequences are the best choice. We also showed that some benefits can be achieved with NOMA when channel estimates of lower quality are available at the mMIMO BS. A PPENDIX A P ROOF OF L EMMA 1 The MMSE estimate of hjli is obtained as [28] o n o−1 n
j j p
H p
p
H b hli = E hli vec Yj E vec Yj vec Yj vec Yjp .
By using vec (ABC) = CT ⊗ A vec (B) we obtain n

H o √ √ E hjli vec Yjp = pli Rjli (φHli ⊗ IM ) = pli φHli ⊗ Rjli (44) (45) since the channels are independent. Similarly, one gets L X K n

H o X E vec Yjp vec Yjp = pl′ i′ (φl′ i′ ⊗ IM ) Rjl′ i′ (φHl′ i′ ⊗ IM ) + σ 2 IM τp l′ =1 i′ =1 = L X K X l′ =1 i′ =1 = L X K X l′ =1 i′ =1
pl′ i′ (φl′ i′ ⊗ IM ) φHl′ i′ ⊗ Rjl′ i′ + σ 2 IM τp pl′ i′ (φl′ i′ φHl′ i′ ) ⊗ Rjl′ i′ + σ 2 IM τp . (46) By substituting (45) and (46) into (44) yields (13). ACKNOWLEDGMENT The authors would like to acknowledge Jakob Hoydis for useful discussions in the development of Algorithm 1. 29 R EFERENCES [1] M. T. P. Le, L. Sanguinetti, E. Björnson, and M. D. Benedetto, “What is the Benefit of Code-domain NOMA in Massive MIMO?” in IEEE 30th Annual International Symposium on Personal, Indoor and Mobile Radio Communications (PIMRC), Sep. 2019, pp. 1–5. [2] E. Björnson, J. Hoydis, and L. Sanguinetti, “Massive MIMO networks: Spectral, energy, and hardware efficiency,” Foundations and Trends R in Signal Processing, vol. 11, no. 3-4, pp. 154–655, 2017. [Online]. Available: http://dx.doi.org/10.1561/2000000093 [3] L. Dai, B. Wang, Z. Ding, Z. Wang, S. Chen, and L. Hanzo, “A survey of non-orthogonal multiple access for 5G,” IEEE Communications Surveys Tutorials, vol. 20, no. 3, pp. 2294–2323, third quarter 2018. [4] M. T. P. Le, G. C. Ferrante, T. Q. S. Quek, and M. Di Benedetto, “Fundamental limits of low-density spreading NOMA with fading,” IEEE Transactions on Wireless Communications, vol. 17, no. 7, pp. 4648–4659, July 2018. [5] S. Parkvall, E. Dahlman, A. Furuskär, and M. Frenne, “NR: The new 5G radio access technology,” IEEE Communications Standards Magazine, vol. 1, no. 4, pp. 24–30, Dec 2017. [6] S. R. Islam, N. Avazov, O. A. Dobre, and K.-S. Kwak, “Power-domain non-orthogonal multiple access (NOMA) in 5G systems: Potentials and challenges,” IEEE Communications Surveys & Tutorials, vol. 19, no. 2, pp. 721–742, 2016. [7] M. T. P. Le, G. C. Ferrante, G. Caso, L. De Nardis, and M. Di Benedetto, “On information-theoretic limits of code-domain NOMA for 5G,” IET Communications, vol. 12, no. 15, pp. 1864–1871, 2018. [8] M. T. Le, G. Caso, L. De Nardis, A. Mohammadpour, G. Tucciarone, and M.-G. Di Benedetto, “Capacity bounds of LowDense NOMA over Rayleigh fading channels without CSI,” in IEEE 25th International Conference on Telecommunications (ICT), 2018, pp. 428–432. [9] K. Senel, H. V. Cheng, E. Björnson, and E. G. Larsson, “What role can NOMA play in massive MIMO?” IEEE Journal of Selected Topics in Signal Processing, vol. 13, no. 3, pp. 597–611, June 2019. [10] D. Kudathanthirige and G. A. A. Baduge, “NOMA-aided multicell downlink massive MIMO,” IEEE Journal of Selected Topics in Signal Processing, vol. 13, no. 3, pp. 612–627, June 2019. [11] D. Zhang, Z. Zhou, C. Xu, Y. Zhang, J. Rodriguez, and T. Sato, “Capacity analysis of NOMA with mmWave massive MIMO systems,” IEEE Journal on Selected Areas in Communications, vol. 35, no. 7, pp. 1606–1618, July 2017. [12] X. Chen, Z. Zhang, C. Zhong, R. Jia, and D. W. K. Ng, “Fully non-orthogonal communication for massive access,” IEEE Transactions on Communications, vol. 66, no. 4, pp. 1717–1731, 2017. [13] G. Yu, X. Chen, and D. W. K. Ng, “Low-cost design of massive access for cellular internet of things,” IEEE Transactions on Communications, vol. 67, no. 11, pp. 8008–8020, 2019. [14] A. S. de Sena, D. B. da Costa, Z. Ding, and P. H. J. Nardelli, “Massive MIMO–NOMA Networks With Multi-Polarized Antennas,” IEEE Transactions on Wireless Communications, vol. 18, no. 12, pp. 5630–5642, 2019. [15] J. Ma, C. Liang, C. Xu, and L. Ping, “On orthogonal and superimposed pilot schemes in massive MIMO NOMA systems,” IEEE Journal on Selected Areas in Communications, vol. 35, no. 12, pp. 2696–2707, Dec 2017. [16] C. Xu, Y. Hu, C. Liang, J. Ma, and L. Ping, “Massive MIMO, non-orthogonal multiple access and interleave division multiple access,” IEEE Access, vol. 5, pp. 14 728–14 748, 2017. [17] T. Wang, L. Shi, K. Cai, L. Tian, and S. Zhang, “Non-coherent NOMA with massive MIMO,” IEEE Wireless Communications Letters, 2019. [18] L. Liu, C. Yuen, Y. L. Guan, Y. Li, and C. Huang, “Gaussian message passing for overloaded massive MIMO-NOMA,” IEEE Transactions on Wireless Communications, vol. 18, no. 1, pp. 210–226, Jan 2019. 30 [19] E. Björnson, J. Hoydis, and L. Sanguinetti, “Massive MIMO has unlimited capacity,” IEEE Trans. Wireless Commun., vol. 17, no. 1, pp. 574–590, Jan. 2018. [20] L. Sanguinetti, E. Björnson, and J. Hoydis, “Towards massive MIMO 2.0: Understanding spatial correlation, interference suppression, and pilot contamination,” IEEE Trans. Commun., vol. 68, no. 1, pp. 232–257, 2020. [21] H. Huh, G. Caire, H. Papadopoulos, and S. Ramprashad, “Achieving “Massive MIMO” spectral efficiency with a not-solarge number of antennas,” IEEE Trans. Wireless Commun., vol. 11, no. 9, pp. 3226–3239, 2012. [22] T. L. Marzetta, E. G. Larsson, H. Yang, and H. Q. Ngo, Fundamentals of Massive MIMO. Cambridge University Press, 2016. [23] H. Yin, D. Gesbert, M. Filippou, and Y. Liu, “A coordinated approach to channel estimation in large-scale multiple-antenna systems,” IEEE J. Sel. Areas Commun., vol. 31, no. 2, pp. 264–273, Feb. 2013. [24] E. Björnson, L. Sanguinetti, H. Wymeersch, J. Hoydis, and T. L. Marzetta, “Massive MIMO is a reality—What is next? Five promising research directions for antenna arrays,” Digital Signal Processing, 2019. [25] X. Gao, O. Edfors, F. Rusek, and F. Tufvesson, “Massive MIMO performance evaluation based on measured propagation data,” IEEE Trans. Wireless Commun., vol. 14, no. 7, pp. 3899–3911, July 2015. [26] J. Flordelis, F. Rusek, F. Tufvesson, E. G. Larsson, and O. Edfors, “Massive MIMO performance—TDD versus FDD: What do measurements say?” IEEE Trans. Wireless Commun., vol. 17, no. 4, pp. 2247–2261, April 2018. [27] H. W. Kuhn, “The Hungarian method for the assignment problem,” Naval research logistics quarterly, vol. 2, no. 1-2, pp. 83–97, 1955. [28] S. M. Kay, Fundamentals of Statistical Signal Processing: Estimation Theory. Upper Saddle River, NJ, USA: Prentice-Hall, Inc., 1993.