On the Correlation Between Incident Power Density and Temperature Increase for Exposures at Frequencies Above 6 GHz

International guidelines/standards for human exposure to electromagnetic fields have recently been revised to update the dosimetric reference limits (DRLs or basic restrictions) and exposure reference levels (ERLs), specifically for frequencies above 6 GHz. At such frequencies, the ERL is defined in terms of incident power density (IPD) and used as a practical quantity to assess compliance with DRLs (absorbed or epithelial power density) and therefore appropriately limits temperature elevation at the body surface. In the exposure standards, IPD is spatially averaged over an area of 4 cm2 below 30 GHz and 1 cm2 above 30 GHz, however the definition of IPD is given in a theoretical manner. With the progress in the development of product safety compliance assessment standards, one concern has been how to define the IPD considering practical measurement procedures. Two definitions or averaging methods were considered: using IPD vectors normal to the averaging surface and using magnitude (norm) of IPD vectors. As the exposure guidelines are intended to prevent excessive tissue heating, statistical analysis was therefore undertaken to investigate which IPD metric better correlates with the temperature increase. To this end, a large data set for several exposure scenarios was collected by different research institutions. The analysis of the obtained results is presented and shows that both definitions have high correlation with temperature rise, with slightly better correlation (0.9 vs. 0.8) for the definition using the magnitude of IPD vectors.


I. INTRODUCTION
With the advent of fifth-generation (5G) mobile devices, assessment of electromagnetic (EM) field exposure at frequencies above 6 GHz is receiving much consideration [1]. At these frequencies, often referred as millimeter-waves (mm-Waves), the established adverse effect on biological tissues is of thermal nature, related to a temperature increase of the superficial tissues, mainly skin and eye [2]- [5].
The associate editor coordinating the review of this manuscript and approving it for publication was Dost Muhammad Khan .
International exposure guidelines/standards for human protection from high-frequency electromagnetic fields have been recently revised by the International Commission on Non-Ionizing Radiation Protection (ICNIRP) [6] and the Institute of Electrical and Electronics Engineers (IEEE), with the International Commission on Electromagnetic Safety (ICES) under the Technical Committee (TC) 95 [7].
The use of the absorbed/epithelial power density (A/E-PD) as the basic restriction (BR) [6] or dosimetric reference limit (DRL) [7] for the frequencies above 6 GHz is one of the main changes in the guidelines/standards. Furthermore, the incident power density (IPD) in free space is defined as the reference level (RL) [6] or exposure reference level (ERL) [7]. It should be noted that at mm-Waves, IPD is a more practical and straightforward measurand than the A/E-PD.
Based on the exposure guidelines/standards, the IPD should be averaged over a squared area of 4 cm 2 for frequencies from 6 to 300 GHz. However, for frequencies higher than 30 GHz, additional criteria of IPD averaged over 1 cm 2 are given with a relaxation of RL/ERL by a factor of 2 for local beam-like exposures [6], [7]. Since no technical specifications have been provided on how to determine these spatial averaging, some product compliance assessment standards have been established by the International Electrotechnical Commission (IEC) TC106 and IEEE ICES TC34 committees. These dual-logo standards, recently published, provide experimental and numerical procedures to assess exposure to IPD from 6 GHz to 300 GHz [8], [9].
Beside the averaging technique, one of the goals in the development of these assessment standards has been the definition of a robust IPD metric considering practical aspects of laboratory measurement procedures. Two definitions have been proposed so far: • the first definition, where only normal components of the Poynting vector crossing the surface are used (sPD n ); • the second definition, where all total components (i.e., magnitude) of the Poynting vector are considered (sPD tot ). More recently, a third definition of IPD (i.e., including the reactive components of the Poynting vector) was proposed in [10], without however considering its relationship with tissue temperature rise. Since the exposure safety guidelines/ standards are established to limit tissue heating, the main criterion for evaluating suitability of a particular IPD metric for compliance should be based on the correlation of that metric with temperature rise and therefore the latter definition has not been hereby considered. To summarize, the above two definitions were used to derive different IPD quantities for compliance by spatial averaging over 1 cm 2 as well as 4 cm 2 of the exposure evaluation surface, thus producing four different metrics.
The effect of the IPD definition as well as the IPD averaging method on the relationship between the resulting compliance quantity and the temperature rise have been studied in several publications [11]- [23]. However, the influence of these IPD definitions on the temperature correlation was not considered. Moreover, in the near-field and for oblique incidence, these definitions lead to different values with significant differences in some exposure conditions [24]- [26]. Therefore, it becomes necessary to study both definitions in terms of correlation with surface temperature elevation.
A new working group (WG) 5 under Subcommittee (SC) 6 of IEEE ICES TC95 has been established to clarify these aspects. A large amount of data from different independent institutions were collected to study the IPD metric correlation with the surface temperature elevation in the frequency range from 6 to 300 GHz, which is essential for compliance assessment of new 5G wireless products. This paper shows the results of the statistical analysis performed by the WG on the collected data, as well as some discussions and conclusions.

II. DATA COLLECTION
Six different organizations collaborated to perform the intercomparison study [27]: the National Institute of Information and Communications Technology (NICT), Nagoya Institute of Technology (NITech), South China Agricultural University (SCAU), Dassault Systèmes SIMULIA (3DS), Foundation for Research on Information Technologies in Society (IT'IS), and the University of Split (UniSplit). The evaluation of the results started with the data collection and its structured aggregation. The data consisted of a total of 227 samples, obtained by combining the results from each participating institution who performed computer simulations of different exposure scenarios. Those simulations included the stratified tissue block models with various EM sources like dipole antennas, patch antennas and different antenna arrays at frequencies f = 10, 30, 60 and 90 GHz as well as a number of separation distances d between the source and body model set to 2, 5, 10, 50 and 150 mm. Each sample, corresponding to one particular simulation scenario, contained the computed peak values of spatially averaged PD (i.e., psPD n and psPD tot averaged over 1 cm 2 and 4 cm 2 ), the peak temperature increase p T , as well as the corresponding heating factors (HFs). The latter is defined as in [28]: All participants of the study have provided the results for the dipole antennas, while results for EM sources made of directional antennas were provided by some of the groups. Therefore, the statistical analysis has been performed on two macro-groups or data sets made of: • 115 samples for the dipole sources at all simulated distances and frequencies, referred to as dipole data; • 112 samples for directional antennas sources at all simulated distances and frequencies, referred to as directional antennas data. Among the latter, 46 samples come from dipole or patch array antennas with a beam shift in azimuth or elevation.

A. SCATTER PLOTS
The first analysis has been conducted to visualize the scatter plots of the peak spatial-averaged power densities, psPD, and the respective relative values of peak temperature increase, p T . Specifically, the following pairs of data series have VOLUME 10, 2022 been analyzed, first for dipole data, then for directional antennas data: • (psPD n avg.1, p T ); • (psPD tot avg.1, p T ); • (psPD n avg.4, p T ); • (psPD tot avg.4, p T ). To evaluate possible influence of near-field exposure conditions, a distinction has been made by excluding the data coming from the simulations with distance d = 2 mm (i.e., d ≥ 5 mm) and the data from all other simulation scenarios, i.e. d ≥ 2 mm. The scatter plots were aimed at visually establishing the type of correlation (linear or not) and its strength (i.e., strong, weak or no correlation).

B. STANDARD DEVIATIONS OF HF
After generating the scatter plots for all data series, the HFs have been analyzed grouping the results by institution/ group. To this end, the average value (µ HF ) and standard deviation (σ HF ) for every group have been computed for each of the four metrics, namely HF n avg.1, HF tot avg.1, HF n avg.4 and HF tot avg.4. The same computations were also performed considering the data from all groups together. The standard deviation of each HFsample was of particular interest, since a smaller value of σ HF would suggest stronger correlation with T . Finally, for each PD metric, the histograms of the HFdata (computed for all the data sets) have been plotted to examine their distribution around the average value µ HF .

C. PEARSON CORRELATION COEFFICIENTS
The analysis of the Pearson correlation coefficients has been performed to identify the PD metric that best correlates with the temperature increase T . Specifically, the following Pearson correlation coefficients were evaluated: • r n1 =Pearson coefficient of the pair (psPD n avg.1, p T ); • r tot1 =Pearson coefficient of the pair (psPD tot avg.1,p T ); • r n4 =Pearson coefficient of the pair (psPD n avg.4, p T ); • r tot4 =Pearson coefficient of the pair (psPD tot avg.4,p T ). Equation (2) was employed for the computation of the correlation coefficients [29]: where Cov(psPD j avg.k, p T ) is the covariance of psPD j avg.k and p T , while σ (psPD j avg.k) and σ (p T ) are the standard deviations of psPD j avg.k and p T , respectively. Each set of four different r jk was computed considering the distinction between data obtained without including the separation distance d = 2 mm and all other data, namely: • r n1 , r tot1 , r n4 , r tot4 for dipole and d ≥ 5 mm; • r n1 , r tot1 , r n4 , r tot4 for all dipole data; • r n1 ,r tot1 ,r n4 , r tot4 for directional antennas and d ≥ 5 mm; • r n1 , r tot1 , r n4 , r tot4 for all directional antennas data. Each of the four analyses above has been performed on data groups distinguished by the working frequency of the sources. Specifically, correlation coefficients have been computed first for all the available frequencies (f = 10, 30, 60, 90 GHz), then for the f = 10, 30 GHz set and finally for the f = 60, 90 GHz set, still distinguishing between dipole and directional antennas data. A clear picture of the analysed data groups is given in Tables 3 and 4 in Sec.IV.C.

D. STATISTICAL SIGNIFICANCE TEST
Finally, a statistical significance test was carried out to quantify the difference between the correlation coefficients and analyze values of r jk , j = {n, tot}, k = {1, 4} in a consistent and objective manner [30]. Such a procedure is also helpful to check whether the noise in the analyzed data had significant effect on the computed values of r jk . In particular, it could be determined if the probability that the difference between r ap and r bq , a = b, is merely due to chance (i.e., noise in the data) with a predefined significance level α (e.g., α = 5% if a confidence level of 95% is chosen). A simplified scheme of the statistical hypotheses for the testing procedure herby adopted is described as follows: 1. Define the null and alternative hypotheses. In this analysis, the null hypothesis is defined as H 0 : ρ ap = ρ bq , a = b or p = q, where ρ ap and ρ bq are the true correlation coefficients of the pairs of data samples (sPD a avg.p, T ) and (sPD b avg.q, T ), while r ap and r bq are their estimators. The alternative hypothesis is defined as H 1 : ρ ap > ρ bq , a, b, p, q such that ρ ap > ρ bq . 2. Employ an appropriate significance test, together with a relevant test statistic. Since the compared correlations have one variable in common (i.e., T ), the adopted test statistics for dependent correlations are those provided by the cocor package [30]. 3. Choose the significance level α. In this analysis, a value of α = 5% has been adopted. 4. Compute the values of z-scores or t-scores (i.e., the test statistics) associated to the difference r apr bq based on the observed data. 5. Reject the null hypothesis H 0 in favour of the alternative hypothesis H 1 if the observed value r apr bq is in the critical region; ''fail to reject'' (retain) H 0 otherwise.
The last two steps of the testing process are automatically performed by the tool provided in the cocor package [30]. Results in terms of z-scores or t-scores for each of the nine different applied tests are reported in Sec. IV.D.
The procedure above was applied to compare the performance of the normal and total component definitions in terms of correlation with T , distinguishing the compared coefficients by the averaging area.

A. SCATTER PLOTS
The scatter plots for the dipole and directional antennas are reported in Figs. 1-4 and 5-8, respectively.     d ≥ 2 mm, have been included in these plots. They provide insight into the collected data spread amongst the several institutions and simulation case studies, noting sporadic high values of HFs for some of the directional antennas, with particular reference to a pair of values with p T between 2 and 2.5 • C (see Figs. [5][6][7][8]. These values correspond to patch array antennas with a beam shift in azimuth or elevation, which have already been reported to exhibit a large spread of HFs [24]- [26]. A quantification of possible effects of these data on the Pearson correlation coefficients is given in Sec. IV.C (see Tables 3-6).

B. RELATIVE STANDARD DEVIATIONS OF HF
The relative standard deviations of the HF obtained for the several PD definitions and metrics from the several research groups are summarized in Table 1. It can be observed from this Table that the psPD tot definition always gives the HF distributions with the least relative standard deviations, except for the case of data provided by NITech group, where for the psPD n definition it is slightly smaller.
The histograms of the HF distributions are shown in Figs. 9-12. The average value of each distribution is demarked with a red dotted vertical line. The occurrence distributions show a noticeable positive skewness in all four cases. Thus, the ratio of the 95 th to the 5 th percentiles of the HF distributions has been performed to quantify the skewness and the tail length [27]. The computed ratios for each of the four HF definitions are reported in Table 2. Considering that the psPD tot values are in general higher than the psPD n while the p T are the same, it can be expected that the associated HF distributions show less dispersed values towards the high end and, thus, a shorter tail. In fact, Table 2 shows that percentile ratios associated with HF tot definitions (bolded highlighted) are smaller than those associated with HF n .     Finally, the plots of the HF rearranged as a function of electric distance, i.e. the distance d over the wavelength λ, distinguished by the PD definitions (HF n vs. HF tot ), are reported in Figs. 13 and 14 for the 1 cm 2 and 4 cm 2 averaging areas, respectively. The graphical evaluation of the plots shows that the average value of HF n is consistently greater than the one of HF tot , but this difference diminishes with electrical distance because the tangential components of the Poynting vector tend to vanish in the far-field. Instead, a greater dispersion is observed at distances d ≤ 0.6 λ (see Figs. 13(a) and 14(a)) due to the near-field zone characteristics. Tables 3 and 4 show the computed values of Pearson correlation coefficients for the dipole source and directional antennas, respectively. In the first case, the PD metric with highest r jk is the one related to the definition of sPD tot1 , while in the second case, coefficients related to sPD n1 are always higher, with the exception of the 10, 30 GHz data group, where the highest coefficients are r totl and r tot4 , respectively. This may be due to the increased capability to focus the beam for the directional antennas.

C. PEARSON CORRELATION COEFFICIENTS
In addition, when grouping data by separation distance, all the highest r jk values are obtained for d ≥ 5 mm, except for r tot4 in the case of directional antennas at the highest frequencies, where including all distances in the analyses seems to yield higher correlation. In general, the increased correlation, obtained excluding the smallest distance, may be due to the reduced level of tangential components of the Poynting vector in the near-field compared to the reactive components. Finally, when grouping data by the frequency range, it can be observed, as expected, that both PD definitions with 1 cm 2 averaging area correlate better than the respective definitions with 4 cm 2 for the 60, 90 GHz data set, while this is not always true for the 10, 30 GHz data set. Tables 6 and 7 show the results of the z-and t-scores analyses considering the data related to 1 cm 2 and 4 cm 2 averaging areas, respectively. The input values provided to the cocor package in order to perform the significance test are reported in Table 5. In particular, it shows the correlation coefficients compared in pairs, i.e., r n1 , r tot1 and r n4 , r tot4 , and the associated coefficients r n1,tot1 and r n4,tot4 , obtained computing the correlation of sPD n avg.1 with sPD tot avg.1 and sPD n avg.4 with sPD tot avg.4. A one-tailed test for the hypothesis that r n < r tot has been performed.

D. STATISTICAL SIGNIFICANCE TEST
The z-and t-scores reported in Tables 6 and 7 show that the results have strong statistical significance, with scores well below the critical values related to 95% confidence level. This result confirms that the correlation of the total component   with the peak temperature rise is (statistically) significantly higher than that of the normal component. In fact, all the employed tests led to the rejection of the null-hypothesis r n = r tot . The only exception to this trend comes from the comparison between r n1 and r tot1 for the directional antennas data, where the two correlations could not be distinguished at the 95% confidence level (see Table 5).

A. EXPLANATION OF THE OBTAINED RESULTS
The results of Table 3 show that, for the simple dipole source, the definition with sPD tot presents highest values of correlation coefficients with temperature rise compared to the TABLE 6. z-scores and t-scores of the psPD avg.1 and p T correlation comparison between data belonging to the dipole source and to the directional antennas. In the Null-Hypothesis column, rej. stands for rejected and ret. stands for retained. TABLE 7. z-scores and t-scores of the psPD avg.4 and p T correlation comparison between data belonging to the dipole source and to the directional antennas. In the Null-Hypothesis column, rej. stands for rejected and ret. stands for retained. definition sPD n , both when grouping the data by separation distance and frequency. The small difference between the values of correlation coefficients corresponding to the two definitions indicates that both resulting quantities correlate strongly with temperature rise (correlation coefficients >0.7).
One possible explanation for this is that the scenarios considered in this study are mainly for the normal incidence of the field on the body, thus the differences are relatively small. However, this slightly better correlation with the total components may be attributed to the near-field exposure conditions, where some of the components are obliquely incident to the model. Furthermore, excluding the data for the shortest distance of 2 mm, which is considered as near-field exposure (at least below 30 GHz), the correlation is generally improved. In other words, the correlation at 2 mm is very poor highlighting the limitation of IPD as RL/ERL metric at close distances.
Differently from the dipole, exposure to more complex sources yields to a less clear picture. In fact, while the power density and heating values corresponding to the dipole sources are consistent and do not show any particular clustering of values or alarming variability, the data coming from the simulations of directional antennas are obviously characterized by larger variations in the correlation coefficients (see Table 4) and by the presence of some outliers in their distribution. The latter can be found also in the scatter plots (see Figs. 5-8) with some p T values between 2 and 2.5 • C. Considerations about possible health risks from such high temperature increases and conservative safety margins lie outside the scope of this work.

B. METHOD AND DATA COLLECTION LIMITATIONS
When conducting a statistical analysis, a crucial aspect is the analysis of the available set of data, both in terms of sample size and nature of the collected data. In particular, for the former aspect, it is a good practice to establish whether the size of the analyzed samples is sufficient to guarantee a certain significance of the statistics. In this study, it has been verified that all the chosen samples for both dipole data (115) and directional antennas data (112) were large enough to properly represent their respective statistical population of values. Specifically, an analysis on the tolerance limits conducted through t-statistics has shown that, at the 95% confidence level, the errors obtained by estimating the population means through the samples average falls between ∼11.8% and ∼16.4%.
Regarding the nature of the collected data, a critical aspect may be represented by the fact that they come from different organizations. However, a study of body model and thermal parameter variability impact for frequencies below 30 GHz was performed in [27] showing that the difference in HFs were globally below 20%. The authors of [27] concluded that, although the slight dependency on the body model, thermal parameters, and antenna models, the deviation of HFs is insignificant when considering the numerical methods used by different organizations.

VI. CONCLUSION
The observation of correlation coefficients, computed for the different data samples, shows that, in general, the spatial averaged power density metrics based on total components (i.e., magnitude) provide a slightly higher correlation than those based on the normal components. However, the difference is marginal (0.9 vs. 0.8) and the close values of correlation coefficients show that both definitions correlate with temperature elevation.
Moreover, the analysis of the heating factor distributions for different normal exposure configurations shows smaller relative standard deviation values when applying an sPD tot definition. Consequently, this definition is expected to yield a slightly better estimate of the induced temperature increase than sPD n . The observed difference is however not large and can mainly be attributed to near-field conditions.
Finally, the analysis of the statistical significance test confirmed that the PD definition using all total components correlates with temperature rise slightly better than the PD definition using the normal components (see Table 5). On the basis of the obtained results, we conclude that the sPD tot is the definition that should be recommended as IPD metric in International exposure guidelines/standards for human protection.
ANTONIO DI FRANCESCO received the Laurea degree (Hons.) in mechanical engineering and the Laurea Magistrale degree (Hons.) in mathematical engineering from the University of L'Aquila, L'Aquila, Italy, in 2017 and 2020, respectively, where he is currently pursuing the Ph.D. degree in mathematics and models.
His current research interests include the physical and mathematical modeling and computer simulation of bioheat conduction and thermal effects of EM fields on biological tissues, nonequilibrium statistical physics, and in stochastic modeling of physical and biological phenomena.  Dr. El Hajj is participating and leading several standardization efforts in the human exposure and product safety domain. He is mandated as an Expert in French Standardization Association (AFNOR). Since 2017, he has been participating in the development of several IEEE/IEC standards on human exposure computational and measurement assessments. He is a member of IEC TC 106 and IEEE ICES TC95. He is the Convener of Working Group five under Subcommittee six of IEEE ICES TC95 established to study the different aspects of incident power density definition in correlation with temperature elevation. He is also a member of CMC TF Radio Group in IECEE. VOLUME 10, 2022