*Q*helical antennae and are dominant in shaping the spectral response in the sub-terahertz region. We found that water absorption, dispersion and multiple interference effects play the major role in shaping the spectrum without the need for the assumption of the sweat ducts acting as low-

*Q*helical antennae. High sensitivities to the water content are found particularly in the ellipsometric parameters at large incidence angles. Hence a new methodology is proposed to detect skin cancer using variable angle ellipsometry or polarized reflectometry. The parameter found with the highest sensitivity to water content is cos Δpp with Δpp being the phase of the on-diagonal reflection matrix ratio between

*p*-to-

*p*polarization.

## 1.

## Introduction

Skin cancer, particularly malignant melanoma (MM) with over 60000 new cases and over 7000 deaths in the US estimated^{1} in 2007, can be cured if detected at early stages. To the total of 59940 new cases of melanoma, one can add an estimated 50000 cases of melanoma *in situ*, a number growing at 15% per/yr in some countries. Unfortunately, the accuracy of physicians in diagnosing MM in the clinic is not high particularly at its early stages because early MM, those having a Breslow thickness less than 1 mm, share many diagnostic features with benign lesions, such as dysplastic nevi, atypical melanocytic nevi or even pigmented basal cell carcinoma (BCC) and keratoacantoma. In a recent study,^{2} general practitioners had a sensitivity and specificity for detection of melanoma of 62%, while dermatologists had a corresponding sensitivity and specificity of 80% and 60%. Although dermatoscopic examination with its better illumination and magnification can improve the accuracy of melanoma diagnostics,^{3} there is still a need for a robust and reliable technique to specify and diagnose melanoma and other skin cancers at an early stage. Dermatoscopy has improved the diagnosis accuracy recently by the incorporation of digital image processing,^{4} and multispectral polarized operation.^{5} However, most of the works published so far use two polarization states, parallel and crossed, and the number of wavelengths is limited and does not cover both the visible and the near-infrared range.^{6} Using multispectral imaging alone has been shown recently^{7} to give limited information. The same is true when using polarimetric imaging alone, while the combination certainly enhances the information on the pictured lesion.^{8, 9} Hence, there is still a need for novel approaches for skin cancer imaging.

Terahertz radiation (or T-rays)—with frequencies between 100 GHz and 10 THz – does not suffer from the scattering that usually circumvents visible and infrared radiation to penetrate the skin tissue. Molecular rotations and dipolar resonances do exist in this electromagnetic range and therefore signatures of water and other biomaterials can be detected using terahertz radiation. For example, T-rays have been demonstrated to effectively image skin burn severity, tooth cavities, and skin cancer. The high dielectric constant of water at these frequencies means that terahertz imaging can detect slight variations in the water content of biological materials.^{10, 11} Together with the nonionizing nature of this radiation and its ability to penetrate clothing with minimal scattering, terahertz imaging appears particularly suited to examine skin abnormalities such as burns, scars, wounds, and skin cancers. Due to the long wavelength compared to the local inhomogenieties in the skin tissue, terahertz waves propagation properties in the skin maybe modelled using some homogenization of the skin layers.^{12} However some anisotropic entities do exist in the skin such as the collagen fibers, the sweat glands, and their helical ducts. The effect of the helical sweat glands on the terahertz spectrum at normal incidence was shown by us recently^{13} not to be significant contrary to some other results^{14, 15} which claimed that they act as low *Q* helical antennae. In this article an anisotropic homogenization procedure is applied and the reflectometric and ellipsometric spectra from the skin versus angle and versus frequency are calculated.

Scattering and absorption of biological tissue depend strongly on the radiation frequency and the polarization. While a more accurate model including scattering mechanisms will be addressed in future research, we discuss here frequencies at which only low scattering occurs. Since dealing with wavelengths in a much larger scale than the region where scattering mechanisms take place—long compared to non homogeneity entities (sweat ducts) and the intralayer skin interface corrugation (tens of micrometers—according to standard histology skin slides), and short compared to the external skin corrugation for relatively young healthy specimens. This simplistic approach for skin modeling has been proven to provide relatively exact results at the specific wavelength region when compared to experimental measurements.^{14, 15} In the low scattering regime, as is discussed in the following, polarization dependence can result from structural properties of the tissue such as the layered nature of the skin, a well known phenomenon of interference in multilayer structures particularly for oblique incidence. In the terahertz range the water in the biological tissue is the main absorber. Because cancerous tissue is overhydrated as compared to the normal tissue, a high contrast reflectance image is obtained in particular when certain terahertz frequencies are selected at which the water absorption is high. Similarly contrast ellipsometric images can be obtained due to polarization dependent localized interference changes due to water content induced dielectric alterations in the skin layers. In addition, reflection from normal tissue is polarization dependent due to the existence of ordered anisotropic entities in the skin such as collagen fibers and sweat glands, while this order is damaged in the cancerous tissue. The scattering, although expected to be weak since the terahertz wavelength is much larger than the inhomogenieties within the skin layers, can give information on tumors in their early stages of sizes of the order of 0.1 mm. Although we are not dealing with scattering here, in Fig. 1 we have drawn a scheme (the vertical imaging column) that will allow measuring the scattering in addition to the specular reflection. Hence the importance of grabbing images at different polarizations or grabbing ellipsometric images. Terahertz ellipsometers capable of measuring spectroscopic Mueller matrix data at terahertz frequencies have already been presented^{16} and applied for several nonmedical uses,^{17, 18} as well as time domain terahertz ellipsometry.^{19} An example of a setup configuration for reflectometric ellipsometric scatterometric type skin imaging system is shown in Fig. 1. The terahertz beam is polarized and incident at an oblique angle and the specularly reflected beam is collected. In Fig. 1 two imaging paths are shown, one for the scattered waves and one for the specularly reflected waves. As mentioned before, although the scattering is not dealt with here and considered very weak, the vertical column in Fig. 1 is drawn for possible scatterometric measurement. A 2nd polarizer (analyzer) in the imaging paths allows measuring the Stokes images, the Mueller matrix components, or the ellipsometric images which we will discuss later.

In this article we are interested mainly in the specular reflection images either by measuring the intensity or the ellipsometric parameters. The specular reflection images are expected to be strong and should be easily modeled as multilayered thin films. Our starting point is a characterization of the skin tissue in an effort to create a reliable model for investigation, discussing the electrical behavior of biological matter in Sec. 2.1 and structural considerations, skin layered nature in Sec. 2.2.1, and the electrical effect of sweat ducts in Sec. 2.2.2 A homogenization-based approach for modeling the sweat ducts in Sec. 2.2.3, accompanied with a short description of the basics of homogenization theory, is presented in the Appendix, Sec. 6. A short description of the 4×4 matrix formalism is given and its application for an electromagnetic analysis of the skin model in Sec. 2.2.4. The results are reviewed and compared to experimental data from literature. Accompanied by a detailed discussion on the true effect of sweat ducts on the reflectivity at normal incidence in the millimeter and lower frequency region we suggest that the sweat ducts do not act as helical antenna in contrast to previous claims^{14, 15} in Sec. 3.1. The preliminary investigation of absorption by water at terahertz frequencies and its influence on both reflectance intensity and generalized ellipsometry parameters are examined, as a contrast mechanism in water content, is presented with the effectiveness of its application for early skin cancer detection.

## 2.

## Skin Structure, Electrical Properties and Modeling

The human skin, the largest organ of the integumentary system, functions as a mediator between human body and the surrounding environment.^{20} Because of its interface with the surroundings, the skin plays a role in protecting the body against pathogens, excessive water loss, radiation, and physical harms, and is responsible for insulation, sensation, and other roles including thermoregulation which is managed mainly by the perspiration system. In light of skin variety of skin functions, it is not surprising that the morphology of skin is a complex one. The skin`s highly complex multilayered structure is interlaced with several subsystems, among them: blood vessels at different sizes, glands, nervous cells, hair and their follicles, fat, sweat glands and ducts. In electromagnetic terms, the skin complexity manifests itself in two main domains: skin`s structural complexity and the unique electrical behavior of the biological matter. In order to examine skin electromagnetic behavior some simplifications are to be introduced in both domains.

## 2.1.

### Dielectric Response of Biological Matter

The past decade has seen a dramatic increase in the exposure of people to electromagnetic fields, sparking many research programs on biological influence resulting from the exposure and its quantification. In addition, fields of science where electromagnetic fields impinge or are used to probe matter, or biomedical fields such as electrophysiology where endogenous bioelectric sources provide fields that are sensed through various body tissues, are affected by the tissue dielectric properties. Any such research must rely on a deep understanding of the dielectric properties of the tissue.^{20} The research on the dielectric properties of biological matter goes back 50 years and is still under intensive work at the present time.

Applying an electrical field on biological matter containing free and bound charges causes them to drift and displace, thus inducing conduction and polarization currents. These dielectric phenomena are determined much by the structure and composition of matter. When approaching the analysis and simplification of human skin electromagnetic behavior, an understanding of its basic components and mechanisms is in order. We assume throughout this work a time dependence of the electromagnetic field of *e*
^{−iω · t}, where ω takes values in the millimeter and terahertz waves range. We rely also on that for most biological matter the magnetic permeability is close to that of free space component of electromagnetic fields at low field strengths.

Formulation of a simplified model for the skin requires defining the dielectric properties of the tissue material that is to characterize its dielectric permittivity and conductivity. Skin tissue is a heterogeneous mixture of lipids, fats, proteins, water etc. Its dielectric spectrum is shaped by the sum of independent dispersion mechanisms, also known as α, β, γ, and δ dispersions, generated by each composite material characterized with specific discrete relaxation parameters (ε_{∞} is the high frequency permittivity, Δε_{n} is the dispersion magnitude, α_{n}-the distribution factor, and τ_{n} is relaxation time constant) as can be seen in Eq. 1 and as is described in more detail in the literature.^{20} An accurate description of the effective complex permittivity of a heterogeneous system can be adequately modeled with a four Cole–Coles and a static conductivity term (σ_{s}):

## Eq. 1

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \varepsilon \left(\omega \right) = \varepsilon _\infty + \sum\limits_{n = 1}^4 {\frac{{\Delta \varepsilon _n }}{{1 + \left({j\omega \tau _n } \right)^{\left({1 - \alpha _n } \right)} }} + \frac{{\sigma _s }}{{j\omega \varepsilon _0 }}}.\end{equation}\end{document} $$\varepsilon \left(\omega \right)={\varepsilon}_{\infty}+\sum _{n=1}^{4}\frac{\Delta {\varepsilon}_{n}}{1+{\left(j\omega {\tau}_{n}\right)}^{\left(1-{\alpha}_{n}\right)}}+\frac{{\sigma}_{s}}{j\omega {\varepsilon}_{0}}.$$^{21}Since we discuss only the millimeter and terahertz wave ranges a simpler model of the complex permittivity than the one presented in Eq. 1 can be applied, thus we will not dwell on its terms. While the other tissue composites are dispersedly inert, water is the main cause for dielectric dispersion at the wavelengths region of interest, causing the γ dispersion.

^{14, 20, 22, 23, 24}A binary mixture model can be applied to describe the complex permittivity of skin,

^{14}where the two constituent materials of the mixture are water [ε

^{w}(

*ω*)] and all other dry skin components (ε

^{bm}). Consequently, the complex permittivity of the mixture can be expressed by the Maxwell-Wagner equation

^{14, 20}[Eq. 2], where ϕ stands for water mixture concentration and ε

^{bm}= 2.2.

^{22}

## Eq. 2

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \varepsilon _{{\rm mix}} = \varepsilon ^{{\rm bm}} \cdot \frac{{2\varepsilon ^{{\rm bm}} + \varepsilon ^w - 2\phi ({\varepsilon ^{{\rm bm}} + \varepsilon ^w })}}{{2\varepsilon ^{{\rm bm}} + \varepsilon ^w + \phi ({\varepsilon ^{{\rm bm}} + \varepsilon ^w })}}. \end{equation}\end{document} $${\varepsilon}_{\mathrm{mix}}={\varepsilon}^{\mathrm{bm}}\xb7\frac{2{\varepsilon}^{\mathrm{bm}}+{\varepsilon}^{w}-2\phi \left({\varepsilon}^{\mathrm{bm}}+{\varepsilon}^{w}\right)}{2{\varepsilon}^{\mathrm{bm}}+{\varepsilon}^{w}+\phi \left({\varepsilon}^{\mathrm{bm}}+{\varepsilon}^{w}\right)}.$$^{25}describing the permittivity of water up to

*f*= 2THz:

## Eq. 3

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \varepsilon ^w \left(\omega \right) = \frac{{\varepsilon _0 - \varepsilon _1 }}{{1 - i\displaystyle \frac{\omega }{{2\pi \cdot \gamma _1 }}}} + \frac{{\varepsilon _1 - \varepsilon _2 }}{{1 - i\displaystyle \frac{\omega }{{2\pi \cdot \gamma _2 }}}} + \varepsilon _2 \end{equation}\end{document} $${\varepsilon}^{w}\left(\omega \right)=\frac{{\varepsilon}_{0}-{\varepsilon}_{1}}{{\displaystyle 1-i\frac{\omega}{2\pi \xb7{\gamma}_{1}}}}+\frac{{\varepsilon}_{1}-{\varepsilon}_{2}}{{\displaystyle 1-i\frac{\omega}{2\pi \xb7{\gamma}_{2}}}}+{\varepsilon}_{2}$$## Eq. 4

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \varepsilon _{{\rm skin}} \left(\omega \right) = \varepsilon _{{\rm mix}} \left({\omega,\phi } \right) + j\frac{{\sigma _s }}{{\omega \cdot \varepsilon _0 }}.\end{equation}\end{document} $${\varepsilon}_{\mathrm{skin}}\left(\omega \right)={\varepsilon}_{\mathrm{mix}}\left(\omega ,\phi \right)+j\frac{{\sigma}_{s}}{\omega \xb7{\varepsilon}_{0}}.$$^{20}results in the conclusion that the presence of water alone cannot account for skin's higher conductivity values. It was recently proposed

^{14}that blood content is the dominant factor in the determination of skin conductivity, ignoring the conductivity contribution of water. This approach however results in too low conductivity values for the upper skin layers, so we suggest a combined model considering the contribution of both water and blood. The following approximation for the complex permittivity of skin is suggested:

## Eq. 5

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \varepsilon _{{\rm skin}} \left(\omega \right) = \varepsilon _{{\rm mix}} \left({\omega,\phi } \right) + j\frac{{\sigma _s }}{{\omega \cdot \varepsilon _0 }} + \rho \cdot \left[ {j\frac{{\sigma ^{{\rm blood}} \left(\omega \right)}}{{\omega \cdot \varepsilon _0 }}} \right]. \end{equation}\end{document} $${\varepsilon}_{\mathrm{skin}}\left(\omega \right)={\varepsilon}_{\mathrm{mix}}\left(\omega ,\phi \right)+j\frac{{\sigma}_{s}}{\omega \xb7{\varepsilon}_{0}}+\rho \xb7\left[j\frac{{\sigma}^{\mathrm{blood}}\left(\omega \right)}{\omega \xb7{\varepsilon}_{0}}\right].$$Blood conductivity values were calculated based on experimentally verified model results^{21} with 0 ⩽ ρ < 1.

## 2.2.

### Skin Structure

After setting up a model for the dielectric properties of the human skin, we turn to deal with the challenges presented by the skin structural complexity.

## 2.2.1.

#### Skin layered nature

Human skin is a multilayer structure composed of three main layers—from the outermost to the innermost: the stratum corneum (SC), the epidermis and the dermis. The layered nature of skin allows us to design an appropriate multilayer model with specific structural and electrical properties of each isotropic subdomain. Thus, for simplicity, each skin layer can be defined by its thickness and its complex permittivity (Table 1). Different body regions are characterized by different skin layer widths and water content values, yet changes within a specific region are slow in normal tissue. Since the skin abnormalities of interest demonstrating water content changes (and morphological changes) are localized and have defined boundaries, they are expected to present high gradient in water concentration and thus contrast in reflectivity or generalized elliposmetry (GE) based images as will be discussed later. Throughout the paper the skin layer widths and water concentration values used are characteristic to the human palm, for comparison with available experimental measurements. Yet the specific values or skin region are of lesser importance than the ability to differentiate between normal water content in healthy tissue and the elevated water content in abnormal skin sites. For illustrative purposes and for results comparison the skin parameters used by Feldman
^{14} will be discussed. In dielectric terms, different skin layers are characterized by different bulk water contents reflecting on the 1st term in Eq. 5, different static conductivity reflecting on the 2nd term, and different blood content values reflecting on the 3rd term. The model described here is in close agreement with the complex permittivity values for the dermis and epidermis layers presented by Gabriel
^{21} and Feldman
^{14} for the frequency region 75–110 GHz. However it differs by a higher conductivity value for the SC layer in comparison to the value of 10^{−5} s m^{−1} derived by Feldman`s model, suggesting a very severe case of dry skin. In addition, our model allows validity for a wider frequency range (up to the terahertz region) and it takes into account the dispersive nature of the effective conductivity which was considered negligible by Feldman
^{14} when applying the coarse consideration of blood content in each layer. The values of the variables of Eq. 5 used for each layer are presented in Table 2 and the complex permittivity for each skin layer in Fig. 2.

## Table 1

Skin layered-model parameters: permittivity, ac conductivity, and thickness.

Layer | εr at 100 GHz | σ at 100 GHz (S m−1) | Thickness (μm) |
---|---|---|---|

SC | 2.4 | 10^{−5} | 30 |

Epidermis | 3.2 | 1 | 350 |

Dermis | 3.9 | 30 | 1000 |

Subcutanous | 3.11 | 12.31 | Semi-infinite |

HSD | 4 | 3500 | 350 |

## Table 2

Conductivity and water and blood content values for each skin layer used in the dielectric model.

Layer | SC | Epidermis | Dermis | Subcutaneous tissue |
---|---|---|---|---|

ϕ(%) | 3.7 | 16.6 | 26 | 21 |

ρ(%) | 0 | 2 | 40 | 16 |

σ_{S} | 10^{−5} | 0.025 | 0.2 | 0.08 |

## 2.2.2.

#### Sweat ducts

Among other functions, the skin is appointed to manage the body`s thermoregulation utilizing the perspiration system incorporated within it. The main building blocks of the perspiration system are the sweat glands and the sweat ducts. The sweat ducts run through the epidermis to the skin surface and are responsible for the delivery of the sweat from the sweat glands embedded in the dermis.^{26} It is well known^{27} since the middle of the last century that the end part of the sweat ducts has a helical shape and it has been confirmed recently^{28, 29} using high resolution images of optical coherence tomography. The combination of the geometric shape and the high conductivity of the ducts are intriguing since they may possibly present a strong electromagnetic response that is sensitive to structural changes. This sensitivity might be utilized for the detection of skin cancer related changes.^{10} Based on the dielectric properties of the human sweat ducts described by^{14} we embed sweat ducts into our simplified layered skin model, thus breaking the isotropy of the Epidermis layer.

## 2.2.3.

#### Modeling sweat ducts using homogenization theory

After establishing a structurally and dielectrically simplified model of the human skin we continue to the examination of its electromagnetic response at a variety of incidence angles, frequencies and polarizations. A rigorous analysis of the model, including the epidermis layer embedded with sweat ducts, coerces one to analytically extract the electromagnetic (EM) fields using the Maxwell's equations in a geometrically complex structure, a very difficult and non elegant task. The homogenization process is proven to be useful in simplifying the investigation of the electromagnetic behavior of complex composite structures for wavelengths much longer than the scale of the inhomogenieties. Implementing the Bruggeman formalism for homogenization provides a prediction for the behavior of the *effective* electromagnetic response of an orderly textured particulate composite media (made up from a host material and inclusions in a relatively small concentration). Having an average density of 150 to 340 per square centimeter of skin,^{30} and being homogeneously spread throughout the epidermis layer, the overall structure of sweat ducts embedded in the epidermis satisfy the conditions for the implementation of the homogenization algorithm. In the following paragraph a more general road-map is provided describing the use of homogenization while a more elaborated description of the mathematics behind homogenization can be found in the Appendix, Sec. 6.

In order to model the helical sweat duct, we treat it as if it is comprised of a collection of very small sequential finite cylinders placed along the center line of the helix [Fig. 3a] since the depolarization dyadic (required for homogenization) for cylinders are easier to calculate.^{32, 33} The epidermis layer is divided into sub-layers, each containing cylinders pointing at an appropriate orientation to their location along the helix. The Bruggeman formalism is used to deduce the effective permittivity dyadic of the single planar sublayer containing cylinders pointing at an orientation defined by the polar and azimuth angles θ and, ϕ respectively. The geometric nature of the cylindrical inclusions determine the nature of the anisotropy of each effective sub-layer, in essence the orientation angles define the direction of the principal axis of the sublayer's dielectric tensor. In our case ε_{3} is the principal dielectric constant along the cylinder axis while along its radius we have: ε_{1} = ε_{2}. A simulation for the electrical behavior of the entire helix is established by stacking a set of these planes one on top of the other. The permittivity dyadic for each plane in the xyz frame is determined by orientation of the appropriate cylindrical slice.^{33} By implementing this approach, the structure embedded with the helical inclusions is simplified to a stack of dielectric anisotropic (uniaxial) layers, a problem well investigated^{34, 35} and easy to approach similar to the simulation of the optics of helical liquid crystals^{36} and chiral sculptured thin films.^{37, 38}

## 2.2.4.

#### Application of the 4×4 matrix formalism for EM analysis of skin model

An accurate analysis of the skin model electromagnetic response at a variety of incidence angles, frequencies and polarizations can be achieved by applying the 4×4 matrix formalism first described by Berreman,^{39} and later further developed by Abdulhalim.^{34} For simplicity, we define
[TeX:]
$\hat z$
$\widehat{z}$
as the direction of propagation at normal incidence and
[TeX:]
$\vec \psi ^T = (\!{\begin{array}{*{20}c}{\sqrt {\varepsilon _0 } E_x,} & {\sqrt {\mu _0 } H_y,} & {\sqrt {\varepsilon _0 } E_y,} & { - \sqrt {\mu _0 } H_x }\\\end{array}}\!)$
${\stackrel{\u20d7}{\psi}}^{T}=\left(\begin{array}{cccc}\sqrt{{\varepsilon}_{0}}{E}_{x},& \sqrt{{\mu}_{0}}{H}_{y},& \sqrt{{\varepsilon}_{0}}{E}_{y},& -\sqrt{{\mu}_{0}}{H}_{x}\end{array}\right)$
, where E and H represent the electric and magnetic fields. The 4×4 matrix formalism uses an algebraic first order differential equation representation of Maxwell's equations [Eq. 6] to describe the propagation of the electromagnetic field in a homogeneous anisotropic media:

## Eq. 6

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \frac{{\partial \vec \psi }}{{\partial z}} = ik_0 \underline{\underline \Delta } \vec \psi.\end{equation}\end{document} $$\frac{\partial \stackrel{\u20d7}{\psi}}{\partial z}=i{k}_{0}\underline{\underline{\Delta}}\stackrel{\u20d7}{\psi}.$$The elements of
[TeX:]
$\underline{\underline \Delta} $
$\underline{\underline{\Delta}}$
were described in detail in Ref. 34. The general formal solution to Eq. 6 that finds the field components after a travel of distance *h* inside a homogeneous medium is given by

## Eq. 7

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \vec \psi \left({z + h} \right) = e^{ik_0 h\underline{\underline \Delta } } \cdot \vec \psi \left(z \right) = \underline{\underline P} \left(h \right) \cdot \vec \psi \left(z \right). \end{equation}\end{document} $$\stackrel{\u20d7}{\psi}\left(z+h\right)={e}^{i{k}_{0}h\underline{\underline{\Delta}}}\xb7\stackrel{\u20d7}{\psi}\left(z\right)=\underline{\underline{P}}\left(h\right)\xb7\stackrel{\u20d7}{\psi}\left(z\right).$$*i*of width

*h*

_{i}one can derive the equivalent propagation matrix for the entire stratified model. The application of the 4×4 matrix formalism as described by Abdulhalim

^{34}provides a simple method for extracting the reflection matrix [Eq. 8] from the equivalent propagation matrix relating the amplitudes of the incident and reflected electric field for

*s*and

*p*polarizations for the multilayered anisotropic skin model.

## Eq. 8

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \left({\begin{array}{*{20}c} {A_{{\rm rp}} } \\ {A_{{\rm rs}} } \\ \end{array}} \right) = \left({\begin{array}{*{20}c} {r_{{\rm pp}} } &\quad {r_{{\rm ps}} } \\ {r_{{\rm sp}} } &\quad {r_{{\rm ss}} } \\ \end{array}} \right) \cdot \left({\begin{array}{*{20}c} {A_{{\rm ip}} } \\ {A_{{\rm is}} } \\ \end{array}} \right). \end{equation}\end{document} $$\left(\begin{array}{c}{A}_{\mathrm{rp}}\\ {A}_{\mathrm{rs}}\end{array}\right)=\left(\begin{array}{cc}{r}_{\mathrm{pp}}& \phantom{\rule{1em}{0ex}}{r}_{\mathrm{ps}}\\ {r}_{\mathrm{sp}}& \phantom{\rule{1em}{0ex}}{r}_{\mathrm{ss}}\end{array}\right)\xb7\left(\begin{array}{c}{A}_{\mathrm{ip}}\\ {A}_{\mathrm{is}}\end{array}\right).$$*R*

_{ij}= |

*r*

_{ij}|

^{2}with

*i*,

*j*=

*s*,

*p*represent the reflectivity when the incident polarization is along

*i*and the analyzer axis is along

*j*. A 4×4 matrix based Matlab program was used to calculate and cascade the propagation matrices for each of the skin layers and epidermis sublayers for the sweat duct embedded skin model. Simulating the propagation and reflection of millimeter and terahertz plane waves impacting on the skin at different incidence angles allowed deriving the reflection matrix. Using the values of calculated reflectivity amplitudes for

*p*and

*s*polarizations and their conjugations upon reflection (

*r*

_{pp},

*r*

_{ss},

*r*

_{ps},

*r*

_{sp}) further studies of the spectral and incidence angle dependence of reflectivity magnitude and GE variables

^{40}have been made as described in Secs. 3.2, 4. A similar approach based on Berreman formulation was presented previously by Schubert

^{41}to describe GE from helical anisotropic materials. A laboratory measurement setup suited to experimentally verify the calculated reflectivity amplitude values can be in principle similar to the one presented in Fig. 1, although in a true system additional sophistications are embedded.

As with any multilayered system, phase retardations and polarization changes play a dominant role in determining the layers’ parameters such as using reflectometry, polarimetry, or ellipsometry. Therefore, we have decided to perform calculations of the ellipsometric parameters as well. The GE variables can be calculated from the reflectivity matrix based on the following definitions:^{42}

## Eq. 9

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \begin{array}{*{20}c} {\displaystyle\frac{{r_{{\rm pp}} }}{{r_{{\rm ss}} }} = \tan \psi _{{\rm pp}} \cdot e^{i\Delta _{{\rm pp}} } } \\[10pt] {\displaystyle\frac{{r_{{\rm ps}} }}{{r_{{\rm pp}} }} = \tan \psi _{{\rm ps}} \cdot e^{i\Delta _{{\rm ps}} } }. \\[10pt] {\displaystyle\frac{{r_{{\rm sp}} }}{{r_{{\rm ss}} }} = \tan \psi _{{\rm sp}} \cdot e^{i\Delta _{{\rm sp}} } } \\ \end{array} \end{equation}\end{document} $$\begin{array}{c}{\displaystyle \frac{{r}_{\mathrm{pp}}}{{r}_{\mathrm{ss}}}=\mathrm{tan}{\psi}_{\mathrm{pp}}\xb7{e}^{i{\Delta}_{\mathrm{pp}}}}\\ {\displaystyle \frac{{r}_{\mathrm{ps}}}{{r}_{\mathrm{pp}}}=\mathrm{tan}{\psi}_{\mathrm{ps}}\xb7{e}^{i{\Delta}_{\mathrm{ps}}}}.\\ {\displaystyle \frac{{r}_{\mathrm{sp}}}{{r}_{\mathrm{ss}}}=\mathrm{tan}{\psi}_{\mathrm{sp}}\xb7{e}^{i{\Delta}_{\mathrm{sp}}}}\end{array}$$## 3.

## Results

## 3.1.

### Effect of the Helical Sweat Glands on the Terahertz and Millimeter Waves Response

Recently, some studies^{14, 15} have claimed that the skin helical sweat ducts act as an array of low-*Q* helical antennae and are dominant in shaping the spectral response at the sub-terahertz region and have used experimental results to try and support their theory. Based on the above skin model and the simulation results, we suggest and demonstrate that interference due to multiple reflections (etaloning effect) from the skin layered structure is in fact the dominant player in shaping the spectral response of the skin and not the sweat ducts as claimed. Particularly, there is no need for strange assumptions such as having the sweat ducts acting as low-*Q* helical antennae. The importance of this comparison is in the apparent similarity in the spectral behavior of the experimental measurements presented by Feldman and the results of our simulation, giving it an experimental validation. In order to present a one-to-one comparison, we use the structural and dielectric parameters (Table 1) of the skin layers and sweat ducts evaluated and described in more details by the same previous works.^{14, 15} Although the approach presented there for deriving the dielectric properties of skin layers is simplified in comparison to the method described by us in Sec. 2.1, as have been explained earlier, the results of both methods still do not differ significantly in the narrow spectral region of 75 to 110 GHz. However our method allows us to extend the spectral range further.

Observing the structural and electrical characteristics of the skin, the assumption that the helical sweat ducts act as low-Q helical antennae^{14, 15} has been based on the following arguments: 1. The sweat ducts are embedded in the epidermis, having smaller permittivity than the Dermis; 2. High ac-conductivity in the extremely high frequency (EHF) domain due to the mechanism of fast proton-hopping through the distributed H-bond networks along the duct surface; 3. A potential drop caused by the difference in the *p*H value between the skin surface and the dermis. The reflection coefficient`s spectral response of skin tissue was claimed^{14, 15} to be dominated by the electromagnetic response of the low-*Q* helical antennae-like behavior of the sweat ducts. These experiments showed that the spectral response of the reflectance coefficient from the skin varies in agreement with stimuli variable changes with the amount of sweat conducted by the sweat ducts (in essence, variations in their conductivity—resulting in changes in the effective conductivity of the epidermis and affecting absorption and dissipation in the entire skin).

In order to demonstrate the dominance of the effect of interference on the spectral response within the simplified multilayered model of the skin, our results are presented in three stages as shown in Fig. 4a, 4b, 4c: (a) isotropic nonconducting skin layers, (b) isotropic conducting layers with increasing conductivity, and (c) the skin model embedded with the helical sweat ducts, with the embodiment of the helical sweat ducts in the epidermis resulting in breaking the isotropy of the epidermis layer. Normal incidence is assumed with air as an ambient material, and as a substrate (representing the tissue under the dermis) a homogeneous tissue with dielectric properties as presented in Table 1 were used.

Starting from the case of lossless media one can clearly see in Fig. 4a the interference pattern of the reflected signal. When introducing conductivity into the model, the electrical field decays along its propagation path due to the frequency dependent losses. By gradually increasing the conductivity from 10% to 100% of the conductivity (σ) values of each layer presented at Table 1, we can see the gradually increasing damping of the oscillations of the interference pattern and the convergence to a single lower minimum at about ƒ_{0} = 92.5 GHz caused by the increasing losses [Fig. 4b]. The sweat duct coils were then introduced to the Epidermis. In Fig. 4c one can see that the spectral response of the skin model with the sweat ducts presents a similar behavior as the model without the sweat ducts. It is thus obvious that the sweat ducts play no significant role in shaping the spectral response. The only main observable difference is the shift in the location of the minimum point in the reflectance that coincides with the one presented by Feldman,^{14, 15} ƒ_{0} ∼ 92 GHz. Hence we propose that multiple interferences are the mechanism responsible for the origin of the skin spectral response seen in Fig. 4c.

The question that may arise, however, is whether the location of the minimum is dominated by the interference effect or by the response of the sweat ducts characterized by their parameters (size and substance). To further help clarify this point the thickness of one of the layers was changed resulting in the relocation of the minimum of the spectral response function, thus proving again that the spectral response shape is interference dominated and that the effect of the sweat ducts is secondary. Figure 5 illustrates the spectral response of the skin model having different SC thickness, first (a) without the presence of sweat ducts and (b) with the sweat ducts present. Observing the spectral response in Fig. 5a, our expectation is approved that the spectral response of the layered structure, without the sweat ducts, is interference dominated.

Observing the spectral response in Fig. 5b we can see that the overall spectral response is negatively-shifted in frequency due to the introduction of the sweat glands. However when changing the SC thickness, the minimum location shifts in resemblance to the case in Fig. 5a. This shift suggests that the effect caused by the introduction of the sweat glands is a secondary player in determining the location of the minimum in the spectral response and that this minimum is not a resonance of a conducting helical antenna.

To summarize, the dominant effect in shaping the spectral response in the millimeter and sub-terahertz region seems to be in fact the interference, rather than an effect caused by the presence of the helical sweat ducts as low-Q helical antennae. Though sweat ducts do not appear to have major effect on the intensity of reflected field, the introduction of sweat ducts to the skin model results in breaking the isotropy of the model. Nonisotropic media results in coupling between EM field polarizations, suggesting an ellipsometric analysis of the skin model. In addition, we have used the experimental measurements presented by Feldman to validate the results of our skin model.

## 3.2.

### Reflectivity and Generalized Ellipsometry Variables Dependence on Frequency and Incidence Angle

The dependence of reflectance coefficients *R*
_{pp}, *R*
_{ss}, *R*
_{sp}, and *R*
_{ps} on frequency and incidence angle are presented in Figure 6. One can notice the identical behavior of *R*
_{pp} and *R*
_{ss} at normal incidence and the reflectance dip, due to interference, at about 92 GHz as was experimentally demonstrated by Feldman
^{14} The crossed polarization reflectivities *R*
_{sp}, *R*
_{ps} are a result of the anisotropic helical sweat ducts but they are too small to be practical for any diagnostic technique. Hence the conclusion is that the helical sweat glands, although, they produce some anisotropy and crossed reflectivity signal, their concentration is so small to give any noticeable effect even at oblique incidence.

Figure 7 presents the dependence of the generalized ellipsometry variables on frequency and incidence angle. The spectral and angular behavior of field reflectance intensity and GE variables are not very informative for the purpose of diagnosing skin cancer when examining them in any specific state of tissue water content. However an examination of the relative changes in each of the variables as a function of water content variation can be useful for diagnostic purposes.

Our goal is to examine and present the potential of possible novel concept for skin cancer detection using terahertz ellipsometric imaging. Since we discuss the theoretical ideal model only, in this early phase of the research we do not deal yet with measurement errors and it's nonuniformity along the spectrum or angular region of interest. These system considerations, although important for fitting, will be considered once experimental measurements will be performed.

## 4.

## Absorption by Water at Terahertz Frequencies as a Contrast Enhancement Mechanism

Recent studies of terahertz imaging have revealed a significant difference between skin cancer and healthy tissue.^{44, 45, 46} Terahertz radiation sensitivity to polar molecules results in strong absorption by water.^{47} Since tumors tend to have different water content than normal tissue,^{48, 49} variation in water content between tumor tissue and the surrounding normal tissue can serve as a likely “contrast” mechanism suitable for early skin cancer detection. Our research goal is to examine the sensitivity to water content changes on not only reflectance intensity, as previous studies on terahertz imaging in reflectance geometry, but also the sensitivity of the polarization changes as can be deduced by an investigation of the GE variables. In addition, in order to determine the ability of a possible future cancer diagnosis imaging setup to detect tumor tissue at an earliest possible stage, we examine the limit of detection of water content changes by both reflectance intensity and GE variables measurements. Since the presented imaging setup can vary both the frequency and the incidence angle of the impinging terahertz waves we move to determine the frequency and incidence angle at which the limit of detection is the lowest.

## 4.1.

### Preface—Skin Cancer

The development of skin cancer begins in the epidermis causing an abnormal cell division in one of the epidermis composite cell types: Squamous cells, basal cells or melanocytes resulting in squamous cell carcinoma, BCC or Melanoma, accordingly. All cancer types tend to spread to neighbor skin layers and may potentially spread to other organs at a different rate for each cancer type. With the disruption of normal cell division at malignant tumor sites, localized elevated water and blood contents are felt due to the production of irregular cell masses. We assume that the water content changes uniformly throughout the skin layers, and examine the changes in the reflected field intensity and polarization at low water content changes (at the early stages of cancer).

## 4.2.

### Results

Figure 8 demonstrates the changes in *R*
_{pp} and *R*
_{ss} for tissue with elevated water content (+4%, +20%, and +45%) relative to the results at normal tissue water content (Fig. 6).

Figure 9 demonstrates the gradual change in tan ψ_{pp}, tan ψ_{ps}, and tan ψ_{sp}, relative to the results for normal tissue water contents (Fig. 7) at the elevated water content tissue (+4%,+20% and +35%) and Fig. 10 demonstrates the gradual change in cos Δ_{pp}, cos Δ_{ps} and cos Δ_{sp}.

Observing Figs.
8, 9, 10 it seems that the amplitude of the change in *R*
_{pp}, *R*
_{ss}, tan ψ_{pp}, tan ψ_{ps}, tan ψ_{sp}, cos Δ_{pp}, cos Δ_{ps}, and cos Δ_{sp} is monotonous with the rise in water content. That is, a tissue suffering from elevated water content can be distinguished from normal tissue by a comparison based on one of the measured values previously mentioned. Since “normal” water content varies slightly from one person to the other, normal water content can be defined as the one compatible with the behavior presented by the majority of the skin area examined. In other words, non-normal tissue will be prominent in its presented measured values relative to the majority of surrounding tissue in such a manner that localized and finite in size abnormalities would be easily detected based on contrast in the terahertz image of the tissue examined.

One can also conclude by observing Figs.
8, 9, 10 that each variable presents sensitivity to changes in water content at a specific range of frequencies and incidence angles. For example, observing Fig. 8, one can notice that the increase in water content effect on changes both in *R*
_{pp} and *R*
_{ss} is mainly noticeable at the lower frequency region 10 to 500 GHz and at all incidence angles. Observing Fig. 9 and Fig. 10, the increase in water content effect on changes in GE variables is significant at different frequency and incidence angle for each variable.

## 4.3.

### Spectrum Sensitivity and Lower Detection Limit of Water Content Variation

In practice, an imaging setup is designed to analyze a suspicious tissue by either performing a scan at different incidence angles at a fixed frequency or by performing a spectral scan at a fixed angle. Fixing the frequency or angle at each of the methods must be efficient in the sense of satisfying the optimized conditions for performing the measurement. Optimized conditions are met by positioning the fixed frequency or angle at a value demonstrating the highest sensitivity—which enables detection of the smallest possible variation in water content, although systems for variable angle spectroscopic ellipsometry exist these days.^{50} We examine the sensitivity in frequency and incidence angle for each of the variables (reflection intensity and GE) by calculating root mean square error (RMSE) values separately for each [Eq. 10]. For example, RMSE(Δ*C*, θ)[*X*(θ, *f*)] is the RMSE for one of the frequency, angle or water content dependent variables examined in Figs.
8, 9, 10, here noted by *X*, for a variation in water concentration Δ*C* at a given incidence angle θ. *N*stands for the number of samples in frequency or angular spectra. By performing summation on angles instead of the frequencies, we can similarly define RMSE(Δ*C*, *f*)[*X*(θ, *f*)] for a variation in water concentration Δ*C* at a given frequency*f*.

## Eq. 10

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{eqnarray} &&\hspace*{-10pt} {\rm RMSE}\left({\Delta C,\theta } \right)\left[ {X\left({\theta,f} \right)} \right]\nonumber\\ && = \sqrt {\frac{{\sum\limits_{i = 1}^N {\left[ {X\left({\theta,f_i,\Delta C = 0} \right) - X\left({\theta,f_i,\Delta C} \right)} \right]^2 } }}{{N - 1}}}. \end{eqnarray}\end{document} $$\begin{array}{ccc}& & \mathrm{RMSE}\left(\Delta C,\theta \right)\left[X\left(\theta ,f\right)\right]\hfill \\ & & =\sqrt{\frac{\sum _{i=1}^{N}{\left[X\left(\theta ,{f}_{i},\Delta C=0\right)-X\left(\theta ,{f}_{i},\Delta C\right)\right]}^{2}}{N-1}}.\hfill \end{array}$$The RMSE(Δ*C*, θ)[*X*(θ, *f*)] is calculated for a set of water concentration variations and for all angles. Similarly the RMSE(Δ*C*, *f*)[*X*(θ, *f*)] is calculated for a set of water concentration variations and for all frequencies. The frequency or incidence angle that presents the highest sensitivity to water content variations is defined as the one that presents the highest gradient in water concentration values around the normal tissue water concentration valueΔ*C* = 0. The results are presented in Fig. 11.

The sensitivity of each variable to small water content variations is defined as the slope around Δ*C* = 0 of the RMSE (frequency or angle) plot for different water content variations (Fig. 11). Note that the variables (*R*
_{pp}, *R*
_{ss}, tan ψ_{pp}, tan ψ_{ps}, tan ψ_{sp}, cos Δ_{pp} and cos Δ_{sp}) present a similar behavior for small water content variations: the RMSE for frequency scan at a fixed angle (marked with circles) demonstrates a slower inclination then RMSE for incidence angle scan at a fixed frequency (marked with triangles). In other words, each variable is expected to present higher sensitivity to small water content variations when examined by an imaging system at a fixed frequency and performing an angular scan. The variable cos Δ_{ps} however demonstrates an opposite behavior, displaying a higher sensitivity for spectral scan. Figure 11 also illustrates an almost linear behavior of the frequency RMSE plots at small water content variations.

The lowest detectable value of water content variation is an important parameter for quantifying the ability of a reflectance or ellipsometry-based imaging setup to detect delicate changes in water content—hence early skin cancer detection. The limit of detection (LOD), for a specific variable is defined as follows:

## Eq. 11

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} {\rm LOD} = \left\langle {\delta C} \right\rangle = \frac{{\left\langle {\delta \left({{\rm RMSE}} \right)} \right\rangle }}{S}, \end{equation}\end{document} $$\mathrm{LOD}=\u2329\delta C\u232a=\frac{\u2329\delta \left(\mathrm{RMSE}\right)\u232a}{S},$$*S*= Δ(RMSE)/Δ

*C*is the sensitivity.

Standard ellipsometric systems are capable of measurements with a repeatability of 〈δ(RMSE)_{R}〉 = 0.001 in reflectance intensity as well as a repeatability^{51} of 〈δ(RMSE)_{α}〉 = 0.001, where α stands for each of the ellipsometric variables tan ψ_{pp}, tan ψ_{ps}, tan ψ_{sp}, cos Δ_{pp}, cos Δ_{ps}or cos Δ_{sp}. LOD values for all variables are presented in Fig. 11 and Table 3. Optimal frequency values demonstrating the highest sensitivity for water content variations are summarized at Table 3. Optimal incidence angle values for imaging by performing a spectral scan at a fixed incidence angle and the corresponding LOD values can be viewed in Table 3 as well.

## Table 3

LOD in percentage estimated using spectral (LODs) or angular (LODa) scans for R pp, R ss, tan ψpp, tan ψps,tan ψsp,cos Δpp, cos Δps, and cos Δsp.

Variable | Rpp | Rss | tan ψpp | cos Δpp | tan ψps | cos Δps | tan ψsp | cos Δsp |
---|---|---|---|---|---|---|---|---|

f_{optimal} (GHz) | 30 | 120 | 90 | 90 | 1830 | 1600 | 1980 | 80 |

LODa (%) | 1.54 | 0.8 | 0.24 | 0.063 | 0.077 | 0.36 | 4.33 | 0.11 |

θ_{optimal} | 80.95 deg | 72.81 deg | 63.77 deg | 55.18 deg | 55.63 deg | 57.44 deg | 0 deg | 80.52 deg |

LODs (%) | 4.22 | 2.11 | 0.99 | 0.096 | 0.1 | 0.12 | 8.457 | 0.39 |

## 4.4.

### Discussion

An examination of Table 3 illustrates that LOD has consistently lower values when performing an angular scan (except for cos Δ_{ps}), in agreement with the initial examination of sensitivity. We, therefore, suggest that the preferred geometry for the imaging system is reflectance geometry, setting a fixed frequency, and performing an angular scan. One can also notice that an optimal imaging setup configuration for both spectral and angular scan takes different values for each variable suggesting that an imaging system designed for skin cancer diagnosis should be robust by allowing the acquisition of measurements at a wide spread of angles or, if preferred, a wide frequency spectrum.

Examining the LOD values in Table 3 one learns that the GE variables present lower LOD values than reflectance intensity variables *R*
_{pp} and *R*
_{ss} even up to an order of magnitude. This is due to two reasons: first their higher sensitivity as they contain the phase information and second because ellipsometry is a self referenced technique and so its precision is much higher than reflectivity measurement. Malignant tissue presents water content elevated up to ∼5%, over that of normal tissue based on nuclear magnetic resonance (NMR) imaging,^{52} refractometry,^{45} infrared spectroscopy,^{53} terahertz pulse imaging^{54} and Raman scattering^{55} observations. This overhydration is higher than most of the detection limits presented in Table 3, hence we conclude that both reflectometry and generalized ellipsometry are suitable for detection of skin cancer based on delicate variations in the water content, however GE is more recommended. Note that according to Table 3 the best parameter to be used is cos Δ_{pp} as it exhibits the lowest LOD at the optimum frequency of 90 GHz in the angular scan and 55 deg incidence angle in the frequency scan. Experimentally however measurement accuracy can depend on the frequency or on the angle, a fact that can reduce the sensitivity predicted above. However, it is believed that even with typical ellipsometric systems, high enough sensitivity will be obtained to allow detecting skin cancer in relatively early stages.

## 5.

## Summary

Using the 4×4 matrix formalism and homogenization theory we have established a simplified model for skin tissue designed to integrate skin's layered and anisotropic nature and tailored for the millimeter to terahertz frequency domain. Presenting a comparison to experimental measurements^{14} we have verified the validity of our model and suggested that sweat ducts do not serve as the dominant player in shaping skin's spectral response but its layered structure.

Using the reflectance matrix components the polarized reflectometry and GE parameters were examined. By relating water content changes associated with malignant skin cancer tumors a potential was demonstrated for differentiating normal tissue from cancerous tissue by a simple measurement of the resulting variations in reflectance intensity and GE variables.

Assessing the LOD for each variable by spectral scan or angular scan we have illustrated the potential of detecting even delicate water content variations in tissue well below the variations caused by advanced tumors. An ellipsometric measurement setup is recommended performing an angular scan at a fixed frequency with the parameter that gives the lowest LOD iscos Δ_{pp}. The expected LOD for water content variations for each variable was presented both by spectral scan, and angular scan and the optimal incidence angles and measurement frequencies suited for their optimal measurement were calculated.

We have demonstrated the potential of a novel technique for skin cancer detection in the millimeter and terahertz regions, based on a simplistic and ideal skin model, with some comparison to experimental measurements. Experimental results are expected to vary slightly from the calculated theoretical values presented here, and summarized in Table 3. For a better fit to the realistic case, future research should examine an improved skin model coping with scattering mechanisms and more realistic measurement system challenges.

## 6.

## Appendix: Homogenization—Mathematical Description

## 6.1.

### Preface

The homogenization theory is described here only shortly; while for a more thorough description of the theory the reader is referred to the work by Lakhtakia
^{32} Dealing with Maxwell's equations in dielectric composite materials of more than one different dielectric properties can be a mathematically nontrivial task. Dielectric particulate composite materials in different fields of research, from thin films and liquid crystals to biological matter, can present high complexity in both dielectric properties (dispersive and nonisotropic materials), the geometrical shapes of each component material in the mixture and their spatial distribution—all of which should be taken into account when calculating the EM properties of the material. The common classical electromagnetism provides us with mathematical tools for calculating Maxwell's equations in media with a small number of inclusions usually in the form of Green functions theory. However, its application for the case of media seeded with a large number of inclusions, a large number of different types of inclusions or inclusions with complex shapes is nontrivial and can be excruciating. In short, the homogenization theory, applicable under some conditions (that will be discussed later), transforms complex particulate composite media to dielectrically equivalent homogeneous material–homogenized composite material (HCM). The homogenization process is based on mathematical manipulation of the differential Maxwell's equations and the EM field continuity constraints taking into account the geometrical characteristics of the composite material. We shall discuss only the Bruggeman formalism for homogenization since it has been shown to provide more accurate results.^{31}

## 6.2.

### General Homogenization Approach

Consider a medium composed of two dielectrically different materials—the host material and the inclusions. The inclusions, having distinct geometrical shapes and dielectric tensor
[TeX:]
$\underline{\underline {\varepsilon _{{\rm inc}} }}$
$\underline{\underline{{\varepsilon}_{\mathrm{inc}}}}$
are uniformly scattered throughout the media with a relative concentration denoted by*f*
_{inc}. The host material, filling the gap between inclusions, possesses a dielectric tensor
[TeX:]
$\underline{\underline {\varepsilon _h }}$
$\underline{\underline{{\varepsilon}_{h}}}$
and a relative concentration of *f*
_{h} = 1 − *f*
_{inc}. In our case, homogenizing sweat duct cylindrical slices (inclusions) in the Epidermis layer (host), both inclusions (uniaxial) and host (isotropic) material are dispersive. Since sweat ducts are scattered throughout the skin with an average density of 150 to 340 per square centimeter of skin^{33} combining the geometrical parameters of the sweat duct (*r*
_{2} = 20 μm, Fig. 3) one obtains*f*
_{inc}≅0.003. Considering the geometrical symmetry in the approached particulate composite media, the only unique direction is the axis of the cylinders. Since the geometrical shape of the inclusions play a role in the HCM dielectric characteristics, as will also be seen in Sec. 6.3, one can expect the HCM to be a uniaxial material due to the uniaxiality of the cylinders. We denote the cylindrical axis direction with
[TeX:]
$\hat c$
$\widehat{c}$
and define
[TeX:]
$\underline{\underline {I_t }} = \underline{\underline I} - \hat c\hat c$
$\underline{\underline{{I}_{t}}}=\underline{\underline{I}}-\widehat{c}\widehat{c}$
, where
[TeX:]
$\underline{\underline I}$
$\underline{\underline{I}}$
is the idempotent. The dielectric tensor of the host and inclusion materials can also be written in the following form:
[TeX:]
$\underline{\underline {\varepsilon _h }} = \varepsilon _t^h \underline{\underline {I_t }} + \varepsilon _c^h \hat c\hat c{\rm },{\rm }\underline{\underline {\varepsilon _{{\rm inc}} }} = \varepsilon _t^{{\rm inc}} \underline{\underline {I_t }} + \varepsilon _c^{{\rm inc}} \hat c\hat c$
$\underline{\underline{{\varepsilon}_{h}}}={\varepsilon}_{t}^{h}\underline{\underline{{I}_{t}}}+{\varepsilon}_{c}^{h}\widehat{c}\widehat{c},\underline{\underline{{\varepsilon}_{\mathrm{inc}}}}={\varepsilon}_{t}^{\mathrm{inc}}\underline{\underline{{I}_{t}}}+{\varepsilon}_{c}^{\mathrm{inc}}\widehat{c}\widehat{c}$
and the electric displacement field inside the host [Eq. 12] or inside the inclusions [Eq. 13] as:

## Eq. 12

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \vec D ({\vec r }) = \varepsilon _0 \underline{\underline {\varepsilon _r^h }} \cdot \vec E \left({\vec r } \right) = \varepsilon _0 \big[ {\varepsilon _t^h \underline{\underline {I_t }} + \varepsilon _c^h \hat c\hat c} \big] \cdot \vec E ({\vec r }), \end{equation}\end{document} $$\stackrel{\u20d7}{D}\left(\stackrel{\u20d7}{r}\right)={\varepsilon}_{0}\underline{\underline{{\varepsilon}_{r}^{h}}}\xb7\stackrel{\u20d7}{E}\left(\stackrel{\u20d7}{r}\right)={\varepsilon}_{0}\left[{\varepsilon}_{t}^{h}\underline{\underline{{I}_{t}}}+{\varepsilon}_{c}^{h}\widehat{c}\widehat{c}\right]\xb7\stackrel{\u20d7}{E}\left(\stackrel{\u20d7}{r}\right),$$## Eq. 13

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \vec D ({\vec r }) = \varepsilon _0 \underline{\underline {\varepsilon _r^{{\rm inc}} }} \cdot \vec E ({\vec r }) = \varepsilon _0 \big[ {\varepsilon _t^{{\rm inc}} \underline{\underline {I_t }} + \varepsilon _c^{{\rm inc}} \hat c\hat c} \big] \cdot \vec E ({\vec r }). \end{equation}\end{document} $$\stackrel{\u20d7}{D}\left(\stackrel{\u20d7}{r}\right)={\varepsilon}_{0}\underline{\underline{{\varepsilon}_{r}^{\mathrm{inc}}}}\xb7\stackrel{\u20d7}{E}\left(\stackrel{\u20d7}{r}\right)={\varepsilon}_{0}\left[{\varepsilon}_{t}^{\mathrm{inc}}\underline{\underline{{I}_{t}}}+{\varepsilon}_{c}^{\mathrm{inc}}\widehat{c}\widehat{c}\right]\xb7\stackrel{\u20d7}{E}\left(\stackrel{\u20d7}{r}\right).$$## Eq. 14

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \begin{array}{*{20}c} {\vec D ({\vec r }) = \varepsilon _0 \underline{\underline {\varepsilon _r^{{\rm MG}} }} \cdot \vec E ({\vec r }) = \varepsilon _0 \big[ {\varepsilon _t^{{\rm MG}} \underline{\underline {I_t }} + \varepsilon _c^{{\rm MG}} \hat c\hat c} \big] \cdot \vec E ({\vec r })} \\[12pt] \hspace{-5pt}{\vec D ({\vec r }) = \varepsilon _0 \underline{\underline {\varepsilon _r^{{\rm BR}} }} \cdot \vec E ({\vec r }) = \varepsilon _0 \big[ {\varepsilon _t^{{\rm BR}} \underline{\underline {I_t }} + \varepsilon _c^{{\rm BR}} \hat c\hat c} \big] \cdot \vec E ({\vec r })} \\ \end{array}. \end{equation}\end{document} $$\begin{array}{c}\stackrel{\u20d7}{D}\left(\stackrel{\u20d7}{r}\right)={\varepsilon}_{0}\underline{\underline{{\varepsilon}_{r}^{\mathrm{MG}}}}\xb7\stackrel{\u20d7}{E}\left(\stackrel{\u20d7}{r}\right)={\varepsilon}_{0}\left[{\varepsilon}_{t}^{\mathrm{MG}}\underline{\underline{{I}_{t}}}+{\varepsilon}_{c}^{\mathrm{MG}}\widehat{c}\widehat{c}\right]\xb7\stackrel{\u20d7}{E}\left(\stackrel{\u20d7}{r}\right)\\ \stackrel{\u20d7}{D}\left(\stackrel{\u20d7}{r}\right)={\varepsilon}_{0}\underline{\underline{{\varepsilon}_{r}^{\mathrm{BR}}}}\xb7\stackrel{\u20d7}{E}\left(\stackrel{\u20d7}{r}\right)={\varepsilon}_{0}\left[{\varepsilon}_{t}^{\mathrm{BR}}\underline{\underline{{I}_{t}}}+{\varepsilon}_{c}^{\mathrm{BR}}\widehat{c}\widehat{c}\right]\xb7\stackrel{\u20d7}{E}\left(\stackrel{\u20d7}{r}\right)\end{array}.$$The permeability tensor and the magnetic field [ [TeX:] $\vec B \left({\vec r } \right)$ $\stackrel{\u20d7}{B}\left(\stackrel{\u20d7}{r}\right)$ ] can be represented similarly.

The mathematical challenge expresses itself in calculating the polarizability dyadic for the inclusions. We first approach the case of a single electrically small inclusion, occupying the geometrical shape defined by the closed geometric shape [TeX:] $\vec r \in \{ {\vec {r}_{{\rm inc}} } \}$ $\stackrel{\u20d7}{r}\in \left\{{\stackrel{\u20d7}{r}}_{\mathrm{inc}}\right\}$ , embedded in the host medium occupying the space [TeX:] $\vec r \notin \{ {\vec {r}_{{\rm inc} } } \}$ $\stackrel{\u20d7}{r}\notin \left\{{\stackrel{\u20d7}{r}}_{\mathrm{inc}}\right\}$ . The time-harmonic version Maxwell curl postulates can be written everywhere as:

## Eq. 15

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \begin{array}{*{20}c} {\vec \nabla \times \vec E ({\vec r }) + i\omega \mu _0 \cdot \vec H ({\vec r }) = - \vec {\rm K} _{{\rm eq}} ({\vec r })} \\[12pt] {\vec \nabla \times \vec H ({\vec r }) - i\omega \mu _0 \cdot \vec E ({\vec r }) = \vec J _{{\rm eq}} ({\vec r })} \\ \end{array}. \end{equation}\end{document} $$\begin{array}{c}\stackrel{\u20d7}{\nabla}\times \stackrel{\u20d7}{E}\left(\stackrel{\u20d7}{r}\right)+i\omega {\mu}_{0}\xb7\stackrel{\u20d7}{H}\left(\stackrel{\u20d7}{r}\right)=-{\stackrel{\u20d7}{\mathrm{K}}}_{\mathrm{eq}}\left(\stackrel{\u20d7}{r}\right)\\ \stackrel{\u20d7}{\nabla}\times \stackrel{\u20d7}{H}\left(\stackrel{\u20d7}{r}\right)-i\omega {\mu}_{0}\xb7\stackrel{\u20d7}{E}\left(\stackrel{\u20d7}{r}\right)={\stackrel{\u20d7}{J}}_{\mathrm{eq}}\left(\stackrel{\u20d7}{r}\right)\end{array}.$$The source current densities K_{eq}(*r*) = 0 and *J*
_{eq}(*r*) = 0 everywhere except inside the inclusions:

## Eq. 16

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \begin{array}{*{20}c} {\vec J _{{\rm eq}} ({\vec r }) = - i\omega \varepsilon _0 \big({\underline{\underline {\varepsilon _r^h }} - \underline{\underline {\varepsilon _r^{{\rm inc}} }} } \big) \cdot \vec E ({\vec r })} \\[12pt] {\vec {\rm K} _{{\rm eq}} ({\vec r }) = i\omega \mu _0 \big({\underline{\underline {\mu _r^h }} - \underline{\underline {\mu _r^{{\rm inc}} }} } \big) \cdot \vec H ({\vec r })} \\ \end{array}. \end{equation}\end{document} $$\begin{array}{c}{\stackrel{\u20d7}{J}}_{\mathrm{eq}}\left(\stackrel{\u20d7}{r}\right)=-i\omega {\varepsilon}_{0}\left(\underline{\underline{{\varepsilon}_{r}^{h}}}-\underline{\underline{{\varepsilon}_{r}^{\mathrm{inc}}}}\right)\xb7\stackrel{\u20d7}{E}\left(\stackrel{\u20d7}{r}\right)\\ {\stackrel{\u20d7}{\mathrm{K}}}_{\mathrm{eq}}\left(\stackrel{\u20d7}{r}\right)=i\omega {\mu}_{0}\left(\underline{\underline{{\mu}_{r}^{h}}}-\underline{\underline{{\mu}_{r}^{\mathrm{inc}}}}\right)\xb7\stackrel{\u20d7}{H}\left(\stackrel{\u20d7}{r}\right)\end{array}.$$As was presented by Ref. 32 and following Ref. 56, the solution of Eq. 15 for both the electric and magnetic fields may be stated as in a similar form to that of EM field in scattering problems; composed of a term describing the field in the absence of inclusions and spatial integrals on the inclusion volume of the products of the source current densities and the matching Green functions as presented in detail in Refs. 32 and 56.

In an effort for estimating the field inside the inclusions, the homogenization theory often makes use of the fact that the inclusions are electrically small (small relative to the fields wavelength), and assume that the field is uniform inside the inclusion. In the Rayleigh approximation, only the terms of order *R*
^{−3} in the dyadic Green functions contribute to the above-mentioned volume integrals and can be written in approximation in a simpler manner, considerably simplifying the calculations. The mathematical process ends with the following approximated results for the field inside the inclusions:

## Eq. 17

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \begin{array}{*{20}c} {\vec E (0) = \vec E ^h (0) + i\omega \mu _0 \displaystyle\frac{1}{{k_0 ^2 }}q_e^h \underline{\underline {L_e^h }} \cdot \big({\underline{\underline {\varepsilon _r^h }} } \big)^{ - 1} \cdot \vec J _{{\rm eq}} \left(0 \right)} \hfill \\[20pt] {\vec H (0) = \vec H ^h (0) + i\omega \varepsilon _0 \displaystyle\frac{1}{{k_0 ^2 }}q_m^h \underline{\underline {L_m^h }} \cdot \big({\underline{\underline {\mu _r^h }} } \big)^{ - 1} \cdot \vec K _{{\rm eq}} (0)} \hfill \\[10pt] \end{array}, \end{equation}\end{document} $$\begin{array}{c}{\displaystyle \stackrel{\u20d7}{E}\left(0\right)={\stackrel{\u20d7}{E}}^{h}\left(0\right)+i\omega {\mu}_{0}\frac{1}{{k}_{0}^{2}}{q}_{e}^{h}\underline{\underline{{L}_{e}^{h}}}\xb7{\left(\underline{\underline{{\varepsilon}_{r}^{h}}}\right)}^{-1}\xb7{\stackrel{\u20d7}{J}}_{\mathrm{eq}}\left(0\right)}\hfill \\ {\displaystyle \stackrel{\u20d7}{H}\left(0\right)={\stackrel{\u20d7}{H}}^{h}\left(0\right)+i\omega {\varepsilon}_{0}\frac{1}{{k}_{0}^{2}}{q}_{m}^{h}\underline{\underline{{L}_{m}^{h}}}\xb7{\left(\underline{\underline{{\mu}_{r}^{h}}}\right)}^{-1}\xb7{\stackrel{\u20d7}{K}}_{\mathrm{eq}}\left(0\right)}\hfill \end{array},$$^{57}

## Eq. 18

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \begin{array}{*{20}c} {\underline{\underline {L_e^h }} = \displaystyle\frac{{\underline{\underline {I_t }} }}{{2q_e^h \sqrt {1 + \displaystyle\frac{{4q_e^h }}{{u^2 }}} }} + \displaystyle\frac{1}{{q_e^h }}\left({1 - \displaystyle\frac{1}{{\sqrt {1 + \displaystyle\frac{{4q_e^h }}{{u^2 }}} }}} \right)\hat c\hat c} \\[12pt] {\underline{\underline {L_m^h }} = \displaystyle\frac{{\underline{\underline {I_t }} }}{{2q_m^h \sqrt {1 + \displaystyle\frac{{4q_m^h }}{{u^2 }}} }} + \displaystyle\frac{1}{{q_m^h }}\left({1 - \displaystyle\frac{1}{{\sqrt {1 + \displaystyle\frac{{4q_m^h }}{{u^2 }}} }}} \right)\hat c\hat c} \end{array} \end{equation}\end{document} $$\begin{array}{c}{\displaystyle \underline{\underline{{L}_{e}^{h}}}=\frac{\underline{\underline{{I}_{t}}}}{2{q}_{e}^{h}\sqrt{{\displaystyle 1+\frac{4{q}_{e}^{h}}{{u}^{2}}}}}+\frac{1}{{q}_{e}^{h}}\left({\displaystyle 1-\frac{1}{\sqrt{{\displaystyle 1+\frac{4{q}_{e}^{h}}{{u}^{2}}}}}}\right)\widehat{c}\widehat{c}}\\ {\displaystyle \underline{\underline{{L}_{m}^{h}}}=\frac{\underline{\underline{{I}_{t}}}}{2{q}_{m}^{h}\sqrt{{\displaystyle 1+\frac{4{q}_{m}^{h}}{{u}^{2}}}}}+\frac{1}{{q}_{m}^{h}}\left({\displaystyle 1-\frac{1}{\sqrt{{\displaystyle 1+\frac{4{q}_{m}^{h}}{{u}^{2}}}}}}\right)\widehat{c}\widehat{c}}\end{array}$$*u*stands for the relation between the cylinders` height and radius. Using Eqs. 16, 17 we finally obtain:

## Eq. 19

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \begin{array}{*{20}c} {\vec J _{{\rm eq}} (0) = - i\omega \varepsilon _0 \big[ {\underline{\underline {\varepsilon _r^h }} - \underline{\underline {\varepsilon _r^{{\rm inc}} }} } \big] \cdot \big\{ {\underline{\underline I} - q_e^h \underline{\underline {L_e^h }} \cdot \big({\underline{\underline {\varepsilon _r^h }} } \big)^{ - 1} \cdot \big[ {\underline{\underline {\varepsilon _r^h }} - \underline{\underline {\varepsilon _r^{{\rm inc}} }} } \big]} \big\}^{ - 1} \cdot \vec E ^h (0)} \hfill \\[20pt] {\vec K _{{\rm eq}} (0) = - i\omega \mu _0 \big[ {\underline{\underline {\mu _r^h }} - \underline{\underline {\mu _r^{{\rm inc}} }} } \big] \cdot \big\{ {\underline{\underline I} - q_m^h \underline{\underline {L_m^h }} \cdot \big({\underline{\underline {\mu _r^h }} } \big)^{ - 1} \cdot \big[ {\underline{\underline {\mu _r^h }} - \underline{\underline {\mu _r^{{\rm inc}} }} } \big]} \big\}^{ - 1} \cdot \vec H ^h (0)} \hfill \\ \end{array}. \end{equation}\end{document} $$\begin{array}{c}{\stackrel{\u20d7}{J}}_{\mathrm{eq}}\left(0\right)=-i\omega {\varepsilon}_{0}\left[\underline{\underline{{\varepsilon}_{r}^{h}}}-\underline{\underline{{\varepsilon}_{r}^{\mathrm{inc}}}}\right]\xb7{\left\{\underline{\underline{I}}-{q}_{e}^{h}\underline{\underline{{L}_{e}^{h}}}\xb7{\left(\underline{\underline{{\varepsilon}_{r}^{h}}}\right)}^{-1}\xb7\left[\underline{\underline{{\varepsilon}_{r}^{h}}}-\underline{\underline{{\varepsilon}_{r}^{\mathrm{inc}}}}\right]\right\}}^{-1}\xb7{\stackrel{\u20d7}{E}}^{h}\left(0\right)\hfill \\ {\stackrel{\u20d7}{K}}_{\mathrm{eq}}\left(0\right)=-i\omega {\mu}_{0}\left[\underline{\underline{{\mu}_{r}^{h}}}-\underline{\underline{{\mu}_{r}^{\mathrm{inc}}}}\right]\xb7{\left\{\underline{\underline{I}}-{q}_{m}^{h}\underline{\underline{{L}_{m}^{h}}}\xb7{\left(\underline{\underline{{\mu}_{r}^{h}}}\right)}^{-1}\xb7\left[\underline{\underline{{\mu}_{r}^{h}}}-\underline{\underline{{\mu}_{r}^{\mathrm{inc}}}}\right]\right\}}^{-1}\xb7{\stackrel{\u20d7}{H}}^{h}\left(0\right)\hfill \end{array}.$$The inclusion, with volume *V*
_{inc}, can be thought of as an ensemble of electric dipole and magnetic dipole:

## Eq. 20

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \begin{array}{*{20}c} {\vec p _{{\rm eq}} \left(0 \right) = \displaystyle\frac{i}{\omega }V_{{\rm inc}} \vec J _{{\rm eq}} \left(0 \right) = \underline{\underline \pi } _{ee} \cdot \vec E ^h \left(0 \right)} \hfill \\[8pt] {\vec m _{{\rm eq}} \left(0 \right) = \displaystyle\frac{i}{\omega }V_{{\rm inc}} \vec K _{{\rm eq}} \left(0 \right) = \underline{\underline \pi } _{mm} \cdot \vec H ^h \left(0 \right)} \hfill \\[10pt] \end{array}, \end{equation}\end{document} $$\begin{array}{c}{\displaystyle {\stackrel{\u20d7}{p}}_{\mathrm{eq}}\left(0\right)=\frac{i}{\omega}{V}_{\mathrm{inc}}{\stackrel{\u20d7}{J}}_{\mathrm{eq}}\left(0\right)={\underline{\underline{\pi}}}_{ee}\xb7{\stackrel{\u20d7}{E}}^{h}\left(0\right)}\hfill \\ {\displaystyle {\stackrel{\u20d7}{m}}_{\mathrm{eq}}\left(0\right)=\frac{i}{\omega}{V}_{\mathrm{inc}}{\stackrel{\u20d7}{K}}_{\mathrm{eq}}\left(0\right)={\underline{\underline{\pi}}}_{mm}\xb7{\stackrel{\u20d7}{H}}^{h}\left(0\right)}\hfill \end{array},$$where the desired dielectric and magnetic polarizability dyadics of the cylindrical inclusions in the host medium are expressed by:

## Eq. 21

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \begin{array}{*{20}c} {\underline{\underline \pi } _{ee} = - V_{{\rm inc}} \varepsilon _0 \big[ {\underline{\underline {\varepsilon _r^h }} - \underline{\underline {\varepsilon _r^{{\rm inc}} }} } \big] \cdot \big\{ {\underline{\underline I} - q_e^h \underline{\underline {L_e^h }} \cdot \big({\underline{\underline {\varepsilon _r^h }} } \big)^{ - 1} \cdot \big[ {\underline{\underline {\varepsilon _r^h }} - \underline{\underline {\varepsilon _r^{{\rm inc}} }} } \big]} \big\}^{ - 1} } \hfill \\[24pt] {\underline{\underline \pi } _{mm} = - V_{{\rm inc}} \varepsilon _0 \big[ {\underline{\underline {\mu _r^h }} - \underline{\underline {\mu _r^{{\rm inc}} }} } \big] \cdot \big\{ {\underline{\underline I} - q_m^h \underline{\underline {L_m^h }} \cdot \big({\underline{\underline {\mu _r^h }} } \big)^{ - 1} \cdot \big[ {\underline{\underline {\mu _r^h }} - \underline{\underline {\mu _r^{{\rm inc}} }} } \big]} \big\}^{ - 1} } \hfill \\ \end{array}. \end{equation}\end{document} $$\begin{array}{c}{\underline{\underline{\pi}}}_{ee}=-{V}_{\mathrm{inc}}{\varepsilon}_{0}\left[\underline{\underline{{\varepsilon}_{r}^{h}}}-\underline{\underline{{\varepsilon}_{r}^{\mathrm{inc}}}}\right]\xb7{\left\{\underline{\underline{I}}-{q}_{e}^{h}\underline{\underline{{L}_{e}^{h}}}\xb7{\left(\underline{\underline{{\varepsilon}_{r}^{h}}}\right)}^{-1}\xb7\left[\underline{\underline{{\varepsilon}_{r}^{h}}}-\underline{\underline{{\varepsilon}_{r}^{\mathrm{inc}}}}\right]\right\}}^{-1}\hfill \\ {\underline{\underline{\pi}}}_{mm}=-{V}_{\mathrm{inc}}{\varepsilon}_{0}\left[\underline{\underline{{\mu}_{r}^{h}}}-\underline{\underline{{\mu}_{r}^{\mathrm{inc}}}}\right]\xb7{\left\{\underline{\underline{I}}-{q}_{m}^{h}\underline{\underline{{L}_{m}^{h}}}\xb7{\left(\underline{\underline{{\mu}_{r}^{h}}}\right)}^{-1}\xb7\left[\underline{\underline{{\mu}_{r}^{h}}}-\underline{\underline{{\mu}_{r}^{\mathrm{inc}}}}\right]\right\}}^{-1}\hfill \end{array}.$$## 6.3.

### Bruggeman Formalism

We now turn to implement the Bruggeman homogenization formalism. Assuming nonmagnetic materials, only
[TeX:]
$\underline{\underline \pi } _{ee}$
${\underline{\underline{\pi}}}_{ee}$
needs to be treated. On the contrary to the Maxwell Garnett formalism (the reader is referred to Ref. 32 for a detailed mathematical description) the Bruggeman formalism does not distinguish between the host and inclusion media in the sense of treating both as if composed of electrically small “particles” (in a similar manner described above for the inclusions alone) in appropriate proportions *f*
_{h}, *f*
_{inc}.^{32, 58} The motivation of homogenization, as mentioned above, is to present the electric displacement field (and the magnetic field) in the entire media by a single equation for each in a similar form as a homogeneous media—in essence characterizing an equivalent homogeneous media [Eq. 14]. When treating both the host and inclusions as particulate materials one can achieve this by simply weighting each of the host and inclusions contribution to the polarization of the media via the polarizability density dyadic for each of them:

## Eq. 22

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} &&{f_h \cdot \underline{\underline \pi } ^{h }_{ee, BR} + f_{{\rm inc}} \cdot \underline{\underline \pi } ^{\rm inc} _{ee, BR} = \underline{\underline 0} } \nonumber \\ [12pt] &&\underline{\underline \pi } ^{X }_{ee, BR} = - \big[ {\underline{\underline {\varepsilon _r^X }} - \underline{\underline {\varepsilon _r^{{\rm BR}} }} } \big]\nonumber\\ [12pt] &&\qquad \cdot\, \big\{ {\underline{\underline I} - q_e^X \underline{\underline {L_e^X }} \cdot \big({\underline{\underline {\varepsilon _r^X }} } \big)^{ - 1} \cdot \big[ {\underline{\underline {\varepsilon _r^X }} - \underline{\underline {\varepsilon _r^{{\rm BR}} }} } \big]} \big\}^{ - 1}. \end{eqnarray}\end{document} $$\begin{array}{ccc}& & {f}_{h}\xb7{\underline{\underline{\pi}}}_{ee,BR}^{h}+{f}_{\mathrm{inc}}\xb7{\underline{\underline{\pi}}}_{ee,BR}^{\mathrm{inc}}=\underline{\underline{0}}\hfill \\ & & {\underline{\underline{\pi}}}_{ee,BR}^{X}=-\left[\underline{\underline{{\varepsilon}_{r}^{X}}}-\underline{\underline{{\varepsilon}_{r}^{\mathrm{BR}}}}\right]\hfill \\ & & \phantom{\rule{2.em}{0ex}}\xb7\phantom{\rule{0.16em}{0ex}}{\left\{\underline{\underline{I}}-{q}_{e}^{X}\underline{\underline{{L}_{e}^{X}}}\xb7{\left(\underline{\underline{{\varepsilon}_{r}^{X}}}\right)}^{-1}\xb7\left[\underline{\underline{{\varepsilon}_{r}^{X}}}-\underline{\underline{{\varepsilon}_{r}^{\mathrm{BR}}}}\right]\right\}}^{-1}.\hfill \end{array}$$*X*stands for either “

*h*” for host material and “inc” for inclusion. The host and inclusions “particles” do not necessarily need to be of the same shape, resulting in different form for the polarizability dyadics. For example, we could have chosen to treat the host particles as if being spherical. We set both host and inclusions as cylindrical particles driven by the motivation to “decrease” the gaps between particles in the composite materials. Though some difference can arise from another choice, its overall influence on the entire model should be minimal.

The relative permittivity dyadic of each homogenized epidermis layer
[TeX:]
$\underline{\underline {\varepsilon ^{{\rm BR}} }}$
$\underline{\underline{{\varepsilon}^{\mathrm{BR}}}}$
can be extracted from the nonlinear Eq. 22 by numerical means. We use the iterative scheme
[TeX:]
$\underline{\underline {\varepsilon ^{{\rm BR}} }} ^{\left(j \right)} = \Im \{ {\underline{\underline {\varepsilon ^{{\rm BR}} }} ^{\left({j - 1} \right)} } \},\left({j = 1,2,.\;.} \right)$
${\underline{\underline{{\varepsilon}^{\mathrm{BR}}}}}^{\left(j\right)}=\Im \left\{{\underline{\underline{{\varepsilon}^{\mathrm{BR}}}}}^{\left(j-1\right)}\right\},\left(j=1,2,.\phantom{\rule{0.28em}{0ex}}.\right)$
^{31} with 10 iterations and the initial value:
[TeX:]
$\underline{\underline {\varepsilon ^{{\rm BR}} }} ^{\left(0 \right)} = f_h \underline{\underline \varepsilon } ^h + f_{{\rm inc}} \underline{\underline \varepsilon } ^{{\rm inc}}$
${\underline{\underline{{\varepsilon}^{\mathrm{BR}}}}}^{\left(0\right)}={f}_{h}{\underline{\underline{\varepsilon}}}^{h}+{f}_{\mathrm{inc}}{\underline{\underline{\varepsilon}}}^{\mathrm{inc}}$
, where:

## Eq. 23

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \begin{array}{l} \Im \{ {\underline{\underline {\varepsilon ^{{\rm BR}} }} } \} = \{ {f_h \underline{\underline \varepsilon } ^h [ {\underline{\underline I} + \underline{\underline L} ^h _e \cdot ({\underline{\underline \varepsilon } ^h - \underline{\underline {\varepsilon ^{{\rm BR}} }} })} ]^{ - 1} + f_{{\rm inc}} \underline{\underline \varepsilon } ^{{\rm inc}} [ {\underline{\underline I} + \underline{\underline L} ^{{\rm inc}} _e \cdot ({\underline{\underline \varepsilon } ^{{\rm inc}} - \underline{\underline {\varepsilon ^{{\rm BR}} }} })} ]^{ - 1} } \} \\[12pt] \cdot\, \{ {f_h [ {\underline{\underline I} + \underline{\underline {L_e^h }} \cdot ({\underline{\underline \varepsilon } ^h - \underline{\underline {\varepsilon ^{{\rm BR}} }} })} ]^{ - 1} + f_{{\rm inc}} [ {\underline{\underline I} + \underline{\underline {L_e^{{\rm inc}} }} \cdot ({\underline{\underline \varepsilon } ^{{\rm inc}} - \underline{\underline {\varepsilon ^{{\rm BR}} }} })} ]^{ - 1} } \}. \\ \end{array} \end{equation}\end{document} $$\begin{array}{c}\Im \left\{\underline{\underline{{\varepsilon}^{\mathrm{BR}}}}\right\}=\left\{{f}_{h}{\underline{\underline{\varepsilon}}}^{h}{\left[\underline{\underline{I}}+{\underline{\underline{L}}}_{e}^{h}\xb7\left({\underline{\underline{\varepsilon}}}^{h}-\underline{\underline{{\varepsilon}^{\mathrm{BR}}}}\right)\right]}^{-1}+{f}_{\mathrm{inc}}{\underline{\underline{\varepsilon}}}^{\mathrm{inc}}{\left[\underline{\underline{I}}+{\underline{\underline{L}}}_{e}^{\mathrm{inc}}\xb7\left({\underline{\underline{\varepsilon}}}^{\mathrm{inc}}-\underline{\underline{{\varepsilon}^{\mathrm{BR}}}}\right)\right]}^{-1}\right\}\hfill \\ \xb7\phantom{\rule{0.16em}{0ex}}\left\{{f}_{h}{\left[\underline{\underline{I}}+\underline{\underline{{L}_{e}^{h}}}\xb7\left({\underline{\underline{\varepsilon}}}^{h}-\underline{\underline{{\varepsilon}^{\mathrm{BR}}}}\right)\right]}^{-1}+{f}_{\mathrm{inc}}{\left[\underline{\underline{I}}+\underline{\underline{{L}_{e}^{\mathrm{inc}}}}\xb7\left({\underline{\underline{\varepsilon}}}^{\mathrm{inc}}-\underline{\underline{{\varepsilon}^{\mathrm{BR}}}}\right)\right]}^{-1}\right\}.\hfill \end{array}$$As described in Sec. 2.2.2, the helical sweat ducts are treated as if being composed of a collection of very small sequential finite cylinders placed along the center line of the helix [Fig. 3a]. Each of the cylinders is set to have similar electrical characteristics as the sweat ducts. The epidermis is divided into sublayers, each containing finite cylinders all pointing at a common orientation. Each sub-layer is then homogenized based on the orientation of the cylinders embedded within it, yielding a homogeneous material but one that is described by an anisotropic dielectric tensor. Thus we are left with an anisotropic but homogeneous, multilayer structure of the epidermis [Fig. 3b].

In order to simplify the calculation, we use the Bruggeman formalism to deduce the effective permittivity dyadic of a single planar sublayer containing cylinders at a concentration of *f* = 0.3%, pointing at the *z* direction, resulting in a diagonal–uniaxial permittivity tensor with the principal axis oriented along the *z* direction. The permittivity dyadic of the effective sublayers differ by the appropriate orientation of their embedded cylinders defined by the polar and azimuth angles θ and ϕ. The permittivity dyadic for each sublayer in the *xyz* frame can be determined by rotating the permittivity dyadic of the homogenized diagonal–uniaxial permittivity tensor to the orientation (θ,ϕ) of the appropriate slice. The permittivity dyadic for a sublayer containing cylinders pointing at polar and azimuth angles θ and ϕ is given by the general expression:^{33}

## Eq. 24

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{eqnarray} \varepsilon (\theta,\phi) = \left({ \begin{array}{ccc} {\varepsilon _2 + \delta \cos ^2 \phi } &\quad {0.5\delta \sin 2\phi } &\quad {0.5\left({\varepsilon _3 - \varepsilon _1 } \right)\sin 2\theta \cos \phi } \\ {0.5\delta \sin 2\phi } &\quad {\varepsilon _2 + \delta \sin ^2 \phi } &\quad {0.5\left({\varepsilon _3 - \varepsilon _1 } \right)\sin 2\theta \sin \phi } \\ {0.5\left({\varepsilon _3 - \varepsilon _1 } \right)\sin 2\theta \cos \phi } &\quad {0.5\left({\varepsilon _3 - \varepsilon _1 } \right)\sin 2\theta \sin \phi } &\quad {\varepsilon _1 + \left({\varepsilon _3 - \varepsilon _1 } \right)\cos ^2 \theta } \\ \end{array}} \right). \end{eqnarray} \end{document} $$\begin{array}{c}\hfill \varepsilon (\theta ,\phi )=\left(\begin{array}{ccc}{\varepsilon}_{2}+\delta {\mathrm{cos}}^{2}\phi & \phantom{\rule{1em}{0ex}}0.5\delta \mathrm{sin}2\phi & \phantom{\rule{1em}{0ex}}0.5\left({\varepsilon}_{3}-{\varepsilon}_{1}\right)\mathrm{sin}2\theta \mathrm{cos}\phi \\ 0.5\delta \mathrm{sin}2\phi & \phantom{\rule{1em}{0ex}}{\varepsilon}_{2}+\delta {\mathrm{sin}}^{2}\phi & \phantom{\rule{1em}{0ex}}0.5\left({\varepsilon}_{3}-{\varepsilon}_{1}\right)\mathrm{sin}2\theta \mathrm{sin}\phi \\ 0.5\left({\varepsilon}_{3}-{\varepsilon}_{1}\right)\mathrm{sin}2\theta \mathrm{cos}\phi & \phantom{\rule{1em}{0ex}}0.5\left({\varepsilon}_{3}-{\varepsilon}_{1}\right)\mathrm{sin}2\theta \mathrm{sin}\phi & \phantom{\rule{1em}{0ex}}{\varepsilon}_{1}+\left({\varepsilon}_{3}-{\varepsilon}_{1}\right){\mathrm{cos}}^{2}\theta \end{array}\right).\end{array}$$_{1}, ε

_{2}, ε

_{3}are the resulting principal values of the dielectric tensor after homogenization, and δ = ε

_{1}cos

^{2}θ + ε

_{3}sin

^{2}θ − ε

_{2}is an effective local dielectric anisotropy. In our uniaxial case ε

_{3}is the principal dielectric constant along the cylinder axis while along its radius we have ε

_{1}= ε

_{2}. By implementing this approach, the epidermis embedded with the helical inclusions is simplified to a stack of dielectric anisotropic (uniaxial) layers, a problem well investigated.

^{34}

## Acknowledgments

The authors are grateful to Profofessor Akhlesh Lakhtakia for useful discussions on the homogenization procedure.

## References

*p-p*+silicon homojunction determined by terahertz and midinfrared spectroscopic ellipsometry,” Appl. Phys. Lett., 95 032102 (2009). https://doi.org/10.1063/1.3184567 Google Scholar

*In vivo*study of human skin using pulsed terahertz radiation,” Phys. Med. Biol., 49 1595 –1607 (2004). https://doi.org/10.1088/0031-9155/49/9/001 Google Scholar