INTRODUCTION
Human exposure limits for the millimeter (mm)-wave radiofrequency (RF) band (i.e., above 6 or 10 GHz) transition from limiting whole-body averaged and spatial-peak specific absorption rate (SAR) to specifying limits based upon maximum power density (^{ICNIRP 1998} ; ^{IEEE 2005} ). This is due to the highly superficial nature of mm-wave energy deposition, which is known to produce surface heating of the skin and is not conducive to volumetric averaging (^{Wu et al. 2015} ). The limits are specified both in terms of their temporal steady-state values, under the assumption that the exposure is indefinite, and in terms of a time-averaging period for which higher intensity exposures of shorter durations are allowed.

The steady-state limits pertain to uniform exposure over large areas of the body (whole-body exposure), with separate provisions for localized exposure in terms of spatially averaged power density. Because of the nature of propagation at these frequencies, wireless devices operating in the mm-wave band will likely make use of adaptive, phased-array techniques, producing narrow beams (^{Andrews et al. 2014} ). Coupled with the short wavelengths from mm-waves, exposures from such devices are likely to be highly localized. In addition to the localized nature of the exposure, the effect of tissue curvature takes on less importance as the wavelength decreases. This factor serves to increase the relevance of planar models when characterizing skin heating due to exposure.

The whole-body power density limits are based upon an acceptable surface temperature rise in the skin, which has been extensively modeled (^{Alekseev and Ziskin 2003} ; ^{Kanezaki et al. 2010} ; ^{Foster et al. 2016} ; ^{Hashimoto et al. 2017} ) and to a lesser extent, measured (^{Blick et al. 1997} ; ^{Walters et al. 2000} ; ^{Nelson et al. 2003} ; ^{Alekseev et al. 2005} ; ^{Adibzadeh et al. 2016} ). It is known that as the irradiated area decreases, maximum steady-state temperature rise also decreases for the same spatial-peak power density (^{Foster et al. 2016} ). With characterization of narrow-beam temperature rise response, spatial-averaging schemes can be developed that result in similar surface temperature rise as for plane waves. In addition, criteria for spatial-peak power density as a function of beam width can be developed that maintain maximum temperature rise below a predefined threshold. Thus, a simple and direct method for calculating steady-state temperature rises in response to narrow-beam exposure would be beneficial in the development of such schemes.

Steady-state temperature rise in tissue has been calculated through solution of the Pennes bioheat transfer equation (BHTE) (^{Pennes 1948} ). Response to whole-body or large-area exposure is obtained from solutions of the one-dimensional, planar formulation (^{Anderson et al. 2010} ; ^{Kanezaki et al. 2010} ; ^{Sasaki et al. 2017} ), while responses to narrow-beam exposure have been obtained from solutions of a cylindrically symmetric, two-dimensional formulation of the BHTE (^{Alekseev et al. 2005} ; ^{Zhadobov et al. 2015} ). Many of these works have concentrated on a specific multilayer model (i.e., with fixed layer thicknesses) that may not reflect the true variation of thicknesses among different individuals and body sites. ^{Anderson et al. (2010)} used a Monte Carlo approach that incorporated statistical tissue thickness variations in a one-dimensional formulation up to 10 GHz, while ^{Sasaki et al. (2017)} performed a similar analysis from 10 GHz to 1 THz.

The wide variation of tissue thicknesses suggest the need for information on temperature rise due to narrow-beam exposure, taking into account the statistical variation of tissue layer thicknesses and in a form that is easily calculated and which can be used in the development of mm-wave exposure standards. This work has aimed at developing a relatively simple analytical model for the steady-state temperature rise that embodies the statistical variation of layer thickness at frequencies from 10 to 80 GHz. Monte Carlo techniques are used in conjunction with finite-difference (FD) solutions of the BHTE where localized exposure is treated by solving the two-dimensional, rotationally symmetric version (henceforth referred to as a 3D solution). Finally, the results of the Monte Carlo/FD solutions are used to derive parameters that when substituted into a closed-form formula can reproduce the Monte Carlo/FD temperature rise data.

METHODS
It will be assumed that all SAR distributions and resulting temperature rises in the tissues result from far-field, normally incident fields on planar tissues. Table 1 explains symbols, units, and quantities. For a circularly symmetric beam, the distribution of intensity that would occur on the surface (without the presence of the tissues) is defined as the projected power density S (r ), where r is the radial distance from the center of the beam. For the purpose of this study, the projected power density beam width (m) is taken to be the span over which the intensity falls to one-half the maximum value occurring at the center and is defined as the half-power projected beam width (HPBW; shown in Fig. 1 ). For frequencies from 28 GHz up to 75 GHz, the smallest HPBWs for which far-field conditions may exist lie within the bounds 0.01–0.004 m, respectively (based upon the idealized case of an isotropic point source located one-half wavelength from the projection plane).

Table 1: Symbols, units, and quantities.

Fig. 1: Definitions of half-power projected beam width (HPBW) of the projected power density distribution and full width at half maximum (FWHM) of the transverse SAR distribution from a narrow, circularly symmetric power density beam. Both distributions are presented as Gaussian with width parameters g and g _{s} related to FWHM and HPBW, respectively.

The resulting distribution of SAR in the tissues is also circularly symmetric and through any transverse slice at z = z _{o} , has a bell-shaped intensity distribution, SAR (r, z _{o} ), as shown in Fig. 1 . The transverse distribution of SAR will be assumed to be Gaussian, whose width at half-intensity will be termed the full width at half maximum (FWHM). The FWHM is related to the HPBW but, in general, is not equal to it. The relationship between the two parameters is dependent upon the source that produces the power density.

The tissue surface is assumed to extend to infinity in the directions transverse to the propagation direction of the electromagnetic wave and is usually truncated in depth as shown in Fig. 2 . The truncated or maximum depth is chosen such that the calculated temperature rise approximately decays to zero, which would hold true if the tissue extended to infinity.

Fig. 2: Geometry of the electromagnetic/thermal problem for a 3-tissue configuration. Transverse electric and magnetic (TEM) waves are normally incident on the tissue half-space with incident (unperturbed), plane-wave power density S _{inc} .

Calculation of axial SAR distribution
Calculation of the axial SAR distribution follows the method employed in ^{Drossos et al. (2000)} and ^{Kanezaki et al. (2010)} . This method makes use of the known solution of a plane wave propagating in layered, lossy media where the amplitudes of the forward- and reverse-travelling electric and magnetic fields constitute the unknowns in an algebraic system of equations. The system is solved for the unknown field amplitudes, from which the SAR can be calculated as a function of the axial coordinate (i.e. from the surface into the tissue). The intensity or amplitude of the resulting one-dimensional (1D) SAR distribution, SAR (z ), is directly proportional to the incident power density, S _{inc} , which is used as the normalizing quantity for the 1D temperature rise.

Calculation of steady-state temperature rise
Both 3D (from a Gaussian projected beam) and 1D (from a plane wave) solutions of the BHTE are required to develop the analytical model. The BHTE in cylindrical coordinates (i.e., having variation in the axial z and radial r directions) is:

where SAR (r,z ) is the 3D distribution of SAR in the tissues and U = U (r,z ) = T (r,z ) – T _{blood} represents the difference between the local tissue temperature T (r,z ) and the baseline temperature (taken to be T _{blood} = 37^{o} C).

For layered tissues, the z -dependent thermal/physical parameters are: m _{b} , the blood perfusion (volumetric) rate in units of m^{3} kg^{−1} s^{−1} ; ρ , the mass density in units of kg m^{−3} ; and k , the tissue heat conductivity in units of W m^{−1 o} C^{−1} . Parameters that are uniform throughout all layers are: C _{b} , the specific heat capacity of blood in units of J kg^{−1 o} C^{−1} ; and ρ _{b} , the mass density of blood.

Equation (1) is solved subject to the boundary condition on the surface

where T _{air} is the ambient temperature of the surrounding air (for all subsequent calculations, an air temperature of 22^{o} C is assumed), T (r, 0) is the temperature on the surface of the tissues, and h is a heat transfer coefficient with units of W m^{−2 o} C^{−1} . The term on the left-hand side of eqn (2) represents the rate of heat flow away from the surface. If insulating or adiabatic conditions exist on the surface, then h = 0 in eqn (2).

Three additional boundary conditions are used to solve eqn (1) (^{Alekseev et al. 2005} ). One is ∂U/∂r = 0 along the z -axis, which arises because of the circular symmetry. Another is T (r,∞ ) → T _{blood} (or U (r,∞ ) = 0), which arises because the RF power is absorbed superficially, and the temperature rise, therefore, decays to zero as tissue depth increases. The third is T (∞ ,z ) = T _{null} (z ) or U (∞ ,z ) = U _{null} (z ) where T _{null} (z ) is the temperature distribution in the tissues in the absence of RF exposure. It is found by solving the 1D BHTE (as a function of z only) with SAR (z ) = 0, which is discussed in more detail below.

For a 1D solution of eqn (1) (i.e., for plane-wave exposure), the differential terms with respect to r are dropped and only the axial SAR distribution (i.e., SAR (z )) is used. For this case, the solution U (z ) is a function only of the axial coordinate z . The net temperature rise for the 1D case ∆T _{1D} (z ) is the difference between the temperature calculated at full SAR T _{full} (z ) and that calculated with SAR (z ) = 0, T _{null} (z ), or is given by ∆T _{1D} (z ) = U _{full} (z ) – U _{null} (z ).

For the 3D case, the 1D null distribution represents the temperatures in the tissues prior to exposure (this distribution is independent of r ) and serves as the boundary condition as r → ∞. The net temperature rise for the 3D case ∆T _{3D} (r,z ) is equal to the difference between the 3D solution at full SAR U _{full} (r,z ) and the 1D null solution U _{null} (z ), or is given by ∆T _{3D} (r,z ) = U _{full} (r,z ) – U _{null} (z ).

For the 3D analyses, the 3D SAR distribution (i.e., SAR (r,z )) is obtained from the axial SAR distribution by multiplying the latter by a Gaussian term in r (^{Alekseev et al. 2005} ):

where g is the Gaussian width parameter, which is equal to g = 0.601 FWHM.

The SAR distribution of eqn (3) is assumed to have been the result of a narrow, projected beam of power density with maximum intensity in the center of the beam given by S _{o} .

Analytical model for steady-state temperature rise from a narrow, projected beam
Examples of distributions of steady-state surface temperature rise for a 4-layer tissue (under adiabatic conditions), illuminated with different projected beam widths at the same spatial-peak power density S _{o} , are shown in Fig. 3 . These distributions were obtained from finite-difference solutions of eqn (1), which will be dealt with later. In the figure, the FWHM of the transverse distribution of SAR is used as a surrogate for the projected HPBW, a feature that is used throughout this work.

Fig. 3: Examples of steady-state, surface temperature rise distributions for different FWHM at the same spatial-peak power density S _{o} (424 W m^{−2} ) for an adiabatic, 4-tissue configuration (skin: 1 mm; SAT: 2 mm; skull: 5 mm; and brain: 42 mm) at 10 GHz. The maximum surface temperature rise, ∆T _{o} = ∆T _{3D} (0,0), occurs in the center of the distribution (i.e., at r = 0).

An analytical solution for the temperature rise at r = 0 (i.e., the maximum temperature rise) was obtained by ^{Foster et al. (2016)} using a Green’s function formulation of the BHTE. Its derivation was based on the assumption of purely superficial SAR deposition on a homogeneous tissue and under adiabatic surface boundary conditions. The analytical model proposed here is based on the ^{Foster et al. (2016)} solution extended to the case of multilayer tissues, SAR depositions that penetrate into the tissue, and for convective surface boundary conditions. The maximum temperature rise (in the center of the beam) ΔT _{o} in the multilayer tissue is modeled as:

where the quantity (∆T _{plane} /S _{inc} ) is the maximal steady-state temperature rise per unit incident plane-wave power density (in the same tissue configuration) and R _{1-} _{eff} is the effective diffusion length for the multilayer tissues. The function erfc {X } is the complementary error function of the variable X .

For homogeneous tissues, the effective diffusion length becomes the diffusion length R _{1} given by (^{Foster et al. 2016} ):

where k , ρ , and m _{b} are constant throughout the single tissue. The effect that the parameter R _{1-} _{eff} has on the maximal temperature-rise response vs. beam width is shown in Fig. 4 . As seen from this figure, the value of R _{1-} _{eff} controls the rate at which ΔT _{o} /S _{o} asymptotically approaches ΔT _{plane} /S _{inc} (or for the case where S _{o} = S _{inc} , where the narrow-beam temperature rise approaches that of a plane wave). This parameter is determined by using a procedure that takes values of ΔT _{o} /S _{o} from 3D solutions for a range of values of FWHM and a single 1D solution to determine (∆T _{plane} /S _{inc} ) (for the same tissue configuration) and curve fits them to eqn (4).

Fig. 4: Maximum, normalized, steady-state, temperature rise ∆T _{o} /S _{o} normalized to (∆T _{plane} /S _{inc} ) as a function of the FWHM of a narrow, projected beam with Gaussian SAR distribution (from eqn 4) for three different values of effective diffusion length R _{1-} _{eff} .

The four parameters in the analytical model of eqn (4) can be divided into two categories. Parameters S _{o} and FWHM are deterministic, electromagnetic quantities obtained from the known exposure conditions (e.g., for a given source and SAR distribution). The remaining two parameters, (∆T _{plane} /S _{inc} ) and R _{1-} _{eff} , are thermally based quantities that depend on frequency and tissue configuration (i.e. composition and thicknesses of the different layers). Because of the random variation of layer thicknesses in the population and among different exposure sites, these two parameters are stochastic in nature and can be described by a statistic (e.g., a percentile value p ). When used with the analytical model (eqn 4), the deterministic exposure parameters along with the stochastic, tissue-related ones allow the narrow-beam temperature rise to be calculated without the need for solving the 3D BHTE. The goal of this study was to derive values of the two stochastic parameters as a function of frequency and tissue configuration.

Model parameter determination
Two different layered tissue configurations were used. One was comprised of three tissues consisting of thin layers of skin and subcutaneous adipose tissue (SAT) backed by a bulk layer of muscle. The other was a 4-tissue configuration comprised of thin layers of skin, SAT, and skull (cortical bone) backed by a bulk layer of brain (grey matter). The latter configuration was intended to represent head tissues (excluding the eye) while the former represents most other body surface locations. An additional configuration consisting only of bulk skin or muscle was used to the check the numerical algorithms against analytical solutions applicable to a single-tissue structure.

Solutions to the 1D and 3D BHTE were carried out using an FD formulation (^{Zhadobov et al. 2015} ). A scheme for gridding the tissues was employed that successively increased the grid spacing with increasing depth and radial distance. To ensure less than 1% variation in the 1D result as the grid spacing is varied, the minimum grid step in the z direction was set to 0.05 mm. Grid spacing could be relaxed for the 3D solutions to achieve the same goal because the determination of R _{1-} _{eff} is a process involving ratios of the temperature rise solutions of the 3D and 1D problems. Provided the same grid in the z direction was used for both, effects of the z -grid resolution largely cancelled out. In this case, the smallest grid step was set to 0.1 mm in the z direction and 0.2 mm in the r direction. Sparse matrix techniques were employed to solve these systems.

Curve fitting to determine R _{1-eff} involved the minimization of the sum of squared differences between FD-calculated and modeled temperature rise data using a nonlinear least-squares function in the programming software. For each tissue configuration, a single determination of (∆T _{plane} /S _{inc} ) (from the 1D BHTE) and four determinations of ΔT _{o} /S _{o} (from the 3D BHTE), one for each value of FWHM, were carried out. Values of FWHM equal to 0.005 m, 0.017 m, 0.035 m, and 0.060 m were used for this purpose.

Monte Carlo simulations were performed with the 1D and 3D solvers to extract statistics for (∆T _{plane} /S _{inc} ) and R _{1-} _{eff} , respectively. Each iteration of the 1D Monte Carlo simulation required two solutions of the 1D BHTE to obtain the null and full temperature distributions and subsequently, the maximum net temperature rise ∆T _{plane} . An iteration of the 3D Monte Carlo simulation included the previous 1D procedure (to find both U _{null} (z ) and ∆T _{plane} ) as well as four solutions of the 3D BHTE at the four values of FWHM given above. This was then followed by a curve fitting to obtain R _{1-} _{eff} .

In each iteration of either the 1D or 3D Monte Carlo simulation, the thickness of skin and SAT (3-tissue) or skin, SAT, and skull (4-tissue) were chosen randomly according to distributions derived from thickness data in the published literature (^{Hwang et al. 1999} ; ^{Mahinda and Murty 2009} ; ^{Anderson et al. 2010} ; ^{Kakasheva-Mazenkovska et al. 2011} ). Different percentile values of (∆T _{plane} /S _{inc} ) were extracted from the ensuing 1D solution distributions while the mean and standard deviation (SD) of R _{1-} _{eff} were extracted from the 3D ones.

All calculations used the electrical and physical/thermal parameters for the various tissues obtained from the IT’IS database (^{Hasgall et al. 2018} ).

Validation
For a single tissue and for the 1D or plane-wave case, eqn (1) is an ordinary differential equation that can be solved using well-known methods. The complementary solutions are given by the functions A _{1,2} exp {±z/R _{1} } where A _{1} and A _{2} are constants. The form of the particular solution is obtained from the 1D SAR distribution given by (^{Foster et al. 2016} ):

where δ is the power penetration depth (the depth at which SAR (z ) drops to exp {−1 }) and T _{R} is the power transmission coefficient (the fraction of the incident power density that is transmitted across the interface).

From eqn (6), the particular solution is the function B exp {−z/δ } where B is a constant. After applying standard methods, all unknown constants can be found, resulting in the solution U (z ) (the details are omitted for brevity). Recalling that the absolute temperature in the tissues is T (z ) = U (z ) + T _{blood} , it can be written as:

The absolute temperature in eqn (7) has two terms, one is the net temperature rise (the term proportional to S _{inc} ) and the other is the null-exposure temperature distribution (independent of power density). Eqn (7) illustrates the significance of the diffusion length R _{1} , which also governs the rate of decay of temperature rise (for δ < R _{1} ) in the deep tissues.

For the adiabatic case (h = 0), the maximum temperature rise occurs at the surface while for the convective heat-loss case, it occurs a very short distance inside the surface of the tissue. It can be shown that the single-tissue temperature rise maximum ΔT _{single-max} at this location is given by:

while the temperature rise on the surface ΔT _{single} (0) is given by:

The difference between ΔT _{single} (0) and ΔT _{single-max} is small; for example, being on the order of 0.2%, 0.6%, and 1.1% for values of h = 10, 20, and 30 W m^{−2 o} C^{−1} , respectively, and for skin tissue at 20 GHz. From eqn (9) it is seen that the convective heat-loss condition produces smaller temperature rises than the adiabatic condition.

Validation calculations were carried out for skin (R _{1} = 0.00670 m) and muscle (R _{1} = 0.0131 m) using the 1D FD solver with a maximum tissue depth of 50 mm. For muscle, this depth represents about 4 times the diffusion length R _{1} . Since muscle has the larger value of R _{1} between the two bulk tissues (muscle vs. brain), the temperature rise decays faster in brain. The differences between the analytical values of ΔT _{plane-max} given by eqn (8) and the equivalent 1D FD results ranged from less than 0.05% at 10 GHz to 1% at 80 GHz for both skin and muscle configurations and for either adiabatic or convective heat-loss boundary conditions.

Additional validation was obtained by observing the 3D FD-derived value of R _{1-} _{eff} for a single tissue under adiabatic conditions. As noted earlier, eqn (4) is exact under the condition where the penetration depth approaches zero. Thus, the curve-fitted value of R _{1-} _{eff} for a single tissue, obtained from a set of 3D FD solutions at different FWHM, should converge to the value of R _{1} given in eqn (5) as frequency increases. Fig. 5 shows results of the curve-fitted value of R _{1-} _{eff} vs. frequency for a skin-only configuration from 10 to 80 GHz. The value of R _{1-} _{eff} obtained at 80 GHz is 0.00678 m, which is a 1.2% difference from the exact value R _{1} = 0.00670 m.

Fig. 5: Curve-fitted values of R _{1-} _{eff} vs. frequency, obtained from narrow-beam, 3D FD solutions for a skin-only tissue.

A further validation of the 3D FD solver was carried out using the software package Sim4Life (ZMT Zurich MedTech AG, Zurich, Switzerland), which possesses both voxel-based finite-difference time-domain (FDTD) electromagnetic and thermodynamics solvers. Fig. 6 shows the results of 3D FD and Sim4Life calculations for a 3-tissue configuration (skin: 1.5 mm; SAT: 1.5 mm; and muscle: 67 mm) at 30 GHz under adiabatic conditions. All electrical and thermal/physical parameters in the Sim4Life simulation were set to be identical to the ones found in ^{Hasgall et al. (2018)} and used in the 1D and 3D FD solutions.

Fig. 6: Comparison of Sim4Life results with 3D FD for the 3-tissue configuration (skin: 1.5 mm; SAT: 1.5 mm; muscle: 67 mm) at 30 GHz under adiabatic conditions. Continuous lines represent curve-fitted representations of the above numerical data using the model of eqn (4).

To implement the electromagnetic source in the Sim4Life simulations, two orthogonal edge sources were placed a distance from the 3-tissue block. The distance was made variable to obtain different projected power density beam widths and values of FWHM (in Sim4Life, an edge source is the electric field defined on one of the edges of a voxel). The two edge sources were excited in time-phase quadrature to produce circular polarization, which resulted in a circular, projected power density distribution on the surface of the tissue block. For each source distance, two simulations were performed, one without the tissue block to obtain the spatial distribution of the unperturbed projected power density and its peak value S _{o} . The second, with the block present, was to obtain the SAR distribution width (FWHM) and the resulting temperature rise distribution.

In Fig. 6 , the Sim4Life data (circles) appear to correspond well with the 3D FD results (triangles). Because of limitations on computing capacity, a simulation large enough to closely approach the 1D results could not be carried out, and thus the parameter (∆T _{plane} /S _{inc} ), had to be fitted from the Sim4Life data. Nevertheless, its value is in good agreement (approximately 2% difference) with the calculated 1D FD values shown in the figure. Comparison of R _{1-} _{eff} between the two methods gives values within 1% of each other.

From the Sim4Life results, the transverse distribution of the SAR can be examined in detail. This allows testing of the assumption of a Gaussian distribution of the SAR for the two edge sources. The transverse SAR distributions for two selected values of source-to-tissue spacing (and for which the FWHM values were measured from the Sim4Life SAR data) are shown in Fig. 7 . These are compared to analytical Gaussian distributions that give the best fit. From the figure, it is seen that the transverse SAR distributions from the edge sources are very close to Gaussian in shape, especially in the center of the distribution where most of the power is deposited.

Fig. 7: Transverse SAR distributions through a z slice just under the surface of the tissue block for FWHM = 0.068 m (source spacing: 50 cm in Sim4Life) and FWHM = 0.013 m (source spacing: 10 cm in Sim4Life).

Size of the hot spot in relation to the projected beam width
In addition to modeling the maximum temperature rise for narrow-beam exposures, the spatial distribution of temperature rise on the surface was characterized. It is seen in Fig. 3 that the spatial distribution is somewhat Gaussian shaped, which can be characterized by a half-maximum width parameter (in this case, the 50% isotherm diameter) D _{50} . The closeness of the fit of the computed surface temperature rise distribution to a true Gaussian was found to be excellent for relatively wide projected beam widths (and corresponding FWHM). For narrower beams (e.g., FWHM under 10 mm), the fit was found to be good for a radius up D _{50} /2. The value of D _{50} can be extracted from each solution of the 3D FD solver. Once ∆T _{o} and D _{50} are known for a given FWHM, the surface temperature rise distribution is fully characterized, and quantities such as spatially averaged temperature rise or other isothermic diameters can be computed.

From 3D solutions of selected tissue configurations and at different frequencies, it was observed that D _{50} was found to be approximately linearly related to the FWHM with a slope that is close to unity (for large FWHM, D _{50} asymptotically approaches the FWHM). In view of this, a model of the form D _{50} = FWHM + D _{o} was adopted so that only one parameter, D _{o} , is required to model the width of the surface temperature rise distribution. The intercept D _{o} is easily computed from the 3D FD output data using linear regression, and it will be seen that its value changes minimally with frequency. While simplistic in nature, the linear model proposed for D _{50} is sufficiently accurate for most purposes.

Statistical modeling of skin, SAT, and skull thicknesses
The statistical model of skin thickness was derived from data contained in ^{Anderson et al. (2010)} and ^{Kakasheva et al. (2011)} . ^{Anderson et al. (2010)} report average skin thicknesses for males and females in three age groups and for eight different body sites resulting in a total of 48 thickness values. ^{Kakasheva et al. (2011)} report mean thicknesses for five age groups and 15 body sites with no distinction between sexes. Some of the body sites in ^{Kakasheva et al. (2011)} were omitted from inclusion because of their location (e.g., soles, lower leg, etc.). There were a combined number of 93 data points used to estimate the distribution. A lognormal distribution with geometric mean (GM) of GM = 0.00166 m and geometric standard deviation (GSD) of GSD = 1.518 was used to represent the distribution of skin thickness in the Monte Carlo simulations.

The data in ^{Anderson et al. (2010)} were used for modeling the thickness of SAT. Mean SAT thicknesses were tabulated for both sexes in three age groups and for eight body sites. A lognormal distribution with GM = 0.00652 m and GSD = 1.781 was used to represent the distribution of SAT thickness in the Monte Carlo simulations.

Skull thickness data was obtained from ^{Mahinda and Murty (2009)} and ^{Hwang et al. (1999)} . The former study reported mean thicknesses at six sites for males and females from 76 individuals. The latter study reported mean thickness at 13 sites from 51 adult Korean individuals. Skull thickness was modeled with a uniform distribution with minimum and maximum limits of 2 mm and 9 mm, respectively.

RESULTS
1D Monte Carlo results—(∆T _{plane} /S _{inc} )
Distributions of (∆T _{plane} /S _{inc} ) for the 10–80 GHz range (in 10 GHz increments) were obtained as the result of 10,000 Monte Carlo iterations. At 30 GHz and above, the distributions were bell shaped with slight positive skewness and with median and mean values within 0.8%. At 20 GHz, there was small negative skewness with median and mean values within 0.3%. Distributions at 10 GHz were slightly peaked (small positive kurtosis) with small positive skewness and with median and mean values within 1.7%.

Percentile statistics at selected levels were obtained from the distributions at each frequency. Fig. 8 shows the 95th, 80th, and 50th percentile values for both the 3- and 4-tissue configurations under adiabatic conditions plotted vs. frequency. For all percentile levels and for both adiabatic and convective heat-loss conditions, differences between the 3- and 4-tissue statistics are quite small (typically within 2% of each other).

Fig. 8: Maximum plane-wave temperature rise per unit incident power density at different percentile levels vs. frequency from 1D Monte Carlo simulations on 3- and 4-tissue (tiss) (adiabatic) models from 10 GHz to 80 GHz.

The percentile values of (∆T _{plane} /S _{inc} ) display an almost linear behavior with frequency above 30 GHz. Below this frequency, they have characteristics of a high-pass function. Over the 10 to 80 GHz range, it was found that they can be modeled as a function of frequency by the product of a linear term and a high-pass function given by:

where the subscript p denotes the percentile level, f _{c} is the high-pass corner frequency, while A _{p} and B _{p} are the slope and intercept of the linear term at percentile level p , respectively.

Since the results for the 3- and 4-tissue configurations were so close in value, they were averaged together, and those averages were used in the fitting of the data to eqn (10). Table 2 gives curve-fitted values of the parameters in eqn (10) for a range of percentile levels and for the adiabatic and convective heat-loss (h = 10 W m^{−2 o} C^{−1} , T _{air} = 22^{o} C) conditions. The curve fitting was carried out by nonlinear least-squares over the frequency range 10–80 GHz. The difference between the modeled values using the parameters in Table 2 and the 3- and 4-tissue FD data is on average +1% and −1%, respectively, for the adiabatic condition and +0.7% and −0.7%, respectively, for the convective heat-loss condition.

Table 2: Tabulated values of the parameters A _{p} and B _{p} and f _{c} used in eqn (10) for calculating (∆T _{plane} /S _{inc} )_{p} as a function of frequency f (in GHz). Values are given for the cases of adiabatic (h = 0 W m^{−2 o} C^{−1} ) and convective (h = 10 W m^{−2 o} C^{−1} , T _{air} = 22^{o} C) boundary conditions.

At the same percentile level, the adiabatic results are approximately 30% larger than the corresponding ones for the case of convective heat loss (h = 10 W m^{−2 o} C^{−1} ).

3D Monte Carlo results—R _{1-} _{eff}
A total of 500 iterations were performed at each frequency for the 3-tissue configuration while for the 4-tissue one, 700 iterations were carried out. This was due to the fact that three parameters were varied for the latter (skin, SAT, and skull thickness). The distributions were bell shaped with slight negative skewness at all frequencies. Differences between the mean and median values were under 1% for both tissue configurations and boundary conditions and at all frequencies. Mean values of R _{1-} _{eff} ranged from 0.011 m (3-tissue, adiabatic, 10 GHz) to 0.0076 m (4-tissue, convective, h = 10 W m^{−2 o} C^{−1} , 80 GHz). Ratios of SD/mean (coefficient of variation) were found to be on the order of 7–8% and 2–3% for 4-tissue and 3-tissue configurations respectively, over the 10–80 GHz range and for both adiabatic and convective heat-loss conditions.

Fig. 9 presents the mean R _{1-} _{eff} vs. frequency for the 3- and 4-tissue (adiabatic and convective, h = 10 W m^{−2 o} C^{−1} ) Monte Carlo simulations, respectively. It was found that the mean could be modeled with respect to frequency by the high-pass function

Fig. 9: Mean R _{1-} _{eff} vs. frequency from 3- and 4-tissue 3D Monte Carlo (MC) simulations under adiabatic and convective heat-loss (h = 10 W m^{−2 o} C^{−1} , T _{air} = 22^{o} C) conditions.

where F _{m} is a different constant depending on whether the configuration is 3- or 4-tissue and f _{c} is the corner frequency.

Table 3 gives values of these parameters for the four separate cases consisting of 3-tissue/4-tissue configurations and with adiabatic/convective heat loss, as well as two cases where the 3- and 4-tissue data are averaged. The fitting of the mean R _{1-} _{eff} data to eqn (11) was carried out using nonlinear least-squares over the frequency range 10–80 GHz. The SD obtained from the Monte Carlo distributions of R _{1-} _{eff} are presented in Fig. 10 vs. frequency.

Table 3: Constants used for calculating the mean value of R _{1-} _{eff} using eqn (11), where the result, R _{1-} _{eff} , is in m and f is in GHz. Convective boundary condition has parameters h = 10 W m^{−2 o} C^{−1} and T _{air} = 22^{o} C.

Fig. 10: SD of R _{1-} _{eff} vs. frequency for the 3- and 4-tissue configurations under adiabatic and convective heat-loss (h = 10 W m^{−2 o} C^{−1} , T _{air} = 22^{o} C) conditions. Obtained from the output distributions of 3D Monte Carlo simulations.

3D Monte Carlo results—50% isotherm diameter intercept D _{o}
The distributions of D _{o} were found to be quite narrow, with a worst-case coefficient of variation (SD/mean) of under 8% for the 4-tissue, convective heat-loss condition. Generally, SD/mean of the 3-tissue configurations were half those of the 4-tissue configurations for both types of boundary condition. The results suggest that the value of the intercept is relatively insensitive to the skin and SAT thickness. Differences of the mean D _{o} between 3- and 4-tissue configurations, for the same boundary conditions, were under 7%. For simplicity, the mean values for the 3- and 4-tissue configurations were averaged and high-pass functions, similar to eqn (11), were fitted. Fig. 11 shows the average of the 3- and 4-tissue mean D _{o} vs. frequency for the adiabatic and convective heat-loss (h = 10 W m^{−2 o} C^{−1} , T _{air} = 22^{o} C) Monte Carlo simulations along with their fitted high-pass functions.

Fig. 11: Mean values of D _{o} vs. frequency from 3D Monte Carlo simulation. Curve-fitting was performed on the average of 3- and 4-tissue (tiss) results for the same boundary condition. Convective heat-loss boundary condition has parameters h = 10 W m^{−2 o} C^{−1} and T _{air} = 22^{o} C.

Applications
The development of a closed-form model for predicting narrow-beam temperature rise permits straightforward calculation of any one of the model parameters from knowledge of the others. For instance, the spatial-peak power density in a narrow beam that produces a fixed temperature rise can be computed as a function of the FWHM. This is easily carried out by solving eqn (4) for S _{o} and inputting the desired value of ∆T _{o} and percentile level of (∆T _{plane} /S _{inc} )_{p} and R _{1-} _{eff} . These latter two parameters are calculated using the data in Tables 2 and 3 and eqns (10) and (11), respectively. An example is shown in Fig. 12 a for the spatial-peak power density that produces 1^{o} C vs. FWHM for the 50th percentile (∆T _{plane} /S _{inc} )_{p} and mean R _{1-} _{eff} , both under adiabatic conditions.

Fig. 12: Spatial-peak power density of a narrow beam at four selected frequencies that produces 1^{o} C maximum surface temperature rise for (a) the 50th percentile, adiabatic (∆T _{plane} /S _{inc} ) and mean R _{1-} _{eff} (average of 3- and 4-tissue, adiabatic) and (b) the 95th percentile, convective (h = 10 W m^{−2 o} C^{−1} , T _{air} = 22^{o} C) (∆T _{plane} /S _{inc} ) and mean R _{1-} _{eff} (average of 3- and 4-tissue, convective).

A similar calculation was performed for the case of the 95th percentile (∆T _{plane} /S _{inc} ) and mean R _{1-} _{eff} (from the averaged 3- and 4-tissue data) for the convective heat-loss boundary condition: h = 10 W m^{−2 o} C^{−1} and T _{air} = 22^{o} C. The results are shown in Fig. 12 b. The corresponding curves (at the same frequencies) in Fig. 12 a and b are similar due to the fact that the 95th percentile, convective (∆T _{plane} /S _{inc} ) is close in value to the 50th percentile, adiabatic value. Both sets of curves show that for relatively wide beam widths, the 1^{o} C spatial-peak power density asymptotically approaches a constant value that is somewhat larger (depending on frequency) than the occupational or controlled, plane-wave limit of 50 W m^{−2} used in the International Commission on Non-Ionizing Radiation Protection (ICNIRP) exposure standard (^{ICNIRP 1998} ).

Another application of the temperature rise model could be in the calculation of spatially averaged surface temperatures using the 50% isotherm diameter D _{50} . Additionally, the simple nature of the model lends itself to investigations of spatial-averaging areas for power density in the development of exposure standards. For instance, ^{Hashimoto et al. (2017)} propose an averaging area of 400 mm^{2} , stating that nonuniform power densities averaged over this area and meeting the plane-wave limit produce temperature rise close to that of a plane wave at the limit (^{Hashimoto et al. 2017} ). This conclusion is tested using the model of eqn (4).

To begin, the projected beam in Fig. 1 , characterized by the width parameter g _{s} = 0.601 × HPBW and peak value S _{o} , is spatially averaged over a circular area of radius R _{a} by integrating the distribution over this radius (details are given in ^{Neufeld et al. 2018} ). If the resulting spatial average is made equal to the plane-wave limit, S _{Limit} , then the relationship between S _{Limit} and S _{o} is:

The temperature rise for this power density distribution, whose peak power density is defined by eqn (12), will be denoted as ∆T _{o,ave} and can be found by solving for S _{o} from eqn (12) and substituting the result into eqn (4). The resulting temperature rise is a function of both HPBW and FWHM, and thus the relationship between these two parameters must be characterized. For simplicity, a simple linear relationship can be assumed whereby FWHM/HPBW is a constant K (K < 1). In addition, to complete the calculation using eqn (4), a percentile value of (∆T _{plane} /S _{inc} )_{p} is chosen for the frequency of interest:

The averaging area that results in a fixed ∆T _{o,ave} (e.g., ∆T _{o,ave} = 1^{o} C) can be found by numerically solving eqn (13) for R _{a} given the assumed HPBW and values of (∆T _{plane} /S _{inc} )_{p} and R _{1-} _{eff} at the chosen frequency.

A slightly different approach can be taken whereby the efficacy of the averaging area is evaluated by observing how close the narrow-beam temperature rise ∆T _{o,ave} (for spatially averaged power density) is to the temperature rise from a plane wave at the plane-wave limit ∆T _{limit} . The latter can be calculated (for the percentile level p ) as: ∆T _{limit} = S _{Limit} (∆T _{plane} /S _{inc} )_{p} . A test ratio can be formed by the quotient ∆T _{o,ave} /∆T _{limit} , whereby a value of unity is ideal. This ratio is given by:

where it is seen that the test ratio is independent of (∆T _{plane} /S _{inc} )_{p} and is only a function of HPBW and R _{1-} _{eff.}

As an example, the test ratio for 28 GHz is plotted in Fig. 13 for several averaging areas (area = πR _{a} ^{2} ), that include the ^{ICNIRP (1998)} recommendation of 2,000 mm^{2} and the 400 mm^{2} area proposed in ^{Hashimoto et al. (2017)} . For this plot, the 3- and 4-tissue-averaged R _{1-} _{eff} (adiabatic) was selected, and it is assumed the power density and SAR distribution widths are related by FWHM/HPBW = K = 0.8. From the figure, it is clear that the 2,000 mm^{2} averaging area is too large since its curve never crosses unity. The 400 mm^{2} averaging area gives values of the test ratio both above and below unity over a wide range of HPBWs, showing that it is a reasonable choice. The figure demonstrates the trade-offs that are incurred by perturbing the size of the area about the proposed 400 mm^{2} area.

Fig. 13: Plot of the test ratio

∆T _{o,ave} /

∆T _{plane} vs. HPBW for different spatial-averaging areas where

∆T _{o,ave} is the temperature rise for a narrow beam with spatially averaged power density at the plane-wave limit and

∆T _{plane} is the temperature rise for a plane wave at the limit. An assumed value of

R _{1-} _{eff} = 0.00935 m was used corresponding to the mean value of the 3- and 4-tissue average under adiabatic conditions at 28 GHz (shown in

Fig. 9 ). In addition, a relationship between power density and SAR distributions, given by FWHM/HPBW =

K = 0.8, is assumed.

DISCUSSION
One limitation of this work is the assumption that exposure is from the far field. At frequencies in the mm-wave range, this is not too serious a limitation since wavelengths and antenna lengths are short, as are source separations that would qualify as far field. Expanding the model to near-field exposures would likely need to be done on a case-by-case basis and would probably require a much more detailed electromagnetic analysis of the source.

An additional limitation is that effects of local vasodilation were ignored (i.e., the blood perfusion rate at the site of exposure is assumed to be the same in the postexposure steady-state as in the preexposure state). Thus, the temperature rise model is pertinent to skin temperature elevations below 39^{o} C (^{Margerl and Treede 1996} ; ^{Walters et al. 2004} ). Since exposure standards are usually designed to restrict temperature rise at or below 1^{o} C, the model developed in the current study is still relevant to this purpose.

The uppermost frequency of 80 GHz was chosen since the highest frequency in the 5th generation wireless band is in the region of 75 GHz. Yet, this choice does not appear to limit the applicability of the analytical model for higher frequencies; as it can be seen from Figs. 8–11 , the model parameters approach asymptotic limits (a straight line for (∆T _{plane} /S _{inc} )_{p} and a constant for mean R _{1-} _{eff} ) at frequencies approaching 80 GHz and beyond.

Throughout this study, the SAR distribution width, FWHM, has been used as a surrogate for the projected beam width of the incident power density, HPBW. For practical purposes, it would be desirable to characterize temperature rise directly in terms of the HPBW of a device since that is more easily measured. For a specific source or antenna, the relationship between the HPBW and the resulting FWHM can be characterized through the use of electromagnetic simulation methods such as FDTD. This approach gives the model developed in this work broader applicability than analyzing specific antennas and spacings.

As an example, the FDTD edge sources used in the Sim4Life validation described above were evaluated for both the HPBW and resulting FWHM as a function of distance. It was found that there was an approximate linear relationship (with zero intercept) between the two with a ratio FWHM/HPBW = 0.82 at 30 GHz. Thus, for this source and frequency, the temperature rise model of eqn (4) could easily be modified by substituting 0.82 HPBW for the parameter FWHM, making it a function of the projected power density beam width.

Comparison of the 1D results in this work can be made to that of ^{Sasaki et al. (2017)} who performed a 1D Monte Carlo analysis of a 4-tissue configuration consisting of epidermis, dermis, SAT, and muscle from 10 GHz to 1 THz. Their solution of the 1D BHTE was carried out for a heat transfer coefficient of h = 10 W m^{−2 o} C^{−1} and an air temperature of 20^{o} C. Comparison of the mean values of (∆T _{plane} /S _{inc} ) in Fig. 5 of their manuscript and the convective (h = 10 W m^{−2 o} C^{−1} ), 50th percentile values computed from Table 2 of this work give differences ranging from 1.2% to 4.9% over the 10–80 GHz range (with the values obtained in the current study consistently larger). ^{Sasaki et al. (2017)} also note an increase in the mean (∆T _{plane} /S _{inc} ) of 25% when going from a heat transfer coefficient of 1 to 10 W m^{−2 o} C^{−1} . By comparison, an increase in the heat transfer coefficient from 0 (adiabatic) to 10 W m^{−2 o} C^{−1} in the current study (with T _{air} = 22^{o} C) produces an increase of 29% for the 50th percentile (∆T _{plane} /S _{inc} ) over the 10–80 GHz range.

It should be noted that ^{Sasaki et al. (2017)} also independently measured the dielectric parameters of dermis, SAT, and muscle and found large differences between the measured SAT parameters and those reported in the database of ^{Hasgall et al. (2018)} . However, they conclude from a Monte Carlo analysis of the effect of dielectric parameter values that values of (∆T _{plane} /S _{inc} ) are more strongly affected by tissue thickness variations than by variations in the dielectric parameters. They further noted that the results were more strongly affected by epidermis/dermis thickness than by SAT and muscle thickness. Finally, they observed that the dielectric parameters of SAT had little effect on the temperature rise above 30 GHz.

^{Neufeld et al. (2018)} analyzed the relationship between circular averaging areas and steady-state temperature rises for spatially averaged power densities of 10 W m^{−2} for several antenna types and beams. Their criterion for temperature rise was to be below 1^{o} C in all cases. They conclude that the maximal averaging area is frequency dependent with areas ranging from 300 mm^{2} at 10 GHz to 190 mm^{2} at 100 GHz.

The approach to arrive at eqn (13) in the current study is similar to that in ^{Neufeld et al. (2018)} except that the distribution widths of SAR and projected power density are treated as separate (but related) entities. Whereas ^{Neufeld et al. (2018)} modeled the relationship between temperature rise and power density beam width or source distance for a multitude of antenna types and spacings, they did not explicitly relate narrow-beam temperature rise to SAR distribution width as done in the current study.

CONCLUSION
An analytical model consisting of eqns (4), (10), and (11) has been developed that enables maximum, steady-state surface temperature rise due to a narrow beam of power density to be calculated using basic computational tools such as a spreadsheet. The input parameters related to the exposure are the spatial-peak power density S _{o} and a surrogate parameter for the projected width of the beam on the tissue surface HPBW. This surrogate parameter is the transverse width of the resulting SAR distribution FWHM. While not equal to the HPBW, the relationship between the two can be ascertained from numerical simulation of the actual source of the power density.

The model was developed for two tissue configurations representing head (4-tissue, excluding the eye) and general body sites (3-tissue) and for two surface boundary conditions, adiabatic and convective heat loss (heat-loss coefficient of 10 W m^{−2 o} C^{−1} and 22^{o} C air temperature). Monte Carlo techniques, using human tissue layer thickness data obtained from the published literature, were applied to solutions of the 1D and 3D BHTE. From the resulting Monte Carlo output distributions, statistics of the two additional input parameters in the analytical model were extracted. These (stochastic) parameters were the temperature rise per unit incident plane-wave power density at different percentile levels and the mean effective diffusion length. The analytic model was developed over the 10–80 GHz frequency range, and the two stochastic input parameters were themselves modeled using simple closed-form expressions.

For selected tissue configurations, the numerical solutions to the 1D and 3D BHTE were compared against known analytical solutions and results from a commercial software package (Sim4Life). The differences were typically no greater than 1–2%. Comparison of the 50th percentile 1D results with those in the literature showed correspondences within 5% despite the use of different dielectric data for SAT and skin. This corroborates the conclusions of other authors that temperature rise is more strongly dependent on tissue thickness, especially that of skin and SAT, than on the value of the dielectric parameter.

An example of the use of the model was given for the calculation of spatial-peak power density vs. beam width that produces a predefined maximum surface temperature rise. It was shown that the current ICNIRP occupational limit for plane waves or large-area exposure produces temperature rise below 1^{o} C. A second example evaluated the effect of spatial-averaging area on the maximum steady-state temperature rise in relation to that for a plane wave at the plane-wave power density limit. Quantitative comparison between different averaging areas as a function of beam width was demonstrated. It was subsequently shown that the current ^{ICNIRP (1998)} averaging area of 2,000 mm^{2} is far too large, while averaging areas on the order of 300–500 mm^{2} were more appropriate for consideration in future exposure standard development.

The result of this work is an analytical tool that can be used to compute temperature rise responses at any percentile level with respect to the variability of the population and at various body sites as a result of exposure to narrow-beam mm-wave RF. This information may be useful for the development of future exposure standards.