Secondary Logo

Journal Logo


Model of Steady-state Temperature Rise in Multilayer Tissues Due to Narrow-beam Millimeter-wave Radiofrequency Field Exposure

Gajda, Gregory B.; Lemay, Eric; Paradis, Jonathan1

Author Information
doi: 10.1097/HP.0000000000001036
  • Open



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.


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
Table 1:
Symbols, units, and quantities.
Fig. 1
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 = zo, has a bell-shaped intensity distribution, SAR(r, zo), 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
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, Sinc, 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) – Tblood represents the difference between the local tissue temperature T(r,z) and the baseline temperature (taken to be Tblood = 37oC).

For layered tissues, the z-dependent thermal/physical parameters are: mb, the blood perfusion (volumetric) rate in units of m3 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 oC−1. Parameters that are uniform throughout all layers are: Cb, the specific heat capacity of blood in units of J kg−1 oC−1; and ρb, the mass density of blood.

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

where Tair is the ambient temperature of the surrounding air (for all subsequent calculations, an air temperature of 22oC 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 oC−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,∞) → Tblood (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) = Tnull(z) or U(∞ ,z) = Unull(z) where Tnull(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 ∆T1D(z) is the difference between the temperature calculated at full SAR Tfull(z) and that calculated with SAR(z) = 0, Tnull(z), or is given by ∆T1D(z) = Ufull(z) – Unull(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 ∆T3D(r,z) is equal to the difference between the 3D solution at full SAR Ufull(r,z) and the 1D null solution Unull(z), or is given by ∆T3D(r,z) = Ufull(r,z) – Unull(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 So.

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 So, 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
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) ΔTo in the multilayer tissue is modeled as:

where the quantity (∆Tplane/Sinc) is the maximal steady-state temperature rise per unit incident plane-wave power density (in the same tissue configuration) and R1-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 R1 given by (Foster et al. 2016):

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

Fig. 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 So 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, (∆Tplane/Sinc) and R1-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 R1-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 R1-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 (∆Tplane/Sinc) (from the 1D BHTE) and four determinations of ΔTo/So (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 (∆Tplane/Sinc) and R1-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 ∆Tplane. An iteration of the 3D Monte Carlo simulation included the previous 1D procedure (to find both Unull(z) and ∆Tplane) 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 R1-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 (∆Tplane/Sinc) were extracted from the ensuing 1D solution distributions while the mean and standard deviation (SD) of R1-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).


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 A1,2expz/R1} where A1 and A2 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 TR 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) + Tblood, it can be written as:

The absolute temperature in eqn (7) has two terms, one is the net temperature rise (the term proportional to Sinc) and the other is the null-exposure temperature distribution (independent of power density). Eqn (7) illustrates the significance of the diffusion length R1, which also governs the rate of decay of temperature rise (for δ < R1) 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 ΔTsingle-max at this location is given by:

while the temperature rise on the surface ΔTsingle(0) is given by:

The difference between ΔTsingle(0) and ΔTsingle-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 oC−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 (R1 = 0.00670 m) and muscle (R1 = 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 R1. Since muscle has the larger value of R1 between the two bulk tissues (muscle vs. brain), the temperature rise decays faster in brain. The differences between the analytical values of ΔTplane-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 R1-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 R1-eff for a single tissue, obtained from a set of 3D FD solutions at different FWHM, should converge to the value of R1 given in eqn (5) as frequency increases. Fig. 5 shows results of the curve-fitted value of R1-eff vs. frequency for a skin-only configuration from 10 to 80 GHz. The value of R1-eff obtained at 80 GHz is 0.00678 m, which is a 1.2% difference from the exact value R1 = 0.00670 m.

Fig. 5
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
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 So. 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 (∆Tplane/Sinc), 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 R1-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
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) D50. 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 D50/2. The value of D50 can be extracted from each solution of the 3D FD solver. Once ∆To and D50 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 D50 was found to be approximately linearly related to the FWHM with a slope that is close to unity (for large FWHM, D50 asymptotically approaches the FWHM). In view of this, a model of the form D50 = FWHM + Do was adopted so that only one parameter, Do, is required to model the width of the surface temperature rise distribution. The intercept Do 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 D50 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.


1D Monte Carlo results—(∆Tplane/Sinc)

Distributions of (∆Tplane/Sinc) 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
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 (∆Tplane/Sinc) 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, fc is the high-pass corner frequency, while Ap and Bp 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 oC−1, Tair = 22oC) 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
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 oC−1) and convective (h = 10 W m−2 oC−1, T air = 22oC) 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 oC−1).

3D Monte Carlo results—R1-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 R1-eff ranged from 0.011 m (3-tissue, adiabatic, 10 GHz) to 0.0076 m (4-tissue, convective, h = 10 W m−2 oC−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 R1-eff vs. frequency for the 3- and 4-tissue (adiabatic and convective, h = 10 W m−2 oC−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
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 oC−1, T air = 22oC) conditions.

where Fm is a different constant depending on whether the configuration is 3- or 4-tissue and fc 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 R1-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 R1-eff are presented in Fig. 10 vs. frequency.

Table 3
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 oC−1 and T air = 22oC.
Fig. 10
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 oC−1, T air = 22oC) conditions. Obtained from the output distributions of 3D Monte Carlo simulations.

3D Monte Carlo results—50% isotherm diameter intercept Do

The distributions of Do 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 Do 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 Do vs. frequency for the adiabatic and convective heat-loss (h = 10 W m−2 oC−1, Tair = 22oC) Monte Carlo simulations along with their fitted high-pass functions.

Fig. 11
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 oC−1 and T air = 22oC.


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 So and inputting the desired value of ∆To and percentile level of (∆Tplane/Sinc)p and R1-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. 12a for the spatial-peak power density that produces 1oC vs. FWHM for the 50th percentile (∆Tplane/Sinc)p and mean R1-eff, both under adiabatic conditions.

Fig. 12
Fig. 12:
Spatial-peak power density of a narrow beam at four selected frequencies that produces 1oC 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 oC−1, T air = 22oC) (∆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 (∆Tplane/Sinc) and mean R1-eff (from the averaged 3- and 4-tissue data) for the convective heat-loss boundary condition: h = 10 W m−2 oC−1 and Tair = 22oC. The results are shown in Fig. 12b. The corresponding curves (at the same frequencies) in Fig. 12a and b are similar due to the fact that the 95th percentile, convective (∆Tplane/Sinc) is close in value to the 50th percentile, adiabatic value. Both sets of curves show that for relatively wide beam widths, the 1oC 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 D50. 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 mm2, 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 gs = 0.601 × HPBW and peak value So, is spatially averaged over a circular area of radius Ra 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, SLimit, then the relationship between SLimit and So is:

The temperature rise for this power density distribution, whose peak power density is defined by eqn (12), will be denoted as ∆To,ave and can be found by solving for So 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 (∆Tplane/Sinc)p is chosen for the frequency of interest:

The averaging area that results in a fixed ∆To,ave (e.g., ∆To,ave = 1oC) can be found by numerically solving eqn (13) for Ra given the assumed HPBW and values of (∆Tplane/Sinc)p and R1-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 ∆To,ave (for spatially averaged power density) is to the temperature rise from a plane wave at the plane-wave limit ∆Tlimit. The latter can be calculated (for the percentile level p) as: ∆Tlimit = SLimit(∆Tplane/Sinc)p. A test ratio can be formed by the quotient ∆To,ave/∆Tlimit, whereby a value of unity is ideal. This ratio is given by:

where it is seen that the test ratio is independent of (∆Tplane/Sinc)p and is only a function of HPBW and R1-eff.

As an example, the test ratio for 28 GHz is plotted in Fig. 13 for several averaging areas (area = πRa2), that include the ICNIRP (1998) recommendation of 2,000 mm2 and the 400 mm2 area proposed in Hashimoto et al. (2017). For this plot, the 3- and 4-tissue-averaged R1-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 mm2 averaging area is too large since its curve never crosses unity. The 400 mm2 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 mm2 area.

Fig. 13
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.


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 39oC (Margerl and Treede 1996; Walters et al. 2004). Since exposure standards are usually designed to restrict temperature rise at or below 1oC, 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 (∆Tplane/Sinc)p and a constant for mean R1-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 oC−1 and an air temperature of 20oC. Comparison of the mean values of (∆Tplane/Sinc) in Fig. 5 of their manuscript and the convective (h = 10 W m−2 oC−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 (∆Tplane/Sinc) of 25% when going from a heat transfer coefficient of 1 to 10 W m−2 oC−1. By comparison, an increase in the heat transfer coefficient from 0 (adiabatic) to 10 W m−2 oC−1 in the current study (with Tair = 22oC) produces an increase of 29% for the 50th percentile (∆Tplane/Sinc) 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 (∆Tplane/Sinc) 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 1oC in all cases. They conclude that the maximal averaging area is frequency dependent with areas ranging from 300 mm2 at 10 GHz to 190 mm2 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.


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 So 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 oC−1 and 22oC 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 1oC. 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 mm2 is far too large, while averaging areas on the order of 300–500 mm2 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.


The authors are grateful for the helpful assistance of James McNamee and Lindsay Beaton in the preparation of this manuscript.


Adibzadeh F, van Rhoon GC, Verduijn GM, Naus-Postema NC, Paulides MM. Absence of acute ocular damage in humans after prolonged exposure to intense RF EMF. Phys Med Biol 61:488–503; 2016. DOI 10.1088/0031-9155/61/2/488.
Alekseev SI, Radzievsky AA, Szabo I, Ziskin MC. Local heating of human skin by millimeter waves: effect of blood flow. Bioelectromagnet 26:489–501; 2005.
Alekseev SI, Ziskin MC. Local heating of human skin by millimeter waves: a kinetics study. Bioelectromagnet 24:571–581; 2003.
Anderson V, Croft R, McIntosh RL. SAR versus Sinc: what is the appropriate RF exposure metric in the range 1–10 GHz? Part 1: using planar body models. Bioelectromagnet 31:454–466; 2010. DOI 10.1002/bem.20578.
Andrews JG, Buzzi S, Choi W, Hanly SV, Lozano A, Soong ACK, Zhang JC. What will 5G be? IEEE J Sel Areas Comm 32:1065–1082; 2014. DOI 10.1109/JSAC.2014.2328098.
Blick DW, Adair ER, Hurt WD, Sherry CJ, Walters TJ, Merritt JH. Thresholds of microwave-evoked warmth sensations in human skin. Bioelectromagnet 18:403–409; 1997.
Drossos A, Santomaa V, Kuster N. The dependence of electromagnetic energy absorption upon human head tissue composition in the frequency range of 300–3000 MHz. IEEE Trans Microw Theory Tech 48:1988–1995; 2000.
Foster KR, Ziskin MC, Balzano Q. Thermal response of human skin to microwave energy: a critical review. Health Phys 111:528–541; 2016.
Hasgall PA, Di Gennaro F, Baumgartner C, Neufeld E, Lloyd B, Gosselin MC, Payne D, Klingenböck A, Kuster N. IT’IS Database for thermal and electromagnetic parameters of biological tissues, version 4.0 [online]. 2018. DOI 10.13099/VIP21000-04-0. Available at Accessed 31 July 2018.
Hashimoto Y, Hirata A, Morimoto R, Aonuma S, Laakso I, Jokela K, Foster KR. On the averaging area for incident power density for human exposure limits at frequencies over 6 GHz. Phys Med Biol 62:3124–3138; 2017.
Hwang K, Kim JH, Baik SH. The thickness of the skull in Korean adults. J Craniofac Surg 10:395–399; 1999.
International Commission on Non-Ionizing Radiation Protection. Guidelines for limiting exposure to time-varying electric, magnetic, and electromagnetic fields (up to 300 GHz). Health Phys 74:494–522; 1998.
Institute of Electrical and Electronics Engineers. Standard for safety levels with respect to human exposure to radio frequency electromagnetic fields, 3 kHz to 300 GHz. Piscataway, NJ: IEEE; IEEE C95.1; 2005.
Kakasheva-Mazenkovska L, Milenkova L, Gjokik G, Janevska V. Variations of the histomorphological characteristics of human skin of different body regions in subjects of different age. Prilozi 32:119–128; 2011.
Kanezaki A, Hirata A, Watanabe S, Shirai H. Parameter variation effects on temperature elevation in a steady-state, one-dimensional thermal model for millimeter wave exposure of one- and three-layer human tissue. Phys Med Biol 55:4647–4659; 2010. DOI 10.1088/0031-9155/55/16/003.
Margerl W, Treede RD. Heat-evoked vasodilatation in human hairy skin: axon reflexes due to low-level activity of nociceptive afferents. J Physiol 497:837–848; 1996.
Mahinda HAM, Murty OP. Variability in thickness of human skull bones and sternum—an autopsy experience. J Forensic Medicine and Toxicology 26:26–31; 2009.
Nelson DA, Walters TJ, Ryan KL, Emerton KB, Hurt WD, Ziriax JM, Johnson LR, Mason PA. Inter-species extrapolation of skin heating resulting from millimeter wave irradiation: modeling and experimental results. Health Phys 84:608–615; 2003.
Neufeld E, Carrasco E, Murbach M, Balzano Q, Christ A, Kuster N. Theoretical and numerical assessment of maximally allowable power-density averaging area for conservative electromagnetic exposure assessment above 6 GHz. Bioelectromagnet 39:617–630; 2018.
Pennes HH. Analysis of tissue and arterial blood temperatures in the resting human forearm. J Appl Physiol 1:93–122; 1948.
Sasaki K, Mizuno M, Wake K, Watanabe S. Monte Carlo simulations of skin exposure to electromagnetic field from 10 GHz to 1 THz. Phys Med Biol 62:6993–7010; 2017.
Walters TJ, Blick DW, Johnson LR, Adair ER, Foster KR. Heating and pain sensation produced in human skin by millimeter waves: comparison to a simple thermal model. Health Phys 78:259–267; 2000.
Walters TJ, Ryan KL, Nelson DA, Blick DW, Mason PA. Effects of blood flow on skin heating induced by millimeter wave irradiation in humans. Health Phys 86:115–120; 2004.
Wu T, Rappaport TS, Collins CM. Safe for generations to come. IEEE Microw Mag 16:65–84; 2015.
Zhadobov M, Alekseev SI, Le Dréan Y, Sauleau R, Fesenko EE. Millimeter waves as a source of selective heating of skin. Bioelectromagnet 36:464–475; 2015. DOI 10.1002/bem.21929.

exposure, radiofrequency; modeling, dose assessment; Monte Carlo; tissue, body

Copyright © 2019 The Author(s). Published by Wolters Kluwer Health, Inc., on behalf of the Health Physics Society.