Realistic Human Exposure at 3.5 GHz and 28 GHz for Distributed and Collocated MaMIMO in Indoor Environments Using Hybrid Ray-Tracing and FDTD

Realistic human downlink exposure at 3.5 and 28 GHz to electromagnetic fields is evaluated for distributed and collocated base stations using a hybrid ray-tracing/finite-difference time-domain method. For the first time, the absorbed power density is computed for distributed massive multiple-input multiple-output (DMaMIMO) 6G base stations at 28 GHz. The results are compared with 3.5 GHz 5G base stations. Computational costs are drastically increased at 28 GHz. A large analysis is realized by speed improvements and using two configurations. In the first, exposure distributions of DMaMIMO BS show clusters of low and high exposure. These clusters disappear when results are normalized with respect to the incoming power at the user. In the second, the influence of BS to user distance in line-of-sight (LOS) and non-line-of-sight (NLOS) scenarios shows expected results. This includes a power law relationship in LOS and shadowing in NLOS. The vast majority of exposure quantities are less than 4% of the limits of the International Commission for Non-Ionizing Radiation. Basic restrictions are respected when reference quantities are set to their limits. With equal power, distributed base stations contribute 2 to 3 times less to exposure than collocated base stations. Expressed as a ratio to their limits set by ICNIRP, the basic quantities are 5 to 10 dB lower than the reference quantities.


I. INTRODUCTION
The demand for more performant wireless networks is continuously increasing. Applications such as autonomous vehicles, Internet of Things (IoT) and Industry 4.0 require extremely reliable, ultra-low latency and very high throughput wireless connectivity [1]. The fifth generation (5G) of cellular networks has addressed these needs with a number of solutions [2], [3]. One such technology is Massive The associate editor coordinating the review of this manuscript and approving it for publication was Hassan Tariq Chattha .
Multiple-Input-Multiple-Output (MaMIMO) [4]: a large number of antennas in the transmitter (Tx) enables precoding (or beamforming [5]) to exploit the environment's spatial diversity, thereby increasing the network capacity [6]. Commercial MaMIMO deployment is ongoing [7,Sec. 5], [8]. This has prompted the research literature to focus on the question ''What is next?'', proposing Distributed MaMIMO proposed as a prospective sixth generation (6G) technology [9]. A MaMIMO base station (BS) is distributed (DMaMIMO) when its antenna elements are located in different geographic positions [10], i.e., the separation is much larger than the wavelength. When the separation is electrically small, channel capacity is hindered by unfavorable propagation conditions (reduced ''mutual orthogonality among the vector-valued channels to the terminals'' [11]) and pilot contamination [12]. When the separation is electrically large, the channel orthogonality increases with antenna density [13]. Also, the broader paradigm envisions a cell-free MaMIMO network of Access Points (AP), effectively removing pilot contamination [13], [14]. Another way 5G delivers on its ultra-low latency and exceptional data throughput is by increasing the operating frequency. The first frequency band (FR1) around 3.5 GHz is employed by the majority of commercially sold 5G BSs in Europe. The second frequency band (FR2) around 28 GHz is also known as the mm-Wave 5G band. In ideal conditions, the larger bandwidths enables record-breaking data rates [15], [16]. However, the higher free-space path loss limits the use of mm-Waves to short UEto-BS distance, such as near an AP. Therefore, DMaMIMO and high mm-Waves frequencies are two 6G directions that work constructively.
The assessment of human exposure to electromagnetic fields (EMFs) is of great importance for public health [17]. The International Commission on Non-Ionizing Radiation Protection (ICNIRP) provides internationally accepted safety guidelines that limit certain exposure metrics [18]. While worst-case exposure has been the subject of much research guaranteeing safety [19], [20], less is known about realistic exposure. To achieve this, the propagation and exposure steps need to be rigorously integrated into one method, for example by hybridizing ray tracing (RT) with the finitedifference time-domain (FDTD) method [21] or a spherical near-field transformation with FDTD [23]. Recently, the influence of DMaMIMO on the behavior of hotspots and their resulting exposure has been examined at 3.5 GHz [24]. The latest update of the ICNIRP guidelines in 2020 [18] introduces the absorbed power density S ab in the basic restriction above 6 GHz. The extension of realistic exposure of DMaMIMO BSs at 28 GHz is, to the best knowledge of the authors, an open question. This paper is first to provide answers in this unexplored research domain. The main novelty is proposing for the first time the hybrid RT/FDTD method at mm-Wave frequencies. Furthermore, the use of realistic state-of-the-art BSs is considered here. These enable two important novel results. In a realistic exposure scenario at 28 GHz, we present the first 1) comparison between distributed and collocated MaMIMO BS exposure; 2) quantification of the new basic quantity S ab for new MaMIMO technologies. We also provide the first comparison between these results at 3.5 GHz and 28 GHz.

A. CONFIGURATIONS
Figures 1a and 1b each show a realization of the two industrial environments studied. In both, a closed concrete (ε r = 7, σ = 1.5 × 10 −2 S/m) room containing metal cuboids models a Factory of the Future (FotF). The cuboids have a fixed footprint and variable height as in [21], and model, e.g., metal machinery. Their position is sampled from two spatial distributions covering the room, detailed below. These scatterers reflect and diffract the EMFs transmitted by the BS's antenna elements, which aims to replicate the propagation diversity of a realistic indoor environment. This collocated or distributed MaMIMO BS models a 5G BS or extremely large aperture array (ELAA), respectively. An array of patch antennas is suspended 1 and 2 m from the ceiling in configuration 1 and 2, respectively. They emit a signal at 3.5 GHz or 28 GHz (mm-Waves). This is shown in blue in Fig. 1. Isotropic receiver antennas (Rxs) receive this signal. Their height is 1.75 m, their polarization is vertical and their positions are on a straight line parallel to the room's walls. They model the user equipment (UE) held next to the head of humans (shown in red) walking in the FotF. Two pairs of channel types (configurations) are studied, namely (i) collocated and distributed MaMIMO channels and (ii) LOS and NLOS channels.

1) COLLOCATED (COL) AND DISTRIBUTED MaMIMO (DIS) CHANNELS
The room's footprint is square and has dimensions 50 × 50 × 10 m 3 . The distribution of the scatterers is as in [24], VOLUME 10, 2022 i.e., by removing randomly 30% of the scatterers on a uniform grid.
Three UEs receive a signal. The number of Rxs is chosen low to reduce computational costs. The number of channel realizations is chosen high (100) to enable a statistical analysis of the exposure distributions. The first UE is placed at a horizontal separation from the BS (henceforth UE separation) of 2.5 m. The second and third have a UE separation of 7.5 and 12.5 m.

2) LINE-OF-SIGHT (LOS) AND NON-LINE-OF-SIGHT (NLOS) CHANNELS
The room's footprint is rectangular with dimensions 40 × 20 × 5 m 3 . The distribution of the scatterers is as in [21], i.e., using the Poisson Disk sampling algorithm. Results are given with and without a 2 × 0.5 × 5 m 3 blocker. Nine UEs receive a signal to enable an analysis as function of receiver position. The first and last Rx have a UE separation of 6 and 22 m, respectively. The number of channel realizations is chosen low (5) to reduce computational costs.

3) PARAMETERS
The n77 and n261 frequency bands are used throughout the numerical pipeline. In the following, we refer to these as the 3.5 and 28 GHz bands, respectively. Their spectrum parameters are listed in Table 1. Two realistic BSs are modeled for 5G coverage at these frequencies. Their physical parameters are listed in Table 3. The parameters used for the RT simulations are shown in Table 2. These parameters are sufficient to accurately predict the channel in the modeled indoor environment [22]. The MaMIMO BS is an array consisting of equally spaced dual-polarized antenna elements. Their parameters are also listed in Table 3. The inter-element spacing in the collocated MaMIMO BS is 0.5λ for both the vertical and horizontal direction. In the distributed case, the inter-element spacing is 10 m, such that it covers the whole room.

B. RT-FDTD
The hybrid RT-FDTD tool developed in [21] is extended to mm-Wave frequencies. The fundamental method is the same: the scattered EMFs by plane waves (PWs) are first stored and then combined based on the full channel matrix. The schematic specifies how this fundamental method is looped over all parameters of interest, for one configuration. The RT and FDTD initialization followed by their hybridization are detailed schematically in Fig. 2.  For the FDTD initialization and exposure evaluation, a number of adaptations are made. These are shown visually in Fig. 3. A Multimodal Imaging-Based Detailed Anatomical  Model of the Human Head and Neck (MIDA) [25] is used as a phantom. The phantom is highly detailed and features a spatial isotropic resolution of 500 µm [26], shown in the upper front head. The model is finely voxelized due the small wavelength used, shown in the lower head. The head faces the BS, i.e., in the positive y 1 -direction or negative x 2 -direction for configuration 1 and 2, respectively. Maximum ratio transmission (MRT) precoding is used in a singleuser down-link (DL) transmission scenario. The incoming rays focus the z-component of the electric field to the location of a vertical dipole, representing the UE. A vertical and horizontal slice of the domain indicates by means of yellow and white colors that the volume around the focus point contains increased EMFs compared to the background. This hotspot is a result of focused rays interfering constructively. The focus point is located 2 cm behind the right ear, within the bounding box of the phantom. This location results in a worst-case exposure scenario because the most common region of maximum exposure is the skin of the ears [24]. The surface absorbed power density S ab averaged over a 4-cm 2 surface is shown on the back of the head. The indicated red colors reveal a region near the hotspot where it is increased.

C. COMPUTATIONAL CONSIDERATIONS
The MaMIMO BS requires the knowledge of the wireless channel to perform the precoding. We calculate the channel matrix coefficient as the sum of the signals with different directions of arrival (DoAs) [27]. A signal is the electric field's z-component received at the location of the UE.
The field is computed with an FDTD simulation of a unit PW incident on the human head phantom.
The higher frequency results in prohibitively large computational costs (due to, e.g., the Courant limit [28]). Therefore, the execution time of each FDTD simulation at 28 GHz is reduced by a combination of techniques. First, the workflow of Fig. 2, including pre-and post-processing, is parallelized on a high-performance computing cluster with 40 CPUs and a performant GPU. Second, the ray approximation is more justified at higher frequencies, reducing the number of required rays in the RT simulations. As a result, the number of reflections and diffractions is decreased, as shown in Table 2. Finally, at 27.9 GHz the millimeter-sized voxels overresolve the model's head features. This gives us the freedom to reduce the voxel size to λ/14 = 1.75 mm. This is not the case for 3.75 GHz, where the voxel size must be at least 1.75 mm = λ/46. The number of voxels at 3.5 and 28 GHz is 1.54 and 22.12 million, respectively. However, the computing and memory costs remain high. More importantly, the substantially increased memory load of the electromagnetic field data restricts the number of DoAs, channel realizations, and UE separations because, for each of these, the data is exported, added and weighted. Therefore, the number of FDTD simulations must be reduced. The number of DoAs is set to 20, which still enables the modeling of hotspots. The error is estimated in [21].
In configuration 1, the number of channel realizations is chosen large (100). Therefore, we study the exposure distributions in this configuration. In configuration 2, the number of receivers is chosen large (9). Thus, we study the influence of position in this configuration.

D. EXPOSURE ASSESSMENT
The reference levels defined by the ICNIRP for human exposure at both 3.5 and 28 GHz is the power flux density S inc in free-space. Here, a free-space computational domain is simulated without phantom. The MIDA phantom is added when simulating for the basic restrictions. The basic restriction at 3.5 GHz is given by the peak-spatial specific absorption rate averaged over a 10-g cube (psSAR 10g in W/kg) [18], [29]. At 28 GHz, the maximum surface absorbed power density S ab (in W/m 2 ) averaged over a 4-cm 2 surface is used to better account for the superficial nature of the exposure [18], [30]. We define η as the percentual ratio between each simulated exposure quantity γ (|ReS inc |, psSAR 10g or S ab ) and its maximum allowed value for the general public (GP) γ limit , as set by the ICNIRP guidelines: The values of γ limit are listed in Table 4. Note that all exposure quantities are local to the head and time-averaged. When the reference quantity (S inc ) is set to the reference level (10 W/m 2 ), the basic quantities (psSAR 10g and S ab ) should be below their exposure limit. Dividing these by η r is equivalent to setting the reference quantities to the reference levels. This results in the definition of α, where η b is the η of the basic restrictions. A value of α above 100% means that the basic restrictions would be exceeded if incident power densities are set to the reference levels. This is called hotspot normalization [31].

A. GENERAL RESULTS
In Figures 4 through 7, results from the first configuration (COL/DIS) are presented in a matrix of plots, comparing  Fig. 1a. The ''All rx'' line combines all 300 channel realizations into one CDF. A vertical line denotes the quantity that is greater than 95% of all these realizations' quantities. In Fig. 8 and 9, results from the second configuration compare the 3.5 and 28 GHz bands side by side. The main interest is how the exposure values depend on position to the BS and presence of a blocker. The left and right vertical axis indicates the relevant exposure quantity and its corresponding η value at 320 W BS output power. The horizontal axis reads the UE separation. Two lines are associated with the absence (LOS) or presence (NLOS, shown in Fig. 1b) of the blocker. Each result is represented by a dot. The number of channel realizations is 5, which estimates the mean and standard deviation of the exposure in LOS or NLOS. Lines connect the average of the 5 channel realizations across UE separation.

B. CONFIGURATION 1: COL VS DIS
In general, the Rxs closest to the BS have the highest exposure, in particular for the first Rx. However, the focus of this configuration's analysis is not on position, but on the distributions of exposure quantities. One hundred channel realizations were necessary to sufficiently resolve the exposure distributions.

1) REFERENCE QUANTITIES
The results of the reference quantities are shown in Figure 4. Realistic values are at most 4% of their limits in 95% of cases for a BS radiating 320 W. All exposure metrics are 2 to 3 times less with a distributed BS compared to a collocated BS. For the distributed BS, more elements are distant from the receivers, reducing the total incoming power. The exposure distribution with the collocated BS can be well modeled by a Ricean fading channel. This suggests that a LOS component or strong reflection path dominates the channel. The distributed BS does not generate a Ricean fading channel. Instead, two regions are identified of low and high exposure (Fig. 4, right), delimited around 0.19 mW/m 2 . We call these exposure clusters. In the first cluster, a similar exposure distribution to the collocated BS is seen. The standard deviation of the second cluster is 4.5 times broader than the first. These results can be explained as follows. When a UE is surrounded by scatterers, the probability that exposure belongs to the first cluster is highest. The channel resembles that of the collocated BS, because the interaction happens mainly by one or a few nearby antenna elements. When no nearby scatterers shadow the UE, the probability that exposure belongs to the second cluster is highest. The interaction can now occur with a larger number of antenna elements farther from the UE. Therefore, the number of reflections and diffractions taken by the rays is higher. The variability in exposure is thus higher. As seen in Fig. 1a, the scatterers belong to a regular grid, the antenna elements are directly above this grid, and the receivers are located at the center of each grid cell. This set-up makes it possible for receivers to be completely encircled or not by scatterers. Similar exposure distributions are observed when comparing the 3.5 and 28 GHz bands. With the collocated BS, the 95 th percentile is 12% lower at mm-Wave frequencies. With the distributed BS, the reduction is 55%. In this case, the standard deviation of the second exposure cluster is much lower. This is likely due to the increased path loss at these frequencies reducing the number of interactions with distant antenna elements.

2) BASIC QUANTITIES
The results of the basic quantities are shown in Fig. 5. Realistic values are at most 0.85% of their limits in 95% of cases. Note that these η values are generally 5 to 10 dB lower compared to those of the reference quantities. As a result, the α value is less than one and the basic quantities stay 5 to 10 dB below their limits when the reference quantities are set to their limits. The same conclusions as above can be drawn about the values and exposure distributions. The incoming electric fields are identical because the propagation of the rays onto the phantom is the same (the propagation step). Differences between Fig. 4 and 5 can only be attributed to the exposure metric used and the presence of the phantom (the exposure step, shown in Fig. 3). Because similar results are seen, the exposure step is mostly independent of the propagation step. However, in the 28 GHz band with the distributed BS, three exposure clusters are seen, delimited around 0.05 and 0.15 mW/m 2 . This deviation can only be explained by the interplay of the propagation and exposure step. Because high exposure clusters are more prominent in NLOS propagation, the formation of hotspots is more likely in the high exposure clusters. When the location of the hotspot is within a certain distance of the skin, the hotspot causes increased exposure. This is registered only by the superficial S ab exposure metric on the phantom's skin. The aforementioned α value is highest in this case because of this effect.

3) NORMALIZED RESULTS TO INCOMING RX POWER
Before the precoding step, the total power is split evenly among the antenna elements. In the distributed MaMIMO BS case, the receiver is surrounded by a small number of antenna elements each having a fraction of the total output power. In the collocated MaMIMO BS case, the receiver is near all antenna elements radiating the total output power. The exposure quantities are therefore substantially lower with distributed BSs compared to collocated ones. Commercial distributed MaMIMO BSs do not exist yet; neither do corresponding regulations. It is therefore unknown what realistic total output powers will be of these systems. The output power is presumably higher, and likely dependent on the number of antenna elements and the area covered. Therefore, the above conclusion may not be realistic due to the equal power assumption. There are a number of ways to normalize the exposure quantities such that a fair comparison becomes possible. One of them is dividing the result by the total incoming power P Rx on the Rx. The total electric field is where i denotes a DoA angle, w i is the corresponding weight from the precoding matrix andÊ i is the electric field for a unit PW with DoA angle i. The total incoming power is defined as The results are now normalized w.r.t. both 1 W total power from the Tx and 1 W incoming power on the Rx. The simulated exposure value γ is retrieved as follows: γ = P BS · P Rx γ 1 W 1 W . In this equation, P BS is the BS's total emitted power and γ 1 W 1 W is the doubly normalized exposure quantity. The latter is shown in Figures 6 and 7 for the reference and basic quantities, respectively. As the propagation step is effectively removed, there is no statistically significant dependence on VOLUME 10, 2022 receiver position for the reference quantities. For the basic quantities, there is a pronounced dependence on position in the mm-Wave distributed case. This validates the previous conclusion that a strong skin coupling with the hotspots in the exposure step is present. For Rx 2 and 3, two clusters are seen delimited by approximately 5 µW/m 2 . This causes the three clusters observed earlier, which were only present for Rx 2 and 3. In all other cases, exposure distributions are well modeled by a normal distribution.

C. CONFIGURATION 2: LOS VS NLOS
Figures 8a and 8b compare the reference quantities as a function of distance at 3.5 and 28 GHz, respectively. Figure 9 does so for the basic quantities. All distance profiles are highly correlated across frequency and exposure quantity. The maximum η values in LOS are highest for the reference quantities. They are lowest for the basic quantities in the FIGURE 8. Reference quantities normalized for a 1 W radiating BS in configuration 2. Clarification of this figure is provided in Section III-A. The exposure is closer to the limit for 3.5 GHz than for 28 GHz. 28 GHz band. Therefore we find again that the α values are lower than 1. The results for the 3.5 GHz band are similar to those found in [31]. Differences are attributed to the different configurations used. We thus find the reference levels to be conservative, guaranteeing basic quantities to be within their limits. When comparing the LOS case with the NLOS case, expected results are obtained. First, the LOS exposure quantities decrease according to a power law; the NLOS exposure quantities do not. Second, the standard deviation is lower for LOS than in the NLOS case. Both of these can be explained by the free-space LOS rays dominating the channel. The results' dependence on the environment is then lower for LOS than in the NLOS case. In the NLOS case, shadowing causes the exposure quantities to be highest 12 to 17 meters separated from the BS. This is a consequence of the NLOS blocker seen in Fig. 1b, separated 5 meters from the BS.

IV. CONCLUSION
The 5G and 6G advances in mm-Wave frequencies and distributed MaMIMO demand realistic exposure assesment. A state-of-the-art hybrid RT-FDTD tool is extended to 28 GHz by including the new absorbed power density as basic quantity set forth by the ICNIRP 2020 guidelines [18]. A multi-dimensional comparison is enabled by increasing the computational efficiency and budget. In a first configuration, clusters are observed in the distribution of exposure quantities for distributed MaMIMO BSs. We discuss the need to normalize the exposure to the incoming power on the Rx. It is shown that exposure clusters are mainly caused by the propagation step. The highest 95 th percentile among the exposure quantities expressed w.r.t. to their maximum allowed by ICNIRP is found for collocated BSs at 3.5 GHz, at 4%. With equal power, distributed BSs contribute 2 to 3 times less to all exposure metric than collocated BSs. Exposure is lower at 28 GHz compared to 3.5 GHz, due to the increased path loss. In a second configuration, the influence of UE separation and presence of a NLOS blocker is analyzed. It is found that LOS exposure follows a power law and NLOS exposure shadows the UE. In both configurations, with respect to their limits, basic quantities are 5 to 10 dB lower than reference quantities, guaranteeing ICNIRP's assumption. In future work, more comparisons with different important factors could be considered. For example, the influence of using a different precoding technique on realistic exposure with distributed MaMIMO at mm-Waves is unknown. Furthermore, the realism of industrial indoor environments could be improved. LIDAR-based models are seen as a good candidate due to their high accuracy. Lastly, the exposure of other 6G technologies such as ELAAs and holographic MaMIMO could be investigated, ideally at sub-THz frequencies.