{'fraction_non_alphanumeric': 0.03705332600329852, 'fraction_numerical': 0.010637713029136889, 'mean_word_length': 4.137108161536289, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 14, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Code-domain NOMA in Massive MIMO: When is it Needed?

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.


I. INTRODUCTION
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.
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 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 [x i ] and [X] i,j the ith element of the vector x and (i, j)th element of the matrix X, respectively. x 2 denotes the L 2 -norm of vector x, i.e. x 2 = i |[x] i | 2 , whereas the Frobenius norm of matrix X is denoted by X F = i,j |[X i,j ]| 2 . X T , X * , X H , 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. C M ×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 N C (0, R).

SINGLE-CELL DEPLOYMENT WITH TWO UES 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 h k ∈ C M 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 response [2, Sec. 1.3.2]: h k = √ β k 1, e jπ sin(φ k ) , . . . , e jπ(M −1) sin(φ k ) T where β k is the large-scale fading attenuation and φ k ∈ [0, 2π) is the angle-of-arrival (AoA) from UE k, measured from the 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 u k ∈ C N the spreading signature randomly assigned to UE k and assume that u k 2 = N. The N × 2 matrix U = [u 1 , u 2 ] ∈ C N ×2 denotes the signature matrix. The received signal Y ∈ C M ×N for the duration of the spreading signatures is where s i ∼ N C (0, p) is the data signal from UE i and N ∈ C M ×N is thermal noise with i.i.d.
elements distributed as N C (0, σ 2 ul ). We assume that β 1 = β 2 = β and define the average received signal-to-noise ratio (SNR) as SNR ul = βp/σ 2 ul . Note that, in the absence of spreading signatures, (1) reduces to the classical mMIMO signal model for the UL.
To detect s 1 from Y in (1), the BS uses the combining vector v 1 ∈ C M N , multiplied by the where g k = vec h k u T k = u k ⊗ h k ∈ C M N for k = 1, 2 is the effective channel vector. By treating the interference as noise, the achievable SE for UE 1 is where γ 1 is the signal-to-interference-and-noise ratio (SINR) 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 v 1 in (2), which must be selected as a function of {g 1 , g 2 }, rather than {h 1 , h 2 } as would be the case in classical mMIMO. We begin by considering the popular MR combining with perfect channel knowledge, defined as v 1 = g 1 , leading to with accounts for the interference generated by UE 2 and MNSNR ul 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., u H 1 u 2 = 0), the interference vanishes irrespective of the similarity of the AoAs, and the SE grows without limit as SNR ul → ∞. On the contrary, it saturates to log 2 (1 + 1/| 1 M h H 1 h 2 | 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 v 1 and thus is maximized by the minimum mean square error (MMSE) leading to where (a) follows from the matrix inversion lemma. The above SINR contains the same terms as (5), but has a different structure. In (5), | 1 M h H 1 h 2 | 2 | 1 N u H 1 u 2 | 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     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 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. SYSTEM MODEL
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 h j lk ∈ C M the channel between UE k in cell l and BS j. In each coherence block, an independent correlated Rayleigh fading channel realization is drawn: where R j lk ∈ C M ×M is the spatial correlation matrix. The normalized trace β j lk = tr(R j lk )/M is the average channel gain from BS j to UE k in cell l.

A. Channel Modeling
The spatial correlation matrix R j lk describes both the array geometry and the multipath propagation environment. Models for generation of R j lk 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 and average path loss β j 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 [ϕ j lk −∆, ϕ j lk +∆], where ϕ j lk is the nominal geographical angle-of-arrival (AoA) and ∆ is the angular spread. This makes the (m 1 , m 2 )th element of R j lk equal to 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 (m 1 , m 2 )th element of R j lk is given by Vertical correlation where f (ϕ, θ) is the joint probability density function of the azimuth ϕ and elevation θ angles. 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 φ jk 2 = τ p . The elements of φ jk are scaled by the square-root of the pilot power √ p jk and transmitted over τ p channel uses, giving the received signal Y p j ∈ C M ×τp at BS j: where N p j ∈ C M ×τp is noise with i.i.d. elements distributed as N C (0, σ 2 ul ). Note that we are not assuming mutually orthogonal pilot signatures, but arbitrary spreading signatures. Hence, the MMSE estimator of h j jk takes a more complicated form than in prior works, as described in [2,Sec. 3.2].
where Y p j is given in (12) and The estimation errorh j Proof: The proof follows from standard results and is provided in Appendix A.
Note that the MMSE estimate in (13) where and P li 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 u jk ∈ C N the spreading signature assigned to UE k in cell j and assume that u jk 2 = N. As for pilot transmission, the spreading signatures {u jk } are also selected from an arbitrary set and different options will be compared below. The received signal Y j ∈ C M ×N at BS j for the duration of a spreading signature is given by where s li ∼ N C (0, p li ) is the data signal from UE i in cell l with p li being the transmit power and N j ∈ C M ×N is thermal noise with i.i.d. elements distributed as N C (0, σ 2 ul ). In the DL, the transmitted signal X j ∈ C M ×N is given by where ς jk ∼ N C (0, ρ jk ) is the data signal intended for UE k in cell j and W ji ∈ C M ×N is the corresponding precoding matrix that determines the spatial directivity of the signal. The received signal y jk ∈ C N ×1 at UE k in cell j, during the transmission of a spreading signature, is where n jk ∈ C N ×1 is thermal noise with i.i.d. elements distributed as N C (0, σ 2 dl ). No a priori assumption is made on the structure of the precoding matrices {W ji }. 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. SPECTRAL EFFICIENCY
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 s jk from Y j in (18), BS j selects the combining vector v jk ∈ C M N , where g j li = vec h j li u H li ∈ C M N or, equivalently, is the effective channel vector with correlation matrix E{g j Note that (21) is mathematically equivalent to the signal model of a classical mMIMO system where the effective channel vectors are distributed as (23) and the effective channel estimates are distributed as with Φ j lk 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
where the effective instantaneous SINR γ ul jk is given in 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 N τu τ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 in the data transmission. MR combining with v jk = g j jk 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.
The combining vector v jk in (27) is a function of the effective MMSE estimates { g j jk }, rather than { h j jk } as would be the case in classical mMIMO. We call it NOMA MMSE (N-MMSE) combining since it also minimizes the mean-squared error (MSE) , that represents the conditional MSE between the data signal s jk and the received signal v H jk vec (Y j ), after receive combining. So far, we have not taken into account the structure of the spreading signatures {u jk }, 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 s jk at BS j is obtained by first correlating Y j with the spreading signature u jk and then by multiplying the processed data signal 2 Y j u jk ∈ C M by the combining vectorv jk ∈ C M . We let C jk 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 bȳ This leads to the maximum SINR value

B. DL Spectral Efficiency
We assume that, to detect the data signal ς ji from y jk in (20), UE k in cell j correlates y jk with its associated spreading signature u jk to obtain We denote the vectorized version of W li as w li = vec(W li ) ∈ C M N and observe that Hence, z jk reduces to 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 should best estimate the effective precoded channel (g j jk ) H w jk 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
where the effective SINR γ dl jk is given as The expectations are with respect to the realizations of the effective channels.
where the scaling factor is chosen to satisfy the precoding normalization constraint E{||w jk || 2 } = 1. By selecting v jk according to one of the UL combining schemes described earlier, the corresponding precoding scheme is obtained. and If the spreading signatures {u jk } are selected from a set of mutually orthogonal vectors, then we can choose W jk =w jk u H jk wherew jk ∈ C M is the precoding vector associated to UE k in cell j. Therefore, (31) reduces to from which the effective SINR in (35) reads as If MR precoding is used withw jk = h jk / E{|| h jk || 2 }, then (40) can be computed in closed form: 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 R l li with (l, i) ∈ C jk and (l, i) ∈ {P jk ∩ C jk \ (j, k)} are to R l jk . By assigning orthogonal spreading signatures to the UEs with similar channel conditions, the SE can be higher than with mMIMO.

UES
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 channel gain per antenna stays the same, i.e., β 1 11 = β 1 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. 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 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  1,12 = 0 if the correlation matrices R 1 11 and R 1 12 are orthogonal, i.e., tr (R 1 11 R 1 12 ) = 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 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 GROUPING ALGORITHM
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.
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 {R j jk ; 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 ∈ C N ×p : g = 1, . . . , G}, representing the center (or mean) of each group, and the sets {C g : g = 1, . . . , G}, where C g 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  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. for k = 1 : . . . , K do Find the row r with least value d C for r = 1 : . . . , K do r k ← arg min r D H [r, k] ∪ (r k = r ′ k ) Update the row index corresponding to group slot

VII. PERFORMANCE 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 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. and grouping-based assignment. UL (Fig. 6a) and DL (Fig. 6b)

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 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 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- 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. CONCLUSIONS
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 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.
APPENDIX A PROOF OF LEMMA 1 The MMSE estimate of h j li is obtained as [28] h By using vec (ABC) = C T ⊗ A vec (B) we obtain since the channels are independent. Similarly, one gets 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.