A Statistical Estimation of 5G Massive MIMO Networks’ Exposure Using Stochastic Geometry in mmWave Bands

: This paper aims to derive an analytical modelling of the downlink exposure in 5G massive Multiple Input Multiple Output (MIMO) antenna networks using stochastic geometry. The Poisson point process (PPP) is assumed for base station (BS) distribution. The power received at the transmitter is modeled as a shot ‐ noise process with a modified power law. The distributions of 5G massive MIMO antenna gain and channel gain were obtained by fitting simulation results from the NYUSIM channel simulator. The fitted distributions, e.g., exponential and gamma distribution for antenna and channel gain respectively, were then implemented into an analytical framework. In this paper, we obtained the closed ‐ form expression of the moment ‐ generating function (MGF) for the total exposure in the network. The framework is then validated by numerical simulations. The sensitivity analysis is carried out to investigate the impact of key parameters, e.g., BS density, path loss exponent, and transmission probability. We then proved and quantified the significant impact the transmission probability on global exposure, which indicates the importance of considering the network usage in 5G exposure estimations.


Introduction
The explosive growth in communication devices in the last decade pushed the industry towards exploiting the existing communication spectrum. The recent mass-deployment of Internet of Things (IoT) devices paired with increasing demand for quality in the streaming and entertainment industries made satisfying these requirements using the same spectrum particularly challenging. To ensure that the huge spectrum capacity required to maintain a stable communication system is met, the industry started shifting with the design of 5G networks to millimeter-wave (mmWave) frequency bands. For indoor communication, mmWave bands are already in use, e.g., IEEE 802.11ad [1]. However, mmWave channels are fundamentally different from Ultra High Frequency channels used for existing cellular communication protocols, they have higher propagation loss in air, and low penetration for construction materials and foliage [2]. The high path loss in mmWave frequencies makes it necessary to implement methods to increase the received power at the mobile terminal (MT) in order to maintain an acceptable signal-to-noise ratio (SNR) to decode the signal. One way of increasing the SNR is by increasing the gain at the base station (BS). Therefore, beamforming became one of the essential technologies for 5G mmWave communication. The usage of a large number of radiating elements at the antennas will allow them to emit signals in fast-varying, high-gain beams focused on the desired MT, while simultaneously decreasing the interference to other MTs.
Exposure compliance assessments are performed by manufacturers and operators to ensure that emitting devices are compliant with safety regulations concerning human exposure to electromagnetic fields. Safety regulations concerning human exposure to EM fields have been specified by International Commission on Non-Ionizing Radiation Protection (ICNIRP) [3] and IEEE [4]. Multiple studies have been made over the years to characterize the exposure in current and legacy communication systems which are based on assuming the maximum theoretical power emitted for the antenna. However, these characterizations struggle to stay accurate for 5G and future technologies due to the fundamental differences between the systems. In [5,6], it has been shown that the actual power contributing to the electromagnetic field exposure is much lower than the theoretical maximum. This is due to the uncertainty in the values of the emitted and received power from a massive MIMO antenna deploying beamforming. The spatio-temporal variation of the antenna pattern, and its dependency on the MT distribution and channel characteristics make the transition to statistical methods for the exposure estimation a necessity if accurate estimations are desired. Statistical analysis of cellular networks' exposure has been a case study for the last couple of years, using simulations to model the behavior of the channel and the transmission e.g., in [7], but accurate analytical analysis is yet to be conducted.

Current Approaches in Exposure Estiamation
The assessment of the exposure is often realized either by in-situ measurements [8] or systemlevel simulations [6]. For 5G specifically, in-situ measurements are generally unconvincing in studying the exposure of a 5G system due to the variability of the transmission beams. The beam serving a specific area or MT changes in a period much smaller than the sweep time of most commercially available measuring equipment which means that during the time the equipment is measuring the whole band, the beam would have already changed. Moreover, the high dependence of the received power on the channel and the MT distribution make the analysis of the measurements especially challenging. Thus, it is why recent studies suggest the extrapolation of the received power from a single constant beam, ignoring the system usage and variability. The difficulty of generalizing measurements to study the exposure make theoretical studies and simulations more attractive. These studies allow freedom in assuming the network to emulate different scenarios. While system simulations are straightforward in design, they are very time and resource heavy, and they do not give direct insight on the different system parameters, therefore analytical analysis becomes more attractive to consider. Analytical estimations of the exposure have been previously considered, e.g., [5] although with a simplified antenna beamforming and system design.
An area in analytical modeling that was widely used in the domain of cellular networks is stochastic geometry. The assumption that the BSs are distributed following a point process gives this approach a relatively accurate basis and provides the mathematical tools necessary for the modeling of the system. Obtaining the closed-form expression will also give a direct way to study different parameters affecting the exposure, hence we chose it as the basis of our study. In this study, we also use the active antenna model developed by 3GPP [9,10], which assumes a uniform rectangular array antenna, for communication in bands above 6 GHz. As for the channel model, we use the NYUSIM statistical mmWave channel simulator [11], which is backed by extensive measurements in the mmWave band, and was shown to be better suited for above 6 GHz frequencies [12]. Multiple largescale channel measurements have been conducted in the mmWave band over the last few years [13][14][15][16], and multiple studies compared the different channel models used by the different standardization bodies as in [12,17].

Our Approach and Contributions
The mathematical approach to analyzing cellular systems is often considered as a promising solution [18]. However, for cellular networks, analytical tractability has been proven by relying on the methodology of stochastic geometry while modeling the BSs using point processes. Stochastic geometry has been widely applied in the field of communication theory, and it has proven that it is a powerful approach in conducting performance analysis and solving optimization problems. Examples of the usage of stochastic geometry in modeling cellular networks can be found in [19][20][21][22][23][24]. In this study, we model the locations of the base stations in the network as identically and independently distributed following a homogeneous Poisson point process (PPP).
In this paper, we fitted the channel gain and the antenna gain into statistical distributions based on channel simulations from NYUSIM and the 3GPP antenna array model. We determined the closed-form expression of the MGF for the exposure in the network, which is equivalent to the total power received at a certain MT location, with which we compute its cumulative distribution function (CDF). We analyzed the behavior of the 90th percentile of the exposure for different network parameters. We showed the comparison between this model and a simplified, less accurate, model with simple antenna and channel models in [25]. We performed a sensitivity analysis on the exposure expression to show the impact of each parameter on the global exposure. We also proved that ignoring the network usage when analyzing exposure will lead to large inaccuracies. The list of the symbols used in the paper is presented in Table 1 and the list of all abbreviation used can be found in in Appendix B at the end of the article.

System Model
The aim of this study is to determine the closed-form equation of the electromagnetic field exposure in a 5G massive MIMO cell emitting in the mmWave band. The base stations are modeled as points of a stationary and isotropic Poisson point process Ψ with density . The locations of said points are denoted by ∈ Ψ ⊆ ℝ . The MTs are distributed uniformly and independently in a bounded region of ℝ and thus, in addition to the assumption on the BS distribution, we can generalize the distribution of the exposure over the whole cell area and the framework can be developed for a typical MT denoted by MT0 located at the origin without losing generality. To investigate the effect of the network load on the exposure, we assume that the BSs emit with probability , i.e., at a certain time-frame the number of active MTs is less than the number of base stations , which means that the BSs will have an effective density of . Since line-of-sight (LOS) paths are the highest contributors to the exposure, we assume that all the paths between the BS and the MT are LOS paths. The BSs will emit with a constant power denoted by . We assume a 5G massive MIMO network where the BSs are composed of antenna arrays with identical structure, and thus identical radiation patterns. The BS is assumed as single user MIMO, and it emits towards one MT in a single time slot. We assume no downlink power control exist in the network, which seems to be the case for 5G networks [26], so every MT is being allocated the whole power resource for each transmission slot. We implement the channel model NYUSIM, developed by New York University Wireless Group [11] which divides the channel power into clusters and paths between the BS and the MT. The received power at each MT is dependent on the transmitted power, the propagation distance, the channel gain, and the antenna gain. The global exposure can then be defined as ∑ ∈ and it can be rewritten as where , , and L are the link's channel gain, antenna gain, and path loss, for the link between the MT to the ith BS, respectively. To better analyze the distribution of total received power , we focus on studying the CDF of total received power, which is defined as In the following subsections, we describe the system model, including the antenna array gain, the channel gain, and derive the expression to compute the CDF of the exposure of the assumed network.

Path Loss Model
In the current paper, a two-state path loss model is considered. It is modeled as a shot noise process with a modified power law [27]. We divide path loss into two regions, with and . The explicit expression for path loss for the link between UE and ith BS is denoted as , , where is the distance between and , and is the path loss exponent assumed constant in the network.
is an independently distributed random variable drawn from a common distribution and independent of . It accounts for the transmitted power, the channel gain and the antenna gain for the i th link, defined as . The details of and are presented in the following subsections.

Remark 1.
In the present paper, we adopt this two-state path loss model, which has already included the simpler case by setting 0 used in the authors' previous work [25]. In the region where , we approximate the path loss with a constant which can be interpreted as a rough average of the attenuation in the near-field region while for , the power law decay model is used.

Antenna Model
We implement the 3GPP active antenna model used to model the channels above 6 GHz [9]. The model can be used on frequencies up to 100 GHz. We have implemented the 3D model (2D steering) with steering in both elevation and azimuth since both scanning directions will be exploited in mmWave antennas [28]. The total antenna radiation pattern is the combination of the radiating element's radiation pattern and the array factor. The modeling of the antenna array is explained in the following subsections.

Element Pattern
By denoting and as the zenith and azimuth angles in the local coordinate system centered at the center of the antenna array, respectively, and defined as 0° 180° and 180° 180°, we can express the radiation pattern of an antenna element as where , and , are the vertical and horizontal radiation patterns of the array element, expressed as , min 12 where, and are the zenith and azimuth angles respectively in the local coordinate system centered at the center of the antenna array and defined as 0° 180° and 180° 180°. 65° are the 3 dB vertical and horizontal beamwidths of the antenna, 30 dB are the vertical and horizontal sidelobe attenuation levels, respectively, and 30 dB is the antenna's front-to-back ratio.

Array Pattern
Implementing a large number of antenna elements in the array will allow high-gain transmission of power towards the served MT by concentrating the field intensity in the direction of transmission [29]. The interference of the transmitted electromagnetic fields will form points of high intensity in certain directions, while producing low intensity transmissions towards unwanted directions. The array pattern is determined according to the physical implementation of the antenna elements and the electrical steering of the main beam. The array's radiation pattern is determined using the element radiation pattern in Equation (6), and assuming a uniformly spaced rectangular array (URA) in the horizontal and vertical dimensions with identical element radiation pattern as in where , is the phase shift due to the antenna element placement, , is the weighting factor which provides the electrical down-tilt and attenuation of the antenna's sidelobes, and and are the number of horizontal and vertical antenna elements respectively. As it is shown in Equation (9), is dependent on the antenna steering angles, the electrical down-tilt steering angle and the electrical horizontal steering angle , and since we assume that the main beam of the antenna is centered towards the receiving MT, the antenna array will have a different radiation pattern towards each MT depending on their location. The 2D slice of the antenna pattern on the azimuth plane for a 16 × 16 array can be presented as in Figure 1.

Channel Model
In our model, we separate the two gains obtained from the channel simulation into channel gain obtained directly from the NYU channel simulation, and antenna gain obtained by determining the antenna radiation pattern for every scanning angle after each channel simulation. Using the NYUSIM channel simulator, we determine, in each simulation, the LOS link characteristics between the serving antenna and the receiving MT. From that link, we can determine the channel gain, and the angles of departure (AoD) from the BS towards the MT. The AoDs are then used to determine the array gain using the 3GPP active array model [9], modeled in Section 2.1.2. The antenna gain for each subpath can be determined from this model, according to its AoD, as shown in the realization in Figure 1. We assume that the antennas are deployed in an urban microcell (UMI) scenario. The channel model is presented in Figure 2. The parameters used for the channel simulation are presented in Table 2. Running a large number of simulations, we empirically fitted the distributions of the gain components into statistical models. We chose the distributions for the channel and antenna gain by determining the best fit out of the common distributions using MATLAB. The distributions' parameters summarized in Table 3 were determined accordingly.

. Array Gain
Using the AoDs from the simulated paths in NYUSIM we have determined, for each path, the array gain towards each specific MT. The distributions of the array gain using the 3GPP active antenna array model, presented in Figure 3, has been fit into an exponential distribution with probability density function (PDF) defined as .

Channel Gain
Channel gain ℎ is the gain produced by the exploitation of the channel between the Tx and Rx such as the diversity gain, multipath, hardening, etc. [11]. We fitted the channel gain into a gamma distribution, as shown in Figure 4, with PDF given in Equation (11). The parameters for the fit gain distributions are summarized in Table 3.

Exposure Estimation
In this section, we derive the analytical framework based on the system model of the 5G network given in Section 2.
Since the distribution of is unknown, we adopt the Gil-Pelaez inversion theorem [30] to compute the CDF defined in Equation (2), which is based on characteristic function of random variable . The CDF is given by where Φ represents the characteristic function of , which can be derived from momentgenerating function as Φ . The moment generating function (MGF) of a random variable, defined as ∶ , ∀ ∈ ℝ, is an alternative representation of the random variable other than its probability distribution. It can be used to determine the distribution's moments as , ∀ ∈ ℝ . Since the expected value of the distribution is determined by derivation of the MGF in comparison to integrating the probability density, using the MGF becomes more attractive for complicated random variables thanks to the simpler operation. Following the path loss model defined in Equation (3), the total power at the can then be represented by: Where ℬ is a set of random variables such that ℬ 1 1 ℬ 0 . Since we divided the space into two regions, the total power received in the whole region can be written as where is the total power received from transmitters within distance from the MT, and is the total received power from transmitters at distances greater than . To account for transmission gain, we assume that is drawn from the distribution having the PDF: . ℎ . Where, G and H are the antenna gain and the channel gain distributions, respectively. As mentioned in 2.2, we fit the antenna gain into an exponential distribution with PDF , and the channel gain into a gamma distribution with PDF ℎ ℎ ⁄ . The MGF of this model is presented in Theorem 1.

Theorem 1.
The MGF of the exposure in the network modeled in Section 2, given the array gain and channel gain in (10) and (11) is formulated as: where and where, , Γ , is the upper incomplete gamma function, , is the lower incomplete gamma function [31], is the generalized hypergeometric function defined as ∑ ∏ ∏ ! , and Γ is the standard gamma function.
Proof See Appendix A Just by observing the exposure definition in Equation (1), we notice some of the effects the parameters will have on the total exposure. However, we are also interested in estimating the importance of each parameter on the total exposure and analyzing the way these parameters affect the total exposure. In the next section, we first verify our model with a Monte-Carlo simulation, then we investigate the variation of the 90 ℎ percentile of the exposure, we perform a sensitivity analysis to quantify the importance of each parameter on the total exposure, and then we compare the model we developed in this article with the old model developed in [25].

Numerical Results
In this section, we validate the model against numerical simulations, and we present some numerical results of the variation of the 90th percentile of the exposure against different network parameters.
In Figure 5, we validate the CDF from with given by Theorem 1, with Monte-Carlo simulations of the same scenario, to verify the model we developed while the parameters used in the validation are presented in Table 4.
.  The expression derived above has an indeterminate expectation, but from the determined CDF we can easily determine the 90th percentile of the exposure, which is a percentile that is often used to express the exposure in a network. We determine the 90th percentile as function of the base station density for different values of the path loss exponent , the results are presented in Figure 6. We also determined the variation of the 90th percentile of the exposure as a function of the system utilization in Figure 7. The parameters used to simulate the results of Figures 6 and 7 are presented in Table  5.   We also perform a comparison between the newly created model and the model in our previous work [25] versus the exposure simulation derived directly from the NYUSIM data in Figure 8. This comparison shows that the old model overestimates the exposure especially at the lower percentiles. This overestimation can be attributed to the fact that the old model assumes uniform array pattern, which is equivalent of having an isotropic antenna gain at the transmitter, and gain values towards the MT. It is also apparent, in the model derived in this study compared to the simulations, the small error between the simulation and our model which can be attributed to the error in fitting finite gain data into infinite distributions alongside the errors from the imperfect fitting. The former error can be calculated depending on the maximum gain the antenna can produce as in Equation (16) where is the probability density function of the fitted distribution, and is the maximum gain the antenna array can produce.
(16) Figure 8. Comparison between the current model developed in this paper, its verification using Monte-Carlo simulations assuming the gain distributions from Section 2.2, and the old model from [25], versus a Monte-Carlo simulation of the exposure using the gain simulated by NYUSIM.
Since massive MIMO will be implemented in really diverse scenarios, it is important to investigate the effect of each variable on the exposure. For this, we perform a variance-based sensitivity analysis on the calculated model which estimates the effect of the variance of the inputs on the variance of the output as form of indices called Sobol indices [32]. Sensitivity analysis is often performed to determine the importance of the input variables on the output of the model. nth degree Sobol indices , … determine the fraction of the output variance attributed to the set of inputs of degree . We are interested here in determining the total sobol indices which can be defined as in Equation (17). The total Sobol indices determine the variance of the output due to the variance of an input, in addition to the interactions between the specified input variable and the other input variables. It is an alternative to computing the higher-order sobol indices for every variable.
Since obtaining the exposure from the CDF requires solving an inverse function, it takes a long time to determine the sobol indices and it may introduce inaccuracies at extreme input values. To avoid this, we use polynomial chaos expantion (PCE) using Latin Hypercube samples to estimate a metamodel to represent the 90th percentile of the exposure in the cell. PCE approximates the relation between the model's output to its inputs by expanding it in an orthogonal polynomial basis [32]. The metamodel can be denoted as where Ψ are miltivariate polynomials defined as the product of univariate polynomials of degrees , … , and are the polynomial's coefficients. The fit model has a leave-one-out error of 0.14%. Using this metamodel we estimate the Sobol indices for each one of the variables and we obtain the results in Table 6. We can see that the sum of the total Sobol indices is greater than one, showing the interactions between the input variables in determining the percentile.

Discussion
In this section, we discuss the results of the numerical outcomes of our model. We discuss the effect of some of the network parameters considered throughout this study on the exposure, and we discuss the result of the sensitivity analysis considering the importance of each parameter on the global exposure.
From the results presented in Figure 5, we can show the accuracy of the developed model in estimating the total exposure in the cell. The result shows a good overlap between the analytical framework and the numerical validation. Here, the blue solid line is obtained by considering the fitted distribution of the channel and antenna gain, e.g., Equations (10) and (11). The purpose of Figure 5 is to verify the correctness of the analytical framework. It should be noted that the execution time in obtaining results from the analytical framework is much quicker than the numerical simulations. The sensitivity analysis shows the big effect the path loss exponent has in comparison to the other variables. In terms of exposure, this effect is desirable especially knowing that the path loss exponent is relatively high in mmWave channels leading to lower power in the cell. Although, since cellular networks are usually designed to maintain a sufficient SNR, it would give clearer insight analyzing the exposure in relation with the SNR. On the other hand, our study sheds the light on a usually ignored, but increasingly discussed, aspect of cellular network design and analysis which is the total electromagnetic wave power present in the cell area.
As previously mentioned in Section 1.1, the 5G NR architecture make it difficult to accurately measure the exposure because of its beam behavior. In terms of in-situ measurements of the exposure, they are being conducted without consideration of the variability in the system utilization in order to have constant beam behavior. This assumption will ignore the effect the system usage will have on the actual exposure in the cell. As we can see from our analysis, the transmission probability, which is directly related to the system utilization, it is the second-most important variable affecting the exposure since it accounts for 0.44 of the total exposure's variance, and simply measuring at constant full transmission may lead to major overestimation of the exposure.

Conclusions
In this study, we have determined a closed-form analytical representation of the downlink exposure of a 5G massive MIMO network distributed following a PPP with realistic transmission gain and channel representation using statistical distributions instead of approximated value. We have analyzed the distribution of the exposure for different implementation scenarios and we have shown the impact that the network characteristics have on the exposure in the cell. This approach allows the accurate study of the massive MIMO network without the need for costly simulations. We have also shown the significance of using a realistic antenna model as compared to a simple one. We also studied the significance of the key parameters in the network, and we showed the importance the network usage has on the total exposure and the importance of considering it when conducting exposure analysis.
This approach can be extended in future work to analyze different deployment scenarios of massive MIMO (e.g., cell-free, IoT … etc.). It can also benefit from realistic data representing the transmission gain after the deployment of the 5G massive MIMO antennas into the 5G network.
The 2 th ⁄ moment of can then be expressed by the close form equation by solving Equation (A14) knowing that . Since we assume that is constant, the integral can be solved to give Equation (A15) is the moment generating function of . The MGF can be determined by solving Equation (A16) using the PDF expression from Equation (A13).
Step (a) follows from a change of variables as giving us the Equation (A17). To determine the closed-form expression of the exponential integral in Equation (A18) we perform, in step (b), another change of variables as and we obtain in Equation (A20) the incomplete upper gamma function defined as Γ , ≔ . The expectation ⁄ Γ 1 2⁄ , can be determined using the law of the unconscious statistician [33] as in Equation (A21). ⁄ Γ 1 2⁄ , where, is the generalized hypergeometric function defined as : ∑ ∏ ∏ !