Massive MIMO Propagation Modeling With User-Induced Coupling Effects Using Ray-Tracing and FDTD

The article describes a numerical approach for Massive MIMO channel modeling that accounts for the effects of electromagnetic coupling between a user and the receiving device. The modeling is performed by a combination of the Finite-Difference Time-Domain and the Ray-Tracing methods, supplemented with a stochastic geometry model of the propagation environment. The influence of user-coupling on the channel properties was studied statistically using the singular value spread and matrix power ratio metrics of the channel correlation matrix. The time-averaged Poynting vector distribution in the near-field of the receiver was evaluated using a realistic human phantom model and the Maximum Ratio Transmission precoding scheme in the downlink. The average enhancement of the time-average Poynting vector magnitude at the receiver location, compared to the surrounding area, was found to be around 10 dB when using 36 antenna elements at the base station. The electromagnetic field exposure of the phantom was assessed in terms of the 10g-average peak-spatial Specific Absorption Rate and compared with the existing public guidelines. Comparison of the EMF and exposure results provides a new perspective on the future regulatory procedures.


I. INTRODUCTION
R ECENT years brought many advances to the development of next generation wireless networking. Significantly increased throughput, capacity and connection density are some of the requirements that 5G wireless technology must meet. Massive Multiple-Input Multiple-Output (MIMO) is one of the promising candidates to fulfill these requirements. Its operation is based on simultaneous utilization of a very large number of base-station (BS) antennas, (linear) transmission precoding and reception decoding schemes. An accurate channel model is needed to facilitate its design and deployment. Both deterministic [1] and probabilistic [2], [3] approaches to Massive MIMO channel modeling are currently being developed.
One of the key features of realistic Massive MIMO channels is the scattering cluster spatial non-stationarity observed at the BS side [4], i.e. channel variability across the BS antenna elements. Another important effect is the variability of the receiver's radiation pattern due to its orientation in space and the effects of the near-field coupling with the user body.
In [5], the cluster non-stationarity is modeled introducing cluster visibility regions. However, effects on the user side are left unaccounted.
This contribution, to the best of the authors' knowledge, for the first time describes the numerical approach that takes both aforementioned effects of cluster non-stationarity and user body coupling into account. These effects will be utilized by the Massive MIMO downlink precoding schemes to produce compact space regions of an elevated electromagnetic field (EMF) around the receiver antennas. Therefore, we apply the proposed approach to calculate the EMF distribution in proximity of a user on a small scale, and use it to estimate the human EMF-exposure. The same method can also be adapted for a 5G near-user antenna design.

II. METHODS
The proposed approach is based on a hybrid Ray-Tracing (RT) method which is used to model large-scale propagation and the Finite-Difference Time-Domain (FDTD) method, used for refinement of the RT results. In the following sections each step of the approach is discussed in detail and the connection between them is explained.
A. Ray-Tracing 1) Method Description: The wireless channel modeling was performed using a ray-launching [6] variant of the RT method [7]. A commercially available REMCOM Wireless InSite 3.2 software package was used [8]. A general RT procedure relies on a ray-optics approximation of the Maxwell equations. A transmitter is modeled by launching rays from its center in a finite set of directions, distributed over the complete  unit sphere. A ray is propagated through the environment predefined with its geometry and material properties, undergoing reflections, refractions, and diffractions, until its power reaches a predefined threshold. If a ray passes in the vicinity of a receiver, it contributes to the total field at that receiver's location.
2) Environment Model: The RT environment model is shown in Fig. 1. It consists of a rectangular room, cuboid scatterers, a massive MIMO BS and receiver locations positioned on a straight line. The room has dimensions of 40 m×20 m×5 m. Dielectric properties ε r = 7, σ = 1.5 · 10 −2 S/m were assigned to its floor, ceiling and walls to model concrete material. The BS was modeled with a 6-by-6 planar rectangular array of uniform 1-wavelength (λ) interelement spacing. The BS elements were excited individually with a continuous sinusoidal signal at f c = 3.5 GHz frequency (λ 86 mm). The array plane was normal to the x-axis and its center was set at at x = 7 m, y = 10 m, z = 4 m, as indicated in Fig. 1.
The 19 UE locations were arranged in a straight line along the x-axis separated by 1 m, in a range from 15 m to 33 m at the height of 1 m.
A perfect electrical conductor (PEC) cuboid of size 2 m×0.2 m×4 m (around 20λ×2λ×40λ at 3.5 GHz) was placed between the BS and the UEs (x = 10 m, y = 10 m), thus blocking the line-of-sight (LOS) propagation. Other PEC cuboid scatterers (2 m×0.5 m size and height from 2 m to 3 m) were randomly distributed along the perimeter of the room and had uniform random rotation along the vertical axis. This model generates realistic industrial environments of a fixed layout (e.g. a warehouse, assembly line) [9].
A 0.02 • ray-spacing was set for the ray-launching at the Tx. Each antenna element had vertical polarization and used an isotropic radiation pattern. Other RT solver parameters are summarized in Table I and were chosen as suggested in [10].
A single realization of the random scatterers produces an environment sample. By successively generating environment samples with random arrangement of scatterers and tracing rays in these samples, we obtained statistics of the EM-field that was incident at the NLOS locations in the described environment. The EM-field properties crucial for the operation of the massive MIMO (e.g. correlation between signals emitted from different antennas) are calculated based on a realistic geometry.
For each Tx-Rx pair (n, k) in the simulation, the RT solver returned a collection of rays {r} n,k . A ray holds properties of the propagation path it models, such as directionof-arrival (DOA), time-of-flight (TOF), Path Loss (PL), EMfield strength, etc. We used these properties to calculate the wireless channel in Section II-C.

B. EMF Coupling Effects
This section explains how the radiation pattern of the UE close to a phantom's head can be calculated using FDTD simulations.
In realistic scenarios, such as a mobile phone call, a user penetrates the near-field of his UE. This results in EMcoupling which affects the radiation pattern of the UE. We accounted for this by simulating the UE antenna in a usage mode together with a realistic human body model using FDTD method.
To model the UE antenna, we used a generic halfwavelength dipole antenna model, designed for the central frequency of 3.5 GHz. The overall dipole antenna length was 37 mm, arm diameter and the feed gap were 2 mm. The dipole arms were modeled as perfect electric conductors.
We used the ViP v.3.1 heterogeneous Duke human phantom [11] (further referred as the phantom) as a human model. The Sim4Life v4.4 (Zürich, Switzerland) software package was used for FDTD simulations throughout the paper. Fig. 2 shows the FDTD computational domain. The antenna was positioned at the center of the global coordinate system and oriented vertically. The phantom's head was positioned near the dipole, such that its bounding box center was coincident with the y-axis and the distance from the left ear to the dipole center was 20 mm. This setup aimed to reproduce the effects observed in a real-life phone call scenario, such as signal blockage by the user head.
The dipole was fed with a sinusoidal signal at 3.5 GHz and 1 W of total input power. The output of the simulation was a complex amplitude of the vertically polarized far-field E-field A(θ, φ).Â(θ, φ) was sampled on a surface of a 1 m radius sphere for elevation angle θ in range [0, π] and azimuth angle ϕ in range [0, 2π) with a 2 • step. Its ratio to the far-field E-field strength of an isotropic radiator fed the same input power ( 7.74 V/m) is the normalized far-field A(θ, φ), which will be used throughout the following sections.
The computational domain center was coincident with the center of the phantom's head bounding box and its dimensions were set to 256 mm×256 mm×250 mm. It fully enclosed the head and the antenna, at the same time significantly reducing the computational demands needed for simulations described in the next section. A maximum discrepancy of around 8% was observed in the antenna directivity if compared to the full-body simulation, which was considered acceptable, taking into account a simplified model of the UE.

C. Channel Matrix
This section explains how the channel matrix was calculated using rays at the Rx and the radiation pattern from the FDTD simulations.
First, a free-space channel matrix H fs was constructed from the collection of all rays {r} = ∪{r} n,k . A channel between the Tx antenna with index n and the Rx antenna with index k is found as Here s(n, k) are the indices of the rays in {r} n,k , E θ r and τ r are the vertical polarization component of the E-field and the TOF of the r th ray respectively.
Second, a reciprocity of an antenna in transmit-receive was utilized to introduce the DOA dependence into (1). As a ray holds the information about its DOA (θ r , ϕ r ), by weighting its contribution to the channel in (1) with the corresponding radiation pattern value A(θ r , ϕ r ), coupling effects are introduced into the channel matrix elements For an arbitrary incident direction (θ r , ϕ r ), A(θ r , ϕ r ) is calculated by the bilinear interpolation of A(θ, ϕ). Evaluating (3) for every Tx-Rx pair, we obtain a full massive MIMO channel matrix H nf .
We calculated the Singular Value Spread (SVS) κ(G) and Matrix Power Ratio (MPR) γ(G) of the channel correlation matrix G to quantify how the inclusion of the UE radiation pattern into the model affects the channel. The channel correlation matrix is obtained as where (·) H denotes the conjugate transpose of a matrix. The SVS of the channel correlation matrix is widely used in MIMO analysis [12]; it is a ratio between the largest and the smallest singular values of G.
The MPR is the ratio between the sum of the squared magnitudes of diagonal elements of G and all its elements [13] where g i,j are elements of G and K is the number of receivers.

D. EMF Distribution
To determine the EMF distribution in proximity of the user and the UE, FDTD simulation based on the RT results were performed. A ray with index r was modelled as a planewave source p r = (k r , a r , φ r ) that spans across the complete computational domain. It is described by its wave-vector k r , amplitude a r and phase delay at the domain center φ r .
To reduce the required computational resources, DOAs of the rays in r were substituted by the outer normal vectors of the faces of an icosahedral sphere (ico-sphere) {n ico (m)} [14], with m being the ico-sphere frequency. The wave-vector of a plane wave that corresponds to the i th face of the ico-sphere is then the inner normal of that face The amplitude and phase of the plane wave with index i were obtained by taking a sum of the complex amplitudes of the rays, DOAs of which have the smallest angular distance to the i th ico-sphere outer normal. Each ray's amplitude was weighted with the precoding matrix element w n,k , where (n, k) is the Tx-Rx pair indices for which the ray was calculated.
In this contribution we investigated the Maximum Ratio Transmission (MRT) precoding scheme [12]. MRT matrix W is proportional to the complex conjugate transpose of the channel matrix and its elements are given by where the normalization coefficient α is chosen such that W has unit Frobenius norm. The amplitude and phase of the plane wave with index i at the Rx with index k is then given by the magnitude and argument respectively of Here the outer sum is taken over the BS elements, with N being their overall number. Using (7) to calculate the precoding with (3) as the channel matrix, we obtained a set of the plane wave sources that models the EMF incidence at a massive MIMO user One FDTD simulation was required for each user in the environment. The result of a simulation was a 3-dimensional distribution of the electromagnetic field in free-space around the phantom and inside its tissues. Knowing the dielectric properties of the phantom's tissues we were able to calculate its EMF-exposure in terms of the Specific Absorption Rate (SAR). We used the peak-spatial SAR averaged over a 10-g cube (psSAR 10g ), calculated according to the IEEE/IEC 62704-1 standard [15]. psSAR 10g captures local exposure peaks, that highly-focused fields distinctive to the massive MIMO technology, are expected to produce.
International Comission on Non-Ionizing Radiation (ICNIRP) specifies basic restrictions for psSAR 10g [16] in the head for the general public (2 W/kg) and workers (10 W/kg). Based on these limits, reference levels on time-averaged power density are established using FDTD simulations with a single incident plane wave (10 W/m 2 and 50 W/m 2 for the general public and workers respectively). We will compare massive MIMO exposure and free-space power density with the ICNIRP guidelines in the following section.
For a more detailed description of the FDTD simulation setup, its sensitivity and error analysis, see [9].

III. RESULTS
In this section, the results are given in the order of the subsections of the previous section.

A. Far-Field Pattern
The normalized radiation pattern A(θ, ϕ) was calculated as described in Section II-B. Fig. 3 shows its magnitude and phase in spherical coordinates. Azimuth angle ϕ = π/2 corresponded to the incident radiation that is not obstructed by the user's head (see the coordinate system at Fig. 2).
A global maximum of |A| was observed around this azimuth angle in the horizontal plane (θ = π/2). This was expected, as the UE dipole is vertically oriented, it favours propagation in the horizontal plane. The phase response shown at the bottom of Fig. 3 is relatively flat around the magnitude maximum, varying for no more than π/4 in the π/2 neighbourhood of (θ = π/2, ϕ = π/2).
A global minimum of A magnitude was found near φ = 3π/2. This was also expected, as at these angles the incident radiation is attenuated by the user's head. The amount of attenuation is significant; it reaches around -18 dB (a factor of 63) if compared to the global maximum.
The phase of A oscillated rapidly around φ = 3π/2. This can be explained by a superposition of multiple diffraction paths that become dominant when the LOS is blocked by the head. However, this rapid phase variation does not make any noticeable contribution to the channel due to a low relative power of the propagation paths associated with it.

B. Channel Correlation Matrix
100 environment samples were simulated using the RT as described in Section II-A2.
First, a free-space channel correlation matrix G fs was calculated. This was done by evaluating (4) with H given by (1) for rays traced in each environment sample. At the top of Fig. 4 the magnitude of a sample average of G fs is shown.
Second, the channel correlation matrix G nf , that accounts for the UE radiation pattern was calculated. For each UE in the environment, the radiation pattern A(θ, ϕ) (see Section II-B) was rotated around the z-axis to an angle sampled randomly in [0, 2π). This models a random positioning of users with respect to the BS, when all user orientation directions were equally probable. The average G nf taken over all 100 environment samples is shown at the bottom of Fig. 4.
The diagonal elements of G were proportional to the power received in the downlink (DL) by users, if the MRT precoding was applied. Both G fs and G nf were diagonally dominant. This means that a significant portion of the DL transmitted power reached the intended receiver. However, in case of G fs there was an apparent increase in a relative magnitude of the super-and sub-diagonal elements. This corresponds to the increased interference of users with their closest neighbours in the simulated environment. Such effect was not observed in G nf . An average ratio between the elements on the super-and sub-diagonal and the diagonal elements is around 12% for G fs and less that 8% for G nf .
This difference (around 30% decrease) can be explained by an elevated correlation in the incident field as a function of DOA between closely separated locations. Presence of a highly-directional individual radiation patterns helps to resolve users that do not have enough spatial separation, which will have a positive effect on a massive MIMO system performance.
To study statistical properties of the channels, we calculated MPR and SVS for a different numbers of active UEs in the environment and Tx elements used at the BS. First, two and five UEs were selected randomly from the full simulated set and channel matrices (1) and (3) were calculated with only those UEs (index k) and a full 36-element array at the BS. Second, the same calculation was performed with five and ten randomly selected BS antenna elements and five randomly selected UEs. This was repeated 100 times in 100 environment samples, resulting in 10 4 samples per dataset. Additionally, κ(G) and γ(G) of the full channel (19 UEs) were calculated. Fig. 5(a) depicts histograms of κ(G) and γ(G) for 2, 5 and 19 UEs in black, blue and red, respectively. As the number of users in the channel increases, MPR (left column) decreases for both G fs (top) and G nf (bottom). This is in line with definition (5), as the number of the off-diagonal elements of a square matrix is roughly proportional to the square of the number of its diagonal elements. However, γ(G nf ) has a Higher γ(G) corresponds to channels with less correlation between UEs, being unity for purely orthogonal channels (with the correlation matrix G being diagonal). This means that the near-field user-UE coupling decreases (on average) correlation between massive MIMO users.
The right column in Fig. 5(a) shows κ(G) (in dB) for G fs (top) and G nf (bottom). For both G fs (top) and G nf larger user count yields larger SVS values. This was expected, as κ(G) is determined by only two most correlated rows of G. The more UEs are included in the channel (the greater dim(G) is), the higher chance there is for any two users to have correlated channels. Fig. 5(b) presents histograms of κ(G) and γ(G) for a fixed number of 5 UEs and a varying number of the BS antenna elements. With an increase of the element count, the mean of γ(G) monotonically grows for both G fs and G nf . The mean of γ(G nf ) remains larger than γ(G fs ) for all (equal) studies BS antenna counts.
Increasing the number of the BS antennas leads to a decrease of κ(G), as shown in the right column of Fig. 5(b). Moreover, the amount of κ(G) variation strongly depends on the ratio N/K between the number of the BS antennas and  the UEs. Increasing N from 5 to 10 (M/K from 1 to 2) leads to around 6 dB (4 times) decrease of mean values of both κ(G nf ) and κ(G fs ). Increasing N from 10 to 36 (N/K from 2 to 6.2) decreases the mean value of κ(G fs ) around 1 dB. This observation agrees well with the results of measurements in NLOS conditions [17].
Mean values of the histograms shown at Fig. 5 are compiled in Table II. In summary, the near-field coupling leads to the increase in the average of both γ(G) and κG.

C. Hot-Spot
In this section we present the results of the FDTD simulations described in Section III-C.
We performed FDTD simulations for all UEs in 10 environment samples. To evaluate the focusing of the EMF in proximity of the UE antenna we calculated the time-averaged Poynting vector at the receiver k as where electric and magnetic field vectors E k and H k were interpolated on a rectilinear grid. A horizontal slice coincident with the phantom's head center (z = 0) of S k , averaged over all 19 UEs in 10 environment samples is shown at Fig. 6. A strong focusing of S k is present near the dipole center, shown with a black circle in Fig. 6. This was the compound effect of the MRT precoding and the propagation environment. The precoding forced the signals from different BS elements to arrive coherently at the antenna terminal and a sufficiently To quantify the focusing effect we calculated the ratio between time-averaged power density at the antenna center r Rx = (0, 100, 0) (shown with a circle) and r sym = (0, −100, 0) that was symmetric to it with respect to the xz-plane (marked with a black cross at Fig. 6) As both the environment model and the phantom body were symmetric with respect to the xz-plane, it is reasonable to assume that for DL transmission with no precoding, the average of S k would be 1. Fig. 7 shows a histogram of S k (dB scale). Average S k is around 10 dB, and 5 th -95 th percentile range spans from just above 3 to slightly over 16 dB. This means that that the average hot-spot values of the power flux density are 10 times higher than in the surrounding space.

D. Specific Absorption Rate
The psSAR 10g was calculated in the phantom's head using built-in algorithms available in the Sim4Life software. We evaluated SAR 10g in the same environment samples and UE locations as were analyzed in Section III-C, resulting in 190 overall number of exposure samples. Fig. 8 depicts SAR 10g distribution in the horizontal slice at z = 0 and projections of the peak-cubes on that plane.
The largest portion of the peak-cubes ( 78%) were found near the left ear (see Fig. 8, red circle). These exposure peaks are produced by the EMF hot-spot at the UE antenna, located close to the ear.
In contrast to that, only around 9% of all peak-cubes is found at the right ear (indicated with the green circle at Fig. 8). These peak-cubes are spatially more concentrated than the cubes at the left ear, which suggests that they were mainly caused by the ear geometry. High curvature of the ear tends to concentrate EMF in the surrounding tissues even if exposed to a plane-wave [18]. This effect can be observed in Fig. 6 by noting the in-tissue local maximas of the power flux density. The most of the remaining peak-cubes were located at the back of the head, in the direction to the BS. In three samples the psSAR 10g was found in the eyes.
This allows to conclude that most of the Massive MIMO DL exposure was produced in the vicinity of the UE. To perform the exposure analysis, we introduced η -the psSAR 10g normalized to the time-averaged Poynting vector magnitude. Using the normalized exposure aids the analysis in two aspects. First, the influence of the BS transmit power and the PL is eliminated. Thus, values assessed at the UEs located at different distances to the BS can be correctly compared. Second, η is an indicator of how adequate the existing safety regulatory limits are. Current regulations require the EMF measured in free-space to comply with a predefined threshold. The threshold value was established based on the results of extensive direct measurements, simulations and additional safety factors, in which conventional single-antenna transmitters were used.
The definition of η depends on the location and conditions at which the EMF is assessed for normalization. Two normalization strategies were compared. First, η hs was calculated as psSAR 10g normalized to the maximum S k in the domain (hot-spot normalization). This is the time-averaged power flux density, that would occur close to the UE in operation.
Second, η pw was calculated by normalizing psSAR 10g to the average of the power flux densities in the incident planewave sources defined in (9). This is the value that would be observed in free-space when no user is present, as currently followed by the ICNIRP exposure guidelines [16].
As a reference we used η ref = 0.2 m 2 /kg, calculated from the ICNIRP general public restrictions (psSAR ICNIRP 10g = 2 W/kg, S ICNIRP = 10 W/m 2 ). Fig. 9 depicts mean values and 5 th -95 th percentiles of η hs and η pw calculated for 190 exposure samples, plotted versus the distance to the BS. For all studied distances η hs was nearly constant and around 3 times lower than η ref . It also had a relatively low variance: 90% of all samples fell between -3 and -7 dB. Additionally, psSAR 10g was highly correlated with the hot-spot power flux density (Pearson correlation coefficient ρ = 0.96). This was in agreement with the analysis of the peak-cube locations: most of the exposure was observed near the hot-spot. Hence, power flux density close to the Fig. 9.
psSAR 10g normalized to the time-averaged Poynting vector (dB scale). η hs (blue)-normalization using values at the hot-spot. ηpw (green)-normalization with the average of the incident plane waves. For all UE locations η hs is lower and ηpw is higher than the value calculated from the ICNIRP basic restrictions (used for the dB reference). receiver could be used to estimate the EMF-exposure. In addition, applying existing reference levels would overestimate psSAR 10g at least by a factor of 2.
In contrast, the average η pw is about 15 dB higher than η ref and has larger sample variation compared to η hs . This indicates that free space EMF measurements cannot reliably estimate the massive MIMO exposure. The existing free-space reference levels underestimate psSAR 10g by a factor ranging from 10 to 100.

IV. CONCLUSION AND FUTURE WORK
A novel numerical approach to massive MIMO channel modeling based on the RT and the FDTD methods was presented. For the first time, to the authors' knowledge, massive MIMO channels were simulated with the effects of near-field coupling between the receiver antenna and the user body. The importance of this effect was studied with models of a generic dipole antenna, a realistic human phantom and an industrial indoor environment. It was shown that including the coupling effects decreases correlation between closely spaced UEs by around 30%. The time-averaged Poynting vector magnitude enhancement at the UE was around 10 dB. psSAR 10g in the phantom's head was found to be directly proportional to the hot-spot power flux density. Normalized psSAR 10g complies with the ICNIRP guidelines, but assessing the reference levels in free space (with no active receiver) would lead to its underestimation by at least a factor of 10.
Although the presented study was carried out at a sub-6 GHz frequency, the same approach can be applied to a mm-Wave system analysis. The RT method is more accurate at higher frequencies with no performance loss. The FDTD memory grows as the third power, and time as the first power of frequency, which can limit the overall performance.
Future work will focus on studying the coupling between the BS antenna elements on the one hand and the user and the UE on the other hand, as these effects could have a significant impact on the channel [19].