## Materials and Methods
### MIMO Channel Matrix, Beamforming and Precoding  
By \(s_{nk},n\in[1,2,\ldots,N],k\in[1,2,\ldots,K]\) we denote a set of rays between \(n\)th Tx and \(k\)th Rx points, with \(N\) and \(K\) being the total number of the Tx and Rx antennas, respectively. Then, the channel coefficient \(\hat{h}_{kn}\) is given by (see, e.g., in [17])  
\[\hat{h}_{kn}=\sum_{r\in s_{nk}}u_{r}\exp(-2\pi if\tau_{r}), \tag{1}\]  
where \(u_{r}\) is the voltage amplitude induced by the ray \(r\) at the Rx antenna terminal and \(\tau_{r}\) is the time-of-flight of the ray \(r\). The MIMO channel matrix \(\hat{\mathbf{H}}\) is obtained by evaluating (1) for each Tx-Rx pair in the simulation. Columns \(\hat{\mathbf{h}}_{k}\) of \(\hat{\mathbf{H}}\) are channel vectors to the Rx with index \(k\). The distance from different UEs to the BS may vary significantly. Therefore, UEs will experience differences in PL in comparison to one another. Thus, channel equalization is performed by normalizing \(\hat{\mathbf{H}}\) column-wise [2]  
\[\mathbf{H}=[\frac{\mathbf{h}_{1}}{||\mathbf{h}_{1}||},\frac{\mathbf{h}_{2}}{|| \mathbf{h}_{2}||},\ldots,\frac{\mathbf{h}_{K}}{||\mathbf{h}_{K}||}], \tag{2}\]  
where \(||\cdot||\) denotes the Frobenius norm. Having the normalized channel matrix, the Massive MIMO precoding matrices \(\mathbf{W}\) are given by [8]  
\[\mathbf{W}= \begin{cases}\alpha\mathbf{H}^{H},&\text{for MRT},\\ \alpha\mathbf{H}^{H}(\mathbf{H}\mathbf{H}^{H})^{-1},&\text{for ZF},\end{cases}  \tag{3}\]  
where \((\cdot)^{H}\) denotes the Hermitian transpose and \(\alpha\) is a real-valued normalization coefficient, chosen such that \(\mathbf{W}\) has unit Frobenius norm. Furthermore, the codebook steering matrix \(\mathbf{W}_{CB}\) is constructed from the steering column-vectors \(\mathbf{b}_{k}\in\mathbb{C}^{N}\) as [18]  
\[\mathbf{b}_{k}=[\exp(2\pi i(\mathbf{d}_{1},\mathbf{c}_{k})/ \lambda),\exp(2\pi i(\mathbf{d}_{2},\mathbf{c}_{k})/\lambda),\ldots,\exp(2\pi i (\mathbf{d}_{N},\mathbf{c}_{k})/\lambda)]^{T}, \tag{4}\] \[\mathbf{W}_{CB}=\alpha[\mathbf{b}_{1},\mathbf{b}_{2},\ldots, \mathbf{b}_{K}]. \tag{5}\]  
In (4), \(\mathbf{d}_{n}\) is a distance vector from the BS array center to the \(n\)th element, \(\mathbf{c}_{k}\) is a unit vector in the codebook direction assigned to the \(k\)th UE, and \((\cdot,\cdot)\) denotes the dot product. Knowing the channel vector \(\hat{\mathbf{h}}_{k}\), the beamforming direction \(\mathbf{c}_{k}\) is chosen as  
\[\mathbf{c}_{k}=\operatorname*{arg\,max}_{\{\mathbf{c}_{i}\}}[(\hat{\mathbf{h} }_{k},\mathbf{b}_{i})]. \tag{6}\]  
Figure 2: User equipment (UE) direction-of-departure (DOD) density averaged over 25 environment samples, viewed from the center of the BS array. The dashed line marks the cell limits. In (6), the maximization is carried out over the set of all beamforming directions \(\mathbf{c}_{i}\) supported by the BS. The beamforming direction vector \(\mathbf{c}_{i}\) is a unit vector in the direction \((\theta_{i},\varphi_{i})\) of the \(i\)th beam center  
\[\mathbf{c}_{i}=\mathbf{c}(\theta_{i},\varphi_{i}), \tag{7}\] \[\mathbf{c}(\theta,\varphi)=[\sin\theta\sin\varphi,\sin\theta\sin \varphi,\cos\theta]^{T}, \tag{8}\]  
where \(\theta\) and \(\varphi\) are the polar and azimuthal angles, respectively, in the spherical coordinate system depicted at Figure 1. The modeled system differentiates 32 beam directions in azimuth and 8 in elevation. The transmit vector \(\mathbf{t}\in\mathbb{C}^{N}\) is obtained by multiplying the precoding or steering matrix by the vector of transmitted symbols \(\mathbf{s}\)  
\[\mathbf{t}=\mathbf{Ws}. \tag{9}\]  
As the EMF is further assessed in terms of the time-average root mean square (RMS) values, with no loss of generality we set all transmitted symbols to be real-valued. In addition, we assume that no per-user power management is implemented at the BS and equal share of transmit power is directed to each UE. Therefore, we define \(\mathbf{s}=[\sqrt{1/K},\sqrt{1/K},\ldots,\sqrt{1/K}]^{T}\).