Exact BER Analysis of NOMA with Arbitrary Number of Users and Modulation Orders

—Non-orthogonal multiple access (NOMA) is a promising candidate for future mobile networks as it enables improved spectral-efﬁciency, massive connectivity and low latency. This paper derives exact and asymptotic bit error rate (BER) expressions under Rayleigh fading channels for NOMA systems with arbitrary number of users and arbitrary number of receiving antennas and modulation orders, including binary phase-shift keying and rectangular/square quadrature amplitude modulation. Furthermore, the power coefﬁcients’ bounds, which ensure users’ fairness, and solve the constellation ambiguity problem, are derived for N = 2 and 3 users cases with any modulation orders. In addition, this paper determines the optimal power assignment that minimizes the system’s average BER. These results provide valuable insight into the system’s BER performance and power assignment granularity. For instance, it is shown that the feasible power coefﬁcients range becomes signiﬁcantly small as the modulation order, or N , increases, where the BER performance degrades due to the increased inter-user interference. Hence, the derived expressions can be crucial for the system scheduler in allowing it to make accurate decisions of selecting appropriate N , modulation orders, and power coefﬁcients to satisfy the users’ requirements. The presented expressions are corroborated via Monte Carlo simulations.


I. INTRODUCTION
F UTURE WIRELESS networks are envisioned to provide ubiquitous and unlimited wireless coverage, which require integrating space, air, ground, and underwater networks into one large multidimensional network architecture [1].However, spectrum scarcity is one of the main challenges for realizing such ultra-wide wireless networks with massive connectivity.To this end, non-orthogonal multiple access (NOMA) has attracted tremendous attention as a promising candidate for future mobile networks because of its ability to provide high spectral efficiency, massive connectivity and low latency [2]- [7].Hence, much research was focused on the integration of NOMA in various applications, such as the Internet of Things (IoT), satellite communication, unmanned aerial vehicle (UAV) communications, and underwater communication [8]- [13].For example, Perez et al. [8] studied NOMA for IoT networks to provide reliable secure short packet communication for downlink and uplink.The work in [9], [10] investigated the application of NOMA in the forward link of multibeam satellite, whereas [13] studied its performance in underwater channels.Furthermore, a framework for UAVs serving ground users using NOMA is studied in [11], while the integration of NOMA with visible light communication (VLC) systems for indoor environments is discussed in [12].

A. Related Work
The widely considered power-domain NOMA (PD-NOMA), denoted as NOMA for short, is based on utilizing the powerdomain to multiplex different users' signals through superposition coding (SC), where distinct power coefficients are allocated to the users before combining their signals [14].The absence of orthogonality between users' signals introduces inter-user interference (IUI) which causes performance degradation to all users [15].Therefore, bit error rate (BER) and symbol error rate (SER) analysis of NOMA has received increased attention [13], [16]- [40].For example, Cejudo et al.
[18] attempted to approximate the BER using SER expressions for a two-user downlink NOMA, where each user may use a different modulation scheme.The considered modulations schemes are binary phase-shift keying (BPSK), quadrature phase-shift keying (QPSK), and quadrature amplitude modulation (QAM).Nonetheless, the SER analysis is limited to additive white Gaussian noise (AWGN) channels.On the other hand, the authors of [24] derived the exact closedform BER expressions for uplink two-user NOMA with QPSK modulation over AWGN channels.They assumed that this model is perfectly synchronized while considering imperfect successive interference cancellation (SIC) at the base station.In [28], exact BER expressions are derived for VLC-NOMA system with an arbitrary number of users employing on-off keying (OOK).
On the other hand, the BER for a single-input-single-output (SISO) Rayleigh fading wireless channel is considered in [29], where exact closed-form BER expressions are derived for the downlink while approximate expressions are derived for the uplink.However, these expressions are limited to a two-user NOMA considering QPSK for the near user and BPSK for the far user.The authors in [41] derived closedform expressions for the union bound on the BER of downlink NOMA with imperfect SIC over Nakagami-m fading channels.The tightness of the derived bounds varies based on various system parameters, and the gap between the bound and exact BER may exceed 3 dB.Furthermore, analytical expressions of the pairwise error probability (PEP) are given in [42] for an arbitrary number of users and modulation orders while considering imperfect SIC over Nakagami-m fading channels.
The presented analytical and simulation results show that the gap between the exact BER and PEP can be substantial, particularly for low and moderate signal-to-noise ratios (SNRs).Assaf et al. [31] derived exact BER closed-form expressions for downlink NOMA over SISO Nakagami-m fading channels for the two and three-user scenarios with QPSK.In [32], closed-form BER expressions are derived for a two-user downlink NOMA-VLC system while considering a limited set of modulation orders for phase-shift keying (PSK), pulse amplitude modulation (PAM), and QAM.In addition, Alqahtani and Alsusa [37] derived exact closed-form BER expressions for a two-user case employing BPSK in flat fading channels that are modelled by α-η-µ fading distribution to study the significance of different fading parameters on the BER performance.
Aldababsa et al. [36] presented closed-form BER expressions for an arbitrary number of users employing BPSK in Rayleigh flat fading, assuming perfect SIC.They also derived the range of proper power assignment for each user to ensure reliable BER performance.In [35], Assaf et al. extended the work in [31] for a two-user NOMA, where each user may use square QAM with arbitrary modulation.Additionally, proper power assignment was formulated to ensure fairness between the users and to avoid constellation points overlap.Besides the fact that the work is limited to the two-user scenario, modulation schemes such as BPSK and 8-QAM modulation orders are not considered.In [13] exact closed-form BER expressions are derived for VLC-NOMA system consisting of two users with OOK modulation in underwater environments.
Analytical SER expressions for NOMA are given in [18]- [23].The authors of [19] considered the two-user case in downlink NOMA using arbitrary QAM with imperfect SIC.In addition, the condition for proper power assignment is considered for the two users case.The approximated BER using SER is found to be inaccurate for high modulation orders, or at low SNR values [35].Moreover, the authors of [23] considered a threshold detector instead of SIC.It is found that the analytical performance of the proposed detector is very close to the SIC detector.Nonetheless, the SIC detector outperforms the threshold detector at low SNRs.A comprehensive survey of work that considers BER and SER of NOMA is given in Table I.In this table, SIC refers to imperfect SIC, while Approx.SIC represents perfect SIC.Additionally, DL, UL and SIMO stand for downlink, uplink and singleinput-multiple-output, respectively.

B. Motivation and Contribution
As can be noted from the surveyed literature, and to the best of the authors' knowledge, the existing work that considers BER analysis for downlink NOMA has one or more constraints in terms of the number of users, modulation order or accuracy.However, the availability of analytical BER analysis tools is indispensable for efficient system design and optimization.Hence, this paper considers the BER performance analysis of NOMA with an arbitrary number of users, where each user employs an arbitrary modulation order.The considered modulation schemes are BPSK and M -QAM with square and rectangular constellations.The obtained BER analysis is crucial for various applications such as adaptive modulation, resource allocation, user pairing, optimal power allocation and QoS requirements' satisfaction.The main contributions of this paper can be summarized as follows: • Derived closed-form BER expressions for downlink NOMA with arbitrary number of users, where each user may use BPSK, or M -QAM with square and rectangular constellations.The analysis is applicable to NOMA systems with receiver diversity as well.
• Derived the asymptotic BER to simplify the BER calculation at high SNRs.
• Evaluated the BER for different power assignments and provided insights into the error performance of large number of users and high modulation orders.• Derived closed-form expressions for the power coefficients' bounds (PCBs) to solve the constellation points ambiguity problem for N ∈ {2, 3}, where arbitrary modulation orders are also considered.• Evaluated the impact of changing the modulation order of certain user on the BER of other NOMA users, which is necessary for adaptive modulation and resource allocation operations.• Computed the optimal power assignments that minimize the system's average BER for N = 2 and 3 cases while considering the PCBs as linear and non-linear constraints.

C. Paper Organization
The rest of the paper is organized as follows.In Sec.II, the system and channel models are introduced.Then, with the aid of an example, the generalized BER expressions are derived in Sec.III for N NOMA users with arbitrary modulation orders while considering SISO and SIMO systems.Sec.IV demonstrates the analysis of the PCBs, while Sec.V presents analytical and Monte Carlo simulation results, as well as the optimal power assignments.Finally, Sec.VI concludes the paper with a summary of the main findings.

D. Notations
The notations used throughout the paper are as follows.Boldface uppercase and lowercase symbols, such as such as X and such as x, will denote matrices and row/column vectors, respectively.The transpose is denoted by (•) T , the Hermitian transpose is denoted by (•) H , and the denotes the Hadamard element-wise product.The real, complex, integer domains are denoted by R, C and Z, respectively.Moreover, B represents the set of binary numbers.Pr( The identity a × a matrix is denoted as I a , and the complex Gaussian random variable with a zero mean and σ 2 variance is denoted as CN (0, σ 2 ).

II. SYSTEM AND CHANNEL MODELS
In downlink NOMA, the base station multiplexes the information symbols of N users using the same radio resources by assigning each user a distinct power coefficient based on its channel conditions.Without loss of generality, we assume that N users are ordered in ascending order based on their average channel gain, i.e.
, where h n is the channel frequency response of the link between the base station and the nth user, i.e.U n , n ∈ {1, 2, . . ., N }.Therefore, the power assignment is performed such that a user with severe fading conditions is assigned higher power than a user with good channel conditions [31], [41].Consequently, the power coefficients α = [α 1 , α 2 , . . ., α N ] are assigned such that α 1 < α 2 < • • • < α N , where N n=1 α n = 1.Fig. 1 shows an illustrative diagram of the system model for a single cell with joint-multi-user maximum likelihood detector (JMLD) receivers.Therefore, the NOMA symbol is described by where x n is the information symbol of the nth user, which is drawn uniformly from a BPSK or M -QAM constellation χ n .The nth user modulation order is For QAM signals, the information symbols typically have   for each user is presented after scaling with its respective power coefficient α n .Fig. 2b shows the resultant constellation after superposition of U 1 and U 2 symbols.In Fig. 2c, the overall NOMA constellation is presented showing the real and imaginary components for each constellation point.The maximum imaginary component of the NOMA symbol is 2 + α2 2 as U 3 symbol does not have imaginary components.On the other hand, the maximum real component is The other real and imaginary components can be found by considering the different combinations of ν 1 , ν 2 , . . ., ν N .The bit-to-symbol mapping considered in this work follows the widely used model, where only the individual user bit mapping is based on Gray coding as shown in Fig. 2a.Therefore, the NOMA constellation will not be Gray coded.The constellation diagrams in Fig. 2 (3) For the example in Fig. 2c, It is worth noting that using nonlinear mapping for the NOMA constellation may provide some error rate performance improvement, however, the gained improvement is generally small and increases the receiver complexity [43], [44].The BER analysis for Gray coded NOMA generally follows the same approach used for linear mapping.At the receiver side, the received baseband signal in flat fading channels is written as where w n ∼ CN (0, σ 2 wn ) is the AWGN, and wn ∼ N (0, 0.5σ 2 wn ).In channels with small scale Rayleigh fading and large scale pathloss, the channel gain can be decomposed as n , Υ n is the distance between the base station and U n , and λ is the pathloss exponent.The coefficients 1 , 2 , . . ., N are mutually independent and identically distributed (i.i.d) random variables.
The most common schemes for multi-user detection of NOMA signals are the SIC and the JMLD detectors [34], [35].The main difference between SIC and JMLD is that the former attempts to cancel the interference of other users, while the latter detects the users' signals jointly without interference cancellation.Nevertheless, Assaf et al. [34] proved that the BER performance of SIC and JMLD is identical for the downlink NOMA under perfect knowledge of channel state information (CSI).Consequently, JMLD is considered in this work to enable a compact systematic analysis.Given that CSI is known perfectly at the receiver, the information symbols can be recovered using JMLD as follows, where { x 1 , x 2 , ..., x N } are the jointly detected N users' symbols, and x i represents the trail symbols for the ith user taken from the symbol alphabet χ i .

III. GENERALIZED BER ANALYSIS
Generally speaking, to evaluate the BER, the error events for all possible transmitted symbols should be considered.Nonetheless, when equally probable symbols are assumed, due to symbols' symmetry, the BER can be calculated by considering only the symbols in the first quadrant of the constellation diagram.

A. Decision Regions' Boundaries
To evaluate the BER, the decision regions' boundaries (DRBs) for each bit should be specified.This can be achieved by segmenting the NOMA constellation into q constellation diagrams each of which corresponds to a particular bit.For example, the binary representation of the top-left symbol in Fig. 2 is 01011 since q = 5.Consequently, the top-left Input: ḃk , p, P Output: Bk , ∀k bit in each of the 5 constellations will be 0, 1, 0, 1, 1, respectively.The remaining points in all bit constellations can be obtained by following the same approach.To generate the constellations more systematically, we define p ∈ C 1×2 q as a vector that contains all possible NOMA symbols.Then vector ḃ is generated by converting the symbols in p into binary, and ḃk ∈ B 1×2 q , k ∈ {1, 2, . . ., q}, is obtained by segmenting ḃ into q vectors ḃ1 , ḃ2 , . . ., ḃq , where ḃ1 contains MSB of all symbols.
The next step is to generate the scatter matrix P ∈ C u×v , which can be performed using Algorithm 1.The elements of matrix P are the NOMA symbols arranged exactly according to the constellation diagram.For the example in Fig. 3, P 1,1 = A 111 + A 110 .The values of u and v are given in Algorithm 1.Similarly, each kth bit vector, ḃk , will have a scatter matrix Bk ∈ B u×v .The algorithms to produce the scatter matrices P and Bk , ∀k are given in Algorithms 1 and 2, respectively.Algorithm 1 is based on finding the real and imaginary amplitude levels and creating the scatter matrix P by moving from left to right and top to bottom.Furthermore, Algorithm 2 is mainly based on the results from Algorithm 1, where symbol to binary mapping is done to find the kth bit for each constellation point in P. Note that the functions used in Algorithm 1 can be found in advanced mathematical software packages such as Matlab, where length(•) finds the length of a vector, unique(•) returns the unique elements of a vector, sort(•) orders the elements of a vector, and reshape(•) transforms the array size.
Once P is calculated, the real and imaginary primary DRBs, d ∈ R 1×u−1 and d ∈ R 1×v−1 , can be found by using a sliding window averaging filter with a window of size 2 whose output can be written as and where U ∈ {1, 2, . . ., u−1} and V ∈ {1, 2, . . ., v−1}.The kth bit DRBs can be computed using the scatter matrix Bk , where the DRBs appear if there is a bit flip.Therefore, exclusive-OR operator can be used over one row or one column of Bk to find the bit flip location.This can be expressed as follows, and where t k ∈ B 1×u−1 and t k ∈ B 1×v−1 .Using the constellation diagrams in Fig. 2, it can be noted that the kth bit flips in either t k or t k , but not in both at the same time.The total number of DRBs for the kth bit is given by where i .Therefore, the kth bit's DRBs, ḋk ∈ R 1×ϑ k , can be found by considering the indices of t k or t k where the entry is 1.These indices are stored in z k ∈ R 1×ϑ k , and hence ḋk elements can be found by Note that dk sorts the elements of ḋk in a descending order, which is required to find the bit error probability as well as the coefficients matrix which will be explained in the following subsection.The DRBs of all NOMA bits for m = [4, 4, 2] are shown in Fig. 4 using different color shading.

B. Euclidean Distance Computation
The BER expressions can be obtained by computing the Euclidean distance between the constellation points and the DRBs in dk ∀k.However, due to constellation diagram symmetry, it is sufficient to consider only the first quadrant.Thus, we are interested in p + which contains the first quadrant symbols of the scatter matrix P, where p + ∈ C 1×2 q−2 .Therefore, the displacement matrix for the kth bit, ∆ k ∈ R ϑ k ×2 q−2 , computes the displacement between the first quadrant constellation points and the DRBs which can be expressed by (12).Note that each column in ∆ k corresponds to a specific constellation point in p+ k , which is defined in (13), while each row corresponds to a specific DRB in dk .

C. Conditional BER Analysis
The conditional BER can be derived by considering all transmitted and received bit combinations and their relation to the DRBs.After exhaustive manipulations, the conditional BER per bit can be expressed as where and the squared Euclidean distance matrix E k can be calculated using ∆ k , where , . . ., N }, and σ 2 w = 0.5σ 2 w .Note that the user index is dropped for notational simplicity.The coefficients matrix is defined as . For ϑ k = 1, C k will reduce to a row vector of length 2 q−2 where c where g and ϕ i, g j : {≥ i, odd} or {< i, even} (−1) i , Otherwise .
(18) The expression in ( 14) considers all the constellation points in the first quadrant of the space diagram.Thus, the weighting factor of 1/2 q−2 is considered as these constellation points are equally probable.For the special case of identical BPSK modulation orders, this weighting factor becomes 1/2 N −1 .Furthermore, each column in Γ k corresponds to a specific constellation point where the sum over that column gives the probability of error for that constellation point.
To demonstrate ( 12)-(18), the example shown in Figs.2-4 is considered and C k , ∀k is found and shown in Fig. 5.For brevity, P 2 does not flip by moving vertically, i.e., it flips only by moving horizontally.The displacement between the NOMA word b E and d2 can be calculated as 2,1 . . ., ∆ The squared Euclidean distance for (19) can be calculated by squaring the vector elements, i.e. e (2) 2,1 , . . ., E . Thus, the first column of Γ k in (15) can be written as 2,1 , . . ., Γ .
The probability that b 2 is detected erroneously, given that the NOMA word b E is transmitted, can be calculated as By noting that wn ∼ N (0, σ 2 wn ), then it is straightforward to show that The coefficients c i,1 , ∀i can be calculated by noting that g (2) 1 = 3 for this case.Therefore, using (16) gives c T .Consequently, the conditional BER is calculated for b 2 given that the NOMA word b E is transmitted, and the same approach should be repeated for all constellation points in p+ k to compute the overall conditional BER per bit (14).Finally, the nth user conditional BER can be found by averaging its P (k) B |Γ k expressions ( 14).This can be written as

D. BER Analysis without Receiver Diversity
For i.i.d Rayleigh fading channels, the PDF f Γ is exponentially distributed [40], [45], hence, the nth user average BER can be computed by noting that, where . Substituting ( 14) into ( 23) and using the result of ( 24) yield It is worth noting that using the binomial series expansion, the asymptotic BER can be obtained by substituting [46, p. 185].
In the case of adaptive power assignment based on the channel statistics, the conditional BER in (23) can be used to derive the average BER for any channel model by using the corresponding PDF of Γ (k) i,j in (15).It is worth noting that the power adaptation process is performed after avenging over the PDF of Γ (k) i,j .The analysis for the Rayleigh channel model are given by (24).For the third type, the BER should be conditioned on both, the channel instantaneous fading coefficients and the used power coefficients.In such scenarios, the power coefficients depend on the fading coefficients.Hence, they can be written as α = F(h 1 , h 2 , . . ., h N ), where F(•) is a general function that depends on the adopted optimization criterion.Therefore, the BER can be evaluated by replacing α with F(h 1 , h 2 , . . ., h N ) in (23), and then averaging over the PDF of Γ (k) i,j .However, deriving the function F(•) and evaluating the average BER for this power assignment strategy is highly challenging because the relation between the power and channel coefficients is highly nonlinear [20], [21].In such scenarios, the desired power coefficients can be obtained using a particular numerical search method using (23).

F. BER Analysis with Receiver Diversity
In this work, we consider that the nth user receiver is equipped with L n receiving antennas, and the channels between the base station and all receiving antennas are i.i.d.Therefore, the received signal for the nth user is where {y n , h n , n , w n } ∈ C Ln×1 , and the entries of these vectors are defined similar to those of (4).The detector in this case can be realized as a maximal ratio combiner (MRC) followed by the JMLD.Therefore, (27) However, the BER for the SIMO case is similar to the SISO one except that the matrix Γ k will be replaced by k is equal to Γ k per receiving antenna.Therefore, the PDF of Ω f where Based on ( 28) and [33, Eq. ( 23)], the BER per bit can be evaluated as The nth user average BER can be computed by averaging its P (k) B expressions similar to (23).Moreover, a tight asymptotic BER per bit can be calculated as follows [45, pp. 326-327],

IV. POWER COEFFICIENTS' BOUNDS ANALYSIS
To enable a reliable detection of NOMA symbols using SIC, the power coefficient for each user should be selected such that the constellation for each user does not overlap with the other users' constellations [16], [19], [35], [36], [47].The PCB can be generally derived by noting that the nearest constellation point to the origin (NCO) in the first quadrant of the x-y plane should not cross the x or y axes.However, due to axes symmetry, it is sufficient to consider the real part of the NCO.Therefore, for N = 2, we have Consequently, the PCB can be written as α1 α2 < κ1 κ2 and the maximum possible power coefficient for U 1 can be obtained The PCB for N = 3 can be evaluated recursively by considering first the N = 2 case where U 1 and U 2 will be combined.Then the constellation that resulted from combining U 1 and U 2 constellations is considered as one constellation that will be combined with that of U 3 .Therefore, the first constraint that should be satisfied is identical to N = 2 case which is The second step is to ensure that the NCO in the combined constellation does not cross the y-axis to the negative side.Therefore, the constraint can be written as Solving the inequality in (32) for either α 1 or α 2 results in two solutions, ε (α n ) and (α n ) for n ∈ {1, 2}, but one of them is not applicable.The desired solution of ( 32) with respect to α 2 will be denoted as ε(α 1 ).Therefore, the second constraint becomes α 2 < ε(α 1 ).As an example, consider the case of m = [4,4,2] in Fig. 2, where the obtained solution is  32) which makes it easier to infer the PCBs.The intersection between (31) and ε (α 1 ) reflects the maximum possible power coefficient for U 1 , i.e., α 1,max .Therefore, the pair (α 1 , α 2 ) that satisfies the PCBs must be inside the region having the bounds of α 2 > α 1 and α Generally, the condition in (31) could intersect with the desired solution in (32) more than once.However, the desired α 1,max is the intersection that gives the minimum value of  For the arbitrary modulation orders cases, the PCBs are given by (31) and Table II summarizes the PCBs for the identical modulation schemes.It can be seen that the power coefficients space becomes smaller as the modulation order increases.Also, it is worth mentioning that the space of the power coefficients for N = 3 becomes narrower compared to N = 2 because of the extra PCB introduced.Following the same approach, the PCBs for other values of N can be derived utilizing the conditions in [47].
These closed-form PCBs expressions can be used as linear and non-linear constraints while solving minimization or maximization optimization problems.For example, by noting that the average BER of an N users NOMA system is given by P (N ) B,Avg.= 1 N N n=1 P Bn .Then the optimal power assignment that minimizes P (N =2) B,Avg. is formulated as subject to, where (35b) satisfies the PCB condition and (35c) is considered to ensure that the transmitted power is normalized to unity.The objective function in (35a) is non-linear, and hence, it is difficult to find a closed-form analytical solution for this problem.Therefore, the problem can be solved using interior-point optimization (IPO), which provides near-optimal solutions [31].Similarly, the optimization problem can be extended to N = 3 case.As such, the optimization problem is formulated as subject to, where (36b) is linear, and it is similar to N = 2 case.However, (36c) is an additional inequality constraint that is introduced for N = 3, where its upper bound is non-linear, whereas the lower bound is linear.In addition, (36d) is to ensure normalized transmission power.

V. NUMERICAL RESULTS AND DISCUSSIONS
This section presents the exact and asymptotic BER results of a downlink NOMA system using various number of users and modulation orders.The BER and asymptotic BER are computed analytically using the expressions in Sec.III, while the BER is validated by Monte Carlo simulations.In addition, the optimal power assignments that minimize the system's average BER for N = 2 and 3 cases are computed, where the PCBs, derived in Sec.IV, are used as constraints to solve the non-linear optimization problem using the IPO algorithm [31].It is worth noting that the PCB constraints increase from one constraint for N = 2, to three constraints for N = 3.Moreover, the small scale fading is considered to be flat and it follows the Rayleigh distribution with σ 2 n = 1.The large scale fading is considered as fixed pathloss with an exponent of λ = 2.7, where the users are at a normalized distance of Υ n = 10 3 5λ (n−1) from the base station.The AWGN variance σ wn is assumed to be common for all users, which corresponds to the transmit SNR 1/2σ 2 w [29], [31], [35], [40], [41].The base station and all users are assumed to be equipped with a single antenna unless stated otherwise.The power assignment for N = 2 is performed such that ( 31) is satisfied for all the considered modulation orders, where the worst-case scenario is when m = [8,64], and thus, α = [1 × 10 −2 , 0.99].For N = 3, the power assignment should satisfy (31) and (32) simultaneously for all the considered modulation orders, thus, α = [1 × 10 −4 , 1 × 10 −2 , 0.9899].The legends in all figures, where 2 ≤ N ≤ 4 represent the modulation order vectors m = [M 1 , M 2 , . . ., M n ].Furthermore, for N > 3 , the power coefficients are selected using linear search such that the constellation points do not overlap and the power coefficients order is maintained.Fig. 7 shows the analytical and simulation BER results for N = 2, where both users adopt identical modulation orders.The figure shows that the analytical and simulation results match very well for all the considered scenarios.Moreover, the asymptotic BER can be considered as an accurate approximation at high SNR values.As expected, increasing the modulation order degrades the BER performance, and the degradation generally follows the case of orthogonal multiple access.For example, QPSK modulation for both users requires about 3 dB additional power as compared to BPSK to achieve a BER of 10 −2 .Moreover, when comparing BPSK and 64-QAM, the latter requires about 14.4 and 17.5 dB for U 1 and U 2 , respectively, to achieve BERs of 10 −2 .
Fig. 8 shows the BER for U 1 and U 2 , where the modulation orders are not necessarily identical.It can be noted that at high SNR values, the BER for each user does not generally depend on the modulation order of the other user, except for m = [8,64].It is worth noting that the selected power coefficients for this case are close to the PCB (31), while for the rest of modulation orders combinations the selected power coefficients are relatively close to the center of the PCB.Therefore, the IUI can cause significant BER variations based on the assigned power for each user.The exact BER change that results by changing the modulation orders is given in Fig. 9, which is computed as the percentage of BER change relative to the identical modulation order case.For example, the BER percentage of change for U 1 with m = [8,4] is calculated as (P B1 |m = [8,4]−P B1 |m = [8,8])/P B1 |m = [8,8].However, for U 2 with m = [8,4], the percentage of change is calculated as (P B2 |m = [8,4] − P B2 |m = [4,4])/P B2 |m = [4,4].As can be noted from the figure, the percentage of change converges to a constant value at high SNR values.In addition, the percentage of change for U 1 at high SNR values is within ±2% for all modulation orders, except for m = [8,64], which saturates at +30%.Similarly, the percentage of change for U 2 is within ±2%, except for m = [8,64], and it saturates at +50%.Fig. 10 shows the BER using various modulation orders for N = 3.Similar to the N = 2 case, the BER increase when higher modulation orders are used.However, the overall BER performance degrades when compared to N = 2 because the power budget is shared by three users, and additional interference is introduced by U 3 .Due to the equal power assignment for all modulation orders, U 1 is assigned a very low power coefficient, which results in poor BER performance.The power coefficient assigned to U 2 is higher than U 1 by 20 dB, however, the BER advantage is reduced by 6 dB, which corresponds to the large scale fading.The same observations apply to U 3 since it has 20 dB power advantage over U 2 .Nonetheless, the power advantage is reduced by the 6 dB relative pathloss difference.The figure also shows the impact of changing the modulation order for certain users on the BER of other users.More specifically, it can be noted that changing the modulation orders for any two users will have a negligible effect on the BER of the other users.
Fig. 11 quantifies the BER variation due to the modulation orders change.As can be noted from the figure, the change at high SNR values is roughly within ±2% for all users.At low and moderate SNR values, the change roughly less than ±5%.
Fig. 12 presents the BER for the case of N = 4.The power coefficients used are α 1 = 10 −6 , α 2 = 10 −4 , α 3 = 10 −2 and α 4 = 0.989899.Consequently, such power coefficients require extremely accurate power control at the base station.Similar to the N = 2 and 3 values, the BER for a particular user is roughly independent of other users' modulation orders.As can be noted from the figure, the small power coefficients, IUI and the large scale fading lead to a degraded BER when compared to cases with a lower number of users.As depicted in Fig. 13, the impact of changing the modulation orders of other users on a particular user's BER generally follows the other considered cases, where the change of the BER with respect to the identical modulation order is bounded by ±2%, at high SNR values for all users and modulation orders, except for m = [4,2,16,8] for U 2 , which saturates in the range of +10%.At low and moderate SNR values, the change is within ±5% for all users and modulation orders, except for m = [4, 2, 16, 8] for U 2 , where it is about +19%.
Fig. 14 shows the BER for N = 2, 3, . . ., 7, where QPSK modulation is adopted for all users.The power coefficients for each case are given by: • The power coefficients are selected to minimize the system's average BER given that SNR = 80 dB.Therefore, the system's average BER can not be considered minimum at low SNR values.Nevertheless, by noting that the power coefficients remain approximately unchanged for a wide range of SNR values, then the presented average BER can be considered near-optimum.Tables III and IV present the optimal power coefficients that minimize the system's average BER for various modulation orders and SNR conditions for N = 2 and 3, respectively.It can be noted from the selected power coefficients that the difference between α 1 and α N becomes significant as N increases.For example, the power difference between α 1 and α 7 is about 42.9 dB.Such a high power difference is due to the necessity of compensating the high attenuation caused by the large scale fading.Moreover, it is noted from Fig. 14 that the BER for all users and values of N approaches the system's average BER at high SNR values.It is also worth noting the trade-off between N and the BER, where adding one additional user to the network causes degradation in the system's average BER between 6.65 and 8.24 dB at BER of 10 −3 .As can be noted from the figure, the analytical and simulation results match very well for all the considered scenarios.Moreover, the figure presents the asymptotic BER, which approaches the exact BER high SNR values.As compared to SISO, the diversity gain for L n = 2 is about 12.3 and 10.2 dB for U 1 and U 2 , respectively, at BER of 10 −3 .The gain increases to about 18.9 and 15.6 dB, respectively for U 1 and U 2 , for L n = 4.

VI. CONCLUSIONS
This paper has derived asymptotic and exact analytical BER expressions for NOMA over i.i.d Rayleigh flat fading channels.The derived expressions are applicable for any number of users, where each user has an arbitrary modulation order including BPSK and rectangular/square QAM.The results were corroborated via Monte Carlo simulation results.The derived expressions were used to provide insights about the BER performance in various conditions and system configurations, including some extreme scenarios in terms of the number of users and modulation orders.For example, the BER results revealed that the BER of all users converge to the system's average BER at high SNR values when the optimal power coefficients are adopted.Moreover, when the power coefficients are roughly in the middle of the PCBs, the BER of each user becomes almost independent of the modulation orders of other users, which might be necessary for adaptive modulation.The closed-form PCBs were derived for the N = 2 and 3 cases with arbitrary modulation orders.These expressions were utilised as linear and non-linear constraints to compute the optimal power assignment that minimizes the system's average BER.Interestingly, using high modulation orders and large number of users make the power control process very critical, where extremely fine power tuning is required.Fig. 14.BER using optimal power coefficients for various number of users, where Mn = 4 ∀n.

Fig. 2
Fig.2shows an example for N = 3, where m =[4,4,2].In Fig.2a, the constellation diagram for each user is shown separately and the x SC real and imaginary components are annotated accordingly.It is worth noting that the constellation

Fig. 2 .
Fig. 2. Constellation points of: (a) All users without superposition coding.(b) U 2 and the superposition coding of U 1 and U 2 .(c) U 3 and the superposition coding of U 1 , U 2 and U 3 .
are labeled using integer numbers that represent the symbol values, which are obtained by converting the binary bits of each constellation point into an integer using linear mapping.The NOMA bit-word is denoted by b = [b 1 , b 2 , . . ., b q ], where q = N i=1 M i and b 1 is the most significant bit (MSB).The nth individual user bits can be expressed as b n = [b On , b On+1 , . . ., b On+Mn−1 ], where b On is user's MSB, and

Fig. 9 .
Fig. 9. BER percentage of change with respect to the identical modulation case for N = 2.

Fig. 15
Fig. 15 considers the BER with receiver diversity.The results are obtained for N = 2 considering L n = 1, 2 and 4 receiving antennas, and modulation orders of m = [64, 64].As can be noted from the figure, the analytical and simulation results match very well for all the considered scenarios.Moreover, the figure presents the asymptotic BER, which approaches the exact BER high SNR values.As compared to SISO, the diversity gain for L n = 2 is about 12.3 and 10.2 dB for U 1 and U 2 , respectively, at BER of 10 −3 .The gain

Fig. 11 .
Fig. 11.BER percentage of change with respect to the identical modulation case for N = 3.

Fig. 13 .
Fig. 13.BER percentage of change with respect to the identical modulation case for N = 4.

TABLE II THE
PCBS FOR SELECTED MODULATION ORDERS,

TABLE III OPTIMAL
POWER ASSIGNMENT INCLUDING THE SYSTEM'S AVERAGE BER FOR DIFFERENT m AND SNR, N = 2