Mouhamadou Fall, Serigne Bira Gueye
Department of Physics, Faculty of Sciences and Techniques Cheikh Anta Diop University, Dakar, Senegal.
DOI: 10.4236/jemaa.2026.189009   PDF    HTML   XML   2 Downloads   29 Views   Citations

Abstract

This paper presents a two-dimensional finite-difference time-domain (FDTD) investigation of the local specific absorption rate (SAR) distribution in a heterogeneous human brain model derived from a segmented magnetic resonance imaging (MRI) dataset. The electromagnetic problem is solved at 900 MHz for a normally incident plane wave using a total-field/scattered-field (TF/SF) formulation. A uniaxial perfectly matched layer (UPML) is employed to truncate the computational domain and reduce artificial boundary reflections. The anatomical model is derived from the BrainWeb database and includes gray matter, white matter, cerebrospinal fluid (CSF), and skull. Published dielectric-property data at 900 MHz are assigned to the different tissue classes. The numerical model uses a spatial resolution of 1 mm, a 280 × 290 computational grid, a 40-cell absorbing layer, and 15,666 time steps. The calculated SAR distribution is strongly non-uniform as a consequence of the dielectric contrasts between tissues, multiple reflections at tissue interfaces, and local electromagnetic-field interference. The CSF exhibits the largest mean and maximum SAR values among the considered tissues, whereas the skull exhibits the lowest mean SAR. The tissue-averaged absorption is quantitatively interpreted through the decomposition $\text{SAR}=\left(\sigma /\rho \right)\langle {\left|E\right|}^2\rangle$ , which separates a material factor determined by the tissue properties from a field factor governed by the local electromagnetic environment. Tissue-wise statistical indicators are reported in terms of the mean, standard deviation, and maximum SAR. A systematic verification strategy is established using one-dimensional benchmarks for numerical dispersion, dielectric interfaces, and parallel dielectric slabs, together with a two-dimensional benchmark based on the Mie-series solution for a lossy dielectric cylinder. These tests provide a quantitative assessment of the numerical accuracy of the FDTD solver and establish the reliability of the computed electromagnetic fields used for the anatomical configuration. The results demonstrate the relevance of anatomically heterogeneous models for investigating local electromagnetic absorption in brain tissues. The present two-dimensional framework provides a controlled and quantitatively verified basis for future three-dimensional electromagnetic and electromagnetic-thermal dosimetric studies.

Keywords

FDTD, Specific Absorption Rate, SAR, Electromagnetic Dosimetry, Human Brain, MRI, BrainWeb, Biological Tissues, 900 MHz, UPML

Share and Cite:

1. Introduction

The widespread deployment of wireless communication systems has increased human exposure to radiofrequency (RF) electromagnetic fields, with the head being a particularly relevant anatomical region. This exposure has motivated numerous numerical and experimental studies aimed at characterizing the electromagnetic energy absorbed by biological tissues.

The specific absorption rate (SAR) is one of the principal physical quantities used in electromagnetic dosimetry. It characterizes the rate at which electromagnetic energy is absorbed per unit mass of biological tissue. Accurate prediction of SAR is therefore essential for assessing RF exposure and for investigating spatially localized electromagnetic absorption, and provides the physical basis for regulatory exposure limits [1] [2].

The finite-difference time-domain (FDTD) method is particularly suitable for this class of problems because Maxwell’s equations are solved directly in the time domain, while heterogeneous and lossy media can be represented naturally [3] [4]. The electric field is obtained throughout the computational domain, allowing the local SAR to be evaluated directly from the simulated electromagnetic fields.

The accuracy and physical relevance of numerical dosimetry depend strongly on the anatomical representation. Homogeneous and multilayer head models are useful for studying fundamental electromagnetic phenomena, but they cannot fully reproduce the geometrical complexity and dielectric heterogeneity of anatomical structures. Several studies have therefore investigated SAR in simplified and heterogeneous head models using FDTD techniques [5]-[8].

The influence of tissue dielectric properties on RF absorption has also been investigated. Husni et al. [9] examined the influence of dielectric properties on RF absorption in human-head models, while Wiart et al. [10] investigated differences in RF exposure between children and adults. Dimbylow and Gandhi [11] demonstrated the importance of heterogeneous anatomical models in FDTD calculations of SAR. These studies highlight the need to account for tissue-dependent electromagnetic properties when analyzing local absorption.

The objective of the present work is to develop and investigate a two-dimensional FDTD model of a heterogeneous human brain based on a segmented MRI-derived anatomical dataset. The model considers gray matter, white matter, cerebrospinal fluid (CSF), and skull and is exposed to a normally incident plane wave at 900 MHz. Particular attention is given to the spatial distribution and statistical characteristics of the local SAR, as well as to the differences in absorption among the considered tissue classes.

A further objective is to provide a physical interpretation of the tissue-averaged absorption through the decomposition

$\text{SAR}=\frac{\sigma }{\rho }{\left|E\right|}^2,$

which separates the contribution associated with the intrinsic material properties of each tissue from that associated with the local electromagnetic field. This decomposition provides a direct framework for distinguishing material-dependent and field-dependent effects in heterogeneous anatomical media.

Finally, a systematic verification procedure is established before interpreting the anatomical results. The numerical implementation is assessed using one-dimensional benchmarks for numerical dispersion, dielectric interfaces, and parallel dielectric slabs, together with a two-dimensional benchmark based on the Mie-series solution for a lossy dielectric cylinder. The purpose of this verification strategy is to quantify the numerical accuracy of the FDTD solver and to establish a controlled basis for the subsequent analysis of SAR in the heterogeneous brain model.

2. Anatomical Model and Electromagnetic Properties

2.1. BrainWeb Anatomical Model

The anatomical model used in this study is derived from the BrainWeb Normal Anatomical Brain Model (Version 1.0) database, which provides simulated three-dimensional MRI brain datasets with an isotropic spatial resolution of 1 mm [12]. The volume considered in the present study has dimensions 181 × 217 × 181 mm³.

A single coronal slice extracted at the exact middle index $Y=100$ along the anterior-posterior dimension (corresponding to an anatomical position of approximately 100 mm with respect to the posterior boundary of the volume) is used to construct the two-dimensional computational model. The original BrainWeb segmentation contains several anatomical and non-anatomical tissue labels. To construct the present electromagnetic model, a crisp segmentation thresholding criterion was applied, and four principal biological tissue classes were retained:

• gray matter (GM),

• white matter (WM),

• cerebrospinal fluid (CSF),

• skull.

All other anatomical labels present in the original dataset (such as skin, fat, muscle, and glial tissue), as well as internal cavities or nasal sinuses, were excluded from the computational domain and treated as background air.

The spatial step of the FDTD grid is chosen as $\Delta x=\Delta y=1$ mm, matching the native isotropic resolution of the BrainWeb dataset. This choice avoids introducing additional spatial interpolation into the segmented anatomical map and provides a direct correspondence between the anatomical labels and the FDTD computational cells. A finer FDTD grid could improve the numerical representation of electromagnetic field variations and tissue interfaces, but would not introduce additional anatomical information beyond that contained in the original 1-mm segmented dataset.

The brain is surrounded by an air region to provide sufficient space between the anatomical structure and the absorbing boundary. An absorbing layer is applied at the outer boundary to emulate an open electromagnetic domain and to minimize artificial reflections. Specifically, a minimum air-buffer distance of 15 cells (15 mm) is maintained between the outermost skull boundary and the total-field/scattered-field (TF/SF) framework. The computational arrangement is schematically represented in Figure 1, showing the brain region, surrounding air, the location of the TF/SF boundary, and the bounding uniaxial perfectly matched layer (UPML).

Figure 1. Schematic representation of the computational domain.

2.2. Electromagnetic Properties of the Tissues

The electromagnetic properties of biological tissues are frequency dependent. At 900 MHz, the relative permittivity, electrical conductivity, and mass density adopted in this study are listed in Table 1. The dielectric data are taken from the published database of Gabriel et al. [13].

Table 1. Electromagnetic properties of the tissues considered in the model at 900 MHz.

Tissue	${\epsilon }_r$	$\sigma$ (S/m)	$\rho$ (kg/m3)
White matter	38.9	0.59	1041
Gray matter	52.7	0.94	1045
CSF	68.6	2.41	1007
Skull	12.45	0.14	1908

The strong differences in conductivity and relative permittivity between the tissue classes are expected to affect both the propagation and absorption of the electromagnetic field.

3. FDTD Formulation and Numerical Parameters

3.1. Maxwell Equations

The electromagnetic propagation is simulated using the FDTD method introduced by Yee [3]. The two-dimensional formulation adopted here corresponds to the transverse magnetic (TM_z) configuration, for which the electric field has one non-zero component, $E_z$ , while the magnetic field has components $H_x$ and $H_y$ .

The governing equations are written as

$\frac{\partial D_z}{\partial t}=\frac{1}{\sqrt{{\mu }_0{\epsilon }_0}}\left(\frac{\partial H_y}{\partial x}-\frac{\partial H_x}{\partial y}\right),$ (1)

$\frac{\partial H_x}{\partial t}=-\frac{1}{\sqrt{{\mu }_0{\epsilon }_0}}\frac{\partial E_z}{\partial y},$ (2)

and

$\frac{\partial H_y}{\partial t}=-\frac{1}{\sqrt{{\mu }_0{\epsilon }_0}}\frac{\partial E_z}{\partial x}.$ (3)

For a lossy dielectric, the frequency-domain constitutive relation can be expressed in terms of the complex relative permittivity

${\epsilon }_r^{\text{*}}\left(\omega \right)={\epsilon }_r+\frac{\sigma }{j\omega {\epsilon }_0},$ (4)

so that

$D_z\left(\omega \right)={\epsilon }_r^{\text{*}}\left(\omega \right) E_z\left(\omega \right).$ (5)

3.2. Propagation Direction and Polarization

The two-dimensional computational domain corresponds to the selected coronal anatomical slice and is represented in the $\left(x,y\right)$ plane. The electromagnetic excitation is restricted to the TM_z configuration. Accordingly, the electric field has only an out-of-plane component,

$E=E_z\left(x,y,t\right) \widehat{z},$

whereas the magnetic field lies in the plane of the coronal slice,

$H=H_x\left(x,y,t\right)\widehat{x}+H_y\left(x,y,t\right)\widehat{y}.$

The incident plane wave propagates along the $+x$ direction. Consequently, the present configuration represents a single, highly specific two-dimensional exposure geometry with a fixed propagation vector and a uniform out-of-plane polarization. The dosimetric results and conclusions reported in this work are strictly bounded by this 2D framework and TM_z polarization state, and must not be interpreted as polarization-independent characteristics of electromagnetic absorption in a real three-dimensional human head.

3.3. UPML Absorbing Boundary

A uniaxial perfectly matched layer (UPML) is used to emulate an open electromagnetic domain and reduce artificial reflections at the boundaries [14] [15]. The absorbing boundary condition is implemented using a standard polynomial grading profile of order $m=3$ for the conductivity scaling, with a target theoretical reflection coefficient at normal incidence set to $R_0={10}^{-6}$ (−120 dB).

The implementation uses the standard auxiliary-variable formulation. The electric-displacement update is written in the form

$\begin{array}{c}{D}_z^{n+1/2}\left(i,j\right)=g_i^{\left(3\right)}\left(i\right)g_j^{\left(3\right)}\left(j\right)D_z^{n-1/2}\left(i,j\right)\\ +g_i^{\left(2\right)}\left(i\right)g_j^{\left(2\right)}\left(j\right)\frac{\Delta t}{\sqrt{{\mu }_0{\epsilon }_0}\Delta x}[H_y^n\left(i+\frac{1}{2},j\right)-H_y^n\left(i-\frac{1}{2},j\right)\\ -H_x^n\left(i,j+\frac{1}{2}\right)+H_x^n\left(i,j-\frac{1}{2}\right)].\end{array}$ (6)

The magnetic-field updates are

$\begin{array}{c}{H}_x^{n+1}\left(i,j+\frac{1}{2}\right)=f_j^{\left(3\right)}\left(j+\frac{1}{2}\right)H_x^n\left(i,j+\frac{1}{2}\right)\\ +f_j^{\left(2\right)}\left(j+\frac{1}{2}\right)\frac{\Delta t}{\sqrt{{\mu }_0{\epsilon }_0}\Delta x}\left[E_z^{n+1/2}\left(i,j\right)-E_z^{n+1/2}\left(i,j+1\right)\right]\\ +f_i^{\left(1\right)}\left(i\right)I_{H_x}^{n+1/2}\left(i,j+\frac{1}{2}\right),\end{array}$ (7)

and

$\begin{array}{c}{H}_y^{n+1}\left(i+\frac{1}{2},j\right)=f_i^{\left(3\right)}\left(i+\frac{1}{2}\right)H_y^n\left(i+\frac{1}{2},j\right)\\ +f_i^{\left(2\right)}\left(i+\frac{1}{2}\right)\frac{\text{Δ}t}{\sqrt{{\mu }_0{\epsilon }_0}\text{Δ}x}\left[E_z^{n+1/2}\left(i+1,j\right)-E_z^{n+1/2}\left(i,j\right)\right]\\ +f_j^{\left(1\right)}\left(j\right)I_{H_y}^{n+1/2}\left(i+\frac{1}{2},j\right).\end{array}$ (8)

The auxiliary variables and dissipative coefficients are determined by the selected UPML profile.

3.4. Electric-Field Update in Lossy Tissues

After updating the electric displacement, the electric field is recovered according to

$E_z^{n+1/2}\left(i,j\right)=\frac{D_z^{n+1/2}\left(i,j\right)-\frac{\sigma \left(i,j\right)\Delta t}{{\epsilon }_0}{E}_z^{n-1/2}\left(i,j\right)}{{\epsilon }_r\left(i,j\right)+\frac{\sigma \left(i,j\right)\Delta t}{{\epsilon }_0}}.$ (9)

This formulation allows each FDTD cell to be associated with the electromagnetic properties of the corresponding segmented tissue.

3.5. Incident-Field Excitation

A normally incident plane wave is introduced using a total-field/scattered-field formulation, following the two-dimensional prescription of Sullivan [15]. The TF/SF boundary is positioned exactly 10 cells away from the inner interface of the UPML absorbing layer. The incident electric field is defined by

$E_{\text{inc}}\left(x,t\right)=E_0\sin \left(\omega t-\tilde{k}x\right)\left(1-{\text{e}}^{-t/\tau }\right),$(10)

where

$E_0=27.45\text{V}/\text{m},\omega =2\pi f_0.$ (11)

The amplitude $E_0$ is chosen such that the incident time-averaged power density $S_{\text{inc}}=E_0^2/\left(2Z_0\right)$ equals 1 W/m², with $Z_0\simeq 377\Omega$ the free-space impedance. Under this normalization, the SAR values reported in Section 5 scale linearly with the actual incident power density and can therefore be transposed to any exposure scenario.

The factor $1-{\text{e}}^{-t/\tau }$ provides a progressive establishment of the field and reduces the influence of the initial transient. The numerical wavenumber is calculated from the FDTD dispersion relation:

$\tilde{k}=\frac{2}{\Delta x}\arcsin \left(\frac{\Delta x}{c\Delta t}\sin \left(\frac{\omega \Delta t}{2}\right)\right).$(12)

3.6. Numerical Parameters

Table 2. Main numerical parameters of the FDTD simulation.

Parameter	Value
Frequency	900 MHz
Spatial step	$\Delta x=\Delta y=1\text{mm}$
Temporal constant	$\tau =3T_0$
Time step	$\Delta t=\Delta x/\left(2c\right)$
Computational domain	280 × 290 cells
UPML thickness	40 cells
UPML Profile order ($m$ )	3
UPML Target reflection ($R_0$ )	10−6
TF/SF boundary offset	10 cells from UPML
Minimum Air padding	15 cells
Number of time steps	15,666
Incident electric-field amplitude	27.45 V/m

The principal numerical parameters are summarized in Table 2. At the operating frequency of 900 MHz, the optical period is $T_0=1/f_0\simeq 1.111$ ns. With the chosen time step $\Delta t=\Delta x/\left(2c\right)\simeq 1.6678$ ps, a single harmonic period is discretized into exactly 666 time steps. The total number of time steps (15,666) corresponds to 47.04 full optical periods. This window consists of a transient phase of 15,000 steps (equivalent to 22.5 periods), after which the ramp envelope $1-{\text{e}}^{-t/\left(3T_0\right)}$ exceeds 0.9994 and the SAR bias associated with an unestablished regime remains below 0.1%, followed by a single full period (666 steps) used to compute the root-mean-square electric field. Given the lossy nature of the biological tissues, internal multi-reflections are rapidly attenuated, ensuring that the computational domain has rigorously reached a harmonic steady state prior to the final sampling period.

4. SAR Calculation

The local SAR in tissue $k$ is evaluated from the RMS electric field as

${\text{SAR}}_k=\frac{{\sigma }_k{\left|E_{k,\text{rms}}\right|}^2}{{\rho }_k}.$ (13)

Here, ${\sigma }_k$ is the electrical conductivity and ${\rho }_k$ is the mass density of tissue $k$ .

The RMS electric field is calculated after the transient regime has decayed and the harmonic steady state has been reached. The SAR is evaluated cell by cell and then statistically analyzed separately for each tissue class.

For each tissue, the following quantities are calculated:

• mean SAR,

• standard deviation,

• maximum local SAR.

The resulting distribution provides a direct measure of the spatial heterogeneity of electromagnetic energy absorption.

5. Results and Discussion

5.1. Spatial Distribution of Dielectric Properties

The segmented anatomical model produces a strongly heterogeneous distribution of electromagnetic properties as shown in Figure 2. In particular, the CSF presents a much larger conductivity than gray matter and white matter, whereas the skull exhibits comparatively low conductivity.

The heterogeneous distribution of permittivity and conductivity determines the local propagation, reflection, transmission, and absorption of the electromagnetic field.

5.2. SAR Distribution

The spatial distribution of SAR at 900 MHz is shown in Figure 3.

The calculated SAR distribution is strongly non-uniform. Local variations occur near tissue interfaces and in regions where the electric field is enhanced by electromagnetic interference.

This behavior results from the differences in permittivity and conductivity between the tissue classes. These differences modify the local electric field and consequently the absorbed electromagnetic power.

Figure 2. Spatial distribution of the dielectric parameters in the segmented brain model.

Figure 3. Spatial distribution of SAR in the heterogeneous brain model at 900 MHz.

5.3. Statistical Analysis of SAR

Table 3 summarizes the mean, standard deviation, and maximum SAR obtained for the four tissue classes.

The average SAR in gray matter is higher than that in white matter. This is consistent with the higher conductivity of gray matter.

Table 3. SAR statistics for the different tissue types at 900 MHz.

Tissue	Mean SAR (10−3 W/kg)	Standard deviation (10−3 W/kg)	Maximum (10−3 W/kg)
White matter	3.08	3.69	25.27
Gray matter	5.73	7.40	41.36
CSF	12.15	16.20	109.64
Skull	1.82	2.54	12.92

The CSF exhibits the largest average and maximum SAR values among the considered tissues. Its comparatively high electrical conductivity facilitates conversion of electromagnetic-field energy into dissipated power.

The dependence of the SAR on the local tissue properties can be made explicit by rewriting Equation (14) as ${\text{SAR}}_k=\left({\sigma }_k/{\rho }_k\right)\langle {\left|E_k\right|}^2\rangle$ . This decomposition separates a material factor ${\sigma }_k/{\rho }_k$ , fixed by the tissue class, from a field factor $\langle {\left|E_k\right|}^2\rangle$ , which depends on the local electromagnetic environment. The ratios obtained from Table 3 illustrate the balance between these two contributions. The CSF-to-gray-matter SAR ratio is 2.12, whereas the corresponding $\sigma /\rho$ ratio is 2.66; the mean square field is therefore smaller in the CSF by a factor of about 0.80, consistent with the partial shielding associated with the high permittivity of this tissue. Conversely, the skull-to-white-matter SAR ratio is 0.59 for a $\sigma /\rho$ ratio of only 0.13, indicating a mean square field about 4.6 times larger in the skull. This is consistent with the peripheral location of the skull, where the field is enhanced by the incident wave before it is progressively absorbed by the underlying tissues. Conductivity alone therefore does not determine the spatial distribution of the SAR: the local field structure plays a comparable role.

The standard deviations are relatively large compared with the corresponding mean values, indicating substantial spatial variability within individual tissue classes. This variability is associated with the non-uniform electric field resulting from the anatomical interfaces and multiple electromagnetic interactions.

The skull exhibits the lowest average SAR among the considered tissues, consistent with its relatively low conductivity.

The absolute SAR values reported here correspond to a normalized incident power density $S_{\text{inc}}=1$ W/m². The largest local value (about 0.11 W/kg in the CSF) is therefore approximately eighteen times smaller than the ICNIRP local basic restriction of 2 W/kg for the head and trunk averaged over 10 g of tissue [1], and correspondingly below the IEEE C95.1-2019 basic restriction [2]. The comparison is illustrative rather than regulatory, since a full compliance assessment would require a three-dimensional anatomical model, a 10 g spatial averaging procedure, and a realistic exposure configuration.

5.4. SAR Attenuation along the Propagation Direction

The mean SAR profile along the direction of propagation is shown in Figure 4.

Figure 4. SAR attenuation along the direction of electromagnetic-wave propagation.

The SAR generally decreases from the incident side toward the opposite side of the brain, reflecting the progressive absorption of electromagnetic energy.

The decrease is not strictly monotonic because the wave interacts with successive tissue interfaces. Local increases in SAR may occur in regions associated with tissues having higher conductivity, particularly the CSF.

These observations illustrate the importance of representing the heterogeneous anatomical structure of the brain when evaluating local electromagnetic absorption.

6. Discussion

The numerical results demonstrate that the SAR distribution in a heterogeneous brain model cannot be described adequately by a uniform absorption model. The electromagnetic response depends simultaneously on the local tissue properties and their anatomical arrangement.

The CSF exhibits the largest mean and maximum SAR values in the present simulation. This behavior is consistent with its comparatively high electrical conductivity. However, conductivity alone does not completely determine the spatial distribution of SAR, because the local electric field is also controlled by reflection, transmission, interference, and geometrical effects.

The decomposition ${\text{SAR}}_k=\left({\sigma }_k/{\rho }_k\right)\langle {\left|E_k\right|}^2\rangle$ provides a compact way to quantify this statement. In the present model, the material factor ${\sigma }_k/{\rho }_k$ spans nearly two orders of magnitude between the skull and the CSF, whereas the tissue-averaged squared field spans roughly a factor of five in the opposite direction. The final SAR ordering is set by the competition between these two factors: the CSF combines a large $\sigma /\rho$ with a moderately reduced field to yield the largest absorption, while the skull combines the smallest $\sigma /\rho$ with the largest peripheral field to yield an intermediate absorption level. This structure supports the general observation that dielectric contrasts drive both the local field redistribution and the resulting energy deposition, reinforcing the relevance of anatomically heterogeneous head models compared with homogeneous or few-layered approximations.

The quantitative verification against the analytical Mie solution for a lossy dielectric cylinder, detailed in Appendix A.7, yields a tissue-averaged SAR error of 2.7% at the working resolution of $\Delta x=1$ mm. However, in accordance with the well-documented behavior of staircased FDTD meshes on curved lossy boundaries, this error does not scale monotonically with sub-millimeter refinement due to localized phase dispersion mismatches. Therefore, this 2.7% metric must be interpreted strictly as a performance indicator for the canonical cylinder geometry, rather than a universal accuracy floor or direct uncertainty bound for the heterogeneous brain anatomy. Nevertheless, the macroscopically resolved tissue boundaries at 1 mm, coupled with the massive physical contrasts in tissue conductivities, ensure that the reported statistical SAR ordering (CSF > Gray Matter > White Matter > Skull) remains robust and invariant to grid-sensitivity variations.

The substantial standard deviations observed within each tissue class further demonstrate that a single average value does not fully describe the local dosimetric behavior.

The present results are qualitatively consistent with previous FDTD investigations of heterogeneous head models [7] [10] [11]. They also reinforce the importance of realistic anatomical representations in numerical dosimetry.

The two-dimensional nature of the present model must nevertheless be emphasized. The simulation represents one selected anatomical section and does not reproduce the complete three-dimensional geometry of the human brain. Consequently, out-of-plane scattering, axial wave diffraction, cross-polarization variables, and complex near-field effects typical of real-world exposure scenarios—such as mobile phone proximity or ambient oblique fields—are entirely ignored. The present numerical results are strict mathematical descriptions of a transverse-magnetic 2D exposure and must not be used as polarization-independent compliance limits for a real 3D human head. Instead, this 2D formulation serves as a highly controlled, verified environment for examining the intrinsic coupling between dielectric tissue property factors and field distribution profiles.

In addition, the present study considers one frequency and one incidence configuration. A systematic investigation of frequency, polarization, incidence angle, and anatomical orientation is therefore required before general conclusions concerning RF exposure can be drawn.

Finally, the present SAR values describe electromagnetic energy absorption and do not directly provide the associated temperature rise. Such an analysis requires coupling the electromagnetic model to a thermal model, for example through the Pennes bioheat equation.

7. Conclusions

A two-dimensional FDTD model was developed to investigate the spatial distribution of the specific absorption rate in a heterogeneous human brain model derived from an MRI-based anatomical dataset.

The electromagnetic problem was solved at 900 MHz using a TF/SF formulation for the normally incident plane wave and a UPML to truncate the computational domain. The anatomical model includes gray matter, white matter, cerebrospinal fluid, and skull, with electromagnetic properties assigned from published tissue-property data.

The calculated SAR distribution is strongly heterogeneous. The CSF exhibits the highest mean and maximum SAR values, whereas the skull presents the lowest mean absorption among the considered tissues. The large spatial variability within each tissue class demonstrates the importance of local electromagnetic-field distributions in determining dosimetric quantities. The tissue ordering is quantitatively interpreted through the decomposition of the SAR into a material factor and a field factor, showing that both factors contribute at comparable levels to the observed absorption pattern.

A systematic verification strategy was also established for the numerical implementation. It includes one-dimensional benchmarks for numerical dispersion, dielectric interfaces, and parallel dielectric slabs, together with a two-dimensional comparison against the analytical Mie solution for a lossy dielectric cylinder. These tests provide quantitative verification of the numerical solver for the corresponding canonical configurations. The Mie-cylinder benchmark, in particular, characterizes the performance of the complete solver at the working spatial resolution, but does not by itself constitute a direct uncertainty bound for the SAR values obtained in the anatomically heterogeneous model.

The primary limitation of this work remains its simplified two-dimensional planar configuration, which omits out-of-plane scattering and cannot reproduce the complex field patterns induced by 3D head features or portable communication devices. The quantitative framework established here constitutes an essential validation benchmark rather than a complete standard compliance test. Future works will extend this verified foundation toward complete three-dimensional anatomical head models, heterogeneous mobile phone source placements, wide-spectrum parametric frequency sweeps, and fully coupled electromagnetic-thermal solvers using the Pennes bioheat equation to assess local temperature elevations.

Appendix A: Verification

Before applying the FDTD model to the heterogeneous anatomical brain, the numerical implementation was tested using canonical configurations for which analytical solutions are available. The verification strategy combines one-dimensional benchmarks that characterize the intrinsic properties of the Yee scheme (numerical dispersion, dielectric interface, parallel dielectric slab, mesh convergence in a homogeneous medium) with a two-dimensional benchmark against the analytical Mie solution for a lossy dielectric cylinder. The latter test exercises the complete solver chain - UPML, TF/SF, staircased curved interfaces, lossy medium - on a configuration whose dielectric properties are representative of the anatomical tissues used in the anatomical simulation. Table A1 summarizes the principal quantitative results.

Table A1. Quantitative results of the numerical verification tests.

Test	Quantity	Result
Numerical dispersion (1D)	Maximum relative error, $\text{Δ}x\in \left[0.5,2\right]$ mm	$4.44\times {10}^{-5}$
Numerical dispersion (2D)	Anisotropy at $\text{Δ}x=1$ mm	$<{10}^{-5}$
UPML reflection (2D)	$\left|R\right|$ for 40 cells	$4.5\times {10}^{-6}$ (−107 dB)
Dielectric interface	Relative error in $R$	$4.03\times {10}^{-4}$
Dielectric interface	Relative error in $T$	$9.48\times {10}^{-5}$
Parallel slab	Relative error in $R$ (frequency-domain extraction)	$3.9\times {10}^{-4}$
Parallel slab	Relative error in $T$ (frequency-domain extraction)	$1.1\times {10}^{-4}$
Mesh convergence (homogeneous)	Observed order	1.999
Mie cylinder (2D)	$L_2$ error on the field	$3.1\times {10}^{-2}$
Mie cylinder (2D)	Error on the tissue-averaged SAR	$2.7\times {10}^{-2}$

The detailed configurations, analytical reference solutions, numerical results, error metrics, and additional verification figures are provided in the following subsections.

A.1. Free-Space Propagation

The first test considers the propagation of a modulated plane-wave packet in a homogeneous lossless medium. The analytical solution is

$E_z\left(x,t\right)=E_0\exp \left[-{\left(\frac{x-x_0-ct}{w}\right)}^{\text{}2}\right]\cos \left[k\left(x-x_0-ct\right)\right],k=\frac{\omega }{c}$ (14)

The numerical and analytical spatial profiles are compared after a prescribed propagation time. The relative $L^2$ error is defined by

${ϵ}_{\text{FS}}=\frac{{\Vert E_{\text{FDTD}}-E_{\text{analytical}}\Vert }_2}{{\Vert E_{\text{analytical}}\Vert }_2}$ (15)

For $\text{Δ}x=1$ mm, the measured relative error is 8.10 × 10⁻³. The numerical and analytical profiles are shown in Figure A1.

Figure A1. Free-space propagation test: FDTD field compared with the analytical solution.

The error is sufficiently small for the intended anatomical simulation, but it also shows that the finite grid and time step produce a measurable numerical phase and amplitude error. This error is further quantified independently through the dispersion test.

A.2. Numerical Dispersion

For the one-dimensional Yee scheme, the numerical wavenumber satisfies

$\sin \left(\frac{\omega \Delta t}{2}\right)=\frac{c\Delta t}{\Delta x}\sin \left(\frac{\tilde{k}\Delta x}{2}\right).$(16)

Hence,

$\tilde{k}=\frac{2}{\Delta x}\arcsin \left(\frac{\Delta x}{c\Delta t}\sin \left(\frac{\omega \Delta t}{2}\right)\right).$(17)

The relative dispersion error is

${ϵ}_k=\frac{\left|\tilde{k}-k\right|}{k}.$ (18)

At 900 MHz, with the Courant factor $S=c\Delta t/\Delta x=0.5$ used in the anatomical simulation, the maximum relative error over $\Delta x\in \left[0.5,2\right]$ mm is 4.44 × 10⁻⁵, and the error at $\Delta x=1$ mm is 1.11 × 10⁻⁵. Figure A2 shows the rapid reduction of the dispersion error with mesh refinement.

Figure A2. Numerical-dispersion error as a function of the spatial step at 900 MHz.

A.3. UPML Reflection (Two-Dimensional Benchmark)

Figure A3. Residual reflection of the 2D UPML at 900 MHz.

The residual reflection of the two-dimensional UPML was characterized directly on the production solver by propagating a plane wave in a homogeneous medium (free space) using the TF/SF formulation, and measuring the residual field in the scattered-field region after the transient regime has decayed. In this configuration, the field in the scattered-field region should vanish identically, so any non-zero value measures the combined residual of UPML reflection and TF/SF leakage. The measured reflection coefficient is

$\left|R\right|=\frac{\underset{\text{SF}}{\text{max}}\left|E_z\right|}{E_0}.$ (19)

For UPML thicknesses of 10, 20, 30 and 40 cells, the measured values are 2.6 × 10⁻⁵, 6.4 × 10⁻⁶, 5.5 × 10⁻⁶ and 4.5 × 10⁻⁶, respectively, corresponding to reflection levels of −92, −104, −105 and −107 dB. The residual reflection saturates at approximately −107 dB for thicknesses of 20 cells and above. Figure A3 shows the resulting reflection coefficient as a function of the absorbing-layer thickness.

A.4. Verification at a Dielectric Interface

The material-interface treatment was verified for normal incidence on a lossless interface between two media with ${\epsilon }_{r1}=1$ and ${\epsilon }_{r2}=4$ . The analytical reflection and transmission coefficients are

$R={\left|\frac{n_1-n_2}{n_1+n_2}\right|}^2,T=\frac{n_2}{n_1}{\left|\frac{2n_1}{n_1+n_2}\right|}^2.$ (20)

The analytical values are $R=0.111111$ and $T=0.888889$ . The FDTD calculation gives $R=0.111156$ and $T=0.888973$ , corresponding to relative errors of 4.03 × 10⁻⁴ and 9.48 × 10⁻⁵, respectively. The results illustrated in Figure A4 demonstrate that the material-interface treatment is accurately reproduced by the FDTD discretization.

Figure A4. FDTD reflected and transmitted fields for a normally incident wave at a dielectric interface.

A.5. Verification Using a Parallel Dielectric Slab

A finite dielectric slab with ${\epsilon }_r=4$ and thickness $d=20$ mm was considered as a second interface benchmark. For a lossless slab, the analytical reflection and transmission coefficients were calculated using the Fabry-Pérot formulation. Energy conservation imposes

$A=1-R-T,$ (21)

where $A$ is the absorption coefficient, which vanishes identically for a lossless medium. The analytical values at 900 MHz are $R=0.208788$ , $T=0.791212$ , and $A=0$ . A frequency-domain extraction of the FDTD reflected and transmitted fields at 900 MHz gives $R=0.208685$ , $T=0.791308$ , and $A\simeq 7\times {10}^{-6}$ , the latter being a numerical round-off-level residual of the energy balance. The relative errors in R and T are 4.9 × 10^-4 and 1.2 × 10⁻⁴ respectively, consistent with those obtained for the single interface (Section A.4). The small energy-balance residual, several orders of magnitude smaller than the individual errors in R and T, indicates that the numerical scheme conserves electromagnetic energy: the errors on R and T compensate one another rather than accumulating.

Figure A5 compares the analytical and FDTD reflection, transmission, and absorption coefficients for the dielectric slab benchmark.

Figure A5. Comparison of analytical and FDTD reflection, transmission and absorption coefficients for the dielectric slab benchmark.

A.6. Mesh-Convergence Study in a Homogeneous Medium

The spatial convergence of the Yee scheme was evaluated using a periodic plane-wave problem in a homogeneous, lossless medium. The spatial step was varied from 2 to 0.5 mm, while the Courant factor was kept constant. The relative $L_2$ error was calculated with respect to the analytical plane-wave solution.

The measured errors are 1.26 × 10⁻³, 7.09 × 10⁻⁴, 3.15 × 10⁻⁴, 1.77 × 10⁻⁴ and 7.86 × 10⁻⁵ for $\text{Δ}x=2$ , 1.5, 1, 0.75 and 0.5 mm, respectively. A least-squares fit in logarithmic coordinates gives an observed convergence order of 1.999, in excellent agreement with the expected second-order convergence of the Yee scheme.

This result certifies the correctness of the implementation in the regime where the scheme admits an exact dispersion analysis. The behavior of the solver in the presence of curved staircased interfaces, which is the regime relevant to the anatomical model, is characterized independently in Section A.7. The measured convergence is shown in Figure A6.

Figure A6. Spatial convergence of the FDTD solution. The observed order is approximately two.

A.7. Two-Dimensional Benchmark against the Mie Series

The complete solver - Yee scheme, staircased curved interfaces, UPML, TF/SF - was benchmarked against the analytical Mie series for a two-dimensional lossy dielectric cylinder illuminated by a plane wave. The cylinder radius is 40 mm, and its dielectric properties are those of gray matter at 900 MHz (${\epsilon }_r=52.7$ , $\sigma =0.94$ S/m, $\rho =1045$ kg/m³). The complex Mie series was truncated at $N_{\text{max}}=80$ multipoles, well above the convergence radius $ka\simeq 0.75$ that governs the truncation error.

At the working spatial step $\text{Δ}x=1$ mm, the relative $L_2$ error between the FDTD phasor and the Mie phasor is 3.1% inside the cylinder and 2.5% on a 5-mm-thick boundary annulus. The corresponding error on the tissue-averaged SAR is 2.7% on the mean value and 4.3% on the maximum value. When the spatial step is refined to 0.75 and 0.5 mm, the $L_2$ error on the field reaches 3.2% and 4.0%, respectively. The error does not decrease with refinement, which is consistent with the theoretical behavior of a Yee-based TF/SF formulation on a lossy scatterer: the residual is dominated by a small dispersion mismatch between the one-dimensional numerical wavenumber used to inject the incident field and the two-dimensional numerical dispersion of the surrounding grid. Consequently, these benchmark errors characterize the baseline performance of the mathematical formulation on a staircased canonical object, but they do not constitute a direct baseline or geometric limit for the real anatomical brain configuration. The observed error level remains fully comparable to that reported in the FDTD dosimetric literature for staircased voxel models [11] [16].

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	International Commission on Non-Ionizing Radiation Protection (2020) Guidelines for Limiting Exposure to ELECTROMAGNETIC Fields (100 kHz to 300 GHz). Health Physics, 118, 483-524. [Google Scholar] [CrossRef]
[2]	IEEE (2019) IEEE Standard for Safety Levels with Respect to Human Exposure to Electric, Magnetic, and Electromagnetic Fields, 0 Hz to 300 GHz. IEEE Std C95.1-2019, 1-312.[CrossRef]
[3]	Yee, K. (1966) Numerical Solution of Initial Boundary Value Problems Involving Maxwell’s Equations in Isotropic Media. IEEE Transactions on Antennas and Propagation, 14, 302-307.[CrossRef]
[4]	Taflove, A., Hagness, S.C. and Piket-May, M. (2005) Computational Electromagnetics: The Finite-Difference Time-Domain Method. In: Chen, W.K., Ed., The Electrical Engineering Handbook, Elsevier, 629-670.[CrossRef]
[5]	Lwin, Z.M. and Yokota, M. (2019) Numerical Analysis of SAR and Temperature Distribution in Two Dimensional Human Head Model Based on FDTD Parameters and the Polarization of Electromagnetic Wave. AEU—International Journal of Electronics and Communications, 104, 91-98.[CrossRef]
[6]	Sabbah, A.I., Dib, N.I. and Al-Nimr, M.A. (2011) Evaluation of Specific Absorption Rate and Temperature Elevation in a Multi-Layered Human Head Model Exposed to Radio Frequency Radiation Using the Finite-Difference Time Domain Method. IET Microwaves, Antennas & Propagation, 5, 1073-1080.[CrossRef]
[7]	Gandhi, O.P., Lazzi, G. and Furse, C.M. (1996) Electromagnetic Absorption in the Human Head and Neck for Mobile Telephones at 835 and 1900 MHz. IEEE Transactions on Microwave Theory and Techniques, 44, 1884-1897.[CrossRef]
[8]	Bielli, P., Richiardi, G., Bertotto, P. and Schiavoni, A. (2000) SAR Generated by Commercial Cellular Phones-Phone Modeling, Head Modeling, and Measurements. IEEE Transactions on Microwave Theory and Techniques, 48, 2064-2071.[CrossRef]
[9]	Husni, N.A., Islam, M.T., Faruque, M.R.I. and Misran, N. (2013) Effects of Electromagnetic Absorption towards Human Head Due to Variation of Its Dielectric Properties at 900, 1800 and 1900 MHz with Different Antenna Substrates. Progress In Electromagnetics Research, 138, 367-388.[CrossRef]
[10]	Wiart, J., Hadjem, A., Wong, M.F. and Bloch, I. (2008) Analysis of RF Exposure in the Head Tissues of Children and Adults. Physics in Medicine and Biology, 53, 3681-3695.[CrossRef] [PubMed]
[11]	Dimbylow, P.J. and Gandhi, O.P. (1991) Finite-Difference Time-Domain Calculations of SAR in a Realistic Heterogeneous Model of the Head for Plane-Wave Exposure from 600 MHz to 3 GHz. Physics in Medicine and Biology, 36, 1075-1089.[CrossRef] [PubMed]
[12]	Cocosco, C.A., Kollokian, V., Kwan, R.K.S. and Evans, A.C. (1997) BrainWeb: Online Interface to a 3D MRI Simulated Brain Database. NeuroImage, 5, S425.
[13]	Gabriel, S., Lau, R.W. and Gabriel, C. (1996) The Dielectric Properties of Biological Tissues: III. Parametric Models for the Dielectric Spectrum of Tissues. Physics in Medicine and Biology, 41, 2271-2293.[CrossRef] [PubMed]
[14]	Gedney, S.D. (1996) An Anisotropic Perfectly Matched Layer-Absorbing Medium for the Truncation of FDTD Lattices. IEEE Transactions on Antennas and Propagation, 44, 1630-1639.[CrossRef]
[15]	Sullivan, D.M. andHoule, J.E. (2020) Electromagnetic Simulation Using the FDTD Method with Python. 3rd Edition, Wiley.
[16]	Wang, J., Fujiwara, O., Kodera, S. and Watanabe, S. (2006) FDTD Calculation of Whole-Body Average SAR in Adult and Child Models for Frequencies from 30 MHz to 3 GHz. Physics in Medicine and Biology, 51, 4119-4127.[CrossRef] [PubMed]