where \({\bf R}_t\)  and \({\bf R}_r\)  denote transmit- and receive-side covariance matrices, respectively, the elements of matrix \({\bf H}_w\)  are independent and identically distributed (i.i.d.) complex zero-mean unit-variance Gaussian random variables. Thus, the distribution of \({\bf H}_k\)  follows as \({\rm vec}\left({\bf H}_k\right) \sim {\cal C N}({\bf 0}_{N_tN_r}, {{\bf R}_t}^{\rm T}\otimes {\bf R}_{{{r} }})\) . It is worthwhile to point out that \({\bf R}_t\)  and \({\bf R}_r\)  are positive semi-definite Hermitian matrices. The validity of Kronecker model has been corroborated regardless of antenna configurations and intra-array spacing, provided that the transmit and receive angular power profiles are separable [1]}. Besides, in order to conserve energy, relieve signaling overhead and avoid frequent CSI feedback in large-scale energy- and bandwidth-limited IoT networks, we assume that only statistical CSI is available to the BS [2]}, [3]}.
Since NOMA is employed to simultaneously accommodate multiple clusters by sharing the same physical resources, the information bearing vector of \(M\)  data streams can be constructed by using superposition coding as follows
\({\bf {s}} = \sqrt{P} \sum \limits _{k = 1}^K {{\rm diag}\left\lbrace \sqrt{\zeta _{1,k}},\cdots ,\sqrt{\zeta _{M,k}}\right\rbrace {{\bf {x}}_k}} ,\) 
