D ow nloaded from http://journals.lw w .com /health-physics by B hD M f5eP H K av1zE oum 1tQ fN 4a+ kJLhE Z gbsIH o4X M i 0hC yw C X 1A W nY Q p/IlQ rH D 3i3D 0O dR yi7T vS F l4C f3V C 4/O A V pD D a8K K G K V 0Y m y+ 78= on 03/15/2024 P 1Non-ionizing Radiation Physical Sciences Division, Consumer a Clinical Radiation Protection Bureau, Health Canada. The authors declare no conflicts of interest. For correspondence contact Gregory Gajda, Consumer and Clinic Radiation Protection Bureau, Health Canada, 775 Brookfield Road, Ottaw Ontario, Canada, K1A 1C1, or email at greggajda99@gmail.com. (Manuscript accepted 21 October 2018) 0017-9078/19/0 Copyright © 2019 The Author(s). Published by Wolters Kluwer Heal Inc., on behalf of the Health Physics Society. This is an open-access article d tributed under the terms of the Creative Commons Attribution-Non Commerci No Derivatives License 4.0 (CCBY-NC-ND), where it is permissible to downlo and share the work provided it is properly cited. The work cannot be chang in any way or used commercially without permission from the Journal. DOI: 10.1097/HP.0000000000001036 254 aper nd al a, th, is- al- ad ed MODEL OF STEADY-STATE TEMPERATURE RISE IN MULTILAYER TISSUES DUE TO NARROW-BEAM MILLIMETER-WAVE RADIOFREQUENCY FIELD EXPOSURE Gregory B. Gajda, Eric Lemay, and Jonathan Paradis1 Abstract—The assessment of health effects due to localized expo- sures from radiofrequency fields is facilitated by characterizing the steady-state, surface temperature rise in tissue. A closed-form analytical model was developed that relates the steady-state, sur- face temperature rise in multilayer planar tissues as a function of the spatial-peak power density and beam dimensions of an incident millimeter wave. Model data was derived from finite-difference so- lutions of the Pennes bioheat transfer equation for both normal- incidence plane waves and for narrow, circularly symmetric beams with Gaussian intensity distribution on the surface. Monte Carlo techniques were employed by representing tissue layer thicknesses at different body sites as statistical distributions compiled from hu- man data found in the literature. The finite-difference solutions were validated against analytical solutions of the bioheat equation for the planewave case andagainst anarrow-beamsolutionperformed using a commercial multiphysics simulation package. In both cases, agreement was within 1–2%. For a given frequency, the resulting ana- lyticalmodel has four input parameters, two of which are determin- istic, describing the level of exposure (i.e., the spatial-peak power density and beamwidth). The remaining two are stochastic quanti- ties, extracted from theMonteCarlo analyses. The analyticalmodel is composed of relatively simple functions that can be programmed in a spreadsheet. Demonstration of the analytical model is provided in two examples: the calculation of spatial-peak power density vs. beam width that produces a predefined maximum steady-state surface temperature, and the performance evaluation of various proposed spatial-averaging areas for the incident power density. Health Phys. 117(3):254–266; 2019 Key words: exposure, radiofrequency; modeling, dose assessment; Monte Carlo; tissue, body INTRODUCTION HUMAN EXPOSURE limits for the millimeter (mm)-wave radio- frequency (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 maxi- mum power density (ICNIRP 1998; IEEE 2005). This is due to the highly superficial nature of mm-wave energy deposi- tion, 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 tempo- ral steady-state values, under the assumption that the expo- sure 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 aver- aged 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, expo- sures from such devices are likely to be highly localized. In ad- dition to the localized nature of the exposure, the effect of tissue curvature takes on less importance as thewavelength 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 de- creases 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 www.health-physics.com http://creativecommons.org/licenses/by-nc-nd/4.0/ http://creativecommons.org/licenses/by-nc-nd/4.0/ 255Model of steady-state temperature rise in tissue c G.B. GAJDA ET AL. D ow nloaded from http://journals.lw w .com /health-physics by B hD M f5eP H K av1zE oum 1tQ fN 4a+ kJLhE Z gbsIH o4X M i 0hC yw C X 1A W nY Q p/IlQ rH D 3i3D 0O dR yi7T vS F l4C f3V C 4/O A V pD D a8K K G K V 0Y m y+ 78= on 03/15/2024 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 calcu- lated through solution of the Pennes bioheat transfer equa- tion (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 solu- tions of a cylindrically symmetric, two-dimensional formu- lation 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 dif- ferent individuals and body sites. Anderson et al. (2010) used a Monte Carlo approach that incorporated statistical tissue thickness variations in a one-dimensional formula- tion 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 calcu- lated 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 tem- perature 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 symmet- ric version (henceforth referred to as a 3D solution). Finally, the results of theMonte 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 www.health-phy 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). 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 cho- sen such that the calculated temperature rise approximately decays to zero, which would hold true if the tissue extended to infinity. 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 pro- portional 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 coordi- nates (i.e., having variation in the axial z and radial r direc- tions) is: k 1 r ∂U ∂r þ ∂2U ∂r2 � � þ ∂U ∂z � � ∂k ∂z � � þ k ∂2U ∂z2 −rbCb rmbU þ rSAR r; zð Þ ¼ 0; ð1Þ 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 be- tween 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; r, the mass density in units of sics.com Table 1. Symbols, units, and quantities. Symbol Units Quantity S(r) W m−2 Circularly symmetric projected power density distribution HPBW m Span over which the projected power density falls to one-half the maximum value SAR(r, z) W kg−1 Circularly symmetric SAR distribution FWHM m Span over which the transverse SAR distribution falls to one-half the peak intensity Sinc W m−2 Plane-wave power density incident normally on the tissue surface SAR(z) W kg−1 Axial SAR distribution due to a plane wave T(r,z) oC Distribution of temperature in the tissue with respect to radial and axial coordinates Tblood oC Temperature of blood—fixed at 37oC U(r,z) oC Temperature difference between T(r,z) and Tblood mb m3 kg−1 s−1 Blood perfusion rate specific to each tissue r kg m−3 Mass density specific to each tissue k W m−1 oC−1 Heat conductivity specific to each tissue Cb J kg−1 oC−1 Specific heat capacity of blood rb kg m−3 Mass density of blood Tair oC Ambient temperature of the air on the surface of the tissue h W m−2 oC−1 Surface convection heat transfer coefficient Tnull(z) oC Distribution of temperature in the absence of RF exposure (1D and 3D cases) Unull(z) oC Difference between Tnull(z) and Tblood Tfull(z) oC Distribution of temperature with full SAR due to a plane wave (1D case) Ufull(z) oC Difference between Tfull(z) and Tblood ΔT1D(z) oC Distribution of temperature rise due to a plane wave (1D case) (= Ufull(z) − Unull(z)) ΔT3D(r,z) oC Distribution of temperature rise due to a narrow beam (3D case) (= Ufull(r,z) − Unull(z)) g m Gaussian width parameter of the 3D SAR distribution (= 0.601 FWHM) So W m−2 Maximum intensity in the center of a narrow, projected beam of power density DTo oC Maximum temperature rise in multilayer tissues due to a narrow beam (ΔTplane/Sinc) oC m2 W−1 Maximal steady-state temperature rise per unit incident plane-wave power density R1−eff m Effective diffusion length for multilayer tissues R1 m Diffusion length for a single (homogeneous) tissue d m Depth at which SAR(z) drops to exp{�1} TR Fraction of incident power density that is transmitted across the air-tissue interface DTsingle-max oC Temperature rise maximum for a single tissue under plane-wave exposure DTsingle(0) oC Temperature rise on the surface of a single tissue under plane-wave exposure D50 m Diameter of the 50% surface temperature rise isotherm Do m Intercept in the linear model of D50 vs. FWHM fc GHz Corner frequency in high-pass models of (ΔTplane/Sinc) and R1−eff vs. frequency Ap and Bp Slope and intercept of the linear term of (ΔTplane/Sinc) at percentile level p, respectively Fm m Constant multiplier in the model of R1−eff vs. frequency gs m Gaussian width parameter of the projected power density distribution (= 0.601 HPBW) SLimit W m−2 Plane-wave power density guideline limit ΔTo,ave oC Maximum temperature rise for a spatially averaged narrow-beam at the plane-wave limit ΔTlimit oC Maximum temperature rise from a plane wave at the plane-wave limit 256 Health Physics September 2019, Volume 117, Number 3 D ow nloaded from http://journals.lw w .com /health-physics by B hD M f5eP H K av1zE oum 1tQ fN 4a+ kJLhE Z gbsIH o4X M i 0hC yw C X 1A W nY Q p/IlQ rH D 3i3D 0O dR yi7T vS F l4C f3V C 4/O A V pD D a8K K G K V 0Y m y+ 78= on 03/15/2024 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 rb, the mass density of blood. Equation (1) is solved subject to the boundary condi- tion on the surface −k dT dz jz¼0 ¼ −h T r; 0ð Þ−Tair½ �; ð2Þ where Tair is the ambient temperature of the surrounding air www.health-phy (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 Wm−2 oC−1. The term on the left-hand side of eqn (2) repre- sents the rate of heat flow away from the surface. If insulat- ing 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. sics.com Fig. 1. Definitions of half-power projected beam width (HPBW) of the projected power density distribution and fullwidth 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 gs related to FWHM and HPBW, respectively. Fig. 2. Geometry of the electromagnetic/thermal problem for a 3-tissue configuration. Transverse electric and magnetic (TEM) waves are nor- mally incident on the tissue half-space with incident (unperturbed), plane-wave power density Sinc. 257Model of steady-state temperature rise in tissue c G.B. GAJDA ET AL. D ow nloaded from http://journals.lw w .com /health-physics by B hD M f5eP H K av1zE oum 1tQ fN 4a+ kJLhE Z gbsIH o4X M i 0hC yw C X 1A W nY Q p/IlQ rH D 3i3D 0O dR yi7T vS F l4C f3V C 4/O A V pD D a8K K G K V 0Y m y+ 78= on 03/15/2024 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 expo- sure), the differential termswith respect to r are dropped and only the axial SAR distribution (i.e., SAR(z)) is used. For this case, the solutionU(z) is a function only of the axial co- ordinate 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 distribu- tion is independent of r) and serves as the boundary condi- tion as r → ∞. The net temperature rise for the 3D case ΔT3D(r,z) is equal to the difference between the 3D solution at full SARUfull(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): SAR r; zð Þ ¼ SAR zð Þ exp �r2=g2 � � ; ð3Þ 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. www.health-phy Analytical model for steady-state temperature rise from a narrow, projected beam Examples of distributions of steady-state surface tem- perature rise for a 4-layer tissue (under adiabatic condi- tions), 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 so- lutions of eqn (1), which will be dealt with later. In the fig- ure, 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. 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) solu- tion extended to the case of multilayer tissues, SAR deposi- tions that penetrate into the tissue, and for convective surface boundary conditions. The maximum temperature rise (in the center of the beam)DTo in the multilayer tissue is modeled as: DTo ¼ So DTplane Sinc � � ffiffiffiffi p p X exp X 2 � � erfc Xf g ; X ¼ 0:601FWHM 2R1−eff ; ð4Þ 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): R1 ¼ ffiffiffi k pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi rmb rb Cb p ; ð5Þ where k, r, andmb are constant throughout the single tissue. sics.com Fig. 4. Maximum, normalized, steady-state, temperature rise ΔTo/So normalized to (ΔTplane/Sinc) as a function of the FWHM of a narrow, projected beamwith Gaussian SAR distribution (from eqn 4) for three different values of effective diffusion length R1-eff. Fig. 3. Examples of steady-state, surface temperature rise distribu- tions for different FWHM at the same spatial-peak power density So (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 maxi- mum surface temperature rise, ΔTo = ΔT3D(0,0), occurs in the center of the distribution (i.e., at r = 0). 258 Health Physics September 2019, Volume 117, Number 3 D ow nloaded from http://journals.lw w .com /health-physics by B hD M f5eP H K av1zE oum 1tQ fN 4a+ kJLhE Z gbsIH o4X M i 0hC yw C X 1A W nY Q p/IlQ rH D 3i3D 0O dR yi7T vS F l4C f3V C 4/O A V pD D a8K K G K V 0Y m y+ 78= on 03/15/2024 The effect that the parameter R1-eff has on the maximal temperature-rise response vs. beamwidth is shown in Fig. 4. As seen from this figure, the value of R1-eff controls the rate at which DTo/So asymptotically approaches DTplane/Sinc (or for the casewhere So = Sinc, where the narrow-beam temper- ature rise approaches that of a plane wave). This parameter is determined by using a procedure that takes values of DTo/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). The four parameters in the analytical model of eqn (4) can be divided into two categories. Parameters So and FWHM are deterministic, electromagnetic quantities ob- tained from the known exposure conditions (e.g., for a given source and SAR distribution). The remaining two parame- ters, (ΔTplane/Sinc) and R1-eff, are thermally based quantities that depend on frequency and tissue configuration (i.e. com- position 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 de- rive 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) www.health-phy backed by a bulk layer of brain (grey matter). The latter con- figuration was intended to represent head tissues (excluding the eye) while the former represents most other body sur- face 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 in- creased 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 re- laxed 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 direc- tion 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 minimiza- tion 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 tis- sue configuration, a single determination of (ΔTplane/Sinc) (from the 1D BHTE) and four determinations of DTo/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 1DBHTE to obtain the null and full temperature distributions and subsequently, the maximum net temperature rise ΔTplane. An iteration of sics.com 259Model of steady-state temperature rise in tissue c G.B. GAJDA ET AL. D ow nloaded from http://journals.lw w .com /health-physics by B hD M f5eP H K av1zE oum 1tQ fN 4a+ kJLhE Z gbsIH o4X M i 0hC yw C X 1A W nY Q p/IlQ rH D 3i3D 0O dR yi7T vS F l4C f3V C 4/O A V pD D a8K K G K V 0Y m y+ 78= on 03/15/2024 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. Thiswas then followed by a curve fitting to ob- tain 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 stan- dard 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). 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 A1,2 exp{±z/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): SAR zð Þ ¼ SARoexp �z=df g SincTR rd exp �z=df g; ð6Þ where d 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/d} 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: T zð Þ ¼ SincTRR1 k 1−d2=R1 2 ½ k=R1 þ hd=R1 k=R1 þ h � � exp ( −z=R1g − d R1 � � exp ( −z=dg� þ Tblood−h Tblood−Tairð Þ k=R1 þ hð Þ exp ( −z=R1g : ð7Þ 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 (in- dependent of power density). Eqn (7) illustrates the signifi- cance of the diffusion length R1, which also governs the rate of decay of temperature rise (for d < R1) in the deep tissues. For the adiabatic case (h = 0), the maximum tempera- ture rise occurs at the surface while for the convective heat-loss case, it occurs a very short distance inside the www.health-phy surface of the tissue. It can be shown that the single-tissue temperature rise maximum DTsingle-max at this location is given by: DTsingle‐max ¼ SincTRR1 k � � 1 ½1þ d=R1ð Þ� k þ hd k þ hR1 � �a ; a ¼ 1−d=R1ð Þ−1 ð8Þ while the temperature rise on the surfaceDTsingle(0) is given by: DTsingle 0ð Þ ¼ SincTRR1 k � � 1 ½1þðd R1= Þ� 1þ h R1=kð Þ½ � : ð9Þ The difference between DTsingle(0) and DTsingle-max is small; for example, being on the order of 0.2%, 0.6%, and 1.1% for values of h = 10, 20, and 30Wm−2 oC−1, respectively, and for skin tissue at 20GHz. From eqn (9) it is seen that the convective heat-loss condition produces smaller tempera- ture 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 ofDTplane-max given by eqn (8) and the equivalent 1D FD re- sults ranged from less than 0.05% at 10 GHz to 1% at 80 GHz for both skin and muscle configurations and for ei- ther 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 adia- batic 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, ob- tained from a set of 3D FD solutions at different FWHM, should converge to the value of R1 given in eqn (5) as fre- quency 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. 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) electro- magnetic and thermodynamics solvers. Fig. 6 shows the re- sults 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. sics.com Fig. 5. Curve-fitted values of R1-eff vs. frequency, obtained from nar- row-beam, 3D FD solutions for a skin-only tissue. 260 Health Physics September 2019, Volume 117, Number 3 D ow nloaded from http://journals.lw w .com /health-physics by B hD M f5eP H K av1zE oum 1tQ fN 4a+ kJLhE Z gbsIH o4X M i 0hC yw C X 1A W nY Q p/IlQ rH D 3i3D 0O dR yi7T vS F l4C f3V C 4/O A V pD D a8K K G K V 0Y m y+ 78= on 03/15/2024 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 den- sity 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 re- sulted 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 corre- spond 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 agree- ment (approximately 2% difference) with the calculated 1D 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). www.health-phy 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 SARdata) are shown inFig. 7. These are compared to analytical Gaussian distributions that give the best fit. From the figure, it is seen that the transverse SARdistributions 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. 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 tem- perature 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 tem- perature 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 dis- tribution 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 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). sics.com 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. 261Model of steady-state temperature rise in tissue c G.B. GAJDA ET AL. D ow nloaded from http://journals.lw w .com /health-physics by B hD M f5eP H K av1zE oum 1tQ fN 4a+ kJLhE Z gbsIH o4X M i 0hC yw C X 1A W nY Q p/IlQ rH D 3i3D 0O dR yi7T vS F l4C f3V C 4/O A V pD D a8K K G K V 0Y m y+ 78= on 03/15/2024 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 tomodel thewidth of the surface temperature rise distribution. The intercept Do is easily computed from the 3D FD output data using linear re- gression, 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 thick- nesses 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 be- tween 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 lognor- mal 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 re- ported mean thicknesses at six sites for males and females from 76 individuals. The latter study reportedmean 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—(Δ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 skew- ness 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 www.health-phy 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). 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 lin- ear term and a high-pass function given by: DTplane=Sinc p ¼ Ap f þ Bp ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1þ fc=fð Þ2 q ; ð10Þ 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 dif- ference 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. 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). sics.com Table 2. Tabulated values of the parameters Ap and Bp and fc used in eqn (10) for calculating (ΔTplane/Sinc)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, Tair = 22oC) boundary conditions. Adiabatic Convective Percentile level Ap (oC W−1 m2 GHz−1) Bp (oC W−1 m2) fc (GHz) Ap (oC W−1 m2 GHz−1) Bp (oC W−1 m2) fc (GHz) 95 7.50 � 10−5 0.0170 10.8 5.85 � 10−5 0.0124 9.5 90 7.01 � 10−5 0.0165 10.8 5.60 � 10−5 0.0121 9.5 80 6.59 � 10−5 0.0157 10.8 5.36 � 10−5 0.0117 9.8 70 6.41 � 10−5 0.0151 10.8 5.08 � 10−5 0.0114 10.3 60 6.22 � 10−5 0.0145 10.8 4.88 � 10−5 0.0112 10.6 50 6.04 � 10−5 0.0141 10.8 4.72 � 10−5 0.0109 10.8 262 Health Physics September 2019, Volume 117, Number 3 D ow nloaded from http://journals.lw w .com /health-physics by B hD M f5eP H K av1zE oum 1tQ fN 4a+ kJLhE Z gbsIH o4X M i 0hC yw C X 1A W nY Q p/IlQ rH D 3i3D 0O dR yi7T vS F l4C f3V C 4/O A V pD D a8K K G K V 0Y m y+ 78= on 03/15/2024 3D Monte Carlo results—R1-eff A total of 500 iterations were performed at each fre- quency 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. Dif- ferences 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 configu- rations 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 R1−eff mean ¼ Fm ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1þ fc=fð Þ2 q ; ð11Þ where Fm is a different constant depending on whether the configuration is 3- or 4-tissue and fc is the corner frequency. Fig. 9. Mean R1-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, Tair = 22oC) conditions. www.health-phy 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 meanR1-eff data to eqn (11) was carried out using nonlin- ear 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. 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 un- der 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 bound- ary condition. The results suggest that the value of the inter- cept is relatively insensitive to the skin and SAT thickness. Differences of the mean Do between 3- and 4-tissue config- urations, for the same boundary conditions, were under 7%. For simplicity, the mean values for the 3- and 4-tissue con- figurations 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. Table 3. Constants used for calculating the mean value of R1-eff using eqn (11), where the result, R1-eff, is in m and f is in GHz. Convective boundary condition has parameters h = 10 W m�2 oC�1 and Tair = 22oC. Tissue configuration Boundary condition Fm (m) fc (GHz) 3-tissue Adiabatic 9.40 � 10−3 6.4 4-tissue Adiabatic 8.87 � 10−3 5.6 Average of 3- and 4-tissue R1−eff data Adiabatic 9.14 � 10−3 6.0 3-tissue Convective 8.04 � 10−3 7.1 4-tissue Convective 7.69 � 10−3 6.3 Average of 3- and 4-tissue R1−eff data Convective 7.86 � 10−3 6.8 sics.com Fig. 10. SD of R1-eff vs. frequency for the 3- and 4-tissue configura- tions under adiabatic and convective heat-loss (h = 10 W m−2 oC−1, Tair = 22oC) conditions. Obtained from the output distributions of 3D Monte Carlo simulations. Fig. 11. Mean values of Do vs. frequency from 3DMonte Carlo simu- lation. 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 = 10Wm−2 oC−1 and Tair = 22oC. 263Model of steady-state temperature rise in tissue c G.B. GAJDA ET AL. D ow nloaded from http://journals.lw w .com /health-physics by B hD M f5eP H K av1zE oum 1tQ fN 4a+ kJLhE Z gbsIH o4X M i 0hC yw C X 1A W nY Q p/IlQ rH D 3i3D 0O dR yi7T vS F l4C f3V C 4/O A V pD D a8K K G K V 0Y m y+ 78= on 03/15/2024 Applications The development of a closed-form model for predic- ting narrow-beam temperature rise permits straightforward calculation of any one of the model parameters from knowl- edge of the others. For instance, the spatial-peak power den- sity 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 de- sired 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), respec- tively. 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. A similar calculation was performed for the case of the 95th percentile (ΔTplane/Sinc) and mean R1-eff (from the aver- aged 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 sim- ilar due to the fact that the 95th percentile, convective (ΔTplane/Sinc) is close in value to the 50th percentile, adi- abatic value. Both sets of curves show that for relatively wide beam widths, the 1oC spatial-peak power density as- ymptotically 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 Inter- national 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. Addi- tionally, the simple nature of the model lends itself to inves- tigations of spatial-averaging areas for power density in the development of exposure standards. For instance, Hashi- moto et al. (2017) propose an averaging area of 400 mm2, www.health-phy 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�HPBWand peak value So, is spatially averaged over a circular area of radius Ra by in- tegrating 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 relation- ship between SLimit and So is: SLimit ¼ So gs2 Ra 2 � � 1−exp −Ra 2=gs 2 � �� � ; gs ¼ 0:601HPBW : ð12Þ The temperature rise for this power density distribution, whose peak power density is defined by eqn (12), will be de- noted 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 rela- tionship can be assumed whereby FWHM/HPBW is a con- stant 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: DTo;ave ¼ SLimit Ra 2=gs2 1−exp −Ra 2=gs2 � �� � DTplane Sinc � � p ffiffiffiffi p p X exp X 2 � � erfc Xf g; X ¼ 0:601KHPBW 2R1−eff : ð13Þ 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. sics.com 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 (ΔTplane/Sinc) and mean R1-eff (average of 3- and 4-tissue, adiabatic) and (b) the 95th percentile, convective (h = 10 W m−2 oC−1, Tair = 22oC) (ΔTplane/Sinc) and mean R1-eff (average of 3- and 4-tissue, convective). 264 Health Physics September 2019, Volume 117, Number 3 D ow nloaded from http://journals.lw w .com /health-physics by B hD M f5eP H K av1zE oum 1tQ fN 4a+ kJLhE Z gbsIH o4X M i 0hC yw C X 1A W nY Q p/IlQ rH D 3i3D 0O dR yi7T vS F l4C f3V C 4/O A V pD D a8K K G K V 0Y m y+ 78= on 03/15/2024 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 spa- tially 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: DTo;ave DTlimit ¼ Ra 2 gs2 � � ffiffiffiffi p p X exp X 2 � � erfc Xf g 1−exp −Ra 2=gs2 � �� � ; X ¼ 0:601KHPBW 2R1−eff ; ð14Þ 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 =pRa 2), 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 se- lected, and it is assumed the power density and SAR distri- bution 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, show- ing 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. Plot of the test ratio ΔTo,ave/ΔTplane vs. HPBW for different spatial-averaging areas where ΔTo,ave is the temperature rise for a nar- row beam with spatially averaged power density at the plane-wave limit and ΔTplane is the temperature rise for a plane wave at the limit. An assumed value of R1-eff = 0.00935 m was used corresponding to the mean value of the 3- and 4-tissue average under adiabatic condi- tions 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 expo- sure 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 www.health-phy 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 vasodi- lation 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 tempera- ture 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 re- strict 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 sics.com 265Model of steady-state temperature rise in tissue c G.B. GAJDA ET AL. D ow nloaded from http://journals.lw w .com /health-physics by B hD M f5eP H K av1zE oum 1tQ fN 4a+ kJLhE Z gbsIH o4X M i 0hC yw C X 1A W nY Q p/IlQ rH D 3i3D 0O dR yi7T vS F l4C f3V C 4/O A V pD D a8K K G K V 0Y m y+ 78= on 03/15/2024 limit the applicability of the analytical model for higher fre- quencies; 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 HPBWof 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 dis- tance. It was found that there was an approximate linear re- lationship (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 eas- ily be modified by substituting 0.82 HPBW for the param- eter 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 epi- dermis, 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 tem- perature 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 10Wm−2 oC−1. By comparison, an increase in the heat trans- fer coefficient from 0 (adiabatic) to 10Wm−2 oC−1 in the cur- rent 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 inde- pendently measured the dielectric parameters of dermis, SAT, and muscle and found large differences between the measured SAT parameters and those reported in the data- base 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 af- fected by tissue thickness variations than by variations in the dielectric parameters. They further noted that the results www.health-phy 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 sev- eral 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 dis- tribution 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 ex- plicitly relate narrow-beam temperature rise to SAR distri- bution 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 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 be- tween the two can be ascertained from numerical simulation of the actual source of the power density. The modelwas developed for two tissue configurations representing head (4-tissue, excluding the eye) and general body sites (3-tissue) and for two surface boundary condi- tions, adiabatic and convective heat loss (heat-loss coeffi- cient 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 solu- tions of the 1D and 3D BHTE. From the resulting Monte Carlo output distributions, statistics of the two additional in- put parameters in the analytical modelwere 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 solu- tions to the 1D and 3D BHTE were compared against known sics.com 266 Health Physics September 2019, Volume 117, Number 3 D ow nloaded from http://journals.lw w .com /health-physics by B hD M f5eP H K av1zE oum 1tQ fN 4a+ kJLhE Z gbsIH o4X M i 0hC yw C X 1A W nY Q p/IlQ rH D 3i3D 0O dR yi7T vS F l4C f3V C 4/O A V pD D a8K K G K V 0Y m y+ 78= on 03/15/2024 analytical solutions and results from a commercial software package (Sim4Life). The differenceswere typically no greater than 1–2%.Comparison of the 50th percentile 1D resultswith 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 tem- perature 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. beamwidth that produces a predefined maximum surface temperature rise. It was shown that the current ICNIRP occupational limit for planewaves or large-area exposure produces temperature rise below 1oC. A second example evaluated the effect of spatial- averaging area on themaximum steady-state temperature rise in relation to that for a plane wave at the plane-wave power density limit. Quantitative comparison between different av- eraging 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 devel- opment of future exposure standards. Acknowledgments—The authors are grateful for the helpful assistance of James McNamee and Lindsay Beaton in the preparation of this manuscript. REFERENCES Adibzadeh F, van Rhoon GC, Verduijn GM, Naus-Postema NC, Paulides MM. Absence of acute ocular damage in humans af- ter 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, ZiskinMC. Local heating of human skin bymillimeter waves: a kinetics study. Bioelectromagnet 24:571–581; 2003. Anderson V, Croft R, McIntosh RL. SAR versus Sinc: what is the appropriate RFexposure metric in the range 1–10GHz? Part 1: using planar body models. Bioelectromagnet 31:454–466; 2010. DOI 10.1002/bem.20578. Andrews JG, Buzzi S, ChoiW, 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 electromag- netic energy absorption upon human head tissue composition in the frequency range of 300–3000MHz. IEEE TransMicrow Theory Tech 48:1988–1995; 2000. www.health-phy 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 Data- base for thermal and electromagnetic parameters of biological tissues, version 4.0 [online]. 2018. DOI 10.13099/VIP21000- 04-0. Available at https://itis.swiss/virtual-population/tissue- properties/database/. 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, mag- netic, 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 fre- quency 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 nocicep- tive 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 Medi- cine 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 frommillimeter wave irradiation: model- ing and experimental results. Health Phys 84:608–615; 2003. Neufeld E, Carrasco E,MurbachM,BalzanoQ, 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 simula- tions 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 irra- diation 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. ■■ sics.com https://itis.swiss/virtual-population/tissue-properties/database/ https://itis.swiss/virtual-population/tissue-properties/database/